{'fraction_non_alphanumeric': 0.036600805678112414, 'fraction_numerical': 0.01304431229618262, 'mean_word_length': 4.081984792357185, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Entanglement entropy of two disjoint blocks in XY chains

We study the Renyi entanglement entropies of two disjoint intervals in XY chains. We exploit the exact solution of the model in terms of free Majorana fermions and we show how to construct the reduced density matrix in the spin variables by taking properly into account the Jordan-Wigner string between the two blocks. From this we can evaluate any Renyi entropy of finite integer order. We study in details critical XX and Ising chains and we show that the asymptotic results for large blocks agree with recent conformal field theory predictions if corrections to the scaling are included in the analysis correctly. We also report results in the gapped phase and after a quantum quench.


Introduction
The bipartite entanglement for a given division of the Hilbert space into a part A and its complement B, can be measured in terms of the Rényi entropies [1] S (α) where ρ A = Tr B ρ is the reduced density matrix of the subsystem A, and ρ = |Ψ Ψ| is the density matrix of the whole system in a pure state |Ψ . The knowledge of the S (α) A for different α characterizes the full spectrum of non-zero eigenvalues of ρ A [2], and gives more information about the entanglement than the more commonly studied von Neumann entropy S (1) A . It also gives a fundamental insight into understanding the convergence and scaling of algorithms based on matrix product states [3].
In [4,5,6] it has been shown that for a one-dimensional critical system whose scaling limit is described by a conformal field theory (CFT), in the case where A is an interval of length in an infinite system, the asymptotic behavior of the quantities determining the Rényi entropies is where c is the central charge of the underlying CFT. Thus the Rényi entropies (and in particular the von Neumann one for α = 1) give one of the best way of detecting the value of the central charge. Less attention has been devoted until now to the entanglement of two disjoint intervals in a CFT (and also in massive theories). It turned out that entanglement of disjoint intervals is sensitive to universal details of the CFT that are not encoded in the central charge and it is connected with the full spectrum of operators of the CFT underlying the lattice model.
We consider here the case of two disjoint intervals A = A 1 ∪ A 2 = [u 1 , v 1 ] ∪ [u 2 , v 2 ] depicted in Fig. 1. By global conformal invariance, in the thermodynamic limit, Tr ρ α A can be always written as where x is the four-point ratio (for real u j and v j , x is real) Normalizing such that F α (0) = 1, we have that c α is the same non-universal constant appearing in Eq. (2). The universal function F α (x) depends explicitly on the full operator content of the theory and must be calculated case by case. Originally it was proposed that F α (x) = 1 identically [5]. This erroneous prediction has been tested in free fermion theories [7,8,9], and only the calculations in more complex theories allowed to detect this error [10,11,12,13]. Furukawa,Pasquier,and Shiraishi [11] calculated F 2 (x) for a free boson compactified on a circle of radius R where τ is pure-imaginary, and is related to x via x = [θ 2 (τ )/θ 3 (τ )] 4 . θ ν are Jacobi theta functions. η is a universal critical exponent dependent on the compactification radius R (in Luttinger liquid literature η = 1/(2K)). This has been extended to general integer α ≥ 2 in Ref. [12] F α (x) = Θ 0|ηΓ Θ 0|Γ/η where Γ is an (α − 1) × (α − 1) matrix obtained in [12], η is the same as above, while Θ is the Riemann-Siegel theta function The XX model considered in the following is described by the compactified boson with η = 1/2. We mention that the analytic continuation of this result to real α for general values of η and x (in order to obtain the entanglement entropy) is still an open problem, but the results for x 1 and η 1 are analytically known [12]. Nowadays, the only other example of non-trivial F α (x) known exactly is F 2 (x) for the Ising model that has the rather simple expression [13] . (8) We mention that all these new results for entanglement entropy were derived by using some old 'classical' results about CFT on orbifolds [14]. Figure 1. Typical bipartition we consider in this paper. The subset A is the union of two disjoint intervals A 1 and A 2 of length 1 and 2 respectively. The block separating them is denoted with B 1 of length r. The 'environment' is B = B 1 ∪ B 2 . The thermodynamic limit is obtained by sending the total length L → ∞, while 1 , 2 , r remain finite (i.e. the length of B 2 goes to ∞).
It is worth to recall that in the case of many intervals, the entanglement entropy measures only the entanglement of the disjoint intervals with the rest of the system. It is not a measure of the entanglement of one interval with respect another, that instead requires the introduction of more complicated quantities because A 1 ∪ A 2 is in a mixed state (see e.g. Refs. [15,16,17] for a discussion of this and examples).
In this manuscript we report the exact evaluation of the Rènyi entropies for integer α in XY spin chains with periodic boundary conditions, whose Hamiltonian is given by where σ α l are the Pauli matrices at the site l, h is the magnetic field and γ is the socalled anisotropy parameter. For γ = 1 the Hamiltonian reduces to the Ising model, while for γ = 0 to the XX model. Hamiltonian (9) can be diagonalized in terms of free-fermions (see below) but the fundamental difficulty is that the reduced density matrix of spin and fermion variables are not the same [13,18]. While calculating the fermion matrix is straightforward, we derive here a general method to deal with spin reduced density matrices and present actual calculations for the Hamiltonian (9).
The paper is organized as follows. In the next section 2 we give a general description of the density matrix mapping and we introduce our formalism. In section 3 we derive the expressions for the Rényi entropies of integer order. These two sections are rather technical and are intended only for those readers that want to repeat the calculations. The readers interested only in the results can skip directly to section 4 where we report our results for the XX chain and compare them with CFT results by properly detecting corrections to the scaling. In section 5 we present the same analysis for the Ising model and we check the universality while changing γ. In section 6 we consider non-critical chains. In section 7 we consider the Rényi entropies after a quantum quench. Finally in section 8 we summarize our main results and we discuss topics deserving further investigations. Four appendices contain technical details of the calculations and some background material.

The reduced density matrix
Generically a spin 1/2 chain can be described by means of the Pauli matrices σ µ i with µ = 0, 1, 2, 3 (σ 0 = 1 and 1 ≡ x, 2 ≡ y, 3 ≡ z) and i labels the lattice sites. We are interested in the entanglement of a subsystem A = N j=1 A j consisting of N disjoint spin blocks A j = [u j , v j ]. The reduced density matrix (RDM) can be written as a sum over all operators belonging to the blocks forming A where |A| = j j is the sum of the lengths of the blocks j . The expectation value is taken over the state we need to calculate the entanglement. The knowledge of multipoint correlation functions l∈A σ µ l l is however a very hard task even for eigenstates of integrable models (see e.g. [19,20] for the XXZ chains) and the above formula is particularly useful only when the spin-chain admits a representation in terms of free fermions, for which the Wick theorem suffices to calculate all of them.
The mapping between spin variables and fermion ones is obtained via the Jordan-Wigner transformation (that for the XY chains makes the Hamiltonians quadratic) It is also useful to define the Majorana fermions [21] satisfying anti-commutation relations {a l , a n } = 2δ ln . For the case of a single interval, a fundamental observation [21] is that the Jordan-Wigner string m<l σ z m maps the space of the first spins in the space of the first fermions and thus the RDMs of the first spins and fermions are the same, i.e.
with a 0 l = 1 and a 1 l = a l . At this point, it is worth to remind that the above equation gives directly all the eigenvalues of ρ A and so all Rènyi entropies (see e.g. the reviews [22] for more details). Indeed, one recognizes ρ A as an exponential form and using the Wick theorem (see also Appendix A) The 2 eigenvalues of the correlation matrix Γ have the form ±ν i with ν i ∈ [−1, 1], so one has all the 2 eigenvalues of ρ A as for the 2 possible choices of the variables s i = ±1. The Rènyi entropies for any complex α are then given by The calculation of ν i only requires the diagonalization of a 2 × 2 matrix, that is an extremely easy task compared to the diagonalization of the full 2 × 2 density matrix. This expression can be further analytically massaged to obtain exactly the non-universal constant c α [23] and also universal corrections to the scaling [24]. This is no longer true in the case of more disjoint spin-blocks. In the fermionic space, the subspace of two disjoint spin-blocks is not a local quantity: the non-locality of Jordan-Wigner transformation causes the fermions between the two blocks (i.e. those in B 1 ) to contribute to expectation values. For example, considering the product σ y l σ y n with l ∈ A 1 and n ∈ A 2 and calling r the distance between the blocks (see Fig.  1 that depends on the 'unpleasant' string of Majorana operators The string contribution cannot be neglected (as also shown in [13,18]) and so the spin representation is not equivalent to the fermionic one. For practical reasons we limit ourselves to the case of two blocks |A 1 | = 1 and |A 2 | = 2 . We denote with |B 1 | = r the space between the two blocks, see Fig. 1 for a graphical representation. We focus on spin-chains with the ground-state property that only even numbers of fermionic operators have non-vanishing mean values, as e.g. those with an Hamiltonian commuting with j σ z j . In this case, in Eq. (10) the product of operators σ µ l l does not depend on the string if we have an even number of Majorana operators (i.e. an even number of σ x,y ) in each block and it does depend if in each block we have an odd number of Majorana fermions (odd-even configurations give an odd total number of Majorana fermions that have vanishing expectation value by hypothesis). Thus, the density matrix splits in two sums extended only over even and odd number of Majorana fermions in each block, i.e.
where S is the complete string of σ z belonging to the region between the blocks (19). We repeat that the two sums are intended over all possible products of Majorana fermions belonging to each interval and even/odd refers to the number of µ l = 1, i.e. signaling the presence of a Majorana operator. We introduced also the 'short' O 1,2 for a general product of Majorana operators belonging to A 1,2 . Eq. (20) is an immediate consequence of the fermionic representation of spin variables and of the considerations above: the string S appears only in operators with odd number of σ x(y) in each block. Spin variables on different sites commute, and so [S, O 1,2 ] = 0. Thus the spin RDM can be written as where ρ ± are the fermionic RDMs We will be interested in the Rènyi entropies and in writing them in terms of the matrices ρ ± . The orthogonal projectors (1 ± S)/2 commute with ρ ± and so the n-th power of the spin RDM is the following combination of the n-th powers of fermionic RDMs We note that ρ + is unitary equivalent to ρ − , (the complete string S 1 (or S 2 ) of σ z belonging to A 1 (A 2 ) is the unitary operator mapping one into the other, i.e. S 1 ρ + S 1 = ρ − ). As a consequence we have where we used Tr S = 0. Still the matrices ρ ± are not of the form (14) that are the only objects we are able to deal with. Some further manipulations are still needed to bring them in a useful form. However, before continuing on the main road, we notice that ρ ± are the RDMs in fermionic variables of the state where |Ψ 0 is the ground state. These are the ground states of the Hamiltonians We remark this unitary transformation is nontrivial only for operators on the contact surface between the blocks and the remaining chain. Summarizing we can write a fermionic RDM equivalent to the spin RDM at the price of adding a finite number of non local terms to the Hamiltonian H. Now we write the fermionic RDM ρ ± in Eq. (23) as linear combinations of quadratic RDMs of the form (23) that are the ones we are able to deal with. The basic intuitive reason why this is indeed possible is that the σ z string is the exponential of a quadratic form and, in the case of a quadratic Hamiltonian, it can be interpreted as a term added to the exponential factor W nm in Eq. (14). These observations show that Wick theorem applies in some form for Eq. (23), but the 2-point correlations of ρ odd are different from the 2-point correlations of ρ even . The RDM's part ρ odd includes multipoint correlations in which the string S is inserted. To deal properly with this kind of objects we introduce the operator (that can be thought as a fake density matrix) where |Ψ Ψ| is the density matrix of the full system and the denominator is added to ensure the normalization Tr A ρ S A = 1. This fake density matrix could be introduced for any operator, but in the following we need only to consider the string S. If O is an operator acting on A, ρ S A satisfies the property Tr A [ρ S A O] = OS / S . If the operator S commutes with any operator with support in A, then ρ S A is hermitian. In Eq. (23) the sums are done only on even and odd numbers of Majorana operators in each block. We need to write these sums in terms of all Majorana fermions. A useful observation is that if we change sign to all the a l with l ∈ A 1 (or l ∈ A 2 ) in the multipoint correlators, then O 1 (O 2 ) changes sign in the odd sum, but not in the even one. By taking the appropriate combination of them we are left with a sum over all possible Majorana fermions and not only over the even/odd ones.
The operator that makes this useful change of sign is the Jordan-Wigner string S 1 (S 2 ) restricted to the first (second) interval, in fact using also S −1 1 = S 1 , we have Thus we finally arrive to (1 is the identity matrix) and so Notice that ρ 1 A and S 1 ρ 1 A S 1 have the same spectrum but different eigenvectors related by the matrix S 1 . The same is true for the matrices ρ S A and S 1 ρ S A S 1 , but they can have eigenvalues smaller than 0 and so in no way can be seen as true density matrices. In Appendix B we generalized this form to the case of more disjoint intervals.
Eq. (32) is the main result of this section: the rewriting of spin RDM as a linear combination of four fermionic RDM, i.e. exponential of a quadratic form as in Eq. (14). These four matrices do not commute and so they cannot be diagonalized simultaneously to find all eigenvalues of the spin RDM. However, if we are interested in Rènyi entropies with integer α, we can handle this problem in a constructive way: we determine the product rules between RDMs and then we construct recursively any finite order Rènyi entropy. Note that the RDM of the free-fermions is just given by ρ 1 A in Eq. (32). This matrix can be simply diagonalized in the same way explained above for the single interval as already done in Ref. [8].

Product rule
In this subsection we analyze the algebra of RDMs generated by a quadratic form, i.e.
where we do not assume W to be hermitian (to include the non-hermitian contribution of the string). Anti-commutation relations of Majorana operators make always W a complex skew-symmetric matrix, i.e. W T = −W . The constant ensures the normalization Tr ρ W = 1. In Appendix C we show that this normalization for complex diagonalizable skew-symmetric matrices is where {w}/ ± is the set of eigenvalues of W with halved degeneration (W is a skewsymmetric matrix so any even function of W has eigenvalues with even degeneracy). We assumed Z(W ) = 0, pathological cases can be cured as explained in Appendix D.
The product of fermionic RDMs of the form (33) is Indeed the commutator of operators in the exponent of (33) l,n,j,k is the essential ingredient in the Baker-Campbell-Hausdorff formula for the product of exponential of operators. Fermionic RDMs are specified by correlation matrices which can be written as (see Appendix C) Clearly the second equation is true only when 1 − Γ is an invertible matrix. At this point, let us briefly summarize the logic of the following derivation. We can easily calculate/manipulate the correlation matrix Γ that via Eq. (39) gives the exponential factor W that defines the quadratic density matrix in Eq. (33). We need to find what are the consequences of the product rule of RDMs for the correlation matrices, i.e. we need to find what is the correlation matrix corresponding to the product of two RDMs. While, through the chains of equations above, any W defines a single ρ A , the opposite is not true and there are several possible W 's for each ρ A . Nevertheless, we can give a unique recipe for the composition of correlation matrices. We indicate this matrix operation as Γ × Γ (notice it is not the product of the matrices) and it is formally defined by Eq. (36) as To specify this operation we still need two ingredients: (i) an usable expression for the correlation matrix associated to the product ρ W ρ W ≡ ρ[Γ]ρ[Γ ]; (ii) an expression for the trace of two fermionic RDMs in terms of the correlation matrices Γ and Γ .
The first requirement is easily obtained if we assume 1 − Γ and 1 − Γ invertible. Indeed, if we make explicit the exponential products after simple algebra we obtain Γ × Γ is a skew-symmetric matrix, even if it is not obvious from the above formula. We checked that this relation remains true also if 1 − Γ is not invertible (at least for the kind of matrices we are interested in), but a complete rigorous proof of Eq. (44) is beyond the goal of this paper.
The second request is less trivial because the correlation matrix Γ does not determine univocally the matrix W , and the sign of Z(W ) remains ambiguous. However, {Γ, Γ } is a functional of Γ and Γ , i.e. it is the product of the eigenvalues of (1 + ΓΓ )/2 with halved degeneration (the spectrum of ΓΓ is double degenerate [25]) The unspecified ± sign in front is the ambiguity that (as shown in Appendix C) can be solved by rewriting the composition rule as that does not depend on the curve γ.
It is evident that the operation × is associative and so we are in position to make any product of fermionic RDMs: where If, for some i and j, the matrix (1 + Γ i Γ j )/2 is not invertible, Eq. (47) cannot be applied. It is still possible to isolate the pathological parts and we report the details of this procedure in Appendix D.

Rényi entropies
From equations (25), (32) and (47), the Rènyi entropies for integer α can be written as where we defined the variables ζ i = 1, 2, 3, 4 (that label which of the terms in Eq. (49) is taken in the particular product) and we defined the shorts for the 2-point correlation matrices and In the case of α = 2 there are no problems with singular terms (see Appendix C) and the above expression can be re-written as where here and in the following equation we leave the sign ambiguity unspecified. Taking into account the trace's invariance under cyclic permutations, S 2 becomes the logarithm of a sum of 10 terms These formulas are already usable for a direct computation of Rènyi entropies. There are however some simplifications that occur by using the property of the correlation matrices Γ's. The first matrix Γ 1 is the standard fermionic correlation matrix (i.e. the one corresponding to free fermions in the absence of the Jordan-Wigner string, already considered in Ref. [8]), Γ 2 can be obtained from Γ 1 as Γ 2 = P 1 Γ 1 P 1 , with P 1 the Hermitian unitary matrix The same relation occurs between the third and the forth matrix Γ 4 = P 1 Γ 3 P 1 . Instead Γ 3 is not trivially related to Γ 1 . In Appendix A we prove the following identity where the double subscripts take into account restrictions to rectangular correlation matrices, i.e. the first (second) subscript identifies the region where the row (column) index runs. We show in Appendix B that similar properties are valid also for an arbitrary number of intervals.
Using these relations after some algebraic manipulations one can write down the full sums for the Rènyi entropies. It is important to notice that S = Pf(Γ B1 ) and so S 2 = det(Γ B1 ). Furthermore to short the notations we denote Finally S 2 can be written in a rather simple way But, increasing the order α, the explicit expressions become soon long. For example, we give the (simplified) formulae for S 3 and S 4

Critical XX model
In this section we report the explicit results for the XX chain in zero magnetic field (i.e. γ = 0 and h = 0 in Eq. (9)). XX chains have been previously analyzed by Furukawa et al [11] by means of exact diagonalization techniques that allowed to explore relatively small chains with at most 30 spins. The asymptotic results from CFT (cf. Eq. (5)) were obscured in this previous analysis by large oscillating corrections to the scaling. The smallness of the systems and the lack of a precise knowledge of the form of the correction to the scaling made impossible any finite-size scaling analysis to check the CFT predictions (5). In fact these oscillating corrections to the scaling have been widely observed in last years [26,27,28], but a precise theory about their origin and their exact form is only recently available [24,29]. By exploiting the exact solution, we can avoid these problems and explore large enough values of allowing a finite-scaling analysis similar to the one for a single block [24]. We start from the infinite volume limit. We consider two intervals both of length = 1 = 2 at distance r. In this case, the four-point ratio x in Eq. (4) is given by and the scaling of the quantity determining the Rényi entropies is (we use c = 1) that is an implicit definition of F lat α (x, ) encoding the fact that for finite we expect corrections to the scaling of the form where F CFT α (x) is the universal quantity appearing in Eq. (6). The exponent δ α governs the leading correction to the scaling. It has been shown with CFT in Ref. [29] that this exponent is equal to 2x/α independently of the number of intervals and x is the scaling dimension of an operator introduced by the conical singularity necessary to describe the reduced density matrix. For a single interval, it has been shown by exact analytic calculation that δ α = 2/α [24] and thus we expect the same exponent for the double interval case (the presence of a relevant operator with dimension x = 1 has been justified in Ref. [29]). In the following analysis we will greatly benefit of the exact knowledge of the quantity c α from the exact solution of the single interval entanglement [23,30] and so our results will be also a non-trivial check of the equivalence of the two non-universal multiplicative constants in the single and double interval case.
We report in Fig. 2 results for the function F lat 2,3,4 (x, ) for various and x. It is evident that irrespective of the value of x, with increasing the results approach the CFT prediction. For odd the asymptotic result is approached from below, while for even ones it is approached from above. These are the already mentioned oscillations that made difficult the analysis based on small chains. These plots are the first direct quantitative tests of the correctness of Eq. (6). These figures do not leave doubts about the correctness of the results of Ref. [12]. It is worth to mention that, as already shown elsewhere [11,12,13], the finite curves do not have the symmetry x → 1 − x that is restored only in the → ∞ limit Eq. (6).
Having established the correctness of the asymptotic form, we can now move to the finite corrections and check whether the prediction δ α = 2/α [29] is correct.  Having precise control of the corrections to the scaling is not an academic task: their analysis is fundamental to provide accurate results when such large system sizes are not available (in previous studies they would have been an important tool) and in cases when the asymptotic form is not known (as in the following section for the Ising model for α > 2). In Fig. 3 we report the function |F α (x) − F lat α (x)| for fixed x = 1/9 and α = 2, 3, 4. We both report results for the spin Rènyi entropies and for the fermionic ones (i.e. without considering the string contribution). We recall that for free fermions we have F α (x) = 1 identically. The results show a power law behavior for large enough with the predicted exponent 2/α as for the single interval. Notice that by increasing α, the values where the leading asymptotic correction can be identified become larger and larger, in analogy with the single block case [24] (as obvious because of the smallness of the exponent δ α ). For the fermionic variables the asymptotic behavior is reached before than for spin degrees of freedom: the string introduces further corrections to the scaling that in the present model are subleading. To show the x independence of this exponent in the last panel of Fig. 3, we report the same kind of plots for different values of x, showing that, at fixed α, the corrections lie on parallel lines.
To conclude this section we present some results for finite systems. It is known that all the formulas above (including the x dependence) at the leading order in finite systems can be described by replacing any distance u ij by the chord length u ij → L π sin πuij L . We checked that this rescaling is indeed correct and that, for large enough and L, the results agree with CFT. These plots give no further information compared to the ones already presented and we do not report them. We only show the results for a rather small system of length L = 39 (that nevertheless is above anything obtainable by exact diagonalization). In Fig 4 we report the resulting F lat 2 for all the possible divisions in four parts of this chain of length L = 39. For such small chain, the results are obviously very unclear since the corrections to the scaling are obscuring the CFT scaling represented by a continuous line that is surrounded by the points, signaling the oscillatory nature of the corrections (for clarity, compare with the analogous plot for the Ising model in next section).

Critical Ising universality class
The Ising model is given by Hamiltonian (9) with γ = 1 and it is critical for h = 1. Results for this model have been already derived numerically for α = 2 in Ref. [13] by using a tree tensor network algorithm and Monte Carlo simulations of the two-dimensional classical problem in the same universality class. These calculations allowed a precise determination of F 2 (x), but the system sizes explored were not  enough to analyze Rènyi entropies with larger values of α. In the course of our analysis we always compared our data for small systems with those in Ref. [13], in order to check the correctness of both methods.
We report our results in the thermodynamic limit for various values of = 1 = 2 at different separations r (resulting in the four-point ratio x given in Eq. (60)) in Fig.  5. As already noticed in Ref. [13], oppositely to the XX chain, in the Ising model we have monotonic finite correction to the scaling. For finite , the results do not show the symmetry x → 1 − x valid for infinite . This is restored only by the extrapolated data at → ∞. To perform this extrapolation in the most accurate way, we have first to determine the correction to the scaling exponent. For a single interval, it is exactly known that the leading corrections are characterized by the exponent δ α = 2/α as for the XX chain (the origin of the equality [24] is due to the equivalence of the two RDMs [31]). However for α = 2, it has been already shown [13] that δ 2 = 1/2, different from the single interval one. This result is quite surprising also in view the recent CFT analysis [29] predicting the same behavior for any number of intervals. Thus as a first step we check the corrections to the scaling exponent for the Ising model in the absence of the Jordan-Wigner string between the two blocks (i.e. we consider only the correlation matrix Γ 1 , as done in Ref. [8]). In this case, the results reported in Fig. 6 give a compelling evidence that the leading corrections to the scaling are given by δ α = 2/α as in the single interval (data not reported here for α = 3, 4 show that this remains true). At this point it is natural that the Jordan-Wigner string produces another operator at the conical singularity, that in the Ising model is the leading one. According to Ref. [29] all the corrections to the scaling should be of the form δ α = 2x/α, thus taking the result δ 2 = 1/2 for granted (see also below for further evidences), we conclude that the Jordan-Wigner string introduces an operator with scaling dimension x = 1/2. Such an operator in the continuum limit of the Ising model exists and it is the Majorana fermion, that has exactly the same features of the Jordan-Wigner string (i.e. same symmetry and same non-local character). Such an operator is clearly not present in the single interval case. These considerations allow us to conclude that the leading corrections to the scaling for the double interval entanglement in the Ising model are described by the exponent Unfortunately already for α = 3, the value of δ 3 = 1/3 is very low and subleading corrections to the scaling going with exponents mδ α (with m integer) are expected to influence the results in a considerable manner. For this reason, in order to have an accurate determination of the asymptotic behavior, at fixed x we consider all corrections to the scaling up to those with exponent 1. The resulting extrapolated data are the top points in Fig. 5. In the case of α = 2, the CFT prediction F 2 (x) (cf. Eq. (8)). perfectly with the extrapolated data agree, giving strong support both for the procedure to account for the subleading corrections terms and on the asymptotic form. For α = 2 no CFT prediction is still available. The extrapolated data in Fig. 5 are the first data for infinite . We now turn to consider the issue of universality. All the critical models (h = 1) for any value of γ = 0 are in the Ising universality class. However, the results at finite show a strong dependency on γ (as obvious from the different correlation matrices). In Fig. 7 we report several data for F lat 2 (x) for different values of γ and at fixed x = 1/4 and x = 3/4 (that for the symmetry x → 1 − x have the same asymptotic value). At finite , all results are evidently different. In the figure we report the extrapolation with two corrections to the scaling (i.e. with δ 2 = 1/2 and 2δ 2 ). For → ∞ all data tend to the same value predicted by CFT F 2 (1/4) = (2 + √ 2(1 + 3 1/4 ))/4, confirming in a single plot many results: (i) universality with respect to γ, (ii) correctness of the correction to the scaling form, (iii) correctness of the CFT prediction Eq. (8).
To conclude this section we report the data for a finite chain. As usual, in all the scaling variables we substitute distances with the chord distances. In Fig. 8, we report the values of F lat 2 (x) for all the possible choices of intervals A 1,2 and B 1,2 in a chain of L = 39 spins. Compared to the analogous plot for the XX chain (Fig. 4), the figure is much clearer, due to the fact that corrections to the scaling are monotonous. Notice however that the data points lie much below the asymptotic value because the exponent δ 2 = 1/2 is small (also compared to the XX case with δ 2 = 1) resulting in very large corrections to the scaling.

Non critical models
The richness of the phase diagram of the XY model allows us also to explore gapped phases that are almost everywhere except on the line |h| = 1 and the segment γ = 0 with h 2 < 1. These systems have not be considered so much in the literature until now, because for non-critical models all the correlations and entanglement between the two blocks fall off exponentially (with a decay rate given by the inverse gap or mass). Thus one would always expect S α ( , r, ) = S α ( ) + S a ( ) + O(e −r/∆ ) .
However, this is not so obvious because of the importance of the connected part in the correlations. In Fig. 9 (left) we report S 2 for the double interval case. While for h > 1, the spin entropy is the same as the double of the single interval and the same as the fermionic one, for h < 1 there is clearly an offset (that we quantify in − log 2, see below). The best way of detecting these unexpected effects is to consider the mutual entropy that gives automatically zero when factorization occurs. From numerical data we deduce that, for large blocks, all terms in Eq. (49) have the same absolute value. However, they can have different signs. Furthermore numerical evidence suggests that terms consisting only of the correlation matrices Γ 1 and Γ 2 have always positive sign. If these observations are generally true, we have where [{ζ}] is the sign associated to the element {Γ ζ1 , · · · , Γ ζα } and # 3(4) is the number of correlation matrices Γ 3 (4) . We expect an eventual discontinuity in ∆S α when crossing a critical line. Thus we study non-critical chains with magnetic field close to h = 1. We found numerically that only terms with odd number of correlation matrices Γ 3(4) display sign changes. Thus (see right panel of Fig. 9 for the explicit plot) we conclude from the numerical evidence, the behavior This result is not completely unexpected: also S A of the single interval for |h| = 0 tends to log 2, independently on [5,23,30]. Thus in the definition of ∆S α , since also the double interval S α tends to the same value, we are left with a single log 2. ‡ We stress that at h = 1 the various ∆ α cross in a single point if they are characterized by the same four point ratio x (see right panel of Fig. 9) that is the CFT prediction. These kind of plots could be used to detect the phase transition points in systems where it is not exactly known.
The α independence of the previous expression, allows us to analytical continue the result to α = 1 and to conjecture the same behavior for the asymptotic von Neumann entanglement entropy in non-critical regions where S ferm 1 is the fermionic entanglement entropy obtained from the correlation matrix Γ 1 . Notice that this does not exclude that exponential corrections to this asymptotic form in the fermionic and spin variables could have different amplitudes. It would be interesting to explore this issue with the continuum theory in the formfactor approach [32].

Entanglement evolution following a quench
Finally we consider the evolution of Rènyi entropies after global quenches. In this problem, the system is prepared in the ground-state |ψ 0 of an Hamiltonian H 0 of the form (9) and then it is let evolve from time t = 0 with a different Hamiltonian H, always of the form (9) but with (h, γ) = (h 0 , γ 0 ). The main feature of this nonequilibrium problem is that the initial state differs globally from the ground state and the excess of energy (compared to the ground-state of H) is extensive. A connected spin block reacts to the quench increasing Rènyi entropies linearly in time up to each spin in the subsystem becomes entangled with the environment. Then entropies saturate. There is no reason why disjoint blocks should behave in different way, except ‡ This is a consequence of the double degeneration of the ground state for |h| < 1. It is easily understood at h = 0, where the ground state is any linear combination of the states all up and all down that we can denote with | ↑ and | ↓ . Since there is no symmetry breaking term in the Hamiltonian, the diagonalization selects a state with zero magnetization, i.e. (| ↑ ± | ↓ )/ √ 2, in which the entanglement of any subsystem, connected or not, is always log 2 as stated in the main text.
for the necessity to subtract the mutual entropy. Our only goal here is to understand the differences between the spin representation and the fermionic one, postponing any more accurate analysis to further studies. We only report results for the Rènyi entropy S 2 , but the value of α is unimportant. We compare the numerical data with the prediction that follows from the interpretation of the entanglement evolution in terms of motion of quasiparticles [33]. According to this physical scenario, because the excess of energy of the system is extensive, any site acts as a source of quasiparticle excitations. Particles emitted from different points (further apart than the correlation length in the initial state) are incoherent, but pairs of particles moving to the left or right from a given point are entangled. Thus S α (t) should just be proportional to the number of coherent particles that emitted from any point reach one a point in A and the other a point in B. Since there is a maximum speed for these excitations v M , for a single interval this implies the linear growth for 2v M t < and saturation for very large times. For a general bipartition the result is [33] where we assumed momentum conservation. For a double interval, this formula predicts a series of linear behavior with different slopes, that finally saturates at late time [33]. The function H α does not depend on the subsystem length and topology. It has been exactly derived for a single interval [34] H α (cos ∆) = 1 where ∆(p) is the difference between the Bogolioubov angles before and after the quench, whose explicit expression in term of the dispersion relation ε(ϕ) is reported in [34]. Thus, plugging everything together, the leading order of Rènyi entropies is , 2),r+ 1+ 2] (v) and [34] Let us now go back to the direct computation of the entanglement in two blocks after a quench, we observe that not all terms contribute to the leading order in blocks lengths and distance: the string expectation value decays exponentially and only terms constructed with Γ 1 and Γ 2 survive. For instance the Renyi entropy S 2 is simply the logarithm of two terms (75) Figure 10. Evolution of the Rènyi entropy S 2 after a global quench from a non critical system to a critical Ising. The asymptotic prediction is given by equation (72) while the fermionic entanglement is obtained neglecting the string contribution. When {Γ 1 , Γ 2 } becomes negligible with respect to {Γ 2 1 }, a jump of about ∼ log 2 is seen. In the graph this happen in the neighborhood of the time t ∼ 11. The asymptotic prediction presents a small offset due to the finite entanglement in the initial state.
Analyzing several data (we report only a single example in Fig. 10), we deduce that in the scaling limit the term constructed with the standard fermionic correlations is never negligible: finite order Renyi entropies are controlled by fermionic correlations. When other terms are comparable with {Γ 1 , · · · , Γ 1 }, they give additive O(1) contribution of the form ∼ log k, where k is the coefficient in front of the factor {· · ·}, and k = 2 for S 2 .
The quasi-particle interpretation [33] perfectly agrees (up to a constant) with the fermionic representation ρ Γ1 (see Fig. 10). In other words the non locality of the Jordan-Wigner transformation does not influence (as one expects) the leading order of time evolution of entropies after a global quench. However, as in the non-critical equilibrium case, we observe an additive log 2 contribution (in Fig. 10 a stick of width log 2 is included to appreciate this difference) from the string and we believe that this is a general feature that will persist also away from criticality. As a final comment, we mention that the O(1) offset of the numerical data at finite compared with the asymptotic form is a consequence of the initial entanglement. This subleading term is present also for the single interval case [34], and can be included in the asymptotic result [35].

Discussions
We presented a general method to calculate Rènyi entropies for integer α in the ground state of the XY Hamiltonian (9) in the case when A consists of two (but the method works for even more) disjoint intervals. We carefully analyze the results in the critical cases of XX and Ising universality classes. We found that the asymptotic results for large intervals and separations are described correctly by conformal field theory when results are known. We also provided results for the asymptotic scaling of the α = 3, 4 Rènyi entropies of the Ising model, for which no analytic prediction is still available. In order to check CFT predictions we had to carefully include corrections to the scaling in the analysis. We found that for the XX model, the leading corrections to the scaling are governed by the exponent δ α = 2/α, the same as the single interval case [24], in agreement with general predictions [29]. Oppositely for the critical Ising model we found δ α = 1/α, that is half of the single interval one. This can be understood in terms of the general theory of the corrections to the scaling [29] by interpreting the Jordan-Wigner string in the lattice model as a non-local fermion in the Ising CFT. We also have briefly considered the entanglement in the gapped phase and the time evolution following a quantum quench. The provided exact representation of the reduced density matrix of spin variables could also be useful for the calculation of the entanglement between A 1 and A 2 measured by the negativity [15].
A major problem remains still open. We are unable to find all the eigenvalues of the reduced density matrix of two intervals and so the Rènyi entropies for noninteger α and in particular for α = 1 that would give the widely studied von Neumann entropy. The same problem is also present in conformal calculations. It is possible to obtain Tr ρ α A only for α integer [12] and the analytic continuation of the result to general complex values remains a big open problem.
The contribution of the Jordan-Wigner string also affects the entanglement entropy of a single interval in a system with one or more boundaries, if the block A does not include the boundary. In this case, the correlation matrix is slightly more complicated and we are currently analyzing the problem. Among the other things, this would be relevant also for some quench problems where indeed only the fermionic entanglement has been calculated [36].
Another interesting open question concerns the behavior of the Rényi entropies in systems with disorder for which in the single block case it is known that the entanglement simply follows from the counting of the singlets [37] between A and B. Whether this remains true in the case of more intervals is not yet known.

Appendix A. Correlation matrices
We report in this appendix some simple properties of the correlation matrices with ρ S defined in Eq. (29). Γ S 's elements are ratio of two expectation values (n = 0) which can be evaluated by means of the Wick theorem. We show that Γ S can be expressed in terms of the system correlation matrix Γ ≡ Γ 1 , where 1 is the identity matrix. By isolating the two-point correlations between a l and the other Majorana operators in the numerator of (A.2) we obtain a l+n j∈B1 (a 2j−1 a 2j )a l = S Γ l l+n + ΓΛΓ l l+n , (A.4) where B * ij denotes the set of all indices {2l} ∪ {2l − 1} with l ∈ B 1 and with i and j removed. The expectation value in (A.4) is the Pfaffian of the system correlation matrix restricted to the region B * ij . Checking the Pfaffian sign we have (A.5) with Γ B1 the correlation matrix restricted to the region B 1 . Substituting this expression into (A.2) we get where the double subscript takes into account restrictions to rectangular correlation matrices, i.e. the first (second) subscript identifies the region where the row (column) index runs. Notice that Γ S is the Schur complement of Γ B1 in Γ A∪B1 . This proves Eq. (55).

Appendix B. Generalization to many disjoint spin blocks
We consider in the appendix how to generalize our approach to the case when the subsystem consists of n disjoint blocks of lengths i , with i = 1, . . . , n. We indicate with O + i a product of an even number of Majorana operators lying in the i-th block A i and with O − i an odd product. S B i is the complete σ z string associated to the i-th region between consecutive blocks B i and S A i is the string corresponding to the block A i (see Fig. 1 for the double block case). Eq. (20) can be easily generalized observing that the reduced density matrix has the following expansion where s i = ±1 and S[{s}] are σ z strings determined by the conditions Indeed, any spin operator that can be written as a product of Majorana operators that are odd in number if restricted to the k-th block is constructed with an odd number of σ x(y) of the k-th block. Spin variables σ x and σ y are non-local in the fermionic space where Λ(A 1 , . . . , A k−1 ) is an operator that depends on the Majorana fermions of the blocks A 1 , . . . A k−1 . When multiplying an odd number of σ x(y) the string k−1 i=1 S i does not cancel (we remind S 2 i = 1) so we obtain the conditions (B.2). The construction of the fermionic equivalent reduced density matrix is analogous to the double block case. In fact the unitary transformation brings the spin reduced density matrix to a fermionic one. In particular we obtain the simple result The term associated to the configuration {s} is sum of 2 n−1 fermionic RDMs. This is because, fixed S[{s}], the further dependence on s i can be handle in the following way (remind where (−1) j∈{i} is equal to −1 if j ∈ {i} and 1 otherwise, and Rènyi entropies can be computed in the same way as when the subsystem consists of just two disjoint blocks. But many more terms contribute. And when the number of blocks is comparable with the chain size we expect an extensive behavior of Renyi entropies, as observed in [18], that is reminiscence of the huge number of fermionic RDMs needed to represent the spin RDM. Finally, as in the double block case, Rènyi entropies do not depend on the sign of Pfaffians because any sign change is the result of a unitary transformation (that is just a σ z string).

Appendix C. Some technical details for the product rules
The reduced density matrix of a given number of disjoint spin blocks in a free chain is unitarily equivalent to a finite sum of operators that well behave under multiplication.
In fact one of our results is that, for some σ z strings S j and constants C j , the relation holds. Like in the single block case each operator ρ Sj turns out to be a normalized exponential of a quadratic form in the Majorana fermions. Up to a constant factor, these matrices are closed under multiplication and we identify any of them by the skew-symmetric matrix W in the exponent in Eq. (33). We recall in this appendix some basic properties of skew-symmetric diagonalizable matrices. If λ is an eigenvalue of W , then also −λ is because If w λ and w µ are two eigenvectors with eigenvalues λ and µ, then The vectors w µ and w λ are orthogonal (with respect to the scalar product w µ · w λ = j (w µ ) j (w µ ) j , that is not the standard one because of the absence of complex conjugation) unless µ = −λ. If the eigenvalue λ (and so −λ) is n-times degenerate, in general one can redefine the eigenvectors in the same n-dimensional eigenspace so that each vector w i λ becomes orthogonal to the other n − 1 ones. Thus the following decomposition holds where {λ}/ ± means that each couple ±λ contributes only once in the sum and the symbol v ⊗ w identifies the matrix with elements v i w j . By introducing the auxiliary matrices which satisfy together with Π − λ the orthogonality relations with f e and f o the even and odd part of f . The decomposition (C.4) allows to express ρ W in a factorized form where we made use of w λ 's orthogonality (C.3) We obtain Z(W ) by imposing the normalization of ρ W Z(W ) = Tr exp aW a 4 = . (C.10) The correlation matrix Γ ij = δ ij − a i a j comes directly from the quadratic part of expression (C.8) The correlation matrix does not determine univocally function Z(W ). Indeed we obtain the expression 13) with an evident phase ambiguity (that is actually a sign ambiguity). In spite of the ambiguity in Z(W ) ∼ Z[Γ], we believe that the trace of two fermionic RDMs is a functional of Γ = tanh(W/2) and Γ = tanh(W /2). In fact, we argue that where the product is extended over all eigenvalues λ of 1+ΓΓ with degeneracy (always even) reduced by half. We have not a rigorous proof of Eq. (C.15), but we give an argument to justify it. After infinitesimal variations of Γ and Γ , function Z(W ) changes in the following way Using Eq. (44)  The variations dΓ and dΓ are skew-symmetric by construction. Our task is to integrate the differential equation above. We consider all possible sets of 2N smooth simple curves joining the ΓΓ eigenvalues λ i from ΓΓ = 0 (when Tr ρ W ρ W = 2 −N ) to our desired values. Despite of the initial ambiguity in the sign of Z(W ), the expression for the logarithm is well defined, in view of the fact that ΓΓ has a double-degenerate spectrum and the ambiguity related to the choice of the curves γ i (avoiding any singular point) gives simply an additive constant proportional to 2πi. Considering the path characterized by d(ΓΓ ) = dλ(ΓΓ − 1)/2, the product rule can be expressed in the following integral form Tr ΓΓ − 1 2 − λ + λΓΓ dλ = exp 1 2 where γ is a smooth simple path in C which joins 0 to 1 and avoids any of the points at which 1 + λΓΓ fails to be invertible. Because of the factor 1/2 in front of the integral we get equation (C.15). For this proof we used several properties about Pfaffian reviewed in Ref. [38]. can be 1 or −1 and numerical data suggest that it is independent of the subsystem and of the chain length. In particular we find s = 1 |h| ≥ 1 or γ = 0, −1 otherwise, (C.22) so that We close this appendix with some speculations and some simplified (wrong) expressions to which one could arrive if treating too naively correlation matrices. If we would not worry about the signs and non-invertibility, we could write directly the result of a generic product like Eq. (48) as This expression is generally wrong. When substituting the chain correlation matrices, (1 − Γ 1,3 ) is usually non-invertible (Fig. C1). Thus, except for S 2 , equation (C.24) is wrong. By using this wrong expression one can try to see if any simplification occurs in the general Rènyi entropy. In general it is not possible to push this approach, unless taking further wrong assumptions. Assuming for example the (wrong) hypothesis of commuting Γ matrices, after trivial algebra we have (the wrong expression) We wrote this expression here, because one could have wonder whether it (despite wrong) would have been defined for any real α (i.e. depending only on the eigenvalues of the matrices Γ) to have an expression analytically continuable to α → 1 (and maybe correct in this limit). However, even this simplified wrong expression is only calculable for integer α and thus we did not insist in this direction and we preferred to exploit the exact solution without assumptions.

Appendix D. Accidental Singularities
When the correlation matrix of the region interposed between the blocks Γ B1B1 is non-invertible, that is the σ z string mean value vanishes S = 0, some problems arise. The density matrix (32) is not well defined and all successive discussions are meaningless. Here we give some tricks to solve these problems. We consider a generic operator can be used to perform any product at the price to have the string S E . We observe that (D.5) where P E is the diagonal operator with elements i = −1 if i ∈ E and 1 otherwise. In fact, exploiting the commutation relation, products of operators like Q W can be carried to the form is free from pathologies.