Renyi Entropy of Interacting Thermal Bosons in Large $N$ Approximation

Using a Wigner function based approach, we study the Renyi entropy of a subsystem $A$ of a system of Bosons interacting with a local repulsive potential. The full system is assumed to be in thermal equilibrium at a temperature $T$ and density $\rho$. For a ${\cal U}(N)$ symmetric model, we show that the Renyi entropy of the system in the large $N$ limit can be understood in terms of an effective non-interacting system with a spatially varying mean field potential, which has to be determined self consistently. The Renyi entropy is the sum of two terms: (a) Renyi entropy of this effective system and (b) the difference in thermal free energy between the effective system and the original translation invariant system, scaled by $T$. We determine the self consistent equation for this effective potential within a saddle point approximation. We use this formalism to look at one and two dimensional Bose gases on a lattice. In both cases, the potential profile is that of a square well, taking one value in the subsystem $A$ and a different value outside it. The potential varies in space near the boundary of the subsystem $A$ on the scale of density-density correlation length. The effect of interaction on the entanglement entropy density is determined by the ratio of the potential barrier to the temperature and peaks at an intermediate temperature, while the high and low temperature regimes are dominated by the non-interacting answer.


I. INTRODUCTION
Quantum entanglement and its various measures like entanglement entropies have become an important theoretical tool in the study of quantum many body systems. Starting from their domain of origin in the area of quantum information [1][2][3] , where entanglement is considered as a resource for quantum computation 4 , these measures have played an important role in the study of systems relevant to a wide variety of fields, from condensed matter theory to quantum gravity 5 .
In condensed matter systems, entanglement entropy has been used to detect phase transitions [6][7][8][9][10] and classify nontrivial topology of ground state of quantum systems [11][12][13][14] . Entanglement entropy 15 and entanglement spread 16,17 has been widely used to understand thermalization of complex many body systems. It has also been used to study many body localization [18][19][20] , where strongly disordered interacting quantum systems do not thermalize at long times. Entanglement entropy of many body system has recently been measured in experiments with ultracold atomic systems 21 .
While there is a large body of literature about entanglement entropy of Fermionic 22,23 and spin systems 24 , relatively less attention has been paid to microscopic calculations of entanglement entropy of interacting Bosons. This is mainly due to the fact that the larger Hilbert space for Bosons do not allow numerical evaluation of entanglement entropy, unless the system is in the Tonk's gas limit 25,26 in one dimension, or deep in the Mott insulator regime 27,28 , where the Hilbert space dimension can be manageable. The analytical approach comes from field theoretic considerations, where one either uses a replica based construction with a complicated manifold 29 , or takes advantage of simplifications from conformal field theories (CFT) if one is interested in a critical point 30,31 . While the CFT based approaches are most effective in one spatial dimensions, recent work has extended these considerations to more than one dimensions 32 for critical Bosons.
In a recent work 33 we have developed a new field theoretic method, to calculate the Renyi entanglement entropy for a generic (interacting) open or closed many body Bosonic system undergoing arbitrary non-equilibrium dynamics 34 . This method calculates the Wigner characteristic function 35 of a density matrix as the Schwinger-Keldysh (SK) partition function 36 of the system in the presence of sources turned on only at the time of measurement. For reduced density matrix of a subsystem, the sources need to be turned on only for degrees of freedom residing in the subsystem. The key new development is that unlike the well known replica based field theoretic approach to calculating entangelement entropy 29 , here one does not need to use complicated manifolds and can work with the standard field theory and its correlation functions. As a result, one is not restricted to particular geometries for the subsystem.
In this paper, we use this technique to calculate Renyi entanglement entropy of a thermal system of lattice Bosons at a fixed density, interacting with a local Hubbard repulsion. We consider N species of Bosons and a U(N ) symmetric interaction between the species 37 . We show that in the large N limit, the Renyi entropy of a subsystem is composed of two terms: (a) The Renyi entropy of an effective non-interacting system which has an additional spatially varying potential and (b) The difference in re-scaled thermal free energy between this effective system and the original translation invariant system within mean field approximation. This spatially dependent effective potential occurs due to the presence of the interactions between the Bosons and has no analogy in the case of non-interacting Bosons. Developing a systematic large N expansion for the Renyi entropy of the interacting theory and interpreting it in terms of an effective non-interacting system with a self consistent spatially dependent potential is the key new result of this paper. We divide a system of Bosons into two mutually exclusive sub-systems A and B. We use field theoretic methods for calculating Wigner functions to construct a functional integral for the 2nd Renyi entropy of the subsystem A. We use saddle point approximation in the large N limit to derive a self-consistent equation for the effective potential imposed by inter-particle interactions between the Bosons. We solve this equations numerically for a one dimensional chain and a two dimensional square lattice with nearest neighbour hopping and hence calculate the Renyi entropy of these systems. For the effective potential profile, we find (i) the potential profile takes the form of a step potential, which is flat in the bulk of the subsystem A and the subsystem B with different values and varies sharply across the entanglement cut (ii) The potential in the subsystem B is just the thermal Hartree shift of the original system. (iii) The potential barrier ∆V increases with temperature and saturates at large temperatures. The effect of interactions on the entanglement entropy is controlled by ∆V /T , which shows a peak at a characteristic temperature. (iv) The potential profile varies near the boundary of the subsystem A on a lengthscale given by the density-density correlation length of the original thermal system. This variation would become important near a phase transition where the correlation length diverges.
Focussing on the Renyi entropy, we find that the entropy scales linearly with the size of the subsystem (as is expected for a thermal system), which allows us to define an entanglement entropy density of the subsystem. The entropy density increases monotonically with temperature and saturates to its non interacting value at high temperatures. The high and low temperature values are governed by the density of the system. The effects of the interaction parameter is largest at an intermediate temperature scale where the ratio of the potential barrier to the temperature is largest. We also define an excess entropy density due to tracing, which measures the additional entropy per site of the subsystem relative to the entropy per site of the thermal density matrix of the whole system. This excess entropy, which indicates the additional randomness generated by tracing over the subsystem B, monotonically decreases with temperature and goes to zero at large temperatures. We find these trends in both one dimensional and two dimensional systems with quantitative differences. We would like to note that a microscopic calculation of Renyi entropy in more than one dimension is less explored in the literature and our field theory formalism, which is agnostic about the dimensionality, is a new step in this direction.
We will now present a roadmap for navigating the paper. In Section II, we review the well known relations between the Wigner characteristic function (WCF) of a density matrix and the second Renyi entropy of a subsystem. In Section III, we show that the WCF of a reduced density matrix of a thermal system can be written as a partition function with a specific set of sources within an imaginary time field theory. We work out the case of a system of non-interacting Bosons as an example in this section. In Section IV, we define the U(N ) symmetric interacting model of N species of Bosons and review its equilibrium field theory in the large N limit. In Section 32, we focus on calculating the Renyi entropy of the U(N ) symmetric model in the large N limit. The key results and the interpretation in terms of an effective potential is developed in this section. In Section VI, we apply this formalism to study a one dimensional Bose gas, while in Section VII, we look at the case of a two dimensional Bose gas. We conclude in Section VIII with a summary and possible future directions to extend this formalism.

II. WIGNER FUNCTIONS AND RENYI ENTROPIES
The Wigner function of a density matrix 38 is the closest equivalent to a phase space distribution function for a quantum system. For a many body bosonic system, it is useful to define the Wigner Characteristics function (WCF), χ W of the density matrixρ 35 , where c † r is the Bosonic creation operator for the spatial coordinate r, andD is the displacement operator which creates coherent states by acting on the ground state. Note that one can use any complete single particle basis to define the Wigner characteristics; we choose the position basis since it is the natural choice for calculation of entanglement entropies. The Wigner function is then defined by where we consider a system with Ω lattice sites. One can follow the same procedure, with the density matrix replaced by the operatorÔ, to obtain the Weyl symbol of the operator The expectation of the operatorÔ is then given by the phase space integral To calculate entanglement measures in a system, one considers a subsystem A consisting of Ω A sites. The reduced density matrix of the sub-system A, i.e.ρ A (t) = T r Bρ (t) is obtained by integrating out the degrees of freedom residing outside A (i.e. in B = A c ). The second Renyi entropy of the subsystem A is then given by The Wigner characteristic function of the reduced density matrix can be obtained from the full density matrix in the following way, where the lattice index j runs only over the Ω A sites in the subsystem A, whereas the trace is over all degrees of freedom. The trace over the degrees of freedom in B creates the reduced density matrix, while the displacement operator with support only in A leads to the Wigner characteristic function of the reduced density matrix. Thus the WCF of the reduced density matrix is obtained from Eq. 1 by restricting the support of the displacement operator to the Hilbert space of the subsystem A. For the sake of simplicity and clarity of the notation, throughout this paper, we denote the lattice sites belonging to the sub-system, A by the indices j, j etc, i.e. j, j ∈ A can run only over the Ω A lattice sites, whereas we will use r, r (or x, y) to indicate co-ordinates that run over the full system with Ω sites.
Finally, using Eq. 5 together with Eq. 4, we can write the second Renyi entropy of the subsystem in terms of the WCF of the reduced density matrix as We note that these results are quite general and are applicable to both non-interacting and interacting systems in and out of thermal equilibrium. In Ref. 33, a general Keldysh theory based formalism for calculating Renyi entropies of systems out of equilibrium has been formulated. However, in this paper, we focus on Renyi entropy of systems in thermal equilibrium and provide a detailed blueprint for calculating it within a field theory formalism. In the next section, we relate the Wigner characteristic of the reduced density matrix of a thermal system to a partition function with appropriate source configurations and apply it to a system of noninteracting Bosons.

III. RENYI ENTROPY OF THERMAL SYSTEMS
For a system in thermal equilibrium, the density matrix in the grand canonical ensemble is given byρ = 1 Z e −β(Ĥ−µN ) , whereĤ is the Hamiltonian of the system, β = 1/T is the inverse temperature,N is the total number operator and µ is the chemical potential. The partition function is given by Z = Tr e −β(Ĥ−µN ) , which ensures Trρ = 1. The partition function can be expanded in a standard functional integral where A is imaginary time action andĤ[c † r , c r ] is the normal ordered Hamiltonian of the system, and φ r (τ ) are the Bosonic fields.
The WCF of the reduced density matrix of the subsystem A is given by where β 0 is an arbitrary number between 0 and β, and we have used the cyclic properties of the trace to write the last line. The numerator can be expanded in a functional integral The insertion of the displacement operator at the imaginary time τ 0 = β 0 is equivalent to turning on a source 39 −iγ j coupled to the field φ * j (β 0 ) and iγ * j coupled to φ j (β 0 ). Note that the sources are turned on only for the fields in the subsystem A. They are turned on at the time β 0 and turned off immediately thereafter, creating a delta function profile in imaginary time. This suggests that the WCF χ A W ({γ j }) can be related to the equilibrium partition function 40 Z in presence of linear delta function sources, J j (τ ), turned on only at τ = β 0 for the modes j ∈ A, i.e.
where the denominator is calculated without any additional sources. This identification matches earlier results presented in Ref. 33, and can be considered as a special case of that more general result defined on real time Keldysh contours. Finally we note that the choice of β 0 is completely arbitrary and hence χ A W ({γ j }) is independent of the value of the imaginary time when the source is turned on. This will play a very important role when we consider an interacting model of Bosons in the next section.
We first apply this technique of calculating WCF and Renyi entropy of Bosonic systems to the simplest case of a generic non-interacting system of Bosons, with Hamiltonian, H 0 = p p c † p c p (See Ref. 41 for an alternate approach). We assume that the single-particle eigenstate p has a wavefunction ψ p (r). We consider Bosons with density ρ = (1/Ω) p n B ( p − µ, T ), where n B (e, T ) = 1/[exp(e/T ) − 1], is the Bose distribution function and the chemical potential µ(T ) is tuned to obtain the correct density at all temperatures. In this case, the action A is quadratic and we get Here, G −1 0 (r, τ ; r , τ ) is the imaginary time inverse Green's function of the non-interacting system. Performing the Gaussian integration over the fields φ r (τ ) in Eq.10, we obtain where, We note that one should be careful in taking the equal time limit τ = τ = β 0 inM A from the time ordered Green's function, G 0 (j, j ; τ − τ ), symmetrically along the imaginary time axis; i.e.
M A (j, j ) = 1 2 lim η→0 + G 0 (j, β 0 ; j , β 0 + η) The symmetric limit is due to the fact that the Wigner function is related to a symmetric rather than a (anti) normal ordering of c and c † operators. This limit also matches with the more general answer derived in Ref. 33 using Keldysh field theory. Note that in this case the answer is independent of the choice of β 0 , as we had predicted before on very general grounds. The system is still time-translation invariant along the imaginary time axis, τ and we can obtain the imaginary time Green's functions by the standard Matsubara summation technique 40 . This yields, We now obtain S (2) by performing the Gaussian integration over the γ j variables in Eq. 7, The entanglement entropy scales linearly with the subsystem size Ω A , as expected for a thermal density matrix. The above expression is further simplified if we calculate the Renyi entropy of the full system i.e. j runs over all the lattice sites. In this case,M A is a Ω × Ω matrix with eigenvalues, m p = [1 + 2n B ( p − µ, T )]/2. Substituting this in Eq. 15, we get The above example of a simple case of non-interacting Bosonic systems gives a consistency check of the new field theoretic method proposed to calculate S (2)41 . In the subse-quent sections, we will use this technique to calculate Renyi entropy of a reduced density matrix in presence of interparticle interactions in the system.

IV. U(N ) MODEL AND LARGE-N FIELD THEORY
In the previous section, we have obtained exact expressions for entanglement entropy of a subsystem of non-interacting Bosons which are in thermal equilibrium. Our aim is to extend this formalism to the case of interacting Bosons using a large N formulation. To this end, we consider a Bose Hubbard model with N species of Bosons. The Bosons hop between nearest-neighbour sites on a lattice and interact with local repulsive inter-species interaction given by, r (17) where r, r denote the nearest-neighbour lattice sites, a, b are species indices which run from 1 to N , t is the hopping matrix element and U is the scale of the local repulsion. We note that the scaling of the interaction, U/N , ensures that the energy of the system is extensive in the number of species N . It can be easily seen that this Hamiltonian is invariant under a global unitary rotation between the species, c where U is an arbitrary N × N unitary matrix. We would study this U(N ) symmetric model in the limit of large N .
Although our final aim is to calculate entanglement entropy of subsystems in this model, in this section, we will focus on the thermodynamics of the model. We will now calculate the partition function and hence free energy of this interacting Bosonic system in thermal equilibrium, using a large N approximation 37 . As we will find in the next section, this will form a part of the calculation of Renyi entropy of the subsystems. This will also allow us to discuss the large N approximation and resulting saddle points in a known and simple example 37 , while setting up notation for the large N calculation of entanglement entropy, which we will take up in the next section.
For the standard thermal field theory, the action A of this U(N ) symmetric model is given by Here, µ = µ(T ) is a chemical potential which ensures the correct number density per species ρ in the system. We will vary µ so that ρ is kept fixed at all temperatures. To calculate the partition function of this system in the large N limit, we perform the standard Hubbard Stratonovich (HS) transformation on the quartic interaction term by introducing the auxiliary fields, λ r (τ ), which couples to the bi-linears of the fields, r (τ ). Then the partition function reduces to where the inverse propagator 19 can be performed exactly to obtain, Note that each term in the exponent of Eq. 21 is multiplied by N . Now if we take the limit, N → ∞, we can evaluate the integral over λ in Eq. 21 in the saddle point approximation.
Choosing a static uniform saddle point solution −iλ r (τ ) = λ th , the saddle point equation can be written as where the dispersion for a hypercubic lattice in d dimension In addition, we have to solve the number equation It is then easy to show that the saddle point solution is λ th = U ρ, which is the standard Hartree energy shift in the system.
Substituting this saddle point value of the auxiliary variable λ th in Eq. 21 we obtain the free energy within the large N approximation, We note that for stability of the system, the effective chemical potential µ − λ th must be lower than the bottom of the band dispersion, −2td. Bose Einstein condensation in this system occurs when µ − λ th reaches the bottom of the band. In one dimension, this happens only at T = 0. In two dimensions the large N mean field theory for Schrodinger Bosons precludes the possibility of Bose Einstein condensation due to a weak logarithmic infrared divergence in the number equation. The physics of vortex binding and unbinding, which is missing from this description, leads to a BKT transition in this case. In three dimension, the large N theory predicts a finite temperature phase transition. In this paper we will confine ourselves to the "high temperature" or non-condensed phases of the Bosons. The case of the ordered phase will be taken up in a future work.
In the next section, we will modify the large N procedure described above for calculating the Wigner characteristic and hence Renyi entropies of a subsystem of interacting thermal Bosons.
In this section, we will devise a large N approximation to calculate the Renyi entropy of a subsystem of interacting Bosons. The standard large N approximation for the thermal field theory of U(N ) model has two parts: a Hubbard Stratanovich transformation with an auxiliary field coupled to a U(N ) invariant bilinear of the fields, and a saddle point evaluation of the integrals over the auxiliary fields, justified by the N → ∞ limit. Since the large N approximation for the partition function has already been shown in the last section, here we focus on the partition function in presence of the sources, i.e. the WCF χ A W ({γ j }). We will use the path integral representation of χ A W ({γ j }) shown in Eq. 10, with the action given by Eq. 18. Note that in this case each species of the Boson fields are coupled to sources. Following the steps in the previous section, we will first decouple the interaction terms in A by a HS transformation to get where we would like to emphasize once again that r, r runs over all lattice sites, while j runs only over the subsystem A. The partition function Z in the denominator has already been calculated in the previous section. The Gaussian functional integration over the fields φ r (τ ) yields where,Ĝ −1 λ is the interacting inverse propagator defined in Eq. 20. Once again, M A λ (j, j ) = G λ (j, β 0 ; j , β 0 ), with the equal time limit taken symmetrically as defined in Eq. 13. In this case, it is useful to substitute χ A W from Eq. 26 in Eq. 7 and perform the Gaussian integration over the source variables, γ (a) j . We then obtain the Renyi entanglement entropy, S (2) from, The exponent in the integral is multiplied by N , and in the limit of large N , we can evaluate the above integral over the two auxiliary fields using saddle point approximation. Note that the saddle point solutions will be different from that of the thermal field theory due to the extra term ∼ T r A ln{M A λ +M A ν }, which comes from integrating over the source variables γ (a) j of the WCF. This term only depends on the Green's functions in the sub-system A. Symmetry indicates that the saddle point profiles λ 0 r (τ ) = ν 0 r (τ ). Considering the saddle point condition, we get The first term in Eq. 28 is the standard answer one gets for the thermodynamics of the system. This term is independent of τ as well as independent of β 0 . The second term comes from integrating out the sources γ (a) j . This is an additional term that occurs when evaluating the Renyi entropy and depends explicitly on τ as well as on β 0 . However, as we had mentioned before, the choice of β 0 is arbitrary and answers should not depend on this variable. Hence, we integrate over the variable β 0 from 0 to β and divide by β. This approximation allows us to consider a spatially varying, but (imaginary) time independent saddle point profile λ 0 r , given by Note that we have restored the (imaginary) time translation invariance of the problem, and hence one can use Matsubara frequency sums to evaluate the integrals above. The saddle point is then equivalent to an effective non-interacting Hamiltonian with a spatially varying external potential V (r) = −iλ 0 r , which has to be determined self consistently. If the eigenval-ues and eigenfunctions of this effective single particle Hamiltonian are E n and ψ n (r) (H 0 is the non-interacting part of the original Hamiltonian given in Eq. 17), the self-consistency equation is given by Substituting this saddle point solution, V (r) in Eq. 27, we obtain the Renyi entanglement entropy in the large N limit, .
We note that the leading order approximation to the entanglement entropy simply scales with the number of Boson components N , and hence it is useful to define the entropy per species of the Bosons. From the above formula we see that the entanglement entropy of the system in large N approximation has two contributions: (a) the first line corresponds to the difference of the re-scaled thermodynamic free energy, βF of two systems: the original translation invariant system that we are interested in, and the "effective" system with a spatially dependent potential V (r), stemming from the entanglement cut in the system, which has to be self-consistently determined from Eq. 31 and (b) a second contribution (second line) related to the one-particle correlation function M within the subsystem A. This correlation function is calculated for the system with the "effective" Hamiltonian. Note that in case of a non-interacting system, there is no effective potential, and hence we simply get the second line as the full answer. Our formalism thus generalizes the non-interacting answers to the case of interacting system within the large N approximation.
In the next two sections, we will use this large N solution to calculate the Renyi entropy of a 1-d and 2-d interacting Bose gas.

VI. RENYI ENTROPY OF 1D BOSE GAS
We consider a one dimensional chain of size L, where Bosons hop with a nearest neighbour tunneling amplitude t and interact locally with a repulsive interaction scale U . Considering N species of Bosons, the U(N ) symmetric Hamiltonian is given by where a, b indicate the species index of the Bosons. The Hamiltonian has periodic boundary conditions to ensure translational invariance. This leads to a dispersion p = −2t cos p, where p is the lattice momentum. We will set t = 1 and the lattice-spacing a = 1 to fix units of energy and length in all subsequent discussion in this paper. For 1D chain, the band energies vary between −2 and +2, giving a bandwidth of 4t = 4. We consider the system at a temperature T with a fixed density ρ, while we allow the chemical potential µ to vary with temperature. This ensures the same number density at all temperatures. We are interested in the Renyi entropy of the reduced density matrix of a contiguous subsystem of size L A . Note that due to the translation invariance of the original Hamiltonian, this subsystem can be placed anywhere in the system; we choose to place it in the left half of the system. In this section, we will work with L = 100, unless otherwise mentioned. For a lattice Boson system, the chemical potential lies below the bottom of the band at high temperatures. If the chemical potential reaches the bottom of the band, Bose Einstein condensation occurs in the system. It is well known that in 1 dimension, the low energy divergence of the single-particle density of states D( ) ∼ −1/2 , together with the ∼ −1 diver-gence of the Bose function, n B ( ) = [e β − 1] −1 , prevents the formation of a BEC in the system. This remains true for the interacting system within the large N approximation. Thus our formalism, which does not take into account possible Bose condensation, works at all temperatures for the 1 dimensional system.
The first task in calculating the Renyi entropy of a subsystem is to obtain the self-consistent profile of the effective potential V (x) from the saddle point equation Eq. 31. The selfconsistent profile obtained for three different temperatures, T = 7.5, T = 2.5 and T = 1 are plotted as a function of the lattice sites in Fig. 1(a). In this case, U = 0.5 and ρ = 0.5. From the figure, we see that the potential has a flat profile in the bulk of subsystem A with a value V A and a similar flat profile with a different value V B in the bulk of subsystem B. It changes near the entanglement cut at x = 50 and x = 0 very rapidly over a scale of a few lattice spacings. This is depicted schematically in Fig. 1(b) as a step potential barrier with a height ∆V (T ) = V B (T ) − V A (T ). Near the entanglement cut, we find that the potential in the subsystem A drops near the edge before rising sharply in the B subsystem. The lengthscale of variation on either side is almost equal to the correlation length ξ(T ), obtained from the exponential decay of the connected density-density correlation function in the thermal system. This is not surprising if we think of the entanglement cut as a localized source of disturbance. One would then expect the local density, and hence the effective potential, to adjust to this perturbation on a scale of ξ(T ). To see this clearly, we plot the potential profile at T = 0.1 in Fig. 1(c), where the variation of the potential occurs on a measurable length scale of ∼ 9. In the inset of Fig. 1(c), we plot ξ(T ) with the temperature of the system. We see that at T = 0.1, ξ(T ) is also ∼ 9. The key takeaway from this is that, while the entanglement calculation inevitably involves solving a translation non-invariant problem of the effective potential (even though the original system is translation invariant), one can get away with solving a much simpler problem of a step potential well in the subsystem A, as long as the correlation length ξ(T ) is much smaller than the subsystem size L A . However, the value of the barrier needs to be solved self-consistently. This can in principle be used to calculate entanglement entropies in large systems in higher dimensions. We note that this potential barrier is imposed by the effect of the inter-particle interactions in the system and vanishes in the limit U → 0.
The bulk value of the potential in the subsystem B, V B is independent of the temperature, as seen from Fig. 1(a). In fact, this value is simply the thermal value of the auxiliary field λ th = U ρ. To see this, we plot V B for systems with different values of U and ρ as a function of U ρ in Fig. 1(d), and obtain a straight line with unit slope. The potential barrier across the entanglement cut increases with temperature, saturating at very high temperatures. This is shown in Fig. 1(e) for three different set of parameters: U = 0.5, ρ = 0.5 (blue solid line), U = 1, ρ = 0.5 (purple dashed line) and U = 1, ρ = 1 (green dashed-dotted line). The high temperature saturation value increases with density of the system as well as with the interaction in the system.
Since we are considering a thermal system with a potential . ∆V (T )/T rises to a peak value at T = Tmax before falling as T −1 at high T . Tmax increases with increasing U and ρ. Since this dimensionless quantity ∆V (T )/T controls the effect of the inter-particle interaction in the thermal system, the effects of the latter in S (2) will also be the largest at T = Tmax. In all the plots, we set t = 1 and a = 1 as units of energy and length respectively.
barrier, the effects of this barrier, and hence of interaction in the system, on S (2) will be controlled by the dimensionless parameter ∆V /T . This is plotted in Fig. 1(f) as a function of temperature for three different set of parameters: systems with U = 0.5, ρ = 0.5 (blue solid line), U = 1, ρ = 0.5 (purple dashed line) and U = 1, ρ = 1 (green dashed-dotted line). The dimensionless parameter peaks at a typical value of temperature T max before going down as ∼ T −1 at large temperatures. T max increases with increasing U and ρ. We should expect the effects of interaction on S (2) to be largest around T max .
We now focus on the Renyi entanglement entropy of the subsystems of this interacting thermal system. The entanglement entropy scales linearly with the size of the subsystem, as expected for a thermal density matrix. To show this, in Fig 2(a), we plot the entanglement entropy of a subsystem with L A for a system of size L = 600 and U = 0.5, ρ = 0.5 for two different temperatures, T = 3 and T = 0.1. This allows us to define an entanglement entropy density (entan-glement entropy per site per species of the Bosons), s (2) and we will now focus on how this entanglement entropy density changes with temperature, density and interaction in the system. In Fig. 2(b), we plot the entanglement entropy density for a subsystem of size L A = 50 as a function of temperature. This is done for a system with L = 100 and U = 1 and two different values of density, ρ = 0.5 and ρ = 1.0. The entropy density increases with temperature and saturates at high temperatures. At large T , s (2) (triangles) closely agrees with the non-interacting answer (= log(coth[|µ|/2t])) plotted by red solid line in the inset of Fig. 2(b) for ρ = 0.5. This can be understood in the following way: the chemical potential |µ(T )| (solid line with squares) increases linearly with T as shown in the inset of Fig. 2(b) and at large T , ∆V (T ) imposed by the inter-particle interaction becomes negligible compared to µ(T ). Hence, at large T the effects of interaction becomes completely suppressed in s (2) . This is expected since ∆V /T goes to 0 in this limit, as seen from Fig. 1(e).
In order to understand the effect of interaction on entangle-

FIG. 2.
Renyi entropy of the reduced density matrixρ A consisting of the DOFs belonging to the sub-system A (0 ≤ x < LA) for 1d Bosons: (a) shows S (2) scales linearly with the sub-system length LA both at high temperature T = 3.0 (squares) as well as in low temperature T = 0.1 (circles). We use L = 600, U = 0.5 and ρ = 0.5 for this plot. The entanglement entropy density (entanglement entropy per site per species of the 1D Bosons) s (2) is plotted as a function of T in (b) for two different value of densities, ρ = 0.5 (triangles) and ρ = 1.0 (squares) for same value of U = 1.0. s (2) increases with T and saturates at high T , with the saturation value increasing with ρ. At large T , s (2) (triangles) closely agrees with the non-interacting answer plotted by red solid line in the inset for ρ = 0.5. This can be understood in the following way: the chemical potential µ(T ) (solid line with squares) increases linearly with T as shown in the inset of (b) and at large T , ∆V (T ) imposed by the interparticle interaction becomes negligible compared to µ(T ). Hence, at large T the effects of interaction becomes completely suppressed in s (2) . To see this clearly, we plot ∆s ( ment entropy density, we consider ∆s (2) = s (2) (U ) − s (2) (0) , the difference between the entropy density of an interacting system and that of a non-interacting system at the same temperature and density. In Fig 2(c), we plot ∆s (2) as a function of temperature for a system with ρ = 0.5 for three different values of U = 0.5, U = 1.5 and U = 2.5. We note that ∆s (2) increases at low temperature and reaches a peak before going down at large temperatures. The large temperature decay ∼ T −1 , reflecting similar decay in the dimensionless potential barrier ∆V /T , shown in Fig. 1(e). In fact, the behaviour of ∆s (2) is very similar to that of ∆V /T over the full temper-ature range considered here. The peak position and the peak value of ∆s (2) increases with increasing U/t. In Fig 2(c), we have also shown the peak positions of ∆V /T with color coded arrows for the same parameters. This shows that the peak of ∆s (2) closely tracks that of ∆V /T .
Every density matrix of a quantum system (or subsystem) can be interpreted in terms of an ensemble of pure quantum states drawn with a particular classical probability distribution. If one diagonalizes the density matrix, the eigenstates are the pure quantum states in question and the eigenvalues are the corresponding probabilities. It is obvious that entropy measures defined on density matrices are actually entropy measures of the corresponding classical probabilities associated with the density matrix. For example, the Renyi entropy of a density matrix S (2) = − ln Trρ 2 = − ln i p 2 i , where the eigenvalues p i correspond to the classical probability distribution.
For entanglement measures, when one constructs the reduced density matrix of a subsystem from a pure quantum state in the full system by tracing over degrees of freedom, the loss of the "quantum" information shows up as a classical probability distribution in the reduced description. Entanglement entropies are various entropy measures of this distribution and measures the classical randomness associated with the loss of the "quantum" information. However, when the state of the full system is described by a density matrix, as is the case for thermal systems at finite temperature, there is a classical randomness associated with the full system as well. In this case, we can define the additional randomness introduced by the loss of information (tracing), i.e. an entropy density of tracing by considering the difference between the entropy density of the subsystem, calculated from the reduced density matrix, and the entropy density of the full system, calculated from the full thermal density matrix, S (2) = s (2) (L A ) − s (2) (L). This is a measure of the additional randomness introduced by tracing over the subsystem B. Note that since the number of degrees of freedom is different in the subsystem and the full system, it is important to divide the entropy by the size of the corresponding system before subtracting them.
In Fig 2(d), we plot the additional entropy density S (2) of a subsystem of size L A = 50 for a system with L = 100, U = 1 and ρ = 0.5 as a function of temperature. We see that the additional entropy density is a monotonically decreasing function of temperature, going to 0 at very large temperatures. At low temperatures, the full system is in a almost pure quantum state and the entropy density of the full system s (2) (L) is very small. Hence, a large entropy density is gained by tracing out the degrees of freedom. On the other hand at very high temperatures, the system behaves like a bunch of classical free particles and hence tracing does not lead to generation of any additional randomness in the system.
In the next section, we will extend this discussion to the case of a 2D interacting Bose gas.

VII. RENYI ENTROPY OF 2D BOSE GAS
In this section we extend our formalism to a 2 dimensional system. We would like to note that a microscopic calculation of Renyi entropy in more than one dimension is less explored in the literature and our field theory formalism which is agnostic about the dimensionality of the system is a new step in this direction. We consider a Bose gas on a square lattice with nearest neighbour hopping t and a local Hubbard repulsion U . Once again, we consider N species of Bosons and a U(N ) symmetric model given by With a periodic boundary condition in both x and y direction, we can work in the lattice momentum basis with a dispersion p = −2t(cos p x +cos p y ), resulting in a bandwidth of 8t. The non-interacting system does not undergo a Bose Einstein condensation at any finite temperature due to an infrared logarithmic divergence in the number equation in 2-d , in accordance with the Mermin-Wagner theorem 40 . While the real system would show a Berezinskii Kosterlitz Thouless type transition at a finite temperature to a disordered phase from a phase with quasi-long range order due to proliferation of vortices, this is not captured within a large N theory, which does not account for the vortex excitations. We consider the large N approximation to the Renyi entropy of a L A × L A subsystem located around the center of a L × L lattice. We will present data for L = 24 and L A = 12 in this section, unless otherwise mentioned. We first consider the effective potential V (x, y) for this system.
The effective potential profile for a system with density ρ = 0.5 and interaction strength U = 0.5 is shown in Fig. 3. Fig. 3 (a) plots the effective potential profile for a high temperature T = 4.0. Similar to the one dimensional case, we find that V (x, y) is almost constant in the subsystem B and sticks to its thermal value U ρ. The potential forms an almost

U S E T d A E 0 5 9 i p H S 0 s i s D A O k p r 5 A s 1 R m D 2 m t c Z 6 V R 2 b V r t t z W K v E K U g V C r R H 5 t d w H O I 4 I F x h h q Q c O H a k 3 B Q J R T E j W X k Y S x I h P E M T M t C U o 4 B I N 5
2 H z 6 w z r Y w t P x T 6 c W X N 1 d 8 b K Q q k T A J P T + Z R 5 b K X i / 9 5 g 1 j 5 l 2 5 K e R Q r w v H i k B 8 z S 4 V W 3 o Q 1 p o J g x R J N E B Z U Z 7 X w F A m E l e 4 r L 8 F Z / v I q 6 T b q T r N + d d e s t q 6 L O k p w A q d Q A w c u o A W 3 0 I Y O Y E j g G V 7 h z X g y X o x 3 4 2 M x u m Y U O x X 4 A + P z B x t s l G 0 = < / l a t e x i t > ρ = 0.5 < l a t e x i t s h a 1 _ b a s e 6 4 = " Q o 5 o i t a F G H G t P N P + p / 0 z c 9 2 r P o A = " > A A A B + H i c b V D L S s N A F L 2 p r 1 o f j b p 0 M 1 i E u i l J K a i 7 o h u X F e 0 D 2 l g m 0 0 k 7 d P J g Z i L U k C 9 x 4 0 I R t 3 6 K O / / G S Z u F t h 6 4 c D j n X u 6 9 x 4 0 4 k 8 q y v o 3 C 2 v r G 5 l Z x u 7 S z u 7 d f N g 8 O O z K M B a F t E v J Q 9 F w s K W c B b S u m O O 1 F g m L f 5 b T r T q 8 z v / t I h W R h c K 9 m E X V 8 P A 6 Y x w h W W h q a 5 W R A M E d 3 6 U N S r Z + l p a F Z s W r W H G i V 2 D m p Q I 7 W 0 P w a j E I S + z R Q h G M p + 7 Y V K S f B Q j H C a V o a x J J G m E z x m P Y 1 D b B P p Z P M D 0 / R q V Z G y A u F r k C h u f p 7 I s G + l D P f 1 Z 0 + V h O 5 7 G X i f 1 4 / V t 6 F k 7 A g i h U N y G K R F 3 O k Q p S l g E Z square well inside the subsystem A with a constant value V A . The potentials match up at the boundary separating the subsystem A from the subsystem B on a scale of the connected density-density correlation length. Fig. 3(b) plots the effective potential profile for the same system at a lower temperature of T = 1.5. In this case, the profile varies around the entanglement cut on a larger lengthscale. In fact, for a subsystem size of 12 × 12, the potential profile within the subsystem A looks almost parabolic and leads to deep pockets at the corners of the entanglement cut. To understand this behaviour, we have plotted the thermal connected density-density correlation length ξ(T ) as a function of temperature in Fig. 3(c). We see that the correlation length is ∼ 5 at T = 1.5, and hence the subsystem size is already of the same scale as the correlation length at this temperature within our finite size calculation. This explains the large spatial variations of the effective potential profile in this case. We consider the potential at the central point of subsystem A to be V A (T ) and the potential at the edge of subsystem B to be V B and plot the potential barrier ∆V = V B − V A (T ) as a function of temperature in Fig. 3(d). The two graphs both correspond to a system with ρ = 0.5. The solid line corresponds to a system with U = 0.5 and the dash-dotted line corresponds to a system with U = 1. The potential barrier increases with temperature and saturates at high temperatures to a value which increases with the interaction strength. As mentioned before, the interaction effects are controlled by ∆V /T , which is shown in Fig. 3(e) as a function of temperature. It shows a peak at a characteristic temperature T max , which increases with the interaction strength. We expect the interaction effects on entanglement entropy to be largest around T max .
We plot the Renyi entanglement entropy of the subsystem with the area of the subsystem for two different temperatures T = 4 and T = 1.5 in Fig 4(a), which clearly shows the linear dependence of S (2) with the area of the subsystem, as expected for a thermal density matrix. The entropy per site, s (2) is plotted as a function of temperature for a system size of 24 × 24 and a subsystem size of 12 × 12 in Fig 4(b). The two plots correspond to densities of ρ = 0.5 and ρ = 0.25 respectively. The interaction strength in this case is U = 0.5. The entropy density increases with temperature and saturates at large temperature to values determined by the density of the system. To understand the effects of interaction, in Fig 4(c) we plot ∆s (2) = s (2) (U ) − s (2) (0), the difference between the entanglement entropy density of an interacting system and a non-interacting system, as a function of T for a system with density ρ = 0.5. The interaction effects peak at an intermediate temperature which closely follows T max where ∆V /T is largest. This is seen from the plots, where the location of T max is plotted as additional arrows. The plot also shows that the peak position increases with increasing U/t. Finally, in Fig 4(d), we plot the difference in entropy density of the subsystem from that of the full system S (2) = s (2) (L A )−s (2) (L) as a function of T for L A = 12 and L = 24. We see that this difference monotonically decreases with temperature.
We note that the finite size effects are more severe in two dimensions than in one dimension, restricting us to a regime of relatively high temperatures. There are two reasons for this: (i) the lengthscale we can access numerically is smaller in two dimension than in one dimension and (ii) the infrared divergence of the number equation, which prevents a Bose Einstein condensation, has a weak logarithmic dependence in two dimensions compared to an inverse square root divergence in one dimension. While this can be mitigated to some extent by going to larger system sizes, the general trend that finite size effects will be more severe in two dimensions will remain an intrinsic factor for this model.

VIII. CONCLUSION
In this paper, we have used a Wigner function based field theoretic approach to calculate the Renyi entropy of a subsystem of interacting Bosons in thermal equilibrium. Using a U(N ) symmetric model of lattice Bosons interacting via a local Hubbard repulsion, we derive a functional integral for the second Renyi entropy of a subsystem of these Bosons, when the full system is in thermal equilibrium at a fixed density.
Using a saddle point approximation in the large N limit, we show that the entanglement entropy can be calculated in terms of an effective system with an externally imposed potential, to be calculated self-consistently. Although the system is translationally invariant, the consideration of a subsystem breaks the translation invariance. Hence the effective potential is spatially varying. We derive the self-consistent equation for this potential and solve it numerically for 1 and 2 dimensional Bose gas. In both cases we find that the potential is flat in the bulk of subsystem A and subsystem B, with different values, thus creating a potential barrier across the entanglement cut. The potential barrier increases with temperature and saturates at high temperature. The potential varies near the boundary between subsystem A and B on a scale of the density density correlation length.
The entanglement entropy scales with the size of the sub-system (volume law) and hence one can define an entanglement entropy density in the subsystem. The effect of interaction on this entropy density is largest at an intermediate temperature where the ratio of the effective barrier to the temperature peaks. This peak temperature increases with increase in interaction strength in the system. There are two directions where the current formalism can be extended. One is to go to larger system sizes, specially for higher dimensional systems. The other is to consider fluctuations around the static saddle point that we have considered. We leave these issues for investigation in a future work.