Trigonal Warping in Bilayer Graphene: Energy versus Entanglement Spectrum
TTrigonal Warping in Bilayer Graphene: Energy versus Entanglement Spectrum
Sonja Predin, Paul Wenk, and John Schliemann
Institute for Theoretical Physics, University of Regensburg, D-93040 Regensburg, Germany (Dated: December 2015)We present a mainly analytical study of the entanglement spectrum of Bernal-stacked graphenebilayers in the presence of trigonal warping in the energy spectrum. Upon tracing out one layer,the entanglement spectrum shows qualitative geometric differences to the energy spectrum of agraphene monolayer. However, topological quantities such as Berry phase type contributions toChern numbers agree. The latter analysis involves not only the eigenvalues of the entanglementHamiltonian but also its eigenvectors.We also discuss the entanglement spectra resulting from tracing out other sublattices. As atechnical basis of our analysis we provide closed analytical expressions for the full eigensystem ofbilayer graphene in the entire Brillouin zone with a trigonally warped spectrum.
I. INTRODUCTION
Although first considered as a source of quantum cor-rections to the entropy of black holes , entanglement en-tropy, in particular the von Neumann entropy, evolvedinto a tool in the field of many-body systems. Thisbrought along connections between seemingly unrelatedresearch areas. In condensed matter, the entanglemententropy serves, e.g., as a geometrical interpretation forthe boundary between local quantum many-body sys-tems. This connection has its origin in the area laws .However, Li and Haldane have shown that the re-lated entanglement spectrum contains more informationthan the single number expressed by the entanglemententropy . This spectrum is determined by the Schmidtdecomposition of the ground state of a bipartite system,and the reduced density matrix obtained by tracing outone of the subsystems can always be formulated as ρ red = e −H ent Z (1)with an entanglement Hamiltonian H ent encoding theentanglement spectrum, and a partition function Z =tr( e −H ent ). Following the Li-Haldane conjecture , in agaped phase, the entanglement spectrum can be directlyrelated to the spectrum of edge excitations as shown forthe Fractional Quantum Hall System . This relationto the edge excitations can be also seen analytically inthe case of non-interacting particles. It can be shownby mapping the free fermionic system H onto a flat-band Hamiltonian H flat7 . Now, the eigenenergies e i of the lat-ter are related to the eigenenergies of the correspondingentanglement energies ε i as e i ∼ tanh ( ε i / / . Thereby, the eigenstates of both H flat and H are thesame. Thus, if H contains topologically protected surfacestates the same holds for the entanglement Hamiltonian.This is why the entanglement spectrum, beyond therelated entropy, is considered a tower of states and usedas a fingerprint for topological order. However, this isnot true in general as shown recently by A. Chandran etal. , Ref. 9.As a result of a multitude of studies, there is a plethoraof revisited effects in the context of entanglement spec- trum like the Kondo effect, many-body localization ordisordered quantum spin systems; for recent reviews seeRefs. 10 and 11.A particular situation arises if the edge comprises theentire remaining subsystem as it is the case for spinladders and various bilayer systems . A typicalobservation in such scenarios is, in the regime of stronglycoupled subsystems, a proportionality between the en-ergy Hamiltonian of the remaining subsystem and the ap-propriately defined entanglement Hamiltonian. We notethat the entanglement Hamiltonian entering the reduceddensity matrix (1) is only determined up to multiples ofthe unit operator which has consequences regarding ther-modynamic relations between the entanglement entropyand the subsystem energy .On the other hand, such a close relation between en-ergy and entanglement Hamiltonian is not truly generalas shown in Ref. 18 where a spin ladder of clearly non-identical legs was studied. In the present work we provideanother counter example given by graphene bilayers inthe presence of trigonal warping . As we shall see inthe following, the geometric properties of the the entan-glement spectrum of an undoped graphene bilayer andthe energy spectrum of a monolayer clearly differ quali-tatively. However, certain topological quantities such asBerry phase type contributions to Chern numbers agree.The latter analysis involves not only the eigenvalues ofthe entanglement Hamiltonian (i.e., the entanglementspectrum) but also its eigenvectors.This paper is organized as follows. In section II wediscuss the full eigensystem of the tight-binding modelof bilayer graphene in the presence of trigonal warping;a full account of the technical details is given in ap-pendices A and B. To enable analytical progress we ne-glect here terms breaking particle-hole symmetry. On theother hand, our calculation considers the entire first Bril-louin zone and avoids the Dirac cone approximation usu-ally employed in studies of trigonal warping in graphenebilayers . We compare our results for the full four-band model with an effective Hamiltonian acting on thetwo central bands . The entanglement spectrumobtained from the ground state of undoped graphene bi-layers is analyzed in section III. We discuss the case of a r X i v : . [ c ond - m a t . m e s - h a ll ] M a r one layer being traced out as well as the situation wherethe trace is performed over two other out of four sublat-tices. We close with a summary and an outlook in sectionIV. II. ENERGY SPECTRUM OF GRAPHENEBILAYERS: TRIGONAL WARPING ANDTOPOLOGICAL INVARIANTS
The standard tight-binding Hamiltonian for graphenebilayers in Bernal stacking can be formulated as H = − t (cid:88) (cid:126)k (cid:16) γ ( (cid:126)k ) a † (cid:126)k b (cid:126)k + γ ( (cid:126)k ) a † (cid:126)k b (cid:126)k + h . c . (cid:17) + t ⊥ (cid:88) (cid:126)k (cid:16) b † (cid:126)k a (cid:126)k + a † (cid:126)k b (cid:126)k (cid:17) − t (cid:88) (cid:126)k (cid:16) γ ( (cid:126)k ) b † (cid:126)k a (cid:126)k + γ ∗ ( (cid:126)k ) a † (cid:126)k b (cid:126)k (cid:17) + t (cid:88) (cid:126)k (cid:16) γ ( (cid:126)k ) (cid:16) a † (cid:126)k a (cid:126)k + b † (cid:126)k b (cid:126)k (cid:17) + h . c . (cid:17) , (2)where a † i(cid:126)k ( a i(cid:126)k ) and b † i(cid:126)k ( b i(cid:126)k ) create (annihilate) electronswith wave vector (cid:126)k in layers i = 1 , γ ( (cid:126)k ) = (cid:80) l =1 exp( i(cid:126)k · (cid:126)δ l )where the (cid:126)δ l are the vectors connecting a given carbonatom with its nearest neighbors on the other sublatticein a graphene monolayer. In what follows we will usecoordinates with (cid:126)δ , = a (cid:16) − , ±√ (cid:17) , (cid:126)δ = a (1 ,
0) (3)where a = 1 . (cid:126)K ± = 2 π √ a (cid:16) √ , ± (cid:17) . (4)The parameter t describes hopping within each layer be-tween the sublattices while t ⊥ parameterizes the verticalhopping between the two sublattices in different layerslying on top of each other. The additional hopping pro-cesses described by the skew parameters t , t lead totrigonal warping of the spectrum and electron-hole asym-metry, respectively. Experimentally established values for these quantities are t = 3 . t ⊥ = 0 . t = 0 . t = 0 . | γ ( (cid:126)k ) | .The presence of all four couplings in the HamiltonianEq. (2) makes its explicit diagonalization in terms of an-alytical expressions a particularly cumbersome task. Asthe present study chiefly relies on analytical calculationsrather than resorts to numerics, we will drop the con-tributions proportional to the smallest parameter t inorder to achieve an analytically manageable situation. FIG. 1. (Color Online) Brillouin zone with a density plot of | γ ( (cid:126)k ) | .FIG. 2. (Color Online) Contour plot of the energy band(+ E ( (cid:126)k )) plotted for t ⊥ = t , t = 0 . t . The contour of thecolored region indicates E = 0 . /t ⊥ . The edge of the firstBrillouin zone is marked by dashed lines. Putting t = 0 the full eigensystem of the Hamilto-nian (2) can be obtained in a closed analytical fashionas detailed in appendix A. The four dispersion branches( ± E ( (cid:126)k )), ( ± E ( (cid:126)k )) form a symmetric spectrum with E / = (cid:118)(cid:117)(cid:117)(cid:116) (cid:32) t ⊥ + t | γ ( (cid:126)k ) | + 2 t | γ ( (cid:126)k ) | ± (cid:114) t | γ ( (cid:126)k ) | (cid:16) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1)(cid:17) + (cid:16) t ⊥ − t | γ ( (cid:126)k ) | (cid:17) (cid:33) (5) and γ ( (cid:126)k ) = | γ ( (cid:126)k ) | e iφ (cid:126)k . The two outer branches ( ± E ( (cid:126)k ))are separated from the inner ones ( ± E ( (cid:126)k )) by gaps de-termined essentially by the hopping parameter t ⊥ . Theresult Eq. (5) generalizes the energy spectrum given inRef. 27 within the Dirac cone approximation to the fullBrillouin zone. Moreover, in appendix A we also give thecomplete data of the corresponding eigenvectors. Fig. 3concentrates on the vicinity of a given K -point using re-alistic parameters.The inner branches ( ± E ( (cid:126)k )) dominate the low-energyphysics of the system near half filling and meet at zeroenergy for γ ( (cid:126)k ) = 0 (6)corresponding to the two inequivalent corners K ± of thefirst Brillouin zone, and forcos (cid:0) φ (cid:126)k (cid:1) = − ∧ | γ ( (cid:126)k ) | = t ⊥ t t . (7)The latter condition defines three additional satelliteDirac cones around each K -point two of which lying onthe edges (faces) of the Brillouin zone connecting K ± .The third satellite Dirac cone lies formally outside theBrillouin zone but is equivalent to a satellite cone on theedge around an equivalent K -point. Indeed, the quan-tity γ ( (cid:126)k ) has a constant phase φ (cid:126)k ∈ {− π/ , π/ , π } oneach face: As an example, consider the edge connectingthe two inequivalent K -points given in Eq. (4) where onefinds γ (cid:18) π a , k y (cid:19) = e − iπ/ (cid:32) (cid:32) √ k y a (cid:33) − (cid:33) (8)with the parenthesis being nonnegative for k y rangingbetween ( ± π/ (3 √ a )). Thus, solving for k y the satelliteDirac cones on that edge lie at (cid:126)k = (cid:18) π a , ± √ a arccos (cid:18) (cid:18) t ⊥ t t (cid:19)(cid:19)(cid:19) , (9)and the other satellite cones are located at positions be-ing equivalent under reciprocal lattice translation and/orhexagonal rotation. Note that for t ⊥ t /t = 1 the satel-lite cones merge in the M -points (centers of the faces) andthey vanish for even larger values of that ratio. In Fig. 2we give a sketch of the situation in the entire Brillouinzone for moderate values of t ⊥ and t . For t = 0 thetwo energy bands ( ± E ( (cid:126)k )) touch only at the K -pointswhere they have a quadratic dispersion. Finite t (cid:54) = 0causes a splitting into in total four Dirac cones with lin-ear dispersion, an effect known as trigonal warping .As a further important property, the eigenvectors cor-responding to ( ± E ( (cid:126)k )) are discontinuous as a func-tion of wave vector at the degeneracy points defined byEq. (7); for more technical details we refer to appendixB. As a simplistic toy model mimicking such an effectone can consider the Hamiltonian H = − kσ z with a FIG. 3. (Color Online) The central energy bands ( ± E ( (cid:126)k ))plotted around a given K -point for t ⊥ = 0 . t , t = 0 . t . Thedispersions show a central Dirac cone accompanied by threesatellites. The components of the wave vector are measuredrelatively to the K -point. one-dimensional wave number k and the Pauli matrix σ z describing some internal degree of freedom: In themany-body ground state of zero Fermi energy all occu-pied states with k > k < k = 0. As weshall see below, in the present case of graphene bilayersthis discontinuity is also reflected in the entanglementspectrum.An effective Hamiltonian providing an approximate de-scription of the central bands ( ± E ( (cid:126)k )) can be given fol-lowing Ref. 27. In up to linear order in 1 /t ⊥ one finds H = − t t ⊥ (cid:16) γ ∗ ( (cid:126)k ) (cid:17) + t γ ( (cid:126)k ) t t ⊥ (cid:16) γ ( (cid:126)k ) (cid:17) + t γ ∗ ( (cid:126)k ) 0 (10) with respect to the basis (cid:16) b † (cid:126)k , a † (cid:126)k (cid:17) | (cid:105) . The eigenstatesread | χ ± (cid:105) = 1 √ (cid:18) ∓ e iψ (cid:126)k (cid:19) (11)with e iψ (cid:126)k = t t ⊥ (cid:16) γ ( (cid:126)k ) (cid:17) + t γ ∗ ( (cid:126)k ) (cid:12)(cid:12)(cid:12)(cid:12) t t ⊥ (cid:16) γ ( (cid:126)k ) (cid:17) + t γ ∗ ( (cid:126)k ) (cid:12)(cid:12)(cid:12)(cid:12) . (12)Note that the Hamiltonian (10) vanishes if and only ifthe conditions (6) or (7) are fulfilled implying that thepositions of the central and satellite Dirac cones are thesame as for the full Hamiltonian (2). Moreover , ψ (cid:126)k isa smooth and well-defined function of the wave vectorexcept for the locations of Dirac cones. Accordingly, theBerry curvature F ( (cid:126)k ) = ∂A y ∂k x − ∂A x ∂k y (13)arising from the Berry connection (cid:126)A ( (cid:126)k ) = i (cid:104) χ ± ( (cid:126)k ) | ∂∂(cid:126)k | χ ± ( (cid:126)k ) (cid:105) = − ∂ψ (cid:126)k ∂(cid:126)k (14)vanishes everywhere outside the Dirac cones where con-tributions in terms of δ -functions arise. Integrating theBerry connection along closed path in (cid:126)k -space leads togeometrical quantities often referred to as Berry phases,although no contact to adiabaticity is made here. More-over, if the Berry curvature has only nonzero contribu-tions in terms of δ -functions (as it is the case here and inthe following) these geometrical phases are indeed topo-logical, i.e. they are invariant under continuous varia-tions of the paths as long as the support of the δ -functionsis not touched.As discussed in Refs. 31, 33, and 34, integrating alonga closed path around the central Dirac cones at K ± yieldsa Berry phase of ( ∓ π ), while each of the accompanyingsatellite cones gives a contribution of ( ± π ). Thus, thetotal Berry phase arising at and around each K -pointis, as in the absence of trigonal warping, ( ± π ), and theintegral over the whole Brillouin zone of the Berry con-nection (i.e. the Chern number) vanishes. Naturally, ourpresent analysis going beyond the Dirac cone approxima-tion confirms these results. III. ENTANGLEMENT SPECTRA
For systems of free fermions as studied here, the en-tanglement Hamiltonian can be formulated as a single-particle operator , H ent = (cid:88) λ ξ λ d † λ d λ . (15)Here the d † λ generate eigenstates of the correlation matrix C αβ = (cid:104) Ψ | c † α c β | Ψ (cid:105) , (16)where | Ψ (cid:105) is the ground state of the composite system,and single-particle operators c α , c β act on its remaining part after tracing out a subsystem. The entanglementlevels ξ λ are related to the eigenvalues η λ of the correla-tion matrix via ξ λ = ln (cid:18) − η λ η λ (cid:19) = 2 artanh (1 − η λ ) . (17)In particular, the entanglement Hamiltonian and the cor-relation matrix share the same system of eigenvectors. A. Tracing out One Layer
We now consider the ground state of the undopedgraphene bilayer such that all states with negative en-ergies ( − E ( (cid:126)k )), ( − E ( (cid:126)k )) are occupied while all othersare empty. Tracing out layer 1 leads to the correlationmatrix C ( (cid:126)k ) = (cid:18) u ( (cid:126)k ) u ∗ ( (cid:126)k ) (cid:19) (18)where an explicit expression for u ( (cid:126)k ) is given in appendixC. The entanglement levels corresponding to the eigen-values η ± ( (cid:126)k ) = 1 / ∓ | u ( (cid:126)k ) | are ξ ± ( (cid:126)k ) = ± (cid:16) | u ( (cid:126)k ) | (cid:17) . (19)The modulus | u | can be formulated as | u | = 1 / (cid:113) d/ ( t | γ ( (cid:126)k ) | )) (cid:115) (cid:18) − (cid:15) (cid:15) + b E E (cid:19) (20)with (cf. Eqs. (A14),(A15)) d = (cid:16) t ⊥ − t | γ ( (cid:126)k ) | (cid:17) / (cid:113) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) , (21) b = t ⊥ t | γ ( (cid:126)k ) || sin (cid:0) φ (cid:126)k (cid:1) | (cid:113) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) , (22)and (cf. Eq. (A21)) (cid:15) , = t | γ ( (cid:126)k ) |± (cid:114)(cid:16) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1)(cid:17) / d (23) implying E , = (cid:113) (cid:15) , + b . (24)The r.h.s of Eq. (20) becomes zero if the radicand van-ishes. According to the discussion in appendices B andC this is the case when cos (cid:0) φ (cid:126)k (cid:1) = − b = 0and E = (cid:15) ≥ E = | (cid:15) | such that | u | ∝ (cid:115) (cid:18) − (cid:15) | (cid:15) | (cid:19) (25) FIG. 4. (Color Online) Contour plot of the entanglementspectrum ξ + ( (cid:126)k ) plotted for t ⊥ = t , t = 0 . t . The contour ofthe colored region indicates ξ = 1 .
5. The dashed line delin-eates the first Brillouin zone.
Now equation (B2) shows that | u ( (cid:126)k ) | = 0 is equivalent tocos (cid:0) φ (cid:126)k (cid:1) = − ∧ | γ ( (cid:126)k ) | ∈ (cid:2) , t ⊥ t /t (cid:3) , (26)where the endpoint of the above interval defines accord-ing to condition (7) the location of the satellite Diraccones. As a result, the entanglement levels (19) vanishalong segments of the faces of the first Brillouin zonebounded by the positions of the central Dirac cones andtheir satellites. At the satellite Dirac cones the entan-glement spectrum is discontinuous as a function of wavevector. In Fig. 4 we plotted the entanglement spectrum ξ + ( (cid:126)k ) for the whole Brillouin zone. For a better visual-ization large hopping parameters have been chosen. Thecontour of the colored region connects all three satelliteDirac cones. As discussed in appendix B, this disconti-nuity is inherited from a discontinuity in the eigenvec-tors of the occupied single-particle states. The entangle-ment spectrum in the entire Brillouin zone is illustratedin Fig. 4, whereas Fig. 5 focuses on a given K -point.Moreover, apart from the eigenvalues of the entan-glement Hamiltonian, let us also consider its eigenvec-tor which coincide with the eigenvectors of the correla-tion matrix (18). As discussed in appendix C, the com-plex function u ( (cid:126)k ) entering the correlation matrix be-comes singular at the K -points and the positions of theaccompanying satellite Dirac cones of the energy spec-trum, leading again to δ -function-type contributions tothe Berry curvature which vanishes otherwise. Combin-ing symbolic computer algebra techniques and numericalcalculations we find here a Berry phase of ( ∓ π/
2) aroundthe corners K ± of the Brillouin zone, and ( ± π/
2) for thecorresponding satellite positions. For the central posi-tions the above calculations can also be done fully ana-
FIG. 5. (Color Online) The entanglement spectrum (19) plot-ted around a given K -point for the same parameters as inFig. 3. The density plot shows the upper entanglement level.Zero eigenvalues of the entanglement Hamiltonian occur alonglines connecting the K -point with the locations of satelliteDirac cones of the energy spectrum (thick black lines). Thecomponents of the wave vector are measured relatively to the K -point. lytically by expanding the eigensystem data around K ± .For the satellite locations such an expansion is not pos-sible due to the discontinuity of the eigenvectors.Thus, the total Berry phase contribution from each K -point K ± is ( ± π ) and agrees with the Berry phase aroundthe Dirac cones in monolayer graphene. As a result, al-though the entanglement spectrum of graphene bilayersgenerated by tracing out one layer shows obvious dif-ferences to the energy spectrum of monolayer grapheneregarding qualitative geometrical properties, the topolog-ical Berry phases obtained from the corresponding eigen-vectors still coincide at each K -point. B. Tracing out other Sublattices
Now, we will consider the entanglement spectrum ob-tained by tracing out sublattices A1 and B2 (or A2 andB1) lying in different layers. In the former case one finds C ( (cid:126)k ) = (cid:18) v ( (cid:126)k ) v ∗ ( (cid:126)k ) (cid:19) (27)where an explicit expression for v ( (cid:126)k ) is given in appendixC. The above correlation matrix has eigenvalues η ± ( (cid:126)k ) =1 / ∓ | v ( (cid:126)k ) | leading to the entanglement levels ξ ± ( (cid:126)k ) = ± (cid:16) | v ( (cid:126)k ) | (cid:17) . (28)In Fig. 6 we plotted the eigenvalues η − ( (cid:126)k ) = 1 / | v ( (cid:126)k ) | of the correlation matrix around a given K -point. Themodulus | v ( (cid:126)k ) | reads more explicitly | v ( (cid:126)k ) | = 12 (cid:115) − t | γ ( (cid:126)k ) | t | γ ( (cid:126)k ) | + d (cid:18) − (cid:15) (cid:15) + b E E (cid:19) (29)= 12 (cid:113) − | u ( (cid:126)k ) | (30)and has a similar structure as | u ( (cid:126)k ) | given in Eq. (20).In particular, | v ( (cid:126)k ) | = 1 / ⇔ | u ( (cid:126)k ) | = 0 if the conditions(26) are fulfilled. In this case η + = 0 and η − = 1 indi-cating that the remaining subsystem is unentangled withthe system traced out.Regarding Berry phases generated from the eigen-states of the correlation matrix (27) we note that theoff-diagonal element v ( (cid:126)k ) nowhere vanishes. As a con-sequence the Berry curvature defined analogously as inEqs. (11)-(14) is zero throughout the Brillouin zone,which in turn holds for all Berry phases. The nonvan-ishing of v ( (cid:126)k ) follows from the fact that | v ( (cid:126)k ) | = 0 wouldrequire | u ( (cid:126)k ) | = 1 / C ( (cid:126)k ) = (cid:18) (cid:19) , (31)indicating that these sublattices are maximally entangledwith the part traced out. IV. CONCLUSIONS AND OUTLOOK
We have studied entanglement properties of the groundstate of Bernal stacked graphene bilayers in the pres-ence of trigonal warping. Our analysis includes boththe eigenvalues of the reduced density matrix (givingrise to the entanglement spectrum) as well as its eigen-vectors. When tracing out one layer, the entanglementspectrum shows qualitative geometric differences to theenergy spectrum of a graphene monolayer while topolog-ical quantities such as Berry phase type contributions toChern numbers agree. The latter finding is in contrastto the reduced density matrix resulting from tracing out other sublattices of the bilayer system. Here, all corre-sponding Berry phase integrals yield trivially zero. Thus,our study provides an example for common topologicalproperties of the eigensystem of the energy Hamiltonian
FIG. 6. (Color Online) Eigenvalues η − ( (cid:126)k ) = 1 / | v ( (cid:126)k ) | of the correlation matrix plotted around a given K -point for t ⊥ = 0 . t , t = 0 . t . The thick black lines correspond tothe one in Fig. 5, and the components of the wave vector areagain measured relatively to the K -point. of a subsystem (here a graphene monolayer) and the en-tanglement Hamiltonian, while the geometrical shape ofboth spectra grossly differs. Our investigations are basedon closed analytical expressions for the full eigensystemof bilayer graphene in the entire Brillouin zone with atrigonally warped spectrum.Future work might address bilayer systems of other ge-ometrical structures such as the Kagome lattice, the in-fluence of a static perpendicular magnetic field , andthe effect of time-periodic in-plane electric fields . Note added.
After this paper was made available asan arXiv preprint and submitted to the journal, we be-came aware of Ref. where also Chern numbers calcu-lated from the eigenstates of entanglement Hamiltoniansare studied. Most recent work building upon this conceptis reported on in Ref. . ACKNOWLEDGMENTS
This work was supported by Deutsche Forschungsge-meinschaft via GRK1570.
Appendix A: Diagonalization of the Bilayer Hamiltonian
Putting t = 0 and fixing a wave vector (cid:126)k the Hamiltonian (2) reads with respect to the basis (cid:16) a † (cid:126)k , b † (cid:126)k , b † (cid:126)k , a † (cid:126)k (cid:17) | (cid:105) H = t ⊥ − tγ ( (cid:126)k ) 0 t ⊥ − tγ ∗ ( (cid:126)k ) − tγ ∗ ( (cid:126)k ) 0 0 − t γ ( (cid:126)k )0 − tγ ( (cid:126)k ) − t γ ∗ ( (cid:126)k ) 0 . (A1)Using γ ( (cid:126)k ) = | γ ( (cid:126)k ) | e iφ (cid:126)k we apply the transformation U = 1 √ e iφ (cid:126)k e − iφ (cid:126)k e iφ (cid:126)k − e − iφ (cid:126)k − (A2)such that in H = U HU † = t ⊥ − t | γ ( (cid:126)k ) | − t | γ ( (cid:126)k ) | − t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) it | γ ( (cid:126)k ) | sin (cid:0) φ (cid:126)k (cid:1) − it | γ ( (cid:126)k ) | sin (cid:0) φ (cid:126)k (cid:1) t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) − t | γ ( (cid:126)k ) | − t | γ ( (cid:126)k ) | − t ⊥ (A3)all information on the phase φ (cid:126)k is contained in the matrix elements being proportional to the skew parameter t .Proceeding now with the transformation U = 1 √ − −
10 0 1 1 (A4)we find H = U H U † = 12 e c − is − isc e is isis − is − e cis − is c − e (A5)with e = 2 t | γ ( (cid:126)k ) | + t ⊥ − t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) , (A6) e = − t | γ ( (cid:126)k ) | + t ⊥ − t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) , (A7) c = t ⊥ + t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) , (A8) s = t | γ ( (cid:126)k ) | sin (cid:0) φ (cid:126)k (cid:1) . (A9)Here it is useful to split the above matrix as H = H (cid:48) + H (cid:48)(cid:48) where H (cid:48) = 12 e − is e isis − e − is − e , H (cid:48)(cid:48) = 12 c − isc is − is cis c . (A10) H (cid:48) is diagonalized by U = α + − iσα − − iσα + α − − iσα − α + α − − iσα + (A11)with σ = sign (cid:16) sin(3 φ ( (cid:126)k )) (cid:17) and α ± = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) ± t ⊥ − t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1)(cid:113) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) (A12)such that H = U H U † = ζ idσ b − idσ ζ b b − ζ idσb − idσ − ζ (A13)where d = (cid:16) t ⊥ − t | γ ( (cid:126)k ) | (cid:17) / (cid:113) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) , (A14) b = t ⊥ t | γ ( (cid:126)k ) || sin (cid:0) φ (cid:126)k (cid:1) | (cid:113) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) , (A15)and ± ζ and ± ζ are eigenvalues of H (cid:48) given by ζ / = 12 (cid:18) ± t | γ ( (cid:126)k ) | + (cid:113) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1)(cid:19) . (A16)Splitting now H in the form H = ζ id − id ζ − ζ id − id − ζ + b b b b (A17)the first part is diagonalized by U = − iστ β + β − β − − iστ β + − iστ β + β − β − − iστ β + (A18)with τ = sign( d ) and β ± = (cid:118)(cid:117)(cid:117)(cid:116) (cid:32) ± ζ − ζ (cid:112) ( ζ − ζ ) + 4 d (cid:33) (A19)while the second part is left unchanged by U resulting in H = U H U † = (cid:15) b (cid:15) b b − (cid:15) b − (cid:15) (A20)with the diagonal elements are given in terms of (cid:15) / = 12 (cid:18) ζ + ζ ± (cid:113) ( ζ − ζ ) + 4 d (cid:19) . (A21)Finally, H is brought into diagonal form via U = γ (1)+ γ (1) − γ (2)+ γ (2) − γ (2) − − γ (2)+ γ (1) − − γ (1)+ (A22)with γ (1) ± = (cid:115) (cid:18) ± (cid:15) E (cid:19) , γ (2) ± = (cid:115) (cid:18) ± (cid:15) E (cid:19) (A23)and E / = (cid:113) (cid:15) , + b (A24)= (cid:118)(cid:117)(cid:117)(cid:116) (cid:32) t ⊥ + t | γ ( (cid:126)k ) | + 2 t | γ ( (cid:126)k ) | ± (cid:114) t | γ ( (cid:126)k ) | (cid:16) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1)(cid:17) + (cid:16) t ⊥ − t | γ ( (cid:126)k ) | (cid:17) (cid:33) . (A25)Thus, U H U † = diag ( E , E , − E , − E ) , (A26)and the matrix elements of the corresponding total transformation U = U U U U U can be expressed as U = 12 ( α − − iσα + ) ( τ β + + β − ) (cid:16) γ (1)+ − iσγ (1) − (cid:17) (A27) U = 12 ( α + − iσα − ) ( τ β + + β − ) (cid:16) γ (1) − − iσγ (1)+ (cid:17) (A28) U = − e iφ (cid:126)k α − − iσα + ) ( τ β + − β − ) (cid:16) γ (1)+ + iσγ (1) − (cid:17) (A29) U = e − iφ (cid:126)k α − + iσα + ) ( τ β + − β − ) (cid:16) γ (1)+ − iσγ (1) − (cid:17) (A30)and U = −
12 ( α + + iσα − ) ( τ β + − β − ) (cid:16) γ (2)+ − iσγ (2) − (cid:17) (A31) U = −
12 ( α + − iσα − ) ( τ β + − β − ) (cid:16) γ (2)+ + iσγ (2) − (cid:17) (A32) U = − e iφ (cid:126)k α + + iσα − ) ( τ β + + β − ) (cid:16) γ (2)+ + iσγ (2) − (cid:17) (A33) U = − e − iφ (cid:126)k α + − iσα − ) ( τ β + + β − ) (cid:16) γ (2)+ − iσγ (2) − (cid:17) (A34)which are the complex conjugates of the components of the eigenvectors of the conduction-band states with positiveenergies E ( (cid:126)k ), E ( (cid:126)k ), while U = 12 ( α − − iσα + ) ( τ β + − β − ) (cid:16) γ (2)+ − iσγ (2) − (cid:17) (A35) U = 12 ( α − + iσα + ) ( τ β + − β − ) (cid:16) γ (2)+ + iσγ (2) − (cid:17) (A36) U = − e iφ (cid:126)k α + + iσα − ) ( τ β + + β − ) (cid:16) γ (2) − − iσγ (2)+ (cid:17) (A37) U = − e − iφ (cid:126)k α + − iσα − ) ( τ β + + β − ) (cid:16) γ (2) − + iσγ (2)+ (cid:17) (A38)0and U = 12 ( α + + iσα − ) ( τ β + + β − ) (cid:16) γ (1)+ − iσγ (1) − (cid:17) (A39) U = −
12 ( α − + iσα + ) ( τ β + + β − ) (cid:16) γ (1) − − iσγ (1)+ (cid:17) (A40) U = e iφ (cid:126)k α + + iσα − ) ( τ β + − β − ) (cid:16) γ (1)+ + iσγ (1) − (cid:17) (A41) U = e − iφ (cid:126)k α − + iσα + ) ( τ β + − β − ) (cid:16) γ (1) − + iσγ (1)+ (cid:17) (A42)correspond to the valence-band states with negative energies ( − E ( (cid:126)k )), ( − E ( (cid:126)k )). Note that all factors involving α ± , γ (1) ± , γ (2) ± in the above expressions have modulus one, i.e. they are phase factors. Appendix B: Continuity Properties
The eigenvectors corresponding to the energy branches ( ± E ( (cid:126)k )) are discontinuous at wave vectors determined bythe condition (7). This comes about as follows: The matrix elements U ,n ( (cid:126)k ), U ,n ( (cid:126)k ), n ∈ { , , , } contain thequantities γ (2) ± defined in Eqs. (A23) whereas the U ,n ( (cid:126)k ), U ,n ( (cid:126)k ) corresponding to ( ± E ( (cid:126)k )) involve γ (1) ± Fixing nowcos (cid:0) φ (cid:126)k (cid:1) = − b = 0 such that E = (cid:15) ≥ E = | (cid:15) | such that γ (1) ± remain continuous while γ (2) ± become γ (2) ± = (cid:115) (cid:18) ± (cid:15) | (cid:15) | (cid:19) . (B1)Inspection of Eq. (A21) now shows that for cos (cid:0) φ (cid:126)k (cid:1) = − (cid:15) ( (cid:126)k ) (cid:26) > | γ ( (cid:126)k ) | < t ⊥ t /t < | γ ( (cid:126)k ) | > t ⊥ t /t (B2)such that (cid:15) ( (cid:126)k ) changes sign for | γ ( (cid:126)k ) | = t ⊥ t /t , i.e. γ (2) ± is discontinuous at wave vectors given by the condition (7).This discontinuity is inherited by the correlation matrix and, in turn, by the entanglement spectrum.The technical reason for this discontinuity in the eigenvectors is of course the fact that the dispersions ( ± E ( (cid:126)k ))become degenerate at wave vectors fulfilling (7). In fact the eigenvectors can also be considered as continuous functionsof the wave vector by appropriately relabeling the dispersion branches. In the ground state of the undoped bilayersystem, however, only the lower branch ( − E ( (cid:126)k )) is occupied, which makes the discontinuity unavoidable.To circumvent this discontinuity one can open an energy gap between the upper and lower central band suchthat the corresponding eigenstates are necessarily continuous for all wave vectors. Among the various mechanismsproducing such a gap only few allow for a still halfway convenient analytical treatment of the Hamiltonian. Theseinclude introducing identical mass terms in both layers, i.e. H (cid:55)→ H + H (cid:48) with H (cid:48) = diag ( m, − m, − m, m ) , (B3)or applying a bias voltage Λ between the layers, H (cid:48) = diag ( − Λ / , Λ / , − Λ / , Λ / . (B4)In the former case the four dispersion branches ( ± E ( (cid:126)k )), ( ± E ( (cid:126)k )) are given by E / ( (cid:126)k ) = (cid:34) m + 12 (cid:16) t ⊥ + t | γ ( (cid:126)k ) | + 2 t | γ ( (cid:126)k ) | (cid:17) ± (cid:114) t | γ ( (cid:126)k ) | (cid:16) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1)(cid:17) + (cid:16) t ⊥ − t | γ ( (cid:126)k ) | (cid:17) (cid:35) / (B5)while for a bias voltage one finds E / ( (cid:126)k ) = (cid:34) Λ (cid:16) t ⊥ + t | γ ( (cid:126)k ) | + 2 t | γ ( (cid:126)k ) | (cid:17) ± (cid:114) t | γ ( (cid:126)k ) | (cid:16) t ⊥ + t | γ ( (cid:126)k ) | − t ⊥ t | γ ( (cid:126)k ) | cos (cid:0) φ (cid:126)k (cid:1) + Λ (cid:17) + (cid:16) t ⊥ − t | γ ( (cid:126)k ) | (cid:17) (cid:35) / . (B6)In both cases the central energy bands ( ± E ( (cid:126)k )) are separated by a gap, and the spectrum can still be given in terms ofcomparably simple closed expressions since the characteristic polynomial of the 4 × .Note that applying a bias voltage as well as introducing a mass term in each layer discriminates the layers againsteach other. The latter circumstance is due to the fact that t ⊥ couples sublattices in different layers for which the massterm has different sign. As a result, when tracing out, say, one layer of an undoped (i.e. half-filled) bilayer system, theremaining layer will not be half-filled, what obscures somewhat the comparison with an undoped graphene monolayer. Appendix C: Correlation Matrices
Upon tracing out layer 1 from the ground state of the undoped bilayer system the correlation matrix reads in thebasis (cid:16) a † (cid:126)k , b † (cid:126)k (cid:17) | (cid:105) C ( (cid:126)k ) = (cid:18) U U ∗ + U U ∗ U U ∗ + U U ∗ U U ∗ + U U ∗ U U ∗ + U U ∗ (cid:19) = (cid:18) u ( (cid:126)k ) u ∗ ( (cid:126)k ) (cid:19) (C1)with u ( (cid:126)k ) = e − iφ (cid:126)k (cid:0) β − β − (cid:1) (cid:18)(cid:16) γ (1)+ − iσγ (1) − (cid:17) − (cid:16) γ (2)+ − iσγ (2) − (cid:17) (cid:19) . (C2)This quantity becomes singular at the corners of the Brillouin zone where γ ( (cid:126)k ) is zero such that its phase is ill-defined,and at the positions of the satellite Dirac cones of the energy spectrum where, as discussed in appendix B, γ (2) ± isdiscontinuous.Tracing out the sublattices A1 and B2 one finds in the basis (cid:16) a † (cid:126)k , b † (cid:126)k (cid:17) | (cid:105) C ( (cid:126)k ) = (cid:18) U U ∗ + U U ∗ U U ∗ + U U ∗ U U ∗ + U U ∗ U U ∗ + U U ∗ (cid:19) = (cid:18) v ( (cid:126)k ) v ∗ ( (cid:126)k ) (cid:19) (C3)with v ( (cid:126)k ) = ( α − − iσα + ) (cid:18) ( τ β + − β − ) (cid:16) γ (2)+ − iσγ (2) − (cid:17) + ( τ β + + β − ) (cid:16) γ (1)+ − iσγ (1) − (cid:17) (cid:19) . (C4)Note that the expressions (C2),(C4) obey the interesting sum rule | u ( (cid:126)k ) | + | v ( (cid:126)k ) | = 14 (C5)which s fulfilled whenever the coefficients involved satisfy α + α − = β + β − = (cid:16) γ (1 / (cid:17) + (cid:16) γ (1 / − (cid:17) = 1 , (C6)which is the case here by construction.Finally, the correlation matrix obtained by tracing out the sublattices A1, A2 is proportional to the unit matrix, C ( (cid:126)k ) = (cid:18) U U ∗ + U U ∗ U U ∗ + U U ∗ U U ∗ + U U ∗ U U ∗ + U U ∗ (cid:19) = (cid:18) (cid:19) (C7)implying that the remaining subsystem is maximally entangled with the subsystem traced out. L. Bombelli, R. K. Koul, J. Lee, and R. D. Sorkin, Phys.Rev. D , 373 (1986). J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. , 277 (2010). H. Li and F. D. M. Haldane, Phys. Rev. Lett. , 010504(2008). A. M. L¨auchli, E. J. Bergholtz, J. Suorsa, and M. Haque,Phys. Rev. Lett. , 156404 (2010). R. Thomale, D. P. Arovas, B. A. Bernevig, Phys. Rev. Lett , 116805 (2010). A. Chandran, M. Hermanns, N. Regnault, and B. A.Bernevig, Phys. Rev. B , 205136 (2011). A. Kitaev, Ann. Phys. , 2 (2006). A. M. Turner, Y. Zhang, A. Vishwanath, Phys. Rev. B, , 241102 (2010). A. Chandran, V. Khemani, and S. L. Sondhi, Phys. Rev.Lett. , 060501 (2014). N. Regnault, arXiv:1510.07670. N. Laflorencie, arXiv:1512.03388. D. Poilblanc, Phys. Rev. Lett. , 077202 (2010). J. I. Cirac, D. Poilblanc, N. Schuch, and F. Verstraete,Phys. Rev. B I. Peschel and M.-C. Chung, Europhys. Lett. , 50006(2011). A. M. L¨auchli and J. Schliemann, Phys. Rev. B , 054403(2012). J. Schliemann and A. M. L¨auchli, J. Stat. Mech. (2012)P11021. S. Tanaka, R. Tamura, and H. Katsura, Phys. Rev. A ,032326 (2012). R. Lundgren, V. Chua, and G. A. Fiete, Phys. Rev. B ,224422 (2012). R. Lundgren, Y. Fuji, S. Furukawa, and M. Oshikawa,Phys. Rev. B , 245137 (2013). X. Chen and E. Fradkin, J. Stat. Mech. (2013) P08013. R. Lundgren, arXiv:1412.8612. J. Schliemann, Phys. Rev. B , 115322 (2011). J. Schliemann, New J. Phys. , 053017 (2013), Corrigen-dum: , 079501 (2013). J. Schliemann, J. Stat. Mech. (2014)
P09011. E. McCann and M. Koshino, Rep. Prog. Phys. , 056503(2013). A. V. Rozhkov, A. O. Sboychakov, A. L. Rakhmanov, andF. Nori, arXiv:1511:06706. E. McCann and V. I. Falko, Phys. Rev. Lett. , 086805(2006). J. Nilsson, A. H. Castro Neto, N. M. R. Peres, and F.Guinea, Phys. Rev. B , 244418 (2006). M. Koshino and T. Ando, Phys. Rev. B , 245403 (2006). K. Kechedzhi, V. I. Falko, E. McCann, and B. L. Altshuler,Phys. Rev. Lett. , 176806 (2007). J. L. Manes, F. Guinea, and M. A. H. Vozmediano, Phys.Rev. B , 155424 (2007). J. Cserti, A. Csorda, and G. David, Phys. Rev. Lett. ,066802 (2007). G. P. Mikitik and Y. V. Sharlai, Phys. Rev. B , 113407(2008). E. Mariani, A. J. Pearce, and F. von Oppen, Phys. Rev. B , 165448 (2012). D. A. Cosma and V. I. Falko, Phys. Rev. B , 165412(2015). A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blakeand K. S. Novoselov, Phys. Rev. B , 165406 (2009). I. Peschel, J. Phys. A: Math. Gen. , L205 (2003). S.-A. Cheong and C. L. Henley, Phys. Rev. B , 075111(2004). N. Nemec, and G. Cuniberti, Phys. Rev. B , 201404(2007). Y.-X. Wang and F. Li, arXiv:1511.01972. S. Predin, P. Wenk, and J. Schliemann, unpublished. T. Fukui and Y. Hatsugai, J. Phys. Soc. Jpn. , 113705(2014).43