Entanglement entropies of inhomogeneous Luttinger liquids

We develop a general framework to compute the scaling of entanglement entropy in inhomogeneous one-dimensional quantum systems belonging to the Luttinger liquid universality class. While much insight has been gained in homogeneous systems by making use of conformal field theory techniques, our focus is on systems for which the Luttinger parameter K depends on position, and conformal invariance is broken. An important point of our analysis is that contributions stemming from the UV cutoff have to be treated very carefully, since they now depend on position. We show that such terms can be removed either by considering regularized entropies specifically designed to do so, or by tabulating numerically the cutoff, and reconstructing its contribution to the entropy through the local density approximation. We check our method numerically in the spin-1/2 XXZ spin chain in a spatially varying magnetic field, and find excellent agreement.


Introduction
Many quantum critical systems in one-dimension (1D), including various quantum gases or spin chains, belong to the universality class of Luttinger liquids [1][2][3][4][5][6][7]. From an effective field theory (FT) perspective, Luttinger liquids are the one-dimensional (1D) systems whose low-energy description is provided by a massless free boson FT in 1 + 1D, or in other words a U (1) conformal FT [8]. They are ubiquitous in the 1d quantum realm because U (1) symmetry itself is ubiquitous-it is associated to particle number conservation, often naturally present in spin chains or quantum gases-, and because of the limited amount of relevant perturbations-often forbidden by symmetry-that could drive the system away from the renormalization group massless fixed point.
One celebrated result about Luttinger liquids is that the entanglement entropy of a subsystem I = [x 1 , x 2 ] in the ground state of an infinite homogeneous (i.e. translation invariant) system grows logarithmically with the size of the interval I [9][10][11][12]. Specifically, for the entanglement entropy with Rényi index α, the result reads [12] S α (x 1 , up to corrections which go to zero as the length |x 1 − x 2 | goes to infinity. Recall that the Rényi entanglement entropy of a pure state |ψ , which will always be the ground state in this paper, is defined by where ρ I = TrĪ|ψ ψ| is the reduced density matrix for subsystem I, and Ī is its complement. The constant a α in (1) is non universal. It cannot be fixed by global conformal invariance and its analytic determination can be extremely cumbersome even in free systems [13][14][15][16]. When comparing formula (1) with numerical results in concrete microscopic systems, a α is often treated as a fitting parameter. Motivated by practical aspects of 1D quantum gases in non-uniform trapping potentials [4,7,[17][18][19][20][21][22][23] or quantum spin chains in non-uniform magnetic fields [24][25][26][27][28], the purpose of this paper is to revisit the entanglement entropy of Luttinger liquids in the context of inhomogeneous 1d quantum systems, which are currently attracting a lot of attention, see e.g. [20,[28][29][30][31][32][33][34][35][36][37][38][39][40][41][42] for recent works in that direction. We focus on inhomogeneous Luttinger liquids, namely on quantum systems whose effective FT description corresponds to a free boson theory in a curved metric and with a spatially varying coupling strength or Luttinger parameter K [4,20,22,23,28,43,44]; the theory will be briefly reviewed in section 2.
In an inhomogeneous quantum system, we expect the entanglement entropy of an interval I = [x 1 , x 2 ]-assumed to be larger than any microscopic length scale such as lattice spacing or interparticle distance-to be of the form where S FT (x, y) is the entanglement entropy calculated in the FT with a certain UV cutoffthat cutoff is discussed below-, and c α (x), c α (y) are the microscopic contributions to the entanglement entropy around the two points x and y . In homogeneous systems with entanglement entropy given by equation (1), c α (x) is independent of x. It can be thought of as 1+1/α 12 log(a FT α /a α ) where a FT α /a α is the ratio of the UV cutoff in the FT and a microscopic length scale in the quantum system.
In writing equation (3), two crucial features must be stressed. First, two microscopic models sharing the same low energy description have the same contribution S FT α (x, y) to their entanglement entropies. This observation is at the root of the method developed in this paper: instead of considering a complicated interacting quantum model, it is possible to directly start from its effective FT description as an inhomogeneous free boson theory. However, in replacing the true microscopic quantum model by its effective low energy description, one is left with the difficult task of finding the position-dependent microscopic contribution c α (x).
On the other hand, the second feature apparent from equation (3) is that the microscopic contribution c α (x) is completely determined by the system in proximity of the points x 1 and x 2 . This offers at least two possibilities, both of which we explore in great detail in the rest of this paper.
The first possibility is to consider a suitable combination of the entanglement entropy such that only universal features remain. For instance, one can define the -regularized entanglement entropy as (for similar definitions which appeared in the literature, see for instance [45,46]). The effect of this regularization is to replace the cutoff-dependent contribution to the entanglement entropy by a function of the length which is determined only by the universal FT part. Indeed, since the length can be fixed arbitrarily, one can chose to be much larger than the microscopic scale, yet much smaller than the length scale of the inhomogeneity. Then S ( ) , as can be seen by plugging equations (3) into (4). The second possibility is to numerically tabulate the values of c α by studying the homogeneous problem for different values of the external parameters. The inhomogeneous system can then be addressed by promoting c α to a function of the position x through the local density approximation (LDA). In this approach, a numerical study of the homogeneous system together with the low energy effective theory gives access to the full entanglement entropy (3).
Below, we provide extensive checks of these two possibilities, thus demonstrating that the inhomogeneous Luttinger liquid theory can efficiently capture the scaling part of the entanglement entropies even in interacting inhomogeneous systems. The lack of conformal symmetry in Luttinger liquids with spatially varying coupling constant [23] makes it impossible to apply standard methods of conformal invariance and to reach nice closed analytical results. Nevertheless the gaussianity of the theory still leads to expressions for entanglement entropies that can be efficiently evaluated numerically. The method we develop here is versatile as it applies to all quantum systems described by (inhomogeneous) Luttinger liquid theory, some of which-including the Lieb-Liniger model-can be very difficult to simulate by other means in general.
Finally, it is important to stress that we are dealing with a free compact boson, and that the compactification has non-trivial consequences in the computation of entanglement entropies. This is well understood in the homogeneous case with periodic boundary conditions [47][48][49], however for open boundary conditions the entanglement entropy was studied only recently by one of us [50]. In this paper we focus on open boundary conditions and make extensive use of the results of [50].
The paper is organized as follows. In section 2 we briefly review the theory of the inhomogeneous free boson in 1d and explain how to calculate entanglement entropies starting from the inhomogeneous Luttinger liquid Hamiltonian. We apply this formalism in sections 3.1 and 3.2 to study the entanglement entropy of an inhomogeneous XXZ chain. We perform DMRG checks in the XXZ chain to validate our approach. In section 3.1 we present our results for the -regularized entanglement entropy with Rényi index α = 2. Section 3.2 is focused on the von Neumann entanglement entropy (α = 1) and Rényi entropy with index α = 2 of a half-chain using numerical tabulation of the cutoff-dependent contribution c α (x). We conclude in section 4. The more technical aspects of our results are presented in full detail in the appendices.

Entanglement entropy in the inhomogeneous free boson theory
The FT that underlies inhomogeneous Luttinger liquids (see e.g. [4,17,20,22,23,28,43,44,51]) is the free massless boson in a metric g with a spatially-dependent coupling constant-which is called the Luttinger parameter K in this paper-, see the action (9) below. In this section we start by briefly reviewing this theory in Hamiltonian formalism (better suited than lagrangian formalism for the purposes of this paper). Then we explain how, discretizing the theory on a finite grid, one can obtain formulas that allow to calculate the entanglement entropy in terms of the Green's function of that free theory.

Definition and Hamiltonian
The Hamiltonian of a homogeneous Luttinger liquid defined on the interval [0, L] is (we follow the notational conventions of the textbook of Giamarchi [5]): where v is the sound velocity in the fluid, K is the Luttinger parameter which encodes the renormalized interaction strength between the constituents of the microscopic model, and φ and πΠ = ∂ xθ are canonically conjugated bosonic fields, Throughout this paper, we focus on Dirichlet boundary conditions for the field φ (as in [4,21,23]),φ Importantly, the fields φ , Π are compactified: θ is identified with θ + 2π, while φ is identified with φ + π. One way of understanding this is to look at the identification of operators in the microscopic theory with local operators in the FT. For instance, in a microscopic model of bosons or fermions, the particle density operator has an expansion of the form n(x) = n 0 − 1 π ∂ xφ (x) + . . ., where n 0 is the mean value of the density. Thus, adding or removing one particle corresponds to a ±π-jump of the field φ . Such ±π jumps are most conveniently interpreted as non-zero windings when φ takes values in the circle R/(πZ). Similarly, the field θ is usually interpreted as a phase, and therefore takes values in the circle R/(2πZ). We will see in section 2.3 that this compactification induces non-trivial technicalities in the computation of entanglement entropies.
Finally, we stress that for microscopic models that can be mapped to free fermions, the Luttinger parameter K is always equal to one (see for instance the discussions of that point in [22,52]), but in the presence of interactions it can be any positive real number. Now let us turn to the inhomogeneous case. In the limit of very slow variations of the mean density n 0 (x), it is still possible to describe the low-energy fluctuations of the system by the two canonically conjugated fields φ (x) and Π (x), if one can identify the correct inhomogeneous generalization of the Hamiltonian (5). Of course, the inhomogeneous generalization must be such that, locally, it reproduces the results of the homogeneous model. The simplest way of ensuring this is to promote the parameters in the Hamiltonian (5) to functions of the position v → v(x), K → K(x), the local values of which are fixed by the underlying microscopic theory through LDA [4,20,22,23,28], Notice however that, as discussed in detail in [22,23,52], promoting the model to inhomogeneous v(x) and K(x) leads to possible ambiguities. Different choices are in principle allowed: for example, substituting 1/K(∂ xφ ) 2 in the homogeneous Hamiltonian (5) by would lead to another inhomogeneous Hamiltonian which would look equally valid. In [23,52] it is properly argued that, in lagrangian formalism, the correct action reads (the field h(x, τ ) in these references where x = (x 1 , x 2 ) = (x, τ ), τ is the imaginary time coordinate, and the two-dimensional metric g ab has euclidean signature and reads in our case In appendix A we show that the Hamiltonian (8) is indeed the one that is obtained from the action (9). The covariant form of the action suggests the change of variables which puts the metric in the form ds 2 ∝ dx 2 + dτ 2 [20] and x ∈ [0, 1]. Hence the positiondependent sound velocity drops from most computations, and we find that it also facilitates numerical calculations (see appendix C). The position-dependent Luttinger parameter K(x), however, cannot be absorbed in a simple change of variables, and the ground state correlations can no longer be calculated analytically, contrary to the case when K is a constant. This is an important complication inherent to interacting inhomogeneous Luttinger liquids. The ground state correlations can all be expressed in terms of the Green's function of a generalized laplacian in 2D, ∇ x 1 K(x) ∇ x , and the problem is formally analogous to calculating the Coulomb energy of distributions of electric and magnetic charges in a 2d electrostatic problem with a spatially varying dielectric constant 1/K(x). We refer to [23] for a more detailed discussion of ground state correlations in the inhomogeneous Luttinger liquid.
In this paper, we rather rely on a numerical evaluation of the ground state correlations which is based on a discretization of the system on a finite grid. We refer to appendix C for further details on our numerical procedure.

Entanglement entropy of an interval at the edge
Let us consider a bipartition I ∪Ī of the system, with I = [0,x 1 ] and Ī = [x 1 , L]. For this choice of bipartition, one does not need to be particularly careful about the compactification of the field φ (x) (see figure 1). The calculation can be done straightforwardly, following for instance the review [53]. Here we write the calculation for continuous fields, without paying attention to problems of UV regularization. In practice, the formulas we present here will be used in sections 3.2 and 3.1 by discretizing the system on a finite grid, thus naturally providing a UV regularization (see appendix C for details).
We define the bosonic annihilation/creation operatorŝ The mean density n 0 (x) appears because of dimensionality. We also define the correlation operator C(x, y) in block form with the following notational convention for product of operators: Let us call C I the correlator restricted to the subsystem I = [0,x 1 ]. Since the correlation functions of the whole system satisfy Wick's theorem, correlation functions restricted to the subspace I automatically satisfy Wick's theorem as well. This implies that the reduced density matrix ρ I is gaussian in terms of the bosons and unambiguously fixed by the correlator C I . In particular, there exists a canonical transformation to new bosonic fields γ( which diagonalizes the correlation operator: This transformation also diagonalizes the reduced density matrix, Figure 1. Cartoon of two configurations of particles in the microscopic model, and the corresponding expectation value of the field φ (x), which jumps by −π at the position of each particle according to n(x) = n 0 (x) − 1 π ∂ xφ . Notice that the two configurations of particles are globally different, however they are identical inside the subsystem I. In terms of the field φ (x), the two configurations inside the interval I are viewed as different unless one identifies φ (x) with φ (x) + π. In other words, if one simply treats the field φ (x) as real-valued, then one is overcounting the configurations inside the interval I; if instead the field φ (x) is viewed as taking values in R/(πZ) then it correctly encodes particle configurations. It is extremely important to properly take into account the compactification in the calculation of the entanglement entropy of an interval I in the bulk. Notice that, because of the Dirichlet boundary condition satisfied by φ (x), there is no problem of overcounting when the interval I touches the boundary.
The entanglement entropy is then related to the eigenvalues λ(k): We see that the non-trivial part of the calculation is the diagonalization of the correlation operator C I , equations (14) and (15). Let us elaborate on that point. In order to preserve the canonical commutation rules, the change of basis (16) must be symplectic, It is convenient to take advantage of the symplectic transformation to observe that (19) Thus, the spectrum of C I σ gives access to λ(k), which then gives the entanglement entropy through equation (17). Again, we emphasize that, in practice, we discretize the inhomogeneous Luttinger liquid on a finite grid (appendix C), so the correlation matrix C I is a finite-size square matrix and the spectrum of C I σ is obtained numerically.

Entanglement entropy of an interval in the bulk: the role of compactification
Next, we turn to the case of a bipartition I ∪Ī with I = [x 1 , This situation is more complicated than the one in the previous section because of the compactification of the field φ which now plays a role. To see this, consider the two configurations of particles in the microscopic model drawn in figure 1. Although they are identical inside the interval I, the field φ (x) takes different real values inside I. One sees that the configurations are identical only if one identifies φ (x) with φ (x) + π. Not taking this compactification into account would lead to wrong results for the entanglement entropy.
It is well known that dealing with compactification in the context of calculating entanglement entropies is not an easy task. Indeed, the same problem is faced when calculating the entanglement entropy of a multi-interval in an infinitely long system [47][48][49]54], and there only the Rényi entropies with integer α 2 have been accessed through a direct analytical calculation. No analytical continuation to continuous Rényi parameter α has been found so far even in that well-studied case, albeit the analytical continuation can sometimes be performed numerically [55,56]. In [57,58], by means of conformal blocks techniques, the Rényi entropies for disjoint intervals were analyzed in the PBC case, resulting in an expansion in the small distance among the two intervals: such an expression can be analytically continued term by term, giving access to the same expansion for the von Neumann entropy. However, connecting this method with the exact results of [47][48][49]54] is still an open challenge.
Therefore, in this section we also restrict ourselves to the case of Rényi index α integer and 2. The calculation leading to this result was done very recently by one of us in [50], and we refer to that reference for details of the derivation. The result of [50] reads where the constant offset does not depend on x 1 , x 2 , but may depend on α. The summation is over m 1 , m 2 , . . . , m α−1 ∈ Z. The matrix M is a square matrix of size (α − 1) × (α − 1), while Φ and Φ ω are operators: their determinants are computed by means of a suitable discretization of the continuum model on a finite grid, before taking the limit of zero lattice spacing. In doing so, the entanglement entropy is well defined except for a constant offset, which diverges logarithmically in the lattice size [50], in agreement with the homogeneous result (1).
The function I ω and the matrix M ab are defined as Here s ω (x) solves the following integral equation The operator Φ is nothing but the connected correlation function in the ground state, defined on the interval [0, L], where χ I is the characteristic function of the interval I: χ I (x) = 1 if x ∈ I and χ(x) = 0 otherwise. The last additive constant C α was overlooked in [50], where the primary focus was on the coordinate dependence of the expression. The validity of the results in [50] holds true apart from an overall constant, which is left ambiguous in the definition of the multi-sheet partition function defining the α-Rényi entropy. However, determining the exact value of the entropy is essential for the applicability of our method (in particular for the results in section 3.1): the constant ambiguity can be removing simply by imposing that, in the limit x 1 → 0, the results of section 2.1 are recovered. This calculation is carried out in appendix D, resulting in Again, we emphasize that, in practice, we evaluate all these expressions numerically, by discretizing the system on a finite grid (appendix C). Thus, we are always manipulating large but finite matrices, and formula (20) is evaluated through standard linear-algebra routines. While taking the continuum limit, the logarithm of the ratio of determinants splits into a universal part and a UV-divergent part ∝ log , where is the lattice spacing. All the other quantities appearing in equation (20) are well defined in the continuum limit: it is possible to show [50] that s ω (x) acquires power-law singularities at the edges of the interval [x 1 , x 2 ], however these divergences remain integrable ∼1/|x − x 1,2 | max(ω,1−ω) for any 0 < ω < 1. We observe that our choice of considering the α = 2 Rényi entropy guarantees the best convergence of the expression (20).

Entanglement entropy in the inhomogeneous XXZ spin chain
Now that we have presented the general ideas and methods, we turn to a specific microscopic model where we can check our approach numerically. We focus on the XXZ spin chain in a smoothly varying external magnetic field Here Ŝ x,y,z are the spin-1/2 operators, and N is the spin chain's length, assumed to be large. The local magnetic field h(x) is some smooth function, so that h( j/N − 1/2) varies very slowly when N is large (we shift the argument of the magnetic field in such a way the argument of h(x) is zero when j = N/2). We stress that our chain has open boundary conditions. All numerical calculations in the chain were performed using a DMRG algorithm based on the Itensor [59] C++ library. When the magnetic field h is constant, the model is integrable [60,61]: this provides a powerful set of tools which allows an exact determination of its spectrum, together with an exact computation of the Luttinger parameters v, K, as sketched in appendix B. Within the paramagnetic phase |∆| < 1, the ground state of the homogeneous XXZ chain is gapless for |h| < 1 + ∆. In our numerics we focus on the values ∆ = ±1/2 and h > 0, but other choices of the parameters would lead to similar results, as long as one stays in the gapless regime. For details on the bosonization of the XXZ spin chain we refer to the existing literature (e.g. [4,5]). Here we only need to determine the local velocity of gapless excitations v(x) and the Luttinger parameter K(x), which we do through the LDA.
Hereafter, we chose the magnetic field in such a way that only a subregion of the spin chain has non-trivial magnetization −0.5 < m < 0.5. The magnetization is maximal outside that region, and the spins do not fluctuate there. The inhomogeneous Luttinger liquid decribes only the subregion of the spin chain where the spins are critically fluctuating. Thus we stress that the spatial domain in which our Luttinger liquid lives is not the whole chain but only that subregion, which is the only one giving non-zero contribution to the entanglement entropy to the leading order (see figures 2 and 3). The subleading edge behavior of the entropy corre sponding to near-saturating magnetic field can also be obtained exactly [39,62,63], but belongs to a different edge universality class. Our focus is on bulk Luttinger liquid behavior in this paper. In this section the bipartition consists of the two intervals I = {1, . . . , j} and Ī = { j + 1, . . . , N} of spins in the XXZ chain (26).

-regularized (Rényi) entanglement entropy in the inhomogeneous XXZ spin chain
We approach the FT using the following flat discretization (i.e. in the isothermal coordinates equation (11)),Ĥ where the discrete fields satisfy the canonical commutation relations and = M −1 is the lattice spacing is isothermal coordinates. We impose Dirichlet boundary conditionsφ Further details about our discretization of the free boson theory are deferred to appendix C. Of course, both the XXZ spin chain and the (discretized) Luttinger Liquid bear their own cutoffs, however, in the same spirit of equation (4), from the equation (3) we can construct cut-off independent identities such as where, making use the entanglement entropy in the XXZ model S XXZ α and in the discretized Luttinger S LL( ) α , we defined the regularized entropy While equation (30) is expected to hold true for arbitrary Rényi entropies (and therefore also in the limit corresponding to the von Neumann α → 1), computing equation (31) requires the knowledge of the entropies of an interval embedded in the bulk ( j, j + ). In this case, the non trivial compactification of the bosonic field plays an important role and, to the best of our knowledge, equation (20) [50] is the only available result for open boundary conditions: this allows us to address the Rényi entropies of integer indexes, but leaves out the von Neumann entropy from our analysis. In figure 2 we compare the regularized Rényi entropy (31) for α = 2 in the XXZ model and in the discretized Luttinger, for different values of the interaction ∆ and the magnetic trap, finding good agreement. One can observe oscillations, which are much more pronounced for ∆ = 0.5 than for ∆ = −0.5. This observation is consistent with known results for corrections to scaling in homogeneous systems [64], which are controlled by the Luttinger parameter. While they disappear in the thermodynamic limit, the corre sponding exponent becomes smaller as K is smaller, so they significantly affect finitesize calcul ations when ∆ approaches 1.

Results from numerical tabulation of the UV contribution to the entropy
According to the discussion around equation (3) in the introduction, the entanglement entropy in the XXZ chain will be of the form where 'FT' stands for the FT result, which itself depends on a FT cutoff. To evaluate the FT part we rely on the discrete Luttinger liquid Hamiltonian (27): the entanglement entropy in the discretized Luttinger liquid (which we note 'LL( )') satisfies a relation similar to equation (32), with a different non-universal O(1) part c α , Here the O(1) part depends on the microscopic parameters of the discretized Luttinger liquid: the lattice spacing , and the local parameters v(x) and K(x) entering the Hamiltonian. Thus, we have where d should be a function of the local parameters in the model only, so we can fix it by the LDA; we see that it depends only on the local magnetic field h x . Notice that, on the discrete Luttinger liquid side, the local parameters v(x) and K(x) do themselves depend on h(x) through the LDA.
To evaluate d XXZ/LL( ) α (h), we can focus on a homogeneous spin chain. For a fixed magnetic field h, the difference is constant as a function of j , up to terms that go to zero in the thermodynamic limit. We tabulate the measured constants as a function of the field h, and this gives us access to the function d XXZ/LL( ) α (h), in the homogeneous spin chain. There is one minor subtlety that is left to discuss when one goes from the homogeneous to the inhomogeneous case. Indeed, the discrete Luttinger liquid Hamiltonian does not need to rely on the spatial discretization corresponding to the spin chain one. In general there is a change of variables j →x such that the spatial discretization of the discrete Luttinger liquid corresponds to x ∈ Z. Then one needs to keep track of that change of variables. This is done as follows. Recall that in FT the calculation of entanglement entropy maps to the one of a correlation function of twist fields [65] with scaling dimension α−1/α 12 . Then a change of variables generates an additive term in the final result for the entanglement entropy, which is the logarithm of the Jacobian of the change of variables, times the scaling dimensions. Thus, if we want to use a general coordinate x for the discretized Luttinger liquid, equation (36) becomes instead (36) In figure 3 we compare the data of several DMRG simulations of the full inhomogeneous system, with those extracted from the (discretized) Luttinger liquid, using equation (36). We focus on the von Neumann entanglement entropy S 1 and the Rényi entropy S 2 : the agreement is excellent in both cases.

Conclusion
We have developed a method to calculate the entanglement entropy of Luttinger liquids in inhomogeneous situations, where the Luttinger parameter K becomes position dependent. The crucial ingredient in our analysis is a careful treatment of the contribution to the entanglement entropy coming from the UV cutoff, both in the microscopic model one is interested in (in this paper, the XXZ chain) and in the discretized free boson FT one uses to do the calculations.  (34). We focus on the von Neumann entropy S 1 and on the second Rényi S 2 for ∆ = ±0. 5 and several magnetic trap (see insets). The XXZ spin chain has N = 1024 lattice sites. In the case of the discretized Luttinger liquid we considered isothermal coordinates and used 200 sites, which were enough to reach convergence. In contrast with figure 2, the XXZ data show smaller oscillations around the Luttinger liquid prediction.
Contrary to homogeneous situations, where this contribution is simply a constant often treated as a fitting parameter, here that constant depends on position. It is possible to construct combinations of entanglement entropies, like the -regularized entropy we considered (see also [45,46] for similar definitions of regularized entanglement entropies), such that the non-universal position-dependent constant disappears. Alternatively, we showed that this position-dependent contribution to the entanglement entropy can be tabulated, and then used in numerics to arrive at concrete predictions for the standard (i.e. not regularized) entanglement entropy. In addition, our analysis provides a numerical benchmark for the results of [50] in the interacting case (K = 1), which was not tested in the original reference.
It would be very interesting to extend our results to deal with out-of-equilibrium inhomogeneous situations in critical spin chains or 1D quantum gases. One paradigmatic example of such situations corresponds to the Riemann problem in hydrodynamics, where two semiinfinite systems with different thermodynamic characteristics are suddenly joined together [31,[66][67][68][69][70]. In the particular case of two initial half-systems at zero temperature (but, say, at different filling fractions), the system must be described by a time-dependent inhomogeneous Luttinger liquid. The free fermion case is mostly understood [20], however in the truly interacting case where the Luttinger parameter K is expected to depend on x and t, there is no theory allowing a calculation of the entanglement entropy, to our knowledge. A number of numerical observations on entanglement entropies have been made [71][72][73][74], and it would be very nice to adapt our results to the time-dependent case and make contact with those works. Here again, the non-universal position-dependent (and also time-dependent) contribution to the entanglement entropy will require a careful treatment. This time, to tabulate it one might need to look at homogeneous states that are not just the ground state of the homogeneous problem, but highly excited states. In particular, in some cases it might be necessary to go beyond the simple Luttinger liquid and consider so-called 'split Fermi seas' [42,[75][76][77][78], and calculate entanglement entropies in those, which could be hard. Alternatively, focusing on -regularized entanglement entropies should be easier in those out-of-equilibrium situations; this could be a good starting point to investigate this class of problems.

Acknowledgments
We are grateful to Paola Ruggiero, Jacopo Viti and Pasquale Calabrese for inspiring discussions. We are especially grateful to Paola Ruggiero for pointing out [45,46]. AB acknowledges the support from the European Research Council under ERC Advanced Grant 743032 DYNAMINT.

Appendix A. Inhomogeneous Luttinger liquid: from Lagrangian to Hamiltonian
In this appendix we work with Lorentzian signature. The action of the inhomogeneous Luttinger liquid is ]. The entries of the metric g and the Luttinger parameter K depend on x and t. The canonical momentum associated to φ(x) is and the Hamiltonian derived from the action S is then In the (t, x) basis, the metric is the 2 × 2 matrix where u ± v are the velocities of sound waves in the fluid. In this paper we are working with a time-independent Luttinger liquid with u = 0; but more generally it is useful to consider the case where u does not necessarily vanish. The above metric has determinant −1, so √ −g = 1, Then the final form of the Hamiltonian is where φ and Π are canonically conjugated: [φ(x),Π(y)] = iδ(x − y). When u = 0, this is the Hamiltonian (8) in the main text.

Appendix B. Luttinger parameters v , K in the XXZ spin chain
The Luttinger description captures the low energy sector of interacting 1D systems, but the difficult task of determining the effective parameters v and K is left open. In general, one is forced to use some approximations or numerical analysis, except for an important set of models, namely integrable systems [60,61]. Their thermodynamics can be exactly computed in terms of proper integral equations, thanks to the Thermodynamic Bethe ansatz: then, the exact solution is matched with the Luttinger prediction, providing the expressions of v and K. In the spirit of the LDA, we can compute the local v(x) and K(x) as if the model was homogeneous.
Since in this work we tested our prediction on the XXZ spin chain, we restrict the discussion of the TBA and the subsequent extraction of v and K to this model. We closely follow [79].
We take for granted a basic knowledge of TBA, for which the reader can refer to [60,61], and just give the main formulas. Thus, we consider the Hamiltonian (26) for an homogeneous magnetic field h.
Notice that with a rotation along the x axis one can send Ŝ x →Ŝ x , Ŝ y → −Ŝ y and Ŝ z → −Ŝ z : thus we can change the sign h → −h. Therefore, we limit ourselves to study the case h > 0 and focus on |∆| < 1, which we parametrize as ∆ = cos(πγ). Due to integrability, the model can be understood in terms of stable quasiparticle excitations which are distributed accordingly to functions known as root densities. The description of the ground state requires a single root density, which satisfies the following linear integral equation The Fermi sea Λ is identified by the constraint that the dressed energy is equal to zero at the Fermi point ε(Λ) = 0. Above, the bare energy is E(λ) = −π sin(πγ)a(λ). We now consider the problem of fixing the Luttinger parameters: the velocity of sound appearing in the Luttinger is nothing else than the dressed velocity of the excitations at the edge of the fermi sea, namely where, for an arbitrary test function τ , the dressing operation is defined through the following integral equations Above, we need to use that ∂ λ p(λ) = 2πa(λ). In order to find K, we proceed as it follows. Conveniently, the Hamiltonian (26) can be regarded as a fermionic system with the correspondence Ŝ z → 1 2 −ĉ † jĉ j , being ĉ j standard spinless fermions [61]. Let n = ĉ † jĉ j be the density of Fermions, thus from Bethe Ansatz we have We define the compressibility as the system's response to external variations of the magnetic field κ = −∂ h n. In the fermionic language, adding a small magnetic field δh to the Hamiltonian couples to the fermion density, thus in the Luttinger liquid one gets a shift proportional to ∂ xφ One can then redefine ∂ xφ = ∂ xφ + K v δh and Ĥ is back to the Hamiltonian for h = 0, but expressed in terms of the new variables, thus ∂ xφ δh = ∂ xφ δh=0 . Therefore Thus, within the Luttinger Liquid theory we get We now compute the same object in the Bethe Ansatz formalism. Firstly, we need to study how the Fermi sea Λ is affected by a change in h. Thus, we regard the effective energy ε as an independent function of the rapidity, the magnatic field and Λ, i.e. ε(λ, Λ, h). Considering a small variation in the magnetic field h → h + δh and in the Fermi sea Λ → Λ + δΛ, we get where we neglect to write explicitly the Λ, h dependence, since there is no risk of ambiguities. From the above it follows We can now take the variation of the integral equation defining the root density (B.1) (B.14) It is convenient to define the function F(λ), which satisfies the following integral equation (we explicitly use that ρ(Λ) = ρ(−Λ)) Figure B1. Luttinger parameters and local magnetization in the XXZ spin chain obtained from the Bethe Ansatz solution as a function of the local magnetic field h. We focus on the cases ∆ = ±0.5. See also [28], where similar plots are shown for different interactions.
In this way δρ(λ) = δΛρ(Λ)F(λ). Finally, we compute the variation of the density Using the integral equation it is rather simple to show Above, using the definition of Z(λ) (B.11) and the symmetry of the kernel, one gets where we also used Z(λ) = Z(−λ). One then reaches the following expression for the compressibility which can be further simplified. Deriving the integral equation (B.3), using the symmetry of the kernel, integrating by parts and using that ε(±Λ) = 0 one gets from which we can conclude ∂ λ ε(λ) = (∂ λ E(λ)) dr . Using that 2πρ(λ) = (∂ λ p) dr and v = (∂ λ E(λ)) dr /(∂ λ p) dr , we can compactly write Comparing with the Luttinger compressibility (B.9), we finally have Of course, in the non interacting limit ∆ → 0 case there is no dressing and Z(Λ) = 1, i.e. the free fermion result as it should be. In figure B1 we plot the magnetization, the sound velocity v and the Luttinger parameter K as a function of the magnetic field h for the two values of D we consider.

Appendix C. Numerical methods for the discrete inhomogeneous Luttinger liquid
In this appendix we provide a numerical algorithm to compute the correlator φ (x)φ(y) on the ground state of the Hamiltonian (27). Similarly to what we did in section 2.2 in the continuum case, it is useful to define new bosonic fields It is then convenient to collect the operators in an unique vectorΨ Subsequently, the correlation matrix can be written as ΨΨ † , which is a 2M × 2M hermitian matrix. Since we are ultimately interested in the correlation functions of the φ j operators, it is convenient to express the Hamiltonian in this language as well. In particular, in this matrixvector notation we can rephrase Ĥ LL( ) as a bilinear form The matrix H has the following block form Which have now the correct normalization. With these normalized vectors we can construct the correlation matrix as from which the discretized two point correlator of φ j is easily derived, according to equation (C.1).

Appendix D. The constant C α
In this appendix we determine the constant C α (25), which has been overlooked in the original [50]. Equation (20) was computed by means of a direct computation of Trρ α I using path int egral methods: while doing so, constraints on the field configurations where imposed by means of Dirac's δ-functions. However, notice that replacing δ(. . .) → (const.) × δ(. . .) enforces the same constraint. This ambiguity in the normalization, after taking the logarithm, results in an undetermined constant offset in the Rényi entropies, which can be fixed by consistency arguments, as we are going to see. We start by quoting a result from [50] concerning the Rényi entropies for open boundary conditions for a bipartition Above, the definition of Φ ω is the same as in equation (24),the only difference being that we use the characteristic function of the interval attached to the boundary As already mentioned, the constant offset (which does not depend on the point x 2 , but it could depend on the integer α) is undetermined in the method of [50], but it can be computed from the alternative technique presented in section 2.2, which does not suffer from any ambiguity. A direct analytical comparison between equation (D.1) and section 2.2 is difficult. However one can check numerically, to at least 7 significant digits, that for arbitrary discretization step , the constant in equation (D.1) must be zero. In the following we assume that this identity holds exactly.
We now use equation (D.1) to fix C α in equation (20). This can be done noticing that, when considering the entanglement of the bipartition I ∪Ī with I = [x 1 , x 2 ], if we send x 1 → 0 the result (D.1) must be recovered. Comparing the two, one is immediately led to the conclusion Above, I a/α and M are defined in equation (21). In principle, I a/α (and thus M) still bear a non-trivial dependence on x 2 , but the resulting expression must be x 2 -independent.
One can also check numerically to high acccuracy that I a/α diverges in the x 1 → 0 limit. The matrix elements of M are given by Since the I −1 a/α also enters in the definition above, this means the matrix elements M ab go to zero also. Let us now check that (D.3) has a non trivial finite limit when x 1 → 0. To access it we just approximate the sum by an integral: