\section*{Introduction}
\quad \ HOTIs are topological crystalline insulators that have protected  topological  features on the boundaries with co-dimension greater than one. In contrast to (1st order) three-dimensional time-reversal protected topological insulators that show an insulating bulk with conducting surface modes, three-dimensional 2nd order HOTIs show insulating bulk and surfaces, but with one-dimensional conducting channels at hinges (i.e., the intersection of two surface facets) or equivalent crystal defects  \cite{schindler2018higher,song2017d,benalcazar2018prb,Schindlereaat0346,Olsen2020gaplss,zhu2020quantized}.  While several groups have recently succeeded in observing HOTI phases in  electronic circuit, phononic, and photonic systems \cite{serra2018observation,peterson2018quantized,noh2018topological,imhof2018topolectrical}, the experimental realization and characterization of HOTIs in solid state materials has proven more challenging.  Among the  proposed HOTI materials \cite{Schindlereaat0346,choi2020evidence,lee2020two,zhang2020mobius,wang2019higher,zhang2019higher,park2019higher}, bulk Bi has recently emerged as a key candidate. Bi and Bi$_{1-x}$Sb$_x$ alloys are well-known important topological materials with highly tunable electronic properties.  Scanning tunneling microscopy studies \cite{schindler2018higher,drozdov2014one} have observed signatures of one-dimensional modes on the edges of (111) facets of Bi, but the higher order topological origin of these edge modes remains to be unequivocally established. Indeed an earlier study modeled the (111) facet as a free-standing Bi-bilayer \cite{drozdov2014one}, which is predicted to be a quantum spin Hall insulator, and hence attributed the observed hinge modes to edge states of a quantum spin Hall system. In contrast, a recent study \cite{schindler2018higher} ascribed the metallic edge modes to the manifestation of a bulk HOTI protected by crystalline and time-reversal symmetries.  New experimental data is necessary to distinguish between these scenarios. Importantly, the confirmation of Bi as a HOTI  would indicate that the low-energy bands of bismuth must have trivial time-reversal protected $\mathbbm{Z}_{2}$ topology \cite{schindler2018higher,hsu2019topology}. This not only has implications for the topology of bulk Bi, which still remains controversial, but also the family of Bi$_{1-x}$Sb$_{x}$ alloys.

\section*{Results}
\textbf{STM measurements on Bi and Bi$_{0.92}$Sb$_{0.08}$ films.}
Bi naturally cleaves along the (111) plane, which makes STM studies of other planes, e.g, (110), difficult in bulk cleaved crystals.  To circumvent this problem, we use molecular beam epitaxy to grow Bi and Bi$_{0.92}$Sb$_{0.08}$ films on n-doped silicon (111) wafers  \cite{takayama2015one,kawakami2015one,nagao2004nanofilm,PhysRevB.91.075429,2020APExp..13h5506N}. The films were grown in a custom MBE system (described in more detail in the Methods section), and the composition was confirmed using Rutherford Backscattering Spectrometry (RBS) (see Supplementary Note 4). Film thicknesses were chosen to be in the range of 30-90 BLs to ensure that they are representative of the bulk.

To establish a benchmark for comparison, we first reproduce the edge state data previously obtained on Bi(111) crystals using our films. STM topographies on our Bi(111)/Si(111) films reveal triangular islands and a hexagonal lattice indicative of (111) facets. The dI/dV spectra on the bulk and edges of the triangular islands show signatures of the Van Hove singularity associated with the surface states, as well as the emergence of clear 1D modes at type-A edges (see Supplementary Note 1), as was observed in the cleaved bulk Bi crystals \cite{drozdov2014one}. Having confirmed the presence of edge modes in the (111) films, we move on to the (110) films.

By changing the growth conditions (see Methods section), we were able to grow Bi and Bi$_{0.92}$Sb$_{0.08}$ in the (110) orientation on the same Si(111) substrates. STM images of the (110) oriented films of both Bi and Bi$_{0.92}$Sb$_{0.08}$ are shown in Figure \ref{f1}. Large area topographic images of bismuth films (Figure \ref{f1} (b)) show that the surface is uniformly covered with (110) oriented rectangular islands with longer edges parallel to the ($1\bar{1}0$) direction. Figure \ref{f1} (e) shows an STM image with its corresponding Fourier transform in the inset. Height profiles along line cuts (red and blue color lines in Figure \ref{f1} (e)) confirm the lattice constants of the pseudo-cubic crystal structure of Bi(110), i.e., 0.443 nm and 0.471 nm in Figure \ref{f1} (g) and (h), respectively. Representative spectra at the center of the Bi and Bi$_{0.92}$Sb$_{0.08}$ (110) islands, shown in Figure \ref{f1} (i) and (j), reveal a broad increase in the density of states (DOS) above 130 mV in both samples. This increase can be ascribed to the spin-split Rashba surface states known to exist in bismuth and its alloys \cite{takayama2015one,takayama2014rashba}.
\\

\noindent\textbf{Spectroscopic measurements on a rectangular facet of Bi(110).}
We first focus on the edge states observed on rectangular islands of the Bi (110) films. Spectra on the edge of an island show a sharp peak in the DOS (density of states)  as shown in Figure \ref{f2} (c). This singularity corresponds to bound modes which are localized within a few nanometers of the edge, as seen in the spectroscopic measurements across and along the edge shown in Figure \ref{f2} (d) and Figure \ref{f2} (e) respectively. This sharp peak resembles the DOS feature observed on the edges of Bi(111). To further characterize this edge mode, we obtain spectroscopic data on all four edges of a (110) oriented rectangular island as shown in Figure \ref{f3} (b,c,d,e). Intriguingly we find that only three out of four edges show clear edge modes, while the fourth edge shows no peaks. We confirmed this observation by measurements on different rectangular facet islands on multiple samples. In all cases, the edge that shows no edge modes is  the one that is perpendicular to the ($1\bar{1}0$) direction. Crucially, the observation of edge states on three out of the four edges cannot be attributed to different bonding on just one side, which was a key part of the argument in Ref.~\onlinecite{drozdov2014one}.  By looking at the lattice structure of the ($110$) facet, we know that the two short sides are identical from a bonding perspective (see Supplementary Note 6). Yet, for a given island, one short side shows hinge states while the other does not.
\\

\noindent\textbf{Theoretical Justification.}
To understand the novel pattern of rectangular island edge modes on the (110) facets, let us first briefly review the recent developments for the topology of bismuth. Specifically, Refs.~\onlinecite{schindler2018higher} and \onlinecite{hsu2019topology} argued that the bulk topology of bismuth can be described as two copies of a strong 3D TI, i.e., two band inversions at the $T$-point.
  Since there are an even number of strong TI copies, they conclude that the low-energy bulk bands of bismuth have trivial strong $\mathbbm{Z}_{2}$ topology,  but may support two robust Dirac points on high-symmetry surface facets exhibiting additional crystalline symmetries. For example, the (111) facet has $\hat{C}_{3}$ rotation symmetry; the two Dirac points are located at the $T$-point and belong to different $\hat{C}_{3}$ eigenvalue sectors having $\hat{C}_{3}$ eigenvalues $-1$ and $\exp(\pm i\pi/3)$ respectively \cite{schindler2018higher}. If the two Dirac points are pinned at the same energy, then they cannot be gapped out while preserving the $\hat{C}_{3}$ (no inter-cone gap) and  time-reversal symmetries (no intra-cone gap). However, there are no symmetry constraints that fix the energies of the two cones, and generically it is expected that the (111) surface will be gapped \cite{schindler2018higher}. In this context, Ref. \onlinecite{schindler2018higher} considered the case where the (111) facet is the top facet of a hexagonal crystal (as shown in Figure \ref{f1} (a)). For this geometry the $\hat{C}_{3}$ symmetry and time-reversal symmetry are compatible with sign-alternating surface-state masses/gaps on the side facets (this breaks $\hat{C}_{6}$ but preserves $\hat{C}_{3}$) \cite{schindler2018higher}. Since the surface Dirac cones on the top/bottom (111) facets can also be generically gapped while preserving the symmetries, this will lead to sharply localized hinge modes where the side surfaces intersect the top/bottom surfaces, which were first observed for bulk Bi crystals in Ref. \onlinecite{drozdov2014one,schindler2018higher}, and confirmed for Bi films in our measurements above.

We propose that the observed edge states of an island on the (110) surface can be ascribed to  $\hat{C}_{2}$ (of the bulk crystal around the $(1\bar{1}0)$ axis) and time-reversal symmetry protected higher order topology consistent with Ref. \onlinecite{hsu2019topology}. The (110) surfaces of our islands do not have any protective crystal symmetry and thus we generically expect the two surface Dirac cones on those surfaces (as predicted by Ref. \onlinecite{hsu2019topology})  to be gapped. When considering hinge states though, we must also consider the mass gaps for surface states on the side-surfaces of such islands. First, the two long sides/edges of the (110) rectangular islands do not have any protected crystal symmetry, hence we expect the surface states to be gapped. Second, the two short sides/edges of the (110) rectangular islands are truncated ($1\bar{1}0)$ surfaces, which could show remnants of the two gapless Dirac cone surface states.  Indeed, Ref. \onlinecite{hsu2019topology} recently emphasized that the ($1\bar{1}0$) surface has $\hat{C}_{2}$ rotation symmetry that protects the two, generically unpinned, surface Dirac points related by $\hat{C}_2$ rotation. In this case, the ($1\bar{1}0$) facet, cannot be gapped out (perturbatively) as long as the $\hat{C}_{2}$ and time-reversal symmetries are preserved since the Dirac cones do not coincide at a TRIM. However, since we are considering a surface island, translation symmetry is strongly broken on the ($1\bar{1}0$) facets. Hence, we do not expect the surface states to be robust since their momentum space locations become irrelevant.
From these arguments we expect that the surfaces that form the top and sides of the island can be generically gapped, so we can explore our observations in the context of higher order topology in this system.  In contrast to previous systems with crystalline-protected topology, (e.g., the mirror-protected surface Dirac cones in SnTe \cite{tanaka2012experimental}), the $(110)$ surface of Bi does not preserve the $\hat{C}_{2}$ bulk symmetry responsible for the higher order protection\cite{Schindlereaat0346}. However, this does not mean that we cannot observe signatures of higher order topology, it just implies that the possible hinge state patterns may be affected by the details of the surface termination. For a rectangular crystal as shown in Figure \ref{f3} (f), where ($1\bar{1}0$) is the top edge/facet, there are two inequivalent configurations of surface mass gaps that match our observations of hinge modes on the $(110)$ facet; i.e., that helical modes are localized on only three edges out of four. The configuration of helical hinge modes shown in the left panel of Figure \ref{f3} (f) nominally preserves both $\hat{C}_{2}$ (on the ($1\bar{1}0$) surface termination) and time-reversal symmetries, while the configuration shown in the right panel of Figure \ref{f3} (f) preserves only the time-reversal symmetry (details can be found in Supplementary Note 2). In addition, both configurations of helical hinge modes in Figure \ref{f3} (f) that are consistent with our observations require broken inversion symmetry on the surface termination (which is a symmetry of the bulk of a bismuth crystal). However, this is not problematic since inversion is naturally broken for an island on a surface. Thus, we find that our observations corroborate, and are clearly indicative of, higher order topology  in bismuth.
\\

\noindent \textbf{Spectroscopic measurements on Bi$_{0.92}$Sb$_{0.08}$(110).}
We can now use our established signature of higher order topology in Bi to determine the topological nature of Bi alloyed with Sb. Previous theoretical studies have suggested that the Bi$_{1-x}$Sb$_{x}$ alloy undergoes a trivial-to-topological band-inversion transition at approximately $x= 0.04$ \cite{lenoir1996bi,fu2007topological,teo2008surface}. However, similar to Bi, the narrow band gaps and the existence of multiple bulk and surface bands near the Fermi energy have made the theoretical and experimental determination of the topological properties of Bi$_{1-x}$Sb$_{x}$ complicated. From our studies on Bi, we propose that the pattern of the edge modes on the (110) oriented rectangular islands may provide a clear picture of the topology in this compound. Figure \ref{f4r} (b) shows STM spectra at an edge and plateau of the (110) facet of Bi$_{0.92}$Sb$_{0.08}$. A peak centered around 100meV is observed at the edge, similar to that observed on Bi islands. Spectroscopic measurements across and along the edge shown in Figure \ref{f4r} (c) and Figure \ref{f4r} (d) respectively, reveal that this mode is localized to the edge. Importantly, similar to our data on Bi, we observe these edge modes on three out of the four edges of (110) oriented rectangular islands as shown in Figure \ref{f4r} (f,g,h,i). The explanation of this pattern of edge modes requires the system to be higher order, and hence requires that the bulk has trivial strong $\mathbbm{Z}_{2}$ topology. Hence, our data on Bi$_{0.92}$Sb$_{0.08}$ suggests that it has trivial strong $\mathbbm{Z}_{2}$ topology and falls in the HOTI regime. This bears further study using other experimental techniques since it calls into question some long-accepted predictions about BiSb alloys, and indeed may have implications for the celebrated identification of nearby Bi$_{0.9}$Sb$_{0.1}$ as a strong $\mathbbm{Z}_{2}$ topological insulator.

\section*{Discussion}
\quad \ In summary, we have studied (110) oriented Bi and Bi$_{0.92}$Sb$_{0.08}$ films with STM and STS. We observe bound modes on three out of four edges of a rectangular island on (110) facets. Moreover, the edge showing no edge modes is always perpendicular to the ($1\bar{1}0$) direction. Theoretically, the novel pattern of the metallic edge modes observed on the (110) facets can be explained by a higher order topology protected by time reversal ($\hat{T}$) symmetry and bulk two-fold rotation ($\hat{C}_{2}$) symmetry around the ($1\bar{1}0$)-axis \cite{hsu2019topology}. Based on this we confirm that Bi and Bi$_{0.92}$Sb$_{0.08}$ belong to the HOTI class of topological materials. The latter result is particularly interesting since previous theoretical studies indicate that Bi$_{1-x}$Sb$_{x}$ with $x \sim 0.08$ should already be in a strong $\mathbbm{Z}_{2}$ topological insulator phase heralded by a band inversion predicted at $x\sim 0.04$ \cite{lenoir1996bi,fu2007topological,teo2008surface}. Our work brings new insights to the topology of Bi$_{0.92}$Sb$_{0.08}$ which is revealed as a new HOTI.

\section*{Methods}

\textbf{Sample Preparation}: Highly crystalline films of Bi and  Bi$_{0.92}$Sb$_{0.08}$ were grown on n-doped silicon (111) wafers using a custom ultrahigh vacuum molecular beam epitaxy system. Before deposition, Si substrates were cleaned by ultrasonication in acetone and isopropanol. The substrates were then transferred to the MBE system (pressure $<2\times10^{-9}$ Torr) and degassed at $520^{\circ}C$ for 10-12 hrs. After degassing, the substrate was cleaned by flash annealing where they were repeatedly heated from $650^{\circ}C$ to $900^{\circ}C$ in 80 s. Each cycle helps to remove oxide layers or impurities from the silicon surface. The  surface quality was confirmed by the observation of clear $7\times7$ reconstructed RHEED patterns. The substrate was rapidly cooled from $650^{\circ}C$ to room temperature at the rate of $1-2^{\circ}C/s$ after 16 cycles of flash annealing. For film growth, bismuth (purity $\sim 99.9999 \% $) and antimony (purity $\sim 99.9999 \% $) sources were evaporated from standard Knudcell cells at flux rates 0.0769 $\mathring{A}/s$ and 0.0071 $\mathring{A}/s$ respectively as measured by a Quartz Crystal monitor. The substrates were held at room temperature during growth.
\\

\noindent\textbf{STM measurements}: Samples were transferred from MBE chamber to a vacuum suitcase at pressure $<3\times10^{-10}$ and subsequently transferred to a Unisoku UHV-STM without exposure to air. All the measurements were done at 4.3 K. Normal tungsten tips were prepared by etching method and used after cleaning by electron beam heating. All the topographic images were taken in constant-current mode and spectroscopic signals were recorded by standard lock-in amplifier with ac modulation 3-3.5 mV at 983 Hz frequency.

\vspace{24pt}

\section*{Supplementary note 1: 1D hinge modes at type-A edges in Bi(111) film}

\begin{figure}[h]
\centering
\includegraphics[width=0.75\columnwidth]{Bi111.pdf}
\caption{\textbf{STM studies on Bi(111) films grown on  silicon (111)substrates:} (a) A large area (700 nm $\times$700 nm) atomic force microscopy (AFM) image  showing triangular islands which confirms the (111) orientation of the bismuth film with the color bar showing the relative height.(b) 300 nm $\times$300 nm STM image of the film with the color bar showing the relative height. (c) Topography of a 28 nm $\times$28 nm triangular (111) facet with two kind of edges, marked as Type-A and Type-B. (d) dI/dV spectra on the two types of edges and near the center of the island.  The hinge mode density of states peak can be seen at the type-A edge. (e,f) STM dI/dV maps at 220 mV and 180 mV. The 180mV map clearly shows the localization of hinge modes at the three type-A edges.}
\label{fig:hinge_Bi}
\end{figure}

\section*{Supplementary note 2: Theoretical explanation for helical edge modes on the (110) facet}
\subsection*{Orientation of surface facets}

\label{app: orientation}
As shown in Fig. \ref{fig: lattice_BZ} (a), in the Cartesian coordinate system we choose, the lattice vectors can be expressed as
\begin{equation}
\label{eq: latticevec}
\vec{a}_{1}=(-\frac{1}{2} a,-\frac{\sqrt{3}}{6} a, \frac{1}{3} c),\quad \vec{a}_{2}=(\frac{1}{2} a,-\frac{\sqrt{3}}{6} a, \frac{1}{3} c),\quad \vec{a}_{3}=(0, \frac{\sqrt{3}}{3} a, \frac{1}{3} c),
\end{equation}
where $a$ and $c$ are lattice constants as illustrated in Fig. \ref{fig: lattice_BZ} (a). As shown in Fig. \ref{fig: lattice_BZ} (b), the three corresponding reciprocal-lattice vectors defined by $\bm{b}_{i} \cdot \bm{a}_{j}=2 \pi \delta_{i j}$ with $i,j=1,2,3$, are given by
\begin{equation}
\label{eq: rlatticevec}
\vec{b}_{1}=g(-1,-\frac{\sqrt{3}}{3}, \frac{a}{c}),
\quad \vec{b}_{2}=g(1,-\frac{\sqrt{3} }{3} , \frac{a}{c}),
\quad \vec{b}_{3}=g(0,\frac{2 \sqrt{3}}{3} , \frac{a}{c}),
\end{equation}
where $g=2\pi/a$. From Eq. \eqref{eq: rlatticevec}, we can derive that
\begin{equation}
\label{eq: vector}
\begin{array}{l}
(111)=\vec{b}_{1}+\vec{b}_{2}+\vec{b}_{3}=g(0,0 , 3a/c),
\\
(1\bar{1}0)=\vec{b}_{1}-\vec{b}_{2}=g(-2,0 , 0),
\\
(11\bar{2})=\vec{b}_{1}+\vec{b}_{2}-2\vec{b}_{3}=g(0,-2\sqrt{3},0),
\\
(110)=\vec{b}_{1}+\vec{b}_{2}=g(0,-2\sqrt{3}/3 , 2a/c),
\end{array}
\end{equation}
from which we see that $(1\bar{1}0)$, $(11\bar{2})$, and $(111)$ are the Cartesian $x$, $y$, and $z$ directions (up to a sign), respectively. We can also see that $(110)$, $(1\bar{1}0)$ and $\vec{a}_{3}$ are orthogonal to each other.

\begin{figure}[h]
\centering
\includegraphics[width=0.8\columnwidth]{lattice_BZ.pdf}
\caption{(a) Crystal structure of Bi and/or BiSb alloys with primitive cell, indicating rhombohedral vectors  ($\vec{a}_{1}$,$\vec{a}_{2}$,$\vec{a}_{3}$) and $(11\bar{2})$, $(111)$, $(1\bar{1}0)$, $(110)$ planes. (b) First Brillouin zone of Bi and/or BiSb alloys band structure.}
\label{fig: lattice_BZ}
\end{figure}

\subsection*{Possible configurations of hinge modes under different symmetries}
As discussed in the main text, the exotic three-sided hinge modes configuration observed on the rectangular island of a (110) surface facet is ascribed to higher order topology protected by the $\hat{C}_{2}$ symmetry around the  $(1\bar{1}0)$-axis, and the time-reversal symmetry. We claim that it is necessary to break inversion symmetry to observe this hinge modes configuration. To support this claim, we enumerate all possible configurations of edge states when both $\hat{C}_{2}$ and inversion symmetries are preserved, only $\hat{C}_{2}$ is preserved, and only inversion is preserved. We also show two typical configurations which break both inversion and $\hat{C}_{2}$ rotation symmetries, one of which also gives the three-sided configurations. The enumeration shown in Fig. \ref{fig:enumerate} clearly show that the three-sided edge modes configurations can only happen when the inversion symmetry is broken, irrespective of whether the $\hat{C}_{2}$ rotation symmetry is broken or not.
\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{enumerate.pdf}
\caption{All possible configurations of helical hinge modes on the (110) surface of $\rm{Bi}$ and/or $\rm{Bi_{0.92}Sb_{0.08}}$ when (a) both $\hat{C}_{2}$ and inversion symmetries are preserved, (b) inversion is preserved but $\hat{C}_{2}$ is broken, and (c) $\hat{C}_{2}$ is preserved but inversion is broken. (d) shows two typical configurations when both $\hat{C}_{2}$ and inversion symmetries are broken. The red-color hinges are hinges where the surface mass vanish, $i.e.$, where there are gapless hinge modes. Note that $\hat{C}_{2}$ here refers to the two-fold rotation around $(1\bar{1}0)$ axis.}
\label{fig:enumerate}
\end{figure}

\clearpage

\section*{Supplementary note 3: Raw Data on Bi(110) and Bi$_{0.92}$Sb$_{0.08}$(110) Films}
\begin{figure}[h]
\centering
\includegraphics[width=0.85\columnwidth]{Raw_Data.pdf}
\caption{(a) Topographic image of the of same Bi(110) rectangular island shown in Figure 3 (a) in the main manuscript. (b,c,d,e)  Raw spectroscopic data (without slope subtraction) (at ac modulation 3.5 mV and current 60 pA) along the four edges marked Edge A , Edge B, Edge C and Edge D (indicated by purple arrows in Figure S4(a)). (f) Topographic image of the same Bi$_{0.92}$Sb$_{0.08}$(110) rectangular island as shown in Figure 4 (e) in the main manuscript. (g,h,i,j) Raw spectroscopic data (without slope subtraction) (at ac modulation 3 mV and current 120 pA) along the four edges marked Edge A, Edge B, Edge C and Edge D ( shown as blue arrows along the edges in Figure S4(f)).}
\end{figure}

\section*{Supplementary note 4: Rutherford Backscattering Spectrometery (RBS) to confirm Sb composition in alloy}

\begin{figure}[h]
\centering
\includegraphics[width=0.50\columnwidth]{rbs.pdf}
\caption{RBS data (in solid black color) on BiSb alloy to confirm the composition of Sb with nicely fitted simulated data (in dashed red color line). Fitting confirms Sb composition 0.08$\%$. }
\label{fig:rbs}
\end{figure}

\newpage

\section*{Supplementary note 5: Additional Raw Data on Bi(110) and Bi$_{0.92}$Sb$_{0.08}$(110) Films}

\begin{figure}[h]
\centering
\includegraphics[width=0.85\columnwidth]{additional_data.pdf}
\caption{(a)Topographic image of the of same Bi(110) rectangular island. (b,c,d,e)  Raw spectroscopic data (without slope subtraction) (at ac modulation 3.5 mV and current 80 pA) along the four edges marked Edge A , Edge B, Edge C and Edge D (indicated by purple arrows in Figure S6(a)). (f) Topographic image of the same Bi$_{0.92}$Sb$_{0.08}$(110) rectangular island. (g,h,i,j) Raw spectroscopic data (without slope subtraction) (at ac modulation 2.5 mV and current 120 pA) along the four edges marked Edge A, Edge B, Edge C and Edge D ( shown as blue arrows along the edges in Figure S6(f)). The red block arrows indicate the sharp peaks corresponding to the edge modes that are only seen on three of the four edges.}

\end{figure}

\newpage
\section*{Supplementary note 6: Schematic of Bi(110) lattice structure}

The schematic of Bi(110) lattice structure, as shown in Figure S7(a), aligned correctly with the long edge  and short edge of (110) island. The atoms along the two long edges  of the (110) island bond differently to the layer below. However, in our data, both long edges of the (110) islands always show edge states. Put another way, the two short edge sides are identical from a bonding perspective. But, for a given island, one short side shows hinge states while the other does not. Since the appearance hinge modes at three out of four edges has no connection with the geometry of the bonds at the edges, we can rule out bonding effects as an explanation for the pattern of our edge modes. Our data are instead well explained by the Higher order topological model.

\begin{figure}[h]
\centering
\includegraphics[width=0.65\columnwidth]{dagling_bond2.pdf}
\caption{(a) Schematic of the top view of first two layers of  Bi(110) lattice structure, indicating two long edges of island bonded differently. Red atoms shows the first layer of two-dimensional rectangular lattice with lattice constants. The red lines and black lines between atoms represent long and short bonds respectively. The mirror plane (`m') of the structure are also shown as dashed yellow line.The lattice constants along long edge of island is 0.454 nm and short edge of island is 0.472 nm.  (b, c) Side view of the lattice structure perpendicular and parallel to mirror plane (`m') respectively.(d) STM image of Bi(110) island shows the long and short edge sides of island same as in (a). (e) Zoomed image of Bi(110) island shown in (d), showing pseudo-cubic arrangements of atoms. (f,g) Lattice constants along the long edge side (corresponding red line in (e)) of island and short edge side (corresponding blue line in (e)) of the island respectively.}
\end{figure}

\end{document}


\section{Introduction}
Electromagnetic Transient (EMT) analysis plays important roles in modern power systems \cite{kundur2004definition,dommel1997techniques}.
EMT studies consider three-phase voltage and current during simulations and the time-step is in microseconds or even less to simulate power electronics switching. With the development of modern power grids, power electronics are highly penetrated in both generation system and transmission system \cite{tolbert2005power}.
EMT simulators model power system with details introducing long simulation time.

Many attempts have been made to improve the computation efficiency \cite{anderson1995new,sultan1998combined,van2014advanced,gao2009frequency,yang2015dynamic,plumier2014relaxation,jalili2009interfacing,zhang2012development,heffernan1981computation,su2004recent}. But it’s challenging in the fact of that the conventional numerical method is adopted to solve the differential algebraic equations representing power system electromagnetic transient by time domain simulation which consumes unaffordable time for large-scale systems. Physically, power system is a very large-scale circuit consisting resistors, inductors, reactors, and other circuit elements. Using integrated circuit chip to mimic power system with high fidelity in real time could leverage the advanced technology of integrated circuit industry in terms of modelling, simulation, and manufacturing. In work \cite{zhang2018vlsi}, a very large-scale integration (VLSI) architecture for exciters, governors, and power system stabilizers is presented and a fixed time-step approach is used for power system electromechanical transient simulation purpose. Circuit simulation results demonstrate that power system modeled by VLSI circuits can be simulated faster than real time without losing accuracy. It’s worth noting that the power grid can be essentially simplified as a linear network composed by resistors, inductors and capacitors when investigating the node voltages with the node current injections. The authors of \cite{dai2019simplify} take the advantage of the linearity and propose a method to simplify power flow calculation by utilizing terminal circuit model and PMU measurements. This approach is further extended in work \cite{zhang2021dynamic}, which proposed a new method to model power systems as an integrated circuit and simulate the system dynamics by integrated circuit simulator. In this paper, the terminal circuit method is further investigated and applied to transient simulation. The research goal of this paper is to extend the terminal circuit based steady-state and dynamic cases with appropriate power system equipment electromagnetic models, to study electromagnetic transients under normal or abnormal (i.e., faulty) system conditions.

In this paper, an approach to model power system components is introduced which is capable of simulating electromagnetic transients using terminal circuit. The construction of a 5 bus test case is provided to verify the proposed approach.
Adaptive time-stepping method is used to improve computation efficiency of EMT simulation. Balanced three-phase fault is applied to the system for further analysis and the simulation results are compared with commercial software to
verify that the test case built with terminal circuit is able to simulate electromagnetic transients accurately and efficiently.

The remainder of the paper is organized as follows: In Section \uppercase\expandafter{\romannumeral2}, terminal circuits are introduced for power system electromagnetic modeling.
Section \uppercase\expandafter{\romannumeral3} provides details of modeling a 5-bus system, and Section \uppercase\expandafter{\romannumeral4} presents the conclusions.

\section{Electromagnetic Modelling}
Typically, power system is consisted of generators, trans- mission lines, transformers, and series RLC component, non- ideal voltage source in special cases. They can be further represented by basic circuit elements of resistors, reactors, capacitors, voltage source, current source, and their combination in very large-scale integrated circuit. In this section, power system  component electromagnetic models are constructed by using terminal circuit in VLSI circuit simulator.
\subsection{Park’s Transformation }
In EMT studies, generator is modelled by an equivalent current source on the fictitious direct (d) and quadrature (q) axes \cite{kundur1994power} and interfacing with the three-phase network model by Park’s transformation. The Park’s transformation and its reverse transformation are given by (\ref{Park}) (\ref{equ:park_1}), and (\ref{equ:park_2}). And the transform block of Park’s Transformation is shown in Fig. \ref{fig:K_transform}.

\begin{figure}[thpb]
\centerline{\includegraphics[scale=0.45]{K_transform.PNG}}
\caption{Park's Transformation}
\label{fig:K_transform}
\end{figure}

\vspace{-2.5mm}
\begin{align} \label{Park}
K &= \sqrt{\frac{2}{3}}
\begin{bmatrix}
\cos(\theta) & \cos(\theta - \frac{2\pi}{3}) & \cos(\theta + \frac{2\pi}{3})\\
\sin(\theta) & \sin(\theta - \frac{2\pi}{3}) & \sin(\theta + \frac{2\pi}{3})\\
\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}}
\end{bmatrix}
\end{align}
\vspace{-1.5mm}
\begin{equation} \label{equ:park_1}
V_{ri} = K\times V_{abc}
\end{equation}
\vspace{-5.5mm}
\begin{equation} \label{equ:park_2}
I_{abc} = K^{-1}\times I_{dq} \textcolor{white}{11.1}
\end{equation}

\subsection{Power System Components Modeling}
In power system, the $\pi$-type transmission line is a typical three phase RLC circuit. The single phase transmission line representation is shown in Fig. \ref{fig:pi_line}.
\vspace{-3.5mm}
\begin{figure}[h!]
\centerline{\includegraphics[scale=0.35]{Transmission_line_pi.png}}
\caption{Transmission Line RLC Model}
\label{fig:pi_line}
\end{figure}

For Transformer, Fig. \ref{fig:Transformer_ratio} (a) shows the equivalent circuit for a single-phase two-winding transformer and the turns ratio of winding 1 (\textit{$T_{k1}$}) and winding 2 (\textit{$T_{k2}$}) is given by (\ref{turns ratio}).
\begin{figure}[h]
\centerline{\includegraphics[scale=0.60]{Single_phaseTransformer_new.png}}
\caption{Two-winding Transformer}
\label{fig:Transformer_ratio}
\end{figure}
\begin{equation} \label{turns ratio}
T_{k1} (T_{k2}) = \frac{V_{Rated}}{V_{Base}}
\end{equation}
where $V_{Rated}$ is the rated line voltage of winding 1 (2) and $V_{Base}$ is the voltage base of winding 1 (2). In this work, the authors utilize two sets of voltage controlled voltage source to implement the turns ratio  function and two sets of current controlled current source to enforce the current through the mutual impedance, as shown in Fig. \ref{fig:Transformer_ratio} (b). Similarly, the equivalent circuit for a  single-phase three-winding transformer is depicted in Fig. \ref{fig:Transformer_3_winding} (a) and three sets of voltage (current) controlled voltage (current) source are used in this work to implement the voltage and current relationship, as in Fig. \ref{fig:Transformer_3_winding} (b). In addition, a three-phase delta-wye transformer connection is constructed using three sets of single-phase transformers as shown in Fig. \ref{fig:Delta-wye connection}. In this three-phase formulation, the left-hand-side terminals of the transformer show a delta connection with a positive sequence of a-b-c phase while the right-hand-side terminals have a wye connection with all return terminals \textit{$v_{2ao}$}, \textit{$v_{2bo}$}, and \textit{$v_{2co}$} attached to the ground.

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.60]{Three_windingTransformer_new.png}}
\caption{Three-Winding Transformer}
\label{fig:Transformer_3_winding}
\end{figure}

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.54]{Trans_3phase.PNG}}
\caption{Delta-wye Transformer Connection Diagram}
\label{fig:Delta-wye connection}
\end{figure}

In typical power system electromagnetic simulation, three-phase series RLC and the non-ideal voltage source are commonly used, their models are shown in  Fig.\ref{fig:RLC} (a) and (b).
\begin{figure}[h!]
\centerline{\includegraphics[scale=0.5]{V_source_and_Series_RLC.png}}
\caption{(a) Series RLC (b) Non-Ideal Voltage Source}
\label{fig:RLC}
\end{figure}

In electromagnetic transient simulation, 6$^{th}$ order generator model is applied. The associated governor model and exciter model are illustrated as follow. The example transfer function diagrams of a governor and an exciter are shown in  Fig. \ref{fig:IEEEG3} and Fig. \ref{fig:DC1A}, respectively.
\begin{figure}[h!]
\centering
\includegraphics[scale=0.458]{IEEEG3.png}
\caption{IEEEG3 Governor}
\label{fig:IEEEG3}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[scale=0.484]{DC1A.PNG}
\caption{DC1A type exciter}
\label{fig:DC1A}
\end{figure}

The transfer function diagrams are modeled by integrating basic control blocks, such as low pass filter, high pass filter, gain block, adder, integrator etc. Those basic control blocks are modeled by circuit using integrated circuit simulator as shown in Fig. \ref{fig:Control Blocks}.
\vspace{-1.5mm}
\begin{figure}[h!]
\centering
\includegraphics[scale=0.595]{ControlBlocks.PNG}
\caption{Basic Control Block Circuit Model}
\label{fig:Control Blocks}
\end{figure}

While dependent voltage and current source, as shown in Fig. \ref{fig:Control Blocks} (c) and Fig. \ref{fig:Control Blocks} (d), are also known as B-source devices and commonly used as analog behavioral modeling (ABM) in integrated circuit \cite{Xyce}. They are adopted to accommodate complicated nonlinear relations between the input and output signal such as in the integrator model shown in Fig.  \ref{fig:Control Blocks} (g). The integrator circuit is a nonlinear dependent current source in parallel with an ideal capacitor.

\begin{equation}
\frac{dV_0}{dt} = \frac{1}{C}f(V_i)
\end{equation}
where, the nonlinear dependent current source provides a time dependent current $f(V_i)$ that simulates the time derivative of a state variable labeled as $V_o$. The capacitance of the capacitor is set as $\tau$ to make the state variable $V_o$ representing the desired integration.

The lead-lag filter in the IEEEG3 governor and the DC1A type exciter can be modeled by adding the signals of a low pass filter and a high pass filter. And the complicated saturation function $S_E(E_{FD})$ can be modeled by a nonlinear dependent voltage source.

\section{Case Study}
To verify the effectiveness of the proposed method, a simple 5-bus system is modeled by SPICE-level analog circuit and simulated by using an integrated circuit simulator, Xyce \cite{Xyce}.

As shown in Fig. \ref{fig:5buscase}, the 5-bus system contains 1 synchronous generator, 1 non-ideal voltage source, 1 series RLC, 3 $\pi$-type transmission lines and 1 circuit breaker. The fundamental frequency of the system is 50-Hz.
\vspace{-1.5mm}
\begin{figure}[thpb]
\centering
\includegraphics[scale=0.49]{5_bus_case.PNG}
\caption{5 bus case}
\label{fig:5buscase}
\end{figure}
\vspace{-4.5mm}
\begin{figure}[h]
\centering
\includegraphics[scale=0.59]{TimeStepSize_new.png}
\caption{Adaptive time-step size according to different activities: small time-step during the initialization $\&$ fault period}
\label{fig:TimeStep}
\end{figure}

To illustrate the electromagnetic transient behavior of the studied system, a balanced three-phase fault is applied at bus 4 at 10.0s and is cleared at \textcolor{black}{10.06s}. In the transient simulation, adaptive time step is applied with a small initial time-step value as 20us and the maximal time-step as 2ms to achieve better simulation efficiency. When the simulation starts, the step size is adapted to increase till the fault happened as shown in   Fig.\ref{fig:TimeStep} (a) and (b). Upon the time of the fault start and clear, the step sizes are reduced to keep acceptable truncation errors during the large transient as shown in  Fig.\ref{fig:TimeStep} (c). After the transient is settled down within seconds after a fault is cleared, as shown in Fig.\ref{fig:TimeStep} (d), the step size increases to the maximal time-step to achieve high computation efficiency. The adaptive time-step balances the computation accuracy and efficiency. By using the adaptive time-step simulation, it takes 2.66s to simulate 20 seconds transients. In contrast, by using fixed time-step at 1ms, it takes 3.98s to complete the 20 seconds transient simulation.

\begin{figure*}[ht!]
\centering
\includegraphics[scale=0.680]{5BusFaultWithZoom_ver3.PNG}
\caption{Simulation results from the 5 bus case. The fault is applied at 10.00 s and cleared 10.06 s.}
\label{fig:5busVmag}
\end{figure*}

\begin{figure*}[bh!]
\centering
\includegraphics[scale=0.680]{CompareWIthZoom_ver2.PNG}
\caption{Positive sequence voltages of 5 buses collected from the developed case using the Xyce circuit simulator and ADPSS}
\label{fig:5busVmag_compare}
\end{figure*}

To verify the computation accuracy, the simulation results provided in Fig. \ref{fig:5busVmag} are compared with the results simulated by a commercial electromagnetic simulation software Advanced Digital Power System Simulator (ADPSS) \cite{zhongxi1998power} developed by China Electric Power Research Institute (CEPRI) as shown in Fig. \ref{fig:5busVmag_compare}.

As shown in Fig. \ref{fig:5busVmag}, after the fault is applied at 10.0s, the three phase voltages decrease while  the short-circuit current at the faulted bus is increased. Upon the time of the fault removal after three cycles, the voltage spikes caused by phase jump at the clearing time are observed. The fast transients are eliminated within 1 cycle and then rapidly settled down to an equilibrium point. Note that the post-fault three-phase voltage amplitudes are approaching to their pre-fault values.

In Fig. \ref{fig:5busVmag_compare}, the positive sequence voltages at each bus simulated by the proposed method are compared with those collected from ADPSS which demonstrates that the simulation results by the two simulators are consistent.

\section{Conclusion}
In this paper, a new approach is proposed to model power system electromagnetic transients by integrated circuit using SPICE-level analog circuit language. Preliminary studies are carried out to give insights on modeling power system components using terminal circuits.  A Park’s transformation circuit model is developed as an interface between the three-phase network reference frame and the single-phase dq reference frame. A 5-bus case is used as an example to model power system electromagnetic transients by the proposed method. A balanced three-phase fault is applied for the case study. To improve the computation efficiency, adaptive time-step method is used in the simulation. The accuracy of the simulation results are validated by comparing with commercial software.

This paper proved the concept of using integrated circuit to model power system for the purpose of transient simulation. As enabling technology, the approach and the simulation results verified the possibility and feasibility of using a large-scale integrated circuit chip to model power system and eventually replace the time consuming traditional numerical methods for transient simulation. As demonstrated in the initial practice, the method has a great potential to improve the large-scale power system transient study accuracy and efficiency.

\vspace{12pt}

{
\small