Entanglement entropy of two disjoint intervals separated by one spin in a chain of free fermion

We calculate the entanglement entropy of a non-contiguous subsystem of a chain of free fermions. The starting point is a formula suggested by Jin and Korepin, arXiv:1104.1004, for the reduced density of states of two disjoint intervals with lattice sites P = {1, 2, …, m} ∪ {2m + 1, 2m + 2, …, 3m}, which applies to this model. As a first step in the asymptotic analysis of this system, we consider its simplification to two disjoint intervals separated just by one site, and we rigorously calculate the mutual information between these two blocks and the rest of the chain. In order to compute the entropy we need to study the asymptotic behaviour of an inverse Toeplitz matrix with Fisher–Hartwig symbol using the the Riemann–Hilbert method.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Quantum systems that are spatially separated can share information that cannot be accounted for by the relativistic laws of classical physics. This fundamental property of quantum mechanics, which plays a crucial role in quantum information, is known as entanglement and its measurement is still largely an open problem. There is not a unique way of quantifying entanglement; however, in bipartite systems one of the most popular and successful measure of entanglement is the von Neumann entropy [9]. Suppose that the system is in a pure state |ψ . The density matrix is simply the projection operator ρ PQ = |ψ ψ|, where P and Q refer to the two parts and the Hilbert space H = H P ⊗ H Q . The von Neumann entropy is defined as where ρ P = Tr Q ρ PQ , ρ Q = Tr P ρ PQ (1.2) and Tr P and Tr Q denote the partial traces over the degrees of freedom of P and Q, respectively. In this paper we study the entropy of a two-block subsystem is a chain of free fermions. More precisely, we consider the chain The starting point for this analysis is an integral representation for the von Neumann entropy of the subsystem P of fermions on lattice sites P = {1, 2, . . . , m} ∪ {2m + 1, 2m + 2, . . . , 3m}. (1.5) This was derived by Jin and the fourth co-author in [26], and followed on from the success of this approach to calculating the entropy of a contiguous block of spins in the XX model [25]. Our goal is to compute the entanglement entropy between the subsystem (1.5) and the rest of the chain in the limit as m → ∞.
Over the past two decades the entanglement of bipartite systems have been extensively studied in one-dimensional quantum critical systems, in particular quantum spin chains. Consider a spin chain with N spins; at zero temperature the Hamiltonian is in the ground state and in the thermodynamic limit N → ∞ it undergoes a phase transition for some critical value of a parameter, e.g. the magnetic field. This quantum phase transition is characterized by an infinite spin-spin correlation length. Several papers have addressed the problem of computing the entanglement of the first L consecutive spins and the rest of the chain in various contexts [1,21,23,25,27,28,30,33,34,38]. It is well known that the entanglement entropy grows as S (ρ P ) ∝ log L, L → ∞. (1.6) Recently, there has been considerable interest in computing S (ρ P ) in quantum spin chains when P is made of disjoint regions of space. Up to now this problem has received attention within the framework conformal field theory (CFT) [2, 11-13, 17, 18]. One-dimensional quantum critical systems can be described in terms of a massless CFT. More general holographic descriptions are given in [32,35]. When P is one interval, then the coefficient of the logarithm in (1.6) is proportional to the central charge c, which is a characteristic of the theory [10]. If the theory is bosonic, i.e. if c is an integer, then in the two-interval case the von Neumann entropy depends on the compactification radius of the bosonic field [18]. In the papers [11,12] the moments of the density matrix were obtained for two-intervals as ratios of Jacobi theta-functions. Unfortunately, they could not compute the analytic continuation of their formulae in terms of the exponent of the moments, which would have led them to an expression for the von Neumann entropy, except in the asymptotic limit of small intervals [12]. A well established approach to solve quantum spin chains that goes back to Lieb et al [31] is to map the spin operators into Fermi operators using the Jordan-Wigner transformation. For example, the XX chain is mapped into (1.3). This approach works well when computing the von Neuman and Reny entropies of a single contiguous interval, as the entropy of the first L spins coincides with that of the first L fermions in (1.3). However, in the case of disjoint intervals in a spin chain there is the extra complication due to the fact that in the fermionic space the operators between blocks contribute to the entropy, because the Jordan-Wigner transformation is not local. This problem was tackled using CFT by Fagotti and Calabrese [17]. In order to avoid this technicality, our starting point is the fermionic chain (1.3). In the model (1.3), the Fermi operators in between blocks do not appear in the computation of the reduced density of states; therefore, the approach adopted in [26] applies. This simplification allows a rigorous computation of the asymptotic behaviour of the entanglement entropy as m → ∞ while at the same time preserving the physical phenomenon that we want to study. This idea is not new and was adopted by Ares et al [3], who performed a numerical study and conjectured a formula of the entropy of several disjoint blocks in a chain of Fermi operators. In fact, our main result-formula (2.12)-seems to be consistent with Ares-Esteve-Falceto conjecture. We hope to address this issue in all detail in the forthcoming publication. One of the main features of this representation of the von Neumann entropy derived in [26] for the two-blocks (1.5) is that the computation of the entanglement reduces to an integral involving the determinant of a block-matrix, whose two block-diagonal entries are Toeplitz determinants, see formulae (41), (48)-(51) in [26], or (2.1)-(2.4) below. This calculation would be the ultimate goal, but at the moment it is out of our reach-remark 1 in section 3. In this paper, instead, we consider a simplified example of a subsystem consisting of two intervals separated by just one lattice site. The asymptotic analysis of this model is already much more difficult than that of a single block Hamiltonian. Indeed, we not only have to evaluate the asymptotics of the Toeplitz determinant itself, but we also need to extract detailed information on the asymptotic behaviour of the inverse Toeplitz matrix.
It should also be noticed that, besides its intrinsic interest as a physical problem, the study of the asymptotics of Toeplitz determinants has a long history going back to Szegő [36,37] as such matrices are ubiquitous in mathematics and physics. Indeed, starting from the seminal works of Szegő, Kaufman and Onsager, the Toeplitz determinants have been playing a very important role in many areas of analysis and mathematical physics. Moreover, a growing interest has been recently developed to the study of certain generalizations of Toeplitz determinants. The most known among those are the determinants of Toeplitz plus Hankel matrices-see [8,16,19], the bordered Toeplitz determinants [4], and the integrable Fredholm determinants [15,20]. These determinants appear in the study of Ising model in the zig-zag layered half-plane [14], in the spectral analysis of the Hankel matrices, in the study of the next-to-diagonal correlation functions in the Ising model ( [4]), and in the theory of exactly solvable quantum models. In this paper, motivated by the physical model in the context of quantum information, we are concerned with yet another generalization of Toeplitz determinants, which are certain finite rank deformations of the standard Toeplitz matrices. In order to study such deformations, we need to analyse the asymptotic behaviour not only of the Toeplitz determinants per se but of the inverse Toeplitz matrices as well. The evaluation of the asymptotic behaviour of this new class of determinants which is done in this paper is, we believe, an important analytical result in its own right.
To summarize, in this article we compute the mutual information between a two blocks of Fermi operators separated by one lattice site and the rest of the chain in the Hamiltonian (1.3) explicitly. Our approach is based on the Riemann-Hilbert method, which has the additional advantage of being mathematically rigorous.

The main result
Let C denote the unit circle on the complex plane and The Fourier coefficients of g are In general, the m × m Toeplitz matrix and determinant with symbol φ ∈ L ∞ (C) will be denoted by T m [φ] and D m [φ], respectively. As it is well-known, the spectral norm (or operator norm) In particular, as T m [g] is a self-adjoint matrix, we obtain the relation σ(T m [g]) ⊆ [−1, 1] for its spectrum.
Let k, m, n ∈ N. We introduce the following matrix and determinant and Define the quantity The contour Γ ε goes around the [−1, 1] interval once in the positive direction avoiding the cuts (−∞, −1 − ε] ∪ [1 + ε, ∞) of e(1 + ε, ·), see figure 1. For instance Γ ε can be the circle (1 + 1 2 ε)C. For a general k, m, n we interpret the quantity in (2.4) as a measure of entanglement (kind of an entropy) between the subsystem (2.5) and the rest of the chain of free fermions (1.3) in the limit N → ∞. Here is our motivation for this interpretation.
Let H be an Hilbert space spanned by the fermions in the chain (1.3). Decompose H in the direct product H = H P ⊗ H Q , where H P is the space generated by the fermions b j at the lattice sites P indicated in (2.5). Write P = P 1 ∪ P 2 , where P 1 = {1, . . . , m} and P 2 = {m + k + 1, . . . , m + k + n} and denote by P 1 and P 2 the sizes of P 1 and P 2 , respectively. A standard calculation leads to the formula for the reduced density matrix. The angle brackets in this equation denote the expectation value with respect to the ground state. Applying Wickös theorem gives The above subsystem consists of two blocks/intervals of m and n fermions separated by a gap of length k. Using (2.7), it was shown in [26] that in the special case when k = m = n, and in the thermodynamical limit N → ∞, the quantity S(ρ P ) is indeed the von Neumann entropy of (2.5). We refer the reader to [26] for more details.
Our ultimate interest is to analyse S(ρ P ) as k, m, n → ∞, however, at this point the general problem seems to be far too complicated to attack directly (see remark 1 in section 3 below). Therefore we decided to start with the easier case when the gap between the two intervals is fixed to be k = 1, that is, when (2.5) becomes In this case the entries of A 12 in (2.3) become and, taking into account that g 0 = 0, As we shall see, this simplest case already leads to a mathematically very challenging problem. The asymptotic behaviour of the von Neumann entropy S(ρ (n) P ) of the interval {1, 2, . . . , n} was calculated in [25]. In particular, it was shown there that . Therefore the problem of calculating the limiting behaviour of the entropy of (2.8) reduces to the calculation of the mutual information between the two intervals: To analyse the asymptotic behaviour of this quantity as m, n → ∞ is still mathematically very complicated. However, as we expect this quantity to converge to a finite number, it makes sense to consider the following limit instead, where the ε and m, n limits are interchanged: We point out that a similar interchanged limit was considered in [22,25] for the case of one interval. The value of the limit (2.11) is what we shall calculate and interpret as the mutual information between the two intervals. It will turn out that indeed this is a finite number, which is stated in our main theorem. 1]). The limiting mutual information between the two intervals of the subsystem P from (2.8) is (2.12) The main tool in the proof of the above theorem will be an asymptotic analysis of an inner product involving the inverse Toeplitz matrix T m [φ] −1 . We phrase the related statement in the next section as lemma 3.2. We also notice that the asymptotics of the Toeplitz determinants D m [φ] and D n [φ] are described by the Fisher-Hartwig conjecture which, for the class of the symbols where the symbol φ belongs to, was proven by E Basor in [7].

Some preliminary calculations
We introduce the notations Therefore we obtain where we used standard facts about rank-one matrices and the following identity for blockmatrices: In particular, we infer Thus in order to compute the mutual information, we need to deal with the inner products G j , g j . It turns out that it is sufficient to handle the case j = 1.

Remark 1.
Notice that for a general gap of length k, the matrix λI − A whose determinant is needed to be evaluated, can be written as and the scalar coefficients γ and γ dl are certain k × k, independent of n, m, determinants. This shows what are the new technical challenges when one passes from k = 1 to the values k > 1.
The 'principal' determinant is not a block diagonal Toeplitz determinant anymore; indeed, the non-trivial off-diagonal Toeplitz blocks, generated by new symbols, appear. Moreover, the finite rank perturbation is of rank 2k and, therefore, ceases to be 'finite rank' as we consider the most general setting of the problem when all three sizes, m, n, and k become arbitrarily large.
From now on, until the end of section 8, our goal is to prove the following lemma, which then we shall apply in section 9 to prove theorem 2.1.
where the error term is uniform in λ on compact subsets of |λ| > 1.
In order to analyse G 1 , g 1 , we shall express it in terms of a Riemann-Hilbert problem (RHP) that arises in the theory of integrable operators, see [6, section 5.6] (see also [15], or [24]). Define the kernel This kernel defines a very special type of bounded singular integral operators on L 2 (C), namely a so-called (completely) integrable operator in the following way: where the integral is meant in the principal value sense, and we put the function in between []. By well-known properties of this operator, for all where (·, ·) denotes the complex inner product on L 2 (C). In particular, the connection between the two determinants can be shown by repeating the argument of [6, page 123]. In order to obtain (3.3) we observe that (by (5.157)-(5.158) in [6, page 123]) 1 − K has the block-matrix theorem 5.21] (see also [15], or [24]), where Y K is the unique solution of the following RHP. Y K -Riemann-Hilbert problem The unit circle is oriented in the usual positive direction, and the jump condition (3.6) is meant in the L 2 sense, see [6, definition 5.16].
In the next section we will connect the Y K -RHP with another RHP, but for the rest of this section our aim is to express the inner product G 1 , g 1 in terms of Y K .

Expressing the inner product in terms of the R-Riemann-Hilbert problem
Note that for all λ ∈ C\[−1, 1] the function φ possesses Fisher-Hartwig singularities at z 1 = i = e i π 2 and z 2 = −i = e i 3π 2 ; thus, we can apply the results in [16]. To be more precise, using the notation of (1.2) in [16], we can write φ in the following form: Note that throughout this paper, lnz denotes the principal branch of the logarithm, that is, −π < arg z < π. Since λ+1 λ−1 is a fractional linear map, we can easily examine the real-and imaginary parts of β. We have which is the reason why we shall take Γ ε = (1 + 1 2 ε)C in (2.12) in our calculations. Let us also note that β does not vanish on C.
Next, we shall connect the Y K -RHP with the Y-RHP, see e.g. [16] or [5] for details. Y-Riemann-Hilbert problem for orthogonal polynomials on the circle The jump condition (4.7) is meant in the sense that Y is continuous up to C from both sides, except at the points ±i.
It is well-known that this RHP has a unique solution which can be given in terms of orthogonal polynomials. An easy calculation shows the following connection between the unique solutions Y K and Y: where σ 3 = 1 0 0 −1 is the Pauli matrix. We point out that a similar connection was observed in [5]. Note that even though the jump conditions (3.6) and (4.7) are meant in different ways, one verifies easily that indeed the above Y K solves the Y K -RHP in the L 2 sense. The advantage of involving Y in our analysis is that we can use the powerful results of [16], in particular, we can express our inner product in terms of the R-RHP which can be estimated effectively. Let us recall the R-RHP next, whose associated contour Γ R is shown in figure 2. Notice that the circles ∂U 1 and ∂U 2 around ±i are oriented in the negative direction.

R-Riemann-Hilbert problem
14) The jump conditions (4.12)-(4.14) are meant in the sense that R is continuous up to Γ R from each side. The functions N and P j denote the global and local parametrices, respectively, see [16, subsections 4.1-4.2]. Namely, s−z ds stands for the Szegő function. The local parametrices will be discussed in detail in section 7.
From (4.10) we calculate If we trace back the transformations Y → T → S → R performed in [16], we obtain In particular, Notice that Y 12 (z) = O(z −m−1 ) as z → ∞, hence by (3.8) it does not contribute to our inner product. Therefore we have where D(z) = 1 + ∞ j=1 d j z − j (|z| > 1). Thus, from (3.8) we obtain , (4.17) where and notice that for all z ∈ C, Rz = 0, we have , (4.20) To summarise, we have two kinds of contributions to the inner product, one which comes from the Szegő function and another coming from R − I. Next, we compute the contribution coming from D(z).

The contribution from the Szegő function
Here we calculate the asymptotic behaviour of − 2 For that, we need a formula for the Szegő function. By (4.8) (or (4.10)) in [16], a short calculation gives where the right-hand side is analytic outside [−i, i], and Therefore, since D(1/z) = 1 + ∞ j=1 d j z j belongs to the Hardy class H 2 , we obtain the following expression for the limit of the sum: where we substituted u = tan ϑ and used standard residue calculus. We estimate the speed of convergence below.
Before we proceed with computing the contribution coming from R − I, we need some auxiliary calculations about the integral f m,1 defined in (4.19), and the local parametrices appearing in the analysis of R in [16].

Estimation of f m,1
We start with the following proposition. Proposition 6.1. We have Since ln(1 − u) is concave, we have at most two stationary points, and clearly one of them is u = 0. Simple calculations show that for u = 1 m the derivative is positive, and that for u = 3 m it is negative, therefore there is a second stationary point 1 m < u m < 3 m (m ∈ N, m > 5). It is then obvious that for all u ∈ [0, 1] we have We (a) Proof. (a) is obvious. Note that f m,1 is an odd function, therefore it is enough to prove (b) and (c) for t 0. For 0 t 1 − 1 m we have where I 1 (t) and I 2 (t) are the first and second integrals, respectively. By substituting y = ix and keeping in mind that m ∈ {2M + 1, 2M + 2}, we get for all 0 t 1 − 1 m that as m → ∞, which proves (ii).
as m → ∞, uniformly in t. What remains to prove is that the latter integral is as m → ∞, uniformly in t, which follows from the following calculations: Now, we estimate near i.
where the path for the second integral lies in ζ ∈ C\(−∞, 0], c |ζ| C, as shown in figure 3.

Proof.
Notice that for all c < |u| < C we have Therefore, using the substitution We note that one can similarly estimate near −i.

The local parametrices
In this section we shall compute how the local paramterices look like, with paying special attention to those parts that depend on m. As the two cases are very similar, we shall only examine the parametrix P 1 around i in detail. As in (4.12) and (4.23)-(4.24) in [16] we have and where ± = + when |z| < 1, and ± = − when |z| > 1. By equations (4.18)-(4.22) in [16], one easily sees that the auxiliary function F 1 (z) is constant in U 1 , and its value is The function E(z) is analytic in a neighbourhood of U 1 and is defined in (4.47)-(4.50) in [16]. What is important for our considerations is that where E(z) is independent of m and analytic in a neighbourhood of U 1 . Furthermore, and The function Ψ 1 (ζ) is an auxiliary function which is the main ingredient in constructing the local parametrix in [16], and which is given explicitly in terms of the confluent hypergeometric function ψ(a, c; z). We recall the details now. Let the contours Γ 1 , . . . , Γ 8 be defined as in figure 4. In particular, each of them is a half line starting or ending at 0. Furthermore, Γ k ∪ Γ k+4 is a line (k = 1, 2, 3, 4), and when k = 1, 3, these unions are the imaginary and real axes, respectively. These contours divide the complex plane into 8 sectors, denoted by I, II, . . . , VIII as shown in figure 4. The function Ψ 1 (ζ) is analytic in C\∪ 8 k=1 Γ k , and is uniquely defined by where the jump matrices J k are constant and are given by see (4.25)-(4.29) and (4.32) in [16]. Note that the functions ψ(a, c, ζ) and ψ(a, c, e −iπ ζ) are defined on the universal covering of the punctured plane ζ ∈ C\{0}, and that Ψ (I) 1 (ζ) is the analytic continuation of Ψ 1 | I to 0 < arg ζ < 2π.
which is uniform in λ on compact subsets of |λ| > 1.
Proof. We know from the standard analysis in the steepest descent method that Δ(z) = O(m 2|Rβ|−1 ) and hence R(z) as m → ∞, which is uniform in λ on compact subsets of |λ| > 1, and in z on C\Γ R (see e.g. [16], and note that there is some flexibility in choosing the parameters for Γ R , hence the Cauchy integral does not blow up as z gets closer to Γ R ). Also, elementary observations show that z−i which is uniform in λ on compact subsets of |λ| > 1, and in z on γ m . Therefore, combining the above estimates with lemmas 6.2 and 6.3, we conclude hence we obtain N(z) −1 = O(m |Rβ| ) as m → ∞, uniformly in z on C m i and in λ on compact subsets of |λ| > 1. Using (7.1)-(7.4) we also obtain the following: as m → ∞, uniformly in z and in λ on compact subsets of |λ| > 1. Hence Δ(z) = P 1 (z)N(z) −1 − I = O(m 2|Rβ| ), and therefore as m → ∞, uniformly in λ on compact subsets of |λ| > 1.
Next we show that, by the formulae in section 7, we have the following for z ∈ II ∪ III: In the above we have concentrated on the first column only. Denote the radius of U 1 by ε. Then, Finally, by substituting ξ = e −iπ ζ, dz z = − dξ m , and using the large ξ asymptotics of the confluent hypergeometric function and the usual estimates for R and D, we obtain As last step, we consider the integrals over Σ i,m j and Σ i,m j ( j = 1, 2).
Proof. First, we notice that the integral in (8.7) over Σ i,m j ( j = 1, 2) vanishes. Indeed, we have Δ + (z) − Δ − (z) = P 1,+ (z) − P 1,− (z) N(z) −1 and the jump of P 1 is exactly the same as that of S. Therefore, similarly as in (8.4) we obtain that R(z) Δ + (z) − Δ − (z) 11 In the rest of the proof, we shall only deal with the integral over Σ i,m 1 , and note that the other integral over Σ i,m 2 can be handled very similarly, since the local parametrix does not jump along Γ 7 . Note that for z ∈ Σ i,m 1 we have Thus, the integral over Σ i,m 1 has the following form: With the above proof we have finished proving lemma 3.2, which we use in the next section.

Calculating the mutual information
In this section we prove theorem 2.1, that is, we calculate (2.12). We start with a lemma.

Conclusions
Understanding entanglement in bipartite systems is of fundamental importance in quantum information, but at the same time it is often fraught with technical difficulties. One of the major reasons that makes it such a challenging problem is that even in simple systems the calculations are rather involved and often it is impossible to perform a rigorous analysis. Over the past twenty years a lot of research on bipartite entanglement has focussed on one-dimensional quantum lattice models, because they are amenable to a certain degree of mathematical manipulations. Originally, the interest concentrated on the entanglement entropy of a single interval of contiguous spins with the rest of the chain [1, 10, 18, 21-23, 25, 27-30, 33, 38] but more recently physicists have directed their research on the computation of the entanglement entropy of disjoint blocks [2, 3, 11-13, 17, 26, 32].
In this article we study a quadratic form of Fermi operators and compute the von Neumann entropy of two disjoint intervals separated by one lattice site. The major contribution of this paper is a rigorous analysis of the asymptotic limit of the entropy as the size of the intervals tends to infinity. More precisely, let P = P 1 ∪ P 2 , where P 1 = {1, . . . , m} and P 2 = {m + 2, . . . , m + n + 1}. Write S(ρ (m) P ) and S ρ (n) P for the entropies of the blocks of fermions at P 1 and P 2 , respectively; denote by S(ρ P ) the entanglement entropy between P and the rest of the chain. The quantity S(ρ (m) P ) was computed in [25]. We prove that the mutual entropy between P 1 and P 2 S(ρ (m) P ) + S ρ (n) P − S(ρ P ) → 2 ln 2 − 1, (10.1) as m, n → ∞. The proof is based on the Riemann-Hilbert method and involves computing the asymptotics of a Toeplitz determinant as well as extracting precise information on the asymptotic behaviour of the inverse of a Toeplitz matrix. Besides the intrinsic physical importance of the problem, the mathematical result is of interest in its own right. The asymptotic analysis of Toeplitz determinants has a long history and is still an area of active research. The ultimate goal would be to compute rigorously the entanglement entropy of two disjoint gaps separated by an arbitrary number of lattice sites. Unfortunately, this is still beyond our present ability.