Oriented propagation of magnetization due to chiral edge modes in Kitaev-type models

TL;DR

This paper theoretically investigates the dynamics of chiral edge modes in Kitaev-type models, revealing oriented edge magnetization propagation and providing insights into spin and charge transport in topological materials.

Contribution

It introduces a method to analyze chiral edge mode dynamics in exactly solvable Kitaev models, highlighting oriented magnetization propagation without relying on conserved quantities.

Findings

Edge magnetization propagates along the edge in a chiral manner.

The approach applies to both spin transport in magnets and charge transport in superconductors.

The method does not depend on conserved quantities, broadening its applicability.

Abstract

Detecting chiral edge modes in topological materials has been intensively pursued in experiments. However, the phenomena caused by the modes are not yet elucidated theoretically. We study the dynamics of chiral spinon wave packets at the edge in Kitaev-type magnets. More precisely, by relying on the exact solvability of the models, we construct a spinon wave packet, localized edge magnetization, which shows oriented propagation along the edge, whose behavior is expected from the chiral character of the dispersion relation of the chiral edge modes. In general, this approach enables us to study not only spin transport in anisotropic magnets but also charge transport in Bogoliubov-de Gennes-type superconductors because it does not rely on a conserved quantity.

Equations

$$H=H_{\mathrm{K}}+H_{3-\mathrm{spin}},$$

$$H_{\mathrm{K}}=J_x\sum _{\ell =1}^{L_x}\sum _{m=1}^{L_y}{\sigma }_{(2\ell -1,m)}^x{\sigma }_{(2\ell ,m)}^x+J_y\sum _{\ell =1}^{L_x-1}\sum _{m=1}^{L_y}{\sigma }_{(2\ell ,m)}^y{\sigma }_{(2\ell +1,m)}^y+J_z\sum _{\ell =1}^{L_x}\sum _{m=1}^{L_y}{\sigma }_{(2\ell ,m)}^z{\sigma }_{(2\ell -1,m+1)}^z,$$

$$H_{3-\mathrm{spin}}$$

$${\sigma }_{(\ell ,m)}^{+}=2a_{(\ell ,m)}{e}^{i\pi {\hat{\theta }}_{(\ell ,m)}},$$

$${\sigma }_{(\ell ,m)}^{-}=2e^{i\pi {\hat{\theta }}_{(\ell ,m)}}{a}_{(\ell ,m)}^{†},$$

$${\sigma }_{(\ell ,m)}^z=(-1{)}^{\ell }\left[2a_{(\ell ,m)}^{†}{a}_{(\ell ,m)}-1\right],$$

$${\hat{\theta }}_{(\ell ,m)}=\sum _{m^{\prime }<m}\sum _{{\ell }^{\prime }=1}^{2L_x}{a}_{({\ell }^{\prime },m^{\prime })}^{†}{a}_{({\ell }^{\prime },m^{\prime })}+\sum _{{\ell }^{\prime }<\ell }{a}_{({\ell }^{\prime },m)}^{†}{a}_{({\ell }^{\prime },m)}.$$

$$c_{(\ell ,m)}=$$

$$d_{(\ell ,m)}=$$

$$c_{(\ell ,m)}=$$

$$d_{(\ell ,m)}=$$

$$H_{\mathrm{K}}=i{J}_x\sum _{\ell =1}^{L_x}\sum _{m=1}^{L_y}{c}_{(2\ell -1,m)}{c}_{(2\ell ,m)}+i{J}_y\sum _{\ell =1}^{L_x-1}\sum _{m=1}^{L_y}{c}_{(2\ell ,m)}{c}_{(2\ell +1,m)}+J_z\sum _{\ell =1}^{L_x}\sum _{m=1}^{L_y}{c}_{(2\ell ,m)}{c}_{(2\ell -1,m+1)}{d}_{(2\ell ,m)}{d}_{(2\ell -1,m+1)},$$

$$H_{3-\mathrm{spin}}$$

$$H_{\mathrm{K}}=i{J}_x\sum _{\ell =1}^{L_x}\sum _{m=1}^{L_y}{c}_{(2\ell -1,m)}{c}_{(2\ell ,m)}+i{J}_y\sum _{\ell =1}^{L_x-1}\sum _{m=1}^{L_y}{c}_{(2\ell ,m)}{c}_{(2\ell +1,m)}+i{J}_z\sum _{\ell =1}^{L_x}\sum _{m=1}^{L_y}{c}_{(2\ell ,m)}{c}_{(2\ell -1,m+1)},$$

$$H_{3-\mathrm{spin}}$$

$${\alpha }_{(\ell ,m)}=\frac{c_{(2\ell -1,m)}+i{c}_{(2\ell ,m)}}{2},$$

$${\alpha }_{\ell ,k_y}=\frac{1}{\sqrt{L_y}}\sum _{m=1}^{L_y}{e}^{i{k}_{y}m}{\alpha }_{(\ell ,m)}.$$

$$H=\sum _{k_y}\sum _{\ell ,{\ell }^{\prime }}{\mathbf{\Psi }}_{\ell }^{†}(k_y)[\hat{h}(k_y){]}_{\ell ,{\ell }^{\prime }}{\mathbf{\Psi }}_{{\ell }^{\prime }}(k_y),$$

$${\mathbf{\Psi }}_{\ell }(k_y)={\left({\alpha }_{\ell ,k_y},{\alpha }_{\ell ,-k_y}^{†}\right)}^{\mathrm{T}},$$

$$\displaystyle[\hat{h}(k_{y})]_{\ell,\ell^{\prime}}=\left(\begin{array}[]{cc}[h_{0}(k_{y})]_{\ell,\ell^{\prime}}&[\Delta(k_{y})]_{\ell,\ell^{\prime}}\\ \left[\Delta^{\dagger}(k_{y})\right]_{\ell,\ell^{\prime}}&[-h_{0}(-k_{y})]_{\ell,\ell^{\prime}}\\ \end{array}\right).$$

$$[h_0(k_y){]}_{\ell ,{\ell }^{\prime }}=\left(J_x-J_z\cos k_y\right){\delta }_{\ell ,{\ell }^{\prime }}-\frac{J_y}{2}\left({\delta }_{\ell ,{\ell }^{\prime }+1}+{\delta }_{\ell ,{\ell }^{\prime }-1}\right),$$

$$[\Delta (k_y){]}_{\ell ,{\ell }^{\prime }}$$

$$H=$$

$$+$$

$$H=$$

$$+$$

$$H=$$

$$+$$

Full text

I Introduction

The exact ground state of the Kitaev honeycomb model Kitaev2006 is known as a quantum spin liquid Anderson1973 ; Balents2010 ; Zhou2017 ; Knolle2019 that shows short-range spin-spin correlations Baskaran2007 . As it is the first example of the exactly solvable model of a quantum spin liquid in two dimensions, tremendous theoretical studies have been done to reveal fundamental properties of this model such as thermal quantities and spin dynamics Knolle2014 ; Nasu2014 ; Nasu2015 ; Knolle2015 ; Yoshitake2017 ; Yoshitake2017_2 ; Nasu2018 ; Gohlke2018 ; Motome2019 . In parallel with theoretical studies, experimental attempts to seek the Kitaev spin liquids in candidate materials Jackeli2009 ; Singh2010 ; Singh2012 ; Plumb2014 ; Rau2014 ; Trebst2017 are ongoing Sandilands2015 ; Banerjee2016 ; Do2017 ; Banerjee2017 ; Winter2017 ; Hermanns2018 .

In addition to the nature of the quantum spin liquid, the model is a weak topological material that hosts a Majorana edge flat band Thakurathi2014 which yields unidirectional edge magnetization Mizoguchi2019 . When applying an external magnetic field along [111] direction, the ground state of the model is believed to change to a chiral spin liquid which is characterized by a Chern number Kitaev2006 , and exhibits chiral edge modes although the exact solvability is lost. If three spin interactions are added to the Hamiltonian of the Kitaev honeycomb model instead of the external magnetic field, then the exact solvability is retained, and the resulting model shows the character of a chiral spin liquid. A different variant of the model in the same topological class can be realized as an exactly solvable Kitaev-type model on a triangle-honeycomb lattice Yao2007 ; Chung2010 . Thus, several Kitaev-type models were found to have the character of a chiral spin liquid.

Detecting chiral edge modes Katsura2010 ; Diaz2019 in topological materials has been a central issue in experiments as well. Actually, Kitaev-type magnets Singh2010 ; Singh2012 ; Plumb2014 in the external magnetic field are most likely candidates whose ground state is a chiral spin liquid Kasahara2018 ; Kasahara2018_2 . Therefore, a much deeper theoretical understanding of the related phenomena is required. In this paper, we study propagation of edge magnetization due to the chiral edge modes in Kitaev-type magnets. More precisely, by relying on the exact solvability of the models, we construct a spinon wave packet, localized edge magnetization, which is an excited state composed of the chiral edge states. Because of the chirality of the dispersion relation of the chiral edge mode, the wave packet propagates in only one of the two directions of the edge irrespectively of the initial form of the wave packet. In comparison with previous studies addressing spin transport which use spin current Chen2013 ; Chatterjee2015 ; Hirobe2017 ; Carvalho2018 , the advantage of our method is that we need neither to introduce spin current nor to apply an external driving force such as a potential difference. Thus, our method can be applied to not only anisotropic magnets in which non of the spin components commutes with the Hamiltonian but also Bogoliubov-de Gennes (BdG)-type superconductors Sato2016 in which the charge is not conserved. We can expect that such oriented propagation of a local excitation at a sample edge is experimentally detectable.

For the purpose of actually calculating the chiral edge modes, we also introduce a new technique into the transfer matrix method. Although the method is already established as a standard tool Lee1981 ; Hatsugai1993 ; Hatsugai1993_2 ; Dwivedi2016 , known exact solutions of chiral edge states are still rare so far; the examples include the Bernevig-Hughes-Zhang model Mao2010 ; Dwivedi2016 and the modified Haldane and Kane-Mele models Kunst2019 . Even in the case of the well-known Haldane tight-binding model on the honeycomb lattice Haldane1988 ; Huang2012 ; Doh2014 ; Pantaleon2017 ; Kunst2019 , the exact solution of the chiral edge state has not yet been obtained. In this paper, we demonstrate that our new technique enables us to obtain exact solutions of edge states, which is applicable to generic models described by $4\times 4$ transfer matrices, including the Haldane model and the Kitaev model with three-spin interaction which we will deal with in this paper.

The rest of this paper is organized as follows. In Sec. II, we describe the models considered in this paper, namely the Kitaev honeycomb model with three-spin interaction and the Kitaev triangle-honeycomb model. We also explain the formalism of the Jordan-Wigner transformation by which these quantum spin models are mapped to the fermion models whose Hamiltonian have BdG form which is interpreted as the Hamiltonian of a spinless superconductor. In Sec. III, we derive the exact dispersion relation and the wave function of the chiral edge spinons by using the transfer matrix method. In Sec. IV, by relying on the exact solution thus obtained, we construct an excited state which realizes the oriented propagation of the edge magnetization. We also demonstrate that a wave packet constructed by this formulation indeed moves in only one of the two directions of the edge. In Sec. V, the edge dynamical spin susceptibility is obtained as well. This quantity reflects the chiral nature of the edge mode, and is experimentally detectable. A summary and discussion are given in Sec. VI. Appendix A is devoted to the details of the transfer matrix method used in Sec. III. Finally, in Appendix B, we present the details of the derivation of the exact solution of the chiral edge modes in the Kitaev triangle-honeycomb model.

II Model

In order to demonstrate that chiral edge modes in topological materials cause oriented propagation of local edge magnetization in Kitaev-type magnets, we treat two models, namely, the Kitaev honeycomb model with three-spin interactions, and the Kitaev triangle-honeycomb model. In this section, we rewrite these Hamiltonians into the BdG form which is interpreted as the Hamiltonian of a spinless superconductor. If a reader is familiar with the following formulations, it is possible to skip this section. Of course, in order to understand the subsequent calculations, it is necessary to keep the notations in this section in mind.

II.1 Kitaev honeycomb model with three-spin interaction

We consider the following Hamiltonian on a cylindrical geometry, whose sites are denoted by $(\ell,m)$ with $\ell=1,\cdots,2L_{x}$ and $m=1,\cdots,L_{y}$ (Fig. 1):

[TABLE]

where

[TABLE]

and

[TABLE]

We impose the periodic boundary condition, $\bm{\sigma}_{(\ell,m+L_{y})}=\bm{\sigma}_{(\ell,m)}$, in the vertical direction, while we impose the open boundary condition in the horizontal direction. Here, $H_{\mathrm{K}}$ represents the Kitaev honeycomb model Kitaev2006 , where the real parameters, $J_{x}$, $J_{y}$, and $J_{z}$, are the exchange coupling on $x$, $y$, and $z$ bonds, respectively. As to the classification of the bonds into three types, see Fig. 1 and its caption. The Hamiltonian $H_{3\mathrm{-spin}}$ represents the three-spin interaction with real coupling $J^{\prime}$. Although this term looks fairly artificial, it was shown by Kitaev that the interaction Hamiltonian $H_{3\mathrm{-spin}}$ naturally appears when the Zeeman energy of the external magnetic field in the [111] direction is treated within the third-order perturbation in the Kitaev honeycomb model Kitaev2006 . In the case with $J_{x}=J_{y}=J_{z}=J$, the parameter $J^{\prime}$ can be estimated as $J^{\prime}\sim\frac{h^{3}}{J^{2}}$, where $h$ is the strength of the magnetic field. Therefore, the Kitaev honeycomb model with the Zeeman term is believed to belong to the same topological class of the chiral spin liquid as that of the Kitaev honeycomb model with the three-spin interactions. This motivates us to study the chiral edge mode in the Kitaev honeycomb model with the three-spin interactions.

To obtain the exact ground state of the Hamiltonian of Eq. (1), we map the spin model to the fermion model, by using the Jordan-Wigner transformation. Following Refs. Feng2007, ; Yao2007, ; Chen2007, ; Chen2008, , we introduce a complex fermion, $a_{(\ell,m)}$, and rewrite the spin operators such that

[TABLE]

where $\sigma_{i}^{\pm}=\sigma_{i}^{x}\pm i\sigma_{i}^{y}$, and

[TABLE]

Further, we decompose the complex fermion $a_{(\ell,m)}$ into real and imaginary parts, i.e., we introduce two Majorana fermions in the following manner Feng2007 ; Yao2007 ; Chen2007 ; Chen2008 :

[TABLE]

and

[TABLE]

Then, by using $c$ and $d$, we can rewrite $H_{\mathrm{K}}$ and $H_{3\mathrm{-spin}}$ as

[TABLE]

and

[TABLE]

One finds that the pairs of $d_{i}d_{j}$ in Eqs.(8) and (9) commute with $H_{\mathrm{K}}$ and $H_{3\mathrm{-spin}}$, thus they are the conserved quantities. Consequently, we can replace $d_{i}d_{j}$ with the classical number, as $d_{i}d_{j}=\pm i$. In the following, we call $d_{i}d_{j}$ a link variable. In the ground state, the configuration of the link variables is determined so that the energy for the Hamiltonian for $c$-Majorana fermions is minimized; this can be achieved by $d_{(2\ell,m)}d_{(2\ell-1,m+1)}=+i$ for $\ell=1,\cdots,L_{x}$ and $m=1,\cdots,L_{y}$ Kitaev2006 ; Yao2007 . Then, the Hamiltonian thus obtained can be written in the quadratic form of the $c$-Majorana fermions,

[TABLE]

and

[TABLE]

To rewrite the Hamiltonian of Eqs. (10) and (11) in the BdG form, we introduce a complex fermion Mizoguchi2019 ,

[TABLE]

and perform the Fourier transformation in the vertical direction,

[TABLE]

Using $\alpha_{\ell,k_{y}}$, we can write the Hamiltonian as

[TABLE]

with

[TABLE]

and

[TABLE]

$h_{0}(k_{y})$, $\Delta(k_{y})$, and $\Delta^{\dagger}(k_{y})$ are defined as

[TABLE]

[TABLE]

and $\left[\Delta^{\dagger}(k_{y})\right]_{\ell,\ell^{\prime}}=\{[\Delta(k_{y})]_{\ell^{\prime},\ell}\}^{\ast}$.

The topological class of the BdG Hamiltonian of Eq. (14) can be determined by the symmetries. For $J^{\prime}=0$, it belongs to the BDI class, which has time-reversal, particle-hole, and chiral symmetries Mizoguchi2019 . As a result, the Majorana edge flat band appears in the Ay phase Kitaev2006 ; Mizoguchi2019 ; Thakurathi2014 , which is protected by the weak topological nature. If $J^{\prime}$ is finite, then the particle-hole symmetry is kept, but the time-reversal symmetry is broken. Consequently, the chiral symmetry given by the product of the time-reversal and particle-hole symmetries is broken, too. Therefore, the corresponding topological class is D class, whose topological number is given by $\mathbb{Z}$, or the Chern number Kitaev2006 . For this class, the chiral edge mode appears if the Chern number is nonzero.

II.2 Kitaev triangle-honeycomb model

Next, we describe the Kitaev triangle-honeycomb model Yao2007 , which is another example of the Kitaev-type model having the chiral edge modes.

We consider a system with $6L_{x}\times L_{y}$ sites, and impose a periodic boundary condition, $\bm{\sigma}_{(\ell,m+L_{y})}=\bm{\sigma}_{(\ell,m)}$, in the vertical direction. The Kitaev triangle-honeycomb model is given as

[TABLE]

See Fig. 2 for the geometry and the interactions. We perform the Jordan-Wigner transformation of Eqs. (4a)-(4c), and decompose the complex fermion $a_{i}$ into two Majorana fermions $c_{i}$ and $d_{i}$ in the same way as that in Eqs. (6) and (7). Using these transformations for Eq. (21), we obtain

[TABLE]

Again, the pairs $d_{i}d_{j}$ in Eq. (22) commute with the Hamiltonian and thus can be replaced with classical numbers. We set the link variables as $d_{(6\ell-2,m)}d_{(6\ell,m)}=d_{(6\ell-1,m)}d_{(6\ell-4,m)}=+i$, and $d_{(6\ell-5,m)}d_{(6\ell-3,m)}=-i$, which gives us the ground state Yao2007 . This choice of the link variables leads to the time-reversal symmetry breaking. In fact, the time-reversal operation changes the sign of the link variables ($\pm i\rightarrow\mp i$), which gives another state with the same energy. Therefore, the ground state has at least twofold degeneracy due to the time-reversal symmetry breaking Yao2007 . Then, the Hamiltonian for $c$-Majorana fermions is

[TABLE]

To rewrite Eq. (23) in the BdG form, we define three species of complex fermions:

[TABLE]

and their Fourier transformations,

[TABLE]

for $s=A,B,C$. Then, the Hamiltonian of Eq. (23) is rewritten as

[TABLE]

where $\bm{\Psi}_{\ell}(k_{y})=(\alpha_{\ell,k_{y},A},\alpha_{\ell,k_{y},B},\alpha_{\ell,k_{y},C},\alpha^{\dagger}_{\ell,-k_{y},A},\alpha^{\dagger}_{\ell,-k_{y},B},\alpha^{\dagger}_{\ell,-k_{y},C})^{\rm T}$ and

[TABLE]

with

[TABLE]

and

[TABLE]

The topological class of this Hamiltonian is also D class, and hence the model exhibits a chiral edge mode.

III Dispersion relation of edge mode

III.1 Kitaev honeycomb model

In this section, we derive the dispersion relation of the chiral edge mode at the left edge for the Kitaev model. To obtain the exact solution of the chiral edge mode, we employ a transfer matrix method. As we mentioned in Sec. I, our method is applicable to generic models described by $4\times 4$ transfer matrices.

The left edge mode $\gamma^{\mathrm{L}}(k_{y})$ can be expanded in terms of $\bm{\Psi}_{\ell}(k_{y})$ as

[TABLE]

with

[TABLE]

Note that $\gamma^{\mathrm{L}}(k_{y})$ satisfies $\left[\gamma^{\mathrm{L}}(k_{y}),H\right]=\varepsilon^{\mathrm{L}}(k_{y})\gamma^{\mathrm{L}}(k_{y})$. This leads to the eigenvalue equation for $\bm{\varphi}_{\ell}(k_{y})$, which is give by

[TABLE]

where $\hat{A}(k_{y})$ and $\hat{B}(k_{y})$ are $2\times 2$ matrices given as

[TABLE]

and

[TABLE]

Equation (46) can be rewritten as

[TABLE]

in terms of the $4\times 4$ transfer matrix,

[TABLE]

where $\hat{I}_{2}$ stands for the $2\times 2$ identity matrix.

From Eqs. (59) and (63), one finds that the energy eigenvalue $\varepsilon^{\mathrm{L}}(k_{y})$ has to be determined such that $\hat{A}^{-1}(k_{y})[\varepsilon^{\mathrm{L}}(k_{y})-\hat{B}(k_{y})]$ and $-\hat{A}^{-1}(k_{y})\hat{A}^{\dagger}(k_{y})$ have simultaneous eigenstates; see Appendix A for details. After some calculations, we obtain the dispersion relation as

[TABLE]

Notice that, in order that the wave function $\bm{\varphi}_{\ell}$ is localized at the left edge, the eigenvalues of $\hat{T}(k_{y})$, $\lambda(k_{y})$, has to satisfy the condition $|\lambda(k_{y})|<1$.

In Fig. 3, we plot the energy eigenvalues of the Hamiltonian of Eq. (14) (red solid lines), together with the edge-mode dispersion given in Eq. (64) (a blue dashed line), showing quite good agreement with each other. It can also be seen in Fig. 3 that the chiral edge mode has a linear dispersion near $k_{y}=0$. The velocity, which is defined as

[TABLE]

is approximated in the vicinity of $k_{y}=0$ as

[TABLE]

III.2 Kitaev triangle-honeycomb model

Next, let us treat the Kitaev triangle-honeycomb model. Following Ref. Yao2007, , we restrict the parameters to $J_{x}=J_{y}=J_{z}=J$, and $J_{x}^{\prime}=J_{y}^{\prime}=J_{z}^{\prime}=J^{\prime}$. It has been pointed out that the ground state is topologically nontrivial (i.e. having a finite Chern number) for $|J^{\prime}|<\sqrt{3}|J|$, and is trivial (i.e., having zero Chern number) for $|J^{\prime}|>\sqrt{3}|J|$; $|J^{\prime}|=\sqrt{3}|J|$ is a critical point at which the bulk band gap closed.

Although the transfer matrix formalism shown in the previous subsection is also applicable to the present model, it is fairly difficult to obtain the exact solution, since the size of matrices is large. We therefore employ a different method to obtain the edge modes, which is similar to one used in the prior work Mizoguchi2019 ; all the details of that method are described in Appendix B.

Here we look at the numerical results: In Figs. 4(a)-(c), we plot the dispersion relations for the BdG equation of Eq. (26), obtained by the numerical diagonalization of the Hamiltonian. For $(J,J^{\prime})=(1,1)$, there exist chiral edge modes, originating from the nonzero Chern number of the bulk band. Interestingly, even for a nontopological phase with vanishing Chern number [$(J,J^{\prime})=(1,\sqrt{5})$], there still exist the edge states, but they are not chiral, because the velocity around $k_{y}=\pm\pi$ takes the opposite sign to that around $k_{y}=0$. This leads to the localization of the edge modes in the presence of the impurities due to the back scattering, i.e., the edge modes in the nontopological phase are not stable against disorders unlike the chiral edge modes Elbau2002 in the topological phase.

IV Propagation of magnetization at the edge

In this section, we calculate the time evolution of the edge magnetization. To do this, we combine the technique to calculate the edge magnetization in the previous work Mizoguchi2019 with the result of the eigenvalues and eigenvectors of the left edge mode derived in the previous section.

As a initial state of dynamics, we consider the situation where the chiral edge mode are excited with a certain weight function $\rho_{k_{y}}$ which is a function of the vertical momentum $k_{y}$. Physically, this situation can be achieved by applying an external perturbation, e.g., shining laser pulse to the sample magnet whose energy is lower than the bulk energy gap of the magnet. We assume that such an external perturbation is weak enough not to excite the flux sector.

IV.1 Formulation

IV.1.1 Warm-up: A spinless superconductor

Consider first the case that the link variable, $d_{(2\ell,m)}d_{(2\ell-1,m+1)}$, is equal to $+i$ for all $(\ell,m)$. Then, the Hamiltonian (26) can be also interpreted as the Hamiltonian of a spinless superconductor. Therefore, before proceeding to the calculation of the edge magnetization in the Kitaev-type models, we compute the propagation of a Majorana edge charge excitation in the superconductor. The calculation is much simpler than those in the Kitaev-type models.

To begin with, we recall

[TABLE]

Substituting (67) into Eq. (42), we obtain

[TABLE]

with $\zeta_{\ell,k_{y}}=u_{\ell,k_{y}}+v_{\ell,k_{y}}$ and $\xi_{\ell,k_{y}}=u_{\ell,k_{y}}-v_{\ell,k_{y}}$.

The time evolution of $\gamma^{\mathrm{L}}(k_{y})$ is given by

[TABLE]

Using $\gamma^{\mathrm{L}}(k_{y},t)$, let us consider the wave packet created by the aforementioned external perturbation,

[TABLE]

with the weight function $\rho_{k_{y}}$ which satisfies $\rho_{-k_{y}}^{\ast}=\rho_{k_{y}}$. One can easily show $\eta[\rho]^{\dagger}(t)=\eta[\rho](t)$. Noting that the dispersion relation given by Eq. (64) satisfies $\varepsilon^{\mathrm{L}}(k_{y})=-\varepsilon^{\mathrm{L}}(-k_{y})$, one obtains

[TABLE]

by substituting Eq. (68) into Eq. (70). Here we have introduced

[TABLE]

For ${\cal A}_{n,k_{y}}$, one can show the relation,

[TABLE]

We further define

[TABLE]

with ${\tilde{\cal A}}_{n,k_{y}}=\rho_{k_{y}}{\cal A}_{n,k_{y}}$. Then, using the above relation (73), as well as $\rho_{-k_{y}}^{\ast}=\rho_{k_{y}}$ and $\varepsilon^{\rm L}(-k_{y})=-\varepsilon^{\rm L}(k_{y})$, one has $\psi_{n,m}^{\ast}(t)=\psi_{n,m}(t)$, i.e., $\psi_{n,m}(t)$ is a real function. Clearly, one has

[TABLE]

The wave function $\psi_{n,m}(t)$ is localized at $n$. In the limit $L_{y}\rightarrow\infty$, the wave function is written as

[TABLE]

Therefore, if the Fourier component ${\tilde{\cal A}}_{n,k_{y}}$ is a smooth function with respect to $k_{y}$, then $\psi_{n,m}(t)$ is localized at $m$, too. Note that

[TABLE]

Therefore, the expectation value of the position $m$ is given by

[TABLE]

where $v^{\mathrm{L}}_{\rm g}(k_{y})$ is given in Eq. (65). This implies that the wave packet propagates with the group velocity

[TABLE]

at the left edge of the sample. Since $v^{\mathrm{L}}_{\rm g}(k_{y})$ in Eq. (65) is always positive (or negative) irrespective of the wave number $k_{y}$, so is the group velocity ${\cal V}_{\rm g}$. This implies that an oriented propagation is realized for the chiral edge mode.

In order to obtain the oriented propagation of the Majorana edge charge excitation, we consider a trial state,

[TABLE]

where $\mathcal{N}$ is the normalization factor. Note that the wave function of Eq. (80) is chosen such that the states with different fermion parities are mixed. This is the standard technique to obtain the finite expectation value of a single Majorana fermion Mizoguchi2019 ; Willans2010 ; Willans2011 , as shown below. The expectation value of the Majorana fermion $c_{(1,m)}$ can be calculated as

[TABLE]

where we have used $\eta^{\dagger}[\rho](t)=\eta[\rho](t)$. Thus, the profile of the excess charge of the Majorana edge fermion is expressed in terms of the amplitude $\psi_{1,m}(t)$ of the wave function at the edge site (1,m), and the charge moves in one of the two directions of the edge with the group velocity ${\cal V}_{\rm g}$ of (79).

IV.1.2 Kitaev-type models

Now, to realize the propagation of the edge magnetization, we need to consider a slightly generic trial state. As we have shown in the previous work Mizoguchi2019 , the operator $e^{i\pi\hat{\theta}_{(1,m)}}$ plays a crucial role for the calculation of the edge magnetization, since the spin operators $\bm{\sigma}_{(1,m)}$ inevitably contain it. Since the operator $e^{i\pi\hat{\theta}_{(1,m)}}$ commutes with the Hamiltonian $H$, the state $e^{i\pi\hat{\theta}_{(1,m)}}|0\rangle$ is still a ground state for the fermion vacuum state $|0\rangle$. But the operation changes the link variable, $d_{(2\ell,m)}d_{(2\ell-1,m+1)}$, for some $m$. This causes the change of the signs of the hopping amplitudes of $J_{z}$ bonds for $c$-Majorana fermions. Of course, this change of the signs can be removed by using the gauge transformation $c_{(2\ell-1,m)}\rightarrow\pm c_{(2\ell-1,m)}$.

In this case, instead of the trial state of (80), we consider a trial excited state,

[TABLE]

where the operator $W$ is defined by

[TABLE]

with a positive weight function $w(m)$ which satisfies the normalization condition,

[TABLE]

We take the wave function $\psi_{n,m}(t)$ of (76) and the weight function $w(m)$ so that the overlap between them at the initial time $t=0$ is large. As we will show below, it is essential to set $w(m)$ as a slowly-decaying function. For instance, if $w(m)$ is set to be constant, then Eq. (84) requires that $w(m)=\frac{1}{\sqrt{L_{y}}}$, which leads to the vanishing of the edge magnetization for $L_{y}\rightarrow\infty$ [see Eq. (87)].

Using the parity conservation for the fermions, one can check that $|\Psi(t)\rangle$ is normalized as

[TABLE]

Here, we also have used that the operators $e^{i\pi\hat{\theta}_{(1,m)}}$ change the link variables and that two states with different link variables are orthogonal with each other, and we have taken the wave function of (76) to be normalized to 1.

Let us compute the expectation value of the $y$-component spin at the site $(1,m)$ which is given by $\langle\Psi(t)|\sigma_{(1,m)}^{y}|\Psi(t)\rangle$. We note that the other two components ($x$ and $z$) have vanishing expectation values for $\ket{\Psi(t)}$ because they contain a single $d$-Majorana fermion operator except for the string operator. To obtain finite values of them, one may need to set the trial wave function to be a linear combination of wave functions with different fermion parities of the $d$-Majorana fermions. However, it seems to be fairly difficult not only to find such a favorable wave function, but also to compute the expectation values. The $y$ component of the spin is written

[TABLE]

in terms of the Majorana fermion $c_{(1,m)}$. In the same way as the above calculation, by using this expression, the expression (75) of $\eta[\rho](t)$ and $\{\eta[\rho](t)\}^{\dagger}=\eta[\rho](t)$, we obtain

[TABLE]

We emphasize that the string operator $e^{i\pi\hat{\theta}_{(1,m)}}$ in $W$ is essential to obtain the finite magnetization, since the string operator $e^{i\pi\hat{\theta}_{(1,m)}}$ changes the eigenvalues of the link variables, $d_{(2\ell,m)}d_{(2\ell-1,m+1)}$, of the fermion vacuum $|0\rangle$ as mentioned above, and yields the vanishing of the magnetization. In other words, the finite expectation value is obtained by canceling the string operator in $\sigma_{(1,m)}^{y}$ with those in $W$. Thus, the profile of the magnetization is expressed in terms of the weight function $w(m)$ and the amplitude $\psi_{1,m}(t)$ of the wave function at the edge site $(1,m)$. The amplitude slowly decays by the weight $w(m)$ with distance. Actually, in the limit $L_{y}\rightarrow\infty$, the magnetization is vanishing in the region where $w(m)\rightarrow 0$ [see Eq. (89) as an example]. This implies that it is impossible to suppress the fluctuation of the link variables through the whole region in this approach. Although the local magnetization decays and finally vanishes, it exhibits an oriented propagation along the left edge with the group velocity ${\cal V}_{\rm g}$.

IV.2 Numerical demonstration

IV.2.1 Kitaev honeycomb model with three-spin interaction

Following the formulation shown in the previous subsection, we calculate the time evolution of the edge magnetization numerically. We first consider the Kitaev honeycomb model with three-spin interaction. We choose the Gaussian distribution,

[TABLE]

for the weight function $\rho_{k_{y}}$ in (70), and

[TABLE]

for the weight function $w(m)$ in (83) [Fig. 5(a)].

We plot $\mu(m,t)$ for $(J_{x},J_{y},J_{z},J^{\prime})=(1,0.8,0.9,0.3)$ and $L_{y}=32$, with $k_{0}=0.1\pi$ in Fig. 5(b). We clearly see that the dip of the magnetization propagates from the top row $(m=L_{y})$ to the bottom row $(m=1)$ in both cases, indicating the oriented propagation of the magnetization due to the chiral edge mode. The velocity of the propagating magnetization is estimated from the figure as $1.08$, which is compatible with the group velocity given in Eq. (66), that is, $v_{\mathrm{g}}^{\rm L}(0)\sim 1.12$. We also see that the magnetization decays as the center of the wave packet approaches to $m=L_{y}/2$.

IV.2.2 Kitaev triangle-honeycomb model

Next, we consider the Kitaev triangle-honeycomb model. Practically, we need to replace $\zeta_{1,k_{y}}$ with $\zeta_{1,k_{y},A}$ when applying the formulation to the Kitaev triangle-honeycomb model. As for the initial state of dynamics, we consider the case where the left edge mode crossing the zero energy is excited. To choose such a weight function, we need to bear in mind that the edge mode in this model exists around $k_{y}=\pm\pi$ in the momentum space, in contrast to the Kitaev honeycomb model with the three-spin interaction, where the edge mode exists around $k_{y}=0$. We therefore choose the weight function $\rho_{k_{y}}$ as

[TABLE]

which has a large amplitude at $k_{y}=\pm\pi$.

In Fig. 6, we plot the results for $J=1$, $J^{\prime}=1$, $L_{y}=32$, and $k_{0}=0.1\pi$. We see that the magnetization in this system shows the staggered structure in the real space, originating from the fact that the weight function has a large amplitude at $k_{y}=\pm\pi$, in contrast to the Kitaev honeycomb model with the three-spin interaction. Nevertheless, we see that the magnetization propagates from the top to the bottom, due to the existence of the chiral edge spinons. Again, we see that the magnetization decays as approaching to $m=L_{y}/2$, while the amplitude is recovered near $m\sim 1$ ($t=80$) since $w(m)$ takes large values in that region.

V Edge dynamical susceptibility

In this section, we calculate the edge structure factor,

[TABLE]

where $|0\rangle$ stands for the ground state, i.e., the state in which all the eigenstates with negative eigenvalues are occupied. We consider the situation that the frequency $\omega$ is contained in the region of the bulk spectral gap so that only the chiral edge mode is detected. Note that we only consider a pair of $\sigma^{y}$s for the expectation value in Eq. (91), since the action of the other two components, $\sigma^{x}$ and $\sigma^{z}$, causes the change of the flux sector Mizoguchi2019 , which makes the calculation fairly complicated. As we will explain below, the chiral nature of the edge mode is reflected to this quantity.

Substituting Eq. (86) into the expression (91) of the susceptibility, we have

[TABLE]

Using the facts that $e^{i\hat{\theta}_{(1,m)}}$ commutes with $H$ and that for any two states of $c$-Majorana fermions the matrix element of $e^{i\hat{\theta}_{(1,m)}}e^{i\hat{\theta}_{(1,m^{\prime})}}$ is zero for $m\neq m^{\prime}$ Mizoguchi2019 , we find that only the terms with $m=m^{\prime}$ have a finite contribution in the summation of (LABEL:eq:edgeSSF_2). Then, we obtain

[TABLE]

where we have set $H|0\rangle=0$ for simplicity, and have used the translational symmetry along the vertical direction to obtain the second line.

To proceed, we employ a spectral representation:

[TABLE]

where $|\lambda\rangle$ denotes the eigenstate of $H$ that satisfies $H|\lambda\rangle=E_{\lambda}|\lambda\rangle$. To calculate the expectation value $\langle\lambda|c_{(1,1)}|0\rangle$ in Eq. (LABEL:eq:edgeSSF_4), we expand $c_{(1,1)}$ in terms of $\gamma^{\rm L}(k_{y})$. To do this, let us first recall that

[TABLE]

We expand $\alpha_{1,k_{y}}$ as

[TABLE]

Note that the anti-commutator $\{\alpha_{1,k_{y}},[\gamma^{\rm L}(k_{y})]^{\dagger}\}$ is equal to $W^{\ast}_{k_{y}}$, since $\{\gamma^{\rm L}(k_{y}),[\gamma^{\rm L}(k_{y})]^{\dagger}\}=1$. On the other hand, from Eq. (42), we have $\{\gamma^{\rm L}(k_{y}),\alpha^{\dagger}_{1,k_{y}}\}=\left[\{\alpha_{1,k_{y}},[\gamma^{\rm L}(k_{y})]^{\dagger}\}\right]^{\ast}=u_{1,k_{y}}$, since $\{\alpha_{1,k_{y}},\alpha_{1,k_{y}}^{\dagger}\}=1$. Combining these, we obtain $W_{k_{y}}=u_{1,k_{y}}^{\ast}$. Therefore, we have

[TABLE]

The rest except for the first two terms in the summand in the right-hand side of (97) do not contribute to the edge susceptibility because the frequency $\omega$ is assumed to be smaller than the bulk band gap.

Using Eq. (97) and noting $\gamma^{\rm L}(k_{y})|0\rangle=0$ $\{[\gamma^{\rm L}(k_{y})]^{\dagger}|0\rangle=0\}$ for $\varepsilon^{\rm L}(k_{y})>0$ [$\varepsilon^{\rm L}(k_{y})<0$], we have

[TABLE]

The matrix element $\langle 0|\gamma^{\rm L}(k_{y})|\lambda\rangle$ becomes finite when $|\lambda\rangle=\left[\gamma^{\rm L}(k_{y})\right]^{\dagger}|0\rangle$ and $\varepsilon^{\rm L}(k_{y})>0$, or $|\lambda\rangle=\gamma^{\rm L}(k_{y})|0\rangle$ and $\varepsilon^{\rm L}(k_{y})<0$. Therefore, we have

[TABLE]

where we have introduced $\tilde{k}_{\omega}$ such that $E^{\rm L}(\tilde{k}_{\omega})=\omega$.

From Eq. (99), we find that $\chi_{\rm edge}(\omega)$ is proportional to the inverse of the group velocity, meaning that the sign of the $\chi_{\rm edge}(\omega)$ is determined by that of the group velocity. Hence, by measuring the $\omega$ dependence of $\chi_{\rm edge}(\omega)$, one can determine whether the edge mode is chiral or not. To be concrete, if the edge mode is chiral, then the sign is always positive or negative, whereas if it is not chiral, then the sign change occurs on varying $\omega$.

VI Summary and discussion

In summary, we have treated the propagation of the edge magnetization due to the chiral edge modes and the edge dynamical susceptibility in the Kitaev honeycomb model with the three-spin interaction and the Kitaev triangle-honeycomb model. These quantities have been shown to refelect the chirality of the edge modes, whose emergence is a consequence of the topological nature of these models, i.e., the nontrivial Chern number and the bulk-edge correspondence. By relying on the exact solvability of these models, we have first mapped them to the free fermion models, and then applied the transfer matrix method. Further, by using the resulting edge wave functions and their dispersion relations, we have calculated (i) the time evolution of the edge magnetization and (ii) the edge dynamical susceptibility. For the former (i), we have demonstrated that the wave packet of the edge magnetization due to the chiral edge mode which is excited by an external perturbation indeed moves in only one of the two directions of the edge with the group velocity. For the latter (ii), we have shown that the edge dynamical susceptibility is proportional to the inverse of the group velocity with the positive weight. This property enables us to determine whether an edge mode is chiral or not. In fact, if an edge mode is chiral, then the sign of the edge dynamical susceptibility does not change through the whole range of the frequency; otherwise, it takes both of plus and minus signs.

Our method is also applicable to generic topological insulators and superconductors, whose concrete examples are the Haldane model Haldane1988 and the spinless superconductor Sato2016 , which were already treated in the text. In particular, the oriented propagation of the excess charge can be expected to be experimentally detected in a chiral superconductor.

Before closing this paper, we address the experimental observation of the oriented propagation of the edge magnetization. As we have seen, the propagation can be expected to occur by exciting the chiral edge spinons by shining a laser pulse to the Kitaev magnets under the magnetic field in [111] direction. Then, we expect that the flow of the magnetization can be detected by using the sophisticated techniques developed in the field of spintronics Bakun1984 ; Saitoh2006 ; Kimura2007 ; Maekawa2013 . In particular, the propagation of the edge magnetization can be expected to be stable against disorder that is inevitable in experimental situations because the chiral edge modes are known to be robust against perturbations Schulz-Baldes2000 ; Kellendonk2002 ; Elbau2002 .

Acknowledgements.

We wish to thank H. Katsura and M. Udagawa for fruitful discussion. T. M. is supported by the JSPS KAKENHI (Grant No. JP17H06138), MEXT, Japan.

Appendix A Transfer matrix method for the edge modes

In this Appendix, we present the details of the transfer matrix method by which we obtain the dispersion relation of the left edge mode in Sec. III.

To begin with, we recall that the BdG equation of Eq. (14) is expressed by using the transfer matrix $\hat{T}(k_{y})$ as Eq. (59). Let us consider the eigenvalue problem for $\hat{T}(k_{y})$. In the following, we abbreviate generally a function of $f(k_{y})$ of $k_{y}$ to $f$ by dropping the argument $k_{y}$.

From (59) and (63), the eigenvalue equation for $\hat{T}=\hat{T}(k_{y})$ can be written as

[TABLE]

with an eigenvalue $\lambda$. Although the size of $\hat{A}$ and $\hat{B}$ is $2\times 2$ in the present problem, the formalism of the transfer matrix is applicable to the case with any size of $n\times n$ ($n=1,2,\cdots$). Therefore, in what follows, we consider the case for generic $n$, whose eigenvalue equation can be obtained by replacing $\hat{I}_{2}$ with $\hat{I}_{n}$ in Eq. (108). We will get back to the case of $n=2$ after we show the generic formulation.

From Eq. (108), one has

[TABLE]

and

[TABLE]

When $\lambda\neq 0$, this can be written as

[TABLE]

For a nontrivial solution $\bm{\varphi}_{\ell-1}$, the two parameters, $\lambda$ and $\varepsilon^{\mathrm{L}}$, must satisfy

[TABLE]

Further, using the fact that if the matrix $M$ satisfies ${\rm det}M=0$, then ${\rm det}M^{\dagger}=0$, we have

[TABLE]

where we have used that $\hat{B}$ is Hermitian, and that $\varepsilon^{\mathrm{L}}$ is real. This implies the following: If $\lambda$ is a solution of the above equation for a given real $\varepsilon^{\mathrm{L}}$, then $1/\lambda^{\ast}$ is a solution, too Molinari1997 . Thus the equation has a set of the solutions,

[TABLE]

Here, we assume that the eigenvalues, $\lambda_{j}$, satisfy $0<|\lambda_{j}|<1$ for $\ell=1,2,\ldots,n$. In the present approach, we need this assumption in order to obtain a decay solution with the open boundary condition at the edge. For the corresponding $n$-component eigenvector which satisfies (110) with the eigenvalue $\lambda_{j}$, we write $\bm{\chi}^{(j)}$.

Then, the $2n$-component vector,

[TABLE]

is an eigenvector of the transfer matrix.

In order to find the solution which satisfies the open boundary condition at the edge, we set

[TABLE]

The coefficients $c_{\ell}$ are determined by

[TABLE]

because of the open boundary condition $\bm{\varphi}_{0}=0$. This implies

[TABLE]

Since all the eigenvalues $\lambda_{j}$ and all the vectors $\bm{\chi}^{(j)}$ are determined by the energy $\varepsilon^{\mathrm{L}}$, this equation determines the energy spectrum $\varepsilon^{\mathrm{L}}$ of the edge mode. But it is very difficult to determine the energy $\varepsilon^{\mathrm{L}}$ by the above equation (124) for general $n$.

Now, let us consider the case of $n=2$, which corresponds to our original problem for the Kitaev model. In this case, one clearly has

[TABLE]

This implies that the two vectors are the same, i.e, $\bm{\chi}^{(1)}={\rm(const.)}\bm{\chi}^{(2)}$. In addition to this, if $\lambda_{1}=\lambda_{2}$, then $\varphi_{\ell}=0$ for all $\ell$ from (122) and (123). Therefore, in this case, the single vector $\bm{\chi}^{(1)}$ has two different eigenvalues, $\lambda_{1}$ and $\lambda_{2}$, i.e., $\lambda_{1}\neq\lambda_{2}$. Then, the corresponding equations for (110) become

[TABLE]

and

[TABLE]

Subtracting (127) from (126) and using $\bm{\varphi}^{2}=\bm{\varphi}^{1}$, one has

[TABLE]

This implies that $\bm{\varphi}^{1}$ is the eigenvector of the matrix $\hat{A}^{-1}(\hat{B}-\varepsilon^{\mathrm{L}})$. Then, clearly, $\bm{\varphi}^{1}$ is the eigenvector of $\hat{A}^{-1}\hat{A}^{\dagger}$ as well.

Since the matrix $\hat{A}^{-1}\hat{A}^{\dagger}$ is independent of the energy $\varepsilon^{\mathrm{L}}$, one can obtain an eigenvector $u$ of $\hat{A}^{-1}\hat{A}^{\dagger}$ with the eigenvalue $\bar{\mu}$ which satisfies the decay condition. Then, we can choose the energy eigenvalue $\varepsilon^{\mathrm{L}}$ so that $u$ becomes the eigenvector of $\hat{A}^{-1}(\hat{B}-\varepsilon^{\mathrm{L}})$. More precisely, $\varepsilon^{\mathrm{L}}$ can be chosen so that the ratios of the first and the second components of the vectors in both sides of the eigenvalue equation $\hat{A}^{-1}(\hat{B}-\varepsilon^{\mathrm{L}})u=\bar{\eta}u$ with an eigenvalue $\bar{\eta}$ coincide with each other. Indeed, the dispersion relation of (64) is determined in this way. In consequence, we have

[TABLE]

with the eigenvalue $\bar{\eta}$. Combining this, $\hat{A}^{-1}\hat{A}^{\dagger}u=\bar{\mu}u$ and Eq. (110) with $\bm{\varphi}_{\ell-1}=u$, we have

[TABLE]

The two eigenvalues of the transfer matrix are given by

[TABLE]

Consequently, the wave function of the edge mode is given by

[TABLE]

which satisfies the boundary conditions when

[TABLE]

The condition (133) is the criterion of the existence of the edge mode at given $k_{y}$ (see the area shaded in yellow in Fig. 3).

Appendix B Chiral edge mode in the Kitaev triangle-honeycomb model

In this Appendix, we explain how to construct the solution for the edge modes of the Kitaev triangle-honeycomb model. To begin with, we expand the $i$-th edge mode as

[TABLE]

where the coefficients $u^{(i)}_{\ell,k_{y},\eta}$ and $v^{(i)}_{\ell,k_{y},\eta}$ satisfy the BdG equation:

[TABLE]

for $\ell=1,\cdots L_{x}$. To solve Eq. (135), we define

[TABLE]

and

[TABLE]

for $s=A,B,C$. Then, among six equations of (135), the four equations contain only the variables at the site $\ell$. These are given by

[TABLE]

The rest of the two equations are

[TABLE]

We write

[TABLE]

[TABLE]

and

[TABLE]

Then, the equation corresponding to Eq. (46) is given by

[TABLE]

Although the matrix $\hat{A}(k_{y})$ is not invertible Dwivedi2016 , we can find the left edge solutions which satisfy $\xi_{\ell,k_{y},C}^{(i)}=0$ for all $\ell$. Actually, when choosing $\xi_{\ell,k_{y},C}^{(i)}=0$, the rest of the five components satisfy

[TABLE]

Since this eigenvalue equation must have a nontrivial solution, the determinant of the matrix must be vanishing. As a result, the energy eigenvalue $\varepsilon^{L,(i)}$ can be determined by

[TABLE]

Since the above matrix is hermitian, this algebraic equation has five real solutions which are already found in Fig. 4 as edge modes. Further, the nontrivial solution of (170) is determined except for the normalization. Combining this with the equation for the sixth component of (159), we obtain the recursion equation $\zeta_{\ell+1,k_{y},A}^{(i)}=\lambda(k_{y})\zeta_{\ell,k_{y},A}^{(i)}$ with

[TABLE]

Clearly, the left-edge solutions must satisfy the decay condition $|\lambda(k_{y})|<1$ which determines the permissible region of $k_{y}$ for given $J$ and $J^{\prime}$. Figure 7 shows the region for the edge modes across the zero energy. As seen in Fig. 4, the permissible wave number $k_{y}$ for the edge modes is restricted to the neighborhood of $\pi$ in the topological phase, whereas, in the nontopological phase, the permissible wave number can take any value in the whole range. Similarly, the right edge solutions satisfy $\zeta_{\ell,k_{y},A}^{(i)}=0$ for all $\ell$, and we can obtain them in the same way. Although we cannot solve the algebraic equation (171) analytically, we believe that the useful information about the chiral edge modes has been able to be obtained.