text:
\section{Introduction \label{sec.1}}

Understanding quantum spin liquid (QSL) states and identifying them ~\cite{Balents2010,Savary_2016,RevModPhys.89.025003,Takagi2019} in materials, has become a centerpiece of research in modern condensed matter physics. The QSL is a topologically ordered state with long-range entanglement that appears in spin systems with frustrated magnetic interactions.  Although the beginning of the exploration of various QSL states dates back to the seminal work by Anderson on resonating valence bond states in antiferromagnetic Mott insulators~\cite{Anderson1973}, the research in QSLs rapidly intensified with the discovery of high-T$_c$ cuprates~\cite{Bednorz1986} due to its possible connection to the origin of superconductivity~\cite{Baskaran1987}. Over the last four decades, various theoretical proposals~\cite{PhysRevLett.59.2095,PhysRevB.37.3774,PhysRevB.62.7850,PhysRevLett.86.292} have been made to understand this exotic phase of matter. All these works hinge on a universal consensus about the fractional nature of quasiparticles in a typical QSL state. However, an unambiguous material-based realization of a QSL phase is hitherto absent.

Among the zoo of various predicted QSL states, the Kitaev spin liquid~\cite{Kitaev2006} plays a special role. It is an exactly solvable $S = 1/2$ spin model on a two-dimensional honeycomb lattice, in which the spins fractionalize into (i) Majorana fermions and (ii) $\mathbb{Z}_2$ vortices, also known as visons. Recently, substantial interest has emerged in the study of the Kitaev model owing to its potential realization in a few spin-orbit coupled Mott insulators such as complex iridates A$_2$IrO$_3$~\cite{PhysRevLett.108.127203} (A = Na, Li) and $\alpha$-RuCl$_3$~\cite{Banerjee2016}. For real materials, however, non-Kitaev interactions and further spin-exchange interactions are inevitably present~\cite{PhysRevLett.102.256403,PhysRevLett.105.027204,Bhattacharjee_2012,PhysRevLett.112.077204}, which spoils the integrability of the ideal Kitaev model and makes it difficult to assess its predictions. In this regard, various proposals have been put forward to tune the magnetic exchange interactions towards the ideal Kitaev limit by applying strain~\cite{PhysRevB.98.121107}, modifying spin-orbit coupling (SOC)~\cite{PhysRevB.87.041107} or through Floquet engineering~\cite{PhysRevB.103.L100408,PhysRevB.105.085144,PhysRevB.105.L180414,Kumar2022} in the honeycomb Mott insulators.

Recent experiments on $\alpha$-RuCl$_3$~\cite{PhysRevLett.119.037201,PhysRevB.95.180411,PhysRevB.96.041405,PhysRevLett.118.187203,Banerjee2018,PhysRevLett.120.117204,PhysRevB.100.060405,Czajka2021,Tanaka2022,Bruin2022} have provided a strong evidence for the existence of the Kitaev QSL phase in an intermediate magnetic field regime of $B \sim 8 - 14$ T. Namely, nuclear magnetic resonance~\cite{PhysRevLett.119.037201} and neutron diffraction~\cite{PhysRevB.95.180411} experiments observed field-induced spin gap opening around $10$ T consistent with the physics of an ideal Kitaev model in an applied magnetic field. Particular features of such gap opening have also been confirmed by inelastic neutron scattering studies~\cite{Banerjee2018}. Specifically, the observation of half-quantized thermal Hall conductivity~\cite{Kasahara2018,Bruin2022} has been interpreted as a signature of chiral Majorana edge modes in $\alpha$-RuCl$_3$. On the other hand, the observation of thermal conductivity oscillations in $\alpha$-RuCl$_3$ as a function magnetic field, for an intermediate magnetic-field regime, suggests a metallic nature of the underlying quasiparticles, despite that the material is a perfect Mott insulator with a large Mott gap~\cite{Czajka2021}. Nonetheless, direct experimental observation of the Kitaev phase has remained somewhat controversial~\cite{Trebst2021}. The reason is that these experimental probes provide, mostly, an indirect signature of the underlying excitations of the Kitaev model, and cannot unambiguously determine the nature of the quasiparticles. Hence, the current \textit{state-of-the-art} research demands an alternative, more direct, way to detect the Majorana and vison excitations in the Kitaev QSL phase.

Motivated by these fascinating experimental findings, here, we consider a different route to \textit{electrically} detect the signatures of the fractionalized \textit{neutral} excitations in the Kitaev model. A priori, such an approach might seem counter-intuitive as the parent compounds, realizing proximate Kitaev physics, are Mott insulators. However, in the paper by Bulaevskii \textit{et al}., it was shown that certain frustrated Mott insulators can exhibit non-zero \textit{spontaneous circular electrical currents} or \textit{nonuniform charge distribution} in the ground state~\cite{PhysRevB.78.024402}. On the other hand, spontaneous breakdown of injected electrons on a QSL material, into the fractionalized excitations, can be used to detect Majorana fermions~\cite{Barkeshli2014}. A previous theoretical work~\cite{PhysRevLett.125.227202} has recently proposed an experimental set-up utilizing the local charge distribution to detect vison excitations in $\alpha$-RuCl$_3$; whereas another contemporary work~\cite{PhysRevX.10.031014} proposed a more exotic electrical access to the Majorana fermions in Kitaev materials utilizing the transmutation protocol~\cite{PhysRevX.2.031013,PhysRevB.87.045130,Barkeshli2014}.

In this continuing search for various innovative electrical access in Mott insulators, one very exciting aspect still remains unexplored \--- \textit{the induced electrical loop currents}, which is allowed in QSL without time-reversal symmetry. Here, we develop a theoretical formalism to analyze such localized orbital loop current profile in a multiorbital Mott insulator relevant for $\alpha$-RuCl$_3$ and propose experimental set-ups to detect the signatures of Majorana and vison excitations in a Kitaev quantum magnet. For a half-filled single orbital Hubbard model on a 2D lattice with $|t_{ij}| \ll U$, with an on-site Coulomb repulsion $U$ and the inter-site hopping amplitude $t_{ij}$, the corresponding expression for the loop current operator reads as~\cite{PhysRevB.78.024402}
\begin{equation}\label{eq.1}
\bm{\mathcal{I}}_{ij,k} = \frac{24t_{ij} t_{jk} t_{ki}}{U^2} \mathbf{S}_k \cdot (\mathbf{S}_i \times \mathbf{S}_j),
\end{equation}
where $\mathbf{S}_i = (S_i^x, S_i^y, S^z_i)$ is the spin-1/2 operator, and $(ijk)$ labels the site indices on the smallest triangular plaquette within the lattice. Eq.~\eqref{eq.1} can be understood from symmetry consideration. Both the spin chirality $\mathbf{S}_k \cdot (\mathbf{S}_i \times \mathbf{S}_j)$ and current are odd under time reversal operation and spatial inversion (i.e. interchange of the indices $i$ and $j$). The system also has SU(2) spin rotation symmetry. Unless explicitly mentioned, we use natural units $\hbar=e=c = 1$ throughout this paper.

In this work, we investigate the induced localized loop current profile in a multiorbital spin-orbit coupled Hubbard-Kanamori model relevant for Kiteav quantum magnets such as $\alpha$-RuCl$_3$. The key findings of our work are listed as follows: (i) Utilizing the Schrieffer-Wolff transformation (SWT), we derive the expression for the localized loop current operator in the  Mott insulator to the third order in perturbation expansion [see Eq.~\eqref{eq.4.1}]. The breakdown of the SU(2) spin rotation symmetry, due to spin-orbit coupling, leads to a different loop current expression compared to its single-orbital counterpart.

(ii) In the presence of an external magnetic field, the Kitaev ground state hosts non-zero expectation value of such localized loop currents. Adopting microscopic parameters relevant for $\alpha$-RuCl$_3$~\cite{PhysRevB.93.155143,Sinn2016}, we obtain the orbital current profile in a finite 2D honeycomb lattice with localized edge and bulk loop current distributions, as shown in Fig.~\ref{fig:Fig2}. The induced localized currents in the bulk of a 2D honeycomb lattice constitute an emergent super-structure of two inter-penetrating triangular lattices as illustrated in Fig.~\ref{fig:Fig2}(b,b$'$). When an in-plane [$ab$-plane in Fig.~\ref{fig:Fig1}(a)] magnetic field is applied to the Kitaev system, these emergent orbital currents produce an out-of-plane ($c$-axis) magnetic field, which we identify as the \textit{key feature} of a Kitaev system.

(iii) Finally, we notice that vison excitations drastically modify the current distribution profile and the associated orbital magnetization locally, as shown in Fig.~\ref{fig:Fig3}. We provide quantitative estimates for this change and about the possible measurement of such orbital magnetization using various possible experimental techniques. Our work provides an important platform to directly detect the electrical signatures of the concomitant excitations (Majorana and vison) in a Kitaev quantum magnet.

\section{Phenomenological picture \label{sec.pheno}}

Before moving to full microscopic model analysis, here we provide a more intuitive picture of why we expect an orbital magnetization in Mott insulators hosting the Kitaev quantum spin liquid in terms of the parton theory~\cite{PhysRevB.29.3035}. The electron operator can be written as $c_\sigma=b f_\sigma$, where $b$ is a boson operator carrying electron charge quantum number and $f_\sigma$ is the fermionic spinon operator carrying spin quantum number. The parton decomposition leads to an emergent gauge field $\mathbf{a}$, which inherently couple to the spinons as $f_\sigma \rightarrow f_\sigma \exp(i a)$, while the charged boson couple to both the physical gauge field $\mathbf{A}$ and emergent gauge field $\mathbf{a}$, as $b\rightarrow b \exp[i (A-a)]$. The effective low-energy Lagrangian for the $b$ boson has the standard Ginzburg-Landau form~\cite{Chowdhury2018}
\begin{equation}\label{eq.pheno.1}
\mathcal{L}_b=|(i\partial_\mu+a_\mu-A_\mu) b|^2-g |b|^2-\frac{u}{2}|b|^4+\cdots.
\end{equation}
When $b$ boson condenses for $g<0$, the Anderson-Higgs mechanism generates a term proportional to $(a_\mu-A_\mu)^2$, which locks $\mathbf{a}$ to $\mathbf{A}$~\cite{PhysRevB.62.7850}. As a consequence, the emergent gauge field behaves the same way as the physical gauge field, and the system is metal if $f_\sigma$ forms the fermi surface or a superconductor if $f_\sigma$ forms Cooper pairs and condenses.

In the Mott insulator, $b$ boson is gapped with $g>0$. However, there is still diamagnetic response in $\mathbf{A}-\mathbf{a}$ due to the gapped charged boson, similar to the Landau diamagnetism in metal, albeit weaker. When an external magnetic is applied in Mott insulators, local current loops are induced which generate magnetization opposite to the applied field as a diamagnetic response. The diamagnetic response increases with decreasing charge gap, and hence the effect discussed here is more prominent near the Mott transition. The effective Lagrangian of the system after integrating out $b$ boson has the form~\cite{PhysRevB.97.045152}
\begin{equation}\label{eq.pheno.2}
\mathcal{L}=\mathcal{L}_f(f_\sigma, \mathbf{a})-\frac{\chi_b}{2}[\nabla\times(\mathbf{a}-\mathbf{A})]^2-\frac{\chi_{\mathrm{B}}}{2}(\nabla\times \mathbf{A})^2,
\end{equation}
where $\chi_b$ accounts for the diamagnetic susceptibility due to the gapped boson $b$, $\chi_{\mathrm{B}}$ is the susceptibility of the background. $\mathcal{L}_f(f_\sigma, \mathbf{a})$ is the Lagrangian describing the spinons coupled to the emergent gauge field $\mathbf{a}$. Due to the diamagnetic response, an emergent magnetic field induces a physical magnetic field, which is obtained by minimizing $\mathcal{L}$ with respect to $\mathbf{B}\equiv \nabla\times \mathbf{A}$: $\mathbf{B}=(\nabla\times\mathbf{a})\chi_b/(\chi_b+\chi_{\mathrm{B}})$. We remark one subtle point that $\mathbf{a}$ is compact meaning that the system is invariant under $\mathbf{a}\rightarrow \mathbf{a}+2\pi$. Because of the compactness of $\mathbf{a}$, the flux of $\mathbf{a}$ is defined upto modulo of $2\pi$. In the context Kitaev Mott insulator, therefore a vison carrying $\pi$ flux without magnetic field is time-reversal symmetric and does not induce a nonzero $\mathbf{B}$. Here we have considered isotropic response. For the two-dimensional systems we consider below, the diamagnetic susceptibility is only for the field component perpendicular to the plane, and the induced $\mathbf{B}$ is also along the direction perpendicular to the plane.

In the spinon description, the chiral Kitaev quantum spin liquid stabilized by a magnetic field corresponds to the state with $f_\sigma$ fermions being in the $p+ip$ superconducting state, where the time-reversal symmetry is explicitly broken by the magnetic field~\cite{PhysRevB.84.125125}. In this phase, there exists chiral edge current of $f_\sigma$ and the associated magnetization of $\nabla\times \mathbf{a}$. At the same time, the vortex carries flux of $\mathbf{a}$. Because of the diamagnetism due to the gapped $b$ boson, the emergent magnetic field induces a physical magnetic field. This picture will be corroborated through the explicit calculation of a microscopic model below.

\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{Fig1}
\caption{(a) Schematic of a finite size ($L \times L$) honeycomb lattice made out of magnetic ions with two $\mathbb{Z}_2$ vortices (vison excitations) [\textit{shaded region}]. The vison excitations are separated by a string operator  (\textit{dashed blue line}) which flips $u_{ij}$ eigenvalue when it crosses a particular bond $\braket{ij}$. The colored bonds represent three anisotropic bond-dependent Kitaev interactions. The $\pi$-rotation symmetry along the crystalline axes $b$ is illustrated by the dashed orange lines. There is an additional three-fold rotation symmetry along the crystalline $c$-axis. Note that the $\pi$-rotation along the $a$-axis is not a symmetry. (b) A single hexagonal plaquette with the (next)- nearest neighbor hopping amplitudes ($t_2'$) $t_2$ formed by the bond-dependent orbital overlaps. (c) The relative orientation of the honeycomb plane ($abc$-coordinate system) in panel (a) with respect to the spin quantization axes ($xyz$-coordinate system).}\label{fig:Fig1}
\end{figure}

\section{Microscopic model \label{sec.2}}

We start from the microscopic multiorbital description for $\alpha$-RuCl$_3$, however, this formalism can be suitably generalized to other Kitaev candidate materials. The electronic configuration of the transition metal (TM) ion Ru in $\alpha$-RuCl$_3$ is $4d^{5}$. These \textit{five} electrons reside in the $t_{2\mathsf{g}}$ manifold formed by the three orbitals: $d_{xy}$, $d_{yz}$ and $d_{zx}$. The presence of strong atomic SOC~\cite{PhysRevB.91.144420,PhysRevLett.102.017205,PhysRevLett.110.097204,PhysRevB.93.214431,Winter_2017,PhysRevB.95.024426} further splits this manifold into $J_{\mathrm{eff}} = 1/2$ and $J_{\mathrm{eff}} = 3/2$ states. This leads to a completely filled $J_{\mathrm{eff}} = 3/2$ manifold with the four electrons, and the remaining electron resides in the $J_{\mathrm{eff}} = 1/2$ sector. Consequently, it leads to an effective one-hole model, constructed with the following parameters: an on-site Coulomb repulsion of strength $U$, the Hund's coupling $J_{\mathrm{H}}$, and (next-)nearest neighbor hopping ($t_2'$) $t_2$, respectively~\cite{PhysRevLett.102.017205,PhysRevLett.112.077204} [see Fig.~\ref{fig:Fig1}(a,b)].

In the strong-coupling limit ($t_2,t_2' \ll U, J_{\mathrm{H}}$), the low-energy effective Hamiltonian is obtained by the (SWT)~[see supplemental material (SM)\cite{supp} for more details] and leads to the ferromagnetic Kitaev Hamiltonian as
\begin{equation}
\label{eq.2}
\mathcal{H}_0  = - K \sum_{\langle ij \rangle; \alpha} S^{\alpha}_i S^{\alpha}_j, \quad K = \frac{8J_{\mathrm{H}} t_2^2}{3(U^2 - 4UJ_{\mathrm{H}} + 3J_{\mathrm{H}}^2)},
\end{equation}
where $\langle ij \rangle$ denotes the nearest-neighbor sites, and $\alpha \in \{x,y,z\}$ labels the planar orientation of the bond-type as illustrated in Fig.~\ref{fig:Fig1}(a)-(b). The Kitaev coupling $K$ in Eq.~(\ref{eq.2}) is obtained to the lowest order in perturbation, and does not include $t_2'$.

This model can be solved exactly, as shown by Kitaev himself~\cite{Kitaev2006}, by introducing Majorana representation of the spin-1/2 operators as $S_i^{\alpha} = ib_i^{\alpha} c_i/2$ for $\alpha = x,y,z$. Here $b^{\alpha}_i, c_i$ are four mutually anti-commuting Majorana operators. In terms of $b^{\alpha}_i$ and $c_i$, the Hamiltonian in Eq.~(\ref{eq.2}) can be rewritten as $\mathcal{H}_0 = iK/4 \sum_{\langle ij \rangle} u_{ij} c_i c_j$, where $u_{ij} (= ib^{\alpha}_i b^{\alpha}_j)$ are the $\mathbb{Z}_2$ gauge field operators with eigenvalues $\pm 1$. Since $u_{ij}$ commutes with the Hamiltonian, one can pick a particular choice of $u_{ij}$ (a gauge sector $\ket{\mathcal{G}}$) and solve the remaining non-interacting Majorana Hamiltonian. The pinned flux in each honeycomb plaquette is obtained by the gauge-invariant loop operator $\mathcal{W}_p = \prod_{\braket{ik} \in \mathord{\raisebox{-0.5pt}{\tikz{\node[draw,scale=.5,regular polygon, regular polygon sides=6,](){};}}}_p} u_{ik}$. Here $\mathcal{W}_p = +1$ corresponds to zero-flux in a plaquette, whereas $\mathcal{W}_p = -1$ signifies a $\pi$-flux and corresponds to a $\mathbb{Z}_2$ vortex or a vison excitation [see the shaded hexagons in Fig.~\ref{fig:Fig1}(a)]. Following Lieb's theorem~\cite{PhysRevLett.73.2158}, the ground state of the Majorana fermions is obtained with the gauge choice of $\mathcal{W}_p = 1, \forall p$, for all the hexagonal plaquettes. Starting from a uniform gauge configuration, one may create a single vison excitation by changing the eigenvalue of a single plaquette loop operator to $-1$.

To induce orbital magnetization, it is necessary to break the related symmetries. In crystals realizing the Kitaev model, the symmetry is low because of the spin-orbit coupling. As shown in Fig.~\ref{fig:Fig1}(a), the system has translation symmetry, two-fold rotation symmetry along crystalline $b$-axis denoted by $\mathsf{C}_{2b}$ in panel (a), 3-fold rotation symmetry along the $c$ axis, $\mathsf{C}_{3c}$, time-reversal symmetry (TRS) and inversion symmetry. The seemingly two-fold rotation along the $a$-axis Fig.~\ref{fig:Fig1}(a) is absent when one embeds the honeycomb plane into the parent crystal. To allow for orbital magnetization along the $c$ axis, we need to break the TRS and also the $\mathsf{C}_{2b}$ symmetry. This can be achieved by applying a magnetic field with non-zero component perpendicular to the $b$ axis. One particularly interesting case is when the magnetic field is applied in the $a$-$b$ plane, which induces orbital magnetization along the $c$ axis. Generally, the induced magnetization is much weaker than the applied magnetic field. The case with an external magnetic field in the plane may be more convenient for experimental detection of the induced orbital magnetization.

\section{Circulating loop currents \label{sec.3}}

Since the parent compound of a Kitaev spin liquid is a multiorbital spin-orbit coupled Mott insulator, there is remnant charge response and there can exist localized circulating electrical currents. The charge response is governed by the charge excitation gap, which is of the order of $U$. Similar to the derivation of the effective spin Hamiltonian in Eq. \eqref{eq.2} from the Hubbard model, we perform a strong-coupling expansion (SWT) for the current operator defined on the bond $\langle ij \rangle$ [see Fig.~\ref{fig:Fig1}(b)] as
\begin{equation}\label{eq.3}
\bm{\mathcal{I}}_{ij} = \frac{iet_2'\hat{\ensuremath{\mathbf{r}}}_{ij}}{\hbar} \sum_{\alpha,\beta,\sigma} \left( d^{\dagger}_{j\alpha \sigma}d_{i \beta \sigma} - d^{\dagger}_{i\beta \sigma}d_{j \alpha \sigma} \right),
\end{equation}
where $d^{\dagger}_{i\alpha \sigma}$ creates an electron in the $i$-th site with spin $\sigma$ and orbital $\alpha,\beta \in \{ x,y,z\}$, and obtain a generalization of Eq.~(\ref{eq.1}) for the triangular loop current as (see the details of the derivation in Sec.~(II) in the SM~\cite{supp})
\begin{subequations}
\begin{align}
\nonumber
\tilde{\bm{\mathcal{I}}}_{ij,k}
& =
	  \bm{\mathcal{I}}_0  (3S^x_i S^y_j +
	 			        S^y_i S^x_j) S^z_k +	\\
\nonumber
& 	  \quad \quad
	  \bm{\mathcal{I}}_0  (3S^y_i S^z_j -
					    5S^z_i S^y_j ) S^x_k	+	\\
\label{eq.4.1}
& 	  \quad \quad \quad
	  \bm{\mathcal{I}}_0  (3S^z_i S^x_j	-
					    5S^x_i S^z_j) S^y_k , \\
\label{eq.4.2}
\bm{\mathcal{I}}_0 & = \frac{e\hat{\ensuremath{\mathbf{r}}}_{ij}}{\hbar}\frac{8t_2^2t_2'J_{\mathrm{H}}(U - 2J_{\mathrm{H}})}{9(U^2 - 4 UJ_{\mathrm{H}} + 3 J_{\mathrm{H}}^2)^2},
\end{align}
\end{subequations}
where $S^{\alpha}_i$ is the $\alpha$-th component of the spin at site $i$, and $U,J_{\mathrm{H}}$, and $t_2$ are the microscopic parameters as defined earlier. For the current operator defined on the other bonds, we need to modify the loop current expression by permuting both the orbital and site indices on the three spin operators. The expression \eqref{eq.4.1} is consistent with the $\mathsf{C}_{2b}$ symmetry [orange dashed line in Fig. \ref{fig:Fig1} (b)]. Under $\mathsf{C}_{2b}$, $i\leftrightarrow j$, $k\leftrightarrow k$, $S^x \rightarrow -S^y$, $S^y \rightarrow -S^x$, $S^z \rightarrow -S^z$ and $\tilde{\bm{\mathcal{I}}}_{ij,k} =- \tilde{\bm{\mathcal{I}}}_{ji,k}$. However, the $\mathsf{C}_{2b}$ symmetry cannot uniquely determine the form of $\tilde{\bm{\mathcal{I}}}_{ij,k}$. Actually the absence of any terms with repeated index in the spin space (i.e. $S_i^x S_j^x S_k^z$) in Eq.~(\ref{eq.4.1}) is an artifact of considering an ideal situation, where we retain only $t_2$, and $t_2'$ hoppings in the multiorbital tight-binding (TB) description. Indeed, a more realistic TB modeling~\cite{PhysRevLett.112.077204} for Kitaev quantum magnets would lead to additional three-spin terms with possible repeated indices. However, for the subsequent analysis in this work, we focus on the ideal Kitaev limit.

For an estimate of the amplitude of the induced electrical current, we adopt all the parameters entering Eq.~(\ref{eq.4.2}) from the recent \textit{ab initio}~\cite{PhysRevB.93.155143} and \textit{photoemission} reports~\cite{Sinn2016} for $\alpha$-RuCl$_3$ as $U = 3.0$ eV, $J_{\mathrm{H}} = 0.45$ eV, $t_2 = 0.191$ eV, and $t'_2 = -0.058$ eV. The value of $t'_2$ has been adopted from Ref.~\cite{PhysRevB.100.075110}. Plugging in the numbers in Eqs.~(\ref{eq.2}) and (\ref{eq.4.2}), we estimate $|\bm{\mathcal{I}}_0| \sim 30$ nA, with a Kitaev coupling $K \sim 0.023$ eV.

The key quantity now is the expectation value of the induced loop current operator in the ground state of the Kitaev model. The pure Kitaev model can be solved exactly from the Hamiltonian $\mathcal{H}_0 = iK/4 \sum_{ij} c_i c_j$, with all $u_{ij}$ chosen to be $+1$. The ground state eigenfunction can be written as $\ket{\Psi_0} = \mathcal{Z}\ket{\mathcal{G}_0}\otimes \ket{\mathcal{M}}$ where $\ket{\mathcal{G}_0}$ corresponds to $\mathcal{W}_p = 1$ for all hexagonal plaquettes and $\ket{\mathcal{M}}$ is the ground state of \textit{gapless} Majorana fermions~\cite{supp}. However as the Hilbert space is enlarged due to the fractionalization of the spin degrees of freedom, we need to project the eigenstates of $\mathcal{H}_0$ to the physical Hilbert space of the spin Hamiltonian in Eq.~(\ref{eq.2}). The corresponding projection is achieved by the local constraint $\mathcal{D}_i = i c_i b_i^x b_i^y b_i^z = 1$ and the global projection operator $\mathcal{Z} = \Pi_i(1 + \mathcal{D}_i)/2$~\cite{Kitaev2006,PhysRevB.84.165414,PhysRevB.92.014403}.

Due to the static nature of the vison excitations in a pure Kitaev model, only the first term, $S_i^xS_j^yS_k^z$, in Eq.~\eqref{eq.4.1} contributes to the computation of the loop current expectation value in the Kitaev ground state. All the other five terms in Eq.~\eqref{eq.4.1} do not contribute as they have zero expectation value in $\ket{\Psi_0}$ (see the SM~\cite{supp} for more details). It is straightforward to see that the loop current expectation value vanishes in the pure Kitaev model ground state, i.e., $\braket{\Psi_0|\tilde{\bm{\mathcal{I}}}_{ij,k}|\Psi_0} = 0$. This absence is not a surprise, since any finite loop current would lead to an orbital magnetization that breaks TRS, whereas the ground state $\ket{\Psi_0}$ preserves the TRS.

\subsection{External magnetic field \label{sec.4}}

The TRS of the system is broken in the presence of an external magnetic field. In this case, we may expect a finite loop current expectation value in the ground state. We consider the Kitaev model in an external magnetic field and explore its consequences. The integrability of the model is destroyed in the presence of an external magnetic field $\mathbf{h} = (\mathrm{h}_x,\mathrm{h}_y,\mathrm{h}_z)$, defined along spin quantization axes (we call it the $xyz$-coordinate system). For a small $|\mathbf{h}|$, we can treat the effect of $\mathbf{h}$ in a perturbative manner. The lowest order perturbation correction, that breaks the TRS~\cite{Kitaev2006}, reads as $\mathcal{H}_{\mathrm{eff}} = -\kappa \sum_{ijk, \triangle} S^{x}_{i}S^{y}_{j}S^{z}_{k}$, where $\kappa = \mathrm{h}_x \mathrm{h}_y \mathrm{h}_z/K^2$. Note that $\mathcal{H}_{\mathrm{eff}}$ has the same form as the first term in $\tilde{\bm{\mathcal{I}}}_{ij,k}$, and therefore it is natural to have a nonzero current in the QSL. Choosing the same gauge configuration $\ket{\mathcal{G}_0}$ as before, the total Hamiltonian $\mathcal{H}=\mathcal{H}_0 + \mathcal{H}_{\mathrm{eff}}$ can be simplified as
\begin{equation}\label{eq.5}
\mathcal{H} = \frac{iK}{4} \sum_{\langle ik \rangle} c_i c_k + \frac{i\kappa}{8} \sum_{\llangle ij \rrangle} c_i c_j,
\end{equation}
where all the $u_{ij}$'s are replaced by their eigenvalues. The ground state eigenfunction of the above Hamiltonian can be written as $\ket{\Psi_{\mathrm{mag}}} = \mathcal{Z} \ket{\mathcal{G}_0} \otimes \ket{\mathcal{M}_{\mathrm{mag}}}$, where $\ket{\mathcal{M}_{\mathrm{mag}}}$ is the ground state of the \textit{gapped} Majorana fermions~\cite{supp}, and $\mathcal{Z}$ is the global projection operator, as defined earlier. The loop current expectation value becomes non-vanishing in this case and $\braket{\Psi_{\mathrm{mag}}|\tilde{\bm{\mathcal{I}}}_{ij,k}|\Psi_{\mathrm{mag}}} \propto \kappa$ for $\kappa \ll K$ (see the SM~\cite{supp} for more details).

One key feature of the Kitaev system is that even an in-plane magnetic field [applied in the $ab$-plane shown in Fig.~\ref{fig:Fig1}(b)] can have non-zero components along all the three spin quantization axes. Let us write the external magnetic field $\mathbf{h}^{\mathrm{ext}}$ in the coordinate system containing the honeycomb plane (we call it the $abc$-coordinate system) as $\mathbf{h}^{\mathrm{ext}} = \mathrm{h}(\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)$, where $\phi$ and $\theta$ are azimuthal angle measured from the crystalline $a$ axis and polar angle from the crystalline $c$ axis, respectively, and $\mathrm{h}$ is the strength of the applied magnetic field. After a straightforward transformation between the $abc$- and the $xyz$-coordinate system [see Fig.~\ref{fig:Fig1}(c)], we obtain~\cite{supp}
\begin{widetext}
\begin{equation}\label{eq.6}
(\mathrm{h}_x, \mathrm{h}_y, \mathrm{h}_z)
=
\mathrm{h} \left(
\frac{\cos \theta}{\sqrt{3}} + \frac{\sin\theta \cos \phi}{\sqrt{6}} - \frac{\sin \theta \sin \phi}{\sqrt{2}}, \;
\frac{\cos \theta}{\sqrt{3}} + \frac{\sin\theta \cos \phi}{\sqrt{6}} + \frac{\sin \theta \sin \phi}{\sqrt{2}}, \;
\frac{\cos \theta}{\sqrt{3}} - \sqrt{\frac{2}{3}} \sin\theta \cos \phi \right).
\end{equation}
\end{widetext}
It is evident from Eq.~\eqref{eq.6} that even for an in-plane magnetic field (with $\theta = 0$), one can have non-zero $(\mathrm{h}_x,\mathrm{h}_y,\mathrm{h}_z)$, leading to a finite $\kappa$. By rotating the magnetic field in the $ab$-plane, $\kappa$ can change the sign and even become zero for certain angles (see the SM~\cite{supp} for the details about the variation of $\kappa$ as a function of this rotation angle). For the field applied along the $b$ axis, $\theta=\pi/2$ and $\phi=\pi/2$, we have $\kappa=0$ and therefore the induced orbital magnetization vanishes, which is consistent with the symmetry analysis above.

When the magnetic field strength becomes large, the perturbative treatment is inadequate. In addition to the $\kappa$ term, the magnetic field generates dynamics for the visons, which spoil the integrability of the model. To make analytical progress, we take $\kappa$ as a free parameter with the understanding that it relates to the magnetic field for a weak field. The main results of the orbital magnetization obtained below are expected to hold even when visons are dynamical. For a large magnetic field, the chiral Kiteav QSL is destroyed and the description in Eq. \eqref{eq.5} is no longer applicable~\cite{Patel2019,Hickey2019}.

So far we focused on a specific triangular plaquette, however, in a lattice geometry, the situation becomes more interesting. A specific triangle is shared between three honeycomb plaquettes, whereas each nearest-neighbor (NN) bond on the triangle is shared between four triangles. Among these four triangles, two belong to the hexagonal plaquette containing the next-nearest neighbor (NNN) bond $\langle ij \rangle$, where the current operator is defined, and the other two belong to the neighboring hexagonal plaquettes (see Fig.~\ref{fig:Fig1}). Consequently, the net current along the NN bonds forming the honeycomb plaquettes vanishes because of counter-propagating currents from the two honeycomb plaquettes for translationally invariant systems. In a finite size system with open boundary conditions (OBC), however, such perfect cancellations are not possible near the edge terminations, as shown in Fig.~\ref{fig:Fig2}. On the other hand, since the NNN bonds are not shared between any plaquettes, the localized currents are never canceled and they create an interesting loop current profile distributed over the entire system, as will be discussed in the next section.

\begin{figure*}[t]
\centering
\includegraphics[width=1\linewidth]{Fig2}
\caption{(a) Distribution of the localized edge currents around the zigzag type edges for a finite system with OBC of size $L = 21$ (the bulk of the system is illustrated by the dashed bonds). Note that the localized currents on four different edges are displayed only for the bonds lying on the honeycomb plaquettes. The localized currents on the next-nearest neighbor (NNN) bonds are not shown for simplicity. (b) The macroscopic build-up of localized currents along the NNN bonds is distributed over the entire sample size and forms \textit{filament-like} current textures. A magnified view of the details of the current textures in the middle of the system is shown in the top right panel (b$'$). (c) The variation of the localized current as a function of $\kappa$ for the nearest-neighbor (NN) bonds on the middle of the  edges $\mathrm{E}_1 - \mathrm{E}_4$ (\textit{open triangles}) and NNN bonds inside the lattice (\textit{open circles}), respectively. (d) The current profile around a different geometry with both the zigzag and arm-chair edges. The numerical estimates are provided with $\kappa=0.2 K$.}\label{fig:Fig2}
\end{figure*}

\section{Finite system analysis \label{sec.5}}

We now focus our analysis on a finite system to compute the localized current profile. We first consider a finite system (with OBC) of linear size $L$ without any $\mathbb{Z}_2$ vortices or vison excitations. We set the value of $\kappa$ motivated by the recent experimental work on $\alpha$-RuCl$_3$~\cite{Tanaka2022}. Comparing the gap magnitude from Eq.~(\ref{eq.5}) with the observed field-dependent Majorana gap at an applied magnetic field of $\sim 10$ T, we estimate $\kappa \sim 0.77 K$. However, this value of $\kappa$ corresponds to a large magnetic field, when our approximate Hamiltonian in Eq.~\eqref{eq.5} is no longer applicable. Hence, we use a smaller magnitude of $\kappa=0.2K$ for subsequent quantitative estimates, unless otherwise mentioned. As explained previously, the localized current on a particular bond $\braket{ij}_{\alpha}$ on the honeycomb plaquette is the sum of the contributions from the associated triangles. Following the analytical form in Eq.~(\ref{eq.4.1}), it appears as if there are three different average currents $\mathcal{I}_{\alpha}$ in three different bonds, $\alpha = x,y,z$, respectively [see Fig.~\ref{fig:Fig1}(a)]. However, the $C_3$ rotation symmetry of the underlying lattice implies that all these three currents are identical \textit{i.e.}, $\mathcal{I}_{\alpha} = \mathcal{I}$. In addition, neighboring honeycomb plaquettes of a particular bond host counter-propagating current $\mathcal{I}$, leading to vanishing average current distribution at the sides of each honeycomb plaquette in the bulk. However, such perfect cancellations are absent near the edge terminations of the finite system, which consequently leads to non-vanishing localized currents around the edge as shown in Fig.~\ref{fig:Fig2}(a) and \ref{fig:Fig2}(c). The lack of translational invariance around the edges of the system leads to different localized currents along different bonds in each of the triangular plaquettes along the edge. The total current is conserved at each vertex.

On the other hand, the localized currents in each of the NNN bonds $\braket{ij}$ within the triangle $\triangle_{ijk}$ do not vanish as they are not shared between other triangles in the lattice geometry. It leads to a macroscopic build-up of the localized currents which mimics a \textit{filament} like structure distributed over the entire system, as shown in Fig.~\ref{fig:Fig2}(b). For an infinite system, these localized currents on the NNN bonds, form an emergent superlattice of two interpenetrating triangular lattices as illustrated in Fig.~\ref{fig:Fig2}(b$'$). We emphasize that there is no free-flowing current over long distances because of the Mott nature of the system. Currents only flow around a closed loop of atomic size. The seemingly free flow of current (red line in Fig. \ref{fig:Fig2}(b)) is an illusion due to the superposition of localized current loops. The dependence of the current amplitude as a function of $\kappa$ is shown in Fig.~\ref{fig:Fig2}(c). In the subsequent sections Sec.~\ref{sec.6}, and ~\ref{sec.7}, we focus on two types of edge-termination geometry for finite systems with OBC in the absence of vison excitations. In Sec.~\ref{sec.8}, we again consider a finite system with PBC in the presence of two well-separated $\mathbb{Z}_2$ vortices, and illustrate the drastic modification of the bulk loop currents around these vison excitations.

\subsection{Quantitative estimates and analysis \--- case of zigzag geometry \label{sec.6}}

We performed numerical diagonalization for a system of linear size $L = 21$ with zigzag edge terminations to obtain the average localized current distribution without any vison excitations. In this case, the thermodynamic limit is reached for a relative small system size $L \ge 20$ because the system is fully gapped. (see the SM ~\cite{supp} for the technical details). The corresponding results are shown in Fig.~\ref{fig:Fig2}(a-c). In this geometry, the four different zigzag edges $\mathrm{E}_1 - \mathrm{E}_4$ are related to each other by the three-fold rotation along the $c$-axis, and we focus on one edge viz. $\mathrm{E}_1$ for the quantitative analysis.

As shown in Fig.~\ref{fig:Fig1}(a), the edge $\mathrm{E}_1$ is formed by two types of bonds: $x$- and $z$-bond. In the thermodynamic limit, we obtain a localized current that saturates at $\mathcal{I}^{(0)}_{\mathrm{edge}} \sim 3.0$ nA for both the bonds, that lie far away from the corner sites with $\mathrm{E}_2$ and $\mathrm{E}_4$, respectively. The variation of the the thermodynamic edge current $\mathcal{I}^{(0)}_{\mathrm{edge}}$ as a function of $\kappa$ is shown in Fig.~\ref{fig:Fig2}(c) (see the variation with open triangles).

The value of $\mathcal{I}_{\mathrm{edge}}$ deviates toward the corners of the system and differs for each bond type because of the lack of translation invariance along the edge. The amplitude for $z (x)$-bond near the corner with $\mathrm{E}_4$ is $\sim 4.57$ ($3.83$) nA, whereas near the other corner with $\mathrm{E}_2$ it is slightly less $\sim 2.21$ ($1.58$) nA as shown in Fig.~\ref{fig:Fig2}(a). This variation along the edge is illustrated by the color legend below panel (a) in Fig.~\ref{fig:Fig2}. The saturation length (the distance beyond which the bond currents become $\mathcal{I}^{(0)}_{\mathrm{edge}}$) near the edge with $\mathrm{E}_4$ is slightly larger than the saturation length near the edge with $\mathrm{E}_2$. Furthermore, the currents on the NN bonds inside the system almost vanish as mentioned earlier.

On the other hand, the current amplitude on the NNN bonds [$x',y',z'$-type corresponding to $t_2'$ hopping, see Fig.~\ref{fig:Fig1}(b)] inside the lattice far away from the edges and corners saturates at $\mathcal{I}^{(0)}_{\mathrm{bulk}} \sim 7.70$ nA as shown in Fig.~\ref{fig:Fig2}(b). In comparison to the edge current $\mathcal{I}^{(0)}_{\mathrm{edge}}$, the variation of this thermodynamic bulk localized current $\mathcal{I}^{(0)}_{\mathrm{bulk}}$ as a function of $\kappa$ is shown in Fig.~\ref{fig:Fig2}(c) (see the variation with open circles). This value also deviates from the bulk value as the bonds come close to the edges/corners of the system. The $\mathsf{C}_{3c}$ rotation symmetry, which enforces identical localized currents $\mathcal{I}^{(0)}_{\mathrm{bulk}}$ on all $x',y',z'$ bonds in the bulk, is lost as we go away from the middle of the system to its edge. Consequently, the amplitudes of the localized currents for the three bonds gradually differ as they come closer to the edge/corner as illustrated by the color legend below panel (b) in Fig.~\ref{fig:Fig2}.

\subsection{Quantitative estimates and analysis \--- case of armchair geometry \label{sec.7}}

In the previous section, we focused our quantitative analysis on a finite-size honeycomb lattice with zigzag edges as shown in Fig.~\ref{fig:Fig2}(a). However, such edge terminations are not the only possible options. Motivated by the graphene literature~\cite{RevModPhys.81.109,PhysRevB.82.045409,PhysRevB.54.17954}, we consider a different geometry with both the zigzag and the armchair edge terminations, as shown in Fig.~\ref{fig:Fig2}(d). Such geometry can be created by modifying the choice of the unit vectors~\cite{PhysRevB.76.205402} compared to the previous case. We again perform numerical diagonalization on this geometry of linear size $L = 21$ to obtain the current distribution, without any vison excitations. The localized current on the NNN bonds in the middle of the system saturates at the same value $\mathcal{I}^{(0)}_{\mathrm{bulk}}$ as before. There are four edges and two pairs of edges are related by the inversion operation.
Consequently, we focus only on two different edges $\mathrm{E}'_1$ and $\mathrm{E}'_2$ with zigzag and armchair termination, respectively. For $\mathrm{E}'_1$ (with $z$ and $y$-bonds), the current far away from the corners saturates at the same amplitude $\mathcal{I}^{(0)}_{\mathrm{edge}}$ as quoted in the previous section. Of course, this number varies as the bonds come close to the corners as before. This variation is similar to Fig.~\ref{fig:Fig2}(a).

However, the localized current distribution along the armchair edge $\mathrm{E}'_2$ (with all the three bonds \--- $x$, $y$ and $z$) is very different from its zigzag counterpart. The amplitude of the current on the outer lying $z$-bonds in the middle of the edge, away from the corners, saturates at a larger value $\mathcal{I}^{(0)}_{\mathrm{z-ach}} \sim 10.35$ nA, whereas the current on the accompanying $x$ and $y$-bonds saturates at a lower value $\mathcal{I}^{(0)}_{\mathrm{xy-ach}} \sim 4.31$ nA. The variation of the localized currents on the bonds as they come close to the corners for both the zigzag and the armchair edge terminations is illustrated in Fig.~\ref{fig:Fig2}(d). The current profile on the NNN bonds is similar to the previous case and is not shown for simplicity. A more interesting edge current profile can be obtained by modifying the edge termination patterns to others such as the bearded edge ribbons~\cite{PhysRevB.76.205402}.

\subsection{Quantitative estimates and analysis \--- case of two vison excitations \label{sec.8}}

In this section, we analyze the current distributions in the presence of localized vison excitations. For this purpose, we now consider a finite system with PBC and a bigger linear size $L = 32$ compared to the previous case. We chose a bigger system size to reduce the mutual interaction between the two vison excitations.

We introduce two visons and keep them at the largest possible distance $\lfloor \tfrac{L-1}{2} \rfloor$ to minimize their mutual interaction effects. Such configuration can be achieved by flipping all the link variables on the $z$-bonds which are connected by the string operator between them, as illustrated in Fig.~\ref{fig:Fig1}(a). However, as the string operator is gauge-dependent, the specific choice of the string operator connecting two vison excitations does not affect the electrical currents.

The ground state of the \textit{two-vison} Kitaev model, in the presence of an external magnetic field can be similarly obtained, utilizing the Majorana and gauge degrees of freedom, as $\ket{\Psi^{(2)}_{\mathrm{mag}}} = \mathcal{Z} \ket{\mathcal{G}_2} \otimes \ket{\mathcal{M}_{\mathrm{mag}}}$. Here $\ket{\mathcal{G}_2}$ corresponds to the gauge configuration with two visons, and $\mathcal{Z}$ is the projection operator as defined in Sec.~(\ref{sec.3}). We note that the action of the projection operator $\mathcal{Z}$ translates into satisfying a parity constraint for the matter and bond fermions as $(-1)^{n_f + n_\chi} = 1$ (see the SM~\cite{supp} for the technical details)~\cite{PhysRevB.92.014403,PhysRevB.84.165414}. Next, we compute the average localized current in the ground state $\ket{\Psi^{(2)}_{\mathrm{mag}}}$, by numerically diagonalizing the tight-binding Majorana Hamiltonian with $\ket{\mathcal{G}_2}$~\cite{supp}.

\begin{figure}[t]
\centering
\includegraphics[width=0.7\linewidth]{Fig3.png}
\caption{Current profile around a vison excitation (colored hexagon). In contrast to the previous case in Fig.~\ref{fig:Fig2}, the localized currents on the NN bonds around honeycomb plaquettes containing the vison are nonzero.}\label{fig:Fig3}
\end{figure}

In the previous section with the gauge configuration $\ket{\mathcal{G}_0}$, we observed that the localized currents on the NNN bonds, within the bulk, saturated at a uniform value $\mathcal{I}^{(0)}_{\mathrm{bulk}}$. Such current configuration led to vanishing localized currents on the NN bonds forming the hexagonal plaquettes. In the case with the gauge configuration $\ket{\mathcal{G}_2}$, the localized current profile far away from the static vison excitations is identical to $\mathcal{I}^{(0)}_{\mathrm{bulk}}$.

However, this distribution acquires drastic modification inside and around the vison excitations. The localized currents on the three NNN bonds ($x'$, $y'$, and $z'$-bonds) within the hexagon containing a vison excitation become isotropic but differ from the localized currents on the NNN bonds surrounding the honeycomb plaquettes, and consequently, finite non-vanishing currents emerge on the NN bonds forming the hexagon, as shown in Fig.~\ref{fig:Fig3}. The localized current on the NNN bonds inside the vison $\mathcal{I}^{(2)}_{\mathrm{vison}}$ saturate at a smaller value of $\sim 4.03$ nA, whereas the corresponding currents on the surrounding hexagons saturates at almost $\mathcal{I}^{(0)}_{\mathrm{bulk}} \sim 7.70$ nA. Consequently, NN bonds surrounding the hexagon containing the vison excitation host non-zero localized loop currents $\mathcal{I}^{(2)}_{\mathrm{vison-edge}} \sim 7.29$ nA, which is obtained by adding the contributions of the four triangular plaquettes containing the NN bonds. The variation of the amplitudes of these localized currents on the vison plaquettes as a function of the orientation and magnitude of the external magnetic field follows similar dependence as was illustrated in Fig.~\ref{fig:Fig2}(c) in Sec.~(\ref{sec.5}).

\section{Loop currents: Implications \label{sec.9}}

So far, we discussed patterns of the TRS breaking localized loop currents in a Kitaev system in the presence of an external magnetic field. For an infinite system, we noticed that the localized loop currents constitute a superlattice formed by inter-penetrating triangular lattices as shown in Fig.~\ref{fig:Fig2}(b$'$). The direction of the loop currents is determined by the sign of $\kappa$ and therefore depends on the direction of the external magnetic field. Therefore, in the absence of any vison excitations, one should expect them to contribute a uniform local magnetic moment $\mu_{\mathrm{loop}}$ for each of the hexagonal plaquettes. Since the current orientation for both the triangular plaquettes is the same, the magnetic moment of each of the two triangles adds up and leads to a total orbital magnetic moment in each hexagon as $\mu_{\mathrm{hexagon}} = 2 \times \mu_{\mathrm{loop}}$. Adopting realistic atomic lattice constant for $\alpha$-RuCl$_3$ from the recent \textit{ab intio} study~\cite{PhysRevB.93.155143}, we obtain $\mu_{\mathrm{hexagon}} \sim 0.0003 \mu_{\mathrm{B}}$, where $\mu_{\mathrm{B}}$ is the Bohr-magneton. Consequently, we find a weak emergent orbital ferromagnetic magnetic order within the Kitaev QSL phase. Naturally, the vison excitations would create a local change in this emergent long-range orbital ferromagnetic phase as the localized currents are different around the vison excitations. Such orbital ferromagnetic order can be resolved by polarized neutron diffraction~\cite{PhysRevLett.96.197001,Jeong2017}, muon spin spectroscopy~\cite{Zhang2018}, second harmonic generation~\cite{Zhao2016}, optical birefringence measurements~\cite{PhysRevLett.112.147001}, or superconducting quantum interference device magnetometers~\cite{PhysRevB.98.100403}.
The key feature of the Kitaev system is that for an external magnetic field only applied in $ab$-plane would lead to an out-of-plane (along $c$-axis). The direction and the magnitude of such $c$-axis magnetic field can be tuned by rotating the magnetic field in the $ab$-plane. Note that the direction of the orbital magnetic field depends on the sign of $\kappa$ which is dictated by Eq.~\eqref{eq.6}.

The magnitude of the localized moments is weak based on our quantitative estimate for $\alpha$-RuCl$_3$. However, two polymorphs of $\alpha$-RuCl$_3$, by replacing the ligand Cl$^-$-ion with Br$^-$ and I$^-$, have recently become possible candidate materials for realizing Kitaev QSL phase~\cite{Kaib2022}. Among these two polymorphs, RuI$_3$ is still debated to be a metal, although RuBr$_3$ is an insulator with almost half the Mott gap compared to RuCl$_3$~\cite{Kaib2022}. On the other hand, it has a much larger Ru-Ru distance and considerable Ru 4$d$ and Br 4$p$ orbital overlap~\cite{PhysRevB.105.L041112}. Adopting the parameters from the \textit{ab initio}~\cite{Ersan2019} results for RuBr$_3$, we predict an order of magnitude enhancement of such localized magnetic moments for an analogous situation in RuBr$_3$ with $\mu_{\mathrm{hexagon}} \sim 0.003 \mu_{\mathrm{B}}$  at $\kappa \sim 0.14 K$ if the Kitaev quantum spin liquid phase can be realized there. On the other hand, if we adopt an inflated magnitude of $\kappa$ for RuBr$_3$ motivated by the experimentally measured Majorana gap for $\alpha$-RuCl$_3$~\cite{Tanaka2022} as $\kappa \sim 0.77$, the orbital magnetization can go up as much as $0.01\mu_{\mathrm{B}}$.

\section{Discussion and conclusion \label{sec.10}}

In this paper, we developed a theoretical framework for analyzing the remnant electrical current responses in spin-orbit coupled multiorbital Mott insulators, such as iridates or ruthenates, considered to be plausible candidate materials to realize the Kitaev QSL phase. Our analysis and predictions reveal an exciting and alternative detection protocol for the signatures of the fractionalized excitations in a Kitaev QSL through emergent orbital magnetization due to the localized orbital currents. This formalism can be generalized to other QSL hosting Mott insulators with broken time-reversal symmetry. Starting from the microscopic multiorbital Hubbard-Kanamori model, we derived the functional form of the induced localized currents in a triangular loop [see Eq.~(\ref{eq.4.1})] using SWT. The current operator expressed in terms of the spin operator in Mott insulators without SU(2) spin rotation symmetry is an important result, on its own, and forms the basis for the rest of our paper. In the presence of an external magnetic field, we solve the Kitaev model in the Majorana representation and notice that the ground state hosts a non-zero expectation value for the localized loop currents. Because of the underlying honeycomb lattice structure, the average non-zero triangular loop current leads to a rich current profile in the two-dimensional system, as shown in Fig.~\ref{fig:Fig2}. The current on the NNN bonds on the hexagonal plaquettes extends over the entire system and forms a super-lattice of two inter-penetrating triangular lattices as illustrated in Fig.~\ref{fig:Fig2}(b$'$), whereas the current on the NN bonds vanishes identically within the bulk. However, depending on the edge termination profiles (zigzag, armchair, or bearded) of the 2D Kitaev system, the localized currents on the NN bonds along the edge can survive and lead to rich and interesting patterns.

We then revisited the analysis of our current profile in the presence of vison excitations. While the localized current distribution in the absence of any visons (uniform gauge configuration) remains identical on the NNN bonds within the bulk, leading to the vanishing currents on NN bonds, a single vison excitation leads to a drastic modification of the currents on the bonds containing the $\mathbb{Z}_2$ vortex (see Fig.~\ref{fig:Fig3}). Depending on the orientation of the applied magnetic field, the currents on NNN bonds in hexagons containing the $\mathbb{Z}_2$ vortex decrease compared to the surrounding hexagons away from the vison excitation. This leads to a reappearance of non-zero localized currents on the NN bonds around the hexagon with vison excitation. This implies that the visons carry the physical magnetic flux, which makes visons resemble the Abrikosov vortices in superconductors. The visons can couple to an external magnetic field directly and can interact through electromagnetic coupling. This direct electromagnetic coupling is generally weak because of the small magnetic moment associated with visons. Under the circumstances when the vison gap is small, one may induce vison lattice by an external magnetic field similar to the Abrikosov vortex lattice. \\

\section{Acknowledgments}

We would like to thank Philipp Gegenwart, Alexander V. Balatsky, Cristian Batista, Dieter Vollhardt, Krishnendu Sengupta, Ross David Mcdonald, Vivien Zapf, Avadh Saxena, and Umesh Kumar for useful discussions and careful remarks on the manuscript. This work was carried out under the auspices of the US DOE NNSA under Contract No. 89233218CNA000001 through the LDRD Program.

\section{Schrieffer-Wolff transformation \label{sec.1}}

In this section, we lay out the details of the Schrieffer-Wolff transformation (SWT) to derive the low-energy effective Hamiltonian from a generic strongly correlated electronic model~\cite{PhysRevB.105.L180414,Eckstein}. The original correlated Hamiltonian is written as
\begin{equation}\label{eq.1}
\mathcal{H} = \mathcal{H}_0 + \mathcal{H}_1,
\end{equation}
where $\mathcal{H}_0$ denotes the correlated part of the model (with parameters such as onsite Coulomb repulsion $U$, Hund's coupling $J_\mathrm{H}$, etc.), and $\mathcal{H}_1$ denotes the tight-binding (TB) contribution (with parameters including all the hopping parameters $t_1$, $t_2$ or spin-orbit coupling, etc.). Here, various $t_i$'s correspond to the hopping amplitudes between multiple orbitals or sites. In a strongly correlated system, interaction strengths are naturally much larger than the hopping parameters. In this case, we can reduce the full Hamiltonian in Eq.~(\ref{eq.1}) to a low-energy effective description to explain various properties of the parent system. SWT is a very powerful tool in this regard, which is achieved via a unitary transformation $U = e^{i\ensuremath{\mathcal{S}}}$, where $\ensuremath{\mathcal{S}}$ is a hermitian operator. After applying such a unitary transformation, we obtain the rotated Hamiltonian as
\begin{align}\label{eq.2}
\nonumber
\mathcal{H}'
&
=
U^{\dagger} \mathcal{H} U = e^{i\ensuremath{\mathcal{S}}} \mathcal{H}e^{-i\ensuremath{\mathcal{S}}} \\
\nonumber
&
=
\mathcal{H} + \big[i\ensuremath{\mathcal{S}}, \mathcal{H} \big] + \frac{1}{2!} \big[i\ensuremath{\mathcal{S}}, \big[ i\ensuremath{\mathcal{S}},\mathcal{H} \big]\big] + \frac{1}{3!} \big[i\ensuremath{\mathcal{S}}, \big[i\ensuremath{\mathcal{S}}, \big[i\ensuremath{\mathcal{S}}, \mathcal{H} \big] \big] \big] + \cdots \\
&
=
\mathcal{H}_0 + \mathcal{H}_1 +
i\big[\ensuremath{\mathcal{S}}, \mathcal{H}_0] + i\big[\ensuremath{\mathcal{S}}, \mathcal{H}_1] -
\frac{1}{2}\big[\ensuremath{\mathcal{S}}, \big[ \ensuremath{\mathcal{S}},\mathcal{H}_0 \big]\big] - \frac{1}{2} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}},\mathcal{H}_1 \big]\big] -
\frac{i}{3!} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \mathcal{H}_0 \big] \big] \big] - \frac{i}{3!} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \mathcal{H}_1 \big] \big] \big] + \cdots,
\end{align}
where $\ensuremath{\mathcal{S}}$ is to be determined later. The idea is now to evaluate the generating function $\ensuremath{\mathcal{S}}$ in a perturbative manner in terms of the small TB parameter, such that an effective Hamiltonian in each order acts on a truncated Hilbert space avoiding any doubly occupied or empty sites (which corresponds to high-energy configurations). Note that we consider half-filling. For this purpose, we introduce two projection operators $\mathcal{P}$, and $\mathcal{Q}$, where $\mathcal{P}$ projects to singly occupied sites and $\mathcal{Q} = 1 - \mathcal{P}$, correspondingly, projects to either a doubly occupied or an empty site. We now write the generating function $\ensuremath{\mathcal{S}}$ as a perturbative expansion in terms of the TB parameters as
\begin{equation}\label{eq.3}
\ensuremath{\mathcal{S}} = \ensuremath{\mathcal{S}}^{(1)} + \ensuremath{\mathcal{S}}^{(2)} + \ensuremath{\mathcal{S}}^{(3)} + \ensuremath{\mathcal{S}}^{(4)} + \cdots.
\end{equation}
Next, collecting terms of the same order in perturbation, we obtain
\begin{subequations}
\begin{align}
\label{eq.4.1}
&
\mathcal{H}'
=
\mathcal{H}_0 + \mathcal{H}_1 + i \big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_0 \big] \\
\label{eq.4.2}
&
+
i \big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_1 \big] + i\big[ \ensuremath{\mathcal{S}}^{(2)}, \mathcal{H}_0 \big] - \frac{1}{2} [\ensuremath{\mathcal{S}}^{(1)}, \big[ \ensuremath{\mathcal{S}}^{(1)}(t), \mathcal{H}_0 \big]\big]  \\
\label{eq.4.3}
&
+
i \big[ \ensuremath{\mathcal{S}}^{(3)}, \mathcal{H}_0 \big] + i\big[ \ensuremath{\mathcal{S}}^{(2)}, \mathcal{H}_1 \big] - \frac{1}{2} \big[ \ensuremath{\mathcal{S}}^{(1)}, \big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_1 \big] + \big[ \ensuremath{\mathcal{S}}^{(2)} , \mathcal{H}_0 \big]\big] - \frac{1}{2} \big[\ensuremath{\mathcal{S}}^{(2)}, \big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_0 \big]\big] - \frac{i}{3!} \big[ \ensuremath{\mathcal{S}}^{(1)}, \big[ \ensuremath{\mathcal{S}}^{(1)}, \big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_0 \big]\big]\big] \\
\label{eq.4.4}
& + \mathcal{O}(4) \; \mathrm{terms},
\end{align}
\end{subequations}
where, in each line, we arranged terms of the same order in perturbation.  The rotated Hamiltonian in Eq.~(\ref{eq.4.1}-\ref{eq.4.4}) can be written in a compact form as follows
\begin{equation}\label{eq.5}
\mathcal{H}' = \sum_{m=0}^{n} \mathcal{H}^{(m)}_{\mathrm{eff}} + \mathcal{O}(n+1),
\end{equation}
where the effective Hamiltonian in $m$-th perturbative order $\mathcal{H}^{(m)}_{\mathrm{eff}}$ is computed in such a way that it does not have any mixing terms, \textit{i.e.}, $\mathcal{P} \mathcal{H}^{(m)}_{\mathrm{eff}} \mathcal{Q}$ = $\mathcal{Q} \mathcal{H}^{(m)}_{\mathrm{eff}} \mathcal{P}$ = 0. In the subsequent analysis, we provide the key steps to obtain the analytical structure of the effective Hamiltonian up to third-order in the perturbation expansion.

\subsection{Second-order effective Hamiltonian \label{sec.2}}

We now move on to the analysis of the second-order effective Hamiltonian and obtain the equation for the generating function $\ensuremath{\mathcal{S}}^{(1)}$ from Eq.~(\ref{eq.4.1}). Utilizing this relation in Eq.~(\ref{eq.4.2}), we obtain the effective Hamiltonian in the second-order perturbation as
\begin{subequations}
\begin{align}
\label{eq.6.1}
& \big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_0 \big]  = i\mathcal{H}_1, \\
\label{eq.6.2}
& \mathcal{H}^{(2)}_{\mathrm{eff}}  	= \frac{i}{2} \big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_1 \big].
\end{align}
\end{subequations}
To solve for the generating function $\ensuremath{\mathcal{S}}^{(1)}$, we first notice its $2\times 2$ matrix structure in the basis of the projection operators $\mathcal{P}$, and $\mathcal{Q}$ as
\begin{equation}\label{eq.7}
\ensuremath{\mathcal{S}}^{(1)}  = \begin{pmatrix}
\mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{P} & \mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q} \\
\mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{P} & \mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q}
\end{pmatrix},
\end{equation}
where $\mathcal{P}$, and $\mathcal{Q}$ are defined earlier. Since, we need to remove the off-diagonal elements, the diagonal elements [$\mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{P}$, $\mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q}$] are naturally assumed to be zero, and the off-diagonal elements are obtained from Eq.~(\ref{eq.6.1}) as
\begin{subequations}
\begin{align}
& \mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{H}_0\mathcal{Q} -  \mathcal{P} \mathcal{H}_0 \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q} = i\mathcal{P}\mathcal{H}_1 \mathcal{Q} \nonumber \\
\label{eq.8.1}
& \Rightarrow \mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q} \mathcal{Q} \mathcal{H}_0 \mathcal{Q} - \cancel{ \mathcal{P} \mathcal{H}_0 \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q} } = i \mathcal{P}\mathcal{H}_1 \mathcal{Q}, \\
& \mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{H}_0\mathcal{P} - \mathcal{Q} \mathcal{H}_0 \ensuremath{\mathcal{S}}^{(1)} \mathcal{P} = i \mathcal{Q}\mathcal{H}_1 \mathcal{P} \nonumber \\
\label{eq.8.2}
& \Rightarrow \cancel{ \mathcal{Q} S^{(1)}(t) \mathcal{Q} \mathcal{Q} \mathcal{H}_0\mathcal{P}} -  \mathcal{Q} \mathcal{H}_0 \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{P} =  i\mathcal{Q}\mathcal{H}_1 \mathcal{P}.
\end{align}
\end{subequations}

\subsection{Third-order effective Hamiltonian \label{sec.3}}

The third-order effective Hamiltonian can be obtained similarly as outlined in Sec.~(\ref{sec.2}). We first obtain the equation for the generating function $\ensuremath{\mathcal{S}}^{(2)}$ [see Eq.~(\ref{eq.4.2})], and utilize it to derive the effective Hamiltonian in the third-order perturbation. Consequently, we obtain [Note: the last term in Eq.~(\ref{eq.4.2}) can be recasted in the form as in Eq.~(\ref{eq.6.2}) by utilizing the equation of motion for $\ensuremath{\mathcal{S}}^{(1)}$ in Eq.~(\ref{eq.6.1})]
\begin{subequations}
\begin{align}
\label{eq.9.1}
& \big[ \ensuremath{\mathcal{S}}^{(2)}, \mathcal{H}_0 \big]  =  -\big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_1 \big] \\
\label{eq.9.2}
& \mathcal{H}^{(3)}_{\mathrm{eff}} 		= \frac{i}{2}\big[ \ensuremath{\mathcal{S}}^{(2)}, \mathcal{H}_1 \big] + \frac{1}{6} \big[ \ensuremath{\mathcal{S}}^{(1)}, \big[ \ensuremath{\mathcal{S}}^{(1)}, \mathcal{H}_1 \big].
\end{align}
\end{subequations}
Finally, we evaluate the matrix elements for $\ensuremath{\mathcal{S}}^{(2)}$ in the projection operator basis as
\begin{subequations}
\begin{align}
\nonumber
\mathcal{P}\ensuremath{\mathcal{S}}^{(2)} \mathcal{Q} \mathcal{Q} \mathcal{H}_0 \mathcal{Q}
-
\cancel{\mathcal{P} \mathcal{H}_0 \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(2)} \mathcal{Q}}
&
=
\cancel{\mathcal{P} \mathcal{H}_1 \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q}} - \mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q} \mathcal{Q} \mathcal{H}_1 \mathcal{Q} \\
\label{eq.10.1}
\Rightarrow \mathcal{P}\ensuremath{\mathcal{S}}^{(2)} \mathcal{Q} \mathcal{Q} \mathcal{H}_0 \mathcal{Q}
&
=
-
\mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q} \mathcal{Q} \mathcal{H}_1 \mathcal{Q}, \\
\nonumber
\cancel{\mathcal{Q}\ensuremath{\mathcal{S}}^{(2)} \mathcal{Q} \mathcal{Q} \mathcal{H}_0 \mathcal{P}}
-
\mathcal{Q} \mathcal{H}_0 \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(2)} \mathcal{P}
&
=
\mathcal{Q} \mathcal{H}_1 \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{P} - \cancel{\mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q} \mathcal{Q} \mathcal{H}_1 \mathcal{P}} \\
\label{eq.10.2}
\Rightarrow  \mathcal{Q} \mathcal{H}_0 \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(2)} \mathcal{P}
&
=
-
\mathcal{Q} \mathcal{H}_1 \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{P},
\end{align}
\end{subequations}
where $\mathcal{Q} \mathcal{H}_1 \mathcal{Q}$ corresponds to hopping between either two doubly occupied states or two singly occupied sites (always leaving one empty site after the hopping, in the latter case).

\section{Effective analytical structure of a generic operator \label{sec.4}}

In previous Sec.~(\ref{sec.1}), we provided all the details of the SWT to obtain the generic forms of the effective low-energy Hamiltonian. Once the generating functions $\ensuremath{\mathcal{S}}^{(m)}$ are obtained upto $m$-th order in perturbation, any other operators can be rotated in the same manner as in Eq.~(\ref{eq.2}). Consequently, in the rotated frame, any local operator $\mathcal{O}_i$ becomes
\begin{equation}\label{eq.11}
\tilde{\mathcal{O}}_i = e^{i\ensuremath{\mathcal{S}}} \mathcal{O}_i e^{-i\ensuremath{\mathcal{S}}} = \mathcal{O}_i + i\big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] - \frac{1}{2!} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] \big] - \frac{i}{3!} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] \big] \big] + \cdots,
\end{equation}
where $\mathcal{O}_i$ a generic physical operator \textit{viz.}, the local charge imbalance operator at half-filling \textit{i.e.} $\delta \bm{\rho}_{i} = e \left(d^{\dagger}_{i \alpha \sigma}d_{i \alpha \sigma} - 1 \right)$, or the current operator \textit{i.e.} $\bm{\mathcal{I}}_{ij} = \tfrac{iet' \hat{\ensuremath{\mathbf{r}}}_{ij}}{\hbar} \sum_{\alpha \beta}\left( d^{\dagger}_{j\alpha \sigma} d_{i\beta \sigma} - d^{\dagger}_{i\beta\sigma}d_{j\alpha\sigma} \right)$, where the sum over repeated indices are assumed and $(\alpha,\beta)$, and $\sigma$ correspond to the orbital, and spin degrees of freedom, respectively. Utilizing the projection operators $\mathcal{P}$ and $\mathcal{Q}$, as defined earlier, we analyze the low-energy effective structure for $\mathcal{O}_i$ as
\begin{align}
\nonumber
\mathcal{P} \tilde{\mathcal{O}}_i \mathcal{P}
&
=
\cancel{\mathcal{P} \mathcal{O}_i \mathcal{P}}
+
i \cancel{\mathcal{P} \big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] \mathcal{P}}
-
\frac{1}{2!} \mathcal{P} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] \big] \mathcal{P}
-
\frac{i}{3!} \mathcal{P} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] \big] \big] \mathcal{P} + \cdots  \\
\label{eq.12}
&
=
-\frac{1}{2!} \mathcal{P} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] \big] \mathcal{P}
-
\frac{i}{3!} \mathcal{P} \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] \big] \big] \mathcal{P} + \cdots,
\end{align}
where we utilized the relation $\mathcal{P} \mathcal{O}_i\mathcal{P} = 0$ and $\mathcal{P} \big[\ensuremath{\mathcal{S}}, \mathcal{O}_i \big] \mathcal{P} = 0$. This property holds for the case where $\mathcal{O}_i$ is a diagonal operator in the projection operator space. An example for this is the charge imbalance operator $\delta \bm{\rho}_i$~\cite{PhysRevLett.125.227202}. However, for the current operator $\bm{\mathcal{I}}_{ij}$ (off-diagonal in $\mathcal{P}$-$\mathcal{Q}$ space), such a restriction does not hold \textit{i.e.} $\mathcal{P} \big[\ensuremath{\mathcal{S}}, \bm{\mathcal{I}}_{ij} \big] \mathcal{P} \neq 0$. We now focus our analysis on the low-energy effective form for the current operator.

\subsection{Circulating loop current \label{sec.5}}

Since there are no mobile-charged carriers in a Mott insulator, a free-flowing current cannot exist in the system. However, this constraint does not exclude the possibility of having an induced circulating loop current in a closed loop. A minimal closed loop is made out of three neighboring sites forming a triangular plaquette. Here, we derive the localized loop current operator up to third-order in perturbation. Formally it reads as
\begin{align}
\label{eq.13}
\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}
&
=
i\mathcal{P}\big[ \ensuremath{\mathcal{S}}^{(2)}, \bm{\mathcal{I}}_{ij} \big] \mathcal{P}
-
\frac{1}{2} \mathcal{P} \big[ \ensuremath{\mathcal{S}}^{(1)}, \big[ \ensuremath{\mathcal{S}}^{(1)}, \bm{\mathcal{I}}_{ij} \big] \big] \mathcal{P}
=
i\mathcal{P} \ensuremath{\mathcal{S}}^{(2)} \mathcal{Q} \mathcal{Q} \bm{\mathcal{I}}_{ij} \mathcal{P}
+
\mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q} \mathcal{Q} \bm{\mathcal{I}}_{ij} \mathcal{Q} \mathcal{Q} \ensuremath{\mathcal{S}}^{(1)} \mathcal{P}
+
\mathrm{h.c.},
\end{align}
where we utilized the equation for the generating functions $\ensuremath{\mathcal{S}}^{(1)}$ and $\ensuremath{\mathcal{S}}^{(2)}$ as in Eqs.~(\ref{eq.8.1}-\ref{eq.8.2}) and Eqs.~(\ref{eq.10.1}-\ref{eq.10.2}). The form of the current operator in a single-orbital Hubbard model~\cite{PhysRevB.78.024402} can be obtained utilizing the above scheme. Here, we apply this scheme to the multiorbital Hubbard-Kanamori Hamiltonian, relevant for $\alpha$-RuCl$_3$~\cite{Kanamori1957_1,Kanamori1957_2}, to derive the low-energy form of the current operator.

\subsection{Hubbard-Kanamori model \label{sec.6}}

The microscopic multiorbital Hamiltonian for Kiteav magnets such as $\alpha$-RuCl$_3$ is written as (see Refs.~[\onlinecite{PhysRevLett.105.027204}],[\onlinecite{PhysRevLett.112.077204}])
\begin{subequations}
\begin{align}
\label{eq.14.1}
\mathcal{H} & = \mathcal{H}_0 + \mathcal{H}_1 + \mathcal{H}_2, \\
\label{eq.14.2}
\mathcal{H}_0 &  =   U\sum_{i\alpha} n_{i\alpha,\uparrow}n_{i\alpha,\downarrow} + \frac{U'}{2}\sum_{\substack{\alpha\neq \beta, \\ \sigma, \sigma'}} n_{i\alpha \sigma} n_{i\beta \sigma'}  -\frac{J_\mathrm{H}}{2}\sum_{\substack{\alpha\neq \beta, \\ \sigma, \sigma'}} d^{\dagger}_{i\alpha \sigma} d_{i\alpha\sigma'} d^{\dagger}_{i\beta \sigma'}  d_{i\beta \sigma}^{\phantom\dagger} + \frac{\lambda}{2} \sum_{i} d^{\dagger}_i \left( \mathbf{L} \cdot \mathbf{S} \right) d_{i}, \\
\label{eq.14.3}
\mathcal{H}_1 & = \sum_{\langle ij \rangle \sigma} \begin{pmatrix}
d^{\dagger}_{i xz \sigma} & d^{\dagger}_{i {yz} \sigma} & d^{\dagger}_{i {xy} \sigma}
\end{pmatrix} \begin{pmatrix}
0 & t_2 & 0 \\
t_2 & 0 & 0 \\
0 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
d_{j {xz} \sigma} \\
d_{j {yz} \sigma} \\
d_{j {xy} \sigma}
\end{pmatrix}, \quad
\mathcal{H}_2  = \sum_{\llangle ij \rrangle \sigma} \begin{pmatrix}
d^{\dagger}_{i xz \sigma} & d^{\dagger}_{i {yz} \sigma} & d^{\dagger}_{i {xy} \sigma}
\end{pmatrix} \begin{pmatrix}
0 & t'_2 & 0 \\
t'_2 & 0 & 0 \\
0 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
d_{j {xz} \sigma} \\
d_{j {yz} \sigma} \\
d_{j {xy} \sigma}
\end{pmatrix},
\end{align}
\end{subequations}
where ($\mathcal{H}_2$) $\mathcal{H}_1$ in Eq.~(\ref{eq.14.3}) is the (next)-nearest-neighbor hopping matrix between the $d_{xz}$, $d_{yz}$ and $d_{xy}$ orbitals in the $t_{2\mathsf{g}}$ manifold. The specific tight-binding structure of $\mathcal{H}_2$ is adopted from  Ref.~[\onlinecite{PhysRevB.91.241110}]. For the Kanamori part in Eq.~(\ref{eq.14.2}) $U$ is the strength of the onsite Coulomb repulsion, $J_{\mathrm{H}}$ is the Hund's coupling, $\lambda$ is the atomic spin-orbit coupling. For rotational invariant system, we consider $U' = U -2 J_{\mathrm{H}}$~\cite{PhysRevB.105.L180414,Kumar2022}. Since each individual hopping between the neighboring magnetic ions is oriented along the particular type of the bond [\textit{see Fig. 1(b) in the main text}], we need to modify the above Slater-Koster TB model depending on the particular planar orientation of the bonds. In Eq.~(\ref{eq.14.3}), we assumed a generic notation and wrote both the nearest and next-nearest neighbor hopping assuming the corresponding bond to be lying in the $xy$-plane ($\mathbf{z}$-bond). However, in the actual crystalline environment, the three individual hopping matrix elements (along three different bonds-$\mathbf{x}, \mathbf{y}, \mathbf{z}$) are modified as (for the next-neighbor bonds we label the bond types as $\mathbf{x}', \mathbf{y}', \mathbf{z}'$)
\begin{subequations}
\begin{align}
\label{eq.15.1}
\mathcal{H}_1 & = \sum_{\substack{\langle ij \rangle_{\mathbf{z}} \\ \sigma}} \begin{pmatrix}
d^{\dagger}_{i xz \sigma} & d^{\dagger}_{i {yz} \sigma} & d^{\dagger}_{i {xy} \sigma}
\end{pmatrix} \begin{pmatrix}
0 & t_2 & 0 \\
t_2 & 0 & 0 \\
0 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
d_{j {xz} \sigma} \\
d_{j {yz} \sigma} \\
d_{j {xy} \sigma}
\end{pmatrix}, \\
\label{eq.15.2}
\mathcal{H}_1 & = \sum_{\substack{\langle jk \rangle_{\mathbf{y}} \\ \sigma}} \begin{pmatrix}
d^{\dagger}_{j xz \sigma} & d^{\dagger}_{j {yz} \sigma} & d^{\dagger}_{j {xy} \sigma}
\end{pmatrix} \begin{pmatrix}
0 & 0 & 0   \\
0 & 0 & t_2 \\
0 & t_2 & 0
\end{pmatrix}
\begin{pmatrix}
d_{k {xz} \sigma} \\
d_{k {yz} \sigma} \\
d_{k {xy} \sigma}
\end{pmatrix}, \\
\label{eq.15.3}
\mathcal{H}_1 & = \sum_{\substack{\langle ki \rangle_{\mathbf{x}} \\ \sigma}} \begin{pmatrix}
d^{\dagger}_{k xz \sigma} & d^{\dagger}_{k {yz} \sigma} & d^{\dagger}_{k {xy} \sigma}
\end{pmatrix} \begin{pmatrix}
0 & 0 & t_2 \\
0 & 0 & 0 \\
t_2 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
d_{i {xz} \sigma} \\
d_{i {yz} \sigma} \\
d_{i {xy} \sigma}
\end{pmatrix}, \\
\label{eq.15.4}
\mathcal{H}_2 & = \mathcal{H}_1[\langle ij \rangle \rightarrow \llangle ij \rrangle, \;\; t_2 \rightarrow t_2', \;\; (\mathbf{x}, \mathbf{y}, \mathbf{z}) \rightarrow (\mathbf{x}', \mathbf{y}', \mathbf{z}')],
\end{align}
\end{subequations}
where the subscripts for bond-type are chosen according to Fig.~1(b) in the main text. For subsequent analysis, we adopt all the parameters entering Eqs.~(\ref{eq.14.1}-\ref{eq.15.4}) from the recent \textit{ab initio}~\cite{PhysRevB.93.155143} and \textit{photoemission studies}~\cite{Sinn2016} for $\alpha$-RuCl$_3$ as: $U = 3.0$ eV, $J_{\mathrm{H}} = 0.45$ eV, $t_2 = 0.191$ eV, and $t'_2 = -0.058$ eV. For subsequent analysis, we first rewrite the Hamiltonian in Eq.~(\ref{eq.14.1}) in the irreducible representation of the doubly occupied states of the octahedral point group (O$_\mathrm{h}$) as
\begin{equation}\label{eq.16}
\mathcal{H}_0 = \sum_{i} \sum_{\Gamma} \sum_{\mathsf{g}_{\Gamma}} U_{\Gamma} \ket{i;\Gamma, \mathsf{g}_{\Gamma}} \bra{i;\Gamma, \mathsf{g}_{\Gamma}},
\end{equation}
where $\Gamma$ corresponds to a particular irreducible representation and $\mathsf{g}_{\Gamma}$ characterizes its degeneracy. The energy of the three non-degenerate states are given as follows~\cite{PhysRevB.65.064442,PhysRevB.94.174416,PhysRevB.105.L180414}
\begin{subequations}
\begin{align}
\label{eq.17.1}
& U_{\mathrm{A}_1}  = U + 2 J_{\mathrm{H}}, \\
\label{eq.17.2}
& U_{\mathrm{E}}  = U - J_{\mathrm{H}}, \\
\label{eq.17.3}
& U_{\mathrm{T}_1}  = U - 3 J_{\mathrm{H}}, \\
\label{eq.17.4}
& U_{\mathrm{T}_2}  = U -  J_{\mathrm{H}}.
\end{align}
\end{subequations}
There are three orbitals and two spin degrees of freedom, and we have to put two electrons within this manifold. Hence, there are $^6C_2 = 15$ possibility of doubly occupied states. From the character table of O$_\mathrm{h}$, we write these fifteen intermediate doubly occupied states at site $i$ as (spanned by the Hilbert space of the projection operator $\mathcal{Q}$)
\begin{subequations}
\begin{align}
\label{eq.18.1}
& \ket{i;\mathrm{A}_1} = \frac{1}{\sqrt{3}} (d^{\dagger}_{ixz \uparrow}d^{\dagger}_{ixz \downarrow} + d^{\dagger}_{iyz \uparrow}d^{\dagger}_{iyz \downarrow} + d^{\dagger}_{ixy \uparrow}d^{\dagger}_{ixy \downarrow}) \ket{0}, \\
\label{eq.18.2}
& \ket{i;\mathrm{E}, \mathrm{u}}  = \frac{1}{\sqrt{6}} (d^{\dagger}_{iyz \uparrow}d^{\dagger}_{iyz \downarrow} + d^{\dagger}_{ixz \uparrow}d^{\dagger}_{ixz \downarrow} -2 d^{\dagger}_{ixy \uparrow}d^{\dagger}_{ixy \downarrow}) \ket{0}, \\
\label{eq.18.3}
& \ket{i;\mathrm{E}, \mathrm{v}}  = \frac{1}{\sqrt{2}} (d^{\dagger}_{iyz \uparrow}d^{\dagger}_{iyz \downarrow} - d^{\dagger}_{ixz \uparrow}d^{\dagger}_{ixz \downarrow}) \ket{0}, \\
\label{eq.18.4}
& \ket{i;\mathrm{T}_1,\alpha_+} = d^{\dagger}_{iyz \uparrow} d^{\dagger}_{i zx \uparrow} \ket{0}, \\
\label{eq.18.5}
& \ket{i;\mathrm{T}_1,\alpha_-} = d^{\dagger}_{iyz \downarrow} d^{\dagger}_{i zx \downarrow} \ket{0}, \\
\label{eq.18.6}
& \ket{i;\mathrm{T}_1,\alpha} = \frac{1}{\sqrt{2}} ( d^{\dagger}_{iyz \uparrow} d^{\dagger}_{i zx \downarrow} + d^{\dagger}_{iyz \downarrow} d^{\dagger}_{i zx \uparrow} ) \ket{0}, \\
\label{eq.18.7}
& \ket{i;\mathrm{T}_1,\beta_+} = d^{\dagger}_{izx \uparrow} d^{\dagger}_{i xy \uparrow} \ket{0}, \\
\label{eq.18.8}
& \ket{i;\mathrm{T}_1,\beta_-} = d^{\dagger}_{izx \downarrow} d^{\dagger}_{i xy \downarrow} \ket{0}, \\
\label{eq.18.9}
& \ket{i;\mathrm{T}_1,\beta} = \frac{1}{\sqrt{2}} ( d^{\dagger}_{izx \uparrow} d^{\dagger}_{i xy \downarrow} + d^{\dagger}_{izx \downarrow} d^{\dagger}_{i xy \uparrow} ) \ket{0}, \\
\label{eq.18.10}
& \ket{i;\mathrm{T}_1,\gamma_+} = d^{\dagger}_{ixy \uparrow} d^{\dagger}_{i yz \uparrow} \ket{0}, \\
\label{eq.18.11}
& \ket{i;\mathrm{T}_1,\gamma_-} = d^{\dagger}_{ixy \downarrow} d^{\dagger}_{i yz \downarrow} \ket{0}, \\
\label{eq.18.12}
& \ket{i;\mathrm{T}_1,\gamma} = \frac{1}{\sqrt{2}} ( d^{\dagger}_{ixy \uparrow} d^{\dagger}_{i yz \downarrow} + d^{\dagger}_{ixy \downarrow} d^{\dagger}_{i yz \uparrow} ) \ket{0}, \\
\label{eq.18.13}
& \ket{i;\mathrm{T}_2,\alpha} = \frac{1}{\sqrt{2}} ( d^{\dagger}_{iyz \uparrow} d^{\dagger}_{i zx \downarrow} - d^{\dagger}_{iyz \downarrow} d^{\dagger}_{i zx \uparrow} ) \ket{0}, \\
\label{eq.18.14}
& \ket{i;\mathrm{T}_2,\beta} = \frac{1}{\sqrt{2}} ( d^{\dagger}_{izx \uparrow} d^{\dagger}_{i xy \downarrow} - d^{\dagger}_{izx \downarrow} d^{\dagger}_{i xy \uparrow} ) \ket{0}, \\
\label{eq.18.15}
& \ket{i;\mathrm{T}_2,\gamma} = \frac{1}{\sqrt{2}} ( d^{\dagger}_{ixy \uparrow} d^{\dagger}_{i yz \downarrow} - d^{\dagger}_{ixy \downarrow} d^{\dagger}_{i yz \uparrow} ) \ket{0}.
\end{align}
\end{subequations}
The singly occupied states at site $i$ are written as (spanned by the Hilbert space of the projection operator $\mathcal{P}$)
\begin{subequations}
\begin{align}
\label{eq.19.1}
\ket{i,+} & = \frac{1}{\sqrt{3}} \left( id^{\dagger}_{ixz\downarrow} + d^{\dagger}_{iyz\downarrow} + d^{\dagger}_{ixy\uparrow}\right) \ket{0}, \\
\label{eq.19.2}
\ket{i,-} & = \frac{1}{\sqrt{3}} \left( id^{\dagger}_{ixz\uparrow} - d^{\dagger}_{iyz\uparrow} + d^{\dagger}_{ixy\downarrow}\right) \ket{0}.
\end{align}
\end{subequations}
Consequently, the low-energy Hilbert space for three site problem is written in terms of the eight states ($2^3$) as (considering a three-site triangular subsystem within the crystal)
\begin{subequations}
\begin{align}
\label{eq.20.1}
\ket{+,+,+} & = \frac{1}{3\sqrt{3}} \left( id^{\dagger}_{ixz\downarrow} + d^{\dagger}_{iyz\downarrow} + d^{\dagger}_{ixy\uparrow} \right) \left( id^{\dagger}_{jxz\downarrow} + d^{\dagger}_{jyz\downarrow} + d^{\dagger}_{jxy\uparrow} \right) \left( id^{\dagger}_{kxz\downarrow} + d^{\dagger}_{kyz\downarrow} + d^{\dagger}_{kxy\uparrow} \right) \ket{0}, \\
\label{eq.20.2}
\ket{+,+,-} & = \frac{1}{3\sqrt{3}} \left( id^{\dagger}_{ixz\downarrow} + d^{\dagger}_{iyz\downarrow} + d^{\dagger}_{ixy\uparrow} \right) \left( id^{\dagger}_{jxz\uparrow} + d^{\dagger}_{jyz\uparrow} + d^{\dagger}_{jxy\downarrow} \right) \left( id^{\dagger}_{kxz\downarrow} - d^{\dagger}_{kyz\downarrow} + d^{\dagger}_{kxy\uparrow} \right) \ket{0},  \\
\label{eq.20.3}
\ket{+,-,+} & = \frac{1}{3\sqrt{3}} \left( id^{\dagger}_{ixz\uparrow} + d^{\dagger}_{iyz\uparrow} + d^{\dagger}_{ixy\downarrow} \right) \left( id^{\dagger}_{jxz\downarrow} - d^{\dagger}_{jyz\downarrow} + d^{\dagger}_{jxy\uparrow} \right) \left( id^{\dagger}_{kxz\downarrow} + d^{\dagger}_{kyz\downarrow} + d^{\dagger}_{kxy\uparrow} \right) \ket{0},  \\
\label{eq.20.4}
\ket{+,-,-} & = \frac{1}{3\sqrt{3}} \left( id^{\dagger}_{ixz\uparrow} + d^{\dagger}_{iyz\uparrow} + d^{\dagger}_{ixy\downarrow} \right) \left( id^{\dagger}_{jxz\uparrow} - d^{\dagger}_{jyz\uparrow} + d^{\dagger}_{jxy\downarrow} \right) \left( id^{\dagger}_{kxz\downarrow} - d^{\dagger}_{kyz\downarrow} + d^{\dagger}_{kxy\uparrow} \right) \ket{0},  \\
\label{eq.20.5}
\ket{-,+,+} & = \frac{1}{3\sqrt{3}} \left( id^{\dagger}_{ixz\uparrow} - d^{\dagger}_{iyz\uparrow} + d^{\dagger}_{ixy\downarrow} \right) \left( id^{\dagger}_{jxz\uparrow} + d^{\dagger}_{jyz\uparrow} + d^{\dagger}_{jxy\downarrow} \right) \left( id^{\dagger}_{kxz\downarrow} + d^{\dagger}_{kyz\downarrow} + d^{\dagger}_{kxy\uparrow} \right) \ket{0}, \\
\label{eq.20.6}
\ket{-,+,-} & = \frac{1}{3\sqrt{3}} \left( id^{\dagger}_{ixz\uparrow} - d^{\dagger}_{iyz\uparrow} + d^{\dagger}_{ixy\downarrow} \right) \left( id^{\dagger}_{jxz\uparrow} + d^{\dagger}_{jyz\uparrow} + d^{\dagger}_{jxy\downarrow} \right) \left( id^{\dagger}_{kxz\downarrow} - d^{\dagger}_{kyz\downarrow} + d^{\dagger}_{kxy\uparrow} \right) \ket{0}, \\
\label{eq.20.7}
\ket{-,-,+} & = \frac{1}{3\sqrt{3}} \left( id^{\dagger}_{ixz\uparrow} - d^{\dagger}_{iyz\uparrow} + d^{\dagger}_{ixy\downarrow} \right) \left( id^{\dagger}_{jxz\uparrow} - d^{\dagger}_{jyz\uparrow} + d^{\dagger}_{jxy\downarrow} \right) \left( id^{\dagger}_{kxz\downarrow} + d^{\dagger}_{kyz\downarrow} + d^{\dagger}_{kxy\uparrow} \right) \ket{0}, \\
\label{eq.20.8}
\ket{-,-,-} & = \frac{1}{3\sqrt{3}} \left( id^{\dagger}_{ixz\uparrow} - d^{\dagger}_{iyz\uparrow} + d^{\dagger}_{ixy\downarrow} \right) \left( id^{\dagger}_{jxz\uparrow} - d^{\dagger}_{jyz\uparrow} + d^{\dagger}_{jxy\downarrow} \right) \left( id^{\dagger}_{kxz\downarrow} - d^{\dagger}_{kyz\downarrow} + d^{\dagger}_{kxy\uparrow} \right) \ket{0}.
\end{align}
\end{subequations}

\subsection{Third-order effective form: Induced loop current operator \label{sec.7}}

In this section, we outline the derivation the induced localized loop current operator in the third-order perturbation expansion in terms of the TB parameters $t_2,t_2'$ using Eq.~(\ref{eq.13}). Explicitly writing Eq.~(\ref{eq.13}) with the individual hoppings, we obtain the operator expression for the localized loop current as
\begin{align}\label{eq.21}
\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}
=
\frac{\mathcal{I}_0t_2^2t_2'}{U_{\Gamma}U_{\Gamma'}}
\sum_{\{\alpha, \sigma \}}
\bigg[
&		 \mathcal{P} i d^{\dagger}_{i \alpha \sigma''}d_{k \beta \sigma''} \mathcal{Q}_{\Gamma'}\mathcal{Q}_{\Gamma'} d^{\dagger}_{k \gamma \sigma'}d_{j \delta \sigma'} \mathcal{Q}_{\Gamma}\mathcal{Q}_{\Gamma} d^{\dagger}_{j \eta \sigma}d_{i \kappa \sigma}  \mathcal{P} \nonumber \\
+ &      \mathcal{P} i d^{\dagger}_{k \alpha \sigma''}d_{j \beta \sigma''} \mathcal{Q}_{\Gamma'}\mathcal{Q}_{\Gamma'} d^{\dagger}_{i \gamma \sigma'}d_{k \delta \sigma'} \mathcal{Q}_{\Gamma} \mathcal{Q}_{\Gamma} d^{\dagger}_{j \eta \sigma}d_{i \kappa \sigma}  \mathcal{P} \nonumber \\
+ &      \mathcal{P} i d^{\dagger}_{j \alpha \sigma''}d_{i \beta \sigma''} \mathcal{Q}_{\Gamma'}\mathcal{Q}_{\Gamma'} d^{\dagger}_{i \gamma \sigma'}d_{k \delta \sigma'}  \mathcal{Q}_{\Gamma}\mathcal{Q}_{\Gamma} d^{\dagger}_{k \kappa \sigma}d_{j \eta \sigma}  \mathcal{P} \nonumber  \\
+ & 	 \mathcal{P} i d^{\dagger}_{j \alpha \sigma''}d_{i \beta \sigma''} \mathcal{Q}_{\Gamma'}\mathcal{Q}_{\Gamma'} d^{\dagger}_{k \gamma \sigma'}d_{j \delta \sigma'} \mathcal{Q}_{\Gamma}\mathcal{Q}_{\Gamma} d^{\dagger}_{i \kappa \sigma}d_{k \eta  \sigma}  \mathcal{P}  \nonumber \\
+ &      \mathcal{P} i d^{\dagger}_{k \alpha \sigma''}d_{j \beta \sigma''} \mathcal{Q}_{\Gamma}\mathcal{Q}_{\Gamma}  d^{\dagger}_{j \gamma \sigma}d_{i \delta \sigma} \mathcal{Q}_{\Gamma'}\mathcal{Q}_{\Gamma'} d^{\dagger}_{i \eta \sigma'}d_{k \kappa \sigma'} \mathcal{P} \nonumber \\
+ & 	 \mathcal{P} i d^{\dagger}_{i \alpha \sigma''}d_{k \beta \sigma''} \mathcal{Q}_{\Gamma}\mathcal{Q}_{\Gamma}  d^{\dagger}_{j \delta \sigma}d_{i \gamma \sigma} \mathcal{Q}_{\Gamma'}\mathcal{Q}_{\Gamma'}d^{\dagger}_{k \eta \sigma'}d_{j \kappa \sigma'} \mathcal{P} +  \mathrm{h.c.} \bigg],
\end{align}
where we have the projection operator $\mathcal{Q}$ to $\mathcal{Q}_{\Gamma}$ to denote all the fifteen eigenstates as defined in Eq.~(\ref{eq.18.1}-\ref{eq.18.15}) with $U_{\Gamma}$ being the corresponding eigen-energy [see Eq.~(\ref{eq.17.1}-\ref{eq.17.4})], and $\mathcal{I}_0 = et_2'\hat{\ensuremath{\mathbf{r}}}_{ij}/\hbar$ is the amplitude of the current operator defined on the $\mathbf{z}'$ bond [\textit{see Fig.~1(b) in the main text}]. Note that we defined the current operator on the longer $\langle ij \rangle$-bond [see Fig.~1(b) in the main text]. After careful analysis of the six terms in Eq.~(\ref{eq.21}), we notice that there are two nonequivalent classes of hopping processes: (a) intermediate states are two distinct doublets at two different sites, and (b) intermediate state is a single doublet at the same site. Considering our geometry [\textit{see Fig.~1(b) in the main text}], we can write these processes (a), (b) as:
\begin{itemize}
\item (a) $i\rightarrow j \Rightarrow j \rightarrow k \Rightarrow k \rightarrow i$, and $j \rightarrow k \Rightarrow k \rightarrow i \Rightarrow i\rightarrow j$ and $ k \rightarrow i \Rightarrow i \rightarrow j \Rightarrow j \rightarrow k$,
\item (b) $i\rightarrow j \Rightarrow \underline{k \rightarrow i} \Rightarrow j \rightarrow k$, and $j \rightarrow k \Rightarrow \underline{i \rightarrow j} \Rightarrow k\rightarrow i$ and $k\rightarrow i \Rightarrow \underline{j \rightarrow k} \Rightarrow i \rightarrow j$.
\end{itemize}
We further notice that in the process (b) we always have one hopping process that connects two singly occupied sites. As the initial and final configuration are constrained by the eight states with $J_{\mathrm{eff}} = 1/2$ total angular momentum [\textit{see Eqs.~(\ref{eq.20.1}-\ref{eq.20.8})}], the overlap amplitude between such singly occupied states with $t_2$/$t_2'$ [\textit{underlined processes in (b)}] hopping always vanishes. Hence, the only contribution to the induced current operator comes from the three processes listed in item [(a)]. We use DiracQ package in Mathematica~\cite{Shastry2013} to compute the matrix elements in Eq.~(\ref{eq.21}). Adding all the three terms illustrated in item [(a)], we obtain the final expression for the induced localized current operator within a triangular plaquette [see Fig.~1(b) in the main text] as
\begin{subequations}
\begin{align}
\nonumber
\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}
=
&
	  \bm{\mathcal{I}}_0	(3S^x_i S^y_j +
	 				     S^y_i S^x_j) S^z_k +	\\
\nonumber
& \quad
	 \bm{\mathcal{I}}_0	(3S^y_i S^z_j-
					   	 5S^z_i S^y_j) S^x_k	+	\\
\label{eq.22.1}
& \quad \quad
	 \bm{\mathcal{I}}_0   (3S^z_i S^x_j-
					     5S^x_i S^z_j) S^y_k, \\
\label{eq.22.2}
\bm{\mathcal{I}}_0 & =   	\frac{e\hat{\ensuremath{\mathbf{r}}}_{ij}}{\hbar} \frac{8t_2^2t_2'J_{\mathrm{H}}(U - 2J_{\mathrm{H}})}{9(U^2 - 4 UJ_{\mathrm{H}} + 3 J_{\mathrm{H}}^2)^2}.
\end{align}
\end{subequations}
First, we notice that SU(2) symmetry of the spins is absent with the additional three coefficients $1,3, 5$ and differ from the well-known form for the current operator in the SU(2) symmetric single-orbital Hubbard model [\textit{see Eq.~(1) in the main text}]~\cite{PhysRevB.78.024402}. Furthermore, the most important property of this structure is that it does not contain any repeated terms such as $S^{\alpha}_i S^{\beta}_j S^{\gamma}_k$, where $\{\alpha, \beta, \gamma\}$ = $\{ x,y,z \}$ can be equal to each other. It is an artifact of retaining only $t_2$ and $t_2'$ hopping terms in the Hamiltonian, which preserves the integrability of the Kitaev Hamiltonian. Inclusion of other hopping parameters like $t_1$, $t_3$ for more realistic modeling would modify this structure with additional three-spin terms with repeated spin indices~\cite{PhysRevLett.112.077204,Kumar2022}. However, here we skip the analysis of the effect of the non-integrable terms in the Kitaev model and leave it for another future work.

Plugging in the parameter values from Sec.~(\ref{sec.7}), we estimate the overall magnitude of the loop current as $|\bm{\mathcal{I}}_0| \sim 30$ nA. Since this current flows around the sides of a triangle in a honeycomb plaquette, an induced magnetic field would appear in the center of the triangle, as shown in Fig.~\ref{fig:Fig1}. Considering the lattice constant $a_0 = 3.44$ {\AA} for $\alpha$-RuCl$_3$~\cite{PhysRevB.93.155143} and assuming a finite expectation value of the induced circulating current operator in the ground state, we obtain an induced out-of-plane magnetic field at the center of the triangle as
\begin{equation}\label{eq.23}
\mathbf{B}_{\perp} = \frac{\mu_0}{4\pi} \int_{\mathrm{C}} \frac{|\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}| d\mathbf{l}  \times \hat{\ensuremath{\mathbf{r}}}}{\ensuremath{\mathbf{r}}^2},
\end{equation}
where $\mathrm{C}$ is the contour of the triangular plaquette. We consider $\mu_0 = 4\pi \times 10^{-7}$ NA$^{-2}$ as the free space magnetic permeability, and $\ensuremath{\mathbf{r}}$ is the distance from the sides of the triangle. After a straight-forward trigonometric and algebraic analysis, we evaluate $\mathbf{B}_{\perp}$ as
\begin{equation}\label{eq.24}
\mathbf{B}_{\perp} = \frac{\mu}{2\pi} |\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}| \left( \frac{\sin \tfrac{\phi_1}{2}}{d_1} + \frac{2 \sin 2\phi_2 + 2 \sin \tfrac{\phi_2}{2}}{d_2} \right) \hat{\bm{c}},
\end{equation}
where $d_1 = \tfrac{a_0}{6} \tan \tfrac{\phi_2}{2}$, $d_2 = \tfrac{a_0}{2\sqrt{3}}$, $\phi_1 = 150^\circ$, $\phi_2 = 30^\circ$ in Fig.~\ref{fig:Fig1}, and $\hat{\bm{c}}$ is the unit vector along the crystalline $c$-axis. Note that an equilateral triangular loop (with sides $2L$), carrying current of amplitude $\mathcal{I}$, would produce a similar magnetic field at its center as $\mathbf{B}_{\perp} = 9\mu_0 \mathcal{I}/(4\pi L)\hat{\bm{c}}$.

\begin{figure}[t]
\centering
\includegraphics[width=0.4\linewidth]{Sub_Fig1}
\caption{The triangular plaquette carrying the induced circulating current around its three sides (solid line). The induced magnetic field is out-of-plane (along the crystalline $c$-axis) as illustrated in the center of the triangle formed by the sites $i,j,k$. The respective angles $\phi_1$, and $\phi_2$ and the shortest distances from the center of the triangle to its sides $d_1$, and $d_2$ are illustrated, respectively.}\label{fig:Fig1}
\end{figure}

\section{Low-energy effective model \label{sec.8}}

\begin{figure}[t]
\centering
\includegraphics[width=0.5\linewidth]{Sub_Fig2}
\caption{The gapless (a) and gapped (b) band structure of the real Majorana fermions obtained by diagonalizing the Hamiltonian in Eq.~({\ref{eq.29}}). Note that we illustrate the full Brillouin zone (BZ) although only half of the BZ contributes to the physical Hilbert space. The Dirac crossing point is illustrated by the band touching for the isotropic and homogeneous Kitaev model in panel (a). The parameter values for $K$ and $\kappa = 0.5 K$ are inflated here in order to show the gap opening near the Dirac point in (b).}\label{fig:Fig2}
\end{figure}

In this section, we briefly outline the derivation of the second-order effective Hamiltonian [\textit{see Eq.~(2a) in the main text}] and discuss its exact solution in Majorana representation~\cite{Kitaev2006}. From Eq.~(\ref{eq.6.2}) we obtain
\begin{align}\label{eq.25}
\mathcal{H}^{(2)}_{\mathrm{eff}} & = \frac{i}{2} \mathcal{P} \ensuremath{\mathcal{S}}^{(1)} \mathcal{Q}_{\Gamma} \mathcal{Q}_{\Gamma} \mathcal{H}_1 \mathcal{P} + \mathrm{h.c.}
=  - \frac{1}{2} \sum_{\Gamma} \frac{1}{U_{\Gamma}} \mathcal{P} \mathcal{H}_1 \mathcal{Q}_{\Gamma} \mathcal{Q}_{\Gamma} \mathcal{H}_1 \mathcal{P} + \mathrm{h.c.},
\end{align}
where we utilized Eq.~(\ref{eq.8.1}) to recast the generating function $\ensuremath{\mathcal{S}}^{(1)}$ in terms of the hopping elements. Again using DiracQ package~\cite{Shastry2013}, we evaluate the matrix elements between the two sites $\langle ij \rangle$, and consequently obtain the effective Hamiltonian as
\begin{equation}\label{eq.26}
\mathcal{H}_0
=
- K \sum_{\langle ij \rangle_{\mathbf{x}}} S^x_i S^x_j
- K \sum_{\langle ij \rangle_{\mathbf{y}}} S^y_i S^y_j
- K \sum_{\langle ij \rangle_{\mathbf{z}}} S^z_i S^z_j.
\end{equation}
In the presence of an external magnetic field $\bm{\mathrm{h}} = (\mathrm{h}_x, \mathrm{h}_y, \mathrm{h}_z)$, a Zeeman term is introduced to Eq.~(\ref{eq.26}) as $\mathcal{H}_{\mathrm{mag}} = \sum_{i} \mathbf{h} \cdot \mathbf{S}_i$. However, in this case the total system $\mathcal{H}_0 + \mathcal{H}_{\mathrm{mag}}$ becomes non-integrable. For a small magnetic field, we consider a low-energy effective form of $\mathcal{H}_{\mathrm{mag}}$~\cite{Kitaev2006}, in perturbative expansion, as
\begin{equation}\label{eq.27}
\mathcal{H}_{\mathrm{eff}} = \kappa \sum_{\substack{\langle ijk \rangle \\ \bigtriangleup}} S^x_i S^y_j S^z_k,
\end{equation}
where $\kappa = \mathrm{h}_x \mathrm{h}_y \mathrm{h}_z/K^2$, and $\Delta$ corresponds to sites within the enclosed triangle with $\langle ijk \rangle$. Utilizing Majorana representation $S^{\alpha}_i = ib^{\alpha}_i c_i/2$, and link variables $u_{ij}$ as $u_{ij} = i b^{\alpha}_i b^{\alpha}_j$, Eq.~(\ref{eq.26}) can be simplified as
\begin{equation}\label{eq.28}
\mathcal{H} = \mathcal{H}_0 + \mathcal{H}_{\mathrm{eff}}
=
\frac{iK}{4} \sum_{\langle ij \rangle} u_{ij} c_i c_j + \frac{i\kappa}{8} \sum_{\llangle ij \rrangle} u_{ik}u_{kj} c_i c_j
=
\frac{iK}{4} \sum_{\langle ij \rangle} c_i c_j + \frac{i\kappa}{8} \sum_{\llangle ij \rrangle} c_i c_j,
\end{equation}
where the link variables $u_{ij}$ have been set with their eigenvalues $+1$, as all $u_{ij}$'s commute with the Hamiltonian and hence can be considered conserved quantities. This particular choice of the link variables is equivalent to attaching zero-flux in each honeycomb plaquette~\cite{Kitaev2006}. Following Lieb's flux theorem~\cite{PhysRevLett.73.2158}, we identify this to be the ground state configuration for the underlying Majorana fermions~\cite{Kitaev2006}. Fourier transforming into the momentum-space, we obtain
\begin{equation}\label{eq.29}
\mathcal{H} = \sum_{\ensuremath{\mathbf{q}}} \mathcal{M}^{{\mathrm{A},\mathrm{B}}}_{\ensuremath{\mathbf{q}}} c_{\ensuremath{\mathbf{q}},\mathrm{A}} c_{-\ensuremath{\mathbf{q}},\mathrm{B}} + \mathrm{h.c.}, \quad
\mathcal{M}_{\ensuremath{\mathbf{q}}}
=
\begin{pmatrix}
\Delta_\ensuremath{\mathbf{q}}		&	if_\ensuremath{\mathbf{q}} \\
-if^*_{\ensuremath{\mathbf{q}}}		&	-\Delta_\ensuremath{\mathbf{q}}
\end{pmatrix},
\end{equation}
where $\mathrm{A}$, and $\mathrm{B}$ are the sub-lattice degrees of freedom. Two functions defined in the matrix $\mathcal{M}_\ensuremath{\mathbf{q}}$ are
\begin{subequations}
\begin{align}
\label{eq.30.1}
f_{\ensuremath{\mathbf{q}}} 		& = \frac{K}{4} \left(e^{i \ensuremath{\mathbf{q}} \cdot \mathbf{a}_1} + e^{i \ensuremath{\mathbf{q}} \cdot \mathbf{a}_2} + 1 \right)\\
\label{eq.30.2}
\Delta_{\ensuremath{\mathbf{q}}} 	& = \frac{\kappa}{4} \big[ \sin \ensuremath{\mathbf{q}} \cdot \mathbf{a}_1 - \sin \ensuremath{\mathbf{q}} \cdot \mathbf{a}_2 + \sin \ensuremath{\mathbf{q}} \cdot (\mathbf{a}_2 - \mathbf{a}_1) \big],
\end{align}
\end{subequations}
where $\mathbf{a}_i, i=1,2$ are the two basis vectors in the honeycomb lattice, given by $\mathbf{a}_1 = (1/2, \sqrt{3}/2)$ and $\mathbf{a}_2 = (-1/2,\sqrt{3}/2)$. Note that only half of the Brillouin zone (BZ) contributes to the physical Hilbert space as Majorana fermions are real quasiparticles with the property $c^{\dagger}_{\ensuremath{\mathbf{q}},\alpha} = c_{-\ensuremath{\mathbf{q}},\alpha}, \alpha = (\mathrm{A},\mathrm{B})$. The gapless (gapped) Dirac spectrum of the Majorana fermions is illustrated in Fig.~\ref{fig:Fig2}, in the absence (presence) of an external magnetic field. In both the cases, the ground state in the extended Hilbert space is written as $\ket{\Psi_0} = \ket{\mathcal{M}} \otimes \ket{\mathcal{G}_0}$, where $\ket{\mathcal{M}}$ is the Majorana fermion ground state in a uniform gauge configuration $\ket{\mathcal{G}_0}$.

So far, we ignored the effect of the external magnetic field on the Hamiltonian Eq.~(\ref{eq.14.3}). However, in the presence of an external magnetic field, there will be the orbital coupling of the magnetic field term in the TB Hamiltonian through Peierl's substitution. Such a modification was not considered in the derivation of the induced loop current or the effective Kitaev Hamiltonian, as derived in Eq.~(\ref{eq.22.1}), and Eq.~(\ref{eq.26}), respectively. If we consider such an orbital coupling of the external magnetic field, after a straightforward analysis as outlined in Sec.~(\ref{sec.7}), we obtain a third-order contribution to the Hamiltonian in Eq.~(\ref{eq.26}). The corresponding term can be written as
\begin{equation}\label{aeq.1}
\mathcal{H}^{(3)}_{\mathrm{eff}} = \sum_{\alpha,\beta,\gamma} \sin \left( \frac{\phi}{\phi_0} \right) \mathcal{A}_{\alpha,\beta,\gamma} S^{\alpha}_i S^{\beta}_j S^{\gamma}_k,
\end{equation}
where the coefficients $\mathcal{A}_{\alpha,\beta,\gamma}$ can be obtained in a similar fashion as in Sec.~(\ref{sec.7}), $\phi$ is the total flux within the triangular plaquette due to the external magnetic field, and $\phi_0 = \hbar c/e$ is the flux quantum~\cite{PhysRevB.51.1922}. For a small magnetic field ($\sim 10$ T)~\cite{Czajka2021,Tanaka2022,Bruin2022}, this term is extremely small and we ignore it for the subsequent discussions.

Finally, we comment about the relation between the components of the external magnetic field $\mathbf{h} = (\mathrm{h}_x, \mathrm{h}_y, \mathrm{h}_z)$ in the octahedral geometry and $\mathbf{h}^{\mathrm{crys}} = \mathrm{h} (\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)$ in the crystalline geometry. Doing a straightforward coordinate transformation between the crystalline and octahedral geometry, we can write $\mathbf{h}$~\cite{Yokoi2021}
\begin{equation}\label{aeq.2}
(\mathrm{h}_x, \mathrm{h}_y, \mathrm{h}_z)
=
\mathrm{h} \left(
\frac{\cos \theta}{\sqrt{3}} + \frac{\sin\theta \cos \phi}{\sqrt{6}} - \frac{\sin \theta \sin \phi}{\sqrt{2}}, \;
\frac{\cos \theta}{\sqrt{3}} + \frac{\sin\theta \cos \phi}{\sqrt{6}} + \frac{\sin \theta \sin \phi}{\sqrt{2}}, \;
\frac{\cos \theta}{\sqrt{3}} - \sqrt{\frac{2}{3}} \sin\theta \cos \phi \right),
\end{equation}
where $\phi$ and $\theta$ are azimuthal angle measured from the crystalline $a$ axis and polar angle from the crystalline $c$ axis, respectively, and $\mathrm{h}$ is the strength of the applied magnetic field. Plugging it back in the Majorana gap $\kappa$ and considering an in-plane magnetic field ($\theta = 90^{\circ}$), we obtain the variation of the strength of $\kappa$ upon a rotation of the in-plane magnetic field as illustrated in Fig.~\ref{fig:Fig3}(a,b).

\begin{figure}[t]
\centering
\includegraphics[width=1.0\linewidth]{Sub_Fig3}
\caption{(a) Contour plot of $\mathrm{sign}(\kappa)$  $ (\kappa = \mathrm{h}_x \mathrm{h}_y \mathrm{h}_z/K^2)$ from Eq.~(\ref{aeq.2}) in the $\theta$-$\phi$ plane for the different orientations of the external magnetic field. (b) The variation of the magnitude of the Majorana gap function $\kappa$ as a function of the azimuthal angle measured from crystalline axis $a$ for an in-plane external magnetic field with $\theta = 90^{\circ}$.}\label{fig:Fig3}
\end{figure}

\section{Computation: loop current expectation value \label{sec.9}}

Now, we focus on the analysis for the evaluation of the expectation value of the loop current operator in the Kitaev ground state with Majorana representation. We first consider the uniform gauge configuration with all $u_{ij} = 1$. In this case, the analysis is simple with the Majorana fermion in the momentum-space representation. The target quantity is evaluating three-spin correlation functions \textit{i.e.} $\braket{S_i^{\alpha} S_j^{\beta} S_k^{\gamma}}$, where $ \{\alpha,\beta,\gamma \} = \{ x,y,z \}$ with $\alpha \neq \beta \neq \gamma$ [\textit{see Eq.~(\ref{eq.22.2})}]. Along with the link variables, we further define the gauge invariant loop operators $\mathcal{W}_p = \prod_{\braket{ij} \in \mathord{\raisebox{-0.5pt}{\tikz{\node[draw,scale=.5,regular polygon, regular polygon sides=6,](){};}}}_p} u_{ij}$ in each hexagonal plaquette. In the uniform gauge, as mentioned earlier, $\mathcal{W}_p = 1$ for all the honeycomb cells. However, flipping one link variable in a single hexagon leads to $\mathcal{W}_p = -1$ in the two neighboring hexagons. The energy of the Majorana fermions increases with such non-uniform gauge configurations, and leads to the emergence of $\mathbb{Z}_2$ vortices, aka. visons~\cite{Balents2010,Savary_2016}. Each vison excitation creates a $\pi$-flux in the associated hexagon, and is static due to the integrable structure of the Kitaev model.

Hence, the non-zero expectation value of an arbitrary operator $\mathcal{O}_i$ exists only if the operator $\mathcal{O}_i$ conserves the vison occupation number. In this context, the action of three spin-operator $S_i^{\alpha} S_j^{\beta} S_k^{\gamma}$ on a ground state with finite/zero number vison excitations, has to preserve the vison occupation number. We note that a single spin operator $S_i^{\alpha}$ leads to two vison excitations and can be symbolically written as~\cite{PhysRevLett.98.247201}
\begin{equation}\label{eq.31}
S^{\alpha}_i \rightarrow i c_i \hat{\pi}_1{}_{\langle ij \rangle_\alpha} \hat{\pi}_2{}_{\langle ij \rangle_\alpha},
\end{equation}
where $\hat{\pi}_1{}_{\langle ij \rangle_\alpha}$ and $\hat{\pi}_2{}_{\langle ij \rangle_\alpha}$ are operators that introduce $\pi$-fluxes to the plaquettes $1$ and $2$ shared by the bond $\langle ij \rangle_{\alpha}$. Therefore, if the ground state of the Kitaev system does not have any vison excitations, a necessary condition for a finite expectation value of $S_i^{\alpha} S_j^{\beta} S_k^{\gamma}$ in this ground state, is that the two visons created by $S_k^{\gamma}$ should be destroyed by the other two spin operators.

\subsection{Pure Kitaev model in uniform gauge \label{sec.10}}

In this case, the ground state $\ket{\Psi_0}$ was obtained from Eq.~(\ref{eq.29}) where $\ket{\Psi_0} = \ket{\mathcal{M}} \otimes \ket{\mathcal{G}_0}$. Here, $\ket{\mathcal{G}_0}$ corresponds to the uniform gauge configuration without any vison excitations. After a straightforward algebra with Eq.(\ref{eq.31}), we see that only the first term in Eq.(\ref{eq.10.1}) satisfy the constraint as mentioned in Sec.~\ref{sec.9}. Hence, we have [using the Majorana representation of the spin operators]
\begin{equation}\label{eq.32}
\braket{\Psi_0|\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}|\Psi_0} = i \frac{\bm{\mathcal{I}}_0}{8}\braket{\Psi_0|c_i c_j|\Psi_0} = i \frac{\bm{\mathcal{I}}_0}{8}\braket{\mathcal{M}|c_i c_j|\mathcal{M}},
\end{equation}
where $\ket{\mathcal{M}}$ is the ground state of Majorana fermions. Choosing a convention that $i,j$ lie in sub-lattice B, we can rewrite $c_i c_j$ in momentum space as
\begin{align}
\nonumber
&	c_{i,\mathrm{B}} c_{j,\mathrm{B}} = \sum_{\ensuremath{\mathbf{k}},\ensuremath{\mathbf{k}}'} e^{i\ensuremath{\mathbf{k}}.\ensuremath{\mathbf{r}}_i} e^{i\ensuremath{\mathbf{k}}'.\ensuremath{\mathbf{r}}_j} c_{\ensuremath{\mathbf{k}},\mathrm{B}} c_{\ensuremath{\mathbf{k}}',\mathrm{B}} \\
\label{eq.33}
								  & = \sum_{\substack{\ensuremath{\mathbf{k}},\ensuremath{\mathbf{k}}' \\ \in \\ \mathrm{HBZ}}} e^{i\ensuremath{\mathbf{k}}.\ensuremath{\mathbf{r}}_i} e^{i\ensuremath{\mathbf{k}}'.\ensuremath{\mathbf{r}}_j} c_{\ensuremath{\mathbf{k}},\mathrm{B}} c_{\ensuremath{\mathbf{k}}',\mathrm{B}} +
									  \sum_{\substack{\ensuremath{\mathbf{k}},\ensuremath{\mathbf{k}}' \\ \in \\ \mathrm{HBZ}}} e^{-i\ensuremath{\mathbf{k}}.\ensuremath{\mathbf{r}}_i} e^{-i\ensuremath{\mathbf{k}}'.\ensuremath{\mathbf{r}}_j} c^{\dagger}_{\ensuremath{\mathbf{k}},\mathrm{B}} c^{\dagger}_{\ensuremath{\mathbf{k}}',\mathrm{B}} +
									  \sum_{\substack{\ensuremath{\mathbf{k}},\ensuremath{\mathbf{k}}' \\ \in \\ \mathrm{HBZ}}} e^{-i\ensuremath{\mathbf{k}}.\ensuremath{\mathbf{r}}_i} e^{i\ensuremath{\mathbf{k}}'.\ensuremath{\mathbf{r}}_j} c^{\dagger}_{\ensuremath{\mathbf{k}},\mathrm{B}} c_{\ensuremath{\mathbf{k}}',\mathrm{B}} +
  									  \sum_{\substack{\ensuremath{\mathbf{k}},\ensuremath{\mathbf{k}}' \\ \in \\ \mathrm{HBZ}}} e^{i\ensuremath{\mathbf{k}}.\ensuremath{\mathbf{r}}_i} e^{-i\ensuremath{\mathbf{k}}'.\ensuremath{\mathbf{r}}_j} c_{\ensuremath{\mathbf{k}},\mathrm{B}} c^{\dagger}_{\ensuremath{\mathbf{k}}',\mathrm{B}},
\end{align}
where we defined the momentum summation in the first line over the full Brillouin zone (BZ) and reduced it to the half-Brillouin zone (HBZ) in the second line. Next, we represent the sub-lattice operators on the diagonal basis within the HBZ as
\begin{equation}\label{eq.34}
\begin{pmatrix}
c^{\dagger}_{\ensuremath{\mathbf{k}},\mathrm{A}} \\
c^{\dagger}_{\ensuremath{\mathbf{k}},\mathrm{B}}
\end{pmatrix}
=
\frac{1}{\sqrt{2}}
\begin{pmatrix}
-m_{\ensuremath{\mathbf{k}}} &  m_{\ensuremath{\mathbf{k}}} \\
1 		&   1
\end{pmatrix}
\begin{pmatrix}
\alpha^{\dagger}_{\ensuremath{\mathbf{k}}} \\
\beta^{\dagger}_{\ensuremath{\mathbf{k}}}
\end{pmatrix},
\end{equation}
where $m_{\ensuremath{\mathbf{k}}} = i f_{\ensuremath{\mathbf{k}}}/|f_{\ensuremath{\mathbf{k}}}|$ from Eq.~(\ref{eq.29}) (in the \textit{absence} of any external magnetic field). Finally, the ground state is obtained by filling all the negative energy states and is written as
\begin{equation}\label{eq.35}
\ket{\mathcal{M}} = \prod_{\substack{\ensuremath{\mathbf{k}} \\ \in \\ \mathrm{HBZ}}} \beta^{\dagger}_{\ensuremath{\mathbf{k}}} \ket{0}.
\end{equation}
Consequently, a simple numerical integration (in Mathematica) using Eqs.~(\ref{eq.33}-\ref{eq.35}) yields
\begin{equation}\label{eq.36}
\braket{\Psi_0|\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}|\Psi_0} = 0.
\end{equation}
Hence, in the absence of any external magnetic field, the ground state of the pure Kitaev model does not allow any localized current expectation value. This result is consistent with the global time-reversal symmetry of the underlying system. Therefore, to induce a non-zero expectation value for the localized current operator, we break the time-reversal symmetry by applying an external magnetic field.

\begin{figure}[t]
\centering
\includegraphics[width=0.5\linewidth]{Sub_Fig4}
\caption{The variation of average localized current (in unit of $|\bm{\mathcal{I}}_0|$) as a function of $\kappa$. For small $\kappa$ we notice the linear dependence of the current operator. \textit{Note}: We utilized inflated values for $\kappa$ to illustrate the dependence.}\label{fig:Fig4}
\end{figure}

\subsection{Kitaev model in external magnetic field \label{sec.11}}

In the presence of a magnetic field, the Majorana fermions acquire a gap near the Dirac point [\textit{see Fig.~\ref{fig:Fig2}(b)}]. In this case, the ground state is written as $\ket{\Psi_{\mathrm{mag}}} = \ket{\mathcal{M}_{\mathrm{mag}}} \otimes \ket{\mathcal{G}_0}$, where $\ket{\mathcal{M}_{\mathrm{mag}}}$ is ground state of the Hamiltonian in Eq.~(\ref{eq.29}). Writing the sub-lattice operators in terms of the diagonal operators we have
\begin{equation}\label{eq.37}
\begin{pmatrix}
c^{\dagger}_{\ensuremath{\mathbf{k}},\mathrm{A}} \\
c^{\dagger}_{\ensuremath{\mathbf{k}},\mathrm{B}}
\end{pmatrix}
=
\begin{pmatrix}
-i\tfrac{f_\ensuremath{\mathbf{k}} (\varepsilon_\ensuremath{\mathbf{k}} + \Delta_{\ensuremath{\mathbf{k}}})}{|f_\ensuremath{\mathbf{k}}|\sqrt{|f_\ensuremath{\mathbf{k}}|^2 + (\varepsilon_\ensuremath{\mathbf{k}} + \Delta_{\ensuremath{\mathbf{k}}})^2}} &  i\tfrac{f_\ensuremath{\mathbf{k}} (\varepsilon_\ensuremath{\mathbf{k}} - \Delta_{\ensuremath{\mathbf{k}}})}{|f_\ensuremath{\mathbf{k}}|\sqrt{|f_\ensuremath{\mathbf{k}}|^2 + (\varepsilon_\ensuremath{\mathbf{k}} - \Delta_{\ensuremath{\mathbf{k}}})^2}} \\
\tfrac{|f_\ensuremath{\mathbf{k}}|}{\sqrt{|f_\ensuremath{\mathbf{k}}|^2 + (\varepsilon_\ensuremath{\mathbf{k}} + \Delta_{\ensuremath{\mathbf{k}}})^2}} 									  	   &  \tfrac{|f_\ensuremath{\mathbf{k}}|}{\sqrt{|f_\ensuremath{\mathbf{k}}|^2 + (\varepsilon_\ensuremath{\mathbf{k}} - \Delta_{\ensuremath{\mathbf{k}}})^2}}
\end{pmatrix}
\begin{pmatrix}
\alpha^{\dagger}_{\ensuremath{\mathbf{k}}} \\
\beta^{\dagger}_{\ensuremath{\mathbf{k}}}
\end{pmatrix},
\end{equation}
where $\varepsilon_\ensuremath{\mathbf{k}} = \sqrt{|f_\ensuremath{\mathbf{k}}|^2 + \Delta_\ensuremath{\mathbf{k}}^2}$. In a similar fashion to Eq.~(\ref{eq.36}), the Majorana fermion ground state is written in terms of the diagonal operators as
\begin{equation}\label{eq.38}
\ket{\mathcal{M}_{\mathrm{mag}}} = \prod_{\substack{\ensuremath{\mathbf{k}} \\ \in \\ \mathrm{HBZ}}} \beta^{\dagger}_{\ensuremath{\mathbf{k}}} \ket{0}.
\end{equation}
Since the gauge configuration characterized by $\ket{\mathcal{G}_0}$ constraints all the gauge invariant Wilson loops $\mathcal{W}_p = 1$, only one of the six spin combination, in Eq.~(\ref{eq.22.1}), gives non-zero expectation value in the ground state $\ket{\Psi_{\mathrm{mag}}}$. Consequently, we have
\begin{equation}\label{eq.39}
\braket{\Psi_{\mathrm{mag}}|\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}|\Psi_{\mathrm{mag}}} = i \frac{\bm{\mathcal{I}}_0}{8}\braket{\mathcal{M}_{\mathrm{mag}}|c_i c_j|\mathcal{M}_{\mathrm{mag}}}.
\end{equation}
Expanding the Majorana bilinear operator $c_i c_j$ with Eq.~(\ref{eq.37}) in the diagonal basis and after a straight-forward algebra, we obtain
\begin{equation}\label{eq.40}
\braket{\Psi_{\mathrm{mag}}|\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}|\Psi_{\mathrm{mag}}}
=
i \frac{\bm{\mathcal{I}}_0}{8} \mathrm{Re} \, \left(
\sum_{\substack{\ensuremath{\mathbf{k}} \\ \in \\ \mathrm{HBZ}}}
e^{i\ensuremath{\mathbf{k}}.(\ensuremath{\mathbf{r}}_i - \ensuremath{\mathbf{r}}_j)} [a_\ensuremath{\mathbf{k}}^2 + b_\ensuremath{\mathbf{k}}^2]
-2 i b_\ensuremath{\mathbf{k}}^2 \sin [\ensuremath{\mathbf{k}}.(\ensuremath{\mathbf{r}}_i - \ensuremath{\mathbf{r}}_j)]
\right),
\end{equation}
where $a_\ensuremath{\mathbf{k}} = \tfrac{|f_\ensuremath{\mathbf{k}}|}{\sqrt{|f_\ensuremath{\mathbf{k}}|^2 + (\varepsilon_\ensuremath{\mathbf{k}} + \Delta_{\ensuremath{\mathbf{k}}})^2}}$, and $b_\ensuremath{\mathbf{k}} = \tfrac{|f_\ensuremath{\mathbf{k}}|}{\sqrt{|f_\ensuremath{\mathbf{k}}|^2 + (\varepsilon_\ensuremath{\mathbf{k}} - \Delta_{\ensuremath{\mathbf{k}}})^2}}$. We perform the momentum integration in HBZ numerically in Mathematica. The corresponding variation of $\braket{\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}}_{\mathrm{mag}}$  as a function of the Majorana gap parameter $\kappa$ is shown in Fig.~\ref{fig:Fig3} (in unit of $|\bm{\mathcal{I}}_0|$). For small $\kappa$ the linear dependence of the current operator becomes apparent.

\section{Kitaev model in a finite system: \\ Vison excitations and external magnetic field \label{sec.12}}

In the previous two sections, we analyzed the expectation value of the current operator in the absence of any vison excitations in an infinite (periodic boundary condition) system. As mentioned earlier, this corresponds to choosing all the link variables $u_{ij} = 1$, with the constraint that $i$ in each bond $\braket{ij}$ belongs to the sub-lattice A. Flipping a particular bond $u_{ij}$ to $-1$ would, therefore, create two $\pi$-fluxes in the adjacent honeycomb plaquettes as explained in Sec.~(\ref{sec.9}). We can separate these two adjacent $\mathbb{Z}_2$ vortices or visons by a string operator as illustrated in Fig.~\ref{fig:Fig5}(a), to minimize their mutual interactions. Note that all the link variables defined on the bonds that cross the string operator are flipped. However, such vison configurations destroy the translational invariance of the Majorana system, and we cannot simply go to the momentum space to do our analysis.

Consequently, we first consider a spin-system on a 2D honeycomb lattice with $L \times L$ unit cells as shown in Fig.~\ref{fig:Fig5}(a). Since each unit cell contain two sub-lattices (A \& B), there are $2L^2$ Majorana operators $c_i$ in the system. The lattice vectors are chosen as $\mathbf{a}_1 = (1/2,\sqrt{3}/2)$, and $\mathbf{a}_2 = (-1/2,\sqrt{3}/2)$. The Majorana operators at a site $i$ is written as $c_i = c_{\eta} (m,n)$, where $\eta$ corresponds to the sub-lattice index, and $\mathbf{R}(m,n) = m \mathbf{a}_1 + n \mathbf{a}_2$, $m,n = 1,2,\ldots L$. Following Ref.~[\onlinecite{Kitaev2006}], we impose periodic boundary condition (PBC) as $c_{\eta}(m+L,n) = c_{\eta}(m,n)$, and $c_{\eta}(m,n+L) = c_{\eta}(m,n)$. The $2L^2$-dimensional Majorana vector is constructed as $\mathbf{\tilde{c}} = (c_{\mathrm{A}}, c_{\mathrm{B}})^{\mathsf{T}}$ with
\begin{equation}\label{eq.41}
c_{\eta} = (c_\eta(1,1), c_\eta(2,1), \ldots c_\eta(L,1), c_\eta(1,2), c_\eta(2,2), \ldots c_\eta(L,2), \ldots, c_\eta(L,L))^{\mathsf{T}}, \quad \eta \in \{ \mathrm{A}, \mathrm{B} \}.
\end{equation}
In terms of the Majorana vector $\mathbf{\tilde{c}}$, we can write the Kitaev model (in presence of an external magnetic field) as $\mathbf{\tilde{c}}^{\mathsf{T}} \mathsf{H} \mathbf{\tilde{c}}$, where $\mathsf{H}$ is written as
\begin{align}
\nonumber
\mathsf{H} = & 					i \frac{K}{4} \sum_{m,n}
								c_{\mathrm{A}}(m,n) \bigg( u_z(m,n) c_{\mathrm{B}}(m,n) + u_x(m,n) c_{\mathrm{B}}(m+1,n) + u_y(m,n) c_{\mathrm{B}}(m,n+1) \bigg) + \\
\nonumber
			 & 					i\frac{\kappa}{8} \sum_{m,n} \bigg( c_{\mathrm{A}}(m,n) \Big[ u_x(m,n)u_y(m+1,n-1) c_{\mathrm{A}}(m+1,n-1) + u_z(m,n)u_x(m-1,n) c_{\mathrm{A}}(m-1,n) \\
\nonumber
			 & \qquad \quad \;	u_y(m,n)u_z(m,n+1) c_{\mathrm{A}}(m,n+1) \Big] + c_{\mathrm{B}}(m,n) \Big[ u_x(m-1,n)u_y(m-1,n) c_{\mathrm{B}}(m-1,n+1) + \\
\label{eq.42}
			 & \qquad \quad \;	u_z(m,n)u_x(m,n) c_{\mathrm{B}}(m+1,n) + u_y(m,n-1)u_z(m,n-1) c_{\mathrm{B}}(m,n-1) \bigg) + \mathrm{H.c.},
\end{align}
where $u_{\alpha}(m,n) = u_{\langle ij \rangle_{\alpha}}$, with $i$ ($j$) being in the sub-lattice $\mathrm{A}$ ($\mathrm{B}$) as shown in Fig.~\ref{fig:Fig4}(a). We now analyze the Majorana physics in a finite system with both the PBC and the open boundary condition (OBC), by diagonalizing the Hamiltonian in Eq.~(\ref{eq.42}) in the real space. All numerical estimates will be written in the unit of Kitaev coupling $K$ in the subsequent sections.

\begin{figure}[t]
\centering
\includegraphics[width=0.8\linewidth]{Sub_Fig5}
\caption{(a) A schematic of the finite size ($L \times L$, $L = 20$) honeycomb lattice. Two static vison excitations are located at the farthest distance $\lfloor \frac{L-1}{2} \rfloor$ in the middle of the system with periodic boundary condition (PBC). The link variables $u_{\braket{ij}}$ of the bonds crossing the blue dashed line (\textit{string} \textit{operator}) are all flipped to $-1$. (b) \textit{No vison excitations:} Variation of the energy per unit cell $\varepsilon_0$ [\textit{circles}: absence of external magnetic-field; \textit{triangles}: non-vanishing external magnetic-field] as a function of the linear system size $L$ with PBC. Convergence is obtained for $L \sim 50$ ($L \sim 20$) in the absence (presence) of the external magnetic field. The energy is measured in the unit of $4K$. An inflated value of $\kappa = K$ is chosen for illustrative purposes.}\label{fig:Fig5}
\end{figure}

\subsection{Kitaev model: Numerical convergence}

At first, we diagonalize the above Hamiltonian in the \textit{absence} of $\kappa$ and any vison excitations and analyze the ground state as a function of the system size $L$. The eigenmodes of the Hamiltonian in Eq.~(\ref{eq.42}) can be obtained by a canonical transformation to a new set of Majorana operators~\cite{Kitaev2006} as
\begin{equation}\label{eq.43}
(b^{'}_1,b^{''}_1, \ldots b^{'}_{N},b^{''}_{N}) = \mathbf{\tilde{c}}^{\mathsf{T}} \tilde{\mathcal{R}},
\end{equation}
where $\tilde{\mathcal{R}}$ is the canonical matrix. In terms of these operators the Hamiltonian in Eq.~(\ref{eq.41}) can be rewritten as $\mathsf{H} = \tfrac{i}{2} \sum_{m}\varepsilon_m b^{'}_m b^{''}_m$, where $\varepsilon_m$ are the positive eigenvalues of $\mathsf{H}$. Introducing fermions operators as $a_m = \tfrac{1}{2} (b^{'}_m + i b^{''}_m)$, we have $\mathsf{H} = \sum_{m}\varepsilon_m (a^{\dagger}_m a_m - 1/2 )$. Hence, the ground state energy is given by $\mathcal{E}_0 = -\tfrac{1}{2} \sum_m \varepsilon_m$. We now analyze the variation of average ground state energy $\varepsilon_0$ [$=2\mathcal{E}_0/N^2$] per unit cell with the system size. The latter is shown in Fig.~\ref{fig:Fig5}(b) [\textit{empty circles}]. We notice that the average energy converges to 1.5746 (as noted by Kitaev himself~\cite{Kitaev2006}) for a linear system size of around $L \sim 50$. Performing a similar analysis in the presence of an external magnetic field, we notice that the convergence is achieved for a much smaller system with linear size $L \sim 20$. The variation of the energy per unit cell, in the external magnetic field, is illustrated by the (\textit{empty triangles}) in Fig.~\ref{fig:Fig5}(b).

Corresponding results for the convergence with system size in a periodic system and energy distribution in a finite system with PBC are shown in Fig.~\ref{fig:Fig4}(b) and Fig.~\ref{fig:Fig4}(c), respectively. Note that an inflated value of the magnetic field term $\kappa$ has been used in the numerical computation. However, the convergence is obtained at a much smaller system size with linear size $L \sim 20$, because of the gapped Majorana spectrum in Fig.~\ref{fig:Fig2}(b).

Next, we introduce two $\pi$-fluxes in the adjacent honeycomb plaquettes shared by a bond $\langle ij \rangle$. In a periodic system, vison excitations can be only created in pairs since the string operator has to end inside the system as illustrated in Fig.~\ref{fig:Fig5}(a). For this purpose, at first, we consider two adjacent plaquettes with $\pi$-fluxes which is achieved by putting the string operator on only one bond $\langle ij \rangle$. The energy of a single vison for two adjacent $\pi$-flux configuration comes out to be around $\sim 0.1311$. We compute the energy of a single vison by numerically computing the eigen-energy of the Hamiltonian with two visons and measuring the difference with the vison-free ground state as $(\mathcal{E}_2 - \mathcal{E}_0)/2$. Note that our quantitative estimate for the energy of a single vison is slightly less than Kitaev's original estimate ($\sim 0.1536$)~\cite{Kitaev2006}, but is consistent with Ref.~\cite{PhysRevB.84.165414}. The analysis has been performed for a system of linear size $L = 50$ with two visons placed at the center of the system.

\subsection{Kitaev model in magnetic field : Current expectation value}

Now we provide the main part of our analysis. We consider the Kitaev model in the presence of an external magnetic field which is parametrized by $\kappa$ in Eq.~(\ref{eq.42}) for a finite system with both PBC and OBC. We consider two scenarios as mentioned in the main text: (a) the distribution of localized current in the system in absence of vison excitations (gauge sector $\ket{\mathcal{G}_0}$) in a finite system with OBC, and (b) the profile of localized current distribution around two vison excitations placed far away from the each other (gauge sector $\ket{\mathcal{G}_2}$) with PBC. For subsequent analysis it will be easier to transform the Majorana operators in Eq.~(\ref{eq.41}) in terms of the matter fermions~\cite{PhysRevLett.98.247201} as $f_i = c_{\mathrm{A},i} - i c_{\mathrm{B},i}$. On the latter basis, the $2L^2$-Majorana vector is transformed as
\begin{equation}\label{eq.44}
\mathbf{\tilde{c}}
=
\begin{pmatrix}
c_{\mathrm{A}} 	& c_{\mathrm{B}}
\end{pmatrix}^{\mathsf{T}}
=
\mathsf{R}
\begin{pmatrix}
f 	& 	f^{\dagger}
\end{pmatrix}^{\mathsf{T}}, \quad
\mathsf{R}
=
\begin{pmatrix}
\mathbb{I}_{L^2}	&	\mathbb{I}_{L^2} \\
i\mathbb{I}_{L^2}	&	-i\mathbb{I}_{L^2}
\end{pmatrix},
\end{equation}
where $\mathbb{I}_{L^2}$ is an identity matrix of dimension $L^2$, and $f$-vector is defined as before
\begin{equation}\label{eq.45}
f = (f(1,1), f(2,1), \ldots f(L,1), f(1,2), f(2,2), \ldots f(L,2), \ldots, f(L,L))^{\mathsf{T}}.
\end{equation}
In terms of the matter fermion operators, the Hamiltonian in Eq.~(\ref{eq.42}) can be rewritten as
\begin{equation}\label{eq.46}
\mathsf{H} =
\begin{pmatrix}
c_{\mathrm{A}}		& 		c_{\mathrm{B}}
\end{pmatrix}
\begin{pmatrix}
\mathcal{H}_{\mathrm{AA}}		& 		\mathcal{H}_{\mathrm{AB}} \\
\mathcal{H}_{\mathrm{BA}}		& 		\mathcal{H}_{\mathrm{BB}}
\end{pmatrix}
\begin{pmatrix}
c_{\mathrm{A}} \\
c_{\mathrm{B}}
\end{pmatrix}
=
\begin{pmatrix}
f^{\dagger}		& 		f
\end{pmatrix}
\mathsf{R}^{\dagger}
\begin{pmatrix}
\mathcal{H}_{\mathrm{AA}}		& 		\mathcal{H}_{\mathrm{AB}} \\
\mathcal{H}_{\mathrm{BA}}		& 		\mathcal{H}_{\mathrm{BB}}
\end{pmatrix}
\mathsf{R}
\begin{pmatrix}
f	\\
f^{\dagger}
\end{pmatrix}.
\end{equation}
The block Hamiltonian $\mathcal{H}_{\alpha\beta}$ ($\alpha,\beta = \mathrm{A},\mathrm{B}$) is read off from Eq.~(\ref{eq.42}). Instead of the canonical diagonalization as earlier, we now diagonalize the Hamiltonian in Eq.~(\ref{eq.46}) by introducing the normal mode operators $a,\; a^{\dagger}$ as
\begin{equation}\label{eq.47}
\begin{pmatrix}
f \\
f^{\dagger}
\end{pmatrix} =
U
\begin{pmatrix}
a	\\
a^{\dagger}
\end{pmatrix},
\end{equation}
where $U$ is a unitary matrix that diagonalizes the Hamiltonian $\mathsf{H}$ in Eq.~(\ref{eq.46}). The resultant Hamiltonian in the $a$-basis reads as
\begin{equation}\label{eq.48}
\mathsf{H} =
\begin{pmatrix}
f^{\dagger}		& 		f
\end{pmatrix}
\mathsf{R}^{\dagger}
\begin{pmatrix}
\mathcal{H}_{\mathrm{AA}}		& 		\mathcal{H}_{\mathrm{AB}} \\
\mathcal{H}_{\mathrm{BA}}		& 		\mathcal{H}_{\mathrm{BB}}
\end{pmatrix}
\mathsf{R}
\begin{pmatrix}
f	\\
f^{\dagger}
\end{pmatrix}
=
\begin{pmatrix}
a^{\dagger}	&  a
\end{pmatrix}
\underbrace{U^{\dagger}
\mathsf{R}^{\dagger}
\begin{pmatrix}
\mathcal{H}_{\mathrm{AA}}		& 		\mathcal{H}_{\mathrm{AB}} \\
\mathcal{H}_{\mathrm{BA}}		& 		\mathcal{H}_{\mathrm{BB}}
\end{pmatrix}
\mathsf{R}
U}_{V_\mathrm{D}}
\begin{pmatrix}
a	\\
a^{\dagger}
\end{pmatrix},
\end{equation}
where $V_{\mathrm{D}}$ is the diagonal matrix with diagonal entries as the eigenenergies of the Hamiltonian $\mathsf{H}$, and we obtain in terms of the $a$ fermions as
\begin{equation}\label{eq.49}
\mathsf{H} = \sum_{l = 1}^{L^2}\varepsilon_l(2a^{\dagger}_l a_l - 1),
\end{equation}
where $\varepsilon_l (\ge 0)$ are non-negative real eigenvalues of $\mathsf{H}$ ordered as $\varepsilon_1 < \varepsilon_2 < \cdots < \varepsilon_{L^2}$. Consequently, the ground state is defined by $\ket{\mathcal{G}_0} \otimes \ket{\mathcal{M}_{\mathrm{mag}}}$ with the property $a_l \ket{\mathcal{M}_{\mathrm{mag}}} = 0$. Here, gauge sector $\ket{\mathcal{G}_0}$ corresponds to an absence of any vison excitation in the system, and $\ket{\mathcal{M}_{\mathrm{mag}}}$ is the Majorana fermion ground state in the presence of an external magnetic field. The ground state energy is given by $E_0 = -\sum_{l=1}^{L^2} \varepsilon_l$. We also define the bond fermion operator $\chi_{\langle ij \rangle_{\alpha}}$ as
\begin{equation}\label{eq.50}
\chi_{\langle ij \rangle_{\alpha}}
=
\frac{1}{2} \left(b_i^{\alpha} - i b^{\alpha}_{j}\right), \quad
\chi^{\dagger}_{\langle ij \rangle_{\alpha}}
=
\frac{1}{2} \left(b_i^{\alpha} + i b^{\alpha}_{j}\right).
\end{equation}
In terms of the bond fermions, the link variables can be recasted as $u_{ij} = 1 - 2 \chi^{\dagger}_{\langle ij \rangle_{\alpha}} \chi_{\langle ij \rangle_{\alpha}}$. Hence, for a uniform gauge configuration, with all $u_{ij} = 1$, we do not have any bond-fermions to start with. Flipping a specific link variable, hence, creates one bond fermion and associated two Majorana zero modes. The energy distribution for a finite system with PBC of linear size $L = 32$ is shown in Fig.~\ref{fig:Fig6}(a) (\ref{fig:Fig6}(b)) in the absence (presence) of vison excitations. For two vison excitations placed at a maximal distance, we have two Majorana zero modes (MZM) as illustrated in Fig.~\ref{fig:Fig6}(b).

\begin{figure}[t]
\centering
\includegraphics[width=0.8\linewidth]{Sub_Fig6}
\caption{(a) The eigenvalue distribution as a function of the Hilbert space dimension for the Kitaev model in an external magnetic field in the uniform gauge configuration without any vison excitations. (b) The eigenvalue distribution in the same external magnetic field in the presence of two well-separated vison excitations with two Majorana zero modes (MZM) exactly at the zero energy.}\label{fig:Fig6}
\end{figure}

In terms of the new "$a$" fermions, we can now evaluate the expectation value of the loop current operator. The current operator can be written in a compact form as $\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k} = \mathbf{A}_{\alpha\beta\gamma} S^{\alpha}_iS^{\beta}_jS^{\gamma}_k$, where the coefficients $\mathbf{A}_{\alpha\beta\gamma}$ are read off from Eq.~(\ref{eq.22.1}). With a particular choice of the gauge sector $\ket{\mathcal{G}_0}$ or $\ket{\mathcal{G}_2}$, we can recast the above current operator in terms of the Majorana fermions. Since the visons in the gauge sector are static, it leads to a huge simplification in computing the average of the current operator in the ground state. The non-zero expectation value of the terms in the ground state is only the terms in the current operator that preserve the number of the vison excitations. It turns out that only the first term in Eq.~(\ref{eq.22.1}) satisfies this condition.

On the other hand, since each bond in the honeycomb lattice is shared between four triangles, four terms will contribute to the total average current flowing in the bond. Consequently, we obtain the current expression for a particular bond $\alpha$ on the honeycomb plaquette as
\begin{equation}\label{eq.51}
\bm{\mathcal{I}}_{\alpha}
=
\sum_{\triangle} \braket{\Psi|\tilde{\bm{\mathcal{I}}}^{(2)}_{ij,k}|\Psi} = 3\bm{\mathcal{I}}_0 \sum_{k,\triangle} u_{ik}u_{kj} \llangle c_i c_j \rrangle - \text{"neighboring plaquette"},
\end{equation}
where the indices $ij,k$ in each triangle are arranged according to the orientation of the associated triangles, and $\llangle c_i c_j\rrangle$ signifies that expectation value of the operator $c_ic_j$ when the sites $i,j$ are connected by the next nearest-neighbor bonds on the associated triangle. The "neighboring plaquette" term corresponds to the triangles belonging to the other honeycomb plaquette shared by the bond $\alpha$.

\begin{figure}[t]
\centering
\includegraphics[width=0.5\linewidth]{Sub_Fig7}
\caption{The difference between the Majorana fermion ground energy for two vison configuration [with well-separated visons as shown in Fig.~\ref{fig:Fig5}(a)] in the presence and absence of an external magnetic field plotted as a function of $\kappa$. The ground state energy decreases with increasing magnetic field strength.}\label{fig:Fig7}
\end{figure}

Hence, the analysis of the expectation of the current operator eventually boils down to the computation of the expectation value of the product of Majorana operators $\braket{c_{\mathrm{A}}(\mathbf{R}_{mn}) c_{\mathrm{A}}(\mathbf{R}'_{m'n'})}$ and $\braket{c_{\mathrm{B}}(\mathbf{R}_{mn}) c_{\mathrm{B}}(\mathbf{R}'_{m'n'})}$, where $\ensuremath{\mathbf{r}}_{mn}$ denotes the position of the site at $\mathbf{R}_{mn} = m\hat{\mathbf{a}}_1 + n \hat{\mathbf{a}}_2$. Considering a gauge configuration $\ket{\mathcal{G}_0}$ or $\ket{\mathcal{G}_1}$, such expectation can be computed as following
\begin{equation}\label{eq.52}
\braket{c_{\alpha}(\mathbf{R}_{mn}) c_{\alpha}(\mathbf{R}'_{m'n'})}
=
\big \langle \begin{pmatrix}
c_{\mathrm{A}} & c_{\mathrm{B}}
\end{pmatrix}
\mathcal{W}^{(\alpha)}_{\mathbf{R} \mathbf{R}'}
\begin{pmatrix}
c_{\mathrm{A}} \\
c_{\mathrm{B}}
\end{pmatrix}
\big \rangle
=
\sum_{ll'=1}^{2L^2}
\left( U^{\dagger} \mathsf{R}^{\dagger} \mathcal{W}^{(\alpha)}_{\mathbf{R} \mathbf{R}'} \mathsf{R} U \right)_{ll'} \braket{a^{\dagger}_l a_{l'}}
=
\sum_{l=L^2+1}^{2L^2}
\left( U^{\dagger} \mathsf{R}^{\dagger} \mathcal{W}^{(\alpha)}_{\mathbf{R} \mathbf{R}'} \mathsf{R} U \right)_{ll},
\end{equation}
where $\mathcal{W}^{(\alpha)}_{\mathbf{R} \mathbf{R}'}$ is $2L^2 \times 2L^2$ matrix defined as follows
\begin{equation}\label{eq.53}
\mathcal{W}^{(\alpha)}_{\mathbf{R} \mathbf{R}'}
=
\begin{pmatrix}
\delta_{\alpha,\mathrm{A}} \mathcal{B}^{(\alpha)}_{\mathbf{R} \mathbf{R}'}	& 				\mathcal{O}_{L^2}				\\
		\mathcal{O}_{L^2}					&	\delta_{\alpha,\mathrm{B}} \mathcal{B}^{(\alpha)}_{\mathbf{R} \mathbf{R}'}
\end{pmatrix}.
\end{equation}
Here, $\mathcal{B}^{(\alpha)}_{\mathbf{R} \mathbf{R}'}$ are the matrices corresponding to the non-zero connections allowed by the orientations of the triangles, and $\mathcal{O}_{L^2}$ is a null matrix of order $L^2 \times L^2$. We compute the above expectation value for both the gauge configurations $\ket{\mathcal{G}_0}$, and $\ket{\mathcal{G}_2}$. Depending on the number of the bond-fermions, we may have to consider the effects of non-trivial Majorana zero modes and the issue with the bond fermions. Following Ref.~\cite{PhysRevB.92.014403,PhysRevB.84.165414}, then we have to satisfy the parity constraint for the matter fermions as $(-1)^{n_f + n_{\chi}} = 1$, where $n_f$ ($n_{\chi}$) denotes the number of matter (bond) fermions. Depending on the situation, we may have to modify the above Eq.~\ref{eq.52} as follows
\begin{equation}\label{eq.54}
\braket{c_{\alpha}(\mathbf{R}_{mn}) c_{\alpha}(\mathbf{R}'_{m'n'})}
=
\sum_{l=L^2+2}^{2L^2}
\left( U^{\dagger} \mathsf{R}^{\dagger} \mathcal{W}^{(\alpha)}_{\mathbf{R} \mathbf{R}'} \mathsf{R} U \right)_{ll} +
\left( U^{\dagger} \mathsf{R}^{\dagger} \mathcal{W}^{(\alpha)}_{\mathbf{R} \mathbf{R}'} \mathsf{R} U \right)_{L^2L^2},
\end{equation}
to consider one matter fermion excitation. In the presence of MZM with two well-separated visons, we perform singular value decomposition (SVD) to diagonalize the Hamiltonian matrix in Eq.~(\ref{eq.46}). Summing over the triangles as mentioned earlier, we obtain the localized current profiles as illustrated in Fig.~2 and Fig.~3, in the main text. In the thermodynamic limit, when the visons are far away from each other, the current profile becomes identical around each hexagonal plaquette containing the vison excitations. Hence, we focus our analysis around a single vison excitation.

Furthermore, we note that in two-vison configuration $\ket{\mathcal{G}_2}$, the ground state energy of the Majorana fermions decreases in the presence of an external magnetic field. In Fig.~\ref{fig:Fig7}, we show the variation of ground state energy $E^{\kappa \neq 0}_2 - E^{\kappa = 0}_2$ for various strength and orientation of the magnetic field. As the energy difference decreases in the presence of the magnetic field, the vison configuration is favored in the presence of the external magnetic field. Hence, the magnetic field lowers the energy for the vison configurations and may lead to formation of exotic vison crystal phases~\cite{PhysRevLett.123.057201}. We leave such an explicit analysis for future work.

redpajama_set_name: RedPajamaArXiv