Entanglement Entropy of Fermi Liquids via Multi-dimensional Bosonization

The logarithmic violations of the area law, i.e. an"area law"with logarithmic correction of the form $S \sim L^{d-1} \log L$, for entanglement entropy are found in both 1D gapless system and for high dimensional free fermions. The purpose of this work is to show that both violations are of the same origin, and in the presence of Fermi liquid interactions such behavior persists for 2D fermion systems. In this paper we first consider the entanglement entropy of a toy model, namely a set of decoupled 1D chains of free spinless fermions, to relate both violations in an intuitive way. We then use multi-dimensional bosonization to re-derive the formula by Gioev and Klich [Phys. Rev. Lett. 96, 100503 (2006)] for free fermions through a low-energy effective Hamiltonian, and explicitly show the logarithmic corrections to the area law in both cases share the same origin: the discontinuity at the Fermi surface (points). In the presence of Fermi liquid (forward scattering) interactions, the bosonized theory remains quadratic in terms of the original local degrees of freedom, and after regularizing the theory with a mass term we are able to calculate the entanglement entropy perturbatively up to second order in powers of the coupling parameter for a special geometry via the replica trick. We show that these interactions do not change the leading scaling behavior for the entanglement entropy of a Fermi liquid. At higher orders, we argue that this should remain true through a scaling analysis.


I. INTRODUCTION
The study of entanglement, which is one of the most fundamental aspects of quantum mechanics, has lead to and is still leading to much important progress and applications in different fields of modern physics such as quantum information [1], condensed matter physics [2][3][4], etc.. To name a few, it has lead to better undertanding of density matrix renormalization group (DMRG) [5][6][7]; it has also been proposed to be a tool for the characterization of certain topological phases [8][9][10].
Among various ways of quantifying entanglement, in condensed matter or many-body physics efforts have mainly focused on the bipartite block entanglement entropy (von Neumann entropy) and its generalizations (Rényi or Tsallis entropy). It has become increasingly useful in characterizing phases [11] and phase transitions [12,13]. The area law [14] is one of the most important results on entanglement entropy: it states that the entanglement entropy is proportional to the area of the surface separating two subsystems. However, thus far there are two important classes of systems that violate the area law: in gapless one dimensional (1D) systems, a logarithmic divergence [13,15] is found where according to the area law the entanglement entropy should saturate as the size of the subsystem grows; in higher dimensions, for free fermions the area law is found to be corrected by a similar logarithmic factor log L [16][17][18][19][20][21][22][23], where L is the linear dimension of the subsystem.
In this work, we first show that the scaling behavior of the entanglement entropy for systems with a Fermi surface is the same as that of 1D systems with Fermi points [24][25][26][27]. We then seek for a generalization of the latter to interacting fermions in the Fermi liquid phase. We first develop an intuitive understanding via a toy model, showing that in this model the entanglement entropy has the same form as that given in by Gioev and Klich (GK) in Ref. [16]. We then develop a more general and formal treatment using the method of highdimensional bosonization [28][29][30][31][32]. This approach will not only lead to a reproduction of the result for free fermions obtained by GK based on Widom's conjecture [33,34], but will also lend itself to the inclusion and subsequent treatment of Fermi liquid type (forward scattering) interactions.
This paper is organized as follows. In Sec. II, we describe the toy model for which the entanglement entropy can be written in the same form as the GK result. Then in Sec. III, we briefly introduce the tool box of multidimensional bosonization, and apply it to free fermions to reproduce the GK formula. The main results of this work are presented in Sec. IV in which we calculate the entanglement entropy of a Fermi liquid for a special geometry using a combination of multi-dimensional bosonization and the replica trick. We subsequently summarize and discuss our results. Some technical details are discussed in two appendices.

Convex
Concave Fermi surface

II. THE INTUITIVE PICTURE -A TOY MODEL
Consider a set of decoupled parallel 1D chains of noninteracting spinless fermions with spacing a as shown in Fig. 1(a). Here we only consider d = 2 for simplicity, but this toy model is viable in general d dimensions. The asymptotic behavior of entanglement entropy in large L limit of a convex subsystem A of this model can be obtained by simply counting the number of chains that intersect A, and each segment contributes a (1/3) log L where L is the linear dimension of the subsystem [14,35]. Due to the logarithm, different shapes only lead to differences at the area law level. Since each segment must have two intersections, we can count the intersections instead, which also automatically takes care of non-convex geometries. Although there is an additional correction for multiple intervals on a single chain [37][38][39], as long as only the log L behavior is concerned, that contribution is negligible. For L large enough, we can write the number of these intersections as an integral over the surface of A projected onto the direction perpendicular to the chains times one half of the chain density, 1/a. To make contact with the GK result, we note this model also has Fermi surfaces as shown in Fig. 1(b) with a total "area" of 4π/a. This enables us to replace the density of chains by an integral over the Fermi surfaces of the system where Γ indicates the occupied area in momentum space so its boundary ∂Γ is the Fermi surface(s). Therefore we can write the entanglement entropy as wheren is the direction along the chains which is also normal to the Fermi surface, and an overall factor of 1 2 accounts for the double counting of chain segments. In Eq. (1) we recover the GK formula in this special case but written in a slightly different way. In Ref. [16] the entanglement entropy is given as: where the real space surface integral is carried out over the subsystem whose volume is normalized to 1. The surface area is factored out as L d−1 . However, in our formula the surface area ∼ L d−1 is implicitly included in the integral over the surface of the subsystem.
We note that the model discussed in Ref. [26] is equivalent with our toy model, but motivated from a different perspective. In Ref. [26], models are constructed from the momentum space, either with Fermi surfaces as our toy model, or a square Fermi surface, and a boxlike and a spherical geometry are discussed. In contrast, our toy model is constructed from a real space perspective, and general single connected geometries are discussed.
Motivated by the toy model, in this work we extend this intuitive understanding of GK's result to generic free Fermi systems and generalize it to include Fermi liquid interactions in two dimensions (2D) via high dimensional bosonization. Using the method of multi-dimensional bosonization, the Fermi liquid theory can be written as a tensor product of low-energy effective theories of quasi-1D systems similar to this toy model, along all directions. This provides us with a tool to treat the entanglement entropy of fermions in high dimensions, even in the presence of interactions.
At this point, we could also include forward scattering for each chain, and from 1D bosonization we know that for spinless fermions this only leads to renormalization of the Fermi velocity, thus does not change the logarithmic scaling of the entanglement entropy for this toy model. This hints that the same conclusion might hold for Fermi liquids, as we can include Fermi liquid interactions in a similar way via high dimensional bosonization. Although as we show later, this is indeed true at the leading order, the situation is more delicate than it seems to be. The Fermi liquid interactions couple a family of "toy models" aligned along different directions in the language of high dimensional bosonization, and lead to a correction to the entanglement entropy ∼ O(1) × log L.

III. MULTI-DIMENSIONAL BOSONIZATION
The scheme of multi-dimensional bosonization was first introduced by Haldane [28], followed by others [29][30][31][32]. The basic idea is to start with a low energy effective Hamiltonian (obtained through a renormalization group (RG) approach) restricted to within a thin shell of thickness λ around the Fermi surface, k F − λ/2 < |k| < k F +λ/2. Then one divides this thin shell into N patches with dimensionality ∼ Λ d−1 × λ as shown in Fig.(2) in such a way that λ Λ k F and Λ 2 /k F λ, where d = 2, 3 is the space dimension, Λ is the linear dimension of the tangential extent of each patch. The condition λ Λ minimizes inter-patch scattering; Λ k F and Λ 2 /k F λ together makes the curvature of the Fermi surface negligible. In the end we shall take the limit Λ/k F → 0, so that the sum over all patches can be converted to an integral over the Fermi surface. In this work, we treat the free theory in general d dimensions, but shall restrict ourselves to d = 2 when interactions are included. For an arbitrary patch S, labeled by the Fermi momentum k S at the center of the patch, we introduce the patch fermion field operator where ψ p is the usual fermion field in momentum space, θ(S; p) = 1 if p lies in the patch S, 0 if p lies outside patch S.
The effective Fermi liquid Hamiltonian can be written as with m * being the effective mass, V (S, T ; x−y) the effective interaction in the forward scattering channels. Even though this model is restricted to special interactions of this form, forward scattering is known to be the only marginal interaction in RG analysis [40]. As the leading order contribution of the entanglement entropy is dominated by the low energy modes around the Fermi surface, it is sufficient to consider this model. Similar to the 1D case, the bosonic degrees of freedom are the density modes of the system, in this case defined within each patch of the Fermi surface: (5) Though q is not explicitly bounded in the above definition of the patch density operator, its transverse components q S⊥ = (q where is the occupation number of state with momentum k,n S is the outward normal direction of patch S, q S⊥ represents all other component(s) of q that are perpendicular ton S , and L 0 is the linear dimension of the entire system. The appearance of δ q+p,0 is a result of momentum conservation. The calculation of the commutator is reduced to computing the difference of occupied states, i.e. the area difference below the Fermi surface, between the two θ functions (θ(S; k − q) − θ(S; k − p)) as indicated by Eq. (6). This is similar to 1D bosonization. If we consider both k and q to be 1D momenta, Eq. (6) would give us the 1D bosonization commutator. The 2D result Eq. (7) is similar, because the Fermi surface confined within the patch is essentially flat thus the dispersion is 1D. That leads to then S · q dependence of the commutator as that of the 1D case, even for q S⊥ = 0. The difference is that, as illustrated in Fig. (3), due to the patch confinement on the transverse direction(s), when q S⊥ increases k ± q would increasingly find itself outside the patch thus not contributing to the commutator. According to Fig. (3), one can see that this gives rise to the factor θ 2 (q S⊥ ), which diminishes the commutator at large q S⊥ . It is usually neglected in literature because the long wavelength limit is taken [29][30][31][32]. However, as it is important in the present context to correctly count the number of total degrees of freedom, this θ 2 (q S⊥ ) factor cannot be neglected because it comes from counting the transverse degrees of freedom. To simplify things, we replace θ 2 (q S⊥ ) by and we also limit q S⊥ to this range. This approximation makes it easier to do Fourier transform while keeping the total degrees of freedom intact. To see that, it is sufficient to consider one direction, comparing the area enclosed by the two different functions: θ 2 (q ⊥ ) = 1 − q ⊥ /Λ over the range (−Λ, Λ) and θ 2 (q ⊥ ) = 1 over the range (−Λ/2, Λ/2). Both functions enclose the same area thus the same number of states. This approximation can also be interpreted as relaxation of the hard wall cutoff in Eq. (5), softening of the step function θ(S; k). In Eq. (5), q is not bounded while k is bounded by θ(S; k). If we relax the restriction on k on the transverse direction, allowing is denoted by blue, that of θ(S; k − q) is denoted by red, and the overlapping region is denoted by yellow. Subtracting the remaining blue area from the red, we obtain that θ(S; k − q) occupies (Λ − q S⊥ )qn S more states, which gives us the commutator.
k with |k (α) S⊥ | > Λ/2 in the summation, but require q α S⊥ to be bounded within the patch, we would obtain the alternative θ 2 (q ⊥ ). Using from now on the above approximation, we construct the local bosonic degrees of freedom φ(S; x) = φ(S; x S , x S⊥ ) as where J(S; x) = q e iq·x J(S; q), x S = x ·n S , and x S⊥ = x − (x ·n S )n S . The commutation relations for the φ's are then which is the bosonic commutation relation we are looking for. The factor α arising from transverse directions must be treated with care in different circumstances. In most literature, the focus is the physics at large length scale l 1/Λ; therefore, this factor is usually approximated by δ d−1 (x S⊥ − y S⊥ ) which is good in that limit without further discussion. This is also what we shall do for most of the time unless noted otherwise: However, the more accurate expression (11) is useful for us to understand how to count the transverse degrees of freedom correctly. It tells us that the transverse degrees are not independent on the short length scale l < 1/Λ. More importantly, later on we need to consider the limit δ d−1 (x S⊥ − y S⊥ )| y S⊥ →x S⊥ ; without Eq. (11), this limit would be ill-defined.
With the above, the Hamiltonian H[ψ † , ψ] is found to be quadratic in terms of these J(S; q)'s: where . So it is also quadratic in the bosonic fields associated with the J(S; q)'s.

A. Entanglement Entropy of Free Fermions
The kinetic energy part of Eq. (4) or its bosonized version Eq. (13) can be written in terms of the boson fields constructed above as: We see that there is no coupling between different patches. The theory is thus formally a tensor product of many independent theories, one for each patch. We can therefore calculate the entanglement entropy patch by patch and sum up contributions from each patch in the end. Within a single patch there is no dynamics in the perpendicular direction as dictated by the Hamiltonian, and the problem is reduced to a one dimensional problem! Note that transverse degrees of freedom are not completely independent. According to Eq. (11), the commutator is non-vanishing for x S⊥ = y S⊥ up to a length scale ∼ 2π/Λ. This is a consequence of restricting q S⊥ to within the range [−Λ/2, Λ/2]. Physically one can view this as discretization along the transverse direction due to a restricted momentum range, similar to the relation between a lattice and its Brillouin zone. In this view the single patch problem is reduced to a 1D problem with a chain density of (Λ/(2π)) d−1 . Therefore, the Hamiltonian (14) becomes Note that the bosonized theory of a single patch is chiral. To directly make use of our toy model, we need to consider two patches having oppositen S simultaneously. This is because for a 1D fermion model at non-zero filling, there are two Fermi points. Both need to be considered to construct well-defined local degrees of freedom. Once we consider such two patches together, it is more convenient to combine the two chiral theories into a non-chiral theories. This is also what we will do for the rest of this work. Introduce the non-chiral fields where −S indicates the patch with normal direction opposite to that of patch S:n −S = −n S . One finds that χ and ϕ are mutually dual fields with S restricted to one hemisphere, but ∂ x S ϕ and ϕ now commute while χ and ϕ have a non-trivial commutator: Therefore, two patches with oppositen S are equivalent to a set of ordinary 1D boson fields. Throughout the rest of this work, we shall assume this chiral-to-nonchiral transformation is done, and when we refer to patches we always refer to the two companion patches that form a non-chiral patch together. For the non-chiral boson theory, it is known that the entanglement entropy of a single interval (with two end points) is (1/3) log L.
Before we proceed further, we note that the relation between boson fields and the original fermion fields is not completely local. However, the underlying physical quantity that matters is not the fields, but the fermion density, or in other words, the fermion number basis one chooses to expand the Hilbert space of the problem. This physical basis is also what one uses to do the partial trace. It is known that the fermion density operator obeys a locally one-to-one corresponding relation to the boson fields. Thus we argue that in 1D the nonlocal relation between the fermion and boson fields does not affect the partial trace operation, so as the calculation of entanglement entropy.
By referring to our result for the toy model, the contribution from a single patch is readily given where an additional factor of 1/2 has been introduced in order to count only once each pair of patches forming a non-chiral theory. Identifyingn S Λ d−1 as the surface element at the Fermi surface d S k and taking the N → ∞ limit, the total entanglement entropy is So we recover the GK result for generic free fermions.

B. Solution for the Fermi liquid case and non-locality of the Bogoliubov fields
When Fermi liquid interactions (forward scattering) are included, the full Hamiltonian will no longer be diagonal in the patch index S. But it is still quadratic in terms of the patch density operators, i.e. the bosonic degrees of freedom, and can be diagonalized by a Bogoliubov transformation. According to Eq. (7) and ignoring terms of O(λ/Λ), one can define a set of boson creation/annihilation operatorsâ † (q)/â(q) as follows: It can be shown that the full Hamiltonian is diagonal in q, and it can be diagonalized by a Bogoliubov transformation [32] independently for each q sector. In Ref. [32], only a Hubbard-U like interaction is considered for practical reasons. But in principle, such a Bogoliubov transformation also applies to general interactions: where both i and j refer to the patch index, α j and β j are the Bogoliubov bosonic annihilation operators that diagonalize the Hamiltonian. With proper choice of u's and v's, the Hamiltonian is readily diagonalized. Ref. [32] solves the Hubbard-U like interaction and provides a successful description of Fermi liquids, even in the strong U limit. However, even for this simple case in which U has no qdependence, the Bogoliubov transformation still depends on q. To be more precise, u ij and v ij will depend only on the angle between the patch normal directionn S and q, leading to discontinuities in the derivatives at q = 0. Consequently, the real space fields constructed from the Bogoliubov operators α j and β j are no longer local with respect to the original boson fields. The real space Bogoliubov fields are constructed in a manner similar to Eq. (20): (22) Then one can show that the original local degrees of freedom φ(S; x) can be expressed in terms of above Bogoliubov fields as where f (S, l; x − y) is typically long-range, even for the short-range Hubbard-U interaction. For more general cases, with further q-dependence in the interaction, the non-locality would only be enhanced. The loss of locality prevents us from calculating the entanglement entropy directly using those eigen modes, since it is difficult to implement the partial trace using those non-local degrees of freedoms. Therefore, although the Bogoliubov fields have a local core as we would expect for Fermi liquids from adiabaticity, they do acquire a nonlocal dressing due to interaction. Though in principle the partial trace can be done with those Bogoliubov fields, such nonlocality makes it difficult and we have not been able to do it, which further renders calculating the entanglement entropy impossible. This is very different from the 1D theory, where for local interactions the eigen fields remain local, since there are only two Fermi points. There the transformation can never involve such angular q-dependence due to limited dimensionality. Despite these technical difficulties, the non-locality may suggest possible corrections to the entanglement entropy. This is indeed the case as revealed by our later calculation for Fermi liquid interactions, although in this case such extra contributions are only of O(1)×log L which is of O(1/L) comparing to the leading term. This shows that the mode-counting argument in Ref. [26], though correctly suggesting the log L violation to the area law for Fermi liquids, does not always fully account for all sources of entanglement entropy.

IV. ENTANGLEMENT ENTROPY FROM THE GREEN'S FUNCTION
In order to preserve locality, we need to work with the original local degrees of freedom. To do that, we adopt the approach used by Calabrese and Cardy [41] (CC) on calculating the entanglement entropy of a free massive 1D bosonic field theory. The calculation is done in terms of the Green's function by applying the replica trick. In our case, we find that the CC approach can be generalized in a special geometry for solving the interacting theory which is quadratic after bosonization. In this way, we avoid diagonalizing the Hamiltonian and thus the nonlocality issue. However, we do have to regularize the theory by adding a mass term by hand. In the end we shall take the small mass limit, and replace the divergent correlation length ξ ∼ 1/m by the subsystem size L. The regularization procedure facilitates the calculation, but also strictly restricts us to computing the entanglement entropy only at the log L level.
In this section, by using the replica trick we convey the calculation of entanglement entropy into computing the Green's function on an n-sheeted replica manifold. We first demonstrate the method by applying it to free fermion theory in d-dimensional space; then based on it we compute the entanglement entropy perturbatively for a simple Fermi liquid theory in powers of the interaction strength up to the second order.

A. The Replica Trick and Application to 1D Free Bosonic Theory
In this part, we briefly describe the replica trick in (1 + 1) space-time dimensions ((1 + 1)d) so that later on we can straightforwardly generalize it to (2 + 1) spacetime dimensions ((2 + 1)d) accordingly for our problem.
The replica trick makes use of the following identity: To compute tr ρ n A , CC use path integral to express the density matrix ρ in terms of the boson fields where Z = tr e −βH is the partition function, β is the inverse temperature, and {φ(x)} are the corresponding eigenstates ofφ(x):φ(x)|{φ(x )} = φ(x )|{φ(x )} . ρ can be expressed as a (Euclidean) path integral: where S E = β 0 L E dτ , with L E being the Euclidean Lagrangian. The normalization factor Z, i.e. the partition function is found by setting {φ(x) } = {φ(x) } and integrating over these variables. This has the effect of sewing together the edges along τ = 0 and τ = β to form a cylinder of circumference β as illustrated in Fig. (4) The reduced density matrix of an interval A = (x i , x f ) can be obtained by sewing together only those points which are not in the interval A. This has the effect of leaving an open cut along the line τ = 0 which is shown in Fig. 4 (right panel). To compute ρ n A , we make n copies of above set-up labeled by an integer k with 1 ≤ k ≤ n, and sew them together cyclically along the open cut so Fig. 5(a) we show the case n = 2. Let us denote the path integral on this n-sheeted structure (known as n-sheeted Riemann surface) by Z n (A). Then so that If we consider the theory as that of one field living on this complex n−sheeted Riemann surface instead of a theory of n copies, it is possible to remove the replica index n from the fields, and instead consider a problem defined on such an n-sheeted Riemann surface which can be realized by imposing proper boundary conditions. In Ref. [41], CC consider the entanglement entropy be-tween the two semi-infinite 1D system (i.e. cutting an infinite chain into two halves at x = 0) for free massive boson fields. For such geometry, as illustrated in Fig. 5(b), the n−sheeted Riemann surface constraint is realized by imposing a 2nπ periodicity on the angular variable of the polar coordinates of the (1 + 1)d plane instead of the usual 2π one. In this way, the (1 + 1)d variable x = (x, τ ) acquires n branches x n , and each branch corresponds to one copy of φ. Notation-wise this corresponds to and the sewing conditions φ(x) k = φ(x) k+1 simply becomes the continuity condition for φ(x) across its consecutive branches. Here we use a generalized polar coordinate: x = (r, θ) with 0 < r < ∞, and 0 ≤ θ < 2nπ. The massive free boson theory considered by CC is defined by the following action The (1 + 1)d bosonic Green's function G (n) 0,b (r, r ) = φ(r)φ(r ) on the n−sheeted Riemann surface satisfies the differential equation To compute the partition function, one can make use of the identity Note that here the integration is over the entire n−sheeted space. The above is applicable to general quadratic theories of bosons, and will be applied by us later to bosonized theories of interacting fermions. Here we use G (n) (x, x ), a general two point correlation function on the n−sheeted Riemann surface in d-dimensional space for later use, instead of the specific G (n) 0,b defined above. Accordingly, S A is then given as Here and in the following, we will leave it understood that the first term in the integrand is integrated over the nsheeted geometry, whereas the second is integrated over a one-sheeted geometry. There should be no confusion as the superscript of G generally indicates the geometry.
The benefit of the above approach is that the two point correlation function or Green's function, defined in terms of certain differential equation obtained from the equation of motion, can be solved for on the n−sheeted Riemann surface thus enabling us to compute the entanglement entropy. Although CC's work only considers massive (1 + 1)d boson fields, it is also applicable to our case. The price one has to pay is to to introduce a mass term for regularization. At the end of the calculation the inverse mass, which is the correlation length of the system, shall be considered to be on the same scale as L: 1/m ∼ L, where L is the characteristic length scale of the subsystem. The validity of such consideration is well-established in other cases, [35,42] where the correlation length is either set by finite temperature or mass. The only modification necessary to apply the above to a bosonized Fermi surface in higher dimensions is to introduce a sum over the patch index.

B. Geometry and Replica Boundary Conditions
Through the remainder of this work, instead of the general geometry considered before, we work with a special half-cylinder geometry as shown in Fig. 6(a): the system is infinite in thex direction while obeying periodic boundary condition along theŷ direction with length L. The system is cut along theŷ axis so that we are computing the entanglement entropy between the two half planes. We require L to be large so that it can be considered ∼ ∞ unless otherwise noted.
We choose such this simple geometry for the following reasons. Cutting the system straight along theŷ direction, yielding a two half-plane geometry, is a straightforward (2 + 1)d generalization of the semi-infinite chain geometry considered in CC. It makes any straight line intersect the boundary only once, dividing it into two semi-infinite segments, for all patch directions as in the 1D case, except for lines parallel to then S =ŷ patch direction. The degrees of freedom associated with this special patch do not contribute to the entanglement entropy, since they are not coupled (have no dynamics) alongx, and are of measure zero in the large patch number limit anyway.
For this simple geometry, the (2 + 1)d n-sheeted geometry is constructed from n identical copies sliced along "branch cuts" and then appropriately glued together along these cuts. This happens exactly as in 1D, and the y coordinate is so far a mere spectator. This defines an n−sheeted, or in this case more appropriately the n−layered, replica manifold which is a simple enough generalization of the (1 + 1)d case. The n-sheeted Riemann surface, as discussed in Sec. IV A and shown in Fig. 5(b), now acquires an extra directionŷ perpendicular to thex −τ plane. It can still be implemented by imposing the same 2nπ periodicity boundary conditions on θ, the angular variable of the polar coordinates (x, τ ) = (r cos θ, r sin θ) in thê x −τ plane. Therefore, we can safely make use of the The half-cylinder geometry and equivalence of boundary conditions inx −τ andn S −τ planes. The system is infinite in thex direction while obeys periodic boundary condition along theŷ direction with length L. The system is cut along theŷ axis so that we are computing the entanglement entropy between the two half planes. (a): The half-cylinder geometry. (b): The projection of r S onto thex −τ plane. Consider polar coordinates of an arbitrarŷ n S −τ plane (the blue plane). Since the polar coordinates in thex −τ plane satisfies the 2nπ periodic boundary condition, consider the one-to-one projection of the vector r S onto thê x −τ plane. Consider, if we move the vector in thex −τ plane around the origin n times (the red circle). Due to the one-to one mapping, r S should also move around the origin n times (the blue "circle", it is actually a eclipse), thus obeys the 2nπ periodicity as well.
CC result, i.e. the solution to the Green's function on a n−sheeted Riemann surface, to the free fermion theory, and can further use it as a starting point for treating the interacting theory. This is obviously true for the patch withn S =x, but it also holds for generaln S as we shall validate as the following.
For a general patch directionn S , the noninteracting Green's function associated with this patch embodies correlations in the affinen S −τ "planes". The geometry of each such "plane" is that of the n-sheeted Riemann surface of the (1 + 1)d problem, as we will now argue. With each patch direction we thus associate a different foliation of the n-sheeted (2 + 1)d geometry into (1 + 1)d counterparts.
To be more precise, for given patch S, instead of the Cartesian coordinates (x, y, τ ), we consider a parallel/perpendicular decomposition (x S , x S⊥ , τ ) for each sheet via where then S⊥ are the perpendicular unit vectors aligned with the patch S. The natural choice of coordinates for a given patch is to choose polar coordinates within thê n S −τ plane: because these are the coordinates in which then S −τ planes restricted to each sheet are naturally glued together by extending the range of θ S to 2πn, as we will now show. The shift x S⊥n x S⊥ /n x S of x S is necessary as to ensure that x = 0, the location of the onset of the branch cut, corresponds to r S = 0 which is what makes these coordinates convenient. Then S −τ planes are now defined by fixed x S⊥ .
If we can establish that the 2πn periodicity of θ is equivalent to a 2πn periodicity of θ S , then CC's solution would be justified in the above set-up so that the noninteracting Green's function G 0 (S, S; , r S , θ S , r S , θ S ) can be expressed through CC's result. This can be achieved by establishing a one-to-one correspondence (mapping) between θ and θ S . The mapping is intuitively constructed, as shown in Fig. 6(b), as the vertical projection from then S −τ plane (the blue plane in Fig. 6(b)) onto thex−τ plane along theŷ direction. Consider moving the projection of r S in thex −τ plane around the origin n times (the red circle). It is clear that r S (on the blue ellipse) follows its projection while also moving around the origin n times, always being on the same sheet. In particular, the branch cut is always traversed simultaneously for θ = θ S = π mod 2π. Then S −τ planes, the leaves of our foliation, thus have the familiar 1+1d n-sheeted geometry, and θ S obeys the same 2nπ periodicity as θ.
Finally, the periodicity condition of theŷ direction is necessary for the total entanglement entropy to be finite; it also provides the only length scale for the subsystem which is needed for extracting the scaling behavior of entanglement entropy. However, if we are only concerned with the integral form of the entanglement entropy as in Eq. (19), not requiring it to be finite as a whole, but rather requiring only the entanglement entropy per unit length to be finite, we may take the y direction to be infinite. This point of view will be taken here and in the following in order to simplify our calculation.

C. Entanglement entropy of free fermions revisited
In order to treat the interacting theory, in this section we re-derive the free fermion result for the half cylinder geometry via the replica trick. Later we shall generalize the method to include interactions. Rewriting the Hamiltonian Eq. (14) in terms of the non-chiral fields, and adding the mass term by hand, we have For convenience, we use the Lagrangian formalism and work with the ϕ(S; x) representation through the rest of this work. Switching to imaginary time t → iτ , and rescaling the coordinates in the following manner: we obtain the following Lagrangian density in the ϕ representation : (38) Then we can work out the Euler-Lagrangian (E-L) equation of motion. Making use of the E-L equation of motion, we find the Green's function G (n) 0 (S, T ; x, x ) = T ϕ(S; x)ϕ(T ; y) 0 satisfies the following differential equation: . The rescaling makes the Green's function dimensionless, thus easier to handle when it comes to computing d d+1 xG (n) (x, x). The extra factor C generated on the right hand side (rhs) will be canceled by the Jacobian of the integral over the Green's function, leaving only a factor of 1/m 2 . All that needs to be computed is then an integral over the dimensionless G. Therefore, it is legitimate to ignore this factor from now on. The δ-functions originate from the commutator Eq. (11), and are coarse-grained. After we include the patch index, perform the integral over m 2 , and take the n derivative, Eq. (31) becomes where a 0 is an ultraviolet cutoff, and The exponential factor in Eq. (31) becomes one after the n → 1 limit is applied. Note that 1 2 log(m 2 a 2 ) ∼ − log L. Our major task is now computing C G (S; n).
Observing that there is no x S⊥ dependence on the left hand side (lhs) of Eq. (39), we can write and we obtain a (1 + 1)d equation (42) in which r S,x(y) is as defined in Eq. (35). The same equation appears in CC. We shall also suppress the subscript S unless necessary, as it is normally already specified in the notation for G 0,b .
The transverse part of the integral in C G (S; n) can be factored out as Recalling our discussion about Eq. (11), this is a coarsegrained δ-function. At short distances, instead of a divergence, we should use Therefore, the transverse direction integral becomes Identifying Λ d−1n S as the surface element dS k , for a given patch the integration can be rewritten as (2π) −d+1 ∂ A |dS x · dS k |. This leaves us with only an integral over (G The solution for the (1 + 1)d Green's function on the n−sheeted replica manifold is given in CC: where d 0 = 1, d k = 2 for k > 0, C ν (θ) = cos(νθ), g ν (r, r ) = θ(r − r )I ν (r )K ν (r) + θ(r − r)I ν (r)K ν (r ), and I ν (r) and K ν (r) are the modified Bessel functions of the first and second kind respectively. r and θ are again the polar coordinates of then S −τ plane, and we have suppress the index S of r as only one patch direction is involved.
The integral over G The integral is divergent since the integrand r x g k/n (r x , r x )| rx→∞ = 1/4, a consequence of the fact that we are calculating the partition function of an infinite system. But this divergence should be canceled in C G (S; n) − nC G (S; 1). To regularize the divergence, we use the Euler-MacLaurin (E-M) summation formula following CC, and sum over k first: where B 2n are the Bernoulli numbers, Note that the first term, the integral over k, is always canceled by rescaling k/n → k in g k/n . For the remaining terms, which contain derivatives with respect to k, we may add a constant, in this case −1/4, under the derivative, which allows us to pull the derivative outside the integral. The integrand now is well-behaved at infinity. To be more precise, according to Eq. (46), we need to compute So we have Combining the above results into Eq. (40) and converting the sum over S into an integral around the Fermi surface, we obtain Eq. (19) for this geometry.

D. Differential Equations of the Green's Functions and an Iterative Solution
In this part, we derive the differential equations of the Green's functions for the quadratic boson theory with inter-patch coupling, and provide an iterative solution.
Including the Fermi liquid interaction V (S, T ; x − y) = U S,T , the Hamiltonian becomes where g S,T = U S,T Ω 2πv * F is order 1/N . This Hamiltonian can be written in terms of the non-chiral fields as where we have made use of the fact that g S,T = g −S,−T , which is required by time-reversal symmetry. Here the summation over S is restricted to a semicircle. This Hamiltonian contains generalized type kinetic terms (inter-patch coupling due to interaction) which are not diagonal. To obtain the corresponding Lagrangian, one needs to invoke the general Legendre transformation [43], and obtains the following Lagrangian densities, respectively, in terms of ϕ or χ: where and h 1(2) (S, T ) is defined through Here, I is the identity matrix, and [f (h) i (S, T )] is the matrix formed by f (h) i (S, T ), i = 1, 2. Applying this result and making use equations of motion obtained from the Hamiltonian, we obtain the Lagrangians L ϕ or L χ .
Here we arbitrarily choose to work with Lϕ. Then, by making use of the E-L equation of motion, applying the same rescaling Eq. (37), and letting t = iτ , we obtain the differential equations that the interacting Green's function G (n) = ϕ(S; x)ϕ(T ; x ) 's satisfies: Here the Jacobian due to change of variables is the same as in the free fermion case. The entanglement entropy is still given by Eq. (40), but replacing G (n) 0 with G (n) in C G (S; n). In the following, we omit the replica index n in the Green's function unless different values of n are involved in a single equation.
As is well known, differential equations such as the above can be converted to an integral form [44] relating the full Green's function to the noninteracting one. This leads to an iterative (perturbative) definition of the former in terms of the latter. In the present case, this integral equation reads G(S, T ; x, y) = G 0 (S, T ; x, y) + δG(S, T ; x, y).
Given this equation, we can now compute the Green's function and thus the entanglement entropy perturbatively in powers of U .

E. Entanglement Entropy from the Iterative Solution
In Eq. (54), the G 0 term is the same as that of the free fermions, thus yields the same contribution to entanglement entropy. To study how the correction term δG(S, T ; x, y) affects the entanglement entropy, we need to study where δ (M ) G denotes the M th order correction. There are two distinctive types of terms in the perturbative expansion of δG. In general, at order M , we have in total 3(M +1) integrals. Let us examine one of the many terms contributing to the M-th order correction, to be summed over patch indices: Here we only include the τ −derivatives. In general we would also have spatial (n li ) derivative terms, as well as terms with mixed derivatives. Butτ andn S directions are equivalent. Using rotational symmetry, and the fact that the two different derivatives in each term are with respect to independent variable that are each integrated over, one can see that all terms are identical except for S-dependent pre-factors. The two categories of terms are defined by the set {l i }: 1) l i = S ∀ i, i.e. with intra-patch coupling only; and 2) ∃ l i = S containing inter-patch coupling. We shall label the two categories as

Intra-patch coupling and comparison with 1D
Setting l i = S for all i's in Eq. (56), first we consider the transverse direction We can immediately integrate out the transverse component of all z i 's and obtain In the last line we use again the fact that the transverse δ-function is a coarse-grained one (Eq. (43)).
The rest of δ (M ) intra G is obtained by substituting G 0 (S, S; z i , z i+1 ) with the (1 + 1)d Green's function G 0,b (S; z S,i , z S,i+1 ). Although a direct computation is possible, we first give a general argument that for any M the contribution to entanglement entropy from δ M intra G vanishes. We do so by making a comparison with the 1D case where a rigorous solution is available.
For the 1D Luttinger liquid with only forward scattering, the entanglement entropy can be calculated directly via bosonization and the result remains at 1/3 log L in the presence of interactions. The calculation is possible because, in our language, there are only two patches, so the transformation which diagonalizes the Hamiltonian is not plagued by the nonlocality issue we encounter in the 2D theory. However, we can also treat the 1D case with our perturbative approach. The resulting series of integrals turns out to be identical to the one obtained from the intra-patch contributions in the higher dimensional case except for the transverse δ−function. Therefore, we argue that at all orders, the intra-patch coupling terms have vanishing contribution to the entanglement entropy. We shall demonstrate such behavior explicitly up to second order in U later on.

Scaling analysis of inter-patch coupling
For terms with inter-patch coupling, we find that they are of order O(1/L) comparing to the leading term according a scaling argument. The crucial observation here is that, as long as ∃ l i = S, we do not encounter the factor where θ S (θ l ) is the angle betweenn S (n l ) and thexaxis in the x − y plane. Therefore, when we integrate out the (M + 1) transverse δ-functions, the factor δ D−1 (0) = L D−1 (Λ) D−1 would be suppressed by even a single l i = S.
To examine the remaining integral, we can ignore the angular part as it cannot affect the scaling behavior. The asymptotic expansion of K ν (r) and I ν (r) for real r at large value is [45] K ν (r) π 2r e −r 1 + ∞ n=1 (ν,n) (2r) n , where (ν, n) = Γ(1/2+ν+n) n!Γ(1/2+ν−n) . By using the above asymptotic expansion of Bessel functions, the leading term for ∂ 2 All of these terms peak around z i+1 = z i and are otherwise exponentially suppressed. We may therefore again estimate this integral by letting x = z 0 = z 1 = · · · = z M −1 and removing (M + 1) of the integrals. The remaining integrals yield, at the leading order, d M +1 z 1/z M ∼ dzz M z −M . However, at the leading order, there is no ν dependence. According to the formalism in Sec. IV C, such terms have no contribution to the entanglement entropy. Therefore, the term that contributes to the entanglement entropy is the next order which behaves as dz 1 z and is of order O(log L), leading only to a correction ∼ O(log L) × log L to the entanglement entropy.
Next, we shall demonstrate in detail our above analysis, for both inter-patch and intra-patch coupling terms by explicit calculation up to the second order.

First Order Correction
The first order term correction to dxG(S, S; x, x) is As we have pointed, it is sufficient to calculate either piece of the two terms due to the equivalence of the imaginary time direction and the real space direction. The other piece should be just the same except for the coefficient. Here we choose to compute the first term. The transverse degrees of freedom provide an overall factor counting the total degrees of freedom as discussed in the general case. Then we can also integrate out the angular degrees of freedom in thex S −τ plane, both θ x and θ z as defined in Eq. (44), after which one obtains where δ (1) G k/n = dr x dr z r x r z g k/n (r x , r z ) The two k summation is reduced to one due to orthogonality of the angular function C k/n (θ). By employing the E-M formula and properties of the Bessel functions, we show in Appendix A that sum over k-values in Eq. (61) can be converted into an integral, which cancels in Eq. (40) for the same scaling reasons discussed above, following Eq. (46). Therefore, we find that the contribution of Eq. (60) to the entanglement entropy vanishes.

Second Order Correction
The second order correction is • for l = S: According to our general discussion, we only need to consider the following piece: where δ (2) G k/n = dr x dr z dr 1 r x r z r 1 g k/n (r x , r z ) × (∂ 2 rz − k 2 /(r z n) 2 )g k/n (r z , r 1 ) × (∂ 2 r1 − k 2 /(r 1 n) 2 )g k/n (r 1 , r x ).
In the above, we have proceeded as in the first order calculation, integrating out the angular part first to obtain the expression for δ (2) G k/n . After a lengthy but similar calculation as for the first order (see Appendix B), we find, using the E-M formula: where p k/n (r x , r z , r 1 ) is the product of g k/n dependent terms in Eq. (66). The usual scaling argument for the integral shows that the entire expression is proportional to n, and thus cancels the second (n = 1) term in Eq. : Therefore, at the second order level for the l = S piece we still have no correction to the scaling law of entanglement entropy.
• for l = S: The integrand we need to consider is The first thing to notice in Eq. (68) is that we have derivatives along directions different from the patch normal directionn S acting on the non-interacting Green's function. We expand this term as For the first term, we can decompose the derivative ∂ z S into terms that act alongn l and along its transverse direction, respectively. The non-interacting Green's function only depends on the transverse coordinates via G 0 (l, l; x, y) ∼ δ(x (l) . For the second term, we integrate by parts with respect to z S , which leads to (including now the first G 0 factor, which depends on z) Therefore, the overall integrand is proportional to 1⊥ ). Note that the x (S) ⊥ dependence only appears in these δ-functions, we can integrate it out, leaving only δ(z

V. SUMMARY AND CONCLUDING REMARKS
In this paper, we developed an intuitive understanding of the logarithmic correction to the area law for the entanglement entropy of free fermions in one and higher dimensions on equal footing -the criticality associated with the Fermi surface (or points). Then we used the tool of high dimensional bosonization to compute the entanglement entropy, and generalized this procedure to include Fermi liquid interactions. In the presence of such interactions we calculated the entanglement entropy for a special geometry perturbatively in powers of the interaction strength up to the second order, and find no correction to the leading scaling behavior. We also point out that the situation is the same at higher orders. Our results thus strongly suggest that the leading scaling behavior of the block entanglement entropy of a Fermi liquid is the same as that of a free Fermi gas with the same Fermi surface, not only for the special block geometry studied in this paper, but for arbitrary geometries. Explicit demonstration of the latter is an obvious direction for future work.
In the special geometry in which we performed explicit calculations using the replica trick, a mass-like term is introduced to regularize the theory at long distance, as is done in closely related contexts [35,42]. For a Fermi liquid (which is quantum-critical) the corresponding length scale ξ ∼ v F /m must be identified with the block size L, and is thus not an independent length scale. On the other hand, such a mass-like term can also describe a superconducting gap due to pairing. In particular, for a weak-coupling superconductor, ξ, the superconducting coherence length, is much longer than all microscopic length scales, but finite nevertheless. In this case it is independent of L, and the interplay between the two is interesting. For L < ξ, the Fermi liquid result (1) still holds. But for L > ξ, the logarithmic factor in the en-tanglement entropy saturates at log ξ, and we expect: S(ρ A ) = 1 12(2π) 2−1 log ξ × ∂A ∂Γ |dS x · dS k |, (76) which agrees with the conjecture made in Ref. [26].
More generally, Fermi liquids are (perhaps the best understood) examples of quantum critical phases (or points) in high dimensions. Unlike in 1D where conformal symmetry powerfully constrains the behavior of entanglement entropy, our understanding of entanglement properties of such high-dimensional quantum critical phases or points (many of them have Fermi surfaces but are not Fermi liquids) is very limited. Our work can be viewed as a step in that general direction. Furthermore, the formalism developed in this work has potential applicability to systems with composite or emergent fermions with Fermi surfaces as well, or more generally, non-Fermi liquid phases with Fermi surfaces. The system studied in Ref. [22], where there is an emergent spinon Fermi surface, is a potential example.
The IK 1 term in h(r, r 1 ) results from δ-functions (derivatives of δ-functions) coming from the derivative applied on the step function (θ(r)). The remaining part comes from terms involving a product of two θ-functions. Here one needs to distinguish between r > r 1 and r < r 1 , which gives rise to the terms in θ(r − r 1 ) and θ(r 1 − r), respectively. The remaining integral drdr 1 rr 1 (θ(r − r 1 )h(r, r 1 ) + θ(r 1 − r)h(r 1 , r)) can be carried out by applying the identities of Bessel functions I and K used in Appx. A. In applying the E-M formula to the sum over k in (65), we again arrive at a divergent k-integral that can be rescaled and subsequently canceled (see (67) and below), and a sum over derivative terms that are well-behaved. In the latter terms, we always add proper constants under the derivatives to regularize the integrand at infinity, as before. We divide (B2) into two terms. The first is the one containing the δ-function. After integrating out r 1 , this term becomes ∂ j k (− drrh(r, r)) = ∂ j k (− dr(rh(r, r) + 1 4 )) = ∂ j k (− dr(r(−IK + r 2 2 (−K 2 I + I − + I 2 K + K − )) + 1 4 )) = ∂ j k (0).