Entanglement Entropy: Non-Gaussian States and Strong Coupling

In this work we provide a method to study the entanglement entropy for non-Gaussian states that minimize the energy functional of interacting quantum field theories at arbitrary coupling. To this end, we use a systematic method to build non-Gaussian variational trial wavefunctionals. The calculability \emph{bonanza} shown by these variational \emph{ansatze} allows us to compute the entanglement entropy using the prescription for the ground state of free theories. In free theories, the entanglement entropy is determined by the two-point correlation functions. For the interacting case, we show that these two-point correlators can be replaced by their nonperturbatively corrected counterparts. Upon giving some general formulae for general interacting models we calculate the entanglement entropy of half space and compact regions for the $\phi^4$ scalar field theory in 2D. Finally, we analyze the r\^ole played by higher order correlators in our results and show that strong subadditivity is satisfied.


Introduction
Quantum entanglement is a key concept that distinguishes quantum from classical physics.
It amounts to a class of nonlocal correlations between subsystems that are not present in a classical system. Remarkably, a well known measure of entanglement in quantum information theory, entanglement entropy, has found numerous applications in fields such as condensed matter physics, high energy theory and gravitational physics (see [1] and references therein).
Despite a huge amount of work in this direction, the characterization of entanglement entropy in strongly interacting quantum field theories (QFT's) through explicit and systematic computations, has shown to be rather intractable to do with the exception of holographic theories. Indeed, there is a limited number of field theories for which the entanglement entropy can be exactly computed. These are free field theories and strongly coupled conformal field theories (CFT's) with holographic duals. Noteworthily, for the later ones, the entanglement entropy can be computed using the Ryu-Takayanagi formula, which amounts to be one of the central entries in the AdS/CFT dictionary [2,3].
In general terms, when an observer has only access to a subset of the complete set of observables associated with a quantum system, tracing out the degrees of freedom localized on the non-accessible region, results in a reduced density matrix that represents the knowledge about the state of the quantum system that the observer posses. Being more explicit, let us consider that, at zero temperature, the total quantum system is described by the pure ground state |Ψ . Then, the density matrix is given by ρ = |Ψ Ψ|. If one arbitrarily divides the total system into two subsystems A and B, the total Hilbert space is the direct product H = H A ⊗ H B of the two spaces corresponding to subsystems A and B. The expectation value of an operator O A acting non-trivially on A is given by where the trace Tr A is taken only over the degrees of freedom pertaining to the Hilbert space H A and the reduced density matrix ρ A is defined as by tracing out the degrees of freedom associated to the the Hilbert space H B . Therefore, for the observer having access only to the subsystem A, the physical description of the total system is given by the reduced density matrix ρ A . The entanglement entropy of the subsystem A, which measures the amount of missing information about the total system for this observer, is given by the von Neumann entropy of the reduced density matrix ρ A , e.g., In QFT, the physical content of the theory is fully determined by the knowledge of its n-point correlation functions. Realistic operational assumptions on the observables usually impose that an observer cannot have full access to field configurations along the entire spacetime in which the QFT is defined. This implies an incomplete knowledge about correlations between points pertaining to the accessible region A and those lying outside it. Entanglement entropy quantifies the amount of information on these correlations that is loss for the observer. Despite an overwhelming body of work in this context, exact results are known only for the free QFT's and those are limited. A central feature for theories in d spatial dimensions, is that entanglement entropy is a quantity that depends on a short distance cut-off regulator as where |∂A| amounts to the area of the boundary of region A. The leading divergent term is known as the area law for the entanglement entropy. The area law is due to the large number of high energy modes that induce entanglement across the boundary ∂A of the accessible region. In this sense, while for any QFT there is an infinite number of degrees of freedom per unit volume, it is not possible to devise a realistic procedure to allow an observer in A to resolve infinitely small distances, so a sharp distinction between inside and outside of A is solved in terms of a UV-regulator. The cut-off can be viewed as a coarse graining parameter that represents to which extent the observer distinguishes the region A from the rest of the system. As a result, the entanglement entropy is necessarily cutoff dependent and sensitive to the UV structure of the quantum fluctuations. In case one wishes to use the entanglement entropy to analyze UV cutoff independent properties of the theory, then one needs to explicitly establish the sub-leading corrections in (4).
In addition, for any quantum system under consideration, the most important property of the entanglement entropy is known as strong subadditivity, This imposes that the entanglement entropy must be a concave function as the geometric parameters defining A are changed. It is the strongest condition one may set for the von-Neumann entropy. For the case of two dimensional field theories, where any region A amounts to a connected compact interval of length R and S A can be written as a function S(R), it has been shown that the converse is also true, e.g., the concavity of S(R) implies strong subadditivity for intervals [4].
It is worth to note that most of the results commented above have been explicitly computed only in the special case of free field theories. In these theories, it is widely known that their ground states are represented by Gaussian wavefunctionals and those are completely determined by the two-point functions of the theory. For the case of interacting QFT's, it is assumed that their ground states might be described in terms of non-Gaussian wavefunctionals and thus it is expected that the entanglement entropy should also depend on higher n-point functions. The question is to determine which of these correlations contribute and their weight on the structure of the entanglement entropy in specific models. In this context, the differences between free and interacting QFT's theories are not evident in terms of the entanglement entropy and have not yet been elucidated completely. As a consequence, explicit results on entanglement entropy for the case of strongly interacting field theories are scarce and much of our knowledge comes from holographic results [5]. Therefore, it is natural to develop new tools to investigate how higher order correlation functions give rise to the structure of the entanglement entropy in interacting field theories.
In [6], the half space ground state entanglement entropy 3 of the interacting λφ 4 and gφ 3 scalar field theories at weak coupling regime, was studied using the replica trick and position space Green's functions. The author showed that a consistent renormalization can be performed, providing finite contributions to the entanglement entropy at one loop. In [7], in a quantum mechanical setting consisting on quartic perturbative perturbations on the harmonic oscillator free case, the replica trick was used to check that to first order in pertur- 3 That is to say, the entanglement entropy of a region resulting from tracing out the half of the system.
More explicitly, regions A and B represent half spaces and their dividing boundary amounts to a flat space bation theory, the entanglement entropy can be computed by means of the perturbatively corrected version of the two point correlation functions of the system.
In [8], authors used a variational principle to determine a non-perturbative approximation to the half space ground state entanglement entropy of the λφ 4 theory at arbitrary coupling. The variational trial states used in this study were Gaussian wavefunctionals, for which the entanglement entropy can be exactly computed. Despite the Gaussian variational approximation at large values of the coupling is well defined, the approximation is only accurate up to one loop computations [6] and large N approximations in theories such as the self interacting O(N ) vector model.
To shed some light on this problem, in this work we present a non perturbative variational approach to compute the entanglement entropy of the self interacting λφ 4 scalar theory for arbitrary coupling. To this end, we use a class of non-Gaussian variational trial wavefunctionals that are systematically and non-perturbatively built through nonlinear canonical transformations [9][10][11].
In this respect, we briefly comment on some aspects that stem from applying a variational method to QFT: generality, calculability and ultraviolet modes.
Regarding the generality of the trial state, this must be general enough in order to capture the most salient physical features of the ground state of the theory through the variation of its parameters. For this reason, a systematic method to achieve this is desirable but not enough.
That is to say, even possessing a general ansatz for the vacuum wavefunctional for a QFT, we are interested in efficiently evaluating expectation values of operators. The calculability problem refers to our limited knowledge in evaluating non-Gaussian path integrals which forces us to restrict the set of available trial wavefunctionals to Gaussian states. These states (that represent the exact ground state for the case of free QFTs) have provided a great amount of non-trivial results when applied to interacting field theories. Regrettably, their favorable calculability also constrains their applicability to settings in which the relevant nonperturbative physics of the system is dominated by a single condensate that changes the mass of quasiparticles. Finally, due to the interaction between the high and low momentum modes in an interacting QFT, it would be desirable for a general variational ansatz to include variational parameters that optimally integrate out the effects of high energy modes into the low energy physics.
As it will be shown, the generality and calculability bonanza shown by the class of variational trial states used in this work, will allow us to compute the entanglement entropy using the same prescription as in the Gaussian case, which strictly depends on two-point correlators. In our case, these 2-p functions will be replaced by their nonperturbatively corrected counterparts. Concretely, the method shows how nonperturbative quantum corrections on the entanglement entropy of an interacting QFT can be obtained in a systematic manner.
The paper has the following structure. In Section 2 we introduce the NLCTs and explain the main properties of the trial non-Gaussian wavefunctionals that we build upon their application on Gaussian states. In Section 3 we review the obtaining of the entanglement entropy of the half space and compact regions for Gaussian states. Additionally, we explain our method to obtain the entanglement entropy of the same subsystems for non-Gaussian states. Our results are independent of the theory. Then, in Section 4 we pick the (1 + 1)dimensional φ 4 and evaluate such quantities for the trial state that minimizes the energy functional. In addition, we study the strong subadditivity condition for the case of intervals.
Finally, we discuss our results and future directions in Section 5.

Non-Gaussian States through NLCT
In the context of variational methods in QFT, Gaussian states are trial states that exactly represent the ground state of free field theories. The variational method consists of minimizing the expectation value of the Hamiltonian with respect to a set of variational parameters describing the trial wavefunctional Ψ[φ]. In particular, a Gaussian wavefunctional is parameterized by a real-valued functionφ(x) and a real-valued symmetric kernel G(x, x ), where td is the spatial dimension and N = [det(2πG)] −1/4 enforces Ψ to have unit norm. In the second line we have used a more compact notation.
From the Schrödinger picture of field theory it follows immediately that and that G(x, y) is the two-point function (propagator) or Green function [8] G Following [9][10][11], extensive non-Gaussian trial wavefunctionals can be nonperturbatively built asΨ where Ψ ≡ Ψ[φ] is a normalized Gaussian state (wavefunctional) and U ≡ exp(B), with As a result, the calculation of any expectation value with Ψ reduces to the computation of a finite number of Gaussian expectation values. In addition, the exponential nature of U ensures an extensive volume dependence of observables such as the energy of the system.
Moreover, as U is unitary, the normalization of the state is preserved.
The operator B consists of a product of π's and φ's, given by with m ∈ N, p, q i d-dimensional momenta and p = d d p/(2π) d . We denote these operators symbolically as B ≡ π φ m . Here, s is a variational parameter that tracks the deviation of any observable from the Gaussian case. In order to ensure an efficient truncation of the operator expansion given by the nested commutator series appearing in Hadamard's lemma, the variational function h(p, q 1 , . . . , q m ) is introduced in the ansatz and it must be optimized upon energy minimization. This function is symmetric w.r.t. exchange of q i 's and truncation imposes it has to satisfy: A suitable way of accomplishing (12) is taking where g(p, q 1 , . . . , q m ) is scalar function to be determined by the energy minimization and we have imposed that η(p) · κ(p) = 0, i.e., the domains of momenta where η and ζ are different from zero must be disjoint, up to sets of measure zero. Authors in [9,10] provided the useful ansatz for η and κ given by where ∆ 0 and ∆ 1 are variationally optimized, coupling dependent momentum cutoffs and the Heaviside step function. The coupling dependent momentum cutoffs variationally determine the most relevant region in momentum space to improve the estimation of the ground state energy. This fact is relevant in strongly-coupled theories in which the Gaussian quasi-particle picture is no longer valid.
The above constraints ensure that the commutator series terminates after the first nontrivial term. Namely, the action of U on the canonical field operators φ(p) and π(p) is given byφ where the quantities with a bar are defined as the nonlinear field functionals, A central consequence of U being unitary is that the canonical commutation relations (CCR) still hold under the nonlinear transformed fields (15) and (16) giving, For this reason, the transformations above are known as nonlinear canonical transformations (NLCT). As commented above, the consequence of Eq. (15) is that the expectation value of an arbitrary operator O(π, φ) w.r.t. Ψ, reduces to a Gaussian expectation value for the Throughout this paper we will consider two types of expectation values, which we will denote When considering QFT, n-point correlation functions of the form φ(p 1 ) · · · φ(p n ) Ψ are of particular interest. Using (15) and (16), one obtains . . .
This is a remarkable property of NLCT of the form π φ m , as non-Gaussian corrections to Gaussian correlation functions can be obtained in terms of a finite number of Gaussian expectation values. In particular, the terms proportional to s j in the non-Gaussian npoint correlation function correspond to (n + m(j − 1))-point Gaussian correlators, where j = 0, . . . , n.
Finally, we analyze the effect of the transformation U = exp(B) on wavefunctionals representing the probability amplitude for concrete field configurations. When applied to an interacting field theory such as the λφ 4 theory, the half-mean width of a variational Gaussian functional such as, is (k 2 + µ 2 ) −1/4 , where µ is the variational mass given by the gap equation, where I N (µ 2 ) = 1 2 k (k 2 + µ 2 ) N − 1 2 and m and λ are the bare mass and the bare coupling respectively. As the variational mass µ increases with interactions, then the nonclassical configurations accounted by the wavefunctional are much strongly suppressed than in the free case [10].
To illustrate the action of U on these Gaussian wavefunctionals, we choose the transformation B = π φ 2 for clarity. Noting that [10] whereφ(k) corresponds to (16) for m = 2, we obtaiñ where ellipses stand for the expansion of the exponential.

Entanglement Entropy of Gaussian States
For the sake of clarity in exposition, we describe here the two most general methods to calculate the entanglement entropy of Gaussian states, the replica trick and the real-time approach (see [1,4] for excellent reviews).

Half space
When the region of interest amounts to the half-space, the entanglement entropy of Gaussian states can be easily computed through the replica trick [12] S As said, being A the half-space, it can be shown that with Z δ being the partition function of a massive free scalar field of mass m living on an n-sheeted conical Riemann surface with a deficit angle δ = 2π(1 − n). When A is the half plane, d dn = −2π d dδ and (25) reads The result for the class of the Gaussian wavefunctionals under consideration can be symbolically written as [12,13] S where the constant accounts for the normalization of the Gaussian wavefunctional (6). Noting that one may write

Compact arbitrary regions
We are going to study the entanglement entropy of a compact region A for Gaussian states.
While this can done using the replica trick as well [14][15][16], here we will use the real time approach which was the first method used to compute entanglement entropy in free theories [17]. In this method, which has been mainly applied to numerical calculations in the lattice, the idea is to obtain the reduced density matrix ρ A in terms of the two point correlators restricted to the region A. As a quantum field theory can be completely specified in terms of its correlation functions, Eq. (1) implies that the information encoded by all the correlators inside A fully determines the density matrix ρ A . For the case of Gaussian states this is very simple, as the Wick theorem imposes that the only nontrivial correlators are the two point correlation functions. This was exploited in [18,19] to derive explicit expressions of ρ A in free boson and fermion discrete systems in terms of correlators.
Let us briefly review this approach for the case of (discretized) free boson theories [4].
The two point correlators of the field variables inside A are given by φ m φ n ≡ X mn , π m π n ≡ P mn , where X mn and P mn are real Hermitian and positive matrices. It can be shown that the entanglement entropy in region A is given in terms of the (positive) eigenvalues of C = √ XP as S [ρ A ] = Tr (C + 1/2) log(C + 1/2) − (C − 1/2) log(C − 1/2) .
For free quadratic bosonic Hamiltonians given by the ground state correlators are The matrix C = √ XP is constant, C = 1/2, when region A amounts to the total system, and thus, the entanglement entropy in Eq. 33 vanishes. In general, it does not vanish for an arbitrary compact region A specified by i = 1, · · · R where R < N (N , the total number of sites in the system) as far as is not necessarily 1/2 in general.
To give a simple example 5 , let us specify the Hamiltonian and correlators that must be used in a one dimensional lattice to calculate numerically the entanglement entropy for a real massive scalar: where a is the lattice spacing. In this setting, the Fourier modes of the field φ(p) are periodic functions of momentum restricted to the first Brillouin zone −π/a ≤ p ≤ π/a with the propagator [20] given by 5 Again, we emphasize that this prescription can be trivially generalized to any dimension.
With this, the correlators (35) can be written as φ n φ m = π/a −π/a dp 2π e i p a (m−n) 1 2 G(p) π n π m = π/a −π/a dp 2π e i p a (m−n) G(p) 2 . (39) Further details of this method can be found in [21].

Entanglement Entropy of Non-Gaussian States through NLCT
Having discussed the Gaussian case, now we focus in NLCT wavefunctionals which, as discussed above, take the effective Gaussian form Noting that the reduced density matrix is Gaussian w.r.t. the nonlinear canonically transformed fields, i.e., in order to compute the entanglement entropy of a NLCT non-Gaussian state we apply, and this is the central hypothesis of this work, the techniques for Gaussian states being discussed above. That is to say, as for Gaussian states, the entanglement entropy is fully determined through the two-point correlation functions, our proposal is to use the computational structure of the Gaussian free case but with the two-point correlators replaced by the nonperturbatively-corrected counterparts yielded by a nonlinear canonical transformation.

Half space
Regarding the entanglement entropy of the half-space for a Gaussian wavefunctional, Eq.
(30), a precise statement on the hypothesis presented above is that the entanglement entropy of non-Gaussian wavefunctional reads

Then, plugging these correlators into the expression for S[ρ A [Ψ]] (42), we have
where we have assumed that |s|< 1. As it will be made clear later, this assumption encompasses regimes of the physical coupling ranging from small to very large values.
For the case m = 2, the transformation in momentum space is given by: Here, the h function is given in (13).
Again, we emphasize that for an stateΨ, this expression is valid for any d-dimensional interacting theory. Such integrals must be evaluated when considering the optimal parameters that minimize the energy functional, which are precisely the objects that carry the information about one particular theory (see [10]).

Compact arbitrary regions
In terms of the real time approach, the entanglement entropy of Gaussian states for an area A is given by (33). As commented above, for a NLCT of the form π φ m , non-Gaussian corrections to Gaussian correlation functions can be obtained in terms of a finite number of Gaussian expectation values. In particular, as we focus on 2-point functions, the terms proportional to s j in the non-Gaussian 2-point function correspond to (2 + m(j − 1))-point Gaussian correlators, where j = 0, 1, 2. In this sense, that is how our method includes the effects of higher order correlations in the computation of entanglement entropy.
Then, based on the property (24) forΨ[φ], we can apply the entanglement entropy expression (33) for the correlators associated to the non-Gaussian states. That is to say, our proposal amounts to where the tilded correlators are These non-Gaussian 2-point correlators are obtained from the nonlinear transformation (16). For m = 2, they are given bỹ where u i is a unitary vector that locates the site of the lattice and X ij and P ij written in As in the half space case, we have obtained these expressions in full generality without specifying what is the interacting theory nor any dimensionality.
In the following section we are going to apply these formulas to a particular theory, which translates to picking a state that minimizes the energy through a variational method.

Entanglement Entropy in φ Theory
In this section we calculate the entanglement entropy of non-Gaussian states generated by NLCT that minimize the energy functional of the (1 + 1)-dimensional λφ 4 theory, which has the following Hamiltonian density: In d = 1, this theory does not exhibit any issue when renormalization is considered, and thus the method and our results lie on a solid ground. For higher dimensions, however, things are more subtle. Additional divergences arise from the non-Gaussian terms and they should be properly addressed [9,10].
Here we circumscribe to a treatment of the the theory in terms of a π φ 2 transformation.
The equations for the optimal values of the variational parameters s, G(p) and h(p, q 1 , q 2 ) are obtained, for a fixed φ c = φ(x) Ψ , by deriving H Ψ w.r.t. them and then equating to zero (details can be found in [9,10]). This yields a set of nonlinearly coupled equations that greatly simplify for φ c ∼ 0.
Indeed, in this case, the kernel G(p) reduces to the Gaussian case (29) up to order O(φ 2 c ), with a variational mass µ given by the gap equation (22), the optimal s = −4 λ φ c and Finally, with the results, numerical optimization is used in order to find the optimal values of the coupling-dependent cutoffs ∆ 0 and ∆ 1 .
Here, we carry out an optimization in which we have fixed φ c = 10 −2 . This implies that our formulas are valid for λ ≤ 100 in such a way we keep |s|< 1.

Half space
For the non-Gaussian transformation (11) with m = 2, the energy minimization requires the optimal values of the variational parameters to be ∆ 0 = 0.031 and ∆ 1 = 0.97. When fixing such values, we numerically evaluate the finite part of S[ρ A [Ψ]]/|∂A|.
In Figure 1 we plot ∆s N G , with where S G,div denotes the divergent terms of the Gaussian contribution, namely the divergences of the first term in (44) (see [6,22]). That is to say, we plot the finite part of the entanglement entropy per unit area of non-Gaussian states as a function of the coupling λ for states generated by the m = 2 NLCT. For comparison, we also plot this quantity for the Gaussian state. For a more detailed discussion of our results, we refer to Section 5.

Intervals
Connected compact regions in d = 1 are intervals determined, up to a translation, by their length R. For these intervals the entropy can be written as a function of a real variable S(R). We now calculate S(R) for non-Gaussian wavefunctionals generated by the NLCT transformations (11).
Upon picking the same optimization scheme as in the previous section (we will consider m = 2), we find the variational parameters ∆ 0 and ∆ 1 that minimize the energy functional for a fixed λ. Then we use these values to build the non-Gaussian state and thus evaluate the correlators in the entropy formula (46).
In Fig. 2 we plot S(R) for various couplings λ. These curves can be studied through a nonlinear regression analysis. We assume that S(R) is given by the nonlinear function   where c i and α are the fitting parameters. In Table 1 we show the value of c 2 and the exponent α that we obtain through nonlinear regression for different couplings. Let us note that for λ ≥ 1 6 the result α 1 and the positivity of c 2 give evidence of a well defined notion of mean entropy of the system [21]. 6 Let us remark that for λ = 0.1 the exponent does not satisfy this pattern. It would be interesting to explore this case when other NLCT that do not break the symmetric phase (m = 3) or other perturbative contributions are considered.
Finally, we want to check whether the entanglement entropy for the trial non-Gaussian states satisfies some entropic properties. In particular we will study the strong subadditivity property is fulfilled. Precisely, in 1 + 1 dimensions, an equivalence between strong subadditivity and concavity of the entropy can be established [21], Then, when calculating this quantity we observe that the entanglement entropy clearly satisfies such condition for any various couplings that differ in various orders of magnitude. More in detail, in Figure 3 we plot the second derivative of the entropy for the interval range that we have been studied and show that it is strictly negative.

Discussion and Outlook
In this work we have provided a method to study the entanglement entropy of non-Gaussian states generated by NLCT that minimize the energy functional of an interacting theory formulated at any dimension. Let us first discuss the entanglement entropy for the half space and the dependence on the coupling. In Figure 1 our variational approximation to the entanglement entropy shows that, in the perturbative regime, it is monotonically decreasing with respect to the coupling. This result is consistent for a coupling close to zero, where the Gaussian ansatz represents an accurate approximation to the ground state of the system. Namely, in this regime, the main effect of interactions is to increase the variational mass µ through the gap equation. This implies that field fluctuations contributing to entanglement are suppressed much stronger than in the free case. 7 Nevertheless, one would think that this cannot be the behavior of entanglement along all the coupling regimes, thus expecting that entanglement increases as the coupling constant also grows. This is exactly the behavior shown in Figures   1 and 2. Namely, results indicate that there are values of the physical interaction for which the entropy turns out to be monotonically increasing with the coupling. We interpret such turning point as a value at which the higher order correlators involved in the nonlinear canonical transformation under consideration become relevant, as the entanglement entropy for Gaussian states only captures the features of the connected 2-point correlators. We strongly believe that studying this value of the coupling could be of great physical interest since it corresponds to a local minimum of entanglement entropy. This issue must be explored in more detail in subsequent works.
When studying the entanglement entropy for intervals as a function of their size R, we also note its decreasing for Gaussian states as we increase λ. For the non-Gaussian states we have shown that strong subadditivity is satisfied through the concavity of the function (see Fig. 3). On the other hand, fittings of the curves of Fig. 2 suggest that the leading UV divergence remains proportional to |∂A|. Interestingly, our results also show an extensive contribution ∝ R. Thus the limit lim R→∞ S(R)/R ≥ 0 is well defined and there exists a well defined notion of mean entropy in the system. In other words, for big enough sets the entropy is approximately extensive [21]. This kind of contribution typically appears when a number (scaling with volume) of IR modes turn out to be highly entangled with the outside of the considered region and thus contribute to the entropy. By construction, an optimized nonlinear transformation acts by shifting some of the long distance IR modes of φ by a nonlinear polynomial functional of short distance UV modes. This shifting is strongly modulated by the strength of interactions and our results indicate that the transformed IR modes are promoted to be highly entangled with the outside.
Thus, in this work we have given evidence that the entanglement entropy for non-Gaussian states satisfies some required entropic properties. Still, several questions remain to be understood. For example, we would like to study to what extent the entanglement entropy is affected depending on the parameter m that determines the NLCT transformation. Additionally, it would be interesting to increase the dimensionality of the theory under consideration and check the reliability of our method through other entropic properties.
This method also allows for other further applications. For example, it is worth to mention that the approach of NLCT can be applied to fermionic field theories [23]. Despite the method can be formulated in any dimension and/or type of fermion, in [23] we only considered the case of two-dimensional Dirac fermions. This results specially well suited to analyze the Gross-Neveu model (GN) [24]. In view of the results of this paper, it would be interesting to investigate the behavior of entanglement entropy in a model possessing asymptotic freedom, chiral symmetry breaking and dynamical mass generation using techniques that go beyond the Gaussian analysis [25] such as those used in this work.
Finally, in light of our results, it would also be interesting to benchmark our method and its validity by calculating other quantities for Gaussian states and checking their properties (e.g., complexity [26,27] or quantum relative entropy [28]).