Local and non-local properties of the entanglement Hamiltonian for two disjoint intervals

We consider free-fermion chains in the ground state and the entanglement Hamiltonian for a subsystem consisting of two separated intervals. In this case, one has a peculiar long-range hopping between the intervals in addition to the well-known and dominant short-range hopping. We show how the continuum expressions can be recovered from the lattice results for general filling and arbitrary intervals. We also discuss the closely related case of a single interval located at a certain distance from the end of a semi-infinite chain and the continuum limit for this problem. Finally, we show that for the double interval in the continuum a commuting operator exists which can be used to find the eigenstates.


I. INTRODUCTION
In entanglement studies, one divides a quantum system into two parts and determines how they are coupled in the chosen state. This information is encoded in the reduced density matrix ρ of one of the subsystems or, writing ρ = exp(−H)/Z, in the operator H. The latter has therefore been called the entanglement Hamiltonian, while in quantum field theory the name modular Hamiltonian is often used [1][2][3][4]. Its form depends on the quantum state in question as well as on the type of partition but also shows certain universal features and it has been the topic of various studies in recent years. For a review, see [5].
In one dimension, which we consider here, there are two standard partitions: a division of a chain (or line) into two halves or into a segment (or interval) and the remainder. In both cases, there is a simple analytical result for H if the chain is infinite and one is dealing with a continuous critical system in its ground state. Then where the integration runs over the subsystem A, T 00 (x) is the energy density in the physical Hamiltonian and β(x) is a weight function. For a half-infinite subsystem occupying the region x > 0, it varies linearly, β(x) = x, while for critical systems and an interval between x = a and arXiv:2204.03966v1 [cond-mat.stat-mech] 8 Apr 2022 x = b, it is given by the parabola β(x) = (x − a)(b − x)/(b − a). The first result is due to Bisognano and Wichmann [1,2] and the second one can be obtained from it by a conformal mapping, see [3,[6][7][8]. This expression has two characteristic features: H contains only local terms as the physical Hamiltonian and it is inhomogeneous. More precisely, the terms are large in the interior of the subsystem and small near its boundary where the factor β(x) vanishes linearly.
For discrete systems, the situation is somewhat more complicated. For free fermions on a lattice, one finds for the interval a dominant nearest-neighbour hopping in H which does not quite vary parabolically, but also hopping to more distant neighbours with smaller amplitudes [10,11]. Thus the entanglement Hamiltonian is not strictly local. However, it has been shown numerically [12] and also analytically [13] that in the continuum limit one recovers the conformal result for β(x) by properly including the longer-range terms. For free massless bosons in the form of coupled harmonic oscillators, the same was found through a numerical approach [14], and a similar analysis was carried out in higher dimensions for a spherical domain [15].
An intriguing new feature appears if the subsystem consists of two (or more) disjoint intervals.
This was discovered by Casini and Huerta in a study of massless Dirac fermions, i.e. free fermions in the continuum [16] and later investigated further in [12,[18][19][20][21]. One then finds a peculiar long-range coupling between each point in one interval and one single partner in the other interval.
Mathematically, it results from the form of the eigenfunctions of the correlation kernel from which one can construct the entanglement Hamiltonian. They show spatial oscillations exp(ipw(x)) in a variable w(x) which depends logarithmically on x, and as a result a function δ (w(x) − w(y)) appears in the expression for H which gives contributions not only for y = x but also for a point y = x c located in the other interval. The locus of the coupled points is a hyperbola in the (x, y) plane while the amplitude of the coupling varies with x and is given by a functionβ(x) resembling the β(x) for the local terms.
The situation on the lattice is again more complicated. To some extent, it was studied in [12] for slightly off-critical systems. Here we will investigate it in more detail and strictly for the critical case. In H one then finds hopping between a site and many others in the other interval. In the Hamiltonian matrix, this corresponds to characteristic regions of elements which are larger than those in the surrounding and do not scale with the size of the intervals. The shape of the regions is reminiscent of the hyperbolae in the continuum. We show that, in analogy with the short-range hopping, a continuum limit can be taken in which they combine and give the continuum result for β(x). This is done first for a half-filled system and symmetric intervals and then generalized to other fillings and unequal intervals.
We also treat a geometry which corresponds to some extent to one-half of the symmetric double interval, namely a single interval located at some distance from the boundary of a half-infinite chain. Its continuum version for Dirac fermions was studied recently in [22], and it turned out that a long-range coupling also exists there, but now inside the single interval. On the lattice, this is a certain drawback since now the usual longer-range terms and the new ones appear in the same region of the entanglement Hamiltonian matrix. Nevertheless, we were able to separate them by a proper choice of the summations in the continuum limit. Thus we could reobtain the continuum weight factorβ(x) also in this case, although only for half filling.
Finally, we consider an aspect which plays a considerable role in the treatment of single intervals in free-fermion chains. Namely, a simple operator exists in a number of cases which commutes with the correlation matrix (or kernel) and thus has the same eigenfunctions, see [23][24][25][26]. For the hopping model, this operator first led to the logarithmic oscillations of these functions [27] and later to analytical expressions for the matrix elements in H [11]. Here we show that such an operator also exists in the continuum case for both a single and a double interval. In the first case, it is a simple first-order differential operator, and in the second case it contains an additional difference term. A similar operator can also be found for an interval on the half line.
The layout of the paper is as follows. In section 2 we describe the setting and the models and give some known results for the double interval. In section 3 we consider the hopping model and present numerical results for the matrix elements in H for rather large intervals, obtained as usual from high-precision diagonalizations of the reduced correlation matrix. In section 4 we describe the continuum limit for the non-local terms in H and compare with the continuum result. In section 5, after presenting some numerical results, we do the same for a single interval in a half-infinite hopping chain. In section 6 we present the commuting operator. Finally, in section 7, we sum up our findings and give some outlook. In two appendices we construct the eigenfunctions of the commuting differential operator and derive its form for the interval on the half line.

II. SETTING
We describe here the two situations for which we consider the entanglement Hamiltonian.
The case of one-dimensional Dirac fermions with zero mass has been treated repeatedly in the past. One then is dealing with right-and left-moving particles described by field operators ψ R (x) and ψ L (x) and the Hamiltonian can be written with the energy density T 00 (x) given by In the ground state, all levels with negative energy are occupied and the correlation function for the right-movers is and similarly for the left-movers.
The entanglement Hamiltonian can then be obtained from the eigenfunctions of these integral kernels [17] as sketched below for the discrete case. For a double interval A = A 1 ∪ A 2 consisting of two disjoint intervals A 1 = (a 1 , b 1 ) and A 2 = (a 2 , b 2 ), it was found in [16] that with the local term as in (1) and the bi-local term given by where T bi-loc is the following operator The conjugate point x c in (6) is defined in terms of x as where and lies in A 2 if x is in A 1 and vice versa. It has a simple geometrical meaning, namely x c is the reflection of x on a circle with radius R around x 0 plus a reflection with respect to x 0 ∈ (b 1 , a 2 ).
In order to write down the weight functions β(x) andβ(x) occurring in the local term (1) and in the bi-local term (6), respectively, it is convenient to use the function which occurs in the solution of the eigenvalue problem and has the property w(x c ) = w(x), which actually defines x c . In terms of (10) one has In the special case of the symmetric configuration all formulae simplify. Then x 0 = 0 and R 2 = ab, the conjugate point becomes x c = −x wherẽ x = ab/x and the functions in (11) are andβ Note that β(x) is related to the corresponding quantities for the single intervals A 1 and A 2 via 1/β = 1/β 1 + 1/β 2 . In the following, we will compare the lattice results to these expressions and in the course of this also show graphs of them.

B. Lattice model
On the lattice, we will study an infinite fermionic hopping chain with Hamiltonian with fermionic creation/annihilation operators c † n and c n and chemical potential µ. Setting t = 1/2 and µ = cos q F , the ground state is a Fermi sea with occupied momenta q ∈ [−q F , q F ]. Similarly to the continuum case, we consider two disjoint segments A 1 and A 2 containing N 1 and N 2 lattice sites, respectively, and separated by a distance of D sites. The entanglement Hamiltonian for the combined subsystem A = A 1 ∪ A 2 is then a quadratic expression in the fermionic operators [9,10] where the matrix H is given by via the eigenvalues ζ k and eigenvectors φ k (i) of the reduced correlation matrix C A . This is composed of matrix elements C i,j = c † i c j given explicitly by with indices restricted to i, j ∈ A. Due to the two subintervals, C A and H have a block structure, which for the latter can be represented as The diagonal blocks H (1) and H (2) describe hopping within the corresponding segment, whereas the off-diagonal ones contain long-range hopping terms between the two segments and satisfy H (2,1) i,j = H (1,2) j,i . In the following, we present the results of numerical calculations of H for rather large subsystems with values of N 1 and N 2 up to 160. As is well known, this requires an extreme accuracy in the diagonalization procedure and amounts to working with up to several hundreds of decimal places in Mathematica.

III. ENTANGLEMENT HAMILTONIAN ON THE LATTICE
To get a first impression of the detailed structure of H, we visualize its matrix elements using color coded plots in Fig. 1. For simplicity, we consider a half-filled chain and equal intervals of size N 1 = N 2 = 100 separated by a distance 21, 51 and 101 respectively. On the left hand side of Fig. 1, a 2D-plot of the full matrix H is shown, where each dot corresponds to a matrix element and its magnitude is encoded by a color scale shown by the bars. The various blocks in (18) are separated by grey lines. As expected, the dominant (large negative) entries correspond to nearest-neighbour hopping in the diagonal blocks. These are accompanied by subdominant matrix elements corresponding to long-range hopping within each segment, reminiscent to the structure observed for the single interval case [11]. not further studied there. According to the continuum results of Sec. II A, the off-diagonal block should contain only single hopping terms between the conjugate sites, which lie on a hyperbola as indicated by the dashed black lines in Fig. 1. On the lattice one obtains a more complicated structure, however, one can still recognize the location of the hyperbolae from the color-coded plots. To better visualize the structure of the off-diagonal block, we also give 3D-plots on the right of Fig. 1, showing the absolute values of the matrix elements. Here the hyperbolic shape is even more evident, one has, however, an additional structure along the antidiagonal as the distance between the segments decreases. For large distances the hyperbola moves towards the diagonal and the structure becomes increasingly more peaked with decreasing amplitude.
One should emphasise that the dominant amplitudes in the diagonal blocks are two orders of magnitude larger than those in the off-diagonal ones. Furthermore, the former show an extensive scaling with N , as one would expect from the local weight function (12), which has a dimension of length. In contrast, the matrix elements in the off-diagonal block do not scale with N , as the corresponding bi-local weight (13) is dimensionless. However, in order to recover the functions β(x) andβ(x) from the lattice data, one needs some further considerations. Indeed, it was already argued in Ref. [12] that the longer-range hopping terms on the lattice should be included properly to recover the local weight function. Moreover, for a single interval the continuum limit can even be carried out analytically, recovering the CFT expression [13]. In the following we shall demonstrate how the continuum limit works for the double interval, both for the local as well as bi-local weights.

IV. CONTINUUM LIMIT
The continuum limit provides a relation between the entanglement Hamiltonian on the lattice and the one obtained for the Dirac fermion theory. It amounts to introduce a lattice spacing s → 0 and consider the thermodynamic limit N σ → ∞ such that N σ s = σ is fixed, with σ = 1, 2. Due to the block structure (18) for the double interval, one should distinguish between the diagonal and off-diagonal blocks of the entanglement Hamiltonian. As expected, the former one shall reproduce the local kinetic term (1), whereas the latter one corresponds to the bi-local contribution (6).
We first consider the diagonal blocks and proceed similarly as in the case of a single interval [13]. Namely, we rewrite the entanglement Hamiltonian as an inhomogeneous long-range hopping model where the hopping amplitude satisfies Note that the argument of the r-th neighbour amplitude corresponds to the midpoint of the sites involved and we assume it to vary slowly in the variable i. We then introduce the continuous coordinate x = i s, and apply the usual substitution for the fermion operators thus rewriting them in terms of right-and left-moving fields. The phase factors are needed to account for the quick oscillations on the lattice. In turn, the substitution (21) amounts to linearizing the dispersion and shifting the Fermi points to q = 0, thus reproducing the relativistic dispersion.
Carrying out the continuum limit is then rather straightforward. Substituting (21) into (19), Taylor expanding the fields as well as the hopping amplitude t r (i + r/2) → t r (x) + r s t r (x)/2 to first order in s, and dropping oscillatory terms, one arrives at the following expression [13] where T 00 (x) is given by (3) and N (x) is the number operator while the functions µ(x) and v(x) are defined as The result (22) is an inhomogeneous Dirac theory with velocity parameter v(x) and a chemical potential µ(x). Hence, contrary to naive expectations, the Fermi velocity in the continuum does not simply correspond to the nearest-neighbour hopping on the lattice, but is rather modified by the presence of long-range hopping. Furthermore, in order to get a finite result for the velocity v(x), the hopping terms t r (x) should scale with the size N σ of the segment, which is indeed the case as pointed out in sec. III. In contrast to the single interval case [11], however, we have no analytical results on t r (x) and the sums in (24) must be evaluated numerically. In particular, to recover the local piece H loc = H (1) + H (2) one needs to find v(x) = 2πβ(x) and µ(x) = 0.
could define the corresponding piece of the entanglement Hamiltonian as We would like to reproduce the bi-local term (6) by an appropriate continuum limit, we thus fix is = x and look for the contributions around the conjugate site js ≈ x c . Since the matrix elements in the off-diagonal block do not scale with the segment size, it is enough to keep the zeroth order term in the expansion of the fields. This is given by Note that in the above expression each term has an oscillatory factor with a large argument. In fact, it is not clear a priori, which one of these term provides a proper continuum limit, i.e. a smooth function of i. Since for the infinite chain the left-and right-moving fermions are not supposed to mix, we will consider the exponentials with the j − i factors and will provide further justification later. The bi-local piece of the entanglement Hamiltonian then reads where T bi-loc (x, x c ) is the operator defined in (7) and The corresponding weight functions with x i = i s are given by The same procedure can be carried out for the piece H (2,1) with indices i and j in (25) interchanged.
The continuum limit is the same as (27) with the integral running over A 2 , and the sums in (29) running over A 1 .

A. Equal intervals at half filling
We start with the simplest case of equal intervals N 1 = N 2 = N at half filling q F s = π/2, such that the distance between the intervals is 2d + 1. The coordinates of the endpoints in the continuum description must be chosen as a = ds and b = (d + N )s, such that b − a = . Due to particle-hole symmetry, H has a checkerboard structure with H i,i+2p = 0, which immediately yields µ(x) = 0, whereas v(x) has to be determined numerically. Due to the extensive scaling of the diagonal blocks, it is useful to introduce the rescaled hopping Using this in (24) where we introduced a cutoff P < N/2. Note that for a fixed i ∈ [d + 1, d + N − 1] the sums are carried out perpendicular to the main diagonal of the matrix.
The convergence of v(x i ) as a function of the cutoff is shown in Fig. 2 for N = 80 and d = 10, 20.
The case P = 0 is simply the scaled nearest-neighbour hopping and one observes that as the distance between the segments becomes smaller, it deviates more and more from the continuum limit prediction 2πβ(i/N ), shown by the red lines. Indeed, the data for P = 0 lies well above the red line with its maximum shifted towards the center of the chain, and a good overlap is found only close to the boundaries. Increasing the value of P one finds a slow convergence towards the continuum limit and already for P = 10 one has an almost perfect overlap.
In the above examples the ratio r = d/N between the distance and segment size is large enough to ensure a very good convergence already for moderate N . The situation becomes more between the segments is still too small to ensure a reasonable continuum limit. Nevertheless, Fig. 3 convincingly demonstrates that the limit is approached smoothly for N → ∞.
Finally, we consider the off-diagonal scaling functions in (29). Note that here the range of the At half filling the function C(x i ) = 0 vanishes identically, while S(x i ) simplifies to In contrast to the local term, this definition involves an alternating row-wise sum of the matrix elements. The result is shown in Fig. 4 for increasing segment sizes N and fixed ratios r = 1/5 and r = 1/10. In both cases one obtains a very good convergence towards the bi-local weight 2πβ(x), shown by the red lines. This function is more concentrated near the left end of the segment than 2πβ(x) due to the factor 1/x in (13), and the effect increases for smaller distances. Note that the relation S(−x i ) = − S(x i ) follows by symmetry arguments. One should also remark that keeping the terms in (26) with i + j would lead to another sum that is simply related to (32) by (−1) i S(x i ).
Obviously, this is an alternating expression which vanishes in the continuum limit, thus justifying our choice of dropping these terms.

B. Arbitrary filling
We now move to the case of arbitrary filling, choosing equal intervals for simplicity. For a single interval one can actually show analytically, that the continuum limit leads to the exact same expressions as in the half-filled case [13]. This proof uses the analytical expressions for t r (x), which are nonzero for both even and odd r [11], to evaluate the sums in (24). In the present case, where the sums are carried out numerically, this leads to some complications as the hopping profiles (20) on the lattice are associated to integer and half-integer sites for even and odd r, respectively.
Thus, analogously to the odd sums as in (31) for the half-filled case, one could define an even sum, albeit with coordinates assigned as x i = (i − 1/2)s. The even and odd sums can then be added by summing the i-th term of the odd sequence with the average of the i-th and (i + 1)-th terms of the even sequence. Alternatively, one could also define the row-wise sums with which also give only corrections of order s 2 that vanish in the continuum limit. Note that the factor of two as compared to (24) is now missing, since the sums run over an entire row of the matrix and not only in the upper diagonal part.
We observed that the row-wise sums yield a smoother convergence towards the expected weight functions. This is essentially due to the large oscillations as a function of i in the even/odd diagonal sums, which do not cancel perfectly when following the averaging process described above. The row-wise sums (33)  for the intervals of size N 1 and N 2 , separated by D sites. The matrix plot in Fig. 6 shows that the main features are similar to the symmetric case, and the hyperbola x c (x) can again be recognized.
However, the fine structure in the off-diagonal block differs from the one in Fig. 1, with a curved structure appearing also along the antidiagonal.
Clearly, there are now two scales in the problem and by taking the continuum limit one has to fix N σ s = σ , with the boundary coordinates given by a 1 = 0, b 1 = N 1 s, a 2 = (N 1 + D)s and b 2 = (N 1 + N 2 + D)s. Furthermore, the diagonal blocks also scale with the corresponding segment size and the densities can be introduced as However, the quantities β(x)/ σ andβ(x) are dimensionless and scale invariant, i.e. they remain unchanged under a rescaling of each length in the functions. This property can be used to fix either the length scale 1 = 1 for the left interval or 2 = 1 for the right one, leading to the definition The parameters of the corresponding β(x) function must be scaled accordingly, and can be shown to depend only on the two ratios ν = N 1 /N 2 and r = D/N 2 . In particular, for the right interval one has a 1 = 0, b 1 = ν, a 2 = ν + r and b 2 = ν + r + 1, whereas the parameters for the left interval read a 1 = 0, b 1 = 1, a 2 = 1 + r/ν and b 2 = 1 + (1 + r)/ν. The same arguments apply also for the bi-local weight, which is defined analogously to (32), with arguments x σ,i = (i − 1/2)/N σ . For the massless Dirac field, the most general boundary condition ensuring the global energy conservation corresponds to imposing the vanishing of the energy flow, T 10 (x = 0) = 0, through the boundary [28][29][30]. This boundary condition can be satisfied in two inequivalent ways, corresponding to a vector and an axial U (1) symmetry [31]. Here we discuss only the former case, which corresponds to the following boundary condition on the fields with a scattering phase α ∈ [0, 2π). Thus, the boundary condition at x = 0 provides a scale invariant coupling of the fields having different chirality.
The entanglement Hamiltonian of an interval A = (a, b) on the half line x 0 has been found in [22] and it takes the form (5). In particular, the local term is still given by (1) with the weight function (12) and the energy density T 00 defined in (3). Instead, the bi-local term in (5) is where the weight function coincides with (13), but the bi-local operator is defined as

B. Lattice results
Consider now the lattice model given by a semi-infinite fermionic hopping chain, whose Hamiltonian is (14) (17) then follow from the equations where i, j ∈ A 2 and i , j ∈ A 1 and the k-th eigenvector has been decomposed into two vectors corresponding to amplitudes on the left/right interval, respectively: Due to the reflection symmetry of the geometry, the eigenvectors are either even or odd, φ (1) k (−i) = ±φ (2) k (i). Using this property and introducing j = −j, i = −i, the eigenvalue equations (39) can be decomposed into two sets where the N × N correlation matrices C ± with elements i, j ∈ A 2 are defined as In fact, the piece containing i + j originates from the terms C i,−j and C −i,j in (39) after making the above substitutions, whereas C −i, −j = C i,j . Note that the above observations are independent of the distance d.
Remarkably, the correlation matrices C ± in (42) are nothing but the ones corresponding to a half-chain with Dirichlet or Neumann boundary conditions. They can be obtained by replacing the plane-wave basis of the infinite chain with appropriate standing-wave bases as Hence, the eigenvalue equation (41) simply tells us, that the even/odd eigenvectors and corresponding eigenvalues of the double interval problem are, up to normalization, identical to the ones for the half-chain with Dirichlet/Neumann boundary conditions.
The above observations can be used directly in the entanglement Hamiltonians H ± of the chains with corresponding boundary conditions. These are defined analogously to (16) as such that they are composed of only the even/odd eigenvectors φ ± k of the double interval problem, with the factor two resulting from the different normalization. Comparing now to the block matrix (18), one can see immediately that In other words, the entanglement Hamiltonians for the half-chain with Dirichlet/Neumann boundary conditions can be written as in terms of the diagonal and off-diagonal blocks in the symmetric double interval.
In Fig. 8  blocks for the double interval, as dictated by (45). In particular, one can see that the hyperbola is now reflected and appears atx = ab/x, as shown by the black dashed lines in Fig. 8.
The result (45) already makes it clear that the relation between the half-chain problem and the double interval is completely analogous to the continuum case in sec. V A. However, it remains to understand the difference between the bi-local operators (38) and (7), by considering the continuum limit of H ± . We first apply the substitutions (21), where the left/right-moving fields are coupled via the boundary condition (36) at x = 0. This coupling can be made explicit by writing the mode expansion of the fields ψ L (x) and ψ R (x), where the parameter α enters as a relative phase between the left/right propagating Fourier components. In our case here H + and H − correspond to the choice α = 0 and α = π, respectively. The change of the bi-local operator is due to the fact that the corresponding matrix component in (45) is mirrored. As pointed out already below (26), it is not a priori clear which phase factors give a proper continuum limit in this expansion. Due to the appearance of H (2,1) i,−j , it is now clear that for the half-chain one has to choose the terms with i + j. This leads to the following expression where the operators are defined as and the bi-local weights are given by Interchanging j → −j, one clearly arrives at the same definition of the weights as in (29) for the double interval. Hence, in the limit N → ∞, one has S(x i ) → 2πβ(x) as well as C(x i ) → 0, perfectly reproducing the result (37)- (38) for the boundary Dirac theory with α = 0, π.
The case 0 < α < π can also be studied on the lattice by introducing a chemical potential at the first site, i.e. adding the term µ 1 c † 1 c 1 to the Hamiltonian with µ 1 > 0. The eigenvectors are then cosine functions with a phase shift Φ(n) = 2 π cos[qn + δ(q)] , tan δ(q) = 1 + 2µ 1 cos q 2µ 1 sin q .
The parameter of the continuum theory then corresponds to the relative phase at the Fermi level, α = 2δ(q F ). It is easy to see that the correlation matrix has the form and the corresponding H α can be evaluated numerically. However, in contrast to (45), for generic α it is not clear how to separate the local and bi-local contributions that are superimposed in the matrix H α .  Here we propose such a method for the half-filled case. Inserting the full matrix into the definition of the bi-local weights, the local contribution in H α i,j does not yield a proper continuum limit, since the sum contains the phase factors with i + j. Indeed, one can argue that this term has an alternating form (−1) i f (x i ) with some smooth function f (x i ). In order to get rid of this unwanted contribution, we definē and similarly forC α (x i ). In other words, we combine the terms i and i±1, such that the alternating piece cancels up to second order corrections in the lattice spacing. In order to reproduce the continuum result, one expects to hold. For a generic scattering phase α, the checkerboard structure of the matrix C α in (51) as well as of H α is lost, and both of the sums in (52) are nonzero. This is demonstrated in Fig. 9, where we show the results for α = π/3 by choosing the value of µ 1 accordingly. The matrix plot on the left is more blurred due to the presence of nonzero elements H α i,j with even j − i, and on the right one can see a very good convergence towards (54) for both bi-local weight factors. We note that the local weights can be extracted analogously, inserting H α into (33) and then defining the averages which now kill the bi-local contributions. As expected, we findv α (x i ) → 2πβ(x i ) and µ α (x i ) → 0 independently of α.

VI. COMMUTING OPERATORS
In the following we consider the integral operators whose kernels are the correlation functions restricted to a subsystem A and discuss differential operators that commute with them.

A. Single interval
We consider a single interval A = (a, b) when the entire system is either on the line and in its ground state, or on the circle and in its ground state, or on the line and in a thermal state.
The integral operator associated with the correlation function G(x, y) restricted to A is and depends both on G(x, y) and on the spatial domain. Consider now the following first-order differential operator with a function t(x) still to be determined and evaluate the action of the operators (55) and (56) on a generic function f in the two possible orders, namely and with D A,x denoting the operator (56), while D A,y is (56) written in terms of the variabile y instead of x. In (58), an integration by parts has been performed with the assumption that t(x) vanishes at the endpoints of A, i.e.
One then sees that the operators (55) and (56) commute provided that Consider now the case of an infinite line with the system in its ground state. Then the kernel is given by (4) and the condition (61) becomes explicitly This functional equation is satisfied by an arbitrary quadratic function t(x). The boundary condition (59) then gives, up to a prefactor, For this choice of t(x), the differential operator (56) commutes with the integral kernel and thus has the same eigenfunctions. We note that t(x) is proportional to the quantity β(x) appearing in the entanglement Hamiltonian (1).
With this result, one can easily obtain the solution for a system on a circle with circumference L in its ground state. The kernel then is and one can check that the generalization of (64) to sine functions again satisfies (61). Finally, if the system is on a line but in a thermal state with inverse temperature β, one has and (61) holds if one substitutes hyperbolic sine functions in (64) The common eigenfunctions of I A and D A are easy to obtain from D A , since this is a first-order differential operator. In all three cases one has with eigenvalue is and w(x) given by For t(x) given by (64), w(x) is proportional to ln[(x − a)/(b − x)] and the exponential factor describes the logarithmic oscillations found previously in [12,27].
where ψ(x) corresponds to the field operator at τ = 0 and we dropped the index R, L at ψ(x).

The equation of motion then is
and inserting the expression for H given in (1) leads to the partial differential equation [16,22] 1 2π where D A is given by (56) with t(x) ≡ β(x) and ± corresponds to the right and left chirality, respectively. The expression on the right is the same as the one for finding the eigenfunctions φ(x) of the kernel of H which have to be the same as those of the correlation kernel. Therefore, with the knowledge of H, one could have found the commuting operator from (72). However, in analogy to other cases [23][24][25][26][27], where H is not known in closed form, our aim in Sec. VI A was to determine it directly from the correlation function.

C. Two disjoint intervals
It is possible to extend the considerations of Sec. VI A to the case where A is the union of two disjoint intervals on the line. When the entire system is in its ground state, the correlation kernel is again given by (62) and one has to find an operator which takes into account the non-local effects. This is somewhat difficult to guess, but one can invoke the equation of motion (72), since H is again known. In this way, one is lead to consider where the first piece is the operator D A introduced in (56) and the second one is an additional term where f is evaluated at x c , which is a function of x. For the moment, we leave the functions t(x),t(x) and x c (x) open, except for the condition that t(x) must now vanish at the endpoints of the two intervals. To obtain the commutator ofD A and I A , one only has to consider the additional term, which givesD and in reverse order where the integrand is the result of the change of integration variable y → y c (with y c (y) defined analogously to x c (x)) and followed by a renaming of the integration variable. If one now assumes this simplifies to which corresponds to (75) but with the action on y instead of x. Therefore, the additional terms can be combined with the others and the condition for the vanishing of the commutator is completely analogous to (61) It then remains to show that this relation holds for where both quantities are given explicitly in (12) and (13) and does not vanish because t(x) is not a quadratic function here, but it turns out that it is compensated exactly by the contributions fromt. ThusD A , which might be called a functional differential operator, commutes with the integral operator and can be used again to find the common eigenfunctions. This is sketched in Appendix A.
The commuting operator for the interval on the half line is discussed in Appendix B.

VII. DISCUSSION
We studied the entanglement Hamiltonian for two disjoint intervals in the hopping chain. The diagonal part of the matrix, containing the hopping within each segment, has a similar structure as for a single interval, with dominant nearest-neighbour terms. However, their profile deviates increasingly from the expected weight function β(x) as the distance between the intervals becomes smaller. In complete analogy to the single interval case, the continuum result could be obtained by We also considered the problem of a single interval at some distance from the boundary of a half-infinite chain. In the continuum, this setting is closely related to the case of two symmetric intervals, and we find analogous relations for an open chain. Nontrivial boundary conditions were implemented by a local chemical potential, and our continuum limit perfectly reproduced the bi-local terms of the Dirac theory that depend on the scattering phase.
The continuum limit obtained for the bi-local terms could be used in a number of other situations. In particular, one could apply it directly to the so-called negativity Hamiltonian [32], that is the analogous (albeit non-hermitian) object related to the partially time-reversed density matrix for two disjoint intervals [33]. The continuum weight functions then follow by applying a partial transposition as in [34] to the double interval result. One could also study bi-local terms in the entanglement Hamiltonian of inhomogeneous chains, the simplest example being the case of a defect [35]. For slowly varying inhomogeneities, the continuum entanglement Hamiltonian for a single interval can be obtained [36] via the curved-space CFT approach [37]. Recently it has been shown that the continuum limit can be generalized to recover the local terms for the domain-wall melting problem [38]. It would be interesting to check whether this treatment could be extended to the bi-local terms in case of two intervals. Further examples where non-local terms are expected even for a single interval include the case of finite temperature and volume [39][40][41], systems with zero-modes [42,43], finite mass [12,44,45] or systems driven out of equilibrium [46].
One can also ask, how well the exact lattice entanglement Hamiltonian can be approximated by a simpler one that does not contain all the long-range hopping terms. For single intervals, one can simply replace the nearest-neighbour hopping amplitudes by the properly discretized weights β(x i ), which gives a very accurate approximation of the exact reduced density matrix [47]. The question then is how such an approximation would work for two intervals and how to include the bi-local weightsβ(x i ). Studying these aspects, e.g. by variational methods as in [48], could shed some light on the role of these peculiar long-range couplings.
Somewhat related to this is the question of operators commuting with the correlation kernel (or correlation matrix), which we also considered. In the continuum, a simple first-order differential operator could be given for a single interval, whereas for the double interval an additional longrange coupling appeared. It would be quite interesting to know whether such a quantity also exists in the lattice problem.
The equation can then be integrated as such that one has Using the property as well as w(x c ) = w(x), it is easy to verify that the ratio indeed gives one Analogously, for the other ratio in (A5) one can also easily integrate (A7) with the result where the sign function ensures that the ratio φ − (x c )/φ − (x) is indeed negative. The functions m ± (x) are the same ones as found in [16,17,19] after orthogonalization.
The choice made in (A5) is the simplest one that leads to real-valued m ± (x). One can also see what happens if one assigns the ratios symmetrically Note that in order to have linearly independent solutions, we now need to work with complex functions m ± (x). The equation (A7) can again be integrated and for x > x 0 yields Using the identity i arctan z = 1 2 ln one arrives at This is clearly just a simple linear combination of the eigenfunctions found before. It is easy to check that it satisfies m + (x c )/m + (x) = i sign(x − x 0 ) as it should, and m − (x) = m * + (x).