Entanglement and Disordered-Enhanced Topological Phase in the Kitaev Chain ^†

Abstract

In recent years, tools from quantum information theory have become indispensable in characterizing many-body systems. In this work, we employ measures of entanglement to study the interplay between disorder and the topological phase in 1D systems of the Kitaev type, which can host Majorana end modes at their edges. We find that the entanglement entropy may actually increase as a result of disorder, and identify the origin of this behavior in the appearance of an infinite-disorder critical point. We also employ the entanglement spectrum to accurately determine the phase diagram of the system, and find that disorder may enhance the topological phase, and lead to the appearance of Majorana zero modes in systems whose clean version is trivial.

1. Introduction

Entanglement is the quintessential characteristic of the quantum world. Recently, much attention has been devoted to manifestations of entanglement in quantum many-body systems [1,2]. Motivations include the prospect of using many-body ground and excited states as a resource for quantum communication and computation; the amount of entanglement controls the applicability of a powerful tensor-network-based numerical method; entanglement is an indicator of quantum correlations which is independent of the details of the system, allowing, in particular, one to identify quantum phase transitions even without knowing the relevant order parameter. In this work, we will concentrate on the latter property, in the context of disordered 1D quantum many-body systems.

Majorana zero modes [3,4,5,6] and their more complicated relatives [7] have been studied extensively recently, both due to their predicted exotic non-abelian braiding statistics, which gives rise to the prospect of topologically-protected quantum computation [8], and to many concrete proposals for their realization in the laboratory. Leading candidates are semiconductor quantum wires with strong spin orbit interaction, which are rendered effectively spinless by the application of an appropriate magnetic field and gate voltage, and driven into a topological phase by proximity-coupling to a superconductor, leading to the formation of Majorana end modes [9,10]. Indications for these modes have recently been measured experimentally [11,12,13,14,15]. These Majorana zero modes show up not only in the presence of a real edge, but also in the entanglement spectrum, that is, in the spectrum of the reduced density matrix of a subsystem [16].

Disorder naturally occurs in all these systems, and may hamper their topological characteristics, especially the Majorana edge modes of 1D systems [17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49]. In this work, we will study the interplay of disorder and the topological phase from the entanglement point of view. After introducing the model and our method for calculating the entanglement spectrum and entropy in Section 2, we will examine the entanglement entropy in Section 3. We will show that although superconducting proximity and disorder each tend to separately suppress the entanglement, their interplay may enhance the entanglement and cause it to behave in a non-monotonic fashion. We will explain this as the result of a strong-disorder quantum criticality. In Section 4, we will use the entanglement spectrum, and in particular, the presence of entanglement Majorana zero modes, to distinguish between the two phases (topological and trivial). We will show that in certain regimes of the parameter, space disorder can actually enhance or even be the origin of the topological phase in the system. We will summarize our findings and give a future outlook in Section 5.

2. Model and Method

We will study the standard Kitaev chain [3] model with disorder. The model describes spinless fermions hopping on a 1D lattice in the presence of (proximity-induced) pairing and disorder potential,

$$H_{\mathrm{Kitaev}}=\sum _{n=1}^L(V_n-\mu )c_n^{†}{c}_n-t{c}_{n+1}^{†}{c}_n-\Delta c_{n+1}{c}_n+\mathrm{h}.\mathrm{c}.,$$

(1)

where $c_n^{†}$ creates a fermion (with the standard anticommutation relations) at site $n=1\dots L$ (and $L+1\equiv 1$ if the boundary conditions are periodic; for open boundary conditions, terms which refer to site $L+1$ are omitted), $V_n$ is the disorder potential, independently and uniformly distributed in the interval $[-W/2,W/2]$, $\mu$ is the chemical potential, and t and $\Delta$ are, respectively, the hopping matrix element and the pairing amplitude, which one may choose as real without loss of generality.

Since the model is quadratic, it can be solved using a Bogolubov transformation. First, let us define the Nambu–Gorkov Fermi operators, ${\psi }_i$ ($i=1\cdots 2L$) as ${\psi }_n=c_i$ and ${\psi }_{n+L}=c_n^{†}$ for $n=1\cdots L$. Then, the Hamiltonian has the general Bogolubov–de Gennes form

$$H_{\mathrm{BdG}}=\sum _{i,j=1}^{2L}{h}_{ij}{\psi }_i^{†}{\psi }_j,$$

(2)

where in our specific case the only nonzero elements of the $2L\times 2L$ matrix $h_{ij}$ are $h_{nn}=-h_{n+L,n+L}=v_n-\mu$, $h_{n,n+1}=h_{n+1,n}=-h_{n+L,n+L+1}=-h_{n+L+1,n+L}=-t$ and $h_{n,n+L+1}=-h_{n+1,n+L}=-h_{n+L,n+1}=h_{n+L+1,n}=-\Delta$. The problem then reduces to the diagonalization of the “single-particle Hamiltonian” matrix $h_{ij}$: Let us denote by $U_{ij}$ the unitary matrix whose jth column contains the jth normalized eigenvector of the matrix $h_{ij}$, corresponding to the real eigenvalue ${\lambda }_j$ (which are guaranteed to come in pairs of equal magnitude and opposite sign), so that ${\sum }_{k=1}^{2L}{h}_{ik}{U}_{kj}=U_{kj}{\lambda }_j$ ($i,j,k=1\cdots 2L$). Defining ${\tilde{\psi }}_i={\sum }_{j=1}^{2L}{U}_{ji}^{\ast }{\psi }_j$ (so that ${\psi }_i={\sum }_{j=1}^{2L}{U}_{ij}{\tilde{\psi }}_j$), one has $H={\sum }_{i=1}^{2L}{\lambda }_i{\tilde{\psi }}_i^{†}{\tilde{\psi }}_i$.

In the clean case ($V_n=0$), the model is easily solved analytically [3]. It has two phases, depending on the ratio $|\mu /t|$: For $|\mu /t|>2$, it is a topologically trivial insulator (qualitatively similar to the limit $|\mu /t|\to \infty$, where the system is either completely filled or completely empty, depending on the sign of $\mu$). For $|\mu /t|<2$, the system is a topological superconductor, which features a single Majorana zero mode exponentially localized at each physical boundary (hence their hybridization is exponentially small in the system size), with the localization length inversely proportional to the energy gap. At $|\mu /t|=2$, the gap closes and the system features a 1D Majorana mode which allows the boundary zero modes to merge and disappear.

Let us now return to the general case, and discuss the calculation of observables in the many-body ground state of the model. In this state (to be denoted by $|{\Phi }_0\rangle$), all the modes with ${\lambda }_i<0$ are occupied. Therefore, $\langle {\tilde{\psi }}_i^{†}{\tilde{\psi }}_j\rangle ={\delta }_{ij}\theta (-{\lambda }_j)$, where $\theta$ is Heaviside’s step function, and hence $c_{ij}=\langle {\psi }_i^{†}{\psi }_j\rangle ={\sum }_k{U}_{ik}^{\ast }{U}_{kj}\theta (-{\lambda }_k)$. Our main interest is in the entanglement between a spatial subsystem A (defined by some subset of sites $A\subset \{1,\cdots ,L\}$ of size $\left|A\right|$) and its complement $\overline{A}$, when the total system is in its ground state $|{\Phi }_0\rangle$. In this case, the entanglement can be determined from the mixedness of its reduced many-body density matrix, ${\rho }_A={\mathrm{Tr}}_{\overline{A}}|{\Phi }_0\rangle \langle {\Phi }_0|$. Since the system is quadratic, ${\rho }_A$ is Gaussian, i.e., it can be written as ${\rho }_A\propto exp(-H_A)$ with a quadratic “entanglement Hamiltonian” [50]

$$H_A=\sum _{i,j\in A\cup A+L}{h}_{A,ij}{\psi }_i^{†}{\psi }_j,$$

(3)

where the dimensions of the matrix $h_{A,ij}$ are $2\left|A\right|\times 2\left|A\right|$. By Wick’s theorem, the latter is completely characterized by the single-time two-point fermionic correlation function, $c_{A,ij}=\langle {\psi }_i^{†}{\psi }_j\rangle$, which is simply the submatrix of $c_{ij}$ of size $2\left|A\right|\times 2\left|A\right|$, with $i,j\in A\cup A+L$ [50]. The latter is related to the matrix $h_A$ by $h_A=ln[(1-c_A)/c_A]$. This allows one to extract the entire spectrum of $h_A$ or ${\rho }_A$, as well as its moments. In particular, the von Neumann entanglement entropy, $S_A=-\mathrm{Tr}({\rho }_{A}ln{\rho }_A)=S_{\overline{A}}$, can be expressed in terms of the $2\left|A\right|$ eigenvalues $f_{A,i}$ of the matrix $c_A$ as [50]

$$S_A=-\sum _i{f}_{A,i}ln{f}_{A,i}.$$

(4)

3. Entanglement Entropy

We will start by studying the von Neumann entanglement entropy of a subsystem A which consists of the sites $n=1,\cdots ,L_A$ (a single continuous interval of length $L_A<L/2$). Generally in 1D, if the system is critical (gapless) and not localized by disorder, one expects the entanglement entropy to scale logarithmically with the subsystem size [51,52,53,54],

$$S_A\sim \frac{c}{3}ln(L_A)+\mathrm{const}.,$$

(5)

where c is the central charge of the corresponding conformal field theory, which roughly counts the number of gapless 1D modes in the system. For example, $c=1/2$ for a non-chiral Majorana (real) fermion mode and $c=1$ for a non-chiral Dirac (complex) fermion, or, equivalently, a non-chiral gapless bosonic mode. The above expression is valid when the interval in question has two boundaries with the rest of the system; for open boundary conditions with subsystem A located at one end, the coefficient of the logarithm is halved. The same logarithmic scaling holds as a function of the total system size L if the ratio $L_A/L$ is kept fixed. In the presence of a gap, this logarithmic scaling holds up to the correlation length (inversely proportional to the gap), above which the entanglement entropy saturates. Similar behavior holds for a gapless disordered localized system, with the role of correlation length being taken by the localization length.

Let us examine the behavior of the entanglement entropy in our system as a function of the total system size L. We will henceforth consider only periodic boundary conditions; open boundaries were found to give similar effects (up to the above-mentioned factor of 2). In Figure 1, we concentrate on $\mu =0$ and start from the clean case ($W=0$) without superconductivity ($\Delta =0$). Then, we have just a half-filled band of free fermions, and one finds the expected logarithmic scaling, $S_A\sim (c/3)ln(L)$ with $c=1$ (gapless Dirac fermion). Adding either superconductivity ($\Delta \ne 0$) or disorder ($W>0$) separately makes the entanglement saturate as a function of L, due to either the finite correlation length (inversely proportional to the gap $\Delta$) or localization length (which decreases as W increases). One would then expect that if one introduces both superconductivity and disorder together, the entanglement entropy will saturate even more quickly. Interestingly, this is not the case: as can be seen in Figure 2a, if we fix $\Delta$ to a finite value and increase W from zero, the entanglement starts to increase at some value of W, reaching a logarithmic behavior at $W/t\approx 7.82$ (the way to accurately determine this value will be discussed in the next section), before decreasing again as W is further increased.

To understand the origin of this peculiar behavior, it is better to employ a dual point of view. Let us concentrate on the case $\Delta =t$. Using a Jordan–Wigner transformation to define standard spin-1/2 operators at each site n (presented for open boundary conditions, but could be adapted to periodic ones as well),

$$\begin{array}{cc}\hfill S_n^x& =c_n^{†}{c}_n-1/2,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill S_n^z& =(c_n^{†}+c_n)\prod _{n^{\prime }=1}^{n-1}(2c_{n^{\prime }}^{†}{c}_{n^{\prime }}-1),\hfill \end{array}$$

(7)

the Hamiltonian (1) can be mapped into that of the disordered transverse-field 1D Ising model,

$$H_{\mathrm{Ising}}=\sum _{n=1}^L{h}_n{S}_n^x-J_n{S}_n^z{S}_{n+1}^z,$$

(8)

where $h_n=(V_n-\mu )/2$, $J_n=t=\Delta$. This model has a ${\mathbb{Z}}_2$ symmetry $S_n^z\to -S_n^z$, corresponding to the overall fermion parity in the Kitaev chain. In the absence of disorder ($W=0$), this Hamiltonian has two phases: when $\left|h\right|>\left|J\right|$, the Zeeman term dominates, resulting in a paramagnetic phase; when $\left|h\right|<\left|J\right|$, the exchange term dominates, resulting in a broken-symmetry ferromagnetic phase (for $J>0$). The non-local Jordan–Wigner transformation maps these into the (symmetry-unbroken) trivial and topological superconductor phase of the Kitaev chain, respectively. However, now we can understand that the same phases exist even when $J_n$ and/or $h_n$ become random: when the former dominates, a ferromagnet still results (if the sign of $J_n$ is fixed, it is a usual ferromagnet; if not, one may flip $S_{n^{\prime }}^z$ between each pair of negative $J_n$ to make all the $J_n$ positive, so the original system is a ferromagnet in disguise), while when the latter dominates a paramagnet appears [55,56,57,58].

More interesting is the transition between the two phases in the disordered case, which presents a new type of quantum criticality, known as the infinite-disorder critical point [55,56,57,58]. Its physics has been analyzed using the strong disorder renormalization group. In the latter, at each step, one picks the largest coupling in the chain. If this is a magnetic field $h_n$, the corresponding spin is polarized according to the sign of $h_n$; its neighboring spins are coupled by a new exchange term, $\sim J_{n-1}{J}_n/h_n$. If, on the other hand, the largest coupling is an exchange term $J_n$, the spins coupled by it become restricted to the two-dimensional low energy space $|\Uparrow \rangle \equiv |\uparrow \uparrow \rangle$ and $|\Downarrow \rangle \equiv |\downarrow \downarrow \rangle$ (assuming $J_n>0$), which can be thought of as a new effective spin-1/2. This effective spin-1/2 experiences an effective magnetic field $\sim h_n{h}_{n+1}/J_n$. If, at a later stage of the renormalization, this spin becomes polarized in the x direction, a Bell state is formed in terms of the original spins. Counting how many such Bell states straddle the boundary between subsystem A and the rest of the system, one finds that the entanglement entropy features a logarithmic behavior, as in Equation (5), but with an irrational coefficient, $c_{\mathrm{eff}}=ln2/2$ [59,60]. This is exactly the behavior that we observed in Figure 2a at $W/t\approx 7.82$.

This also has a natural interpretation from the original fermion point of view [17,18,24]. In the clean topological phase, Majorana zero modes exist at the ends of the system. Due to the gap, their coupling and hence splitting are exponentially-small in the distance between them, that is, the system size. When disorder is introduced, bulk modes start filling in the gap. At some critical value of the disorder, these gap modes approach zero energy, creating a channel which allows the Majoranas at the ends to couple and eliminate each other, making the system trivial. This transition point corresponds exactly to the strong-disorder critical point in the spin chain language. The fermionic picture also shows that the transition should persist when $\Delta \ne t$, where the equivalent spin chain becomes more complicated than the disordered transverse-field Ising model, Equation (6). This is supported by our results in Figure 2, where $\Delta /t=0.6\ne 1$.

Let us continue examining our data. New intriguing behavior arises when one introduces a finite chemical potential $\mu$. Figure 2b depicts the behavior for $\mu /t=2$. Here, in the clean limit, $W=0$, the system is at the transition from a topological superconductor to a trivial state, featuring a logarithmic behavior described by Equation (5) with $c=1/2$, due to the Majorana nature of the critical (gapless) bulk mode. One would then naïvly expect the disorder (which suppressed the topological phase at $\mu =0$) to immediately drive the system into the trivial phase. However, in reality, while the entanglement entropy initially drops with disorder, it then increases again to a logarithmic behavior with a $c_{\mathrm{eff}}=ln2/2$ coefficient at $W/t\approx 6.21$, before dropping back again. This shows that, in this case, disorder initially drives the system into the topological phase, and only later on makes it trivial, via a strong-disorder critical point. This behavior persists at $\mu /t=2.2$ (Figure 2c), where the clean system is trivial, and increasing disorder causes the entanglement to grow to a logarithmic behavior at $W/t\approx 2.39$, then drop, increase again to a logarithmic behavior at $W/t\approx 5.62$, then drop again. In both cases, $c_{\mathrm{eff}}=ln2/2$. Thus, here, disorder causes two transitions: from trivial to topological and back. Finally, at $\mu /t=2.5$, logarithmic behavior is never obtained, indicating that the system stays trivial for all values of W (Figure 2d). All of this prompts one to study more carefully the phase diagram of the system, which will be the topic of the next section.

4. Phase Diagram from Entanglement

The non-monotonic changes in the subsystem size dependence of the entanglement entropy as disorder is varied, discussed in the previous section, call for a more careful study of the phase diagram of the system. However, the $S_A(L)$ curves are not a very accurate means to probe the phase boundaries, since for reasonable system sizes it is not easy to distinguish numerically between a critical logarithmic behavior, Equation (5), and a saturation at some constant value.

To overcome this difficulty, we employ a new strategy, looking at the entanglement spectrum, that is, the spectrum of the single-particle entanglement Hamiltonian, Equation (3). In the topological phase, it should feature Majorana zero modes at the ends of the subsystem, just like a system with open ends, whereas in the trivial phase there should be no such zero modes. This is exemplified in Figure 3a, where the median of the gap ${\delta }_A$ between the two lowest eigenvalues of the entanglement Hamiltonian (3) is plotted as a function of disorder for different system sizes. One can clearly distinguish between fast decay to zero for small $W/t$ and apparent saturation for large $W/t$, as the subsystem size $L_A$ is increased. To make this quantitative, for each disorder strength we fit the dependence of ${\delta }_A$ on the subsystem size to:

$${\delta }_A(L_A)={\delta }_{\infty }+{\delta }_0{e}^{-L_A/{\xi }_{\delta }}.$$

(9)

We would like to identify whether the constant term ${\delta }_{\infty }$ is finite or zero. For this purpose, it is very important to use the correct form for the second, length-dependent, term in the case where the constant ${\delta }_{\infty }$ is zero, but not when ${\delta }_{\infty }$ is finite. Hence, the choice of exponential dependence of the second term, appropriate for the topological phase, is sufficient for our purposes (compare with Ref. [61]). The resulting ${\delta }_{\infty }$ is plotted as a function of disorder in Figure 3b, and shows nicely the topological-to-trivial transition. The error bars for ${\delta }_{\infty }$ come from the error in determining the median ${\delta }_A$, which in turn was found using the bootstrap method.

We can now determine the phase transition point as the largest W value where ${\delta }_{\infty }$ most probably vanishes. By this, we mean the value of W at which the lower end of the $2\sigma$ confidence interval for ${\delta }_{\infty }$ approaches zero. This is determined by linear interpolation between the two probed values of W between which this lower end of the confidence interval changes from positive to negative, as indicated in the inset to Figure 3b. We can also define the largest W which is certainly in the topological phase (${\delta }_{\infty }=0$) as the one for which the error in ${\delta }_{\infty }$ is one order of magnitude larger than the estimated value of ${\delta }_{\infty }$. The distance between this point and the transition point will serve to set the error bar for the transition point.

Using this procedure, we determined the phase diagram of our system in the disorder strength-chemical potential plane for several values of $\Delta /t$. The results are displayed in Figure 4. At the transition line, the entanglement entropy should display the critical logarithmic behavior, Equation (5) with $c_{\mathrm{eff}}=ln(2)/2$, except at the clean point ($W=0$ and $\mu /t=2$), where $c=1/2$. This is exactly the behavior found earlier in Figure 2 for several values of $\mu$ at $\Delta /t=0.6$. Let us note that the phase diagram of the system for $\Delta =t$ has been previously determined in reference [39], by looking for the presence of long range order in the $S^z$ correlation function in the spin version of the model, Equation (6), and the results agree with ours. However, the possibility to identify the topological phase by appealing to some spin correlation function is very model-specific, while our method applies generally. Let us also note that reference [44] studied the phase diagram of a tight-binding chain with a large unit cell in which pairing centers are distributed in an ordered or a disordered fashion. This allowed us to characterize the system using the Kitaev’s topological index. However, this cannot be easily generalized to either the current case of a fully-disordered chain, or to the interacting regime. In contrast, our entanglement-based approach applies generically.

The most striking feature of the phase diagram is the existence of a topological region at $\mu /t>2$ for intermediate values of W, as anticipated in the previous section from the behavior of the entanglement entropy. In this region, the clean system is non-topological, and the topological phase is brought about by the disorder. Indeed, in the regime of large chemical potential, the clean system is almost completely filled with fermions, hence trivial. With disorder, regions of positive potential fluctuations can emerge, where a high density of holes may exist. These can create a network throughout the system, allowing for a topological phase with Majorana zero modes at its ends, as evidenced by our results.

5. Conclusions

To conclude, in this work we have examined the interplay of a topological phase and disorder in the archetypical example of the random Kitaev chain, from the point of view of entanglement. We have found the length dependence of the entanglement entropy to change non-monotonically with disorder, indicating strong-disorder quantum criticality of the random transverse field Ising type, which persists even for $\Delta \ne t$, where the corresponding spin model is more complicated. We then developed a method to accurately identify the transition by looking at the entanglement spectrum and extracting from it the transition point and its uncertainty. Using this method, we verified the existence of a region in phase space where the topological phase is enabled by disorder.

In the future, it would be interesting to see how our results extend to interacting systems [38,39,62,63,64,65,66] (where the density matrix renormalization group [67,68] and related methods give direct access to the many-body density matrix and its spectrum), as well as to systems in higher dimensions [69,70].