Quantized thermal Hall conductance from edge current calculations in lattice models

TL;DR

This paper introduces a new method to calculate the quantized thermal Hall conductance in lattice models by focusing on edge currents, simplifying the process compared to traditional linear-response approaches.

Contribution

It proposes an edge current-based approach for computing thermal Hall conductance in chiral topological models, overcoming limitations of the Kubo formula in microscopic simulations.

Findings

Edge current calculation effectively captures quantized thermal Hall conductance.

The method is validated on chiral p-wave superconductors and Hofstadter models.

Finite-size effects are analyzed, recommending long strip geometries for simulations.

Abstract

The quantized thermal Hall effect is an important probe for detecting chiral topological order and revealing the nature of chiral gapless edge states. The standard Kubo formula approach for the thermal Hall conductance $\kappa_{xy}$ based on the linear-response theory faces difficulties in practical application due to the lack of a reliable numerical method for calculating dynamical quantities in microscopic models at finite temperature. In this work, we propose an approach for calculating $\kappa_{xy}$ in two-dimensional lattice models displaying chiral topological order. Our approach targets at the edge current localized at the boundary which involves only thermal averages of local operators in equilibrium, thus drastically lowering the barrier for the calculation of $\kappa_{xy}$. We use the chiral $p$-wave superconductor (with and without disorder) and the Hofstadter model as benchmark examples to illustrate several sources of finite-size effects, and we suggest the infinite (or sufficiently long) strip as the best geometry for carrying out numerical simulations.

Equations

$${\kappa }_{xy}=\frac{\mathrm{d}{j}_{\mathrm{edge}}(T)}{\mathrm{d}T}\text{ for }T\ll {\Delta }_{\mathrm{bulk}}.$$

$$j_{\mathrm{edge}}(T)=\frac{\pi c_{-}}{12}{T}^2\text{ for }T\gg v/N_y,$$

$$\frac{{\kappa }_{xy}}{T}=\frac{\pi c_{-}}{6},$$

$$\frac{\mathrm{d}{\hat{H}}_n}{\mathrm{d}t}=-\mathrm{i}[{\hat{H}}_n,\hat{H}]=-\sum _m\mathrm{i}[{\hat{H}}_n,{\hat{H}}_m],$$

$${\hat{J}}_{m\to n}=\mathrm{i}[{\hat{H}}_m,{\hat{H}}_n].$$

$$\hat{j}=\frac{1}{N_y}\sum _n{\hat{J}}_{n\to n+1}=\frac{\mathrm{i}}{N_y}\sum _n[{\hat{H}}_n,{\hat{H}}_{n+1}],$$

$${\hat{J}}_{n\to n+1}$$

$${\hat{j}}^{(r)}=\frac{1}{N_y}\sum _n{\hat{J}}_{n\to n+1}^{(r)}=\frac{\mathrm{i}}{N_y}\sum _i[{\hat{H}}_n^{(r)},{\hat{H}}_{n+1}^{(r)}],$$

$${\kappa }_{xy}=\frac{\mathrm{d}{j}_{\mathrm{edge}}}{\mathrm{d}T}=\frac{\mathrm{d}⟨{\hat{j}}^{(r)}{⟩}_T}{\mathrm{d}T},$$

$${\kappa }_{xy}$$

$$H$$

$$H=-\sum _{m,n}\left(f_{m,n+1}^{†}{f}_{m,n}+e^{\mathrm{i}2\pi \Phi n}{f}_{m+1,n}^{†}{f}_{m,n}+h.c.\right)-\mu \sum _{m,n}{f}_{m,n}^{†}{f}_{m,n},$$

$$r=q{s}_r+p{t}_r,|t_r|<q/2,s_r,t_r\in \mathbb{Z},$$

$$j_Q\equiv ⟨\mathcal{P}{⟩}_T=\frac{v^2}{2\pi }⟨\mathcal{T}-\overline{\mathcal{T}}{⟩}_T,$$

$$j_Q$$

$$(L_{-1}{)}_{\mathrm{cyl}}=\frac{2\pi }{N_y}\left[(L_0{)}_{\mathrm{plane}}-\frac{c}{24}\right].$$

$$j_Q$$

$$j_Q=-\frac{2\pi v^2}{N_y^2}{\partial }_{\xi }\ln {\chi }_a,$$

$$\ln {\chi }_a=\ln \left(\sum _b{S}_{ab}{\chi }_b\left(-\frac{1}{\tau }\right)\right)=\ln (S_{aI}{\mathfrak{q}}^{-c_{-}/24}+O({\mathfrak{q}}^{1-c_{-}/24})),$$

$$\ln {\chi }_a=\ln (S_{aI}{\mathfrak{q}}^{-c_{-}/24})=\frac{c_{-}{\pi }^2}{6\xi }+\ln S_{aI},$$

$$j_Q(T)=\frac{\pi c_{-}}{12}{T}^2.$$

$$j_Q(T=0)=\frac{2\pi v^2}{N_y^2}\left(h_a-\frac{c_{-}}{24}\right),$$

$$\tilde{M}=OMO^{T}=\bigoplus_{l=1}^{N}\left(\begin{array}[]{cc}0&\lambda_{l}\\ -\lambda_{l}&0\end{array}\right),$$

$$\exp (-\beta \hat{H})=\exp \left(-\frac{\mathrm{i}\beta }{4}{\tilde{\gamma }}^T\tilde{M}\tilde{\gamma }\right)=\prod _{l=1}^N\exp \left(-\frac{\mathrm{i}\beta }{2}{\lambda }_l{\tilde{\gamma }}_{l,1}{\tilde{\gamma }}_{l,2}\right),$$

$$\exp(-\beta\hat{H})=\prod_{l=1}^{N}\exp\left[\beta\lambda_{l}\left(a_{l}^{\dagger}a_{l}-\frac{1}{2}\right)\right]\dot{=}\bigotimes_{l=1}^{N}\left(\begin{array}[]{cc}\mathrm{e}^{-\frac{1}{2}\beta\lambda_{l}}&\\ &\mathrm{e}^{\frac{1}{2}\beta\lambda_{l}}\end{array}\right),$$

$$Z=\mathrm{Tr}({\mathrm{e}}^{-\beta \hat{H}})=\prod _{l=1}^{N}2\cosh \left(\frac{1}{2}\beta {\lambda }_l\right).$$

$$K_{ab}=\frac{\mathrm{i}}{2Z}\mathrm{Tr}({\mathrm{e}}^{-\beta \hat{H}}[{\gamma }_a,{\gamma }_b]).$$

$$K=O^T\tilde{K}O,$$

$${\tilde{K}}_{(l,1),(l,2)}$$

$$=$$

$$⟨\hat{Q}⟩=\frac{\mathrm{i}}{4Z}\sum _{a,b}{M}_{ab}^{(Q)}\mathrm{Tr}({\mathrm{e}}^{-\beta \hat{H}}{c}_a{c}_b)=-\frac{1}{4}\mathrm{Tr}(M^{(Q)}K),$$

Full text

††preprint: APS/123-QED

I Introduction

In recent years there has been a prodigious interest in non-Abelian anyons emerging in condensed matter systems, which are of great significance in the contexts of both fundamental science and applied subjects Nayak et al. (2008). These non-Abelian anyons, which appear as quasiparticles in systems with topological order, allow braiding manipulations—i.e., exchanging these quasiparticles drives the system into a distinct state—and upon these braiding processes lays the foundation for fault-tolerant quantum computation Kitaev (2003); Nayak et al. (2008). Originating from the fractional quantum Hall system with filling factor $\nu=5/2$ Moore and Read (1991); Read and Green (2000), the search for non-Abelian anyons also proceeds to other quantum systems, such as superconductor-semiconductor heterostructures Lutchyn et al. (2018) and chiral spin liquids Kitaev (2006). Despite extensive theoretical and experimental effort, a completely confirmed experimental observation of non-Abelian anyons still remains elusive.

In realistic materials, the quantized thermal Hall effect could serve as a strong indicator of the existence of non-Abelian anyons. For chiral topological systems, there exist gapless edge modes at the boundary that, at the low-energy limit, are generally described by chiral conformal field theories (CFTs) in $(1+1)$ dimensions. The energy current carried by the chiral edge mode has a universal form Cappelli et al. (2002) $I(T)=(\pi k_{B}c_{-}/12\hbar)T^{2}$, where $k_{B}$, $\hbar$, and $T$ denote the Boltzmann constant, the reduced Planck constant, and the temperature, respectively. Here $c_{-}$ is the chiral central charge, which characterizes the edge CFT and reflects the topological properties in the bulk, can be directly measured in a thermal Hall experiment, i.e., $\hbar\kappa_{xy}/k_{B}T=\pi c_{-}/6$ Kane and Fisher (1997); Cappelli et al. (2002), where $\kappa_{xy}$ is the thermal Hall conductance. Specifically, if the edge mode includes a Majorana branch, the Majorana edge mode would contribute a chiral central charge $c_{-}=1/2$. In this regard, a half-integer quantized thermal Hall conductance provides strong evidence for the existence of the Majorana edge mode and non-Abelian anyon excitations in the bulk. Very recently, the anticipated half-integer thermal Hall conductance is observed in the $\nu=5/2$ quantum Hall system Banerjee et al. (2018), as well as the field-induced disordered state in the Kitaev material $\alpha$-RuCl3 Kasahara et al. (2018).

In spite of the remarkable experimental progress, methods for a direct calculation of the thermal Hall conductance from microscopic Hamiltonians are circumscribed, which poses an obstacle to the theoretical understanding of the experimental observations. For example, despite the observation of the half-integer quantized thermal Hall conductance in $\alpha$-RuCl3, the corresponding theoretical interpretation is still under debate. The microscopic model describing the magnetic properties of $\alpha$-RuCl3 remains obscure (see, e.g., Trebst ; Janssen and Vojta (2019)), and it is largely unclear whether the experimentally observed field-induced disordered state is adiabatically connected to the non-Abelian spin liquid with dominant Kitaev interaction perturbed by a magnetic field. One of the important reasons for these obscurities is attributed to the lack of a suitable theoretical toolbox for handling these problems—the existing numerical calculations of the thermal Hall conductance rely on the usage of the Kubo formula Kubo (1957); Luttinger (1964), which requires an evaluation of dynamical quantities, which are commonly unattainable in interacting systems. These calculations are limited, therefore, to the cases in which a mapping of the system to a quadratic fermion model is possible Nasu et al. (2017); Metavitsiadis et al. (2017), or to small-sized systems that are within the capability of exact diagonalization. Moreover, the numerical application of the Kubo formula in thermal Hall problems is plagued by the subtle magnetization correction that originates from circulating currents in the system Cooper et al. (1997); Qin et al. (2011); Nomura et al. (2012); Stone (2012); Sumiyoshi and Fujimoto (2013); Vinkler-Aviv (2019), and the existing analytic solution Qin et al. (2011); Tatara (2015) to this correction is rather difficult to use in numerics, especially in the presence of interactions. For these reasons, a direct calculation of the thermal Hall conductance from interacting systems remains elusive.

In this work, we propose a method to directly address the edge modes of chiral topological systems and compute the quantized thermal Hall conductance from the edge currents, as motivated by Kitaev’s work Kitaev (2006). In contrast to previous work using the Kubo formula, our edge current approach only involves thermodynamic calculations of the locally defined operators, which are of much less computational cost. Furthermore, this approach is rather straightforward and thus is free from impediments such as energy magnetization corrections. In these respects, it has the potential to be applied to interacting systems with prevailing numerical algorithms.

As a preliminary study, we investigate the performance of the edge current approach using several non-interacting fermion systems as benchmark examples. By analyzing the finite-size effects in this approach, including the discretization of the edge spectrum and the overlap between edge modes, we suggest the finite-width long strip as an optimal choice of the lattice geometry when applying it to interacting systems, which might be within the scope of the state-of-the-art tensor-network algorithms. The edge current approach using this suggested lattice geometry is further tested in the Hofstadter model and the disordered $p$-wave superconductor.

This paper is organized as follows. In Sec. II, we present our edge current approach to thermal Hall conductance, and we describe the applicability of this approach. In Sec. III, we show preliminary numerical results in the chiral $p$-wave superconductor, discuss two major finite-size effects in this approach, and suggest an optimal lattice geometry for this approach. In Sec. IV, we apply our method to the Hofstadter model and the disordered chiral $p$-wave superconductor. Finally, in Sec. V we summarize our work and give an outlook. In Appendix A, we review the renowned universal $T^{2}$-dependence of the edge energy current and also derive the analytical form of the ground-state contribution to the energy current. In Appendix B, we provide details on how our approach is performed in free-fermionic systems.

II Method

We begin by introducing our edge current approach for the thermal Hall conductance. We put the system upon which we focused in our work—the gapped chiral topological system—on a cylinder, in order to introduce system boundaries and get access to the gapless edge modes. When the temperature of the system is much smaller than the bulk gap $\Delta_{\mathrm{bulk}}$, the energy current carried by the chiral edge mode, from which the quantized thermal Hall effect originates 111In $U(1)$-conserved systems, the heat current $J_{Q}$ and the energy current $J_{E}$ are conceptually different, and they are related to each other by $J_{Q}=J_{E}-\mu J_{N}$, where $J_{N}$ and $\mu$ represent the $U(1)$ current and the chemical potential of the system, respectively. In these systems, however, it is always legitimate to shift the chemical potential to zero by absorbing $-\mu\hat{N}$ into the Hamiltonian $\hat{H}$, i.e., $\hat{H}\rightarrow\hat{H}-\mu\hat{N}$, where $\hat{N}$ is the particle number operator. Along these lines, we hereby identify the energy current and the thermal current in this paper., can naturally be used to compute the thermal Hall conductance [see Fig. 1(a)]

[TABLE]

Here $j_{\mathrm{edge}}(T)$ denotes the edge current density in the thermal equilibrium state at the temperature $T$, which is predicted to be of the universal form Cappelli et al. (2002)

[TABLE]

where we have set $k_{B}=\hbar=1$. Here $c_{-}$, $v$, and $N_{y}$ denote the chiral central charge, the velocity of the edge mode and the circumference of the cylinder, respectively. A straightforward combination of Eqs. (1) and (2) gives the quantized thermal Hall conductance

[TABLE]

where $T$ satisfies $T\ll\Delta_{\mathrm{bulk}}$ and $T\gg v/N_{y}$.

Despite the simplicity of this scheme, it is far from obvious how to obtain $j_{\mathrm{edge}}(T)$ in numerical calculations. For this purpose, we first introduce the definition of local energy current in the lattice system. For a lattice model constituted by local Hamiltonian terms $\hat{H}=\sum_{m}\hat{H}_{m}$, according to the Heisenberg equation

[TABLE]

one may define the local current flowing from $\hat{H}_{m}$ to $\hat{H}_{n}$ as

[TABLE]

In practice, however, there may exist multiple ways to define such local energy currents in a lattice system. First, the choice of partitioning the total Hamiltonian into local Hamiltonian terms is not unique. Moreover, given a choice of local Hamiltonian terms, Eq. (5) may be only one of the possible definitions—in two dimensions, there may exist more than one definitions of the local currents that satisfy the energy conservationQin et al. (2011); Metavitsiadis et al. (2017). These ambiguities in the definition of the local energy current may call into question the reliability of calculations based on this definition, especially when the Hamiltonian is rather complicated—for example, there may exist multiple-site interactions in the system.

Nevertheless, given a physically correct approach, the calculation results of the thermal Hall conductance—a measurable physical quantity—should not depend on the choice of conventions. Indeed, in the following, we introduce a suitable setup, under which the effect of the convention differences is negligible.

To start with, we view the two-dimensional system as a quasi-one-dimensional system by slicing the two-dimensional lattice into layers, as schematically depicted in Fig. 1(b). We presume that only local Hamiltonians at the nearest-neighboring layers have non-vanishing commutators, i.e., $[\hat{H}_{m},\hat{H}_{n}]\neq 0$ only when $|m-n|=1$, where $\hat{H}_{m}$ denotes the local Hamiltonian at the $m$-th layer. This condition can be achieved with a proper choice of Hamiltonian term at each layer. In this simple one-dimensional system, the energy current density along the circumferential direction is apparently given by

[TABLE]

where we have set the length of each layer to be one unit length in the circumferential direction [see Fig. 1(b)].

Next, in order to investigate the edge currents of the system, we spatially divide the layers into two parts, $\hat{H}_{n}=\hat{H}_{n}^{(l)}+\hat{H}_{n}^{(r)}$, where $\hat{H}_{n}^{(l/r)}$ is on the left/right side [see Fig. 1(c)]. The current operator $\hat{J}_{n\rightarrow n+1}$ becomes

[TABLE]

Since the edge modes are exponentially localized at the boundaries, it is legitimate to concentrate on $\hat{J}_{n\rightarrow n+1}^{(r)}\equiv\mathrm{i}[\hat{H}_{n}^{(r)},\hat{H}^{(r)}_{n+1}]$ which covers the edge current contributions at the right boundary [see Fig. 1(c)]. In this regard, we introduce the energy current density operator in the right part

[TABLE]

whose thermal average $\langle\hat{j}^{(r)}\rangle_{T}$ contains all the edge current contributions at the right boundary. Here $\langle\hat{j}^{(r)}\rangle_{T}$ is given by $\langle\hat{j}^{(r)}\rangle_{T}=(1/Z)\mathrm{Tr}(\mathrm{e}^{-\beta\hat{H}}\hat{j}^{(r)})$, where $Z$ represents the partition function, and $\beta=1/T$. It is apparent that $\langle\hat{j}^{(r)}\rangle_{T}$ is not solely contributed by the edge current—it also contains contributions from bulk states. These bulk-state contributions may come from the circulating bulk currents [see Fig. 1(d)]. Furthermore, due to the spatial separation $\hat{H}_{n}=\hat{H}_{n}^{(l)}+\hat{H}_{n}^{(r)}$, $\langle\hat{j}^{(r)}\rangle_{T}$ may also contain bulk-state contributions at the dividing point, and for different ways of dividing the local Hamiltonians in the central region, these bulk contributions are different. Nevertheless, if we restrict the temperature to be much smaller than the bulk gap (as in Eq. (1)), the bulk states with energies higher than the bulk gap are suppressed—in this case, the bulk-state contributions contained in $\langle\hat{j}^{(r)}\rangle_{T}$ mostly come from the ground state and are hence nearly independent of the temperature. Along these lines, and combining Eq. (1), the thermal Hall conductance can be calculated as

[TABLE]

where we have restricted the temperature to $T\ll\Delta_{\mathrm{bulk}}$.

Equation (9) is the central result of our work. In numerical applications, one needs to calculate $\langle\hat{j}^{(r)}\rangle_{T}$ for several temperatures and make use of the numerical differentiation. An alternative recipe is to subtract the ground-state contribution in $\langle\hat{j}^{(r)}\rangle_{T}$, i.e., to make use of the $T^{2}$-dependence of the edge current and calculate $\kappa_{xy}$ as

[TABLE]

where $T\ll\Delta_{\mathrm{bulk}}$, and bulk-state contributions in $\langle\hat{j}^{(r)}\rangle_{T}$ and $\langle\hat{j}^{(r)}\rangle_{T=0}$ cancel each other out. However, we need to point out that Eq. (10) suffers from a finite-size correction that comes from the subtracted ground-state edge current $j_{\mathrm{edge}}(T=0)\sim O(1/N_{y}^{2})$. The analytical form of this finite-size correction is derived in Appendix A.

We emphasize that Eq. (9) only involves the thermal averages over locally defined quantities, and it requires much less numerical cost compared to the conventional approach based on the Kubo formula Luttinger (1964); Qin et al. (2011); Nomura et al. (2012); Sumiyoshi and Fujimoto (2013), which relies on the evaluation of the dynamics of the system. Furthermore, our approach is free from the subtle issue of energy magnetization correction Qin et al. (2011). As a trade-off, this approach is restricted to gapped chiral topological systems, and the temperature must be much lower than the bulk gap. In this regard, the edge current approach cannot be applied to the thermal Hall effect contributed by magnons Katsura et al. (2010); Onose et al. (2010); Matsumoto and Murakami (2011a, b) or phonons Strohm et al. (2005); Sheng et al. (2006), where the bulk states are gapless and the corresponding thermal Hall conductance is generally nonquantized.

III Discussions

Despite the simplicity of Eq. (9), this equation gives the anticipated quantized thermal Hall conductance only when the following conditions are satisfied. First, the temperature $T$ must satisfy $T\ll\Delta_{\mathrm{bulk}}$ and $T\gg v/N_{y}$. In addition, the length of the cylinder, $N_{x}$, should be large enough so that the overlap between edge modes is minimized. In these respects, there are restrictions on the choice of the system size. For non-interacting fermion systems, one can calculate Eq. (9) to reasonably large system sizes (see Appendix B for details of the calculation). Hence, in this section, using the $p+\mathrm{i}p$ superconductor on the square lattice as a platform, we numerically investigate the finite-size effects arising in the edge current approach and look for an optimal lattice geometry for this calculation. The Hamiltonian of the $p+\mathrm{i}p$ model is given by

[TABLE]

where $f_{m,n}$ is the fermion operator, with $m,n$ representing the site index. $t$, $\mu$, and $\Delta$ correspondingly represent the hopping parameter, the chemical potential, and the pairing potential. We set $t>0$ and $\Delta>0$. When $-4t<\mu<0$ and $0<\mu<4t$, the $p+\mathrm{i}p$ model is in the gapped topological phase (also known as the “weak pairing” phase), where the sign of $\mu$ determines the chirality, while for $\mu<-4t$ or $\mu>4t$, the system is in the trivial phase (also known as the “strong pairing” phase) Read and Green (2000). When the system is in the topological phase, a chiral Majorana edge mode exists at the boundary and leads to a chiral central charge $c_{-}=\pm 1/2$, where the sign is determined by the chirality of the edge mode. In the following, we concentrate on the topological phase corresponding to $0<\mu<4t$, and we demonstrate the performance of the edge current approach in the case of various system sizes.

III.1 Restrictions on $N_{y}$: discretization of the edge spectrum and bulk contributions

Let us first explore the requirements on the circumference of the system. As previously mentioned, Eq. (9) gives the anticipated quantized thermal Hall conductance only when the conditions $T\gg v/N_{y}$ and $T\ll\Delta_{\mathrm{bulk}}$ are satisfied. The first condition $T\gg v/N_{y}$ arises due to the discrete edge spectrum—the temperature should be much larger than the edge spectrum discretization $v/N_{y}$; otherwise, the edge spectrum appears as a few separate energy levels and the prediction from the continuous field theory will break down. The second condition $T\ll\Delta_{\mathrm{bulk}}$ is imposed to prevent bulk-state contributions. Apparently, to ensure the existence of a temperature regime that satisfies the two conditions simultaneously, we should choose a sufficiently large circumference $N_{y}$.

In the $p+\mathrm{i}p$ superconductor, whose Hamiltonian is given by Eq. (11), we numerically investigate the constraint on the circumference of the system. In Eq. (11), we fix the hopping parameter and the pairing potential as $t=\Delta=1.0$. We set the length of the system to be $N_{x}=32$, which is large enough to avoid the overlap between edge modes at opposite boundaries. The circumference of the system varies between $N_{y}=20$ and $N_{y}=200$, which leads to different edge spectrum discretizations. Meanwhile, we adjust the bulk gap by changing the chemical potential between $\mu=1.0$ and $\mu=0.5$.

Fig. 2 shows the numerical results of the calculated $\kappa_{xy}/T$ from the edge current approach. As can be observed in Fig. 2, if the system circumference is sufficiently large, the anticipated plateau appears. Apart from the anticipated plateau, deviations from the quantization appear in both the low-temperature and high-temperature regimes. By comparing results calculated with different system circumferences, one can verify that the low-temperature deviations mainly originate from the edge spectrum discretizations, since these deviations are reduced in the cases of larger circumferences. On the other hand, by comparing Fig. 2(a) and (b), we find that for the case with a smaller bulk gap, the plateau shrinks evidently in the high-temperature regime, which is in agreement with our analysis—the edge current approach is affected by bulk-state contributions at high temperatures.

The numerical results in Fig. 2 also suggest that the required circumference of the system is very large, which is of the scale $N_{y}\sim 10^{2}$ (or probably even higher, if the bulk gap is smaller).

As a comparison, we show the numerical result of $\kappa_{xy}/T$ in the gapped trivial phase in Fig. 3, where $t=\Delta=1.0$, $\mu=4.1$, and the bulk gap $\Delta_{\mathrm{gap}}=0.1$. In Fig. 3, the $\kappa_{xy}$ result vanishes at low temperatures and becomes finite when the temperature is high enough. This observation is consistent with our previous analysis: Due to the existence of the bulk gap, when $T\ll\Delta_{\mathrm{bulk}}$, the bulk contributions to $\langle\hat{j}^{(r)}\rangle_{T}$ are nearly independent of $T$, and the result becomes zero due to the absence of edge modes. When the temperature is high, the bulk states are populated, which give rise to nonzero contributions to the numerical results of $\kappa_{xy}/T$. Unlike the edge current, these bulk-state contributions show no apparent dependence on $N_{y}$ for large enough $N_{y}$, which may be used as a criterion to distinguish the edge and bulk currents. This is in agreement with the results in Fig. 2, where $\kappa_{xy}/T$ in the high-temperature regime also becomes independent of $N_{y}$ for $N_{y}\geq 40$, since $\kappa_{xy}/T$ results at high temperatures are mainly contributed by bulk states. The $\kappa_{xy}/T$ data in the trivial phase also shows a specious plateau, which mostly lies at relatively high temperature and corresponds to a nonquantized value. This suggests that essential prior knowledge of the chiral central charge or energy scale of the system is necessary; otherwise, this specious plateau in the trivial phase may lead to fallacious conclusions.

III.2 Restrictions on $N_{x}$: overlap effect

The other restriction on the geometry of the lattice comes from the overlap between edge modes. When deriving Eq. (9), it is assumed that the edge modes are well localized at the boundaries. In practice, however, the edge modes usually have finite decaying length $l_{e}$. By dimensional analysis, one can infer that the decaying length $l_{e}$ scales as $l_{e}\propto v/\Delta_{\mathrm{bulk}}$, where $v$ is the velocity of the edge mode.

When the system is deep in the topological phase, the bulk gap is usually rather large and the edge modes are well localized. As the system approaches the phase transition point, the $\Delta_{\mathrm{bulk}}$ decreases and the edge modes delocalize into the bulk. In this case, if the system length $N_{x}$ is not sufficiently large, the edge modes on opposite boundaries would have a large overlap with each other, which breaks the validity of Eq. (9).

Using the $p+\mathrm{i}p$ superconductor given by Eq. (11), we numerically investigate the overlap effect. By setting $t=\Delta=1.0$, the velocity of the edge mode is fixed as $v=2.0$, and the $\Delta_{\mathrm{bulk}}$ varies linearly with $\mu$ [see Fig. 4(a)]. We adjust the decaying length of the edge mode by varying the chemical potential and thus the bulk gap $\Delta_{\mathrm{bulk}}$, and we test the performance of the edge current approach for different choices of $N_{x}$’s. For each $N_{x}$, the capability of the edge current approach is determined by a search for the minimal $\Delta_{\mathrm{bulk}}$ where the quantized plateau in $\kappa_{xy}/T$ can appear. The results are summarized in Fig. 4.

From Fig. 4, one can see that the required $N_{x}$ diverges as the bulk gap approaches zero, and when $N_{x}\leqslant 10$, the edge current approach can only capture the central region of the topological phase. This indicates that the edge current approach is capable of detecting the existence of a chiral topological phase, but it lacks the power to precisely determine topological phase transition points, where the bulk gap closes.

In the following, we concentrate on the cases when the overlap effect is not negligible but is of weak strength. We show that in this case, although the condition for Eq. (9) is violated, the edge current approach can still approximately give the quantized $\kappa_{xy}/T$.

The overlap between the edge mode induces a small gap $\Delta_{\mathrm{edge}}$ in the edge spectrum and drive the low-energy part of the edge spectrum apart from linearity. In this case, the low-temperature behavior of the numerical result for $\kappa_{xy}/T$ depends on the choice of system circumference $N_{y}$. For systems with finite circumference $N_{y}$, the edge spectrum is discrete with size $\sim v/N_{y}$. For small $N_{y}$ which satisfies $\Delta_{\mathrm{edge}}\ll v/N_{y}$, when the temperature is in the vicinity of $T\sim\Delta_{\mathrm{edge}}$, the higher-energy excitations (with energies satisfying $E\geq\Delta_{\mathrm{edge}}+v/N_{y}\gg T$) are far from populated and only the lowest-energy excitations at $\Delta_{\mathrm{edge}}$ contribute to the thermal Hall conductance. As a result, a peak in $\kappa_{xy}/T$ would arise in the vicinity of $T\sim\Delta_{\mathrm{edge}}$. This large low-temperature peak can be remedied by choosing a larger $N_{y}$, and when the condition $\Delta_{\mathrm{edge}}\ll v/N_{y}$ is no longer satisfied, the low-temperature peak would gradually disappear. However, due to the non-linearity of the low-energy part of the spectrum, there will still be a weak deviation from the quantization value in this temperature regime. Despite these finite-size effects, one may expect that the edge spectrum remains linear at relatively higher temperatures and gives rise to an approximately quantized $\kappa_{xy}/T$.

We numerically verify the preceding analysis using the $p+\mathrm{i}p$ model with weak overlap effects. The overlap effect is gradually enhanced by adjusting the chemical potential $\mu$ towards the phase transition point. We sequentially calculate $\kappa_{xy}/T$ for $\mu=0.9,0.8,0.7$. The length of the system is fixed as $N_{x}=10$, while the system circumference varies among $N_{y}=200,400,600$.

Fig. 5 shows the numerical results and the corresponding edge spectra. As indicated by the numerical results of $\kappa_{xy}/T$ in Figs. 5(a)(c)(e), the overlap effect arises for each choice of $\mu$ but has different strengths in different cases. This is in agreement with the edge spectra shown in Figs. 5(b)(d)(f), where the overlap effects—the small edge gap and the low-energy deviation from linearity—are gradually enhanced as the chemical potential approaches the phase transition point at $\mu=0$. What is more, for relatively small $N_{y}=200$, we observe the low-temperature peak in the $\kappa_{xy}/T$ result, which is attributed to the singular contribution of the lowest excitation state that arises in the temperature regime $T\sim\Delta_{\mathrm{edge}}$, as the condition $\Delta_{\mathrm{edge}}\ll v/N_{y}$ is satisfied. This finite-size artifact is remedied when we increase the system circumference to $N_{y}=600$, and we are left with a small hump in $\kappa_{xy}/T$ originating from the low-energy nonlinearity of the edge spectrum. Apart from the hump in the data, we can still observe a plateau in $\kappa_{xy}/T$ in the relatively higher temperature regime, which corresponds approximately with the quantized value $\pi/12$ [see Figs. 5(a)(b)]. The size of the hump and the length of the approximately quantized plateau depends on the strength of the overlap effect. These numerical observations are in agreement with the foregoing analysis.

III.3 Remarks

From the preceding numerical results, one can come to the following conclusions regarding the choice of the lattice geometry: (i) Due to the discretization effect and bulk-state contributions, the system circumference $N_{y}$ must be sufficiently large, which usually has to reach a scale of $\sim 10^{2}$ unit cells; otherwise, a quantized plateau in $\kappa_{xy}/T$ would not appear. For interacting systems, calculations with such a large $N_{y}$ are generally impractical, unless one uses an infinite algorithm to avoid a finite-size calculation. (ii) Due to the potential overlap effect between edge modes, there is a constraint on the system length $N_{x}$. However, from our numerical results, one may infer that a small $N_{x}\leqslant 10$ is sufficient to capture the central region of the topological phase, while the precise determination of the phase boundaries remains evasive with reachable choices of $N_{x}$. This suggests that an infinite tensor-network algorithm designed for an infinite strip would be suitable for the edge current calculation.

IV Applications

As preliminary applications and benchmark examples, we apply the edge current approach to two relatively challenging non-interacting systems, the Hofstadter model and the disordered $p+\mathrm{i}p$ superconductor. As suggested from the previous section, to mimic the lattice geometry of an infinitely long strip with a finite width, in our non-interacting calculations, we put the systems on a cylinder with a sufficiently large circumference $N_{y}$ and a relatively small length $N_{x}$.

IV.1 Hofstadter model

The Hofstadter model Hofstadter (1976) describes the behavior of charged particles in a magnetic field, in which we intend to test the performance of the edge current approach in the case of higher Chern numbers—in that case, there are more than one edge modes at each boundary. The Hamiltonian of the Hofstadter model is given by

[TABLE]

where we have used the Landau gauge. $\Phi=p/q$ is a rational number, which is determined by the magnetic flux through each plaquette. When $q>1$, the energy band splits into $q$ subbands, each of which has a nonzero Chern number. When the system has boundaries, there will be edge modes between the subbands, and the number of edge modes is determined by the total Chern number $C$ of the valence bands. Since in the Hofstadter model, the edge modes are complex (Dirac) fermions, each of these will contribute a chiral central charge $c_{-}=\pm 1$, where the sign depends on the chirality of the edge mode. In our present convention, a positive chiral central charge corresponds to the counterclockwise flowing edge current, which in turn corresponds to a negative Chern number, i.e., $c_{-}=-C$.

In the following, by choosing different values of the magnetic flux $\Phi$ and the chemical potential $\mu$, we obtain systems with different Chern numbers. In each choice of parameters, the total Chern number of the valence bands is calculated with the Diophantine equation Thouless et al. (1982),

[TABLE]

where the total Chern number of the $r$ valence bands is given by $t_{r}$. Using Eq. (9), we numerically calculate $\kappa_{xy}/T$, which is compared with the theoretical prediction given by Eqs. (3) and (13). The lattice geometry, as we have stated above, is a cylinder with sufficiently large circumference $N_{y}=400$ and relatively small length $N_{x}=12,16,20,24$.

The calculation results are shown in Fig. 6. For clarification, we also correspondingly plot the energy spectrum of the system, using the same parameters $\Phi$ and $\mu$. The energy spectrum is obtained with the lattice on a cylinder with $N_{x}=24$ and $N_{y}=60$.

First, as the simplest case, Fig. 6(a) shows $\kappa_{xy}/T$ calculated with $\mu=-1.5$ and $\Phi=1/3$, where the Chern number equals $C=1$, and there exists only one complex fermion mode at each boundary [see Fig. 6(b)]. In this case, we expect $c_{-}=-1$. As shown in Fig. 6(a), for each of the $N_{x}$’s that we have considered, a plateau corresponding to the expected quantization value $-\pi/6$ is observed. The low-temperature deviations can be attributed to the discretization of edge spectra, which may be in conjunction with the overlap effect, as discussed in Sec. III B. We also note that for $N_{x}=12$, the plateau in the $\kappa_{xy}/T$ data is only approximately quantized, which indicates the presence of weak overlap between edge modes.

Secondly, we show the results in the cases with the higher Chern number $C=2$, where two electron modes flow at each boundary. The chiral central charge becomes $c_{-}=-2$, and the quantization value of $\kappa_{xy}/T$ becomes $-\pi/3$. In Fig. 6(c), we show the results calculated with $\mu=-1.4,\Phi=1/7$. The corresponding energy spectrum is plotted in Fig. 6(d), which shows the existence of two electron edge modes. In the energy spectrum, one can observe that the bulk gap is smaller than the $C=1$ case, which leads to the inference that the finite-size effect would be more severe. Nevertheless, although only approximately quantized, the anticipated plateaus in the $\kappa_{xy}/T$ data are still observed, even for the smallest choice of the system length $N_{x}=12$.

Finally, in Fig. 6(e), we show results with much stronger overlap effects, where we have set $\mu=-0.6,\Phi=1/5$, and the corresponding chiral central charge is $c_{-}=-2$ [see Fig. 6(f)]. Due to the stronger overlap effect, in order to observe the quantized plateau in $\kappa_{xy}$, a larger system length $N_{x}$ is required. As one can observe from the numerical results, only in the data of $N_{x}=24$ does the expected plateau appear, and the calculations for smaller system lengths all break down in this case.

IV.2 Disordered $p+\mathrm{i}p$ model

As the second application, we apply the edge current approach in disordered systems, where translation symmetry is broken. Since the gapless edge modes are robust due to the protection of the bulk gap, the edge current approach is expected to remain effective in the presence of relatively weak disorder.

As a numerical test, we add on-site disorder to the clean $p+\mathrm{i}p$ model. The total Hamiltonian is expressed as $H_{\mathrm{tot}}=H+H^{\prime}$, where $H$ is the original $p+\mathrm{i}p$ Hamiltonian Eq. (11), and $H^{\prime}=-\sum_{m,n}W_{m,n}f_{m,n}^{\dagger}f_{m,n}$ represents the on-site disorder. Here $m,n$ are site indices, $W_{m,n}$’s are site-dependent random potentials that are uniformly distributed within the range $[-W,W]$, with $W$ representing the disorder strength. We fix the parameters of the original $p+\mathrm{i}p$ Hamiltonian to be $t=\Delta=1.0$ and $\mu=2.0$, in which case the bulk gap is $\Delta_{\mathrm{bulk}}=2.0$. The lattice is put on a cylinder with circumference $N_{y}=600$ and length $N_{x}=12$. We gradually increase the disorder strength $W$, and for each $W$, we calculate the thermal Hall conductance $\kappa_{xy}$ by averaging over 200 different disorder configuration. To minimize the computational cost, we make use of a supercell trick, i.e., the system is constituted by supercells of size $N_{x}\times M_{y}$, and the disorder configuration in each supercell is identical, as schematically depicted in Fig. 7. We compare results obtained with different sizes of supercells, $M_{y}=6,10,12,15$, and we expect the result of $\kappa_{xy}$ to converge if $M_{y}$ is sufficiently large. The numerical results are shown in Fig. 8.

First, in Fig. 8(a) we show the results calculated with disorder strength $W=2.5$, from which the anticipated quantized plateau is clearly observed. Furthermore, the error bars in the numerical results are almost negligible, which indicates that the disorder-induced fluctuation is very small. Although, as one may notice, the disorder strength is slightly over the bulk gap $\Delta_{\mathrm{bulk}}=2.0$, which indicates that the disorder may induce domains with different topological properties in the system, these domains are believed to be very small and cannot lead to considerable effects on our results. This result verifies the robustness of the edge current approach in the presence of weak disorder.

Secondly, Fig. 8(b) shows the result calculated with the disorder strength $W=3.5$. As can be observed, the disorder-induced fluctuations substantially increase and the quantized plateau is evidently diminished by the disorder. The reason for these fluctuations is twofold: (i) Since the disorder strength approaches the phase transition point, the domains in the system largely increase, while the system size remains the same. In this regard, the overall topological property of the finite-size system may become dependent on the different disorder configurations. (ii) In our calculation, we need to divide the system into two parts along some central line, and some disorder-induced domains may cross the central line, which gives rise to additional edge current contributions, as schematically shown in Fig. 7(b). When the sizes of the domains are sufficiently large, these additional edge currents can lead to substantial corrections to the numerical results. Moreover, since the domain distributions vary drastically with the disorder configurations, the additional edge contributions also give rise to large fluctuations in the results, even when the overall topological property of the system persists. These additional corrections may pose a limitation on the application of the edge current approach in the presence of strong disorder. Nevertheless, the disorder-induced fluctuations in the plateau regime are still visibly smaller [see Fig. 8(b)] and the plateau is still approximately quantized, which highlights the robustness of the edge current against disorder influences.

Finally, in Fig. 8(c) with the disorder strength $W=4.0$, the plateau totally disappears and large fluctuations dominate in the whole temperature regime. We infer that in this case the system is quite close to phase transition and the disorder-induced fluctuations predominate. We also expect that in this case, the effective overall bulk gap is very close to zero, which leads to delocalization of the edge modes. As a result, the edge current approach is no longer applicable to this very strong disorder strength.

Other than the quantized thermal Hall conductance, we can also obtain some other relevant information from the numerical result of $\kappa_{xy}/T$. In Fig. 9, we plot the $\kappa_{xy}/T$ with respect to the disorder strength in the low-temperature and high-temperature regime respectively. Fig. 9(a) shows that the deviation from the quantization at low temperatures is reduced by the weak disorder. Recall that this deviation originates from the violation of the condition $T\gg v/N_{y}$, thus we infer that the effective velocity of the edge mode is reduced due to the obstruction of the disorder. On the other hand, Fig. 9(b) shows that $\kappa_{xy}/T$ deviates more as one increases the disorder strength in the high-temperature regime. This numerical observation indicates that the overall effective bulk gap is reduced by the disorder.

V Conclusion and outlook

To conclude, we have proposed an approach to calculate the quantized thermal Hall conductance in gapped chiral topological systems—calculating the thermal Hall conductance by utilizing edge thermal currents directly. As we have seen, in contrast to the bulk Kubo formula paradigm, our edge current approach is physically transparent and hence free from subtleties such as energy magnetization corrections. More appealingly, the edge current approach only involves thermodynamic calculations of locally defined operators, which may get into the capability of prevailing numerical methods for interacting systems. As a preliminary study, using several non-interacting fermion systems as platforms, we investigate the major finite-size effects that would arise in our approach, based on which we suggest that an infinite (or sufficiently long) strip would serve as the best lattice geometry in these calculations, and the capability of this approach lies at the maximal width of the strip that is reachable in numerics. According to our numerical results, the edge current approach is competent when the system is deep inside the topological phase, but it lacks the capability to precisely determine phase boundaries. While it still requires more progress to apply the finite-temperature tensor-network algorithms Verstraete et al. (2004); White (2009); Stoudenmire and White (2010); Chen et al. (2018); Kshetrimayum et al. (2019) to interacting systems with the suggested infinite-strip lattice geometry, our work drastically reduces the computational barrier and serves as an important step towards that goal.

As an outlook, our edge current approach can serve as a numerical instrument for the detection of chiral topological phase in microscopic models. Specifically, a possible application lies in the thermal Hall effect in the Kitaev material $\alpha$-RuCl3 Kasahara et al. (2018). In addition to the recent theoretical interest in the interplay between acoustic phonons and the chiral Majorana edge mode Vinkler-Aviv and Rosch (2018); Ye et al. (2018), one may find several other experimental features intriguing—for example, the sensitivity of the field-induced spin liquid phase upon the direction and strength of an external magnetic field Baek et al. (2017); Zheng et al. (2017); Kasahara et al. (2018). From this perspective, it is desirable to directly compute the thermal Hall conductance in the several proposed microscopic models for $\alpha$-RuCl3 and compare with the experimental observations. Moreover, it would also be interesting to numerically examine and investigate the thermal Hall effect using the edge current approach in other theoretically proposed chiral spin liquid systems—e.g., Refs. Kalmeyer and Laughlin (1987); Nielsen et al. (2013); Bauer et al. (2014), as well as fractional quantum Hall systems Kane and Fisher (1997).

Acknowledgments

We thank W. Brenig, H. Jiang, J. R. Shi, Q. F. Sun, Y. Wan, and Y. H. Wu for helpful discussions. We acknowledge financial support from NBRPC under Grants No. 2015CB921102, No. 2017YFA0303301, and NSFC under Grants No. 11534001, No. 11504008 (W.T. and X.C.X.), Ministry of Science and Technology of China under Grant No. 2016YFA0302400 (L.W.) and the DFG through project A06 (H.H.T.) of SFB 1143 (project-id 247310070). W.T. acknowledges Institute of Theoretical Physics at TU Dresden for hospitality and Graduate School of PKU for financial support during part of this work.

Appendix A Derivation of the edge energy current of chiral topological systems

In this appendix, we derive the edge energy current of a chiral topological system using the partition function of the edge CFT. Along with the renowned result of Ref. Cappelli et al. (2002), which reveals the $T^{2}$-dependence of the edge energy current at finite temperature, we also derive the ground-state edge current, which gives rise to the finite-size correction $O(1/N_{y}^{2})$ in Eq. (10).

We consider the chiral topological system on a cylinder, as depicted in Fig. 1(b). At the two open boundaries of the cylinder there are the gapless chiral edge modes, which are described by the edge CFT. In the edge CFT, the energy current density is defined as the thermal average of the momentum density,

[TABLE]

where $\mathcal{T}(\omega)$ and $\overline{\mathcal{T}(\omega)}$ are the two chiral components of the energy-momentum tensor. Here $\omega=v\tau+\mathrm{i}x$ is a complex number, where $v$, $\tau$, and $x$ denote the velocity, the imaginary time, and the spatial coordinate, respectively. Due to the periodic boundary condition in the spatial dimension, $\omega$ is a coordinate on a spacetime cylinder. Since the thermal current density is translationally invariant, it is equal to the spatial average of the momentum density,

[TABLE]

where $L_{n}(\overline{L}_{n})$ ($n\in\mathbb{Z}$) is the Laurent mode of the energy-momentum tensor $\mathcal{T}(\omega)$($\overline{\mathcal{T}(\omega)}$), and $N_{y}$ denotes the circumference of the original system, which corresponds to the length of the edge theory [see Fig. 1(a)]. By $z=\exp(2\pi\omega/N_{y})$, we map the spacetime cylinder to a complex plane. Then

[TABLE]

In the following, we focus on the complex plane and the subscript “plane” will be dropped.

With regard to the chiral edge CFT, we only consider the simplest case, in which the anti-holomorphic part is completely absent—i.e., we assume the absence of counterflowing edge modes (the discussions below can be easily generalized to the case with counterflowing edge modes being present). The edge CFT partition function then only contains the holomorphic part, which is given by the character $\chi_{a}(\tau)=\mathrm{Tr}_{a}(q^{L_{0}-c_{-}/24})$, where $c_{-}$ represents the chiral central charge, $q=\mathrm{exp}(2\pi\mathrm{i}\tau)$, $\tau=\mathrm{i}v\beta/N_{y}$ (not to be confused with the “$\tau$” in $\omega=v\tau+\mathrm{i}x$), and $\beta=1/T$ represents the inversed temperature. Here $a$ labels the primary state of the conformal tower, which is associated with the quasiparticles in the bulk, and the trace is over the conformal tower of states Francesco et al. (1997). By writing $\tau=v(\gamma+\mathrm{i}\beta)/N_{y}$, where $\gamma\rightarrow 0^{+}$, we then have

[TABLE]

By introducing $\xi\equiv-2\pi\mathrm{i}\tau=(-2\pi\mathrm{i}v/N_{y})(\gamma+\mathrm{i}\beta)\xrightarrow{\gamma\rightarrow 0}2\pi v\beta/N_{y}$, Eq. (17) can be simplified as

[TABLE]

With Eq. (18), we can readily compute the thermal current density $j_{Q}$ in the two different limits $T\gg v/N_{y}$ and $T\rightarrow 0^{+}$, respectively.

When $T\gg v/N_{y}$, by making use of the modular transformation $\chi_{a}(\tau)=\sum_{b}S_{ab}\chi_{b}(-1/\tau)$, we get

[TABLE]

where $\mathfrak{q}\equiv\exp(-4\pi^{2}/\xi)$, and $S$ represents the $S$-modular matrix. Under the condition $N_{y}\gg v\beta$, $\xi=2\pi v\beta/N_{y}\rightarrow 0^{+}$, then $\mathfrak{q}\rightarrow 0^{+}$. Since $\mathfrak{q}$ exponentially vanishes as $\xi$ approaches zero, we infer there are no polynomial corrections in the expression of $\ln\chi_{a}$. Therefore,

[TABLE]

up to some exponentially decaying terms as $\xi\rightarrow 0^{+}$. According to Eq. (18), we get

[TABLE]

On the other hand, we can take the zero temperature limit $T\rightarrow 0^{+}$, where $\xi=2\pi v\beta/N_{y}\rightarrow+\infty$, and $q=\exp(-\xi)\rightarrow 0^{+}$ vanishes exponentially as $\xi$ approaches infinity. In this case, the character $\chi_{a}$ is dominated by the contribution of the primary state $a$, i.e., $\chi_{a}=q^{h_{a}-c_{-}/24}+O(q^{h_{a}+1-c_{-}/24})$, where $h_{a}$ denotes the conformal weight of the primary state. By making use of Eq. (18), we get

[TABLE]

up to some exponentially decaying terms.

Appendix B Majorana basis for non-interacting fermion systems

Here, we provide the details of our calculations in non-interacting spinless fermion systems, where we make use of the Majorana basis. The generalization to systems with multicomponent fermions at each site is straightforward. At each lattice site $i$, we introduce two Majorana fermions $f_{l}^{\dagger}=\frac{1}{2}(\gamma_{l1}+\mathrm{i}\gamma_{l2})$ and $f_{l}=\frac{1}{2}(\gamma_{l1}-\mathrm{i}\gamma_{l2})$ for the spinless fermion mode, where $\gamma_{l,\eta}$ $(\eta=1,2)$ is a Majorana operator that satisfies $\gamma_{a}^{\dagger}=\gamma_{a}$ and $\{\gamma_{a},\gamma_{b}\}=2\delta_{ab}$. Generally, under this basis, one can represent any quadratic Hamiltonian with $N$ sites using a $2N\times 2N$ real antisymmetric matrix, which is suitable to describe non-interacting fermion systems.

B.1 Thermodynamic calculation under the Majorana basis

For the quadratic Hamiltonian $\hat{H}=(\mathrm{i}/4)\sum_{a,b}\gamma_{a}M_{ab}\gamma_{b}$, where $M$ is a $2N\times 2N$ real antisymmetric matrix, $\gamma_{a}$ is the Majorana fermion operator, and $a,b$ are the indices of Majorana fermions. Any $2N\times 2N$ antisymmetric matrix $M$ can be block-diagonalized, with an orthogonal transformation, into the following form:

[TABLE]

where $O$ is the orthogonal matrix, and $\lambda_{l}>0$. By introducing the new Majorana fermion basis $\tilde{\gamma}\equiv O\gamma$, the Hamiltonian is expressed as $\hat{H}=(\mathrm{i}/4)\sum_{a,b}\tilde{\gamma}_{a}\tilde{M}_{ab}\tilde{\gamma}_{b}$, and we have

[TABLE]

where we have introduced $\gamma=(\gamma_{1,1},\gamma_{1,2},\gamma_{2,1},\ldots,\gamma_{N,2})^{T}$. By introducing the complex fermion $a_{l}=(\tilde{\gamma}_{l,1}-\mathrm{i}\tilde{\gamma}_{l,2})/2$ and $a_{l}^{\dagger}=(\tilde{\gamma}_{2l-1}+\mathrm{i}\tilde{\gamma}_{2l})/2$, we have $\mathrm{i}\tilde{\gamma}_{2l-1}\tilde{\gamma}_{2l}=-2a^{\dagger}_{l}a_{l}+1$ and then

[TABLE]

from which one can calculate the partition function as

[TABLE]

Under the Majorana basis, for our purposes, it is useful to introduce the correlation matrix $K$ as

[TABLE]

Under the orthogonal transformation $O$, the correlation matrix transforms as

[TABLE]

where $\tilde{K}_{de}=(\mathrm{i}/2Z)\mathrm{Tr}(\mathrm{e}^{-\beta\hat{H}}[\tilde{\gamma}_{d},\tilde{\gamma}_{e}])$. One can easily verify that $\tilde{K}$ is a $2\times 2$-block diagonalized matrix, with nonvanishing elements given by

[TABLE]

and $\tilde{K}_{(l,2),(l,1)}=-\tilde{K}_{(l,1),(l,2)}$, where $l=1,2,\ldots,N$. In the zero temperature limit, $\tilde{K}_{2k-1,2k}=-\tilde{K}_{2k,2k-1}=-1$.

With the correlation matrix $K$, one is able to calculate the thermal average of an observable $\hat{Q}=(\mathrm{i}/4)\sum_{a,b}\gamma_{a}M^{(Q)}_{ab}\gamma_{b}$ as

[TABLE]

where $M^{(Q)}$ is the real antisymmetric matrix corresponding to the quadratic Hermitian operator $\hat{Q}$.

B.2 Energy current operator in the Majorana basis

Generally speaking, the energy current operator defined in Sec. II is of the form of commutators of local Hamiltonian terms, i.e.,

[TABLE]

where $j,k$ are indices of local Hamiltonians. Under the Majorana basis, we introduce $\hat{H}_{j}=(\mathrm{i}/4)\sum_{a,b}\gamma_{a}M_{ab}^{(j)}\gamma_{b}$, where $M^{(j)}$ is the antisymmetric Hamiltonian corresponding to the local Hamiltonian $\hat{H}_{j}$. By a straightforward calculation, we get

[TABLE]

and after rearrangement of terms, we find

[TABLE]

Therefore, we get the antisymmetric matrix representation corresponding to the energy current operator $\hat{J}_{k\rightarrow j}$,

[TABLE]

and the thermal average of $\hat{J}_{k\rightarrow j}$,

[TABLE]

where $K$ is the correlation matrix.