Anatomy of fermionic entanglement and criticality in Kitaev spin liquids

We analyze in detail the effect of nontrivial band topology on the area-law behavior of the entanglement entropy in Kitaev's honeycomb model. By mapping the translationally invariant 2D spin model onto 1D fermionic subsystems, we identify those subsystems responsible for universal entanglement contributions in the gapped phases and those responsible for critical entanglement scaling in the gapless phases. For the gapped phases, we analytically show how the topological edge states contribute to the entanglement entropy and provide a universal lower bound for it. For the gapless semimetallic phases and topological phase transitions, the identification of the critical subsystems shows that they fall always into the Ising or the XY universality classes. As our study concerns the fermionic degrees of freedom in the honeycomb model, qualitatively similar results are expected to apply also to generic topological insulators and superconductors.

Figure 1

2D system defined on a torus of size $L_x\times L_y$ is bi-partitioned in regions $A$ and $B$ such that translational symmetry along $L_y$ is preserved. If the band structure is characterized by a topological invariant, both the edge energy spectrum and the entanglement spectrum exhibit localized states along the cut $\partial A$.

Figure 2

(Top) The energy gap $G$ in the vortex-free $\left(\theta =0\right)$ and full vortex $\left(\theta =1\right)$ sectors. In the vortex-free sector, the phase transitions between the $\nu =0$ and the gapless phase for $K=0$ or the $\nu =1$ phase for $K>0$ occurs at $J_c=1/2$, while at the full vortex sector the transition to the $\nu =2$ phase or the corresponding gapless phase takes place at $J_c=\sqrt{\frac{1-K^2}{2}}$. (Inset) The transition between $\nu =1$ and $\nu =2$ topological phases occurs for staggered couplings corresponding to ${\theta }_c=\frac{3+K^2}{4}$. (Bottom) Divergence of the rate of change in the entropy correlates with the quantum phase transitions. Here, $L_x=L_y=144$ in the main figures and $L_x=L_y=100$ in the insets.

Figure 3

The correspondence between the physical velocities $v_E$ and the virtual velocities $v_{\lambda }$. The analytic expression (solid line) given by (16) for the physical velocities agrees with the numerical values (dots). Via the ansatz (17) (dashed line) it reproduces the numerically obtained virtual velocities (squares).

Figure 4

Entanglement spectrum (Left) and entropy (Right) in the momentum interval $\mathrm{\Delta }p$ (between vertical lines) in the vortex-free $\left(\theta =0\right)$ sector with $J=1$ and $K=0.12$. A linear approximation ${\lambda }^{\pm }$ with velocity $v_{\lambda }$ obtained form (17) contributes almost all of the entropy inside $\mathrm{\Delta }p$ (blue circles). The analytic expression (19) (red dashed) as a quadratic approximation to $S_{\text{edge}}\left(p\right)$ reproduces the entropy within accuracy of $0.1\%$ inside $\mathrm{\Delta }p$.

Figure 5

Area law $S=\alpha L_y$ obeyed by total entropy. The area-law coefficient $\alpha$, theoretically upper bounded by $ln2$, increases as the virtual velocity $v_{\lambda }$ decreases as a function of $J$ (insets). Here, $K=0.12$.

Figure 6

Energy spectrum $E\left(p\right)$ on a cylinder and the corresponding the entanglement spectrum $\lambda \left(p\right)$ and the entanglement dispersion $S\left(p\right)$ for the $\nu =1$ (left) and $\nu =2$ (right) phases. The real-space edge states are indicated in red. Here, $J=1$ and $K=0.15$.

Figure 7

The effective model for the edge entropy $S_{\text{edge}}$. (Top) Partitioning the system on a torus to $A$ and $B$ gives rise to two edges with counter-propagating Majorana edge states in both partitions. We take the edges to correspond to the sites of a four-site tight-binding model, with Majorana modes $a_L$ and $a_R$ living on the edges of $A$ and $b_L$ and $b_R$ living on the edges of $B$. These hybridize according to the Hamiltonian (21), where the solid (dashed) arrows correspond to strong (weak) tunneling of magnitude $1 \left(w_p={\kappa }_{\text{tm}}\left|E_{\pm }\left(p\right)\right|\right)$. The sign of the tunneling is positive along the arrows, which gives $\pi$ flux on the single plaquette. (Bottom) Entropy $S_{\text{tm}}$ (red dots) between the Fermi points obtained from the effective model reproduces accurately the entropy $S$ (black line) of the $\nu =1$ phase of the honeycomb model. The data are for $L_x=24,L_y=60,K=0.125$ and the value of the optimized proportionality constant ${\kappa }_{\text{tm}}$ as the function of $K$ is shown in the inset.

Figure 8

The energy spectrum $E\left(p\right)$, the entanglement spectrum $\lambda \left(p\right)$, and the entanglement dispersion $S\left(p\right)$ at critical points in the (a)–(c) vortex-free sector $\left(\theta =0\right)$ and (d)–(f) in the full vortex sector $\left(\theta =1\right)$. The transitions are (a) $\nu =0$ to Majorana semimetal with two Fermi points at $J=\frac{1}{2},K=0$, (b) $\nu =0$ to $\nu =1$ at $J=\frac{1}{2},K=0.15$, (c) $\nu =1$ to $\nu =-1$ at $J=1,K=0$, (d) $\nu =0$ to Majorana semimetal with four Fermi points at $J=\frac{1}{\sqrt{2}},K=0$, (e) $\nu =0$ to $\nu =2$ at $J=\sqrt{\frac{1-K^2}{2}},K=0.15$, and (f) $\nu =2$ to $\nu =-2$ at $J=1,K=0$.

Figure 9

Scaling of the entanglement entropy $S\left(p_c\right)$ of the gapless chains at the transitions $(\nu \leftrightarrow {\nu }^{\prime })$ between gapped topological phases characterized by Chern numbers $\nu$ and ${\nu }^{\prime }$ as the function of the chord length (24). For $(0\leftrightarrow 1)$ transition, $p_{\text{c}}=\pi$, for $(-1\leftrightarrow 1)$ transition, $p_{\text{c}}=\pm \frac{2\pi }{3}$, for $(0\leftrightarrow 2)$ transition, $p_{\text{c}}=0,\pi$, for $(1\leftrightarrow 2)$ transition, $p_{\text{c}}=0$, and for $(-2\leftrightarrow 2)$ transition, $p_{\text{c}}=\pm \frac{\pi }{6},\pm \frac{5\pi }{6}$. For each critical chain, the scaling implies a central charge of $c=1/2$.

Figure 10

Scaling of the entropy $S$ in the gapless phase in the vortex-free sector $\left(\theta =0,J=1\right)$ as the function of the partition length $x$ for cut lengths $L_y$ commensurate with the critical momenta $p_c=\frac{2\pi }{3}$. (Left) In the thin torus limit $L_y\ll L_x$ the quasi-2D models scale precisely with $\tilde{c}=1/2$, while in the regime $L_y\approx L_x$ one only finds scaling close to this lower bound in the large chord length limit. (Right) This system size dependence of the scaling can be traced back to quasicritical chains. In the thin-torus limit, there are only few chains with large gaps giving constant contributions, while in a large system some quasicritical chains show scaling all the way to the large $x\to L_x$ limit. Here $L_x=210$.

Figure 11

Composite entropy, $S_{\text{comp}}\left(p\right)$, as given by (27) in the gapless phase of the vortex-free sector. The flat band of edge states contributes $S_{\text{edge}}$ between the Fermi cones at $p_c=\pm \frac{2\pi }{3}$, the critical contributions $S_{\text{c}}$ only arise at $p_c$ and quasicritical chains contribute $S_{\text{qc}}$ outside them. The nonuniversal constant in Eqs. (24) and (25) is chosen such that $S_{\text{qc}}\left(0\right)=S\left(0\right)$. Here, $J=1$ and $L_x=L_y=300$. (Inset) Even if the system is critical, for a fixed partition size $x$ the system still obeys the area law. Here, $x=105$.

Figure 12

Anisotropic nearest-neighbor couplings $J_x,J_y,J_z$ (solid arrows) and chiral next-nearest-neighbor couplings $K$ (dashed arrows) of the honeycomb lattice model with lattice vectors $v_x$ and $v_y$ (dash-dotted arrows). The black and white sites denote the two triangular sublattices of the honeycomb lattice. In the Majorana picture, the tunneling is chosen to be positive along the arrows. In red, we denote every other $z$ link along rows of the honeycomb whose value is dictated by $\theta$, enabling the interpolation between the vortex-free and full vortex sectors.

Figure 13

The set of zeros of the denominator $\{z\in D(p\left)\right\}$ (blue) and the set of zeros of the numerator $\{z\in N_{\pm }\left(p\right)\}$ (red) of ${\tilde{G}}_{\pm }\left(z\right)$ on the complex plane (for $J=0.6$ and $K=0.15$). The unit circle $\left|z\right|=1$ is depicted in green. Simplification among zeros in the numerator and denominator of $\tilde{G}$ occurs when the red dots cover the blue lines. End modes are possible only if portions of both of the two connected blue lines not covered with red dots lie completely outside or inside the green circle. Consistently with Fig. 14, the former (latter) situation is possible only for $G_{+} \left(G_{-}\right)$ on the bottom (top) plot. This implies the existence of edge modes with dispersion $E_{+}$ at the boundary $j=1$ and $E_{-}$ at the boundary $j=N$.

Figure 14

Comparison between the absolute value of the zeros of the denominator ($\left|z\right|:z\in D\left(p\right)$, in blue) and the zeros of the numerator ($\left|z\right|:z\in N_{\pm }\left(p\right)$, in red) of ${\tilde{G}}_{\pm }$ as a function of $p$ in different regimes. From top to bottom, each row considers fixed parameters $J$ and $K$. In particular $J=0.6$ and $K=0.15$ for the first row (non-Abelian regime), $J=1/2$ and $K=0.15$ for the second row (phase transition), and $J=0.4$ and $K=0.15$ for the third row (Abelian regime). Within each row, on the left (right), we plot the zeros of the numerator associated with $G_{-} \left(G_{+}\right)$. In green we plot the line at $\left|z\right|=1$. For each plot, the blue lines not covered by red dots represent the poles of $G_{\pm }$ (when the red dots cover the blue line the zeros of the numerator simplify the zeros of the denominator). In the non-Abelian regime, a range of $p$ exists where all poles of $G_{+} \left(G_{-}\right)$ on the right (Left) plot have absolute value bigger (smaller) than 1 (allowing for end modes at $j=1$ (right) and $j=N$ (left), see also Figs. 13 and 15. In the Abelian case there is no $p$ where all poles have absolute value bigger or smaller than one. In this case, end modes do not arise [29].

Figure 15

Existence of end modes and validity of the edge-mode expression in Eq. (C13). Top row: zeros of the denominator ($\left|z\right|:z\in D\left(p\right)$, in blue) and the zeros of the numerator ($\left|z\right|:z\in N_{\pm }\left(p\right)$, in red) of ${\tilde{G}}_{+}$ as a function of $p$ for $J=0.6$ and $K=0.15$ (left) and for $J=1$ and $K=0.15$ (right). Bottom row: dispersion of the edge modes (C13) and their linear approximation (C16) each for the values of $J$ and $K$ corresponding to the upper row (same side). In dark grey, we highlight the momentum space region where edge modes are not present so that the validity of (C13) (bottom row) is limited to the white area. The dispersion with negative (positive) angular coefficient at $p=\pi$ corresponds to modes at the $j=1 \left(j=N\right)$ boundary.

Figure 16

Momentum space range $\mathrm{\Delta }p$ where edge modes are present (C23) as a function of $J$ (left) and $K$ (right). On the left different lines correspond to different values for $K$: dashed line $[K=0$, (C25)], green line $\left(K=0.05\right)$, orange line $\left(K=0.15\right)$, blue line $\left(K=0.3\right)$. On the right different lines correspond to different values of $J$: green line $\left(J=1/2\right)$, orange line $\left(J=0.51\right)$, blue line $\left(J=0.6\right)$. Both plots show that edge modes exist only in the non-Abelian phase, as expected. In fact, on the left we can see that $\mathrm{\Delta }p\to 0$ as $J\to {1/2}^{+}$ and on the right $\mathrm{\Delta }p=0$ for the $J=1/2$ line consistently with the plots in Fig. 14. On the plot on the left we can see that, for $K\to 0$ the behavior becomes that described in Eq. (C25) (dashed line).