Entropy and modular Hamiltonian for a free chiral scalar in two intervals

We calculate the analytic form of the vacuum modular Hamiltonian for a two interval region and the algebra of a current $j(x)=\partial \phi(x)$ corresponding to a chiral free scalar $\phi$ in $d=2$. We also compute explicitly the mutual information between the intervals. This model shows a failure of Haag duality for two intervals that translates into a loss of a symmetry property for the mutual information usually associated with modular invariance. Contrary to the case of a free massless fermion, the modular Hamiltonian turns out to be completely non local. The calculation is done diagonalizing the density matrix by computing the eigensystem of a correlator kernel operator. These eigenvectors are obtained by a novel method that involves solving an equivalent problem for an holomorphic function in the complex plane where multiplicative boundary conditions are imposed on the intervals. Using the same technique we also re-derive the free fermion modular Hamiltonian in a more transparent way.


Introduction
The reduced state to a local algebra of operators in a region in quantum field theory (QFT) can be expressed, in presence of an ultraviolet cutoff, as a density matrix ρ = e −H . The exponent H, called the modular Hamiltonian, conveniently encodes the reduced state. It retains a meaning in the continuum limit as the generator of unitary transformations corresponding to imaginary powers of the density matrix, ρ iτ = e −iHτ , the so called modular flow.
The structural importance of modular flows in the algebraic formulation of QFT has been being recognized since early times [1]. More recently, statistical properties of reduced states in QFT, have been the subject of much interest. In particular, entropy and relative entropy have simple geometric duals for holographic QFT [2,3]. In this context, the modular Hamiltonian in the boundary theory and its bulk dual have been used to elucidate localization properties of degrees of freedom in quantum gravity [4][5][6]. More generaly, knowledge of the modular Hamiltonian is an important step in relative entropy calculations, and it is fundamental to the formulation of entropy bounds [7][8][9][10][11] and the proof of several energy conditions [12][13][14][15][16].
However, most of our knowledge of the explicit form of modular Hamiltonians reduce to some examples where the modular flow is local, and it is primarily determined by spacetime symmetries [17][18][19]. In more generality, it is possible to identify a local part of the modular Hamiltonian which should have a large degree of universality [20,21], while not much is known about non local terms. An example of non local modular Hamiltonian which has been explicitly computed is for the vacuum state of the free massless fermion in d = 2 [22] (see also [23,24]). In this case H for several disjoint intervals has a local term proportional to the energy density and an additional non local part given by a quadratic expression in the fermion field. This last term however, does not contain all possible products of pairs ψ † (x)ψ(y) for fields located in arbitrary points x, y in the region, but only selected points appear: For each x in one interval only one specific "conjugate" point y appears in any of the other intervals.
There has been much progress in understanding local statistical properties of the vacuum reduced to two intervals in more general CFT. The Renyi entropies for integer index n have been explicitly computed for several models of interest [25][26][27][28][29][30]. However, analytic continuation to n → 1 to obtain the entanglement entropy has shown to be a difficult task. The computation of modular Hamiltonians of the vacuum for two intervals has also proved elusive [31].
A natural candidate to apply kernel methods is the free massless scalar. However, in two dimensions the uncompactified scalar field itself is ill defined due to infrared divergences, and one has to restrict the algebra to the derivatives of the scalar field. In this paper we compute the entanglement entropy and the modular Hamiltonian explicitly for two intervals in the theory of a chiral current, that is, the chiral derivative of the scalar field. The model is free and the modular Hamiltonian is quadratic. We diagonalize the kernel in this quadratic expression. The eigenvectors are obtained by an adaptation of the method in [20]. In the present case the problem is maped to one of finding analytic functions in the complex plane with specific multiplicative boundary conditions on the two intervals, located on the real line. In contrast to the case of the fermion, for the chiral scalar, the modular Hamiltonian is completely non local, while it still contains the expected local term proportional to the energy density operator [21].
The mutual information turns out to be a function of the cross ratio of the four end points of the intervals given by an integral over Hypergeometric functions. We check the result with numerical simulations on a lattice. For global pure states it is naturally expected that the entanglement entropy for complementary regions coincide. In this case the mutual information for two intervals would acquire an additional symmetry property relating the cross ratios η ∈ (0, 1) and (1−η) [32]. This property also follows from modular invariance in the replica trick calculation with Euclidean path integrals [33]. In the present model this symmetry is absent. This is explained by failure of Haag duality for two intervals: The algebra corresponding to the complement of the two intervals is smaller than the commutant of the algebra of the two intervals. Kawahigashi, Longo and Müger related this failure of Haag duality for two intervals in chiral conformal models to an algebraic index (µ-index) on inclusion of subalgebras [34]. This index also determines the amount of asymmetry in the mutual information [35]. For the chiral scalar the µ-index is infinity.
Our new method to deal with eigenvectors of the correlator kernels also leads to a better understanding of the derivation of the free fermion modular Hamiltonian. We start by a detailed derivation of this result in the next section. This will also allow us to introduce the main ideas to be used in the more complicated case of the chiral scalar. However, the treatment of the scalar case is self contained and a reader not interested in the discussion of the fermion field can start directly with section 3.
In section 3 we describe the algebra of the chiral scalar current, and the relevant kernels. In section 4 we show for a one interval region the kernel diagonalization, and the calculation of the entropy and modular Hamiltonian. The case of two intervals is dealt with in section 5, where we also compute the mutual information and modular Hamiltonian. Section 6 discuss the reasons for the breaking of the symmetry property of the mutual information. In section 7 we present the numerical calculations in a lattice and compare with the analytic results for the mutual information. Finally we end with a brief summary in section 8.

The massless free fermion revisited
A complete description of the reduced density matrix of a massless fermion field for multi-interval regions was given in [22]. This was achieved by diagonalyzing the correlator kernel in the region, using previous results on the literature about singular kernels of the Cauchy type [36]. Here we are doing this diagonalization transparent by mapping the problem of integral equations in one dimension to one of partial differential equations in two dimensions, following [20]. The form of the eigenvectors as well as its main properties are easily derived using this trick. This will also show how to generalize this calculation to scalars. We will treat the scalar field in the next section.
Since chiralities decouple in the massless case we consider a chiral fermion in d = 2, which only depends on a null coordinate. In our notation the variables x, y, etc, correspond to this null coordinate. We will consider a region A = (a 1 , b 1 ) ∪ (a 2 , b 2 ) ∪ . . . ∪ (a n , b n ) formed by n intervals. The field satisfies the anticommutation relations {ψ(x), ψ † (y)} = δ(x − y). The correlator kernel is 1 C(x − y) = 0|ψ(x)ψ † (y)|0 = 1 2 δ(x − y) + i 2π where the distribution on the left hand side is understood in principal value regularization. This is a projector when acting on the full line, and on the region A is an Hermitian operator with continuous eigenvalues in the range (0, 1) [22]. To obtain the modular Hamiltonian it is important to solve the spectrum of the correlator reduced to the region C A (x, y) ≡ C(x, y)| x,y∈A , because the modular Hamiltonian is given by [37][38][39] H =ˆA dx dy ψ † (x)H(x, y)ψ(y) , (2.2) where the kernel H is This last equation is understood as an operator equation, where the action of the operators is defined through their kernels.

An equivalent problem in the complex plane
In this section we will relate the orginal problem of solving the spectrum of C A as a kernel in A, to a new problem about a function in the complex plane. At the end, we will arrive at the same results of [20], but here we adapt the discussion to the chiral fermion field. For such purpose, we think the n intervals A as included in the real axis of the complex plane. For each λ ∈ R consider the following problem about a function S(z) in the complex plane S(z) analytic in C −Ā , (2.4) where l z,∂A is the distance from z to the boundary ∂A (formed by 2n disjoint points). Thus, S(z) has a cut over A with multiplicative boundary conditions. Consider now the complex integral where we choose an integration contour that encircles both A and z 1 in the positive (anticlockwise) direction. Then the integral vanishes because of (2.6), but writing it as two separated contributions from the pole at z 1 and the integration around the cut A we get where we have used the boundary condition (2.5). We remark there are no contributions from the end points of the intervals due to (2.7). This equation gives the value for S(z) on any point z ∈ C from its values at the cut A. Taking the limit z 1 → x ∈ A from above, and using we getˆA 11) which means that the boundary value of S(z) plays the role of an eigenvector with eigenvalue λ(λ − 1) −1 for the correlator kernel on A. Since the spectrum of C A (x − y) is restricted to (0, 1) (see [37]), we have that λ ∈ (−∞, 0). For later convenience we write λ = −e 2πs , s ∈ R . (2.12) Conversely, suppose we have a solution S + (x) of the eq. (2.11) for some λ ∈ R − with appropriate boundary conditions on the end points of the intervals as in (2.7). 2 Then equation (2.9) gives a complex valued function S(z) satisfying all the properties (2.4-2.7). For the boundary condition (2.5), the function S(z) defined in this way has boundary value S + (x) at the upper side of the cut, and for the lower side of the cut we have to use instead of (2.10), to get the right value S − (x) = −e −2πs S + (x).
In conclusion, the solutions of the problem in the complex plane (2.4-2.7) are in one to one correspondence with the eigenvectors of the correlator kernel (2.1).

Multiplicity and normalization of eigenvectors
Because of conditions (2.5) and (2.7), the function S(z) must have the following asymptotic behaviour when z → ∂A, where V ai , V bi ∈ C are constants. 3 Below, we will show these constants uniquely determine the solution.
In order to see this, for each s ∈ R we define the Green function G(z, w) for the problem (2.4-2.7), i.e. 18) and in addition G(z, w) satisfies the two boundary conditions (2.6) and (2.7) as a function of w for each z ∈ C fixed. 4 For w → ∂A then we have an expansions analogous to (2.14), Then, the combination G(z, w)S(w) does not have any jump singularity at A as a function of w. On the other hand, it has already simple poles at ∂A and at z, but it does not have a pole at infinity. Since the sum of all its residues must vanish, we have This shows there are at most 2n linearly independent solutions to the problem (2.4-2.7) for fixed s, and they can be viewed simply as elements of C 2n . It also shows that any solution which is bounded on ∂A (i.e. V ai = V bi = 0) must vanish. Now, we will show that the degenerancy of the space of solutions for each s fixed is indeed at most n. Let us take two solutions S 1 (z) and S 2 (z) corresponding to the same value s. The functionS 1 (z) = (S 1 (z * )) * is a solution with parameter −s instead of s. The functionS 1 (z)S 2 (z) does not have a cut, only poles at ∂A. The sum of residues must vanish and we get where V 1 ai , V 1 bi are the coefficients corresponding to S 1 and V 2 ai , V 2 bi the ones corresponding to S 2 . This means that any two solutions of (2.4-2.7) for the same s must be orthogonal according to the canonical (non positive) inner product of C n,n , which includes the case when the two solutions are the same. The argument to justify why the space of solutions is at least n is as follows. Suppose that the s-valued subspace of solutions is spanned by {S 1 , . . . S 2n }, where each S k is of the form (2.21). Then after a diagonalization procedure 5 we can get a new set of solutions {S 1 , . . .S 2n } which spans the same subspace but with the property that V k ai = 0 for all i = 1, . . . , n for all k = n+1, . . . , 2n. Automatically, because of (2.22), we must have V k bi = 0 for all i = 1, . . . , n and for all k = n + 1, . . . , 2n and henceS n+1 (z) = · · · =S 2n (z) = 0. In conclusion, the s-valued subspace of solutions has dimension at most n.
We will show in the next subsection that the dimension is exactly n. Now we make a final comment about the normalization of the eigenvectors S + (x, s), where we are writing explicitly the dependence of the eigenvectors through the eigenvalues s. Since any two eigenvectors S + 1 (x, s) and S + 2 (x, s ) must be orthogonal for s = s , we have 6 In order to orthonormalize the eigenvectors, we need to figure out the proportionality constant on the above equation. For this we note the delta function can only come from the integration around the end points of the intervals on the scalar product. Using the asymptotic expansion near these points we arrive at 7 Note the two terms inside the parenthesis in the right hand side are equal according to (2.22).

Construction of the eigenvectors
In this subsection we will explicitly construct the eigenvectors of the correlator C A (x − y) using the relation developed in subsection 2.1. Concretely, we will find the general structure of any solution S(z) of the problem (2.4-2.7) and through them we will obtain the corresponding eigenvectors. In particular, we will show that all eigenspaces for a given eigenvalue have dimension n. In this subsection s ∈ R is fixed. We start defining the complex valued functioñ where log is the principal determination of the complex logarithm which has a branch cut for z ∈ R ≤0 . The functionω is analytic everywhere on the plane except atĀ where it has a jump discontinuity of the form , the function f (z) = S(z) e −(is+ 1 2 )ω(z) (2.28) 5 Equivalent to the Gauss-Jordan algorithm used to diagonalize a finite dimensional matrix. 6 The function (S + 1 (x, s)) * is the complex conjugate of the boundary value of S + 1 (x, s) which is not the same that the boundary value of (S 1 (z, s)) * . These two operations do not commute. 7 More precisely, we should write each eigenvector as where the function R(z) has finite limit when z → ∂A. Then after replacing in the left hand side of (2.23) we get that the only possible delta Dirac contributions are of the form (2.24).
is analytic on C −Ā and it's also continuous on A, and hence 8 it's analytic on C − ∂A. Then f (z) must be some rational function with poles located at the end points of the intervals and also possibly at ∞. 9 Since S(z) satisfies (2.7) and because of (2.29) it follows that f (z) must be of the form with g(z) an entire analytic function. In order to satisfy the last requirement (2.6) for S(z), we have that g(z) can be a polynomial function of degree at most n − 1. Taking all this into account, any solution S(z) for the problem (2.4-2.7), must be of the form where a k ∈ C parametrize n linearly independent functions. Conversely, it's easy to see that any complex valued function of the form (2.31) is a solution for the problem (2.4-2.7). Taking the limit of z → A from the upper side of the cut on expression (2.31), we obtain the eigenvectors .
(2.33) Therefore, there are exactly n degenerate eigenfunctions for the same s. This space of eigenfunctions coincides with the one obtained in [22].

Scalar product
Due to the degeneracy, we have some arbitrariness for the election of the eigenvectors. Such freedom is encoded in the polynomial P (x) = n−1 k=0 a k x k of equation (2.31). In subsection 2.4 we will fix such freedom in order to get an orthonormal basis of eigenvectors. In order to do that, it is useful to have an expression for the scalar product between two eigenvectors in terms of its corresponding polynomials. In the rest of this subsection, we will obtain such expression. In eq. (2.24), using that the scalar product of two eigenfunctions is proportional to a delta function δ(s − s ), we obtained these scalar products in terms of the coefficients of the expansion of the eigenvectors near the end-points of the intervals. We will re-obtain this result here by explicit integration of the product of eigenfunctions.
First we take two solutions S + 1 (x, s) and S + 2 (x, s ) of the form (2.32) corresponding to two polynomials P 1 (x) and P 2 (x). Then, we compute the scalar product where we have changed the integration variable to ω and the sum in (2.34) runs over the disitinct solutions of the equation ω(x) = ω, which are the n simple roots of the polynomial equation In each interval A l = (a l , b l ), ω(x) is a monotone increasing function that goes from −∞ at a l to +∞ at b l . This fact implies that there exists n distinct simple roots x l , each of one belonging to any distinct interval A l . In (2.34) x l is understood as function of ω, i.e. x l (ω).
In order to proceed we will show that the following function of ω is a constant, i.e. K(ω) is independent of ω for any polynomial Q 2n−2 (x) of degree 2n − 2. Replacing the following expression for ω (x) in (2.36), we arrive at . (2.38) Since ω = −∞ implies x l = a l and ω = +∞ implies x l = b l , then we have the following particular limits . (2.40) Now, we will show K(ω) = K(−∞), and hence constant. For this, from equation (2.35) we have the following polynomial identity Evaluating (2.41) on x = a k (for some k = 1, · · · , n) we get and taking the derivative of (2.41) respect to x and evaluating at x = x l (for some l = 1, · · · , n) we have Then, replacing (2.42) on the denominator of (2.39) we get and replacing (2.43) on the denominator of (2.38) it follows . (2.45) Hence, the expected relation K(ω) = K(−∞) follows from where the last equality to zero is a general fact valid for any polynomial Q 2n−2 of degree 2n − 2: evaluating the polynomial in 2n arbitrary points y 1 , · · · , y 2n there is a linear equation that relates the value on the first 2n − 1 points to the value on y 2n . This equation is Eq. (2.46) follows specializing on y i = a i (for i = 1, · · · , n), and y i = x i−n (for i = n + 1, · · · , 2n).
Let us come back to the scalar product (2.34). Since the integrand on the Fourier transform in ω is constant, we getˆA

A complete orthonormal basis
In order to construct a basis of eigenvectors for each eigenspace of fixed s, we choose the following subset {u k s } n k=1 of eigenfunctions with the normalization factor 11 In the rest of this subsection, we will show that the set {u k s } n k=1 is orthonormal and complete. The orthonormality follows immediately from equation (2.49), and hence we havê The completeness is quite less obvious. The general argument of section 2.2 shows that n is the maximal degeneracy and then any n linearly independent vectors should form a complete basis. This fact can be shown explicitly as follows. 11 Note that the expression apparently differs from eq. (36) of [22]. But and then both equations are in agreement.
Using the eigenfunctions (2.50), we have where x l ≡ x l (ω(y)) ∈ A l are the n roots of the polynomial equation (2.35) for ω = ω(y). In particular, when y ∈ A l then x l ≡ y. Using the following algebraic relation 12 where the function is a polynomial in x of degree n − 1 for each fixed y, and its roots are the points x = x l except when x l = y.
Because of that, the only delta function which survives in (2.54) is for x = y and hence In order get a better expression for P (x, x), from (2.56) we define a new function which allows us to compute

Modular Hamiltonian
In this subsection we re-derive the results of [22] about the modular Hamiltonian using the information about the spectral decomposition of the correlator kernel C A (x − y) obtained in the previous subsections. The modular flow corresponding to this modular Hamiltonian and the entanglement entropy for several intervals have been computed in [22]. Recently, the modular Hamiltonian has also been computed using Euclidean path integral methods in [24]. In [23] it was shown that the modular flow satisfies the KMS condition. In [35] the mutual information between several intevals have been computed using the Araki formula without using a cutoff to compute the entanglement entropy. The results coincide with the ones in [22]. From (2.11) and (2.12) the correlator kernel writes where the function . (2.64) The aim for the rest of this subsection, is to obtain a simplified expression for the modular Hamiltonian kernel (2.63). First we have the following identity for the Dirac delta term

54). From this last equation, the modular Hamiltonian splits into the sum of a local and a non-local operators
where H loc (x, y) comes from the term in (2.65) when x l (y) = y and H noloc (x, y) comes from the n − 1 terms in (2.65) when x l (y) = y. We discussed these two contributions separately.

Local part
We recognize the local part for the modular Hamiltonian kernel as (2.67) In order to simplify the above expression, it is compulsory to understand it as a distribution acting over some smooth compactly supported test function ϕ(x, y). Integrating by parts, the derivative of the Dirac delta is converted to which can be simplified after we recognize the following identities which they follow from (2.64) and the algebraic relation 13 The final steps consist on replacing (2.69) and (2.70) into (2.68), and integrating by parts the term containing ∂ x ϕ(x, y), in order to factorize the test function. We get The local part of the modular Hamiltonian comes from (2.2) and (2.72) is the energy density operator.

Non local part
The non local part of the modular Hamiltonian kernel is (2.74) The first term can be simplified by a similar computation as we did around eq. (2.68). Here the situation is simpler because k(x l , y) ≡ 0 for all x l = y as we showed in (2.55). Hence the unique term which survives it is the one proportional to de derivative of k(x, y), Replacing it on (2.74) we arrive to and

Two intervals
For the case of two intervals A = (a 1 , b 1 ) ∪ (a 2 , b 2 ), the modular Hamiltonian operator H = H loc + H noloc reduces to where 3 The free chiral current We are going to study the model determined by a field operator j(x) in the line, which is identified with the chiral derivative of a massless free scalar in d = 2, that is, Then the line we are considering can be thought as a null line in the d = 2 model, and all variables x, y, etc. will be null variables. We will study the structure of the vacuum density matrix reduced to a region. Because of the complexity of the problem we are going to restrict attention to the case of one or two intervals. The commutator is and the model has Hamiltonian The model is Gaussian and all other multipoint correlators follow from this one according to Wick's theorem.
In order to proceed we will need general formulae for Gaussian states in an algebra of canonical commutation relations. We briefly review the derivation of these formulas in the appendix A.
An algebra of canonical commutation relations can be written in the form for hermitian variables f i , with numerical commutator given by the real antisymmetric matrix C ij . We take a Gaussian state with correlator The modular Hamiltonian is then given by where The entropy is [40] The operator V is not symmetric but it has real eigenvalues ν ∈ ±(1/2, ∞). In the present case this spectrum will be continuous. For later convenience we parametrize We name the left and right eigenvectors of V as where k is a possible degeneracy index. We normalize the eigenvectors as It is not difficult to see from (3.8), (3.11), (3.12), that C|u k s is an eigenvector of V † with the same eigenvalue as the vectors |v k s , and then is a linear combination of these later. The orthogonality (3.13) leave us the freedom to redefine the basis |u k s by an arbitrary matrix and the basis |v k s by the inverse adjoint matrix. We can use this freedom to set |v k s proportional to C |u k s . The phase of the proportionality constant is however fixed to be i sign(s), as can be seen from taking the scalar product of (3.11) with v k s | and using (3.13) and the positivity of F . As a result, we can further fix the vectors by taking In terms of these vectors the kernel (3.7) writes simply

The case of a single interval
Let us first consider the simplest case of a single interval A = (a, b). The inverse of the commutator δ (x − y) acting on a test function h(x) has to be a linear combination of´x a dy h(y), and´b a dy l(y) h(y), for some l(y). This last term is linear in h(y) and being independent of x, is annihilated by the δ . In kernel notation we have to combine 14 There is only one antisymmetric inverse, given by On a test function its action is Hence δ · (δ ) −1 = 1 and (δ ) −1 · δ = 1 on test functions that vanish on the boundary of the interval. Hence our first task is to solve the spectrum of .
We proceed as in the case of the fermion. We think the interval A as included in the real axis of the complex plane. We consider an analytic function S(z), with a multiplicative boundary condition on the interval A, as in (2.5). This is

5)
(4.6) 14 Only two of these three kernels are linearly independent.
Thus, S(z) has a cut over A with multiplicative boundary conditions. We further impose the boundary conditions at the infinity and at the end points of the interval, Consider now the complex integral on a closed curve in the complex plane encircling both A and z 1 . This integral is zero because the integrand has zero residue at infinity. We can shrink the integration contour around the point z 1 and around the cut to get The symbol ffl for the integral means here that it is regularized at the end points of the interval by taking the complex integral on a small circle around these end-points (as implied by (4.9)) and then take the limit of zero circle size. We will soon be more specific on how this regularization can be expressed directly for the function of a real variable such as S + (y).
Taking the limit of z 1 → x − i0 + , x ∈ A, and using (4.6), (4.4) and (3.8), we get Since the eigenvalues of |V | are in (1/2, ∞), we have λ > 0, in contrast to the case of the fermion where the factor λ is negative. We write λ = e −2πs , s ∈ R . (4.12) The eigenvalue in (4.11) then coincides with 1/2 coth(πs) as in (3.11). Therefore, for each solution S(z) of the problem in the complex plane we get an eigenvector of the kernel V (x, y) in the interval. The eigenvector modulus, in contrast to the case of the fermion, is bounded at the end points of the interval, (4.8). This is in accordance with boundary conditions for scalars [37]. Conversely, if we have a solution of (4.11) we can use it as boundary data on the interval in (4.10), which gives a solution S(z) satisfying (4.5-4.8). These problems are then mutually equivalent.
For a single interval we can write a solution to this problem as This obeys all boundary conditions. Any other solution divided by S(z) in (4.13) must be an analytic function on the plane except perhaps the end points of the interval. Consequently this other solution would be either proportional to (4.13) or it is not bounded at infinity or at the end points of the interval. Hence (4.13) is in fact the unique solution to the problem. The eigenvector is given by the boundary value on the interval (4.14) Now we explain more precisely the regularization in (4.10), (4.11). Frequently we will encounter integrals on the real line of the formˆb Then the integral will have an oscillatory but bounded term −c e −is log(x−a) /(−is) in the lower boundary as x → a, and it does not converge. The regularization used above just subtracts this oscillatory phase, that is If the oscillatory term appears on the upper end of the integral an analogous subtraction is understood. As mentioned above, this subtraction appears naturally when the integral comes from a limit of a complex integral around the cut as in the transformation from (4.9) to (4.10). The definition of the kernel V has to be understood with this regularization. 15 Now we have to look at the left eigenvectors of V . These are eigenvectors of V † = −iF C −1 − 1/2. For this we can just use the relation (3.14). However, we find instructive to compute them directly from the kernel. From (4.2) this is .
Now we take a new analytic function S(z), and assume the same multiplicative boundary condition (4.6). However, in order to obtain a solution of the eigenvector problem from the complex integral, we are now forced to impose S(z) → |z| −2 at infinity, and that S(z) must have at most pole singularities at a and b. We in a close contour encircling z 1 and the cut A. We get (4.20) The limit The value of λ is the same as above, giving the same multiplicative boundary conditions for S(z) as for the eigenvectors of V . However, the boundary conditions at infinity and at the end points of the interval are different. These now imply that the unique solution is The poles have to have opposite sign in order that the function decays at infinity as |z| −2 . We recognize this function is proportional to the derivative of (4.13), as it must be, given (3.14).
Orthonormalizing the eigenvectors we get These satisfy (3.13) and (3.14). 15 Note that this regularization eliminates from the bare integral an infinitely oscillatory phase in s, which will produce vanishing terms in any finite calculation involving integrals over the eigenvalues.

Modular Hamiltonian and entropy
Replacing these formulas for the eigenvectors into the equation (3.15) and after a simple integration, we get the following expression for the modular Hamiltonian kernel Then the modular Hamiltonian operator has the known form for an interval in a CFT [18] where in the present case the energy density is T = 1 2 j 2 . According to (3.9) the entropy is The full integral over the x coordinate gives a delta function δ(0) and is divergent. This is just the usual divergence of entropy in QFT due to the continuum spectrum of the modular operator.
A convenient way to regularize the entropy is to integrate up to a distance cutoff from the boundary. Then we compute with g 1 (s) = g(s). We get the well-known result (4.31)

The two interval case
To start, we need first to know the expression of the kernel C −1 ≡ δ −1 (x−y) for two intervals. The commutator is block diagonal in each of the intervals, and we get the same result as (4.2) for each of the intervals separately, or equivalently, Notice this last equation manages to be antisymmetric and its derivative is the delta function. Then we have and where Θ 1 (x) and Θ 2 (x) are the characteristic functions of the two intervals, that is, functions with value 1 inside the first interval (respectively second interval) and zero elsewhere. We have to deal now with kernels that contain theta functions and it might seem at first glance that the analytic method used in previous sections is not applicable here. However, we will show how to bypass this issue.
To begin with, let us consider eigenvectors v s (x) of (5.4) satisfying the extra property For such particular eigenvectors the second and third term in (5.4) vanishes. At the end we will show this is true in the general case. Under this assumption we have that v s is an eigenfunction of (x − y − i0 + ) −1 . Then we use the same ideas as for a single interval, trying to obtain v s as a boundary value of an analytic function. We again look for analytic functions S(z) on the complex plane with multiplicative boundary conditions on the two intervals A as in (4.6). The class of eigenfunctions u s of the problem must behave near the end points of the intervals as in the case of a single interval (or the case of the half line). That is, they should behave as pure phase factors of the form u s ∼ e −is log(x−ai) or u s ∼ e is log(bi−x) near the end-points. Their derivative, the v s functions, should at most have single poles (together with a phase factor) at the end-points. Under this condition the general solution is of the form Integrating S(z) along a contours encircling the two intervals and a large circle at infinity, it is not difficult to see that the integral at infinity is equal to the one over the two intervals, which vanish because of (5.5). Then this functions must fall as |z| −2 to cope with (5.5), and impose the condition Calling q = (a 1 , b 1 , a 2 , b 2 ), we have from (5.5) the coefficients α i satisfy in addition where and Only two of the eqs. (5.8), (5.9) and (5.10) are independent. This follows from the fact that This complex integral around the two cuts is zero because it is equal to the integral at infinity which vanishes because the integrand falls fast enough. Therefore we conclude the dimension of the space of these solutions for fixed s is two. The same argument of the previous section shows that these solutions give the eigenvectors of V † once evaluated on A, and any eigenvector of V † with at most simple poles at the end of the intervals, and satisfying (5.5), are of this form.
Now a simple solution is (5.14) In fact, this satisfies (5.8) and (5.9) because it is proportional to a derivative of the phase e −isω and hence the integral on any of the intervals vanish with the regularization we are using. That is, integrating this function on the intervals we have a further relation for the I l q integrals, The eigenvector v 1 (x) follows from taking the boundary limit ofṽ 1 (z) on A. The corresponding u 1 solution is an integral of this function, where in applying (5.2) to v 1 boundary terms that are oscillatory phases are discarded, in accordance with the regularization we are using. We can check more directly this is an eigenfunction of (5.3) by noting that This follows from (5.13), (5.15) and (5.16). Hence, the two terms with characteristic functions in (5.3) are equal and we can eliminate these functions altogether by replacing Θ 1 (x), Θ 2 (x) → 1/2. After this replacement the proof that u 1 is an eigenvector of the kernel (5.3) follows the same steps as the one for a single interval in the previous section by promoting u 1 toũ 1 (z) ∝ e −isω in the complex plane. We choose the second solution,ṽ 2 (z), of the form (5.6), such that its value on A gives a v 2 eigenfunction orthogonal to u 1 . Collecting the coefficients of the would-be delta functions in the scalar product between u 1 and v 2 , which are generated by the integral near the endpoints of the intervals, we have forṽ 2 From this and (5.8) we getṽ and from (5.9) it follows In order to compute u 2 we use the fact that the integrals of v 2 along the two intervals A 1 and A 2 vanish according to (5.5). Therefore using u = −isign(s) C −1 v, and (5.4), we have In order to normalize the solutions, we compute the coefficient of the delta function in the scalar product, that can only come from the singular part of the integrals near the end points of the intervals. We have that the normalized solutions satisfying (3.14) are We have not yet shown that these are the only possible eigenfunctions, since we imposed the additional conditions (5.5) to derive them. We will show at the end of the next subsection these conditions follow from orthogonality with u s 1 and u s 2 in the limit |s| → 0. To do that we first have to get simpler expressions for these eigenvectors.

Dependence of the eigenvectors through the cross ratio
The aim of this section is to obtain simplified expressions for the function α (5.21) and the eigenfunctions (5.23-5.26), which will be useful for the final computation of the modular Hamiltonian and the mutual information. For such purpose, we analyse the behaviour of such expressions in terms of the cross ratio which is the natural geometric parameter of the problem given the conformal invariance of the model. We consider a change of coordinates x = f (x) given by a general Mobiüs transformation where a, b, c, d ∈ R and ad − cb > 0. Such transformation leaves the cross ratio (5.27) invariant. Let us first understand the dependence of the function α(s; a 1 , b 1 , a 2 , b 2 ) with the interval endpoints. We can use (5.28) to make a change in the integration variables on the integrals (5.11) involved on the definition of the function α. After that, an straightforward computation shows that α(s; a 1 , b 1 , a 2 , b 2 ) = α(s; a 1 , b 1 , a 2 , b 2 ) , (5.29) where the two sets of intervals endpoints are related by Since such relation holds for any general Mobiüs transformation, we have that α(s, η) depends on the intervals endpoints only through the cross ratio. Similarly, a direct computation for the eigenfunctions reveals the following covariance properties under the change of variables x = f (x), where u i (x; q i ) and v i (x; q i ) are the eigenfunctions corresponding to the problem with endpoints q i = (a 1 , b 1 , a 2 , b 2 ) (idem for q i ) and the two sets of endpoints are related by (5.30). The real function K(q i ) is given by Not surprisingly, the u functions transform as a scalar wave, and the functions v as their derivatives under conformal transformations. Simpler expressions are obtained when we specially take the Mobiüs transformation f 1 (x) which sends the points (a 1 , b 1 , a 2 With this transformation we get from (5.21) the compact formula where 2 F 1 (a, b; c; x) is the Gaussian or ordinary hypergeometric function. To derive this result we used the integral representation for such function, where the function F 1 (a; β 1 , β 1 ; c; z 1 , z 2 ) in (5.40) is the Appell Hypergeometric function of two variables. Such function has the following integral representation for x, y < 1 and Re(a) > 0 and Re(c − a) > 0. 19 We remark formula (5.40) is only valid for x ∈ A 1 , and in such case x ∈ (0, η) and hence the arguments of the Appell's function belong to its domain of analyticity.
To get an expression for u 2 valid for x ∈ A 2 in terms of Appell functions we must consider a different Mobiüs transformationx := f 2 (x) which sends (a 1 , b 1 , a 2 Then we obtain for x ∈ A 2 , which is the same expression valid for u 2 in the first interval (up to a minus global sign) but evaluated iñ x instead of x . Such expression indicates that for any point x 1 ∈ A 1 exits a point x 2 ∈ A 2 such u 2 (x 1 ) = −u 2 (x 2 ), and viceversa. In the next subsection we show that all the eigenvectors are classified according to such "parity symmetry".

Parity symmetries of the eigenfunctions
In this section we study the behaviour of the eigenfunctions under a conformal transformation that interchanges the two intervals. For that we introduce the Mobiüs transformation where in this case we have K(a 1 , b 1 , a 2 , b 2 ) = 0. On the other hand, from (5.23-5.25) we easily see that Therefore, using additionally (5.40) and (5.45) for u 2 , we conclude we have the following parity symmetries Then, the first set of eigenfunctions is even and the second is odd under taking the conjugate pointx. We remark in these expressions we have the same eigenfunctions for the same endpoints (a 1 , b 1 , a 2 , b 2 ) appear at both sides. Another particular Mobiüs transformation which implies a quite simple symmetry relation, is the transformationx = q (x) which sends the end points (a 1 , b 1 , a 2 Under such transformation, the eigenvectors satisfy where now we must explicitly write the dependence of the eigenfunctions with the parameter s because the above expressions relate a eigenfunction of eigenvalue s with the one of eigenvalue −s. In this case we have a no null phase factor K (q i ) = 2 log b2−b1 a2−a1 .

Completeness of the eigenvector system
Before we compute the mutual information and modular Hamiltonian, we have to clarify if the eigenvector basis is complete. When we explicitly constructed the eigenvectors, we only consider solutions satisfying the equation (5.5) in order to simplify the calculation, but there was no further reason to assume that. Now, we are be able to show that any other possible eigenvector must satisfy (5.5), and hence there are no other eigenvectors than the ones already obtained. This fact follows considering the s = 0 solutions for u 1 and u 2 .
Taking out an irrelevant factor of |s| −1/2 (which is compensated by the inverse factor in the eigenfunctions v) we have The first one is proportional to a constant, the same in the two intervals, while the second one is proportional to two opposite constants in the two different intervals. Hence, any third solution v 3 for any s would be orthogonal to u 1 (s = 0) and u 2 (s = 0), and therefore must satisfy (5.5). Therefore there cannot be any other eigenvectors for two intervals.

Mutual information
In this section we compute the mutual information using the formulas developed in the previous sections. As we did for the one interval case, the equivalent expression of (4.27) to the present case is is the UV-regularized region and g(s) is given by (4.28). The entropy is UV-divergent, and is more convenient to express the result in terms of the mutual information independently. They must correspond to evaluate the integrals (4.27) and (5.65) along the regularized regions as we have already defined, with the same cutoff parameter for all the terms. After that, we take the limit → 0 + and we get the finite desired result for the mutual information.
Using the formulae (5.23-5.24), the first term on (5.65) can be easily calculated Direct treatment of the integral for the second eigenvector is more complicated due to the presence of hypergeometric and Appell functions in u 2 . We will find convenient to use the following trick instead. Since the integral in (5.65) is regularized, keeping fixed, we can replacê In the last step we have used the fact that vectors u s and v s+δs are orthogonal for δs = 0. The advantage is that we do not need now the precise behaviour of the eigenfunctions along the intervals but only in a small region near the end point of the intervals. Then we can just take the leading terms of u s and v s+δs since all other subleading terms in will disappear in the limit → 0. From (5.23-5.26) these leading terms are Plugging this back in (5.70) we get This implies where the first term in (5.74) coincides with the mutual information of the free chiral fermion field [22,41], and, taking into account that −i α(s, η)∂ s α * (s, η) = i∂ s log(α), and integrating by parts, the second term is given by .

(5.75)
We could not express this last integral in terms of standard known functions, and it has to be performed numerically. The result for U (η) is always negative, as it must be, considering that the chiral current is a subalgebra of the chiral fermion algebra, 20 and hence the mutual information has to be smaller (see formula (5.74)). In figure 1 we show a plot of the mutual information while the function U (η) is shown in figure 2. The mutual Renyi entropies I n (η) = S n (A 1 ) + S n (A 2 ) − S n (A 1 ∪ A 2 ) can be computed by replacing g(s) by g n (s), eq. (4.30), in (5.73). Hence we get .

(5.77)
In figure 3 we show I n (η) for some values of n. We will come back to discuss some aspects of these results in the following sections. In the rest of this section we will work out the asymptotic expansions for the mutual information for large and short distances between the intervals.

Asymptotic behaviour for I(η)
Before we continue with the calculation of the modular Hamiltonian, we want to analize the asymptotic behaviour of the function U (η) in the limits η → 0 and η → 1.
The η → 0 limit corresponds to the large distance limit between the intervals. Since the integrand of (5.75) is analytic at η = 0, a simple Taylor expansion reveals the following asymptotic behaviour This gives The first term coicides with the general result for the leading term of the large distance expansion of the mutual information for CFT [28,29,42,43]. For two intervals this is where ∆ is the lowest dimensional operator of the theory. In the present model this is j(x) = ∂φ(x) and has ∆ = 1. The fermion has ∆ = 1/2 and a different behavior I(η) ∼ (1/6) η at large distances, what is quite visible in figure 1. The short distance limit η → 1 is more tricky, since the integrand of (5.75) converges to zero in a non uniform way in such limit. The main contribution to U (η) in this limit comes from s ∼ 0. We have to expand the Hypergeometric functions at the numerator and denominator inside the logarithm in (5.75) for η ∼ 1 and s ∼ 0 to get (1 − η)) .

Modular Hamiltonian
We have now all the necessary elements to compute the modular Hamiltonian. This is H =ˆA dxˆA dy j(x)H(x, y)j(y) , (5.83) were, according to (4.25), the kernel becomes The kernel is real and symmetric. As the expression for this Hamiltonian will turn out to be quite complex, we need to start with some preliminaries about how we are going to express the results. Our first simplification will be to express the kernel in the four sectors A 1 × A 1 , A 1 × A 2 , A 2 × A 1 and A 2 × A 2 , that we call respectively H 11 , H 12 , H 21 , H 22 , in terms of the kernel in the first interval alone, H 11 , using the parity symmetry of the eigenfunctions.
Let us start by computing the contribution of the first eigenvector u 1 to H 11 , (5.85) where we set x, y ∈ A 1 and hence we have summed over only one of the roots of ω(x) = ω(y) in the delta function. We will find convenient to write that is, we have doubled the delta function contribution by u 1 , and the remaining part is, from (5.84) and (5.85), It will turn out that N (x, y) is a regular distribution. 21 Hence it gives the purely non local contribution to the modular Hamiltonian. Now, we recall the parity symmetry of the eigenfunctions studied in section 5.0.2 under the conformal transformationx = p(x) that interchanges the intervals, eq. (5.47). We have u 1 (x) = u 1 (x) and u 2 (x) = −u 2 (x). These relations give the following relations between the kernel of the modular Hamiltonian in the different sectors x ∈ A 1 , y ∈ A 2 , (5.88) These relations, together with (5.86), reduce the problem to the one of finding the form of the kernel N (x, y) in A 1 × A 1 . Unless otherwise stated, in the following we will assume x and y belong to the first interval. A different parity symmetry, (5.58), implies the kernel N (x, y), x, y ∈ A 1 , satisfies where the transformation x →x is given by (5.58). This follows from the corresponding symmetry of the eigenvectors (5.59), (5.60).
Another simplification is that we can relate all two interval cases with cross ratio η between the four endpoints of the intervals to the case where the two interval region is the standard region A η = (0, η) ∪ (1, ∞). This is done using the action of the conformal transformation x = f 1 (x) given by (5.28) on the eigenvectors, eq. (5.31). This simply gives where we wrote explicitly the dependence on the two interval regions. In the following we will call simply N (x , y ) to the kernel on the right hand side of (5.91), and keep x , y ∈ (0, η). That is, we focus on the first interval in the case where the region is A η . An evaluation of N (x , y ) requires the integration over s, which turns out to be the Fourier transform of products of Appell functions contained in u 2 , eq. (5.40). This obscures the structure of the kernel due to the complexity of these functions, and in particular the analysis of the possible singular terms. Instead we proceed in the following way. We first write the vectors u s as Then we make the integral in s that will be more amenable. Therefore we write We split the above kernel in two contributions K 1 and K 2 corresponding to v 1 and v 2 . Using (5.39) we get For the other term we use (5.41) and hence where z =ω (x ) −ω (y ), and and we have also introduced the function (5.98) and (5.99) are respectively the α-independent and α-dependent contributions to the kernel K 2 (x , y ). This gives the kernel N (x , y ) as a double integral over the sum of (5.95), (5.98) and (5.99). The final result depends on the Fourier transform of the function s 2 α (s, η), which has to be computed numerically. This numerical computation can be done after we have extracted the leading terms for s → ∞ from α(s, η). This will also help understanding the structure of singularities of these kernels. In the following we will make a further analysis of their local and non local parts.

Structure of singular terms
The asymptotic behaviour of the hypergeometric functions in α(s, η) for large s can be computed using the integral representation and the saddle point approximation. This is straightforward. The leading term was computed for example in [44]. Extending this calculation to include fluctuations around the saddle point we get the asymptotic expansion where Instead of extracting these asymptotic terms directly, we will write where nowα and the reminder function α r (s, η) is smooth in the parameter s and α r (s, η) ∼ 1 s 4 when |s| → ∞. The Fourier transform of s 2 α(s, η) iŝ whereα r (z, η) is the Fourier transform of s 2 α r (s, η). This is a continuous function vanishing exponentially fast at infinity. This real function is computed numerically, and ii is shown in figure (4) for some values of η.
Putting all together we finally get K (x , y ) = K 1 (x , y ) + K 2 (x , y ) = K 1 (x , y ) + K 2,α (x , y ) + K 2,/ α (x , y ) (5.111) = k δ (x , y ; η) δ (z) + k δ (x , y ; η) δ (z) + k i (x , y ; η) + k r (x , y ; η) , where where we have defined the positive polynomial function The singular terms can be simplified rewriting the Dirac delta distributions in terms of the variable x − y. A careful computation reveals the relations 22 Therefore, upon integration in (5.93) these terms behave as ordinary functions, and they produce a singularity in N (x , y ) which is just a jump in the first derivative for x = y . In a similar way, expanding for x ∼ y the term (5.114) reveals that it behaves as (x − y ) −2 for x → y . Upon integration, this gives a log |x − y | singularity for N (x , y ). Hence, this shows that the non local kernel N (x , y ) is given by a regular distribution.
The behaviour of the kernel K(x , y ) on the extremes of the interval, for example for x → 0 and y fixed, can be seen more clearly from (5.95), (5.98) and (5.99). The only non local contribution comes from the Fourier transform of s 2 α(η, s), Since α(η, s) is an infinite differentiable function of s its Fourier transform fall faster than any power of the variable. Therefore K(x , y ) is integrable on the boundary, and that is the reason we have not needed to use a regularization in (5.93). 23 As a result N (x , y ) fall to zero faster than any power of log(x ) −1 for x ∼ 0. A simplification in the structure of the non local kernel N (x , y ) arises if we take into account that the integration´A 1 dx v k s (x ) = 0. Hence, using (5.94) we could write (5.93) as ỹ)) , x > y . (5.121) In this way we avoid crossing thex =ỹ line in the integration, and therefore the delta functions (5.117) and (5.118) do not contribute. Moreover, the integrals are now completely regular and can be done numerically since we do not have to cross the singular points of the distribution. We checked these expressions coincide with (5.93). 24 A contour plot of N (x , y ) for η = 9/10 is shown in figure 5.

Summary of the modular Hamiltonian for two intervals
Summarizing the results, the modular Hamiltonian contains a local part and a non local part The local part, generated by the delta function in (5.86) and (5.89), gives a contribution to the modular Hamiltonian that writes on the full region A as where T (x) = 1/2 j 2 (x) is the energy density. The quantity β(x) = 2π (ω (x)) −1 acts as the local inverse temperature multiplying the energy density operator and controls the limit of relative entropy between the vacuum and energetic localized excitations around x [21]. This term, written in terms of the energy density, is equal to the local term of the modular Hamiltonian for the free fermion studied in section 2. This result coincides with general expectations for this term to be universal across two dimensional theories [20,21]. The non local part of the modular Hamiltonian is dx dy j(x) N (x,ȳ) j(y) .

Failure of Haag duality for two intervals
For simplicity, in this section we compactify the line and think in a circle S 1 instead. Let us consider two disjoint intervals A 1 , A 2 in the circle, and the complement is formed by two disjoint intervals, that we call A 3 , A 4 . We call A to the complement of a region A, hence (A 1 ∪ A 2 ) = A 3 ∪ A 4 . If the global state is pure one usually assumes This is equivalent to Taking into account the mutual information is a function of the cross ratio and that the entropies for single intervals are equal to the one of complementary intervals, we can express this relation for the model in the line as That is, the symmetry property for the entropy of complementary regions (6.1) gives the symmetry of the function U (η). This symmetry has also been shown as a consequence of modular invariance of two dimensional CFT's. We have seen this symmetry is not present for the free chiral current. This is not a problem of the continuum limit, the same happens in a finite lattice as we will show in the next section. In the rest of this section we explain how it is possible. The essential reason is that we have two basic choices for algebras of the two intervals. Let us call A to the commutant of the algebra A, that is, the set of all bounded operators that commute with all operators of A. If we take a set of operators S, the smallest algebra containing this set is S . This is the generated algebra by S. The first choice of algebra for two intervals is just the algebra generated by j(x) for x in the region, and we call it, Here we are smearing the field inside the region. 25 There is another algebra that can be attached to the two intervals, which consist of adding another operator to the generating set of the previous algebra. This is given by where f (x) is any smooth function such that The second algebra is generated by A (1) (A 1 ∪ A 2 ) and this new operator 26 It is evident that any two operators Φ 12 given by two different functions f (x) with the above properties will differ by an element of A (1) (A 1 ∪ A 2 ), and therefore the choice of f (x) satisfying these properties will not change the algebra A (2) (A 1 ∪ A 2 ). Note we can define A (2) (A 1 ∪ A 2 ) as the algebra of A 1 ∪ A 2 because Φ 12 actually commutes with all j(x) for x ∈ A 3 ∪ A 4 . This is precisely because f (x) is constant in A 3 ∪ A 4 , and the 25 There is an additional technical point in getting algebras of bounded operators, that can be thought as doing the spectral decomposition of the smeared fields and taking the algebra of the projectors, or, equivalently, taking the algebra of the exponentials of these operators. We are assuming this step when writing the generated algebra as in (6.5). 26 More correctly one should add the spectral projectors of this operator or the unitaries e iaΦ 12 for all a ∈ R. commutation relations (3.1). Notice these are two possible different algebra choices for the same underlying theory and for the same region. We call the mode Φ 12 a "long link" joining A 1 with A 2 , since being the integral of ∂ + φ, it is equivalent to some difference of the field φ localized inside the two intervals. However, in this model the field φ does not actually exist, and the only way to write this operator is through an integral of the current as in (6.6).
There is still the possibility of adding a long link crossing the interval A 4 , rather than A 3 as in Φ 12 . However, a sum of these two long links is (modulo operators in A (1) (A 1 ∪ A 2 )) equivalent to the integral of the current on the circle¸dx j(x). This last operator commutes with all the algebra and generates the global center. In a specific Hilbert space representation this operator will take a fixed value and then would be equivalent to a multiple of the identity operator. Hence the two options for the long link are actually equivalent. Now, we can define in analogous way the algebras . A long link crossing A 1 for example would be It follows from the commutation relations that In general, the algebras of complementary regions V andV commute, i.e. A(V ) ⊆ (A(V )) . If this inclusion is instead an equality, i.e. A(V ) = (A(V )) , it is said that the model satisfies Haag duality for the region V . Hence, in order to have Haag duality for a two interval region in this model one should choose the long link for one of the pairs of intervals and not for the complementary one. This prescription necessarily does not treat in an equivalent way all pair of intervals. Another perhaps more disturbing consequence of this choice is that the algebra A (2) containing the long link is not additive, meaning that it is not the generated algebra by the algebras of the two intervals. This is because the long link does not belong to the algebra generated by the single intervals. Hence we can have additivity at the expense of Haag duality, or viceversa, but not both properties together. The natural choice is the A (1) (A 1 ∪ A 2 ) because it can be consistently assigned to any two interval regions, and because of the additivity property, it is the only choice that allows the definition of mutual information.
The relation (6.1) holds for example when the complementary regions correspond to tensor product of full algebras in a lattice and the global state is pure. This situation will always give Haag duality. In the present case, failure of Haag duality indicates one can still think in terms of tensor products (in a regularized model) but where one of the factors does not have an interpretation in terms of the algebra of two intervals; it contains additional "long-link" operators. Hence, for the chiral current eq. (6.1) fails not because of the global state is not pure but because of failure of Haag duality, i.e., the commutant of the algebra of a region is not the algebra of the complementary region.
It is worth to notice that the algebra of the current j(x) is a subalgebra of the free massless chiral fermion. It is precisely the subalgebra generated by the fermion current ψ † ψ. The fermion is an extension of the current algebra, and it is one that satisfies Haag duality for two intervals. 28 This is why for the fermion U (η) = U (1−η) since trivially U (η) ≡ 0. It is by reducing the theory to the current algebra that we run into this particular trouble. The current subalgebras of the free fermion for two intervals are of the form A (1) , since the long link (6.6) measures the charge in A 3 and does not commute with the charged fermion field. This failure of Haag duality for two intervals in general chiral conformal models has been associated to an algebraic index (µ-index) of the inclusion of subalgebras A (1) This index should also determine the amount of asymmetry in the mutual information [35] as 29 In the present model we have seen this is divergent, in accordance with the fact that the µ-index of the current is infinity 30 .

Two chiralities and restoration of Haag duality
From the point of view of CFT in d = 2 it might seem strange that we could have Haag duality violation for two intervals and hence U (η) = U (1 − η). This property can be derived from modular invariance of the twist operators giving place to the Renyi entropies [33]. The reason is that Haag duality can be restored by adequately combining two chiral theories. Let us look at the example of the massless limit of a free massive scalar. The usual local algebra for a massive scalar in d = 2 is Haag dual and additive for two "diamonds", corresponding to two intervals in the spatial line at t = 0. At zero mass the zero mode of the scalar field has divergent mean quadratic variation and has to be removed. Then what remains are spatial and time derivatives of the field. With these we can form ∂ ± φ, the two chiral currents. For a single diamond, the algebra is then equivalent to the one of two decoupled chiral currents in an interval. For two diamonds however, the algebra also contains the difference φ(x 2 ) − φ(x 1 ), with x 1 and x 2 belonging to the two different diamonds. We can take x 1 , x 2 on the two intervals at t = 0. This is Hence, the two diamonds contain the difference of long link operators corresponding to the two chiral algebras. However, they do not contain the sum of these long link operators, and therefore the chiralities do not decouple for two diamonds. This is the reason these algebras for the two diamonds are compatible with the ones of the two complementary diamonds: the chiral long link operators do not commute to each other but their sum does, since commutators come out with opposite sign. Thinking in terms of the field differences (6.15) this commutation is evident. Therefore these algebras have the same form for all pairs of diamonds and Haag duality is retained. 31 However, without the zero mode, additivity is lost in this example. The massless limit of the mutual information of two intervals is divergent as I ∼ 1/2 log(− log(m)) [37].

The chiral current in the lattice
For doing numerical simulations we put the model in a lattice. We take the lattice Hamiltonian and the commutator Let us take a periodic system, and in order for C to be invertible, we take an even number N = 2n of points. The eigenvectors of the commutator are given by phase factors j C jl e ikl = 2i sin(k)e ikj , k = 2πm N , m = −(n − 1), · · · , n . (7.3) From here it follows that, defining the following variables for 0 < k < π (that is, m ∈ (1, n − 1)), they are canonical conjugates, [φ k , π k ] = i δ k,k . There are another two variables that form a global center of the algebra since they commute with all other elements, The inverse relation is The Hamiltonian then writes in these new variables This gives for the vacuum state φ 2 k = π 2 k = 1/2, φ k π k = i/2. The center can take any value and we set ψ 0 = ψ π = 0. Hence we impose these relations as a constraint. In this way we get a pure vacuum state and a global algebra without center. The full system has now n − 1 degrees of freedom: n − 1 coordinates and n − 1 momentum variables.
We see from (7.8) that we have two sets of low energy degrees of freedom, for k ∼ 0 and k ∼ π. Hence the system shows doubling of degree of freedom in the continuum, analogous to the usual fermion doubling. This is also the reason we have two commuting operators ψ 0 and ψ π .
The correlator of the original variables F (i − j) = f i f j is given by In the limit of a large circle N → ∞ we have The entropy of a region follows from (A.14), (A.16). We first check numerically the entropy for an interval. We calculate the matrices (7.2) and (7.11) for intervals of length R = 10k with k = 1, ..., 20. We fit the pairs (R k , S(R k )) with c 0 + c log log k + c −1 1 k + c −2 1 k 2 obtaining the logarithmic coefficient c log = 1/3 with high precision. Notice this coefficient is twice the expected one for the chiral current model. This reflects the doubling on the lattice.
To calculate the mutual information between two intervals of length a and b separated by a distance c, that is I(A, B) = S(A) + S(B) − S(A ∪ B), we need the entropies of the single intervals S(A) and S(B) and the entropy of the two intervals S(A ∪ B). Each of these entropies are calculated using (A.14, A.16). In the continuum limit the mutual information is a function of the cross ratio η, I(A, B) = I(η) where η is defined as , (7.13) in accordance to (5.27). For a given cross ratio, we repeat the calculation for different configurations that differ one from another just by a dilatation with parameter k = 2, 4, ..., 20. We then fit the pairs (k , I k (η)) with c 0 + c −1 1 k 3 and take the constant coefficient c 0 as the continuum limit of the mutual information for the lattice model, which is twice the chiral current model due to doubling. We then take I(η) = c 0 /2. We repeat the same procedure for different values of η obtaining the red points showed in figure  1.
In doing simulations for this model it is important that, if N is finite, we take the total number of points even N = 2n, and, in order not to have a center, the subsystems need to have an even number of points or variables (intervals of even size). This is because half of them are coordinates and half are momentum. The complementary subsystems automatically must have equal entropy because the global state is pure. For example, for an interval of size 2k in a circle of size 2n, the commutant is an interval of size 2n − 2k − 2, because there are two points in the complementary region adjacent to the interval that do not commute with the original interval. The entropies are indeed equal. When we consider two interval regions the commutant algebra contains a long link as explained in the previous section. In the lattice it will contain two long links operators, because of the doubling. More precisely, these commutant algebras with long links for two intervals are of the form: all points in the closed intervals [a 1 , b 1 ] ∪ [a 2 , b 2 ], and two long links, given by the sums f i and (−1) i f i , where the sums are over all the points in the open interval (b 1 , a 2 ) = (b 1 + 1, · · · , a 2 − 1). The long links crossing the other gap between the intervals are related to these by elements of the algebra of the intervals and the global constraints, and hence they do not give additional contributions. The counting of degrees of freedom is as follows: for a circle of 2n points, if the original intervals have 2k 1 and 2k 2 points, the commutant will have 2n − 2k 1 − 2k 2 − 4 points plus two long links. This gives a total of 2n − 2k 1 − 2k 2 − 2 linearly independent operators. This is precisely (twice) the complementary number of degrees of freedom: 2n − 2 is twice the total number of degrees of freedom in the lattice.
We have checked the entropies of complementary algebras of two intervals are equal in the circle. The entropy for the two intervals with the long linksS can also be completed to form a kind of "mutual information", eliminating UV divergences in the continuum, as I(A 1 , A 2 ) = S(A 1 ) + S(A 2 ) −S(A 1 ∪ A 2 ) . (7.14) The equality of the entropies for two intervals and the one of the complementary region including the long links, These two should be symmetric and equal for a model with Haag duality for two intervals but this is not the case in the present model. We have also checked numerically the relation (7.16) in the infinite lattice. For that we calculateĨ(η) and we note that convergence to the continuum limit is much improved for this case using the fitting function as c 0 + c −1/2 1 k 1/2 + c −3/2 1 k 3/2 , instead of integer powers, as we increase the global size k of the region. The continuum limit again correspond to the coefficientĨ(η) = c 0 /2.

Final Remarks
We have diagonalized the vacuum density matrix for a chiral scalar in two intervals. The modular Hamiltonian contains the usual local term given by an integral of the energy density times a position dependent inverse temperature. This term is identical to the free fermion one, and very probably is universal for all theories in d = 2. In addition, there is a non local term. This is given by a quadratic expression in the current with a locally integrable kernel which does not vanish in any open set of A × A. Hence the modular Hamiltonian is completely non local in contrast to the fermion case.
The mutual information does not have the symmetry property (6.3), and the origin of this is the failure of duality for two intervals.
We treated the case of two intervals. More intervals could in principle be treated in a similar fashion, but the expressions will depend on a higher number of cross rations, and besides the Hypergeometric and Appell functions that parametrize the eigenvectors should be replaced by higher dimensional Lauricella functions.
It would be interesting to understand why the fermion modular Hamiltoninan is quasilocal while the one of the current is completely non local. The technical reason is that one eigenvector of the bosonic model (u 2 ) has a dependence of the eigenvalue s that is not simply a phase factor e isω(x) . For the fermion field and any number of intervals, or the current field in the single interval case, this same phase factor determines completely the dependence of all eigenvectors in s. In this case the modular flow has a simple geometrical picture as a translation in the variable ω, ω → ω + 2πτ [22]. Perhaps a reason for the fermion to be special is the multilocal symmetries described by Rehren and Tedesco [46].
We have shown that the current mutual information is smaller than the fermion one because the former model is a subalgebra of the later. It would be interesting to explore other consequences of this inclusion. For example, the difference of modular Hamiltonians H ψ − H j between these models should be a positive operator. We can compute the expectation value of this difference of operators in a state generated from the vacuum by acting with a unitary in A, for example a coherent state e i´dx γ(x)j(x) |0 . The local contribution vanish in the difference of expectation values of the two modular Hamiltonians and we get an inequality involving exclusively the non local parts of H ψ and H j .

A Formulas for Gaussian state
Formulas for the entropy and modular Hamiltonian for Gaussian states in the algebra of canonical commutation relations in terms of coordinate and momentum correlators are described for example in [37]. We need here a slightly more general approach where we describe the algebra by 2n hermitian operators f i , i = 1 · · · 2n, with a general non degenerate numerical commutator where C is antisymmetric and real. This general case treated for example in [40,47,48]. The state is Gaussian with hermitian correlator and we have C ij = 2 Im(F ij ) . is a vector of field and momentum with canonical commutation relations. We write where X and P are the matrices of correlators of the field and momentum respectively. In writing (A.8) we are assuming there is no real part of the off diagonal blocks, and this is consequence of time-inversion invariance of the state, that is an anti-unitary symmetry mapping φ i → φ i , π i → −π i . In terms of the original variables f i it is an anti-unitary symmetry mapping linearly f i → T ij f j , with T real and T CT T = −C.
We can choose the system of eigenvectors of XP and P X, This was first shown in [40]. Analogously, the Renyi entropies defined by The modular Hamiltonian writes [37] H = Φ T g(P X)P 0 0 g(XP )X Φ , (A. 19) where g(y) = 1 2 √ y log √ y + 1/2 √ y − 1/2 . (A.20) In the present notation we have