Renyi entanglement entropy of Fermi liquids and non-Fermi liquids: Sachdev-Ye-Kitaev model and dynamical mean field theories

We present a new method for calculating Renyi entanglement entropies for fermionic field-theories originating from microscopic Hamiltonians. The method builds on an operator identity which we discover for the first time. The identity leads to the representation of traces of operator products, and thus Renyi entropies of a subsystem, in terms of fermionic-displacement operators. This allows for a very transparent path-integral formulation, both in and out-of-equilibrium, having a simple boundary condition on the fermionic fields. The method is validated by reproducing well known expressions for entanglement entropy in terms of the correlation matrix for non-interacting fermions. We demonstrate the effectiveness of the method by explicitly formulating the field theory for Renyi entropy in a few zero and higher-dimensional large-$N$ interacting models based on the Sachdev-Ye-Kitaev (SYK) model, and for the Hubbard model within dynamical mean-field theory (DMFT) approximation. We use the formulation to compute Renyi entanglement entropy of interacting Fermi liquid (FL) and non-Fermi liquid (NFL) states in the large-$N$ models and compare successfully with the results obtained via exact diagonalization for finite $N$. We elucidate the connection between entanglement entropy and residual entropy of the NFL ground state in the SYK model and extract sharp signatures of quantum phase transition in the entanglement entropy across an NFL to FL transition. Furthermore, we employ the method to obtain nontrivial system-size scaling of entanglement in an interacting diffusive metal described by a chain of SYK dots.


I. INTRODUCTION
Quantum entanglement has emerged as a major tool to characterize quantum phases and phase transitions [1][2][3][4][5][6] and to distill fundamental quantum mechanical nature of non-trivial many-body states, e.g.those with topological order that is otherwise hard to quantify [7,8].Recent developments in condensed matter and high energy physics have revealed beautiful connections between entanglement, thermalization and dynamics, leading to classification dynamical phases of interacting quantum systems into thermal and many-body localized (MBL) phases [9][10][11].Quantum entanglement is quantified in terms of properties of the reduced density matrix, ρ A = Tr B ρ, of a system with density matrix ρ and divided, e.g., into subsystems A and B, where Tr B is the partial trace over the degrees of freedom of B. Typical entanglement measures constructed out of ρ A for a pure state are von Neumann and Renyi entanglement entropies.The latter can be used to compute useful measures [12][13][14][15][16] for entanglement even in a thermal mixed state.
However, calculation of entanglement entropy is much more challenging than, e.g. that of thermal entropy or usual correlation functions.Over the last decade a lot of progress has been made to obtain entanglement entropy, both numerically and analytically, for non-interacting bosonic and fermionic systems [3,4,17], and at critical points described by conformal field theories [1,2].The latter rely on field-theoretic techniques using replicas and path integrals, typically in imaginary time, with complicated boundary conditions on fields and associated Green's function along the time direction [1][2][3].Such methods are often hard to implement for interact-ing systems.Hence, computation of entanglement entropy for interacting systems are only limited to small systems using exact diagonalization (ED) or for systems, like those without sign problem accessible via quantum Monte Carlo (QMC) simulations [18][19][20][21].
A promising path-integral approach, that circumvents the use of complicated boundary condition using known relation between reduced density matrix and Wigner characteristic function [22], has been recently proposed for bosonic systems in ref. 23.Motivated by this, here we develop a new field theoretic method to compute Renyi entanglement entropy for fermions.A similar field theory formalism for fermions has been developed independently by Moitra and Sensarma [24].To this end, we derive new representations of an operator and traces of product of operators in terms of fermionic displacement operators [25].The representations allow us to develop a very transparent fermionic coherent-state path-integral method with simple boundary condition to compute Renyi entropies of a sub-region of the system in terms of a fermionic version of the 'Wigner characteristic function'.The formalism can equally be applied to calculate subsystem Renyi entropy for thermal equilibrium state via imaginary-time path integral or non-equilibrium time evolution described via Schwinger-Keldysh field theory.The approach naturally transcends the effect of boundary condition into time-dependent selfenergy, which acts like 'kick' at a particular time.We show that the method immediately reproduces the known expressions for von Neumann and Renyi entanglement entropies for non-interacting fermions.The effect of the time-dependent self-energy can be implemented for interacting systems treated within standard perturbative arXiv:2004.04751v1[cond-mat.str-el]9 Apr 2020 and non-perturbative field-theoretic approximation and diagrammatic continuous-time Monte Carlo simulation [26].We elucidate this by deriving the subsystem second Renyi enetropy within two well-known approaches to treat correlated fermions -(a) strongly interacting large-N fermionic models based on Sachdev-Ye-Kitaev model [27,28], and (b) dynamical mean-field theory (DMFT) [29].
In the other major part of the paper we explicitly demonstrate the utility of the method by computing the second Renyi entropy (S (2) ) for subsystems in several large-N model in thermal equilibrium -(i) zero dimensional SYK model having infinite-range random four-fermion or two-body interaction with a non-Fermi liquid (NFL) ground state, (ii) SYK model with additional quadratic hopping between fermions having a Fermi liquid ground state, (iii) a generalized SYK model, the Banerjee-Altman (BA) model [30], having quantum phase transition (QPT) between SYK NFL and FL, and (iii) an extended system, a chain of SYK dots [31,32], describing an interacting diffusive metal.We compute the subsystem Renyi entropy at the large-N saddle point as a function of subsystem size for the thermal density matrix at a temperature T and extrapolate to T → 0 to obtain ground-state Renyi entanglement entropy.In all the above cases except (iii), we compare and contrast the results for the interacting systems with a corresponding non-interacting system where the SYK interaction is replaced by infinite-range random hopping.For the non-interacting models we numerically calculate the second Renyi entropy by numerical diagonalization for moderately large systems and compare with large-N field theoretic results.Moreover, we also compute subsystem Renyi entropy via exact diagonalization (ED) of manybody Hamiltonian for small N ∼ 8−16 for the interacting zero-dimensional models.We obtain the following important results using our method for the large-N models.
(1) We show that for the SYK model, in the N → ∞ limit, the zero-temperature residual entropy [27,[33][34][35] of the SYK NFL contributes to the Renyi entanglement entropy, thus making it impossible to recover the true quantum entanglement of the NFL ground state.Moreover, consistent with our numerical results, we analytically prove that the SYK model is maximally entangled when the relative size of the subsystem p → 0.
(2) We demonstrate how the T → 0 bipartite Renyi entanglement entropy crosses over from a NFL to a FL as function of the strength of hopping in the SYK model with random quadratic term.The results show that heavy FL are much more entangled than weakly-or noninteracting FL.
(3) In the BA model with a NFL-FL transition, we establish a precise connection between Renyi entanglement entropy and the residual entropy of the NFL and show that the entanglement entropy carries sharp signature of the underlying QPT.
(4) In the extended one-dimensional (1D) model of SYK dots, we compute Renyi entanglement entropy of a subsystem of length l.We find a crossover from S (2) ∼ log l to S (2) ∼ log(1/(l −2 + l −2 0 ) 1/2 ) with increasing l where the log l behavior, expected for gapless fermions, gets cut off by an emergent mean-free path l 0 in an interacting diffusive metal.
The paper is organized as follows.In sec.II A we derive the operator identities that form the basis of our formalism and discuss the connections of these identities with Renyi entropy of a subsytem.The general formulation for the equilibrium and non-equilibrium path integrals to compute the Renyi entropy is discussed in sec.II B. Sec.II C describes the application of the field-theory formalism to derive well-known formulae for Renyi entropy of non-interacting fermions.In secs.II D and II E, we develop the field-theory for Renyi entropy in several interacting large-N models based on SYK model, and in Hubbard model within DMFT approximation, respectively.We describe our analytical and numerical results for Renyi entropy in the large-N models and the comparison of the large-N results with that obtained from numerical exact diagonalization in sec.III.Additional details of the derivations of the operator identities, path integral formulations and their analytical and numerical implementations in various models to compute Renyi entropies are given in the appendices.

II. COHERENT STATE-PATH INTEGRAL FORMALISM FOR FERMIONS IN AND OUT-OF EQUILIBRIUM
A. Subsystem Renyi entropy, displacement operator and trace formula In this section, we consider a system with fermionic degrees of freedom and derive a useful expansion of an arbitrary operator in terms of the so-called fermionic displacment operator [25] and show that the expansion can be used to represent Renyi entropies of a sub-region of the system.A similar representation of Renyi entropy has been independently developed by Moitra and Sensarma [24].The Renyi entropy of a quantum system described by a density matrix ρ is obtained by dividing the system into two parts A and B (not necessarily equal) and defining a reduced density matrix, for region A. If ρ represents a pure state, then a measure of the quantum entanglement of region A with B can be obtained by evaluating the von Neumann entropy by tracing over the degrees of freedom in A. In practice however, the von Neumann entropy is often hard to calculate directly within field theoretic methods, and a more convenient measure of entanglement [1][2][3][4] is the n-th Renyi entropy, S A , defined as where the integer n > 1.The von Neumann entanglement entropy can be obtained by analytically continuing to n → 1.We refer to the above subsystem Renyi entropy throughout the paper simply as Renyi entropy for brevity.

The main difficulty in evaluating the above comes from the representation of Tr
] in a coherent-state path integral, since each factor of Tr B ρ leads to separate replicas which need to be connected via appropriate boundary condition when represented through Grassmann variables [1,3].To circumvent this difficulty, we derive the following operator expansion (see app.A) for an arbitrary operator F , where ξ ≡ { ξi , ξ i } and η ≡ {η i , η i } denote set of Grassmann variables with index i = 1, . . ., N A , e.g., referring to a set of sites that includes the support of the operator F on the lattice; d 2 (ξ, η) = i d ξi dξ i dη i dη i .Here the weight function f N (ξ, η) is given by The basis of above expansion in eqn.( 3) is formed by the fermionic displacement operators [25], much like more familiar bosonic counterparts [22], defined as where c † i , c i are the creation and annihilation operators on site i.The displacement operator generates the coherent state [25], |ξ = D(ξ)|0 by shifting the vacuum state and more convenient to be used in the path integral representation discussed in the next section.Eqn.
(3) offers a way to decompose a general operator F using only the normal-ordered displacement operators D N (ξ) and is one of the key results of this paper.An important corollary to the decomposition identity of eqn.(3) is to express the trace of the product of operators F and G as To derive the above, we have used the identities, 7) is the second key result of this paper and is crucial for deriving a path-integral representation for evaluating entanglement entropy as we discuss in the next section.We could also make the operator expansion, F = d 2 (ξ, η)f (η, ξ)Tr[F D(ξ)]D(η) and obtain the corresponding trace formula similar to eqn.(7), in terms of the displacement operator D(ξ) and weight function Using the trace formula (eqn.( 7)), e.g., the second-Renyi entropy, S A , can be conveniently expressed as In the above, ξ ≡ { ξi , ξ i } i∈A (and similarly for η), i.e. the displacements operators above only involve the fermionic operators in the region A. As a result, we have i.e., the expectation value of the operator D(ξ) evaluated for the reduced system is same as that obtained using the full density matrix.Eqn. ( 9), therefore, eliminates the need to calculate the reduced density matrix ρ A .The evaluation of ρ A is the difficult step for calculating entanglement entropy, as mentioned earlier.To proceed further, we define the 'normal-ordered' fermionic Wigner characteristic function [22,25] for the density matrix, and arrive at the final expression for the second-Renyi entropy for region A e −S (2) We also sometime use an analogous expression written in terms of D(ξ) by replacing D N (ξ) and the function f N (ξ, η) by f (ξ, η).This leads to the usual fermionic characteristic function [25], The higher-order Renyi entropies, S , can be found in a similar manner by repeated application of the trace identity in eqn.(7).In fact, as we show in sec.II C, a hierarchy for higher-order Renyi entropies can be derived that recursively expresses the characteristic function of higher moments of the reduced density matrix ρ A in terms of the lower order ones.

B. Equilibrium and non-equilibrium field theories for Renyi entropy
A path integral representation of the characteristic function, χ N (ξ), will naturally lead to a similar represen-tation for the second Renyi entropy of region A through eqn.(11).Therefore, we first derive the path integral for χ N (ξ) for -(a) thermal density matrix describing a system in equilibrium, and (b) time-evolving density matrix subjected to a general time-dependent Hamiltonian, e.g.describing a quench.

Path integral for thermal equilibrium
The density matrix for a system described by a Hamiltonian H under thermal equilibrium is given by with the partition function and inverse temperature β = 1/T (k B = 1).Inserting the identity operator d 2 c|c c| = I in the coherent state basis, we get with d 2 c n = i dc in dc in .Using eqn.( 6), the matrix element for the displacement operator is easily evaluated as Finally, following the standard methodology [36] of fermionic coherent state, e.g., evaluating − c 0 |e −βH |c 1 via Trotter decomposition β = N τ ∆τ and taking the continuum limit N τ → ∞, ∆τ → 0 we get as the path-integral representation of the characteristic function.The above equation offers an advantage over usual field-theoretic formalism [1][2][3] to compute entanglement entropy.The imaginary-time boundary condition for the Grassmann fields c i (τ ) and ci (τ ) are still antiperiodic, i.e. c i (τ + β) = −c i (τ ), irrespective of whether i belongs to the region A or not.Instead, the distinction between subsystem A and the rest of the system is encoded by the auxiliary fields ξ which only couple with the fields c i (τ ) and ci (τ ) for i ∈ A. Under time discretization, relevant for the numerical implementation discussed later, we have

Path integral for non-equlibrium evolution
The density matrix, at a time t, evolving under a timedependent Hamiltonian H(t) is given by where ρ 0 is the initial density matrix at time t 0 .The operator U (t 1 , t 2 ) is the unitary evolution operator asso-ciated with the Hamiltonian H(t) and defined as where T, T are the time ordering and anti-time ordering operators respectively.We use Schwinger-Keldysh closed time-contour formalism [37,38] to obtain a path integral representation for the Renyi entropy, e.g., S (2) (t), which is now time dependent and given by e −S (2) Here we have rewritten the trace identity (eqn.( 7)) (see app.A) to define the time-dependent characteristic function as for the sake of convenience in constructing the path integral representation.As in the standard Schwinger-Keldysh closed-time contour formalism [37,38], the last line in eqn.( 22) may be interpreted, from right to left, as starting from an initial density matrix ρ 0 , evolving forward in time (represented by + branch in fig. 1 As done often in the Schwinger-Keldysh formalism, the contour is extended to +∞ (see fig. 1).We incorporate this contour extension in our expression for the characteristic function.
Again following standard route [36,37], as in the thermal equilibrium case, we obtain a path integral representation for χ N (ξ, t), in terms of the entanglement Keldysh action such that Here δ C (z, z 1 ) is the Dirac delta function on C and takes into account both the time argument and the branch index (±).The coherent-state matrix element of the initial density matrix, c i (0, +)|ρ 0 | − c i (0, −) in eqn.(25) encodes the information of the initial state at starting time t 0 .For a thermal initial state, one can add an additional branch to C and redefine it as i.e., the contour is now comprised of a vertical imaginarytime branch of length β and the usual + and − branches, as shown in fig.1(b).The form of the action S C remains exactly the same as given in eqn.(24) such that where Z 0 is the partition function for the initial thermal state described by ρ 0 .We emphasize here that the Hamiltonian describing the ρ 0 and the one dictating the unitary evolution (see eqn. (20)) do not need to be same, e.g. in a quench [39][40][41][42].This distinction is incorporated by choosing the appropriate Hamiltonian in eqn.(24) for the imaginary-time and real-time branches.
In the next section we show that the equilibrium and non-equilibrium coherent-state path integrals can be used immediately to obtain well-known expressions for Renyi and von Neumann entanglement entropies in terms of the correlation matrix [3].

C. Renyi entropies for non-interacting fermions
We consider a non-interacting system described by the thermal density matrix where Z = Tr exp −β ij t ij c † i c j is the partition function.Using eqn.(17), the normal-ordered characteristic function in this case is given by The Gaussian structure of the fermionic-integral allows us to integrate the fermions to give where the matrix G, is the Green's function and is defined as the imaginary-time-ordered two point correlator and can be calculated using However, because of the delta functions in eqn.(30), we need to evaluate only G ij (0, 0 + ).This allows us to define a correlation matrix [3] C ij , as follows where C T is transpose of the matrix C. Therefore, the characteristic functions for an arbitrary non-interacting system is given by The Gaussian structure of the characteristic function is a direct consequence of the underlying non-interacting Hamiltonian.Using eqn.( 11) and the characteristic function above, the second Renyi entropy can be immediately evaluated by integrating out the auxiliary Grassmann variables ξ and η to get the well-known expression [3] for non-interacting fermions, where 1 is a N A × N A identity matrix and only the elements C ij of the correlation matrix (see eqn. ( 33)) with i, j ∈ A are involved in the above expression.In fact, the Gaussian form of the characteristic functions allow us to derive the expression for the higher Renyi entropies as well.To this end, we first define the higher-order generalization of the characteristic function (see eqn. ( 10)) A recursion relation of the form with χ n=1 = χ, can then be derived (see app.B 0 a), which obtains higher-order χ from lower ones.As alluded earlier, we use a Gaussian-ansatz of the form where to solve the recursion.The ansatz is consistent with the expression for χ obtained in eqn.(35), provided we set The ansatz lets us simplify the recursion, see app.B 0 b, and derive recursions involving λ n and A n , i.e.
The last equation can be rearranged (see app.B 0 b) into with X n defined as We find the following solves eqn.(43), see app.B 0 c.The n-th Renyi entropy is obtained by evaluating λ n , since and eqn.(2) implies From the recursion for λ n in eqn.( 42), and eqn.( 44), we see that Looking at the structure of X n in eqn.(45), we find that the inverse in X n cancels with the non-inverse term in X n−1 , so on and so forth, leaving at the end from which we find the n-th Renyi entropy to be in agreement with the known results [3].Taking limit n → 1 above, we get the von Neumann entanglement entropy As shown in app.C, in the non-equilibrium case, using the Schwinger-Keldysh path integral, discussed in sec.II B 2, we obtain the same expression for the n-the Renyi entropy S (n) A (t) as given in eqn.( 50) with a timedependent correlation matrix Here ρ(t) is the density matrix at time t.
Having reproduced the results for non-interacting fermionic systems, we now seek to set up the field theory for Renyi entropy in some well-known models for interacting fermions in the next sections.We discuss the formulation for two types of systems, large-N fermionic models based on SYK model [27,28] and repulsive Hubbard model treated within single-site DMFT [29].).Each dot is connected to its nearest neighbors via q hop -body terms, e.g.q hop = 1, as shown in the zoomed view (below).The highlighted part of the chain is chosen as the subsystem A for which entanglement-entropy is calculated.
SYK model for non-Fermi liquid state: The SYK model describes a zero-dimensional system having N fermion flavors or sites (see fig. 2(a)), that interact via an all-to-all or infinite-range Hamiltonian.To keep the discussion general, we use the q-body version [35] of the SYK model, SYK q (typically referred in the literature as SYK 2q ) described by the Hamiltonian where are properly antisymmetrized Gaussian random numbers with variance J 2 /qN 2q−1 (q!) 2 .Setting q = 2 in the above Hamiltonian gives back the original SYK model [27,28].The ground state for the SYK q model, for q ≥ 2, has been shown to be a non-Fermi Liquid (NFL), lacking quasi-particle excitations, with a fermion scaling dimension ∆ = 1/2q [35].Interestingly, the SYK NFL states possess a residual zero-temperature thermodynamic entropy S 0 in the limit N → ∞, e.g S 0 ≈ 0.464 for q = 2, which, as we shall see later, plays a crucial role in the interpretation of the entanglement entropy as well.The q = 1 model is a special case whose ground state is a non-interacting Fermi Liquid described by a single-particle semi-circular density of density of states (DOS).SYK model with quadratic hopping term for Fermi liquid state: The SYK model has also been generalized to describe strongly-interacting or heavy Fermi-liquids (FL) [32] by adding a random all-to-all hopping (q = 1) term to eqn.(52), i.e., where t ij are the Gaussian random numbers with variance t 2 hop /N (see fig. 2(b)).BA model for non Fermi liquid to Fermi liquid transition: We also study Renyi entropy in the model of ref. [30], H = H c + H ψ + H cψ , where The above model (see fig. 2(c)) has two species of fermions -(1) the SYK fermions (c), on sites i = 1, . . ., N c , interacting with random coupling J ijkl (eqn.(54a)) with variance J 2 /(2N c ) 3/2 ; and (2) the peripheral fermions (ψ), on a separate set of sites α = 1, . . ., N ψ connected via random all-to-all hopping t αβ (eqn.(54b)).The SYK and the peripheral fermions are quadratically coupled via V iα ; t αβ and V iα are Gaussian random variables with variances t 2 hop /N c and V 2 / N c N ψ , respectively.
The model is exactly solvable for N c , N ψ → ∞ with a fixed ratio p = N ψ /N c , that is varied to go through the QPT between NFL and FL at a critical value p = p c = 1 [30].The residual entropy density S 0 (p) of the SYK NFL continuously vanishes at the transition [30].
Lattice of SYK dots for interacting diffusive metal: Extension of the above zero-dimensional models to higher dimensions have also been achieved [31,32,47,48,67].These systems typically involve a lattice of SYK q dots, each having N fermion flavors, connected to their nearest neighbors via q hop -body terms and can be described using the general Hamiltonian where the coordinates x, x denote the position of the individual dot in the lattice or chain and xx represents the nearest neighbors (see fig. 2(d)).The amplitude t xx and J x are independent Gaussian random numbers with variances having the same form as that of the coupling J in eqn.(52); q → q hop and J → t chain for the t xx amplitudes.Setting q hop = q = 1 produces a non-interacting diffusive metal, while q hop = 1, q = 2 and q hop = 2, q = 2 results in a diffusive heavy-Fermi liquid [32] or a non-Fermi liquid [43], respectively.
1. Renyi entropy in the thermal state Here we derive the exact equations to evaluate the Renyi entropy for all the models of the preceding section at the large-N limit in thermal equilibrium.The evaluation of Renyi entropy for non-equilibrium evolution can be performed using the Schwinger-Keldysh path integral formalism of sec.II B 2 and is discussed in app.D. In order to access the second Renyi entanglement entropy of the ground state we use the thermal density matrix to perform the analysis and then take the T → 0 limit.The definition for 2-nd Renyi entropy (see eqn. ( 2)) implies, that we need to calculate the disorder averaged quantity where • • • represents disorder average over the amplitudes We also define the operator Z A = Tr B exp(−βH).We evaluate each term in eqn.( 56) separately.Using the trace identity (see eqn. ( 7)) we formally write down Tr A [Z 2 A ] as where as before ξ and η couple only to sites in the region A, and Tr[ZD N (ξ)] is the characteristic function for the partition function Z. Since we need to evaluate the logarithm of Tr A [Z 2 A ] (see eqn. ( 56)), we use the replica-trick to express the log as a product of disorder replicas The replicas, along with eqn.(57) and the integral representation in eqn.(17), allows us to write the path-integral for the SYK q model (eqn.( 52)) The Grassmann fields ciσa (τ ), c iσa and self-energy Σ σ b,σa (τ 2 , τ 1 ) for the fermions.We point out that the above analysis closely follows that of the thermal-field theory of the SYK model.The latter can be found in refs.[30,33,43,48,49], to name a few.The end result of introducing large-N fields is an integral of the form where the effective action S for the second Renyi Entropy is bilinear in the Grassmann-fields c, c, ξ, and η, i.e.
ξia ηia The Gaussian structure, allows us to integrate the ξ and η variables, changing the effective action to another bilinear of c and c where the matrix M is defined as and T .Interestingly, at this point we see explicitly, from eqn. (63), how the fields ξ, η extract the information regarding entanglement.Instead, of introducing complicated imaginary-time boundary conditions for fermions in sub-region A, the fields alter the self-energy by providing a time-dependent "kick" to the fermions, while leaving the fermions outside A untouched.All the fermions can now be integrated to produce Tr[Z 2  A ] r = D(c, c, Σ, G)e −N S[Σ,G] and the final effective action which depends only on the Green's function G and self-energy Σ.Here dτ 1,2 ≡ β 0 dτ 1 dτ 2 .The symbol p denotes the the ratio of subsystem size to total size, i.e.
and can take a value from 0 to 1.The symbols ∂ τ , Σ, appearing above, represents matrices having elements respectively.We evaluate the action S at the (disorder) replica diagonal and replica symmetric saddle point (i.e.Σ, G ∝ δ ab ) by minimizing with respect to G and Σ to get as the saddle point conditions.Here G is the matrix representation of the Green's function G and g, G are additional matrices that we have introduced to simplify the notation.Due to the similarities of the entanglement action S in eqn.(65) with the thermal free energy, we define the entanglement free-energy which allows us to express the 2-nd Renyi entropy(eqn.( 56)) as which depends on p and temperature β −1 .Note that we have used the fact that thermal free energy F (β) at the saddle-point is given by −β −1 ln Z and is equal to (69).The ground state entanglement entropy for an arbitrary subsystem size p can now be calculated by taking the limit The above analysis can also be carried out for the generalizations of the SYK model given in eqns.(53), (55) and (54).The equations to evaluate the Renyi entropy for the strongly-interacting Fermi liquid Hamiltonian in eqn.(53) and BA Hamiltonian in eqn.( 54) can be similarly obtained and are given in app.E. For extended systems of the kind described by eqn.(55), i.e. a ring of SYK dots connected to their nearest neighbors, we choose the subsystem as the l successive dots, e.g.x = 1, . . ., l (see fig. 2(d)) as the subsystem A. The saddle-point conditions for this arrangement are where the matrix M (defined in eqn.( 64)) adds to the site dependent self-energy -Σ (x) only for dots which belong to subsystem A, as indicated by the Kronecker delta function δ x∈A .The Green's function G (x) also become dependent on the site index x.The entanglement free energy can be calculated in terms of the space dependent G (x) and Σ (x) , i.e. x,σ where xx represents nearest neighbors.The presence of δ x∈A in eqn.(71) explicitly breaks translation symmetry and one needs to retain the Green's functions for each site x in order to calculate the entanglement entropy.
We obtain the Renyi entropy in sec.III for all the models discussed in this section by numerically solving the saddle-point equations, such as eqn.(67) and eqn.(71), and using eqn.(69).We discuss the iterative numerical algorithm to solve the saddle-point equations, in app.H 2. In this section we demonstrate the application of the path integral method of sec.II B for the Hubbard model.We develop the formulation to compute Renyi entropy within dynamical mean field theory (DMFT) [29], one of the successful approaches to treat electronic correlation and metal-insulator transition in Hubbard model.We formulate the DMFT for Renyi entropy in the single-site approximation.In the usual single-site DMFT [29] one reduces a lattice problem into an impurity coupled to a bath and the properties of the impurity and the bath are calculated self-consistently using the non-interacting dispersion.We extend the DMFT approximation to evaluate the path integrals developed in sec.II B. The method can be integrated with the continuous-time quantum Monte Carlo (CTQMC) impurity solver [26] and extended to the cluster implementations [68].Here we only discuss the formulation of the problem within DMFT, the numerical implementation will be discussed in a future publication [69].
We consider the Hubbard model with hopping amplitudes t ij between sites on a lattice, µ the chemical potential and U the on-site repulsive interaction between fermions with opposite spins σ =↑, ↓.
Here n iσ = c † iσ c iσ and n i = σ n iσ are the electronic number operators.
We discuss the second Renyi entropy of a subsystem for equilibrium state described by the thermal density matrix (eqn.( 13)).The method can be easily extended to non-equilibrium situation (sec.II B) via non-equlibrium DMFT [70].As in eqn.(56), we compute where Ω (2) = −T ln Z (2) (Z (2) ≡ Tr A Z 2 A ) and Ω = −T ln Z are relevant grand potentials.Again, using the trace identity (eqn.( 57)) we obtain the path integral, Here we have the auxiliary Grassmann fields ξ = { ξiσα , ξ iσα } i∈A with spin (σ) and entanglement replica (α = 1, 2) indices, and As in the case of chain of SYK dots in sec.II D 1, due to the entanglement subdivision, the above action breaks translational invariance.Hence, we use a single-site inhomogeneous DMFT to reduce the lattice problem (eqn.( 73)) into impurity problem for site i, Z The auxiliary Grassmann fields { ξσα , ξ σα } only appears for sites belonging to the A region; The above impurity action could be derived from eqn. (75) using the cavity method [29].To simplify the notations, we have assumed a paramagnetic state with spin symmetry, and hence the dynamical 'Weiss mean field' is independent of spin and is given by A 2×2 matrix in the replica space, where I denotes a (2× 2) unit matrix and ∆ i (τ, τ ) is the hybridization function which carries the information about the lattice and needs to be found self-consistently (see below).For i ∈ A, we integrate out the entangling Grassmann fields ξσα , ξ σα to get Z (2) i = D(c, c)e − Si and the effective action where the 2 × 2 matrix M is given in eqn.(64).The impurity Green's function G i (τ, τ ) can be obtained from the impurity Dyson equation, where Σ i,αγ (τ, τ ) is the impurity self-energy that needs to be computed using an appropriate impurity solver, e.g.iterative perturbation theory (IPT) [29] or CTQMC [26].
The hybridization function can be obtained in terms of cavity Green's function as Here G (i) (τ, τ ) is the Green's function with i-th site removed from the lattice.The cavity Green's function can be obtained from the lattice Green's functions as Where G(τ, τ ) is the full lattice Green's function.To obtain the full lattice Green's function, within single-site DMFT approximation, one assumes that the lattice selfenergy is local and is given by the impurity self-energy.Then, the lattice Green's function is obtained from the Dyson equation The above equation can be used to obtain the hybridization function in eqn.(82) and thus closes the DMFT selfconsistency loop.The numerical solution of the above DMFT equations is rather involved since both space and time translation invariance are broken, the former due to the choice of subregion A and the latter due to the time-dependent self-energy kick used to extract the entanglement.The numerical solution will be discussed in a future work [69].
We conclude this section by briefly mentioning the important step to compute the grand potentials Ω (2) and Ω appearing in eqn.(74) for the second Renyi entropy.The grand potentials can be obtained via coupling constant integration discussed in detail in app.F. In this method, one considers a modified Hamiltonian H λ = H 0 + λH 1 , where H 0 is the non-interacting part in eqn.( 73) and H 1 is the Hubbard term with the coupling constant 0 ≤ λ ≤ 1; λ = 0 is the non-interacting limit and the λ = 1 is the interacting Hamiltonian of interest.The grand potentials could be obtained as Ω (2) = Ω (2) (0) The integrands above can be calculated by computing the Green's function G (λ) via DMFT for each λ.One important point here is that the equal-time Green's function appearing in the expression for Ω (2) needs to be calculated at the imaginary-time instant where the fermionic source fields are inserted in the path integral, in our case τ = 0.In contrast, τ in the expression of usual thermodynamic grand potential Ω is arbitrary.Ω (2) (0) and Ω(0) are the grand potentials for the non-interacting system (λ = 0) which can be calculated directly using the results of sec.II C.

III. ANALYTICAL AND NUMERICAL RESULTS FOR RENYI ENTROPY IN LARGE-N FERMIONIC MODELS OF FERMI LIQUIDS AND NON-FERMI LIQUIDS
With the formalism for calculating entanglement entropy for large-N systems in place, we proceed forward and discuss the solutions obtained for the models introduced in sec.II D. In this paper, we only describe the results for the 2-nd Renyi entropy obtained from the thermal field theory in sec.II D 1.The results for higher-order Renyi entropies and time-evolution of Renyi entropy in non-equilibrium situations will be communicated in a future work [71].

A. Non-interacting large-N model with disorder
To set the stage and bench mark our large-N field theory formalism of sec.II D 1 we first discuss the Renyi entropy in a non-interacting model.For q = 1, the model defined in eqn.( 52) describes a non-interacting system connected with disordered random hoppings having a variance t 2 hop , and hence is amenable to solution in multiple ways.Therefore, we compare the temperature dependent entanglement entropy, defined in eqn.(69), calculated using -(a) the correlation matrix approach applicable for non-interacting system as discussed in sec.II C, where the correlation matrix is evaluated by diagonalizing the non-interacting Hamiltonian, (b) using manybody exact diagonalization (ED) (see app.H 1) and (c) using the large-N thermal field developed in sec.II D 1.For cases (a) and (b), the large-N entropy is calculated by explicitly averaging over multiple disorder realizations.
Fig. 3 shows the variation of 2-nd Renyi entropy with subsystem size p for multiple temperatures calculated using the above methods.We get an excellent match for all three cases.Particularly remarkable is the fact that the large-N results almost exactly match with that obtained with ED for quite small system size (N = 10), implying that 1/N corrections to entanglement entropy is rather small.This feature persists for the interacting large-N models discussed in the next sections.We find that for small subsystem size, i.e. p → 0, the Renyi entropy grows with p with a slope of ln(2) (see dashed line in fig. 3) indicating a maximal entanglement which can be argued based on the counting of degrees of freedom [63,72].In the limit of large subsystem size, i.e. p → 1, the entanglement entropy attains a finite non-zero value that goes to zero as T → 0, see fig. 3  FIG.3. Renyi entropy in non-interacting large-N model: Second Renyi-Entropy, S (2) , for the non-interacting model with random hoppings shown as a function of subsystem fraction p for several values of temperature T .The predictions obtained from large-N saddle-point theory (points with lines) show an excellent agreement with results (represented by color matched squares) obtained using the correlation matrix, as well as exact diagonalisation (ED) calculations (circles) performed for N = 10 sites.All three techniques show a maximal entanglement scaling of Renyi entropy for small subsystem sizes (p → 0) with a ln(2) coefficient.Inset, gives the temperature dependence of the Renyi entropy at p = 1, i.e. S (2) (p = 1), which goes to zero as temperature is decreased indicating that the thermal density matrix, used for calculating entanglement, approaches the ground state for the model.
mixed-state density matrix picks up contribution from usual thermal entropy, which is only of statistical origin.However, as temperature is lowered (β → ∞) the thermal density matrix approaches that of the ground-state, i.e. ρ 2 → ρ.For any pure-state, in this case the ground state, the entanglement entropy goes to zero as the subsystem size approaches the size of the parent system.

B. SYK model: Renyi entropy of a non-Fermi liquid
We now move on to large-N interacting Hamiltonians, starting with the original q = 2 SYK model described in eqn.(52).The model connects fermions, on N sites, via all-to-all two-body interactions with random complex amplitudes having variance J 2 , which we set to unity for this section.There have been several numerical and analytical works on entanglement entropy [43,[63][64][65][66]72] in the SYK model for various different contexts.The analytical approaches have used standard replica pathintegral appraoch for Renyi entropy.Among these works, Gu et al. have studied growth of bi-partite Renyi entanglement entropy in a chain of SYK dots starting from a thermofield double state [43].They made a diagonal

FIG. 4.
Renyi entropy in SYK model: (a) The second-Renyi entropy (S (2) ) for the SYK (q = 2) model as a function of subsystem size p for temperature values T = 0.50, 0.20, 0.05.Predictions from large-N formalism (points with lines) match exactly with ED calculations (represented by circles, squares and triangles) for subsystem sizes p < 0.4 (both showing a volume-law scaling with a ln(2) coefficient), but deviate significantly at larger p values due to the influence of residual-thermal entropy S0 of the SYK model.The S (2)  vs. p curve (black line with crosses) for the ground-state of the model, obtained using ED, is also shown for comparison.(b) Renyi entropy for p = 1, i.e. S (2) (p = 1), approaches S0 linearly when temperature T goes to zero, as seen explicitly from a linear fit (line) to the large-N data (circles) and recovering an intercept equal to S0 ≈ 0.464.(c) The match between large-N prediction (filled circles) and ED (empty circles), for S (2) (p = 1), is recovered at higher temperatures, see text for explanation.All ED calculations were done for N = 12 sites.approximation in the entanglement replica space while treating inter-dot coupling perturbatively.In our work, we find the (entanglement) replica off-diagonal terms to be very important.Ref.65 has used similar methods to look into ground-state entanglement between two quadratically coupled SYK dots.During the preparation of our manuscript, we became aware of a very recent work [66] on subsystem Renyi entropy in the thermal ensemble for the SYK model using standard approach with the boundary conditions.Overall, their results compares well with our results on SYK model obtained using the new field-theoretic method developed in this work.
As in the non-interacting large-N model case, we compare the entanglement entropy obtained via large-N thermal field theory with that obtained via ED.However, the Renyi entropy cannot be calculated in the interacting case using the correlation matrix approach.Fig. 4(a) gives the 2-nd Renyi entropy with subsystem size p for several values of T .We find that, yet again, for all tem-perature values the results from large-N and ED calculations agree really well, when p < 0.4, even for exact diagonalization performed for N = 12 sites.Also as expected, volume law with a ln(2) coefficient still holds as p → 0 for all temperature values.Interestingly though, as T → 0, the large-N result deviates substantially from the ED result when p ≥ 0.4 and especially when p → 1.While the entanglement entropy, obtained by ED, goes to zero as expected, the large-N result seems to converge linearly to a finite value when T is reduced as shown in fig.4(b), where we plot S (2) (p = 1) for low values of T .In fact, when extrapolated to T → 0, we find the intercept to be equal (within numerical accuracy) to the thermal residual entropy S 0 ≈ 0.464 • • • of the SYK model.Remarkably, the residual entropy has found a way to influence the Renyi entropy for the ground-state of the SYK model in the large-N limit!An explanation is provided when we look back at the expression for Renyi entropy in eqn.(69).Since by definition lim p→1 F EE (p, β) = F (2β), expanding the entanglement free energy to first order in T (= β −1 ) we find i.e. S (2) (p = 1) is indeed equal to the residual entropy S 0 .
The above analysis suggests that the limit N → ∞ and T → 0 do not commute.It is established that the residual entropy for the SYK model is a consequence of exponentially small in N level spacings arising from the large-N limit.Taking the large-N limit first also prohibits "direct" access to the ground-state quantum entanglement when temperature is reduced, since any T however small cannot resolve these exponentially small level spacings.We therefore argue that temperatures higher than the typical level spacings obtained from ED should "mimic" this large-N effect and predictions from large-N thermal field theory should match with finite-N ED calculations.To test this hypothesis, we plot the p = 1 value for 2-nd Renyi entropy obtained from large-N and finite-N ED calculations for higher values of temperature in fig.4(c).Indeed, we find that the results match perfectly for T ≥ 0.5, thereby validating our intuition about the role of residual entropy in entanglement at large-N .

C. SYK model with disordered hopping: Renyi entropy of an interacting Fermi liquid
We now discuss the behavior of entanglement in an interacting Fermi liquid, like the one described by the Hamiltonian in eqn.(53).The addition of a q = 1 quadratic term t ij c † i c j has been shown [30,72] to change the SYK NFL ground state to an interacting zerodimensional Fermi-liquid.We set t hop = 1, J = 1 for the variances of t ij , J ij in eqn.(53) and compare the 2-nd (a) Subsystem size (p) dependence of second Renyi entropy, S (2) , in the zero-dimensional Fermi-liquid model (defined in eqn.( 53)) with SYK interactions, for several values of temperature T .The strength of hopping -t hop and interactions -J are both set to one.The predictions from the large-N formalism (points) and results from exact diagonalisation (circles) match for all values of p and T .Both predict a volume-law scaling with a ln(2) coefficient.(b) The Renyi entropy for p = 1 (represented by circles) shows a linear dependence on T (for values less than T = 0.15) and goes to zero as T → 0, indicating an approach to a ground-state with no residual-thermodynamic entropy.(c) The zero-temperature bipartite second-Renyi entropy, S (2) (p = 0.5, T → 0), shown as a function of increasing hopping strength t hop (J = 1), with t hop = 0 being the SYK-NFL limit.As hopping increases, entanglement decrease from the a rather large value (≈ 0.35) and saturates to a value (≈ 0.15) close to that of the non-interacting model (see fig. 3).
Renyi entropy obtained from large-N theory with ED in fig.5(a).Quite encouragingly, we find an excellent match between the two for all values of temperatures!As expected, the volume law for entanglement is shown to hold as p → 0 due to the all-to-all nature of interactions and hoppings.Interestingly this time, unlike the SYK NFL case, the entanglement entropy for p = 1 goes to zero as temperature approaches zero (see fig. 5(b)).Looking back at the analysis in eqn.(87), S (2) (p → 1) → 0 implies that the FL ground state does not possess any residual entropy.Indeed, we find in literature [73] that the presence of a quadratic term in eqn.(53) drastically reduces the exponentially large density of levels near the ground state thereby reducing the residual entropy to zero.
It is interesting to explore how entanglement entropy changes as the relative strength of interactions is increased with respect to the quadratic term, going from a non-interacting system to a NFL, via heavy Fermi liquid states for t hop J.To this end, we set J = 1 and track the bipartite Renyi entropy, i.e. S (2) (p = 1/2), as the value of t hop is increased from 0 to 1.The bipartite Renyi entropy is the maximal value that can be attained for any homogeneous system.Fig. 5(c) shows the result of this exercise.When t hop = 0, the model is described by purely SYK type interactions and therefore bipartite Renyi entropy attains the maximum possible value ≈ 0.5 ln(2).Interestingly, as t is increased the entanglement decreases continuously, indicating that a FL of this type becomes less entangled when the relative strength of interactions is reduced.Finally, for values of t hop ≥ J the bipartite entanglement entropy saturates to around 0.16, a value close to that of non-interacting large-N model (see fig. 3), implying that interactions become irrelevant in this limit.

D. BA model: Renyi entropy across a non-Fermi liquid to Fermi liquid transition
One of the important aspects of the putative duality between SYK model and black holes in quantum gravity [27,[33][34][35] is the residual entropy and its possible connection with the black hole entropy.Hence, it is interesting to ascertain if there is any relation between the residual entropy and the ground-state entanglement entropy of SYK NFL.We saw that the residual entropy inevitably appears in the Renyi entropy when the large-N limit is taken followed by the T → 0 limit in the pure SYK model.Here we discuss the BA model [30] which helps us to make the connection between the residual and entanglement entropy more explicit by tuning a QPT between a SYK NFL and a FL.
The BA model is defined in eqn.(54).A T = 0 transition between the chaotic NFL fixed point and non-chaotic FL fixed point is achieved by tuning the ratio of sites p = N ψ /N c .The quantum-critical point (QPT) occurs at p = 1, below which the c-fermions behave as a SYK-NFL and the whole system has a finite p-dependent residual entropy, S 0 (p) = [(1 − p)/(1 + p)]S 0 , that goes to zero at the critical point [30].Motivated by this, we ask whether any such sharp features exists in the Renyi entropy.We analyze the 2-nd Renyi entropy for the model at the large-N saddle point (see app.E) using the following two choices of subsystem A -(a) when A is made up of all the c-fermions that have the SYK type interactions, and (b) when the A is comprised of the all the non-interaction ψ-fermions.For our calculations we set J = t hop = V = 1 (see eqn. ( 54)).
The subsystem-Renyi entropy for the c-fermions, S ψ .However, yet again, we find that residual entropy S 0 (p) has managed to influence en-tanglement entropy as seen from S (2) c (represented by circles) approaching a value close to S 0 = 0.464 of the SYK model when p → 0. On the other hand, the S (2) ψ (triangles) goes to zero in the same limit in accordance with the volume-law maximal entanglement, p ln(2).To uncover the signature of the underlying QPT, we extrapolated S (2) c and S (2) ψ to T → 0 and plot the results in fig.6(b).Remarkably, we find that around the critical point, p = 1, the Renyi-entropy for the fermions become exactly equal , i.e., S ψ and remains equal for all values p > 1. Motivated by this, to further understand the relationship between the residual entropy S 0 (p) and the Renyi-entropies ψ , we compare S 0 (p) with the difference S (2) ψ , which approaches S 0 for p → 0 and also vanishes continuously at p = 1.Quite encouragingly, we find the match to be excellent, as seen in fig.6(b).Thus the difference between the Renyi entropies of the two subsystems directly corresponds to the residual entropy and carries the signature of the underlying QPT.
We also look into the mutual information (per site) between the two species of fermions, defined as cψ .I(p) has a broad peak around the critical point as shown in fig.6(b) inset.Indeed, the underlying QPT has left an imprint in the entanglement entropy!E. Higher dimensional large-N models: Renyi entropy in an interacting diffusive metal We now discuss one of the main new results obtained using the path-integral formalism developed in this work.We compute the Renyi entropy of a sub-region in an interacting diffusive metal.To this end, we consider the higher-dimensional generalizations of the SYK derived large-N models, which were discussed in sec.II D and can be described by a Hamiltonian with the general form given in eqn.(55).The large-N equations needed to determine entanglement entropy were derived in eqn.(71) for one-dimensional large-N models.Here we quantify how well our the formalism performs in predicting the Renyi entropy for these systems.As mentioned earlier (see the discussion near eqn.( 71)), naturally solving the saddle-point equations for the extended system of the one-dimensional chain is computationally more expensive than the zero-dimensional case since we need to retain the Green's function G (x) σσ (τ 1 , τ 2 ) for each SYK dot in the chain.In particular, the time complexity for solving the equations scales linearly with the number of dots which we denote as N dots .However, this increase in time complexity does not prevent us from accessing entanglement for large chains and we are able to perform calculations for values of N dots which are well beyond the reach of ED, even for the non-interacting large-N models for which we use the correlation matrix approach of sec.II C. Therefore, we adopt the following strategy -first, we compare the results of our large-N field-theoretic approach with that obtained from correlation matrix for non-interacting systems for moderate system sizes, and then explore the behavior of entanglement for large non-interacting and interacting systems.The interacting systems are much larger than that accessible via ED.
We set q hop = q = 1 in eqn.(55), making the model non-interacting, and look at two configurations of hopping amplitudes t chain = J = 1 and t chain = 1, J = 0. Also, we set N dots = 10 and N = 10, i.e. a moderate number of flavors per dot to perform exact numerics using the correlation matrix approach of sec.II C.This allows us to access subsystems sizes for p = 0.0 − 1.0 in steps of 0.1.The result for the above two scenarios are shown in Fig. 7(a) and (b) respectively.We find, like before, an excellent agreement between large-N predictions and exact numerical calculations, with the results for Renyi entropy from both the techniques collapsing on The model is obtained by setting q = q hop = 1 and intra as well as inter-dot hopping strengths t chain = J = 1 in eqn.(55).(b) Renyi entropy vs. l for the non-interacting model when t chain = 1, J = 0.In both (a) and (b), large-N predictions (points joined by lines) are in agreement with the results obtained from exact numerics (color matched circles) performed using correlation matrices.(c) Extrapolated zero-temperature Renyi entropy as a function of subsystem size l for a ring of SYK dots, (fig.2(d)), connected by random hoppings.We set the interaction strength J = 1 within each SYK dot as well as the strength of hopping between dots t chain = 1.The number of dots in the chain, N dots , is varied from 20 − 100 dots (well beyond the scope of any exact numerics), to verify that the growth of entanglement for small subsystem sizes becomes independent of N dots for N dots > 20.(d) A close up, showing the growth of S (2) near small subsystem sizes, i.e. l ≤ 20, in the interacting SYK-chain (represented by squares) and the non-interacting chain (represented by pentagon) discussed in (a), for a large value of chain size N dots = 50.Both curves can be fitted (lines) with the analytical function S (2) ∼ log(1/(l −2 + l −2 0 ) 1/2 ), where l0 is an emergent "mean-free path" length scale which takes a value between ∼ 4 − 5 in this case.
top of each other for multiple temperature values.
Having established the agreement of our formalism with exact numerics, we move on to larger chain sizes involving interacting as well as non-interacting dots.In particular, we keep q hop = 1 and study chains with either q = 2 or q = 1 body interaction for the dots (see eqn. ( 55)).The resultant ground state in these cases is known to describe either a strongly-interacting diffusive Fermi-liquid [32] or a non-interacting (but still diffusive) FL, both of which are gapless phases.The entanglement scaling with subsystem size (l) for a FL have been argued to violate the area law.For e.g., in 1-d translationally invariant gapless Fermi system, conformal field theory (CFT) predicts [2] that the von-Neuman entanglement entropy S (1) (obtained by taking the limit lim n→1 S (n) ), scales with the universal form S (1) = (c/3) ln(l) + const.
where c is the central-charge of the underlying CFT.Following which, an exact expression for S (1) , of the form S (1) ∼ l d−1 log l, was obtained for non-interacting FLs in arbitrary dimensions d [74].Later on, it was argued that the n-th order Renyi entropy for non-interacting as well as interacting FLs should follow S (n) = 1 2 (1+ 1 n )S (1) , and therefore should have the same scaling with system size [75,76].The above scaling form has been studied and confirmed, for systems with and without interactions, using several numerical approaches [77].Therefore, we may naively expect to recover this log l scaling for our model of a one-dimensional chain of connected dots.
To verify this, we evaluate the 2-nd Renyi entropy, in the T → 0 limit, as a function of subsystem size l for increasing values of N dots and plot the result in fig.7(c).We find that the curves converge as the chain size increases and entanglement entropy growth near l = 0 becomes independent of the total number of dots in the chain.Surprisingly, we find that Renyi entropy quickly saturates to a finite value as subsystem size is increased, a behavior inconsistent with log l scaling, which poorly fits the the entanglement growth curve in our calculation.Instead, we find that a "modified" growth function fits our data rather well, see fig.7(d).This indicates the presence of an emergent length scale l 0 in the gapless system!This could be explained by realizing that the FL state is obtained by connecting sites (dots) with random hoping amplitudes (for e.g.t ijxx in eqn.( 55)).The latter will induce an effective mean free path for the quasi-particles in the system and therefore will cut off the growth of entanglement beyond a length scale l 0 .Indeed, if we identify l 0 as the mean-free path and take l 0 → ∞, we recover the expected log l growth of entanglement.However, we note that a definition of mean free path in a system with random inter-site hopping (with zero mean) is somewhat ambiguous due to the difficulty in identifying a 'Fermi velocity'.Interestingly, the growth function in eqn.(88) has also been suggested in ref. 78 based on numerical studies of system size scaling of entanglement in a disordered non-interacting system, and somewhat heuristic hydrodynamic arguments for diffusive FL.On the contrary, our results explicitly demonstrate the scaling law of eqn.(88) in a interacting diffusive metal for system sizes much beyond any other numerical techniques.
The coefficient of the growth function in eqn.( 88) is also expected to contain the information of effective central charge [78] of the underlying CFT.As evident in fig.7(d), the coefficient of the growth function changes with interaction strength.These aspects will be studied in detail in a future work [71].

IV. CONCLUSION
In this work, we have proposed and developed a new equilibrium and non-equilibrium field theory method to compute Renyi entanglement entropies for interacting fermions.The basis of the field-theory formalism relies on a new operator identity that we have derived here.The path-integral techniques is an alternative, and maybe complementary, to the existing path-integral methods that typically requires complicated boundary conditions on the fields for computing Renyi entropy.Our methods rigorously transforms the complex boundary conditions into time-dependent self-energies while keeping the standard anti-periodic boundary conditions on the fermionic fields as in a usual coherent-state path integral.As a result, the path-integral formalism could be easily incorporated within standard diagrammatic and field-theory techniques and approximations, e.g.mean-field theory and random-phase approximation (RPA), in a very transparent manner.We demonstrate this by formulating the field theory for Renyi entropy in Hubbard model within DMFT approximation, and in several interacting large-N fermion models of current interest.Using the methods, we obtain several new results on entanglement entropy of FL and NFL states in zero-dimensional large-N models based on SYK model.Our results reveal intriguing connections between residual entropy of SYK NFL and ground-state entanglement entropy in the large-N limit.We also obtain system-size scaling of entanglement in an interacting diffusive metal.
The DMFT formulation developed here for Renyi entropy would be very useful to understand quantum entanglement in correlated systems.The formulation could be easily integrated with standard impurity solver, e.g.CTQMC [26].Entanglement entropy has been previously studied within CTQMC treating the interaction correction perturbatively for weak interaction [21].Our method is non-perturbative and can be used in the strongly interacting regime.The field theory is also amenable to standard renormalization group (RG) methods that may be advantageous for deriving analytical results on entanglement entropy in interacting fermionic systems.A protocol, inspired by the "kick" interpretation, may also be designed to measure the Renyi-entropy in quantum circuits.
Integrating out ξ and η, we get . Equating the answer to the Gaussian ansatz for n + 1, we get Implying the following recursion relations for λ n and A n In the last line, we have chosen to simplify the expressions, by assuming that A n at the end will be well behaved functions of the matrix A, and therefore commutes with A, allowing us to treat all matrices as scalars.We rewrite the recursion for A n as The last line gives the expression reported in the main text as eqn.(43).

c. Solving the recursion
We find the following solution B10) solves eqn.(43).To prove it, we first define B = C T (and treat everything like scalars as explained earlier) to simplify the notation and then plug the above solution into the R.H.S of the recursion (eqn.( 43)) and simplify as shown below to find the answer to be X n+1 .Therefore, eqn.(B10) solves the recursion in eqn.(43).
Appendix C: Keldysh formulation for non-interacting systems For a non-interacting system undergoing nonequilibrium time evolution starting from an initial Gaussian state, e.g.described by the thermal density matrix of eqn.(28), the characteristic function (eqns.(23), (24)) can be obtained in the following form where G ij (z 1 , z 2 ) is the contour-ordered single-particle Green's function that encodes the details of the nonequilibrium process, e.g. for time-dependent Hamiltonian H(t) = ij t ij (t)c † i c j with time-dependent hopping t ij (t), G −1 ij (z 1 , z 2 ) = ((i∂ z1 + µ)δ ij − t ij (t 1 )) δ C (z 1 − z 2 ) (z 1 = (t 1 , ±)).We can integrate out the fermions to get   as discussed in sec.II C. The above is of the same form as obtained in eqn.(36) for the thermal case.This implies a similar recursion relation for higher Renyi entropies as the one given in eqn.(38) can be derived and solved in exactly the same way to give the final expression reported in eqn.(50) of the main text.
Appendix D: Non-equilibrium field theory for Renyi entropy in SYK model The imaginary-time thermal field theoretic formulation discussed in sec.II D for SYK model and its extensions allows us to access ground-state entanglement in these large-N model.In this section we discuss the large-N Schwinger-Keldysh formulation which will enable us to track the time evolution of entanglement entropy under non-equilibrium situations [54][55][56][57][58][59][60][61][62].We demonstrate this by deriving the Swinger-Keldysh action and saddle-point equations for out-of-equilibrium evolution in the SY K q model (eqn.( 52)) when the strength of interactions is varying with time, such that (D1) Here J(t) is an arbitrary function of time.Also, we take the initial density matrix as a thermal ensemble for the SY K q model prepared for some initial configuration of J ijkl s.We use the path-integral derived in eqn.(27) and a modified form of the contour defined in eqn.(26) to develop the formulation.Instead, of stretching the realtime branches to +∞ (see fig. 1(b)) we stop it at t, the time at which entanglement is to be measured.Therefore, the new contour is now defined as C = [t 0 − iβ, t 0 ) ∪ [t 0 , t] ∪ (t, t 0 ], (D2) where the imaginary-time contour remains same as that shown in fig.1(b).
The steps are similar to the thermal case of sec.II D, but with the imaginary-time τ generalized to the contourvariable z, and we end up with an intermediate replicadiagonal Keldysh-action q (G σ σ (z 1 , z 2 )) q , (D4) where the contour delta function δ C (z − − t) is nonzero when z approaches the measurement time t along the − branch of the contour and δ C (z − t) is nonzero when z approaches t from the + branch of the contour, see eqn.(D2).The large-N field G is upgraded to a contour version and is defined as The rest of the symbols have the same meaning as the ones defined in eqn.(59).We emphasize, that the time dependence is explicitly encoded in the function J(z), which is equal to J(t) when z ∈ [t 0 , t] ∪ (t, t 0 ] part of the contour, and a static quantity J 0 when z ∈ [t 0 −iβ, t 0 ).If we integrate ξ,η first, in eqn.(D3), and then the fermionfields ciσ (z), c iσ (z), in that order, we arrive at the final expression for the entanglement-Keldysh action − J(z 1 )J(z 2 ) 2q G(z 2 , z 1 ) q G(z 1 , z 2 ) q (D6) in the large-N limit.The symbols ∂ z , Σ represents matrices having elements ∂ z1 δ C (z 1 − z 2 )δ σσ , Σ σσ (z 1 , z 2 ) respectively, and the matrix M is defined as Eqn. (D6) is the time-dependent generalization of eqn.( 65) that was derived in the context of thermalfield theory.And in the same manner, the saddle-point easily.H 1 is the interacting part and we have multiplied by a real variable λ, which will be integrated over; The Hamiltonian H λ = H 0 + λH 1 and λ = 0 is the noninteracting limit and λ = 1 is the interacting Hamiltonian we are interested in.n = 2 correspond to the second Renyi potential and n = 1 the usual grand potential.We obtain Tr B e −β(H0+λH1) λH 1 . (F2) For n = 2, using the trace formula (eqn.( 7)) and integrating both sides over λ from 0 to 1 we get Ω (2) (λ) =Ω (2) where Ω (2) (0) is the Renyi grand potential for the noninteracting system.As in eqn.(75a), one can construct a path integral representation for the above using the generating function Z (2) (λ) and obtain Ω (2) =Ω (2) (0) + .

(F4)
Here α is the entanglement replica index.The important point to note here is that the expectation of λH 1 with respect to Z (2) has to be calculated at the same time where D N (ξ) is inserted in the path integral i.e. at τ = 0 in our case.Naively, it seems that the above requires the evaluation of a four-point function.However one can write λU i n i↑α (0)n i↓α (0) in terms of the single-particle Green's function.To show this, we use the Heisenberg equation of motion [79] dc

FIG. 2 .
FIG.2.Large-N models: (a)The q-body SYK model (SYKq) for complex-fermions shown for q = 2 interactions, that scatter a pair of fermions, for e.g., from sites k, l to sites i, j with a random complex amplitude J ijkl .(b) Model for a zero-dimensional interacting-Fermi liquid obtained from the SYKq model by adding quadratic (one-body) hopping terms (represented by dashed lines) on top of the existing interactions.(c) The BA model for NFL-FL transition, obtained by joining a non-interacting dot of peripheral ψ-fermions (blue dots joined by lines) to the SYK-dot of c-fermions via random hopping amplitudes Viα.(d) A chain of SYKq dots with periodic boundary conditions (P.B.C).Each dot is connected to its nearest neighbors via q hop -body terms, e.g.q hop = 1, as shown in the zoomed view (below).The highlighted part of the chain is chosen as the subsystem A for which entanglement-entropy is calculated.
(τ ), appearing in the path-integral above and representing the fermions in the model, are labeled by several indices apart from the site indices i, j = 1, • • • , N -a imaginary-time coordinate τ , a disorder replica index a = 1, • • • r and an entanglement replica index σ = 1, 2.Here ξ i1a = ξ ia and ξ i2a = η ia , indicating the Grassmann-field origi-nates either from the characteristic function Tr[ZD(ξ)]or Tr[ZD(η)], respectively in eqn.(57).The Grassmannvariables ξ and η have the same indices as c,c except for the time-coordinate τ .The subsystem or the region A in the zero-dimensional model is defined as any of the N A sites, e.g.i = 1, . . ., N A , out of the total N sites.Having obtained the path-integral, we now seek to evaluate the integral at the large-N saddle point, i.e.N → ∞.This limit can be accessed by taking the disorder average over J i1•••iq;jq•••j1 s and then introducing the large-N field

E
. Renyi entropy in Hubbard model: Dynamical mean-field theory (DMFT) inset.At finite temperature the Renyi entropy of the subsystem for a thermal

FIG. 5 .
FIG. 5. Renyi entropy of interacting Fermi-liquid:(a) Subsystem size (p) dependence of second Renyi entropy, S(2) , in the zero-dimensional Fermi-liquid model (defined in eqn.(53)) with SYK interactions, for several values of temperature T .The strength of hopping -t hop and interactions -J are both set to one.The predictions from the large-N formalism (points) and results from exact diagonalisation (circles) match for all values of p and T .Both predict a volume-law scaling with a ln(2) coefficient.(b) The Renyi entropy for p = 1 (represented by circles) shows a linear dependence on T (for values less than T = 0.15) and goes to zero as T → 0, indicating an approach to a ground-state with no residual-thermodynamic entropy.(c) The zero-temperature bipartite second-Renyi entropy, S(2) (p = 0.5, T → 0), shown as a function of increasing hopping strength t hop (J = 1), with t hop = 0 being the SYK-NFL limit.As hopping increases, entanglement decrease from the a rather large value (≈ 0.35) and saturates to a value (≈ 0.15) close to that of the non-interacting model (see fig.3).
case (a)), and ψ-fermions, S (2) ψ (case(b)), are shown as a function of site fraction p in fig.6(a) for multiple values of temperature T .Ideally, at T = 0, in a pure state one expects S

FIG. 7 .
FIG. 7.Renyi entropy of diffusive Fermi liquid: (a) Behavior of Renyi entropy, S(2) , vs. subsystem size l in a model of noninteracting dots coupled to nearest neighbors in a ring of N dots = 10 dots (see fig.2(d)), for temperatures T = 0.03, 0.02, 0.01.The model is obtained by setting q = q hop = 1 and intra as well as inter-dot hopping strengths t chain = J = 1 in eqn.(55).(b) Renyi entropy vs. l for the non-interacting model when t chain = 1, J = 0.In both (a) and (b), large-N predictions (points joined by lines) are in agreement with the results obtained from exact numerics (color matched circles) performed using correlation matrices.(c) Extrapolated zero-temperature Renyi entropy as a function of subsystem size l for a ring of SYK dots, (fig.2(d)),connected by random hoppings.We set the interaction strength J = 1 within each SYK dot as well as the strength of hopping between dots t chain = 1.The number of dots in the chain, N dots , is varied from 20 − 100 dots (well beyond the scope of any exact numerics), to verify that the growth of entanglement for small subsystem sizes becomes independent of N dots for N dots > 20.(d) A close up, showing the growth of S(2) near small subsystem sizes, i.e. l ≤ 20, in the interacting SYK-chain (represented by squares) and the non-interacting chain (represented by pentagon) discussed in (a), for a large value of chain size N dots = 50.Both curves can be fitted (lines) with the analytical function S(2) ∼ log(1/(l −2 + l −2 0 ) 1/2 ), where l0 is an emergent "mean-free path" length scale which takes a value between ∼ 4 − 5 in this case.