Dynamical structure factor in the non-Abelian phase of the Kitaev honeycomb model in the presence of quenched disorder

Kitaev's model of spins interacting on a honeycomb lattice describes a quantum spin liquid, where an emergent static ${\mathbb{Z}}_2$ gauge field is coupled to Majorana fermions. In the presence of an external magnetic field and for a range of interaction strengths, the system behaves as a gapped, non-Abelian quantum spin liquid. In this phase, the vortex excitations of the emergent ${\mathbb{Z}}_2$ gauge field have Majorana zero modes bound to them. Motivated by recent experimental progress in measuring and characterizing real materials that could exhibit spin-liquid behavior, we analytically calculate the dynamical spin structure factor in the non-Abelian phase of the Kitaev's honeycomb model. In particular, we treat the case of quenched disorder in the vortex configurations. Our calculations reveal a peak in the low-energy dynamical structure factor that is a signature of the spin-liquid behavior. We map the effective Hamiltonian to that of a chiral $p$-wave superconductor by using the Jordan-Wigner transformation. Subsequently, we analytically calculate the wave functions of the Majorana zero modes, the energy splitting for finite separation of the vortices and, finally, the dynamical structure factor in presence of quenched disorder.

Figure 1

(a) A plaquette of the honeycomb lattice. The letters $x$, $y$, and $z$ label the type of the link between sites and indicate if the corresponding coupling in the Hamiltonian $H_0$ is of a $xx$, $yy$, or $zz$ type. The direction of the links represents an ordering of the operators in the Hamiltonian. They do not have a meaning for the spin Hamiltonian in Eq. (1), but are important for the Majorana hopping Hamiltonian in Eq. (6). The three labels $s$, $s^{\prime }$, and $s^{″}$ represent an example configuration that contributes to the magnetic perturbation $H_K$. The dashed line indicates how the three neighboring sites partly encircle a plaquette as mentioned in the main text. (b) The plot shows the next-to-nearest-neighbor links (NNN links) $x_{\perp }$, $y_{\perp }$, and $z_{\perp }$ for a site on the black and one on the white sublattice. These links are perpendicular to the original links and point in opposite directions for the white and the black sublattices. The NNN links describe the hopping in $H_M$ of Eq. (6). (c) The honeycomb lattice with the choice of the unit cell (gray) and the basis vectors ${\mathit{n}}_1$ and ${\mathit{n}}_2$. Each unit cell contains a white and black dot connected with a vertical $z$ link. The lattice constant between sites is $a$.

Figure 2

The dashed line shows the Jordan-Wigner string (JWS) used in Eq. (3). It runs in zigzag along a horizontal line from the left to the right. It jumps from the right end of the system back to the left where it runs again from left to right on the horizontal line above, until it reaches the site $s$.

Figure 3

A configuration with a vortex in the middle of the shaded plaquette marked with the cross. We choose the following gauge to describe a vortex: all ${\alpha }_l$ on the dashed line are set to $-1$ while all other ${\alpha }_l$ are kept to be 1. In a continuum model, the dashed line will manifest itself as a branch cut. The wave function of bound Majorana modes changes sign across the branch cut.

Figure 4

The plot shows the analytic result for the wave function of the zero mode (dots) on top of the numerical exact wave function (little crosses) obtained by numerical diagonalization of Eq. (6). The wave functions are evaluated along the branch cut on the white sublattice (with $K=0.03J$, for a vortex at the origin). The dots and crosses represent the wave function on lattice sites while the lines connecting the dots are a guide to the eye. The plot shows that the analytic results are valid even for relatively small distances from the vortex core.

Figure 5

The sketch shows the geometry we choose for the calculation of the energy splitting due to hybridization of two MZMs. The gray plaquettes with a cross mark the position of a vortex. The angle ${\phi }_1$ of the vortex further right on the $x$ axis is defined with respect to an axis that runs through the vortex parallel to the $x$ axis such that ${\phi }_1\in [0,2\pi [$. This adds a branch cut for the wave function pointing along the $x$ axis to positive infinity. The vortex further left on the $x$ axis is described by the angle ${\phi }_2$ with respect to an axis that runs through this vortex parallel to the $x$ axis. This time ${\phi }_2\in [-\pi ,\pi [$, so that the branch cut points along the $x$ axis to negative infinity. We can choose to evaluate the matrix element $\langle {\mathbf{\chi }}_2\left|{\mathcal{H}}_2\right|{\mathbf{\chi }}_1\rangle$ on one of the branch cuts. In the main text, we choose to evaluate it on branch cut 1, indicated by the red color.

Figure 6

The upper panel shows the numerical results of Eq. (47) for the dynamical subgap structure factor normalized to the value at $\omega -E_{\text{fl}}=0$. The result is a narrow peak located around the flux gap $E_{\text{fl}}$. The lower panel shows a closeup of the rectangle in the upper panel. The dashed line represents our analytic result ${\overline{S}}_{\mathit{q}=0,\text{tail}}\left(\omega \right)$, while the solid line depicts the numerical result. We find an unconventional decay proportional to ${ln}^2\{{\widehat{\omega }}^2/{[\omega -E_{\text{fl}}]}^{2}ln[{\widehat{\omega }}^2/{(\omega -E_{\text{fl}})}^2]\}$.