Probing Chiral Kitaev Spin Liquids via Dangling Boundary Fermions

In this work, we are only interested in low-energy dynamics below $\Delta_{\rm v}$, which corresponds to the zero-flux sector of the Kitaev model. In this regime, the spin dynamics is completely determined by the spin excitations generated by the “dangling spins” on the boundary that create “dangling fermions”. For that purpose, we first introduce a relevant effective Hamiltonian governing the low-energy dynamics.

Projecting $\hat{\mathscr{H}}$ into the zero-flux sector yields two distinct contributions. The first one originates from the Kitaev Hamiltonian $\hat{\mathscr{H}}_{K}$ and has a quadratic form in terms of the $c$ fermions:

$$\hat{\tilde{\mathscr{H}}}_{K}=i\sum_{\langle jl\rangle}\hat{c}_{j\in A}\hat{c}_{l\in B},$$

(6)

where the sum runs over all non-equivalent nearest neighbor bonds, $A$ and $B$ refer to the two sublattices of the honeycomb lattice, and we are adopting the standard gauge choice for the zero-flux sector [1]. The second contribution arises from the Zeeman coupling to an external field, $\hat{\mathscr{H}}_{Z}$. Projected into the zero-flux sector, this term generates a boundary contribution,

$$\hat{\tilde{\mathscr{H}}}_{Z}=i\sum_{j\in{\rm boundary}}h^{\mu_{j}}\hat{b}_{j}^{\mu_{j}}\hat{c}_{j},$$

(7)

which is linear in the magnetic field and hybridizes each dangling fermion $\hat{b}_{j}^{\mu_{j}}$ (where $\mu_{j}$ is $x$, $y$, or $z$, depending on $j$) with the matter fermion $\hat{c}_{j}$ via the $\mu_{j}$-th component of the Zeeman field. The two terms described in this paragraph are illustrated by the pink solid arrows in Fig. 2(a).

Higher-order processes generate further symmetry-allowed terms, as indicated by blue dashed arrows in Fig. 2(a). Given the dominance of the Kitaev interactions, we restrict our attention to terms arising from leading-order perturbations. Additionally, we assume that these interactions do not alter the quasiparticle picture, retaining only quadratic terms in the effective Hamiltonian. In the bulk, higher-order processes generate the effective three-spin interaction that immediately gaps out any gapless bulk fermions [1]:

$$\hat{\tilde{\mathscr{H}}}_{1}^{\prime}=-\kappa\sum_{\langle jkl\rangle_{\alpha\beta\gamma}}\sigma_{j}^{\alpha}\sigma_{k}^{\beta}\sigma_{l}^{\gamma}=i\kappa\sum_{\langle jkl\rangle_{\alpha\beta\gamma}}u_{kj\parallel\alpha}u_{kl\parallel\gamma}c_{j}c_{l},$$

(8)

where $(\alpha\beta\gamma)$ is a permutation of $(xyz)$ in each term, and $\langle jkl\rangle_{\alpha\beta\gamma}$ denotes a three-site path formed by two consecutive nearest-neighbor bonds $\langle jl\rangle\parallel\alpha$ and $\langle lk\rangle\parallel\beta$. This effective interaction preserves the integrability of the model. The scaling of $\kappa$ with respect to ${\bm{h}}$ is given by $c_{1}h+c_{2}h^{3}+\dots$ [58], where $c_{1}\sim 4\Gamma^{\prime}/\Delta_{\rm v}$ originates from a second-order process involving the $\Gamma^{\prime}$ term the and Zeeman coupling term, while $c_{2}$ arises from a third-order contribution in ${\bm{h}}$ [1]. The factor of $4$ in $c_{1}$ accounts for two different processes, illustrated by the first two lines of Fig. 2(b), which contribute to the same matrix element. On the boundary, higher-order processes endow the dangling fermions with a finite hopping term:

$$\hat{\tilde{\mathscr{H}}}_{2}^{\prime}=\sum_{jl\in{\rm boundary}}t_{jl}b_{j}^{\mu_{j}}b_{l}^{\mu_{l}}.$$

(9)

This hopping arises, for example, from second-order perturbation processes involving the $\Gamma$ term and the Zeeman coupling term, as shown in Fig. 2(c), which give a hopping amplitude $t_{jl}\sim 4h\Gamma/\Delta_{\rm v}$. In the following, we consider only the nearest-neighbor hopping term, denoted as $t_{1}$. Hopping terms beyond nearest neighbors merely affect the dispersion of the boundary mode without introducing any qualitatively new effects.

To summarize, the low-energy physics is described by the effective Hamiltonian

$$\hat{\tilde{\mathscr{H}}}=\hat{\tilde{\mathscr{H}}}_{K}+\hat{\tilde{\mathscr{H}}}_{Z}+\hat{\tilde{\mathscr{H}}}_{1}^{\prime}+\hat{\tilde{\mathscr{H}}}_{2}^{\prime}+...,$$

(10)

which takes a quadratic form,

$$\hat{\tilde{\mathscr{H}}}=i\sum_{mn}A_{mn}\hat{F}_{m}\hat{F}_{n},$$

(11)

in the fermionic representation. Here $A_{mn}=-A_{nm}\in{\cal R}$, and $\hat{F}_{m}$ refers to the components of the Majorana-fermion vector,

$\displaystyle\hat{\bm{F}}=(\hat{c}_{1},...,\hat{c}_{N},\hat{b}_{1}^{\mu_{1}},...,\hat{b}_{M}^{\mu_{M}})^{T},$

(12)

with $N$ being the number of lattice sites, $M$ the number of dangling fermions on the boundary, and $\mu_{j}$ the spin component corresponding to a given dangling fermion (see, e.g., Fig. 2).

The effective Hamiltonian $\hat{\tilde{\mathscr{H}}}$ can be diagonalized through a unitary transformation,

$\displaystyle\hat{F}_{m}=\sqrt{2}\sum_{s=1}^{\cal N}\Phi_{ms}\hat{f}_{s}+{\rm H.c.},$

(13)

where ${\cal N}=N+M$, and ${\bf\Phi}$ is an ${\cal N}\times{\cal N}$ unitary matrix. The $s$-th column of this matrix, ${\bm{\Phi}}_{s}$, satisfies the eigenvalue equation

$\displaystyle i{\bm{A}}{\bm{\Phi}}_{s}=\epsilon_{s}{\bm{\Phi}}_{s},\ \ i{\bm{A}}{\bm{\Phi}}_{s}^{*}=-\epsilon_{s}{\bm{\Phi}}_{s}^{*},$

(14)

where ${\bm{A}}$ is the matrix of the quadratic Hamiltonian $\hat{\tilde{\mathscr{H}}}$ in Eq. (11). The eigenvalues appear in pairs $\pm\epsilon_{s}$ due to particle-hole symmetry. After applying the transformation in Eq. (13), the effective Hamiltonian takes a diagonal form

$$\hat{\tilde{\mathscr{H}}}=\sum_{\epsilon_{s}\geq 0}\epsilon_{s}\hat{f}_{s}^{\dagger}\hat{f}_{s}.$$

(15)

To distinguish the $b$ and $c$ fermions, we introduce the notations $\hat{\bm{c}}\equiv(\hat{c}_{1},...,\hat{c}_{N})^{T}$ and $\hat{\bm{b}}\equiv(\hat{b}_{1}^{\mu_{1}},...,\hat{b}_{M}^{\mu_{M}})^{T}$. Moreover, we define ${\bm{\Phi}}_{s}^{T}\equiv({\bm{\varphi}}_{s}^{T},{\bm{\phi}}_{s}^{T})$, where ${\bm{\varphi}}_{s}$ is a $N\times 1$ vector and ${\bm{\phi}}_{s}$ a $M\times 1$ vector. The transformation in Eq. (13) leads to two equations

	$\displaystyle\hat{c}_{j}$	$\displaystyle=$	$\displaystyle\sqrt{2}\sum_{s}\varphi_{js}\hat{f}_{s}+{\rm H.c.},$
	$\displaystyle\hat{b}_{j}^{\mu_{j}}$	$\displaystyle=$	$\displaystyle\sqrt{2}\sum_{s}\phi_{js}\hat{f}_{s}+{\rm H.c.},$		(16)

where $\varphi_{js}$ and $\phi_{js}$ are the $j$-th components of ${\bm{\varphi}}_{s}$ and ${\bm{\phi}}_{s}$, respectively. The inverse transformation is given by

$\displaystyle\hat{f}_{s}^{\dagger}={1\over\sqrt{2}}\left(\sum_{j=1}^{N}\varphi_{js}\hat{c}_{j}+\sum_{j=1}^{M}\phi_{js}\hat{b}_{j}^{\mu_{j}}\right).$

(17)

Our primary focus is the local dynamic spin structure factor (DSSF) on boundary lattice sites $j$,

$\displaystyle{\cal S}_{j}^{\mu_{j}\mu_{j}}(\omega)={1\over 2\pi}\int_{-\infty}^{\infty}dte^{i\omega t}\langle\sigma_{j}^{\mu_{j}}(t)\sigma_{j}^{\mu_{j}}(0)\rangle,$

(18)

which can be measured by STM. We specifically consider the spin components $\mu_{j}$ corresponding to the dangling fermions at the boundary sites $j$, as they excite low-energy boundary fermion modes without generating $Z_{2}$ gauge fluxes with a finite energy gap $\sim\Delta_{\rm v}$.

In the Majorana fermion representation, the local DSSF reduces to a convolution of two Majorana fermion spectral functions ${\cal A}_{F_{m}F_{m}}(\omega)$:

	$\displaystyle S_{j}^{\mu_{j}\mu_{j}}(\omega)$	$\displaystyle=$	$\displaystyle\int_{0}^{\omega}d{\nu}\left[{\cal A}_{b_{j},b_{j}}(\nu){\cal A}_{c_{j},c_{j}}(\omega-\nu)\right.$		(19)
			$\displaystyle\left.-{\cal A}_{b_{j},c_{j}}(\nu){\cal A}_{c_{j},b_{j}}(\omega-\nu)\right],$

where we have abbreviated the flavor index of the dangling fermions for notation convenience, and

$\displaystyle{\cal A}_{F_{m}F_{n}}(\omega)={1\over 2\pi}\int_{-\infty}^{\infty}dte^{i\omega t}\langle\hat{F}_{m}(t)\hat{F}_{n}(0)\rangle.$

(20)

Since $\hat{F}^{2}_{m}=1$, the spectral function satisfies the sum rule

$$\int_{-\infty}^{\infty}d\omega{\cal A}_{F_{m}F_{m}}(\omega)=1,$$

(21)

while the off-diagonal elements satisfy

$$\int_{-\infty}^{\infty}d\omega{\cal A}_{F_{m}F_{n}}(\omega)=-\int d\omega{\cal A}_{F_{n}F_{m}}(\omega)$$

(22)

due to the anti-commutation relation of the Majorana fermions: $\{\hat{F}_{m},\hat{F}_{n}\}=2\delta_{mn}$. Furthermore, since $\sigma_{j}^{\mu}\sigma_{j}^{\mu}=\sigma_{j}^{0}$, the DSSF satisfies the sum rule

$\displaystyle\int_{-\infty}^{\infty}d\omega S_{j}^{\mu_{j}\mu_{j}}(\omega)=1.$

(23)

The Majorana fermion spectral functions can be computed by expressing the $\hat{F}_{m}$ operators in terms of the quasi-particle operators $\hat{f}_{s}$ that diagonalize $\tilde{\mathscr{H}}$:

	$\displaystyle{\cal A}_{F_{m}F_{n}}(\omega)$	$\displaystyle=$	$\displaystyle\sum_{\epsilon_{s}\geq 0}\left[{\cal M}^{*}_{F_{m}}(s){\cal M}_{F_{n}}(s)\delta(\omega-\epsilon_{s})\right.$		(24)
		$\displaystyle+$	$\displaystyle\left.{\cal M}^{*}_{F_{n}}(s){\cal M}_{F_{m}}(s)\delta(\omega+\epsilon_{s})\right],$

where

$\displaystyle{\cal M}_{F_{m}}(s)=\langle s\rvert\hat{F}_{m}\rvert 0\rangle=\sqrt{2}\Phi_{ms}^{*},$

(25)

and $\rvert s\rangle=\hat{f}_{s}^{\dagger}\rvert 0\rangle$ represents an excited single-quasiparticle state.

We begin by generalizing the analysis in Ref. [53] to arbitrary Kitaev spin exchange. Since we are adopting PBC along the boundary direction, the eigenmodes have a well defined quasi-momentum $k=2\pi n/L$, where $n$ is an integer number that belongs to the interval $-L/2<n\leq L/2$ ($L$ is assumed to be even). For a generic applied magnetic field with $\kappa\neq 0$, the energy spectrum consists of gapped bulk modes, where a gapless Dirac cone at $k=k_{0}$ acquires a finite mass $\propto\kappa$.

In the absence of the hybridization term $\hat{\tilde{\mathscr{H}}}_{Z}$, the boundary modes consist of two branches: one branch features a flat zero-energy spectrum representing the dangling fermions, while the other branch corresponds to a chiral mode formed by the matter fermions. This chiral mode only exists for $\rvert k\rvert>k_{0}$ and connects the negative-energy Dirac cone at $k=-k_{0}$ with the positive-energy Dirac cone at $k=k_{0}$.

The Zeeman term strongly hybridizes the dangling fermion modes at $\rvert k\rvert>k_{0}$ with the chiral mode, while those at $\rvert k\rvert<k_{0}$ retain their dominant dangling-fermion character because they can only hybridize with bulk matter-fermion modes. Fig. 3(d) shows the boundary modes on one edge before (black dashed line) and after (solid colored line) hybridization for the zigzag boundary depicted in Fig. 3(a).

Fig. 3(a,d,g) show the results obtained from numerical diagonalization of the effective Hamiltonian in Eq. (11). While these calculations are performed for a finite value of the magnetic field, the general structure can be understood by analyzing the small-field limit, $\rvert h_{\mu}\rvert\ll\Delta_{\rm v}$, where the spectrum of the boundary modes can be found analytically. For simplicity, we also consider the limit of $\kappa=0$. An analytical calculation reveals that the crossover momentum $k_{0}$ is then given by

$\displaystyle k_{0}=\cos^{-1}\left(\frac{K_{z}^{2}-K_{x}^{2}-K_{y}^{2}}{2K_{x}K_{y}}\right),$

(26)

which recovers $k_{0}=2\pi/3$ in the isotropic limit. The condition $0<k_{0}<\pi$ holds within the entire non-Abelian phase.

For $\rvert k\rvert<k_{0}$, the boundary-mode energy is zero to linear order in $h_{\mu}$ and has dominant dangling-fermion character, as captured by the wave function

$\displaystyle\varphi_{j,k}\simeq 0,\qquad\phi_{j,k}\simeq\frac{1}{\sqrt{L}}e^{ikj}.$

(27)

Note that Majorana fermions are real operators, implying that $b_{-k}^{\dagger}=b_{k}$ describe the same degree of freedom.

For $\rvert k\rvert>k_{0}$, the dangling and matter fermions form a hybridized mode with energy linear in $h_{\mu}$. For the zigzag boundary orientation shown in Fig. 3(a), this hybridization is linear in the $z$-component $h_{z}$ of the external field [53]:

$\displaystyle E_{k}=2h_{z}\sqrt{1-\lambda_{k}^{2}},$

(28)

where

$\displaystyle\lambda_{k}^{2}=\frac{K_{x}^{2}+K_{y}^{2}+2K_{x}K_{y}\cos k}{K_{z}^{2}}$

(29)

is a dimensionless parameter defined by $\lambda_{k}=\exp(-9/4\xi_{k})$, with $\xi_{k}$ being the the localization length of the boundary mode for wave vector $k$. Clearly, the boundary mode only exists for $\lambda_{k}^{2}<1$, and the crossover momentum $k=k_{0}$ given in Eq. (26) precisely corresponds to $\lambda_{k}^{2}=1$.

The minimum of $E_{k}$ is at $k=k_{0}$ where $E_{k}=0$. Around this point, the dispersion reads

$\displaystyle E_{k}\simeq 2h_{z}\left[\left(\frac{K_{x}-K_{y}}{K_{z}}\right)^{2}+1\right]^{1/4}\left[\left(\frac{K_{x}+K_{y}}{K_{z}}\right)^{2}-1\right]^{1/4}\sqrt{\rvert k\rvert-k_{0}}.$

The singular square-root behavior of the dispersion leads to the following density of states (DOS) for this hybridized mode,

$\displaystyle\rho(\omega)$

$\displaystyle=$

$\displaystyle\frac{1}{\pi}\int_{k_{0}}^{\pi}dk\delta(\omega-E_{k})=\frac{1}{\pi}\frac{1}{\rvert{dE_{k}/dk}\rvert}\rvert_{E_{k}=\omega},$

(31)

which is linear in $\omega$ for small energies: $\rho(\omega)\propto\omega$. The maximum of $E_{k}$ is at $k=\pi$ where $E_{\pi}=2h_{z}\sqrt{1-(K_{x}-K_{y})^{2}/K_{z}^{2}}$. Correspondingly, the DOS has an inverse square-root divergence at the top of the band: $\rho(\omega)\propto 1/\sqrt{E_{\pi}-\omega}$.

Next, we investigate the local DSSF, which requires us to compute the matrix elements defined in Eq. (25). For momenta $0<\rvert k\rvert<k_{0}$, we have

$\displaystyle{\cal M}_{b_{j}}(k)\simeq\sqrt{\frac{2}{L}}e^{-ikj},\ \ {\cal M}_{c_{j}}(k)\simeq 0,$

(32)

while for $k_{0}<\rvert k\rvert<\pi$, we obtain

	$\displaystyle{\cal M}_{b_{j}}(k)$	$\displaystyle\simeq$	$\displaystyle\frac{i}{\sqrt{L}}e^{-ikj+i\phi_{k}/2},$		(33)
	$\displaystyle{\cal M}_{c_{j}}(k)$	$\displaystyle\simeq$	$\displaystyle\frac{1}{\sqrt{L}}\frac{E_{k}}{2h_{z}}e^{-ikj+i\phi_{k}/2},$		(34)

where $\phi_{k}$ is a $k$-dependent phase factor. Notably, these matrix elements are explicit functions of the boundary-mode energies $E_{k}$. From these results, we obtain the Majorana-fermion spectral functions in Eq. (24):

	$\displaystyle{\cal A}_{bb}(\omega)$	$\displaystyle=$	$\displaystyle\frac{2}{3}\delta(\omega)+\rho(\omega),$		(35)
	$\displaystyle{\cal A}_{cc}(\omega)$	$\displaystyle=$	$\displaystyle\left(\frac{\omega}{2h_{z}}\right)^{2}\rho(\omega),$		(36)
	$\displaystyle{\cal A}_{bc}(\omega)$	$\displaystyle=$	$\displaystyle-{\cal A}_{cb}(\omega)=i\left(\frac{\omega}{2h_{z}}\right)\rho(\omega).$		(37)

One can verify that $\int_{-\infty}^{\infty}d\omega{\cal A}_{bb}(\omega)=1$, obeying the sum rule. In contrast, the integrated spectral weight of ${\cal A}_{cc}(\omega)$ is less than 1, as some weight is transferred to the bulk modes. The $\delta$-peak in ${\cal A}_{bb}$ arises from the quasi-flat modes at $\rvert k\rvert<k_{0}$, which, in reality, have a finite broadening $\propto\kappa$ that is neglected in the small-field analysis above. The factors $\omega/(2h_{z})$ in ${\cal A}_{cc}$, ${\cal A}_{bc}$, and ${\cal A}_{cb}$ originate from the factor $E_{k}/(2h_{z})$ in the matrix element ${\cal M}_{c_{j}}$ given in Eq. (34).

The local DSSF is finally obtained by performing the convolution described in Eq. (19):

	$\displaystyle S_{j}^{zz}(\omega)$	$\displaystyle=$	$\displaystyle\frac{2}{3}\left(\frac{\omega}{2h_{z}}\right)^{2}\rho(\omega)$		(38)
		$\displaystyle+$	$\displaystyle\frac{1}{2}\int_{0}^{\omega}d\nu\rho(\nu)\rho(\omega-\nu)\left(\frac{\omega-2\nu}{2h_{z}}\right)^{2},$

where the first term arises from the convolution of the $\delta$-peak in ${\cal A}_{bb}$ with ${\cal A}_{cc}$, while the second term is connected to two-particle excitations from the hybridized mode at $\rvert k\rvert>k_{0}$.

Since $\rho(\omega)\propto\omega$ for small energies, the first term in Eq. (38) dominates over the second one, leading to

$$S_{j}^{zz}(\omega)\propto\omega^{3}.$$

(39)

Moreover, the first term gives an inverse-square-root singularity at the top of the band, corresponding to

$$S_{j}^{zz}(\omega)\propto 1/\sqrt{E_{\pi}-\omega}.$$

(40)

In contrast, the second term results in a regular continuum without any singularities in the entire energy range.

In Fig. 3(g), we show the local DSSF for a finite magnetic field. While the low-energy $\omega^{3}$ behavior is rigorously satisfied, the singularity near the band top is smeared out over the interval $\rvert\omega-E_{\pi}\rvert\lesssim\kappa$ because the quasi-flat mode, corresponding to the $\delta$-peak in Eq. (35), has a finite broadening $\propto\kappa$. Nevertheless, the inverse-square-root behavior still approximately applies in a finite energy range below this interval.

Next, we discuss the non-chiral phase of the Kitaev model, where the Kitaev spin exchange is stronger along one of the three bond directions. As shown in Figs. 3(b) and 3(c), there are two possible scenarios depending on the orientation of the strongest bonds relative to the boundary. These two scenarios exhibit distinct boundary excitation spectra, which can be understood by examining the spectrum in the absence of the hybridization term $\hat{\tilde{\mathscr{H}}}_{Z}$. In the first scenario, where the strongest bonds are perpendicular to the boundary, the boundary spectrum consists of one branch of dangling fermions and another branch of matter fermions [black dotted lines in Fig. 3(e)]. Note that the existence of a matter-fermion boundary mode is linked to boundary sites that are not connected to any dominant bonds as creating a matter-fermion excitation on a dominant bond incurs an energy cost comparable to the largest Kitaev spin exchange. In contrast, for the second scenario, where the strongest bonds are at a $30^{\circ}$ angle with respect to the boundary, all lattice sites belong to a dominant bond, and the boundary spectrum thus contains only a single branch of dangling fermions [black dotted line in Fig. 3(f)]. Upon introducing the hybridization term $\hat{\tilde{\mathscr{H}}}_{Z}$, the two boundary modes in the former case become gapped, whereas the boundary mode in the latter case remains gapless, as indicated by the colored lines in Figs. 3(e) and 3(f).

To make contact with the boundary spectrum of the chiral KSL, we approach the non-chiral phase from the chiral phase by following one of the paths shown in Fig. 4. Along path I, $k_{0}$ decreases monotonically and vanishes at the phase boundary into the non-chiral phase. Therefore, the dispersive hybridized mode, which exists at $\rvert k\rvert>k_{0}$ in the chiral KSL phase, spans the entire BZ at the phase boundary. In contrast, if we approach the non-chiral phase from path II-1 or II-2, $k_{0}$ increases to $\pi$ at the phase boundary, and the quasi-flat mode at $\rvert k\rvert<k_{0}$ extends over the whole BZ accordingly.

In the following, we analyze the two scenarios of the non-chiral KSL separately. We start from the scenario in Fig. 3(b) and the corresponding numerical results in Fig. 3(e,h). In this case, the analytical expressions for the small-field limit of the chiral phase remain applicable. However, they reveal different behaviors for both the Majorana-fermion spectral function and the local DSSF. First, the energy spectrum in Eq. (28) becomes gapped due to ${\rm max}(\lambda_{k}^{2})=\lambda_{0}^{2}=(K_{x}+K_{y})^{2}/K_{z}^{2}<1$. In particular, the spectrum is bounded, $E_{0}\leq E_{k}\leq E_{\pi}$, by the two finite energies

	$\displaystyle E_{0}=2h_{z}\sqrt{1-(K_{x}+K_{y})^{2}/K_{z}^{2}},$		(41)
	$\displaystyle E_{\pi}=2h_{z}\sqrt{1-(K_{x}-K_{y})^{2}/K_{z}^{2}}.$		(42)

The gapped nature of the spectrum distinguishes the non-chiral phase from the chiral phase where a gapless chiral mode is guaranteed by topology. Due to the quadratic dispersions near the band bottom $E=E_{0}$ and the band top $E=E_{\pi}$, the density of states (DOS) of the boundary modes exhibits inverse-square-root singularities:

$\displaystyle\rho(\omega)\propto\frac{1}{\sqrt{\omega-E_{0}}},\quad\rho(\omega)\propto\frac{1}{\sqrt{E_{\pi}-\omega}}.$

(43)

The Majorana-fermion spectral function ${\cal A}_{F_{m}F_{n}}(\omega)$ is then confined to the energy interval $E_{0}\leq\omega\leq E_{\pi}$. Moreover, since the matrix elements ${\cal M}_{F_{m}}(k)$ in Eq. (25) remain finite in the limits of $\omega\to E_{0,\pi}$, the spectral function exhibits the same singular behavior as $\rho(\omega)$.

The convolution of the spectral functions in Eq. (19) leads to a two-particle continuum distributed within the energy range $2E_{0}<\omega<2E_{\pi}$. While the convolution produces finite values for $S_{j}^{zz}({2}E_{0})$ and $S_{j}^{zz}({2}E_{\pi})$, a logarithmic singularity appears at $\omega=E_{0}+E_{\pi}$:

$\displaystyle S_{j}^{zz}(\omega)\propto\log\biggr{\rvert}\frac{1}{\omega-(E_{0}+E_{\pi})}\biggr{\rvert}.$

(44)

These features obtained to the leading order in $h_{\mu}$ explain the local DSSF presented in Fig. 3(h).

Next, we consider the alternative scenario in Fig. 3(c), which corresponds to the numerical results in Fig. 3(f,i). As we discussed above, the dangling fermions can only hybridize with each other as well as the gapped bulk modes, resulting in an effective tight-binding model along the boundary. This effective tight-binding spectrum is gapless, resembling the chiral boundary mode. We also note that the effective hopping amplitude is only finite in the absence of time-reversal symmetry and is not captured by the previous small-field analysis, which only gives rise to a zero-energy mode with the wave function dominated by the dangling fermions [see Eq. (27)].

We calculate the local DSSF numerically for a finite-size lattice. As shown in Fig. 3(i), the DSSF exhibits $\omega^{3}$ behavior at small energies, resembling the chiral phase. Additionally, a broad peak appears at the energy corresponding to the top of the boundary band (denoted as $E_{\rm max}$). The spectral intensity is then gradually reduced upon further increasing the energy and vanishes above $2E_{\rm max}$.

To summarize, the analysis of the chiral and non-chiral KSLs based on the effective low-energy Hamiltonian in Eq. (11) reveals distinct spectral features in the local DSSF. The spectrum remains gapless in the chiral KSL due to the presence of a topologically protected chiral boundary mode. In contrast, for the non-chiral KSL, the spectral response can be either gapped or gapless, contingent on the boundary orientation relative to the dominant bond interaction.

While the characteristic lineshapes observed in both cases can serve as useful indicators, similar low-energy spectral distributions can arise in both chiral and non-chiral KSLs. Moreover, when comparing the local DSSF for the three scenarios in Fig. 3(g,h,i), the distinction between gapped and gapless spectra can become obscured under realistic conditions due to the inevitable presence of disorder. This lack of a clear distinction motivates us to investigate the impact of boundary defects, as their presence may help uncover spectral features unique to the chiral boundary mode, distinguishing it from other types of boundary modes in, e.g., the non-chiral KSL.

We first consider the chiral KSL, and start the investigation by examining the Majorana-fermion spectral functions. In the absence of hybridization, the dangling fermions are localized zero-energy modes, manifesting as sharp $\delta$-peaks at $\omega=0$ in the local spectral function ${\cal A}_{b_{j},b_{j}}(\omega)$, with a total spectral weight of 1. If hybridization is then introduced via Eq. (7), the spectral weight of the dangling fermions spreads across the energy range of the chiral mode, as demonstrated for the clean boundary. The resulting spectrum includes a peak near zero energy—a remnant of the zero-energy peak for bare dangling fermions—and a broad continuum extending up to the top of the boundary band dispersion.

In the presence of disorder, ${\cal A}_{b_{j},b_{j}}(\omega)$ develops a series of peaks superimposed on a continuum background, as shown in Fig. 5(a) for disorder strength $p=0.3$. Importantly, a persistent peak always appears around zero energy, originating from the quasi-flat mode at $\rvert k\rvert<k_{0}$ of the clean boundary. In contrast, the finite-energy peaks appear at different energies depending on the location along the boundary, and correspond to resonances with intrinsic broadening due to the linear-in-$h$ hybridization with the chiral boundary mode. Indeed, by fitting these peaks using a Gaussian function, we find that the peak width increases linearly with $h$, as shown in the inset of Fig. 5(a). The real-space wave functions of the corresponding resonance modes are peaked at the location where the given Majorana fermion is excited but extend across the entire lattice [see Fig. 6(a)]. It is important to distinguish the broadening discussed above from the one caused by coupling with the bulk modes. The latter occurs above the bulk gap indicated by the black dashed line in Fig. 5(a).

Next, we discuss the local DSSF at a boundary site, which is the convolution of two Majorana-fermion spectral functions [see Eq. (19)]. As shown in Fig. 5(b), the main features of the local DSSF can still be captured by the convolution involving the narrow $\delta$-peak in ${\cal A}_{b_{j},b_{j}}(\omega)$ around zero energy, similar to the first term in Eq. (38) for a clean boundary. Therefore, the local DSSF is approximately proportional to ${\cal A}_{c_{j},c_{j}}(\omega)$, whose structure is similar to that of ${\cal A}_{b_{j},b_{j}}(\omega)$. This particular convolution gives rise to peaks around the same positions as in the Majorana-fermion spectral function. Notably, the peak width in the local DSSF is well above the energy resolution used in the calculation and increases linearly with $h$. Note that there are additional peaks originating from the convolution of two finite-energy resonance modes, which could be above the bulk gap. For these peaks, however, it is hard to distinguish the broadening due to the chiral boundary mode and that due to the bulk modes since the dangling fermions couple to both of them through the same Zeeman term. Therefore, for the purpose of detecting the chiral boundary mode, one must focus on the structure of the DSSF below the bulk gap.

It is also worth mentioning that a finite-energy resonance mode already appears around a $60^{\circ}$ corner of the lattice, i.e., the intersection of two zigzag boundaries [see Fig. 6(a)] in the absence of disorder. The convolution of the zero-energy peak in ${\cal A}_{b_{j},b_{j}}(\omega)$ with this “corner mode” leads to an extra peak in the local DSSF.

The calculation above assumes that STM has atomic-scale spatial resolution. In a real experiment, however, the effective tip resolution spans a few lattice sites. To account for a finite tip size, we average the local DSSF over five adjacent boundary sites. As shown in Fig. 7(a), we find that the resonance peaks remain resolved for disorder strength $p=0.3$. Moreover, both the peak position and the peak width continue to exhibit linear scaling with the field $h$. We finally note that, by averaging over progressively more boundary sites, one eventually crosses over into the disorder-averaged regime discussed in the companion paper [53], where the resonance peaks corresponding to different boundary positions coalesce into a single peak that continues to exhibit linear energy scaling.

To contrast with the chiral KSL, we next consider the non-chiral KSL. Following the same approach, we first examine the Majorana-fermion spectral functions. Similar to the chiral case, ${\cal A}_{b_{j},b_{j}}(\omega)$ displays several peaks at discrete energy levels, as shown in Fig. 5(c). Unlike in the chiral case, however, these peaks do not exhibit any intrinsic broadening. Instead, the Gaussian fit indicates that the peak width is determined solely by the energy resolution used in the calculation, meaning that the peak width remains constant as a function of $h$, as long as the energy is below the bulk gap. We also observe that the wave functions of the corresponding modes are localized in real space, as illustrated in Fig. 6(b).

Figure 5(c) corresponds to the relative boundary orientation shown in Fig. 3(c), in which case there is a gapless boundary mode in the clean boundary limit. Upon introducing disorder, this gapless boundary mode becomes localized, giving rise to peaks in ${\cal A}_{b_{j},b_{j}}(\omega)$ near zero energy. In contrast, for the alternative boundary considered in Fig. 3(b), where the boundary modes are gapped in the clean boundary limit, the localized modes only appear at finite energy. However, this distinction does not produce observable differences in the local DSSF.

The convolution of these localized boundary modes in the Majorana-fermion spectral function results in a series of discrete $\delta$-peaks in the local DSSF. The energies of these peaks are found by summing the energies of any two distinct $\delta$-peaks in the Majorana-fermion spectral function, in accordance with the Pauli exclusion principle. Therefore, the local DSSF shown in Fig. 5(d) does not exhibit a peak around zero energy.

As the field strength increases, the energies of these modes rise linearly, while their peak widths remain constant, determined solely by the energy resolution $\eta=0.005{\rm max}(K_{\alpha})$ (i.e., the imaginary part of the complex frequency $\omega+i\eta$) used in the calculation. This behavior in the non-chiral case is clearly distinct from the chiral case, where the peak width increases linearly with the field strength. Moreover, we confirm that this distinction remains observable after averaging over multiple lattice sites within the effective tip size [see Fig. 7(b)]. Finally, the “corner mode” mentioned earlier becomes a truly localized mode for the non-chiral KSL, which implies that its convolution with any other mode localized due to disorder results in an additional $\delta$-peak in the local DSSF.