Renyi entanglement entropies for the compactified massless boson with open boundary conditions

We investigate the Renyi entanglement entropies for the one-dimensional massless free boson compactified on a circle, which describes the low energy sector of several interacting many-body 1d systems (Luttinger Liquid). We focus on systems on a finite segment with open boundary conditions and possible inhomogeneities in the couplings. We provide expressions for the Renyi entropies of integer indices in terms of Fredholm determinant-like expressions. Within the homogeneous case, we reduce the problem to the solution of linear integral equations and the computation of Riemann Theta functions. We mainly focus on a single interval in the middle of the system, but results for generic bipartitions are given as well.


I. INTRODUCTION
The challenge of quantifying the entanglement content of extended systems is at the common intersection of several communities, ranging from high energy to quantum information and condensed matter [1][2][3], with a recent emphasis on out-of-equilibrium features (see Ref. [4] and references therein). It has been understood that different phases of manybody systems can be characterized in terms of their entanglement properties [5,6], which also cover an important role in the efficiency of tensor networks algorithms [8][9][10][11][12][13]. Furthermore, recent times have witnessed a direct engineering of entanglement measurements in cold atom experiments [14][15][16][17][18][19][20]. Given a quantum system initialized in a pure state |ψ and bipartited in two halves A ∪ B, arguably the most widespread measurement of the entanglement between A and B is given by the Von Neumann entropy whereρ B is the reduced density matrix of the subsystem B, obtained through a partial trace of the degrees of freedom of the subsystem A, i.e.ρ B ≡ Tr A |ψ ψ|. Together with the Von Neumann entanglement entropies, one can consider the more general Renyi entropies of index N , defined as The knowledge of all the moments of the reduced density matrix gives information about its whole spectrum [21][22][23]: the Renyi entropies are the central objects of this work. The Von Neumann entanglement entropy can be understood as the limit N → 1 of the Renyi entropies. Such a point of view is often used as a useful computational tool, commonly known as replica trick [5]: the Renyi entropies for integer indices N are computed, subsequently the result is analytically continued to real indices and the limit N → 1 is taken. For generic systems, the computation of Renyi entropies is a tremendous task on its own, not to speak of their analytical continuation. A great insight has been achieved at criticality: the system acquires scale invariance and the powerful methods of conformal field theory (CFT) [24,25] can be used to investigate the scaling part of the entanglement entropies. In the one dimensional world, CFT displays its full power and allows one to derive an array of simple analytical results. For example, for an infinite system in the ground state with a bipartition chosen to be an interval of length L and the rest, the following result holds true [5][6][7] S N L = c Above, c is the central charge of the CFT and corrections to the conformal prediction go to zero as L → ∞. The constant a N (but in principle N −dependent) is an ultraviolet cutoff determined by the microscopic features of the model (i.e. not fixed by the CFT). If the microscopic system lives on a lattice, then a N si proportional to the lattice spacing. Several many-body quantum systems can show critical features, in particular the compactified free boson (Luttinger Liquid) [26][27][28][29][30][31][32] is a ubiquitous low-energy description of models with a U (1) conserved charge, such as the total number of particles in cold atoms systems or the total magnetization in spin chains. The compactified boson is arguably one of the simplest CFT, having the central charge equal to one c = 1 [29,30].
While it remains true that the scaling part of the entanglement entropies is determined by the CFT, on the other hand the global conformal invariance alone does not uniquely fix the entanglement of arbitrary bipartitions, which in turn depends on the details of the underlying CFT. For example, in a system with periodic boundary conditions (PBC) and a bipartition made out of two disjoint intervals with the rest, the global conformal symmetry is not enough to determine the entanglement entropy. A lot of efforts have been directed towards this problem [5,, obtaining closed analytical expressions for the compactified boson [57][58][59].
Another situation where one cannot rely just on the global conformal symmetry is finite subsystems on an interval [0, L], with some boundary conditions. Apart from the per se interest, the study of these boundary CFTs [60] has many applications going from Kondo physics [61][62][63] to string theory. Indeed, the entanglement in Kondo problems has been already intensively studied (see e.g. Ref. [3] and references therein). In this paper, we start analyzing this problem from the simplest possible case, which are open boundary conditions (OBC).
Let us consider a bipartition A ∪ B where we choose B = [x 1 , x 2 ] with x 1 < x 2 : if one of the two extrema coincides with the boundaries of the system, i.e. x 1 = 0 or x 2 = L, conformal invariance unambiguously predicts the entanglement entropy [5,6] (albeit the O(L 0 ) contribution can acquire an extra term due to the non trivial boundaries [5,64], which is the boundary entropy of Affleck and Ludwig [65]). On the other hand, if 0 < x 1 < x 2 < L, global conformal symmetry is not sufficient anymore and a detailed analysis of the operator content of the CFT is required. In this case, closed analytical results are known only for free fermionic systems [66][67][68][69].
OBC are of primary importance to be explored, both for the experimental relevance and for their importance in numerical studies. Indeed, the advent of DMRG techniques [70] provided us a powerful numerical tool to study the one-dimensional world and these are best suited to work with OBC. However, we stress that while the entanglement entropies of bipartitions starting at the boundaries (i.e. choosing A ∪ B with A = [0,x]) are easily computable with DMRG methods, studying the case of an interval in the middle of the system in the scaling regime is much harder. In this perspective, alternative analytical and semi-analytical tools are surely helpful.
Going beyond the global conformal invariance of the theory is extremely hard, but dealing with the compactified boson we can take advantage of an important simplification: the model is free. Efficient techniques exist to deal with free (i.e. gaussian) systems [71], and even though extracting the scaling part of the entanglement in terms of simple analytical expressions can be hard, the result is straightforwardly numerically computed diagonalizing matrices which scale proportionally with the system size. Such a task can be easily performed on a laptop for quite large sizes, much faster compared with extracting the same result with DMRG methods. However, the techniques of Ref. [71] cannot be straightforwardly applied to the compact boson because of the non-trivial compactification radius, which does affect the entanglement entropy (see for instance Ref. [57][58][59]). In this work, we address the problem of computing the Renyi entropies for integer indices on the compact free boson: our final goal are simple semi-analytic expressions which can be regarded as a generalization of the free-model computations of Ref. [71], but properly keeping in account the compactification radius.
Recent times witnessed an increasing interest towards inhomogeneous generalizations of the compact boson, due to their capability of describing the low energy sector of inhomogeneous systems [72][73][74][75][76][77][78]. In this perspective, we consider the following more general inhomogeneous Hamiltonian [74] where the phase and density fields (respectivelyθ andφ) are canonically conjugated [θ(x), ∂ yφ (y)] = πiδ(x − y). The model is free, but the target space of the phase field is compactified The density field is also compactified, but modulus π. Choosing in Eq. (4) the sound velocity v(x) > 0 and the Luttinger parameter K(x) > 0 to be constant, the usual homogeneous case is recovered. It must be stressed that while the presence of an inhomogeneous velocity does not spoil conformal invariance [72,73], an inhomogeneous K(x) does: therefore it is of utmost importance to have other efficient methods to access entanglement entropies which do not rely on CFT. The parameters v and K are model-dependent and must be properly fixed by an independent analysis of the microscopic system at hand: for example, v and K can be extracted in integrable models by mean of Bethe Ansatz techniques [29,30,74], or numerically in generic systems. As anticipated, for the sake of simplicity we focus on finite systems with OBC, but we expect that the same techniques can be applied to different boundary conditions as well. OBC on the microscopic system impose Dirichlet boundary conditions on the density field of the Luttinger liquid [29] φ( while the phase fluctuation at the boundaries is left free. The conservation of the underlying U (1) symmetry of the microscopic model imposes φ left − φ right ∈ πZ [29].
We will see that, once the difference in the boundary conditions is pinned to a multiple of π, the Renyi entropies are otherwise independent on the actual choice of the boundary fields.
We are primarily interested in a bipartition A ∪ B with B = [x 1 , x 2 ] and 0 < x 1 < x 2 < L, but we also present results for arbitrary bipartitions. Within the fully inhomogeneous case, we provide general expressions in terms of Fredholm determinants. Then, we further advance within the homogeneous case, reducing the problem to a solution of linear integral equations and to the computation of Riemann Theta functions, analogously to the findings of Ref. [57][58][59]. These operations can be quickly carried out on a laptop.
For the sake of clarity, we anticipate and discuss our results in the forthcoming Section II: we present the general expression valid also within the inhomogeneous case and then consider its simplified version for the homogeneous system. The technical derivation is contained in Section III and then our conclusions are gathered in Section IV. A short description of the numerical methods is given in Appendix A.

II. SUMMARY OF THE RESULTS
In this Section we present and discuss our main result, leaving the most cumbersome derivations to Section III. As already stated in the introduction, we are mostly interested in the case of a bipartition A ∪ B with B = [x 1 , x 2 ] and 0 < x 1 < x 2 < L. Generalizations to arbitrary bipartitions are possible as well and are presented in Section III D. In the case of a single interval in the middle of the system, the Renyi entanglement entropy of integer index N is given by the following expression Above, we introduced the following notation: 1. Φ is a linear operator constructed from the two-point connected correlator of the density field Φ(x, y) = φ (x)φ(y) − φ (x) φ (y) . We then define Φ Ω starting from Φ and adding a twist with a phase 2πΩ at the boundaries of B with χ B (x) the characteristic function of the subsystem B, being 1 if x ∈ B and zero otherwise.

det
is understood as a Fredholm determinant, being the domain of the operator the whole system [0, L]. From a practical point of view, the two-point correlator is discretized on a proper lattice and the determinants are computed, then the limit of infinitesimal lattice spacing is considered (see Appendix A). While taking this limit, the logarithm of the product of determinants is well defined up to a constant, which logarithmically diverges in the limit of zero lattice size and it is due to the twist in Eq. (8). Such a divergence is ubiquitous in entanglement's studies (see e.g. [5,6]). Indeed, the lattice discretization introduces a UV cutoff which affects the entanglement according to Eq. (3).
3. The summation in Eq. (7) is over all the possible integers m a ∈ Z and is an example of a Riemann Theta function. Its appearance is due to the non-trivial compactification radius and similar functions are known to be present in the Renyi entropy with PBC [57][58][59] as well.
4. g is the ground state degeneracy, giving rise to the Affleck-Boundary boundary entanglement [64,65]. Within our approach, g emerges as an ill-defined constant (i.e. independent on x 1 and x 2 ) whose computation requires a proper regularization scheme yet to be devised. For the homogeneous compact boson with open boundary conditions, i.e. Dirichlet boundary conditions on the density field, the ground state degeneracy has been computed in Ref. [79] with other methods, as we further comment in the next section.

The matrix M is defined as
and lastly the real coefficients I /N are obtained from the solution of a linear integral equation The solution s /N (x) displays a power-law singularity near the endpoint x 1 and similarly at x 2 . We were not able to analytically characterize such a divergence, but the numerical tests provided in Appendix A confirm the power-law singularity and set 0.5 ≤ µ( /N ) ≤ 1, with µ being 1 only at = N (indeed, in the homogeneous case and = N , Eq. (10) can be explicitly solved and the power-law divergence is analytically confirmed). Note that away from the single point = N , the power-law singularity (11) is integrable, thus I /N is finite. Therefore, the Renyi entropy is well defined except for the overall UV divergent offset caused by the Fredholm determinants, which we have already discussed. 6. We stress that the connected correlator Φ(x, y) is independent from the precise choice of the boundary fields in Eq. (6), which therefore do not affect the Renyi entropies.
Eq. (7) requires as an input the two-point connected correlation function of the density field. For the general case where both the sound velocity and the Luttinger parameter are space dependent, the correlator can be determined solving a two dimensional electrostatic problem [74]. The two-point function is efficiently numerically computed, but a general analytical solution is unknown. If K(x) = const., a putative space-dependent velocity can be absorbed in a change of coordinates along the lines of Ref. [72]. In the new coordinates, the sound velocity is simply v = 1 and the homogeneous problem is easily solved [29]. As we already commented, Eq. (7) can be generalized to bipartitions in an arbitrary number of intervals: the structure of the expression remains the same, but I /N becomes a matrix of dimension (N I −2)×(N I −2) with N I the number of intervals of the bipartition (note that for N I = 3 we retrieve the previous result of a single interval in the middle: indeed the dimension of the matrix becomes one). The matrix M ab is promoted to a tensor M ab → M ij ab where lower indices live in a (N − 1)−dimensional space and upper indices in a space of dimension (N I − 2). Therefore, the Riemann Theta function requires a summation over (N − 1) × (N I − 2) independent integers. We provide the explicit expression, together with its derivation, in Section III D. We now specialize our findings to the homogeneous case, providing a benchmark with known results and further simplifying the expression.

A. The homogeneous case
The homogeneous case stands out as one of the most interesting applications of our approach, where it can be compared with CFT results. As we have already commented, CFT gives important insight into the analytic form of the entanglement entropies, albeit it is not able to completely determine it. Exact, close results are known only at the free fermion point [66][67][68][69]. Here, we greatly exploit the available information: first of all, the free fermion result provides a non trivial check of the correctness of our computation. Furthermore, by a direct comparison of the general form of the CFT expression with our result Eq. (7), we can greatly simplify the latter.
For the case of a single interval, global conformal invariance fixes the entanglement entropy up to a function of the four-point ratio X, constructed from the boundaries and the extrema of the interval [6,57,58] (13) Above, the constant part is due to the microscopic UV cutoff. The result can be generalized to arbitrary bipartitions. The function F N is universal, in the sense that is fully determined by the underlying CFT and not by the microscopic features of the model, which contribute as a constant UV-divergent offset.
The function F N for OBC is so far known only at the free fermion point K = 1, where F N = 1 [66][67][68][69] (see also Ref. [72]): we can use the free fermion result as a non trivial check of our general expression Eq. (7). Furthermore, In evaluating Eq. (7) on a discrete lattice, we testified the expected UV divergent constant offset (independent from K): in each plot we shift the numerical curves with a constant offset which guarantees the agreement with the free fermion analytical result at K = 1. Note, in contrast with the PBC results [57,58], the break down of the K → 1/(4K) duality due to the open boundary conditions. such a comparison allows us to derive a compact expression for F N for arbitrary values of the Luttinger parameter K. The ground state degeneracy for the compactified boson with Dirichlet boundary conditions has been computed in Ref. [79] (see also Ref. [80]), leading to the simple expression Let us now revert to Eq. (7), more specifically we consider the two-point correlator of the density field. In the homogeneous case we have the following simple expression [29] Φ(x, The Luttinger parameter can be simply factorized out in the two-point correlator (this can actually be seen directly from the Hamiltonian (4), where an homogeneous Luttinger parameter can be absorbed in a rescaling of the fields). Establishing analytically the equivalence between Eq. (7) and the CFT result is extremely hard (see e.g. [67][68][69] for a closely related problem), but numerical tests can be easily performed. In Fig. 1 we plot Eq. (7) for different bipartitions and different Luttinger parameters, finding perfect agreement with the known result at the free fermion point. We can now take advantage of the free fermion limit to further simplify the general result (within the homogeneous case). It is convenient to make K explicit in Eq. (7) and define the rescaled correlator UsingΦ rather than Φ in Eq. (7), we note that K disappears from the ratio of Fredholm determinants, having a non-trivial contribution only in the Riemann Theta function. Comparing with Eq. (13) and imposing that for K = 1 the free fermion result is recovered, we conclude (17) whereM is defined as per M, but replacing Φ →Φ in all the expressions. Once again, we stress that the constant term does not depend either on the position of the interval, or on the value of K. Comparing the above with Eq. (13), we can express the universal function F N as a ratio of Riemann Theta functions We have numerically checked that F N is a function of the four-point ratio comparing the r.h.s. of the above expression for several choices of the interval, but having the same four-point ratio. In Fig. 2 we focus on the universal function F N for Renyi index N = 2.
While the lack of a fully analytical expression for M makes difficult to grasp the behavior in X directly from Eq. (18), the trend in K at fixed X is easily understood. In F N >1 , for large K the Riemann Theta function at the numerator exponentially approaches 1, leaving out the algebraic dependence given by the prefactor K (N −1)/2 . For small K we rather reach a plateau, indeed we can readily estimate The K−divergence above is exactly compensated by the prefactor in Eq. (18), leading to a finite result in the decompactification limit. We note that in the infinite volume limit we recover the CFT result for a single interval in an infinite system. The infinite volume limit can be achieved letting L → ∞, with x 2 − x 1 fixed and x 1 , x 2 being kept at a divergent distance ∝ L from the boundaries. This limit simply corresponds to send to infinite the four-point ratio X → ∞, which implies lim X→∞ (1 − N ) −1 log F N (X) = − log √ K, as we numerically verified (see Fig. 2). This K dependence exactly compensates the Affleck-Ludwig term Eq. (14), leaving out as the only contribution the limit of the first term in Eq. (13) which is indeed the homogeneous CFT result for an infinite system (3).
It worth to be mentioned an important difference with Ref. [57,58], where the universal function F N for two intervals in a system with PBC is presented. Similarly to the OBC case at hand, also the PBC result can be written in terms of Riemann Theta functions, but the expression is symmetric for K → 1/(4K). Indeed, the compact boson with PBC presents such a duality [29], which is then reflected into the entanglement entropy.
The weak-strong coupling duality can be understood to be equivalent to swapping the phase and density fields: periodic boundary conditions are of course invariant under such an operation and the model is dual. On the contrary, in the OBC case this is no longer true, since the density and phase fields have Dirichlet and free boundary conditions respectively, thus the duality is spoiled. Indeed, Eq. (18) is not symmetric for K → 1/(4K). The breakdown of the duality is clearly displayed both in Fig. 1 and Fig. 2.

III. DERIVATION OF THE RESULTS
This section gathers all the technical steps which allowed us to compute Eq. (7). Before going deeper into the technical analysis, it is convenient to sketch the main steps. We compute the Renyi entropies through a brute-force evaluation of Trρ N B in a path integral formalism. In this perspective, we introduce the coherent states for the phasedensity field as the eigenvectors of the operatorθ(x) Above, θ(x) is just a number. The coherent states are complete and we explicitly perform partial traces and products in this basis. Let |GS be the ground state of the compact boson, we define the wavefunction W [θ] as its projection on the coherent state |θ A common representation of this wavefunction is that of a path integral in Euclidean time [5,6], leading to a representation of the Renyi entropies in terms of a partition function on a N −sheeted Riemann surface. This was used, for example, in Ref. [57][58][59] to tackle the multi-interval case with PBC. However, we will not go through this route: in Section III A we rather compute the wavefunction in terms of the correlator of the density field, which should be determined by means of independent methods. Being the theory free, the wavefunction turns out to have a gaussian form: simple enough to be handled in the forthcoming steps. We now revert to the problem of computing the reduced density matrices and then the products. The matrix elements of the ground state density matrix in the coherent state basis are simply the product of the wavefunction where W * is the complex conjugated. Considering the system bipartited in two subsystems A ∪ B, we split the field in two parts θ( and θ B is defined in a similar manner. The reduced density matrix of the subsystem B is obtained tracing away the degrees of freedom in A. However, the non trivial compactification radius plays an essential role: the fields must be identified modulus 2π. Formally we can write Furthermore, while computing the productsρ B =ρ B ×ρ B × ... ×ρ B and finally the trace, the bra-ket fields on the subsystem B must be still identified modulus 2π. The (formal) infinite product of Dirac−δs is actually less nasty of what it could seem at first glimpse: the wavefunction vanishes if the field θ is not continuous, therefore we are forced to keep θ A (x) − (θ A ) (x) to be constant in any connected region of the subsystem A. More precisely, imagine Of course, we must sum over all the possible integers n i ∈ Z. Having roughly sketched the battleplan, we now present an helpful short summary of the notation we use, then we proceed with the first step, namely the determination of the ground state wavefunction.

List of notations
First of all, we omit the variable in the derivatives when no ambiguities arise ∂θ(x) ≡ ∂ x θ(x): in the next section, where we need to use space and euclidean time, we explicitly write the argument of the derivative. Besides this exception, we stick to this notation. As we have already done, we consider a bipartition A ∪ B and use superscripts A and B for fields living in A and B respectively, being zero otherwise. Computing the N th Renyi entropy, we need to introduce a replica space with N copies of the original system. We use subscripts to distinguish among the copies: given an arbitrary quantity v (it will be either a field, a real number or an integer), with v j we indicate the quantity living in the j th copy. Furthermore, we callṽ the Fourier transform in the replica space, defined according to the following convention For the sake of a more compact notation, we introduce a vector-matrix representation for the integrals. For example, given two test functions γ(x) and β(x) and an operator A(x, y), we write We also make repetitive use of the characteristic function χ S (x) for an arbitrary set S, defined as it follows A. The ground state wavefunction As we have already mentioned, the ground state wavefunction can be represented through a path integral in euclidean time [5,6]. Even though we ultimately do not explicitly use such a representation, it is nevertheless an useful starting point in studying W [θ]. Let us consider the euclidean time τ : now the fields live on the euclidean plane (x, τ ). The ground state wavefunction is then the partition function with suitable boundary conditions Above, the boundary conditions are such that All the remaining boundary conditions are set free. The path integral representation is convenient to deal with the non-trivial boundaries and easily spot a few symmetries. Indeed, in the path integral we can shift the field φ as per This implies φ (x = 0, τ ) = φ (x = L, τ ) = 0. With this substitution, we get where we set Above, the path integral describes the ground state wavefunction with trivial Dirichlet boundary conditions in the field φ , which is pinned at 0 in x = 0 and x = L. Rather than aiming for a brute force solution, we can reason as it follows: the path integral is gaussian, therefore the result must be gaussian as well For suitable kernels K(x, y) and C(x) yet to be determined. Parity symmetry holds true, i.e. the result does not change if we replace θ (x, τ ) = −θ (x, τ ) and φ (x, τ ) = −φ (x, τ ), together with a reflection in the boundary conditions θ(x) → −θ(x). Therefore, we can conclude C(x) = 0. Furthermore, the path integral must be real: indeed, taking the complex conjugate is equivalent to set φ (x, τ ) → −φ (x, τ ). This implies that the kernel K(x, y), which is a symmetric function, must also be real. Based solely on these symmetry-based arguments, we reach the conclusion In principle, the kernel K(x, y) could be fixed solving the path integral, but we rather prefer to express it in terms of a suitable two-point correlation function, which can be then computed by other means. From a direct analysis of the gaussian path-integral, we can connect the two-point correlator of θ with the (inverse of the) kernel K: However, it is more convenient to look at the conjugate fieldφ. From the commutation rules, we see that ∂ xφ acts as a functional derivative on the wavefunction Which quickly leads to the following simple identity Above, we defined the connected correlator Φ(x, y) = φ (x)φ(y) − φ (x) φ (y) . We stress once again that the connected correlator is independent from the actual boundary conditions φ left and φ right . Note that, obviously, Φ(x = 0, y) = Φ(x = L, y) = 0. For later use, it is convenient to express the wavefuntion W [θ] in terms of the derivative of the field ∂ x θ and of the field at the origin θ(0). In this case, we get where we defined In position x = 0, there is the additional degree of freedom θ(0). While computingρB the fields in the interval A must be identified modulus a constant offset multiple of 2π, therefore the derivatives of the field are the same (we represent this operation with an arrow). Later, while computingρ 2 B and finally taking the trace, we perform a similar operation identifying the fields θB pairwise.

B. Warm up: Renyi entropy for an interval attached to the boundary
In order to build some insight in the required technical steps, it is useful to start with a simple case, namely we choose a bipartition A ∪ B where B is an interval attached to the boundary B = [x, L]. In this case and in the homogeneous setup, the Renyi entropies are fixed by the underlying CFT and the Luttinger parameter K does not play any role, except for the Affleck-Ludwig boundary term. Nevertheless, tackling at first this simpler setup is an useful gymnastic for the forthcoming, more complicated, computations. As we anticipated, rather than considering path integrals in the fields θ, it is more convenient to change basis as and replace the path integral over the field θ with an integration over its derivative and the field at the origin θ(0). As it is clear from Eq. (41), the wavefunction can be factorized in a term which depends only on θ(0) and another one which accounts for the derivatives of the field. For the sake of simplicity, we introduce the following notation While gluing the fields together modulus 2π, the derivatives coincide. The identification of the fields can be understood with the help of Fig. 3: the difference of two fields connected by two arrows is an integer multiple of 2π. This can be implemented identifying the field derivatives and adding the proper constraints at the edges of the intervals, namely (see also Fig. 3) with n j , p j ∈ Z. Above, we used the vector notation Eq. (28) and the characteristic function defined in Eq. (29). We can conveniently subtract the first constraint from the second and redefine new integers m j = p j − n j , then write the N th power of the reduced density matrix as it follows PBC on the replica space are imposed, in such a way that ∂θ B N +1 = ∂θ B 1 . The partition function Z must be fixed imposing Trρ B = 1, as we will do later on. The summation is over all the possible integers {n j } N j=1 , {m j } N −1 j=1 . We stress that the index of the integers {m j } N −1 j=1 runs from j = 1 to j = N − 1: indeed, the N th constraint can be obtained as a linear combination of the others. First of all, note that the sum over the integers n j imposes otherwise the result vanishes. As we anticipated, this quantization of the boundary conditions is required by the presence of the U (1) conserved charge in the microscopic model [29]. Furthermore, the absence of any θ j (0) dependence in the wavefunction W allows us to use the integration over θ j (0) to get rid of the δ constraints. Having N − 1 constraints and N integration variables θ j (0), we are still left with a variable to be integrated, which we choose to be θ 1 (0). Furthermore, we have N − 1 summations over the integers m j and N summations over n. We are left with the following formal expression The integrals and summations above are on the constant function 1 and are formally divergent quantities. We are not entitled to discard the divergent prefactor, which we will analyze later on. The path integral is now a simple gaussian integration. Note that, because of the PBC on the replica space, the linear part in the exponential of W gets averaged to zero and we are left with the following expression Above, we use the vector-matrix notation (28) for the integrals. The quadratic term in the exponential acquires a simple block-diagonal form if we consider the Fourier transform in the replica space. We use the convention of Eq. (27) and obtain the following expression Above, we defined ∂θ = ∂θ A + ∂θ B and Φ Ω is defined as per Eq. (8). Putting together all the terms we get where det Φ/π 3 + Φ /N /π 3 must be understood as the determinant of the operator Φ(x, y)/π 3 + Φ /N (x, y)/π 3 . Imposing Trρ B = 1 we are forced to require Employing this result in Eq. (51), we arrive at the following expression for the Renyi entropy This concludes the derivation for the interval attached to the boundary. Above, we can recognize two terms: the first one actually depends on the chosen interval and, in the homogeneous case, it contains the CFT scaling. It is easy to check that, in the homogeneous case, this term is completely independent on the choice of the Luttinger parameter K. Even though an analytical extraction of the CFT result from the above expression is a hard task, a simple numerical test with the discretization presented in Appendix A ensures the matching with the CFT scaling. The second term does not depend on the actual interval neither on the index of the Renyi entropy and it is undetermined within our approach: we interpret it as the ground state degeneracy resulting in the Affleck-Ludwig boundary entropy [64,65] Indeed, the very same term appears also in computing other bipartitions, strengthening once again the above identification: we will see in the next section in the case of an interval in the middle of the system and then later for general bipartitions in Section III D. A first principle derivation of the ground state degeneracy g within the present method would require a proper regularization of the integral (and sum) over the target space appearing into Eq. (54). Such a regularization lays beyond our current understanding.

C. An interval in the middle
We now turn to the main case of interest, namely B = [x 1 , x 2 ] with 0 < x 1 < x 2 < L. We still conveniently integrate over the derivative of the field and its value at the boundary, keeping the same notation as Eq. (44). The trace of the reduced density matrix has a similar expression to Eq. (46), however with a further delta function (see Fig. 4). Indeed, the following constraints must be enforced It is convenient to subtract from the second and third constraint the first one, in this way we reach the following expression for the trace of the reduced density matrix As we did in the previous section, we use the integration over θ j (0) to get rid of the first set of constraints. As before, the summation over n j becomes trivial once Eq. (47) is imposed. However, after these operations, a set of non trivial constraints is still left out Above, note the appearance of the same prefactor that, in the previous section, resulted in the ground state degeneracy (54). The normalization Z is of course independent on the chosen bipartition and we already computed it in Eq. (52). Employing this result together with Eq. (54) we can write In order to take care of the δ−functions, we use a proper integral representation for each constraint. Therefore Above, we conveniently added a dummy integration over an additional field λ N that is then pinned to 0 by the Dirac−δ. This makes the integrand more symmetric and easier to be handled. Furthermore, we insert another integral representation for δ(λ N ) Rather than using λ j , we introduce new variables ζ j in such a way Which implies λ j+1 − λ j = ζ j . Obviously, we have λ N = N i=1 ζ i . The real numbers ζ j and integers z j are Fouriertransformed in the replica space according to the usual convention Eq. (27) and we reach the following expression We now change variable in the summation and define new integers m j = N i=j z i . This is a one-to-one map provided we set m N = 0. We use these new integers to perform the summation, since the constraint associated to the ζ variables gets simplified. We now compute the gaussian path integral, leaving out the integration over λ j and ω for a second moment. We get rid of the imaginary part performing an analytical continuation to complex fields In terms of the fields in the Fourier space, this amounts to a transformation The shift must be chosen in such a way to absorb the linear term, i.e. we require After such a shift, we get The last path integral is now performed as we did in the previous section. We then reach the following expression Where we defined the real numbers I /N as We stress that this definition is equivalent to Eq. (10). The last step requires computing the gaussian integral in the variables ζ and than that on ω, which leads to the following result Note that Trρ B = 1, as it should be. Since the sum must be performed on real integers m j , it is more convenient to write the exponent in the real space rather than in the Fourier one. Therefore, let us define the matrix M ab as per Eq. (9) and write Taking the logarithm, we are immediately lead to Eq. (7).

D. Generalization to arbitrary bipartitions
The same analysis of the previous section can be further generalized to arbitrary bipartitions: the main difficulty consists in choosing the right notation. We now imagine a collection of points on the real axis {x i } N I −1 i=1 such that 0 < x i < x i+1 < L and consider the bipartition Above, we set x 0 = 0 and x N I = L. We now focus on the Renyi of integer index N . We keep the same notation for the fields as in the previous section and retain the vector-matrix notation as well. We define ∂θ j = ∂θ A j + ∂θ B j . Then, it is not hard to understand that the non trivial compactification radius induces the following set of non trivial constraints If we choose N I = 3 we are back to the case of a single interval in the middle studied in the previous section. We then need to compute the following path integral Above, the summation is over all the possible integers z i j . As before, we introduce an integral representation of the δ function for each constraint (77) Now, generalizing what we have previously done, we change variables for λ i j and on the integers z i j as Rather than summing over z i j we can sum over the integers m i j . We furthermore introduce additional variables ω i to write an integral representation of δ(λ i N ). Once we consider the Fourier transform in the replica space, we reach a generalization of Eq. (65) Above, Φ Ω is still defined as Eq. (8), with the difference that now B is made of several disconnected intervals. The gaussian path integral is carried out through an analytical continuation in the complex plane, as we did in the single interval case and we then reach the analogue of Eq. (70), but with extra indices. After that, the gaussian integral on the ζ and ω variables can be computed and the following final expression is reached Above, we defined the tensor It worth to mention that, within the homogeneous case, we could generalize the content of Section II A: conformal symmetry fixes the Renyi entropies apart from an additional contribution (1− N ) −1 log F N . F N is a function of all the four-point conformal invariant ratios which can be constructed with the extrema of the intervals (and the boundary points). The Renyi entropies at the free fermion point are known explicitly for arbitrary bipartitions (see e.g. Ref. [66]), thus by comparison with Eq. (80) it is possible to fix F N in terms of Riemann Theta functions, generalizing Eq. (18).

IV. CONCLUSIONS AND OUTLOOK
This work has been dedicated to the computation of the Renyi entanglement entropies of integer index for the 1d compactified massless boson. While this problem has already been posed (and solved [57][58][59]) for periodic boundary conditions, here we rather address open boundary conditions, which are of utmost importance for tensor networks simulations and the low energy physics of confined systems [72]. We address this problem for the ground state of an inhomogeneous generalization of the compactified boson, described by locally-varying sound velocity and Luttinger parameter: in this case, we provide Fredholm determinant-like expressions which generalize previously known methods for free systems [71]. Within the homogeneous case, which is described by conformal symmetry, we pushed further our analytical investigation reducing the problem to solving linear-integral equations and computing Riemann Theta functions, which account for the role of the non-trivial compactification radius. Several open questions and interesting applications are left for future investigations. We presented a non-trivial benchmark of our result on free fermions systems where the analytical expression was known, but the final goal is the description of ground states of interacting models with K = 1: it would be extremely interesting to compare our findings with DMRG simulations. We provided general expressions for the Renyi entropies having as an input the density-density correlation function: its application to spatially inhomogeneous or out-of-equilibrium setups is a natural question to be investigated. For what it concerns possible open issues, there are at least three main points.
The first question concerns the determination of the Affleck-Ludwig boundary term within our approach, which leaves it undetermined in absence of a proper not-yet-devised regularization scheme. Within the homogeneous case, we can match this constant with the existing literature, but the application of our findings to inhomogeneous systems ultimately needs a direct computation of the ground state degeneracy.
Secondarily, even considering the homogeneous case, the computation of the Renyi entropies requires solving some linear integral equations, task that we performed numerically: having a fully analytical expression would be extremely important, especially for a systematic small and large distance expansion.
The third natural question concerns the possibility of analytically continuing our result to real Renyi indices, finally reaching the Von Neumann entanglement entropy. Solving this surely appealing problem seems quite hard, due to the intricate index-dependence of our result: a similar question has already been asked in the closely related PBC problem of disjoint intervals [57][58][59] but, as far as we know, it does not have a direct solution yet, albeit numerical approaches can be used [81,82]. In Ref. [33,34] the PBC problem has been attacked with conformal blocks techniques, resulting in an expansion of the Renyi entropies in the distance among the intervals. Remarkably, each term of the series can be analytically continued and an expansion for the Von Neumann entropy is achieved. However, establishing the connection with the exact result of Ref. [57][58][59] is a difficult question. It is natural to wonder if conformal blocks techniques can be used to study the OBC case as well. Lastly, we would like to mention the intriguing possibility of extending the methods here presented to other measures of entanglement, for example the entanglement negativity [83][84][85][86][87][88][89], using the replica approach in the path integral formalism [90,91]. i| with the decreasing of the lattice spacing. The perfect linear growth testifies the power-law singularity of the solution to the continuum integral equation (10) ∼ |x − x1| −µ (in the current example, x1 = 0.25). The slope of the linear growth allows to numerically determine the exponent µ (Subfig. (b)), which we verify satisfies 0.5 ≤ µ ≤ 1. We considered 0 ≤ Ω ≤ 0.5, since the case 0.5 ≤ Ω ≤ 1 is trivially connected to the previous one through complex conjugation.
Out of this repesentation, we construct a discrete approximation truncating the series Where the index j runs from 0 to 2N d − 1 and are a discretization of the real space. Then, we construct the discrete correlation function [Φ d ] jj as per We now want to consider an interval [x 1 , x 2 ] = [i/N d , i /N d ], we then construct the discrete version of the Φ Ω operator as per The ratio of Fredholm determinants is discretized as per The matrix [Φ d ] −1 has the following explicit expression As expected, we see that the determinant has a well defined continuous limit up to a constant (i.e. it does not depend on the actual position of the interval) which diverges with the lattice spacing. We now discuss the computation of the I /N coefficients. We define the discrete characteristic function χ d [i,i ] as it follows Then define a vector [s d Ω ] j as the solution of the linear equation below The presence of the lattice spacing ensures that s d Ω reproduces the solution of the continuum problem in the limit a → 0, i.e. [s d Ω ] j → s Ω (j/N d ). Finally, the cofficients I Ω are computed as The presence of the sudden phase twist in the definition of Φ Ω causes a power-law singularity in the solution of s Ω (x) around x ∼ x 1 and x ∼ x 2 . In Fig. 5 (Subfig. (a)) we numerically study the singular behavior, testifying a power-law singularity |x − x 1 | −µ . The exponent of the singularity is numerically extracted in Fig. 5 (Subfig. (b)), showing the singularity is always integrable for Ω = 0 mod 1 and ensuring the finitness of I Ω , as wee commented in the main text.