{'fraction_non_alphanumeric': 0.040634291377601585, 'fraction_numerical': 0.014780023268841298, 'mean_word_length': 4.4774604673117775, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The scaling of entanglement entropy in a honeycomb lattice on a torus

The entanglement entropy of a noninteracting fermionic system confined to a two-dimensional honeycomb lattice on a torus is calculated. We find that the entanglement entropy can characterize Lifshitz phase transitions without a local order parameter. In the noncritical phase and critical phase with a nodal Fermi surface, the entanglement entropy satisfies an area law. The leading subarea term is a constant in the gapped phase rather than a logarithmic additive term in the gapless regime. The tuning of chemical potential allows for a nonzero Fermi surface, whose variation along a particular direction determines a logarithmic violation of the area law. We perform the scaling of entanglement entropy numerically and find agreement between the analytic and numerical results. Furthermore, we clearly show that an entanglement spectrum is equivalent to an edge spectrum.


Introduction
Entanglement entropy is a unifying theme in the interdisciplinary theories during the past decade. It is a helpful bridge between quantum gravity, quantum information and condensed matter physics. It is defined as the von Neumann entropy of one of the reduced states for a pure state. Given a normalized wave function |Ψ AB and a partition of the system into subsystems A and B, the associated reduced density matrix of A part is obtained from ρ A = Tr B |Ψ AB Ψ AB | by tracing out the degrees of freedom in B part. The bipartite entanglement entropy is given by, Related generalized quantities are the Rényi entropies defined as whose limit recovers the von Neumann entropy via relation S vN = lim n→1 S n . Once the complete values of S n for all n are known, one can obtain all eigenvalues of ρ A , which are called the entanglement spectrum. The spooky nonlocal nature makes the entanglement entropy hold great advantage in the context of condensed matter physics. On one hand, it plays an essential role in determining the efficiency of numerical algorithms based on matrix product states, as well as an extension to tensor network states [1] and Projected Entangled Pair States (PEPS) [2]. On the other hand, it is proven that the entanglement entropy is able to capture the quantum criticality in a many-body Hamiltonian [3,4]. The constituents become extraordinarily entangled at criticality. Quantum critical points are successfully pinpointed in both free systems [5,6] and interacting systems [7,8,9,10,11], especially for those critical ground states without prior knowledge of order parameters, e.g., disordered systems [12,13,14] and quantum spin liquids [15]. In particular, the discovery of a topological phase transition, where distinct phases of matter have the same symmetry, was in urgent need of basis-independent way of identifying quantum critical points. In this respect, extensive work has utilized the phase sensitivity of the entanglement entropy to successfully characterize the topological order. The universal constant term in the bipartite entanglement entropy has been shown to be a probe of topological contribution [16,17]. The success of entanglement entropy as a diagnostic tool is essentially rooted in its peculiar scaling properties in a quantum state. The growing awareness of the monotonic scaling of von Neumann entropy can be traced back to the study of black hole physics [18,19,20], and then the profound relations to conformal field theory (CFT) in (1+1)-dimensions were established [21,22]. Despite the formal simplicity, it is surprisingly difficult to calculate the entanglement entropy analytically, even for noninteracting models. The entanglement entropy in onedimensional system has been extensively studied and the scalings have been explicitly established by CFT [3,23]. Intensive research in terms of numerical algorithms unveiled a few qualitative and quantitative aspects in higher dimensions [24,25].
In this work, we study the half-space entanglement entropy of free fermion in the anisotropic honeycomb lattice with fine tuning chemical potential. In order to discriminate the relation between various boundaries and entanglement entropy, we place the fermion system on a torus. The scaling behavior of the entanglement entropy in a two-dimensional hexagonal lattice demands a comprehensive investigation. The basic interest in honeycomb lattices comes from the graphene and associated Dirac fermion physics. More interestingly, a recent influx of ideas and toolboxes can create the so-called "artificial honeycomb lattice" [26], such as photonic lattices [27], coupled molecular systems [28] as well as optical honeycomb lattice [29]. The manipulations of Dirac points are possible through these artificial methods [30,31,32]. The asymmetric tunneling amplitudes cause a topological Lifshitz transition from a semimetal to an insulator, which is characterized by the topological change of the Fermi surface in the Brillouin zone [33]. It was found that the Lifshitz points can be also captured by the entanglement entropy [34]. Moreover, based on the bulk-edge correspondence [35,36], the entanglement entropy and the Berry phase are closely related in a topologically nontrivial system through the zero-energy edge states [37]. A fascinating aspect of graphene-like structure is that there are two types of edges, i.e., zigzag and armchair edges. The role of edges in determining the entanglement entropy in honeycomb lattice and associated scaling behavior are not fully understood. In addition, the chemical potential can be well tuned in graphene by a gate voltage, so the dependence of the entanglement entropy on the chemical potential is a nontrivial issue to be discussed.

Tight binding model on honeycomb lattice
The anisotropic fermionic tight-binding model on a two-dimensional honeycomb lattice reads where a r (b r ) is the annihilation operator on the triangular sublattice Λ A • (respectively Λ B •) of the rth dimer, τ i is the nearest neighbor bond vector, and µ is the chemical potential. For simplicity, we confine the hexagonal lattice to the brick lattice without changing lattice topology as shown in Fig.1(a). The nearest-neighbor separations are given by where the bond length is set to unit. Taking Fourier transformation for all lattice sites, the Hamiltonian (3) reduces to and then the eigenspectrum is where Fig.1(b). The bulk band structure in Eq.(6) is invariant under inversion : k → − k. In the parameter space |t a − t b | < t c < t a + t b (containing symmetric case t a = t b = t c ), the valence and conduction bands intersect linearly at two k D -points (known as the Dirac points), given by which will be suited at K-points, i.e., (± 2 3 π, 0), of the first Brillouin zone for symmetric hoppings t 1 = t 2 = t 3 . The spectra and the corresponding Dirac points for a semimetallic phase is displayed in Fig.2(a). The low-energy electronic states at any Dirac point are described by the massless relativistic Dirac fermions, which lead to a number of unconventional and fascinating phenomena, such as anomalous quantum Hall effect [38]. The asymmetric hoppings make the two Dirac points move away from the Kpoints and approach each other [39]. At a critical asymmetry the Dirac points merge at one of the four inequivalent M -points, that is, (±π, 0) and (0, ±π). A band gap opens upon increasing the asymmetry further; see illustration in Fig.2(b). The topological invariant that characterizes the topological properties is given by the Zak phase, which is the integration of Berry connection over in the reduced Brillouin zone [40,41]. If the integral loop encloses a gap closing point, the Zak phase gives where is the eigenvector of lower band and Otherwise it is zero. The merging is a topological transition since the Berry phases ±π Scaling of entanglement entropy in honeycomb lattice on a torus associated with the two Dirac points annihilate each other at the critical point.

Entanglement entropy of noninteracting fermions
The entanglement entropy in the ground state of free fermions can be determined from the correlation matrix as a fictitious Hamiltonian defined as [42,43,44,45] where c i (c † i ) is the annihilation (creation) operator for site i, j ∈ A. There exist an simple relation between correlation matrix of order O(N A ) and the reduced density 46], where N A is the number of sites in A subsystem. Calculating the eigenvalue λ n of matrixĜ requires only the diagonalization of N A × N A matrix, and then the entanglement entropy is given by The eigenvalue λ n (0 ≤ λ n ≤ 1) ofĜ is termed entanglement spectrum, which is an alternative description of Li and Haldane's version [47]. Recently the entanglement spectrum has been proposed as a ground-state property that captures characteristic edge excitations [48,49,50]. In a large subsystem, most of the eigenvalues lie exponentially close to 0 and to 1 [51]. A special role is played by eigenvalues λ n = 1/2 corresponding to a zero-energy mode [52], and it will contribute a maximal value 1 to the entanglement entropy by Eq. (10).
In the following, we consider a torus of size L x × L y , and a subregion A can be designated as Ω= [0, z 1 L x ]×[0, z 2 L y ] with z i ≤ 1. Here we measure the entanglement between a cornerless cylindrical region A and its complement by two virtual cuts along x or y direction; see illustration in Fig.1(a). Zigzag edges emerge following a horizontal bipartition, i.e., z 1 = 1, z 2 = 1/2, and k x labels the momentum running along the zigzag ribbon. A vertical truncation produces armchair edges, i.e., z 1 = 1, z 2 = 1/2, and k y is well defined instead. One can simultaneously classify the entanglement spectrum with respect to the momentum parallel to the cut, and an intimate relation between entanglement spectrum and eigenspectrum of subblock A can be discovered.
In Fig.3 (a), a flat band of zero energy exists for momenta connecting with two Dirac points along x direction near the zigzag edges. Such zero-energy edge states are reflected by the maximally entangled modes in entanglement spectrum, as shown in Fig.3 (b). The fermion doubling in a chiral symmetric system is a universal property of the bulk, and the bulk-edge correspondence guarantees the edge states appear as pairs. In other words, such doubly degenerate entanglement modes are protected by the chiral symmetry. In the semimetallic phase with zigzag edges, the contributions to the entanglement entropy can be broadly separated into those from the edges and those from bulk, while only bulk states contribute to the entanglement entropy in the insulating phase (see Fig.4). However, the situation is quite different for armchair edges. The interedge interaction of armchair edges in the semimetallic phase makes the localized states fade away [53], and thus edge states play a trivial role in the entanglement entropy (see Fig. 5). Nevertheless, as is shown in Fig.6, the anisotropy will induce number of edge states, which contribute significantly to the entanglement entropy. It is explicit that the zero-energy states have one-to-one correspondences with the maximally entangled modes in the entanglement spectrum. In a sense, the entanglement spectrum is a very promising tool to characterize chiral edge states. We show the half-space entanglement entropy as a function of t c along the normalized path t a + t b + t c = 1 for t a = 0.15, µ = 0 in Fig.7. The entanglement entropy changes discontinuously when t c moves across points t 1 c = |t a − t b | and t 2 c = t a + t b . In fact, the abrupt change associated with the topological phase transition can be grasped more clearly by its first derivative for both types of bipartition and any µ within the bulk bands. The system is an insulator when t c < |t a − t b | or t c >t a + t b , and enters the semimetallic phase when |t a − t b |< t c <t a + t b . In the insulating phase, the   entanglement entropy is linear with the block length for both edges, and the deviation seems negligible. However, the finite-size effect is prominent in the gapless phase. There exists a remarkable correction to the linear dependence on the block length. We note that the correction is likely a logarithmic form for armchair edges, whereas an oscillating correction exists for zigzag edges. We attribute this peculiar correction to finite-size effect of edge modes. The trends in the proliferation of edge modes for both edges are converse as we approach a quantum critical point, as displayed in the inset of Fig.7. We find that the existence of edges state is that t a , t b and t c should satisfy the triangle inequality for zigzag edge [54]. Next we investigate the scaling behavior of bipartite entanglement entropy in the presence of different edges. The scaling of the entanglement entropy has attracted much attention by a flurry of recent works [55,56,57,58]. A general form of the entanglement entropy in two-dimensional system is argued to be of the form where c i (i=1,· · ·, 4) are size independent coefficients. Eq.(11) includes a universal logarithmically divergent correction and a nonuniversal area-law term. Which term is the leading term requires further elucidation of the Fermi surface. We will give both analytical and numerical results for the logarithmical scaling of entanglement entropy in one concrete case. Apparently in the insulating phase, the entanglement entropy grows asymptotically as L [59], as supported in Fig.8(a). For µ = 0 in the semimetallic phase, the Fermi surface reduces to zero dimension, and then the entanglement converges to the area law [55,60,61], as is shown in Fig.9(a). We see that the slope of scaling changes for varying parameters in the semimetallic phase, so the area-law term is nonuniversal. Apart from the leading power-law term, the subleading term is of general interest. To this end, we plot the quantity which will cancel the area-law term and hence extract the subleading contribution. It is remarkable that S sub vN behave diversely in different phases. If the entanglement entropy follows a behavior as S vN (L) = c 2 L + c 4 + O(1/L), S sub vN will approach a constant term like as is shown in the inset of Fig.8(a). In contrast, S vN displays a divergent oscillation in the inset of Fig.9 (a). A careful analysis identifies a logarithmic additive term in these critical phases by following S vN (L) = c 2 L + c 3 log 2 L + O(1), and then the subarea law of the entanglement entropy yields [61,62] S sub vN ∼ −c 3 log 2 L.
The subdominant contribution stems from one chiral mode in the zero-dimensional Fermi surface, which infers c 3 = 1/3 by the FisherCHartwig conjecture [63,64]. The oscillation can be interpreted as the following: In case of symmetric hoppings, the existence of zeroenergy localized states is restricted for 1/3 of the total momentum parallel to the zigzag edge, and hence a multiple of 3 in block length favors more edge states. Hence, the oscillation has a period of 3, as is observed in the inset of Fig.7. A significant difference of the scaling in the entanglement entropy is found in the presence of finite density of state at the Fermi surface. In such situation, every lowenergy radial mode is approximately treated as a chiral excitation that contribute a lnL to the total entanglement entropy [65]. Therefore, a multiplicative logarithmic correction to the power-law behavior arises for finite Fermi surface (i.e., line nodes), which are reported in Fig.8(b) and Fig.9(b) [66,67]. Noticeably, we see that the entanglement entropy with armchair edge is larger than that with zigzag edge, and the reason will be accounted soon.
Widom conjecture [68] predicted that the prefactor of the logarithmic term is determined by chemical potential µ via an accurate formula [69] where Ω is the real-space region of A, rescaled by L such that Vol(Ω)=1. Vectors n x and n k denote the normal vectors on the real-space surface ∂Ω and the Fermi surface ∂Γ(µ). From Eq.(15), we see the topology of the subsystem and the Fermi surface has intimately related to the entanglement entropy. Thus the existence of edge states crucially affects the physical properties of graphene system. Besides, it was argued that Lifshitz quantum phase transitions are accompanied with a nonanalyticity in the prefactor of the leading-order term [70]. In two dimensions the calculation of c 1 (µ) can be simplified. The normal unit vectors in real space are ±ŷ for zigzag-type bipartition and  ±x for armchair-type bipartition. Since the Fermi surface is a curve in two dimensions, we can parameterize the curve by k x = k x (θ) and k y = k y (θ) like [34] (t a + t b ) cos k x + t c cos k y ≡ µ cos θ, Then we have Eq.(16) and Eq. (17) indicates c 1 (µ) are proportional to the variation of Fermi surface along x and y, respectively. For |µ| is large, the Fermi surface is an oval, while the disconnected Fermi components appear with the decrease of |µ|. Without loss of generality, we give a simple demonstration of the solution for symmetric case t a = t b = t c = 1. The dispersion has particle-hole symmetry and 90-degree rotational invariance, so we only consider positive µ in the first quadrant of the Brillouin zone, and then the Fermi surface can be parameterized by for 0 ≤ µ ≤ 1 and θ ∈ (−π, 0). Here R(µ, θ) = (− 1 − µ 2 sin 2 θ + µ cos θ)/2. The longitudinal variation ∆k y is 2 arcsin(µ) and the horizontal variation ∆k x is arccos( −µ−1 2 ) − arccos( µ−1 2 ). When µ > √ 5, it is replaced by for θ ∈ [− arccos( µ 2 +3 4µ ), 0]. The longitudinal variation ∆k y is arccos µ 2 −5 4 and the horizontal variation ∆k x is arccos( µ−1 2 ). Otherwise when 1 < µ < √ 5, for θ ∈ [− sin −1 1 µ , 0] and k x (θ) = π − cos −1 R(µ, θ) k y (θ) = π + sin −1 (µ sin θ) for θ ∈ [arccos( µ 2 +3 4µ ) − π, arcsin( 1 µ ) − π]. The longitudinal variation ∆k y is arccos µ 2 −5 4 and the horizontal variation ∆k x is arccos( µ−1 2 ). The resulting prefactors of leading logarithmic terms are: A generalization to anisotropic case and arbitrary edges is straightforward. The coefficient of the leading term in system of open boundary condition is half of the one of periodic boundary condition, due to the fact that the boundary between the bipartite spaces is half of that of stripes [71]. From Fig. 10, we can recognize the entanglement entropy presents a shoulder-like behavior with respect to the chemical potential µ. The finite-size effect becomes less dominated for large size. c 1 (µ) extracted from fits to our numerical data agrees rather well with the exact form Eq. (18) in the inset of Fig.10. c 1 (µ) goes to zero with µ approaches Dirac point, where the logarithmic correction is vanishing. When µ deviates slightly from the critical point, two small Fermi pockets are formed with Fermi vector k F = |µ|/v F , and Eq.(18) implies that c 1 (µ) scales linearly for the low-energy modes when system size L ≫ l * , where l * is a characteristic length l * ∼ k −1 F . The ratios of c 1 (µ) to c 2 (µ) with both edges are above a factor √ 3, because the Brillouin zone of brick-type lattice along y direction is stretched √ 3 times longer than that along x direction comparing with pristine honeycomb lattice. The ratios reach their maxima passing through a Van Hove singularity, i.e., µ = 1, where two Fermi surfaces merge to form a single Fermi surface; see illustration in Fig.10(c). Thus at µ = 1 there is a Lifshitz transition. Besides a logarithmic divergence of the entanglement entropy, the linear-term coefficient is intricately related to ultraviolet cutoff or correlation length, and it is generally nonuniversal. Shown in the lower insets of Fig.10 (a-b) are c 2 (µ), which decrease with the increase of |µ|. Extraordinarily, there is a dip at µ = 1 by a meticulous observation. We attribute such decrease to the enhanced contribution from a negative subarea term at µ = 1 [61].

Spin-orbit interaction
Next, we extend our scope to Kane-Mele model, which is a tight-binding model with second neighbor spin-orbital interaction defined on a honeycomb lattice, Hereτ i = −τ i is the antiparallel vector of τ i in Eq.(4). The Hamiltonian is time reversal invariant and has a wide range of applicability, such as quantum spin-Hall (QSH) insulator [72,73], graphene/graphane interface [74] and silicene [75,76,77,78]. The energy band structure of the brick lattice is given by where Φ k = −t(e iky + e ikx + e −ikx ) and Θ k = 4t 2 cos k y sin k x . Without the secondneighbor hopping, the model's bulk and surface spectra display the non-generic feature that all energy levels come in pairs at each momentum due to parity symmetry [52]. The presence of the spin-orbit interaction, i.e., t 2 = 0, leads to a gap in the bulk spectrum, as shown in Fig.11. More interestingly, in the case of zigzag ribbon, the flat dispersion is replaced by the helical modes intrinsic to the QSH insulator as the spinorbital is introduced [79]. The dominant contributions originates from the helical edge state is displayed in Fig.11(d). The area-law monotonic scaling at µ = 0 is confirmed in Fig.12(a), and the leading logarithmic divergence is also observed in Fig.12(b).

Conclusion
In this paper, we study the von Neumann entropy in free fermionic systems on a twodimensional honeycomb lattice. Distinct bipartitions on the honeycomb lattice will produce different edges, and the edge states will have an effective impact upon the arealaw terms in the entanglement entropy. We exhibit that the zero-energy edge state has an one-to-one correspondence with the maximal entangled mode in the entanglement spectrum, and hence we speculate the main contributions to the entanglement entropy come from the bulk states and edge states. The edge states play an essential role in semimetallic state with zigzag edges and in insulating phase with armchair edges. Furthermore, we find the entanglement entropy obeys area law when the Fermi surface is within the gap or a nodal point. The corrections to the pow-law term behave differently in critical phase from in noncritical phase. A logarithmic violation to the area law emerges in the presence of line-nodes Fermi surface. We also show that the logarithmic scaling highly depends on the topology of the Fermi surface. The prefactor of logarithmic term is determined by the variation of Fermi surface along the direction parallel the edge. The scaling laws are independent of the statistics of the microscopical constituents. The introduction of spin-orbit coupling will cause a gap in the bulk spectrum and a decrease of entanglement entropy.
Since there is an equivalence between the entanglement entropy and particle number fluctuations in free-fermion systems [56,57,58], it is very promising to measure the entanglement entropy in the graphene-like systems. Especially experimentally the edges can be better fabricated in the artificial honeycomb than those of natural electronic graphene, which tend to be very irregular and contaminated with adsorbates. The flat edge modes can be observed in the zero-bias conductance measurements [80]. Our results can be extended to other free fermionic systems, i.e., band topological insulators and superconductors, as well as weakly interacting systems [81].