Entanglement entropy of inhomogeneous XX spin chains with algebraic interactions

We introduce a family of inhomogeneous XX spin chains whose squared couplings are a polynomial of degree at most four in the site index. We show how to obtain an asymptotic approximation for the Rényi entanglement entropy of all such chains in a constant magnetic field at half filling by exploiting their connection with the conformal field theory of a massless Dirac fermion in a suitably curved static background. We study the above approximation for three particular chains in the family, two of them related to well-known quasi-exactly solvable quantum models on the line and the third one to classical Krawtchouk polynomials, finding an excellent agreement with the exact value obtained numerically when the Rényi parameter α is less than one. When α ≥ 1 we find parity oscillations, as expected from the homogeneous case, and show that they are very accurately reproduced by a modification of the Fagotti-Calabrese formula. We have also analyzed the asymptotic behavior of the Rényi entanglement entropy in the non-standard situation of arbitrary filling and/or inhomogeneous magnetic field. Our numerical results show that in this case a block of spins at each end of the chain becomes disentangled from the rest. Moreover, the asymptotic approximation for the case of half filling and constant magnetic field, when suitably rescaled to the region of non-vanishing entropy, provides a rough approximation to the entanglement entropy also in this general case.


Introduction
The entanglement entropy of spin chains of XX type -or, equivalently, systems of free spinless fermions with nearest-neighbors hoppings-has been intensively studied since the seminal work of Jin and Korepin [1] for the homogeneous chain.Indeed, the bipartite entanglement entropy of one-dimensional models is a convenient indicator of their criticality.The reason is that in their critical phase these models are effectively described at low energies by a (1 + 1)-dimensional conformal field theory (CFT), whose entanglement entropy has been shown to scale logarithmically with the block size L [2,3].In fact, a fundamental property of all XX-type spin chains is the fact that their entanglement entropy can be expressed in terms of the eigenvalues of a (truncated) correlation matrix.In the homogeneous case this matrix is Toeplitz (for closed chains) or Toeplitz+Hankel (for open ones), which makes it possible to apply proved instances of the (generalized) Fisher-Hartwig conjecture [4][5][6] to rigorously derive the leading asymptotic behavior of the entanglement entropy.In this way it was shown that the Rényi entanglement entropy S α of the homogeneous XX spin chain is asymptotically proportional to log L in the open, closed and (semi)infinite cases [1,7], even for subsystems of more than one block [8][9][10][11][12].
Moreover, the coefficient of log L in the asymptotic formula for S α confirms that this model has central charge c = 1, as expected.
In fact, corrections to the logarithmic behavior of S α in the limit of large L were exhaustively analyzed by Calabrese and Essler for the closed homogeneous XX chain [7], and by Fagotti and Calabrese for the open one [13].For α > 1, or α 1 in the open case, these terms present an oscillatory behavior which is particularly simple in the open case.Indeed, in this case the leading order correction is proportional to sin((2L + 1)k F )L −1/α , where k F is the Fermi momentum.Moreover, it is argued in Ref. [13] (and earlier in [14]) that this correction -more precisely, the exponent of L in the formula for the amplitudeencodes additional information about the underlying CFT beyond the central charge.
The situation is less straightforward in the non-homogeneous case, since the correlation matrix is in general neither Toeplitz nor Toeplitz+Hankel.However, at half filling and constant magnetic field, the leading behavior of the bipartite entanglement entropy can be derived through the technique first used in Ref. [15] to study the rainbow chain [16].The main idea is that with these assumptions the chain's continuum limit yields the CFT of a massless Dirac fermion in a static curved (1 + 1)-dimensional spacetime, whose metric's conformal factor is proportional to the square of (the continuum limit of) the hopping amplitude.This suggests that in the thermodynamic limit the leading asymptotic behavior of S α can be obtained from the formula for the homogeneous case replacing the chain's and block's lengths by their conformal versions [17].This was actually shown to be the case for the rainbow chain in Ref. [18], and more recently for several other inhomogeneous XX chains in Refs.[19,20].We stress, however, that to the best of our knowledge the latter results only hold in the case of half filling and constant magnetic field.Moreover, the method just outlined has only been applied to the leading term in the asymptotic formula for S α , without addressing the behavior of the subleading corrections.
Another fundamental property of XX spin chains is their close connection with classical orthogonal polynomials.Indeed, the chain's single-particle Hamiltonian is represented in the position basis by a real tridiagonal symmetric matrix (the so-called hopping matrix), whose elements can in turn be used to define a three-term recursion relation determining a finite orthogonal polynomial system (OPS) {P n } N n=0 , where N is the number of spins.This establishes a one-to-one correspondence between XX spin chains and OPSs, that can be used to derive in a simple way many of the chain's properties.Indeed, the single-particle energies are the roots of the critical polynomial P N , and the correlation matrix elements can be computed in closed form (without need of numerical diagonalization) in terms of the polynomials in the OPS evaluated at the latter energies.This turns out to be more efficient than brute force diagonalization of the matrix of the single-particle Hamiltonian as the number of spins grows.This connection has also been exploited in Ref. [21] to construct in some cases a tridiagonal matrix commuting with the hopping matrix of the entanglement Hamiltonian, which can be used to improve the numerical accuracy of the eigenvalues of the latter matrix and hence of the entanglement entropy.
A key property shared by the chains studied in Ref. [21] is the fact that the square of the interaction strength J n is a polynomial of degree at most four in the site index n.This property is very natural from the point of view of the associated orthogonal polynomial family, since −J 2  n−1 coincides with the coefficient of P n−1 in the recursion relation for P n+1 .In fact, the latter property also holds for the inhomogeneous XX chains related to onedimensional quasi-exactly solvable (QES) models [22][23][24][25] introduced in Ref. [19].We shall show in this work that for all XX spin chains for which J 2  n is a polynomial of degree up to four in n it is possible to evaluate in closed form the leading asymptotic approximation to the ground-state entanglement entropy (at half filling and in a constant magnetic field) by the procedure introduced in Ref. [15].This in done essentially by reducing a suitable elliptic integral to Legendre canonical form, a procedure which depends on the root pattern of J 2  n .
We shall analyze in some detail three inhomogeneous XX chains with interactions J n of the algebraic form described in the previous paragraph (algebraic interactions, for short).Two of these chains arise from well-known QES potentials, namely the sextic oscillator and the Lamé periodic potential, while the third one, introduced in Ref. [21], is associated to the Krawtchouk discrete orthogonal polynomial family.We first of all check that the leading term in the asymptotic approximation to the entanglement entropy obtained through the related massless Dirac fermion CFT in a suitably curved background is in excellent agreement with the numerical results in the standard scenario of constant magnetic field and half filling.Moreover, for the sextic and Krawtchouk chains at half filling and in a constant magnetic field we have found strong numerical evidence that the subleading (constant) term in the asymptotic expansion of the entanglement entropy coincides with its counterpart for the homogeneous XX chain for N large enough.To the best of our knowledge, this remarkable coincidence had not been previously noticed in the literature.We stress in this regard that the connection of the models under study with families of orthogonal polynomials makes it possible to determine the eigenvalues and eigenvectors of the single-particle Hamiltonian in a numerically efficient way when the number of spins is very large.Our numerical calculations also indicate that when α 1 the Rényi entanglement entropy S α features parity oscillations which become more marked as α increases, as in the homogeneous case.Remarkably, these oscillations are reproduced with great precision by the heuristic formula proposed by Fagotti and Calabrese for the homogeneous XX chain, replacing the lengths of the block and the whole chain by their values computed with the metric of the ambient space of the associated CFT.More precisely, this formula depends only on two free parameters, whose fitted values are very close to the theoretical ones for the homogeneous XX chain.This underscores the essential similarity of the homogeneous and inhomogeneous cases, and the universality of Fagotti and Calabrese's formula for this class of models.
All of the above results have been obtained under the standard assumptions of constant magnetic field and half filling, for which the connection with the massless Dirac fermion CFT has a theoretical justification via the continuum limit.In this work we have also analyzed the non-standard situations of arbitrary filling and/or inhomogeneous magnetic field, which to the best of our knowledge have not been addressed in the literature.Our numerical calculations clearly indicate that the new feature in both of these scenarios is the vanishing of the entanglement entropy when the length of the block is small or close to the chain's length.In other words, the first few and last spins become disentangled from the rest of the chain.As a consequence, the (leading) asymptotic approximation to the entanglement entropy derived from the associated CFT cannot be expected to hold in this case.Remarkably, we have checked that this approximation roughly reproduces the average behavior of S α if suitably scaled to the region of non-vanishing entropy.On the other hand, the oscillations of S α in these non-standard cases are found to be much more complex than in the usual situation of half filling and constant magnetic field, and in particular are not well reproduced by the conformally modified version of the Fagotti-Calabrese formula.
The paper is organized as follows.In Section 2 we briefly review the connection between inhomogeneous XX spin chains and free fermion systems, and recall how the latter models can be exactly diagonalized.Likewise, in Section 3 we explain how these models are related to a finite orthogonal polynomial family through its three-term recursion relation, and how to exploit this connection to diagonalize the single-particle Hamiltonian.In Section 4 we summarize the main results on the bipartite entanglement entropy of spin chains used throughout the paper.In particular, we review some known results about the asymptotic behavior of the entanglement entropy of the homogeneous XX chain as the number of spins tends to infinity, and briefly outline their recent extension to the non-homogeneous case at half filling in a constant magnetic field.In Section 5 we present a class of inhomogeneous XX spin chains with algebraic interactions for which it is possible to compute in closed form an asymptotic approximation to the block entanglement entropy by the procedure explained above.The next three sections are devoted to the detailed analysis of three models in the previous family, associated to the QES sextic oscillator Hamiltonian (Section 6), the classical Krawtchouk polynomials (Section 7) and the periodic Lamé potential (Section 8).In Section 9 we present our conclusions and discuss several lines for future research.The paper ends with a technical appendix explaining how to reduce to Legendre canonical form the elliptic integral appearing in the asymptotic formula for the entanglement entropy.

Inhomogeneous XX spin chains and free fermion systems
The Hamiltonian of an inhomogeneous XX spin chain with interactions J n in an external magnetic field B n can be taken as where N is the number of sites and σ α n (with α = x, y, z) denotes the Pauli matrix σ α acting on the n-th site.In what follows we shall assume that the interaction strengths J n do not vanish.As remarked in Ref. [19], the model with J n replaced by ε n J n , where ε n ∈ {±1} is a site-dependent sign, is unitarily equivalent to the original one.Hence we can take all the J n 's to be positive without loss of generality.
It is well known that the Jordan-Wigner transformation where σ ± n := (σ x n ± iσ y n )/2, maps the Hamiltonian (2.1) into that of a system of N hopping spinless fermions, Here c † n (resp.c n ) is the operator creating (resp.destroying) a fermion at site n, and the coefficients J n and B n respectively represent the hopping amplitude and the chemical potential of the fermions.In what follows we shall mainly deal with the free fermion system (2.3), our results being easily translated to its spin chain equivalent (2.1).
In the homogeneous case (i.e., when J n and B n are site independent), the Hamiltonian (2.3) commutes with the translation operator along the chain sites and is thus diagonal in momentum space.In the non-homogeneous case this symmetry is lost, but the Hamiltonian can still be diagonalized by introducing suitable modes.More precisely, let denote the matrix of the restriction of H to the single-particle sector with respect to the position basis where |vac is the fermionic vacuum.Since the hopping matrix H is real and symmetric, it can be diagonalized by a real orthogonal matrix Φ = (Φ nk ) N −1 n,k=0 , namely where which satisfy the canonical anticommutation relations (CAR) on account of the unitary (real orthogonal) character of Φ.It is easily shown that the Hamiltonian (2.3) can be written as and is thus diagonal in the basis consisting of the states whose corresponding energy is given by In particular, the one-particle eigenstates c † k |vac (with 0 k N − 1) represent singlefermion excitation modes with energy ε k .
The case in which the magnetic field B n vanishes for all n deserves special attention.Indeed, in this case H is equivalent to −H under the unitary transformation c n → (−1) n c n (which obviously preserves the CAR), so that the spectrum is symmetric about zero: This implies that the system possesses particle-hole symmetry, since if Moreover, from the equivalence of H to −H under c n → (−1) n c n we immediately obtain the relation up to an n-independent sign.Thus when N is even the ground state is the half-filled state with Fermi momentum π/2 and energy (2.12) (When N is odd the ground state is doubly degenerate, since the zero energy mode does not change the total energy.)

Orthogonal polynomials
As explained in the previous section, the diagonalization of the full Hamiltonian (2.3) of a free fermion system is achieved by diagonalizing the hopping matrix H in Eq. (2.4).Since the latter matrix is tridiagonal, it can be used to define a finite orthogonal polynomial system {φ n (E)} N n=0 through the three-term recursion relation (with φ −1 := 0).It is easily shown that the polynomial φ N (E) is proportional to the characteristic polynomial of H, and that the matrix elements of Φ can be taken as Φ nk = φ n (ε k ), where for each k the constant φ 0 (ε k ) is determined up to a sign by the normalization condition Recall that the eigenvalues ε k of H are non-degenerate (i.e., the roots of P N are simple), and hence the orthogonality relations are automatically satisfied.As is customary, we shall work in what follows with the monic polynomial family {P n (E)} N n=0 , where P n is the unique monic polynomial proportional to φ n .From Eq. (3.1) it follows that and that the polynomials P n satisfy the normalized recursion relation with P −1 := 0 and a n = J 2 n−1 > 0 .Conversely, a monic polynomial OPS defined by a recursion relation of the form (3.2) with a n > 0 determines a free fermion system (2.3) with hopping J n = √ a n+1 and chemical potential B n .From the previous argument it follows that the one-particle energies ε k are the roots of the critical polynomial P N .It is also shown in Ref. [19] that the entries Φ nk of the real orthogonal matrix Φ determining the mode creation/annihilation operators through Eq. (2.6) can be taken as where Note that the orthogonality of the (real) matrix Φ follows directly from the fact that the family {P n } N −1 n=0 is orthogonal with respect to the discrete measure N −1 k=0 w k δ(E − ε k ), with square norm P n 2 = γ n (see Ref. [19] for the details).In the particular case in which B n vanishes for all n the recursion relation (3.2) implies that P n has the parity of n, i.e., P n (−E) = (−1) n P n (E).It then follows from the definition of w k that w N −k−1 = w k , and hence in agreement with Eq. (2.10).

Entanglement entropy
A quantitative measure of the entanglement entropy of a block A of spins of the chain (2.1) -or fermions in the system (2.3)-when the whole system is in a pure state |ψ is the entropy of the block's density matrix ρ A := tr A |ψ ψ|, where the subindex in the trace operator indicates that we are tracing over the degrees of freedom of the complementary set A := {0, . . ., N − 1} \ A. More precisely, we shall take A = {0, . . ., L − 1} and work with the Rényi entropy S α (where α > 0 is a real parameter) defined by whose α → 1 limit S 1 is the usual von Neumann-Shannon entropy In general, the system's state |ψ shall be taken as an eigenstate (2.8) with the first M energy modes excited: An important property of the free fermion system (2.3) is the fact that its eigenstates are Gaussian.Thus the analogue of Wick's theorem can be applied to express the Rényi entanglement entropy in terms of the eigenvalues ν n (n = 0, . . ., L − 1) of the L × L correlation matrix C(L, M ) ≡ C with entries through the formula for the von Neumann entropy (see, e.g., [1,26,27]).Using Eq. (2.6) it is straightforward to show that the correlation matrix can be computed from the matrix Φ -or equivalently, by Eq. (3.3), the OPS {P n } N n=0 -as In other words, C = Φ LM Φ T LM , where Φ LM is the matrix obtained by taking the first L rows and M columns of Φ.
Remark 1.The state |ψ M in Eq. (4.1) can always be regarded as the ground state by adding to the Hamiltonian (2.3) a homogeneous term This of course amounts to adding a multiple of the identity to the hopping matrix H, but does not change its eigenvectors (Φ 0k , . . ., Φ N −1,k ) (with k = 0, . . ., N − 1), and thus leaves the matrix Φ and the correlation matrix C invariant.Since, as explained above, the entanglement entropy is determined by the eigenvalues of C, the entanglement entropy is also invariant.
When 0 < M/N < 1 the system (2.3) is gapless, and is thus described at low energy by an effective (1 + 1)-dimensional CFT.The asymptotic behavior of the entanglement entropy of such a theory in Minkowski spacetime was determined in Refs.[2,28].More precisely, when the spatial manifold is a finite interval [0, L] the entanglement entropy of a subinterval [0, ] is given by where c is the central charge of the CFT, c α is a non-universal constant depending only on the Rényi parameter α, and ∆x is an ultraviolet cutoff.In particular, the leading asymptotic behavior of S α is entirely determined by the central charge c, and is thus universal.
For a homogeneous free fermion system (2.3) the correlation matrix C is "Toeplitz plus Hankel", which makes it possible to derive the leading asymptotic behavior of S α in the limit L, N → ∞ with [13] At half filling this result is in agreement with the CFT formula (4.2) with central charge c = 1 taking ∆x as the chain's spacing, so that This confirms the fact that in the homogeneous case the free fermion system (2.3) (which is equivalent to the homogeneous Heisenberg XX chain) is described by the free fermion CFT (in Minkowski spacetime) with c = 1.
In fact, the asymptotic behavior of the entanglement entropy of the homogeneous XX chain is known in much greater detail [13].To begin with, in this case the non-universal constant c α is given by for α = 1 (and its α → 1 limit for α = 1 [1]).Actually, Eq. (4.2) holds for an arbitrary filling (with c α as above) if we add the extra factor sin k F to the argument of the logarithm, where k F := πM/N is the Fermi momentum.For α < 1, the o(1) term is actually of order N −1 , and thus the formula with provides an excellent approximation to the Rényi entanglement entropy of the homogeneous XX chain even for moderately large values of N .On the other hand, for α 1 the o(1) term in Eq. (4.2) features parity oscillations which for large α can even obscure the leading asymptotic behavior (4.5).A remarkable heuristic formula for these parity oscillations was found in Ref. [13] using CFT arguments, namely ) where c α is given by Eq. (4.4) and with µ 1 = lim α→1 µ α = −1/4.This formula reproduces with great precision the parity oscillations of the Rényi entropy of the homogeneous XX chain when α 1.It was also argued in the latter reference that a subleading term proportional to f H (N, λ) −κ/α -though not the coefficient µ α or even the oscillatory term sin ((2L + 1)k F )-is in fact universal, the parameter κ (which is unity for the homogeneous XX chain) providing information on the scaling dimensions of relevant operators in the associated CFT.
In the general (non-homogeneous) case the chain (2.3) is no longer described by a CFT in Minkowski spacetime, and thus the previous considerations do not directly apply.However, when N is even and B n = 0 for all n -i.e., when the system's ground state is the half-filled state (2.11)-it was shown in Refs.[15,18] that the continuum limit of the Hamiltonian (2.3) coincides with the Hamiltonian of a free massless Dirac fermion in the curved spacetime with static metric (4.9) Here J(x) is the continuum limit of J n , obtained by setting n∆x =: x n , taking the limit N → ∞ and ∆x → 0 with L = (N − 1)∆x fixed, and replacing x n by a continuous variable x ∈ [0, L].The metric (4.9) can be expressed in isothermal coordinates as with d x = dx/J(x).It is therefore natural to assume that in the limit N → ∞ with L/N → λ finite the entanglement entropy of the free fermion system (2.3) with even N and B n = 0 for all n -or more generally, by Remark 1, with B n constant at half fillingcan be obtained from Eq. (4.2) with , L and ∆x respectively replaced by the conformal lengths where is the length of the spatial interval [0, x] computed with the metric (4.10).In other words, we should have for a suitable (non-universal) constant c α (not necessarily given by Eq. (4.4)).This was shown to be the case for the rainbow chain (for which J n = J 0 e −h|n/N −1/2| , J(x) = J 0 e −h|x/L−1/2| ) in Refs.[15,18], and more recently for the Lamé [19], Rindler and sine chains [20].
An interesting open problem motivated by the previous considerations is whether the more precise asymptotic approximations (4.5)-(4.7)also hold in the non-homogeneous case after the replacement (∆x, , L) → (∆ x, , L) in Eq. (4.6).In fact, since the equivalence of the continuum limit of the inhomogeneous chain (2.1) with a CFT in curved spacetime has only been established at half filling (and for zero magnetic field), the latter formulas are only expected to apply when k F = π/2 and B n vanishes (or, more generally, is constant).We are thus led to conjecture the following more detailed asymptotic formulas for the Rényi entanglement entropy of the general (non-homogeneous) chain (2.1) at half filling in a constant magnetic field: with In what follows we shall introduce a family of inhomogeneous XX chains for which it shall be checked that the above conjecture holds.We shall also show that the analogous generalization of Eqs.(4.5)-(4.7) to the case of arbitrary filling and/or inhomogeneous magnetic field is not valid for the models considered in this paper.

Spin chains with algebraic interactions
In order to evaluate the right-hand side of Eqs.(4.14)-(4.15) in closed form it is necessary to compute the integral in Eq. (4.12).We shall introduce in this section a large class of inhomogeneous XX chains for which the latter integral can be explicitly evaluated.This class is characterized by the fact that the coefficient a n in the recursion relation (3.2) is a polynomial of degree at most four in n, and thus J n and J(x) are algebraic functions of degree two.As we shall discuss in the sequel, for these chains the RHS of Eq. (4.12) is an elliptic integral which can be evaluated by transforming it to Legendre normal form.
In fact, chains with this type of algebraic interactions have been recently discussed in the literature in two different contexts.Indeed, the inhomogeneous chains associated to the discrete Krawtchouk and dual Hahn polynomials studied in Ref. [21], whose entanglement Hamiltonian admits a commuting tridiagonal operator, both feature interactions of the above form.The same is true for all spin chains related to quasi-exactly solvable quantum models on the line recently constructed and classified in Ref. [19].In particular, it was shown in the latter reference that the entanglement entropy of an inhomogeneous spin chain associated to the quantum Lamé potential is indeed well approximated by the asymptotic formula (4.13) in the limit of large N .Consider, then, the integral (4.12).Since J(x) is dimensionless (in natural units), it must be a function of the dimensionless variable ξ := x/L.Setting we can rewrite (4.12) as We shall assume in what follows that p is a polynomial with real coefficients, with deg p 4 and p(ξ) 0 for 0 ξ 1.The main idea for reducing the last integral to canonical form is the fact that a real projective change of variable transforms it into an integral of the same type.Indeed, p az+b cz+d a polynomial of degree at most four in z.Using a projective change of variable of the form (5.3), the original polynomial p(ξ) can always be transformed into a suitable canonical form p(z), which is completely determined by the root pattern of p(ξ).
To begin with, it is clear that the integral (5.2) can be transformed into an elementary integral (expressible in terms of rational, trigonometric or hyperbolic functions and their inverses) if p(ξ) has a multiple root.Indeed, if p(ξ) has a multiple root at infinity (i.e., if ξ 4 p(1/ξ) has a multiple root at the origin) then deg p 2, and the integral (5.2) is elementary.Otherwise, if ξ = ξ 0 is a multiple (finite) real root of p the projective transformation z = (ξ − ξ 0 ) −1 transforms p(ξ) into a polynomial p(z) with a multiple root at infinity, i.e., a polynomial of degree at most two.Finally, if p(ξ) has a pair of complex conjugate double roots with c > 0, and the integral (5.2) is again elementary.
In view of the above discussion, we need only consider the case in which all the roots of p(ξ) (real or complex) are simple.In this case (5.2) is a genuine elliptic integral, which can be reduced to its standard Legendre form by the general procedure described, e.g., in Ref. [29].We present in the appendix a simplified version of this procedure adapted to the integral (5.2).The conclusion of this analysis is that in all cases the integral (5.2) can be expressed in terms of the incomplete elliptic integral of the first kind In the following sections we shall present several examples of algebraic inhomogeneous XX spin chains, including the Krawtchouk and Lamé chains previously mentioned, for which the integral (5.2) can be computed in closed form by the procedure described above, and thus the RHS of the asymptotic formula (4.14)-(4.15)can be readily evaluated.We shall study the applicability of the latter formula both in the standard situation considered in the literature of constant B n and half filling, and also outside this regime.We shall verify that Eq. (4.14) is an excellent approximation for the Rényi entanglement entropy in the standard situation, but this is not the case for inhomogeneous magnetic fields and/or other fillings.

The sextic chain
As our first example, we shall consider the inhomogeneous XX chain associated with the QES sextic oscillator potential [22,23,30], whose parameters are given (up to irrelevant constants) by [19] with γ = 0 or γ > 1/2.We shall start by considering the case β = 0, for which the magnetic field term vanishes identically and the asymptotic approximation (4.13) to the entanglement entropy should hold.Note that for finite γ the hopping amplitude J n is not symmetric about the chain's midpoint, i.e., J n = J N −2−n .On the other hand, for γ → ∞, or more precisely when γ N , we have To begin with, we write the coefficient J n as We shall suppose that the limit exists.From the restrictions on γ it follows that a 0; we shall first analyze the generic case a > 0. We can then drop the 1/(N − 1) term in the first factor under the radical in Eq. (6.2) and take the continuum limit of J n (after an obvious rescaling) as The integral in Eq. (5.2) is most easily computed through the change of variable ξ = cos 2 θ, which yields where is the complete elliptic integral of the first kind, and the modulus of the elliptic functions We thus have (dropping, for the sake of conciseness, the modulus k) and hence where we have used the fact that Using Eq. (4.15) we finally obtain the following closed-form expression for f (N, λ) when a is positive: Consider next the case in which the limit (6.3) vanishes, so that the term 1/(N − 1) in the first factor under the radical in Eq. (6.2) cannot be neglected.We now write where and the modulus of the elliptic integral is and therefore Proceeding as before we arrive at the following formula for the function f (N, λ) in Eq. (4.14): The latter limit is finite for λ = 0, in which case for large N we can write Thus when λ > 0 in the limit ε 1,2 → 0+ we have which coincides with the a → 0+ limit of Eq. (6.4).It is also straightforward to compute the a → ∞ limit of the function f (N, λ) in (6.4).Indeed, in this limit the modulus k = (1+a) − 1 2 tends to zero, so that K → π/2, F (ϕ) → ϕ, and therefore We thus obtain the asymptotic formula  where the first equality follows from Schmidt's decomposition and the second one is due to the chain's symmetry about its midpoint.On the other hand, for finite a the sextic chain is not symmetric about its midpoint, and thus neither its entanglement entropy nor the asymptotic approximation (4.14) thereof are invariant under λ → 1 − λ.
Our numerical simulations indicate that the asymptotic formula (4.14)-(4.15)does indeed provide an excellent approximation to the Rényi entanglement of the sextic chain for N 1 in the absence of a magnetic field and at half filling.For α < 1 this is illustrated by Fig. 1 (right), where we present the case of N = 400 spins for a = 10 −2 (for which the hopping amplitude is neither homogeneous nor approximately symmetric about the midpoint) and α = 1/4, 1/2, 3/4.Of course, in order to compare the asymptotic formula (4.14)-(4.15)with the exact (numerically computed) value of the entanglement entropy S α it is necessary to first determine the constant part c α in the former equation.We have simply estimated c α as the average value of the difference between S α and the leading term of its asymptotic approximation (4.14).Surprisingly, this value of c α coincides to a remarkable accuracy with the corresponding one for the homogeneous XX chain given  by Eq. (4.4) (see, e.g., Fig. 2 (left) for a = 10 −2 ).In fact, we have checked that this is also the case for several other values of the parameter a.
As mentioned in the previous section, in the homogeneous XX chain the o(1) term in the asymptotic formula (4.2) for the Rényi entanglement entropy with parameter α 1 is oscillatory and of order N −1/α (cf.Eqs.(4.7)-(4.8)).In particular, at half filling this term features parity oscillations with amplitude roughly proportional to [(N/π) sin(πL/N )] −1/α .We have checked that the behavior of the Rényi entanglement entropy of the sextic chain with parameter α 1 at half-filling and zero magnetic field is very similar, and in particular that its parity oscillations are reproduced with great accuracy by Eqs.(4.14)-(4.15).This can be seen, for instance, in Fig. 2 (right), where we compare the Rényi entanglement entropy with parameter α = 2 for a = 10 −2 and N = 400 spins with its asymptotic approximation (4.14).Remarkably, in all the cases we have analyzed the values of the parameters c α and µ α are very close to the corresponding ones for the homogeneous model, given by Eqs.(4.4) and (4.8).This suggests -as shall be further corroborated by the analysis of the Krawtchouk chain in the next section-that at half filling and in a vanishing (or constant) magnetic field the sextic chain is in the same universality class as the homogeneous XX chain.
The situation is markedly different in the presence of an inhomogeneous magnetic field and/or at arbitrary fillings.Indeed, our numerical calculations clearly indicate that in these cases the behavior of the Rényi entanglement entropy is not well described even to leading order by the conformal analogue of Eq. (4.5)-(4.6),namely 1) , (6.8)    with Consider, to begin with, the case β = 0 and k F = π/2, illustrated in Fig. 3 (left) for the Fermi momentum k F = π/4 and a = 10 −2 or a = 1.The fact that k F = π/2 is seen to have two main effects.In the first place, S α is now virtually zero for small L and N − L (for instance, if N = 400 and a = 10 −2 then S 2 is less than 10 −3 for 1 L 64 and 393 L 399, while for a = 1 we have S 2 < 10 −3 if 1 L 31 and 385 L 399).It is then clear that the entanglement entropy cannot be well approximated in this case by the concave function in the RHS of Eqs.(6.8)-(6.9).Moreover, for α 1 the parity oscillations of S α are much less regular than in the case of half filling, and their amplitude is not well reproduced by a simple formula like (4.14) (see, e.g., the main plot in Fig. 3 (left) for the case a = 10 −2 , β = 0, k F = π/4, N = 400 and α = 2).On the other hand, a rough approximation capturing only the average variation of S α with L can be obtained by restricting ourselves to the interval [L 1 + 1, L 2 − 1] in which S α differs significantly from zero, and replacing accordingly L and N respectively by L − L 1 and L 2 − L 1 .With these changes Eq. (6.9) becomes with The rough approximation (6.8)-(6.10) is represented by a red line in Fig. 3.
The situation is qualitatively similar in the presence of an external magnetic field given by Eq. (6.1) with β = 0, even in the case of half filling; see, e.g., Fig. 3 (right).More precisely, as seen in the latter figure, the length of the intervals in which S α is practically zero increases significantly with β.Moreover, when k F differs from π/2 the interval over which S α is appreciably different from zero is translated by an amount depending on β.As before, the pattern of the parity oscillations is much more involved than in the case B n = 0, k F = π/2 analyzed earlier, although the average variation of S α with L is still reproduced to a certain extent by the heuristic formula (6.8)-(6.10).
Note, finally, that the fact that the entropy of the block {0, . . ., L − 1} is negligibly small for L L 1 and L L 2 clearly indicates that the first L 1 and last N − L 2 spins are approximately in a product state, so that the chain's entanglement is almost entirely concentrated in the central block {L 1 , . . ., L 2 − 1}.It would certainly be of interest to understand how exactly this phenomenon arises as the external magnetic field is turned on, or the standard filling M/N = 1/2 is varied.

Definition and entanglement entropy
The Krawtchouk chain was introduced in Ref. [21] in connection with the family of discrete Krawtchouk polynomials.More precisely, the Krawtchouk polynomial K n (x; q, m) is defined as where 0 < q < 1 and m is a nonnegative integer (see, e.g., [31]).Here 2 F 1 denotes the standard hypergeometric function where (a) k := a(a + 1) is the usual (ascending) Pochhammer symbol.These polynomials satisfy the recursion relation The corresponding monic polynomials are given by and satisfy the normalized recursion relation with B n as before and The polynomial K m+1 (x; q, m) cannot be defined through Eq. (7.1), since (−m) m+1 = 0. On the other hand, we can use the previous recursion relation with n = m to define P m+1 by Using the previous formula and the recursion relation, after a long but straightforward calculation we obtain In view of the above, we define the Krawtchouk chain2 through the polynomials P n (x) with m = N − 1.In other words, the chain's parameters J n and B n are given by 2) In particular, the magnetic field is constant if and only if q = 1/2.Using the definition of the Krawtchouk polynomials we readily obtain the explicit formula . Thus in this case the single-particle energies ε k are simply the integers 0, 1, . . ., N − 1.This makes it possible to express the matrix elements Φ nk in closed form using Eqs.(3.3)- (3.4).Indeed, it is readily found that and therefore Dropping the irrelevant overall factor [q(1 − q)] 1/2 (N − 1) in Eq. ( 7.2) we easily obtain the following formula for the continuum limit of J n : Left: Rényi entanglement entropy S 2 for the Krawtchouk chain with q = 1/2 and N = 400 spins at half filling (blue crosses) compared to its asymptotic approximation (4.14)-(7.6)(red squares).The inset shows a similar plot for k F = π/4, compared to its heuristic approximation (6.8)-(7.7).Right: analogous plot for q = 1/4 at half filling (main plot) and for k F = π/4 (inset).
We thus have and hence the function f (N, λ) in Eq. (4.15) is simply given by As expected, this result coincides with the a → ∞ limit of the analogous function for the sextic chain.In other words, the asymptotic behavior of the entanglement entropy of the Krawtchouk chain with q = 1/2 at half filling, given by Eqs.(4.14)-(7.6),should be the same as for the sextic chain with β = 0 and a 1.This is shown in Fig. 4 (left) for N = 400 spins.For the same reason, at arbitrary fillings the heuristic asymptotic approximation to the entanglement entropy is given by Eq. (6.8) with where [L 1 +1, L 2 −1] is the interval over which S α is appreciably nonzero; see, e.g., the inset of Fig. 4 (left).Note also that when q = 1/2 we have J n = J N −n−2 and B n is constant, so that in this case the entanglement entropy is invariant under L → N − L, i.e., satisfies Eq. (6.7), for all fillings.On the other hand, for q = 1/2 the magnetic field strength B n is non-uniform, so that the entanglement entropy behaves much the same as for the sextic chain with a 1 and β = 0; see, e.g., Fig. 4 (right) for the case q = 1/4 and N = 400.The main difference, as seen in the right inset of the latter figure, is that in this case the entropy is close to zero only on an interval of the form [L 2 , N ].
Remark 2. The value of the subleading (constant) term c α in the asymptotic formula (4.14)-(7.6)for the entanglement entropy of the Krawtchouk chain with q = 1/2 at half filling is remarkably close to its counterpart c α,hom for the homogeneous chain (cf.Eq. (4.4)), even more so than in the case of the sextic chain discussed above.For instance, the difference |c α − c α,hom | is of the order of 10 −3 or less for 1/2 α 2. In fact, the values of c α in the latter range were obtained as in the previous example by taking the average of the differences S α − S α,app for all values of the block size L = 1, . . ., N − 1, where S α,app denotes the leading term in Eq. (4.14).The behavior of the latter differences, however, clearly suggests that in this case a more accurate estimate for c α is given by the average of S α − S α,app for the central block sizes L = N/2 and L = N/2 − 1, or more simply by see, e.g., Fig. 5 (left) for α = 1.The previous equation actually yields a value of c α much closer to c α,hom than the average of the differences S α − S α,app for all block sizes.For instance, for α = 1 the difference between c 1,hom and the value of c 1 computed from the previous formula is 1.1 • 10 −6 , compared to 3.6 • 10 −4 when c 1 is estimated by the average of S − S app over all block sizes.We have also studied how the difference c α,hom − c α varies as the number of spins increases from N = 400 to N = 600 (in increments of 10) for several values of α, where in view of the previous remark we have used Eq.(7.8) to estimate c α .As is apparent from Fig. 5 (right) for the case α = 1, the absolute value of this difference steadily decreases with N .We thus conjecture that for the Krawtchouk chain with q = 1/2 at half filling the constant term c α tends to c α,hom as the number of spins tends to infinity.A similar analysis for the sextic chain (with a = 10 −2 ) also shows a decrease in |c α,hom − c α | as N increases in the same range, although the value of this difference is about three orders of magnitude higher than in the case of the Krawtchouk chain.This different behavior could be explained by the fact that the sextic chain is not invariant under J n → J N −n−2 , as are the homogeneous and Krawtchouk chains.In fact, the latter two chains are the only XX spin chains of algebraic type (in the more general sense that the recursion coefficients a n are polynomial in n) with interactions J n invariant under n → N − n − 2 having a finite total conformal length L. Remark 3. The plots in Fig. 4 suggest that the entanglement entropy of the Krawtchouk chain is invariant under L → N − L not only for q = 1/2 (at arbitrary filling), but also for arbitrary q at half filling.That this is indeed the case can be deduced from a general property of the entanglement entropy, stemming from the fact that J n := J N −n−2 = J n and B n is linear in n.Indeed, setting B n = B 0 + bn, with B 0 and b independent of n, we have From the previous equations for J n and B n we immediately obtain the following relation for the corresponding polynomials P n : This is readily seen to imply that which yields the relation where in the first equality we have applied the well-known invariance of the entanglement entropy under complements in position space.On the other hand, from the energy-position duality of the entanglement entropy [12,32,33] it follows that S α is also invariant under complements in energy space.We thus obtain the relation In particular, this implies that at half filling S α is invariant under L → N − L, as claimed.
Of course, for the Krawtchouk chain with q = 1/2 we can combine Eqs.(6.7) and (7.9) to deduce that S α is also invariant under M → N − M (this also follows from a standard duality argument).Note, finally, that since the couplings of the sextic chain with a 1 are approximately symmetric under n → N − n − 2, and its magnetic field term B n is linear in n, Eq. (7.9) is approximately valid also in this case.

Ground state energy
In Ref. [20] it is conjectured that as N → ∞ the ground state energy E 0 (N ) of an inhomogeneous XX chain (2.1) with B n = 0 for all n behaves as where N = L/∆x, and c 0 , c B , v F are three constants (representing the bulk energy per site, the boundary energy and the Fermi velocity) which in the homogeneous case take the respective values 2/π, 4/π − 1, and 2 [34,35].In the absence of a magnetic field the ground state is the half-filled state (2.11), whose energy E 0 (N ) is given by Eq. (2.12).This quantity can be exactly computed for the Krawtchouk chain with q = 1/2, since its single-particle energies (after subtraction of the constant magnetic field B n = (N − 1)/2) are given by the formula We shall next compare this exact value for E 0 (N ) with its conjectured asymptotic expansion (7.10), which by Eqs.(7.2) and (7.5) reads in this case Using the exact value (7.11) for E 0 (N ) we deduce that in this case On the other hand, the leading asymptotic behavior as N → ∞ of the sum can be determined from the Euler-Maclaurin formula [36] 2Σ where g(x) := x(N − x).The remainder R 3 can be estimated as where ζ(s) denotes Riemann's zeta function.We also have and hence We thus conclude that in agreement with the right-hand side of Eq. (7.13) if we take c 0 = 2/π (as in the homogeneous case).
It can be shown that the higher-order corrections in the Euler-Maclaurin formula are all O(N 1/2 ), so that they cannot be used to compute c B and v F in closed form (this is essentially due to the fact that the derivatives of g(x) diverge at x = 0, N ).On the other hand, the parameter c B can be computed through the formula By evaluating the right-hand side for large values of N we have verified that this limit indeed exists, and that c B = 0.1323449 to seven decimal places.Note that this value is about one half of the corresponding one for the homogeneous case.
Using the previous estimate for the constant c B , we have studied the behavior of the difference N 2 /8 − 2Σ N /π − c B √ N − 1 for 1000 N 15000, obtaining very strong numerical evidence that it is of order N −1/2 instead of N −1 , as predicted by Eq. (7.13) (cf.Fig. 6).In other words, our results suggest that for the Krawtchouk chain (with q = 1/2) Eq. (7.12) should be replaced by For 0 < λ < 1 the numerator of the argument of the cosine tends to the finite limit F (arcsin(2λ − 1), 1) = arctanh(2λ − 1) as N → ∞ (i.e., ε → 0+), while the denominator tends to K(1) = +∞.Thus for sufficiently large N we can take the function f in Eq. (4.15) simply as up to lower-order terms in N .Note that f (N, λ) is invariant under λ → 1 − λ, which is consistent with the symmetry of the coupling (8.1) with respect to the chain's midpoint.
It can be shown that in the N → ∞ limit we have with κ = 0.1882264 . . ., so that f (N, λ) ∼ N log N in this limit.Thus when N → ∞ with λ fixed the entanglement entropy of the Lamé chain (8.1) at half filling should behave as It should also be noted that the logarithmic divergence of f (N, λ)/N as N → ∞ with λ fixed is due to the presence of double zeros of J(x) at both endpoints of the interval [0, L] in this limit, and can thus only happen when the polynomial p(ξ) in Eq. (5.1) is of degree four.
We have numerically checked that the asymptotic formula (4.14)-(8.3)still provides a reasonable approximation to the entanglement entropy at half filling in this case, though not as precise as for the sextic and Krawtchouk chains.In particular, the oscillating term proportional to f (N, λ) −1/α in Eq. (4.14) reproduces with acceptable accuracy the parity oscillations in S α present when α 1 (see, e.g., Fig. 7 (left) for α = 2 and N = 400 spins).As before, at Fermi momenta k F = π/2 the Rényi entanglement entropy S α is virtually zero over two intervals of the form [1, L 1 ] and [N − L 1 , N − 1], and the asymptotic formula (4.14)-(8.3)fails.However, proceeding as above we can derive a rough heuristic approximation to S α in the interval [L 1 + 1, N − L 1 − 1] using Eq.(6.8) with with λ eff given by Eq. (6.11) with Since the dependence on the number spins N of the asymptotic formula (4.14)-(8.3)for the entanglement entropy at half filling is nontrivial, it is also of interest to study the growth of S α (N, λ) with N for fixed values of λ.We have checked that the latter formula captures the behavior of S α (N, λ) with great accuracy, and in particular reproduces the parity oscillations that appear when α 1 (see, e.g., Fig. 7 (right) for the case λ = 1/2 and even N in the interval [500, 600]).Remarkably, although the two parameters c α and µ α appearing in Eq. (4.14) are fitted, they turn out to be of the same order of magnitude as their counterparts (4.4)-(4.8)for the homogeneous XX chain.Note, however, that in this case it should not be expected that the constant term c α tend to c α,hom as N → ∞, since the subleading term in the asymptotic expansion of S α is no longer constant but of the order of log(log N ) (cf.Eq. (8.4)).

Conclusions and outlook
In this work we study a large class of inhomogeneous XX spin chains whose squared couplings are a polynomial of degree at most four in the site index.This class includes some previously studied models related to classical Krawtchouk and dual Hahn polynomials [21], as well as the inhomogeneous XX chains related to QES models on the line classified in Ref. [19].We show how to exactly compute the leading term in the asymptotic expansion of the block entanglement entropy of these models (in a constant magnetic field at half filling) from their continuum limit, which coincides with the CFT of a massless Dirac fermion in a curved (1 + 1)-dimensional background [15,17,18].We next focus on three inhomogeneous chains with algebraic interactions, associated to the sextic oscillator QES potential, the Krawtchouk polynomials and the periodic quantum Lamé potential.We exploit the relation of XX chains with finite families of orthogonal polynomials to numerically compute the Rényi entanglement entropy of the latter chains for a large number of spins.When the Rényi parameter α is less than one we find that, as expected, the asymptotic formula reproduces with great accuracy the behavior of the entanglement entropy.On the other hand, for α 1 the Rényi entropy presents parity oscillations whose amplitude increases with α, as is known to be the case for the homogeneous XX chain [13].We show that these oscillations are reproduced with excellent accuracy by the Fagotti-Calabrese formula for the homogeneous chain [13], replacing the block's and the chain's length by their conformal counterparts.In fact, for the sextic and Krawtchouk chains (at half filling and in a constant magnetic field) we have found rather compelling numerical evidence that the subleading non-universal (constant) term in the asymptotic expansion of the entanglement entropy tends to its counterpart for the homogeneous XX chain as the number of spins tends to infinity.We conjecture that this is actually the case for the class of algebraic chains studied in this paper when the chain's conformal length is finite, i.e., when the squared couplings considered as functions of the site index n have no multiple real roots in the interval [0, N ].
All of the above results apply to the case of half filling and constant magnetic field, which is the one usually considered in the literature in the inhomogeneous case.In this work we have also studied in some detail the non-standard situation of arbitrary filling and/or inhomogeneous magnetic field.We have found that the main difference with the standard situation is that the block entanglement entropy vanishes when the block's length is either small or close to the chain's length.Thus at fillings other than one-half, or in an inhomogeneous magnetic field, the first few and last spins in the chain become disentangled from the rest.This is in fact one of the paper's main results, which certainly deserves further theoretical analysis.Another interesting feature of the non-standard case is that the oscillations of the entropy when α 1 are considerably more complex than in the standard one, and in particular are not well described by a modification of the Fagotti-Calabrese formula along the lines mentioned above.
One of the models studied in this paper, namely the Krawtchouk chain (cf.Section 7), has the rather unusual property that its single-particle energies at zero magnetic field can be exactly computed (they are simply the numbers −(N − 1)/2 + k, with 0 k N − 1).This of course makes it trivial to evaluate the ground-state energy in closed form for an arbitrary number of spins N .We have compared this exact result with the recently proposed asymptotic expansion in Ref. [20], finding that they match only to leading order.
The above results clearly suggest several avenues for future research that we shall now briefly outline.To begin with, we would like to find a theoretical justification of the fact that the constant term in the asymptotic expansion of the entanglement entropy of chains with algebraic interactions and finite conformal length (at half filling and in a constant magnetic field) seems to coincide with the analogous term for the homogeneous XX model in the limit of large N .An outstanding open problem of considerable interest is that of deriving an asymptotic formula for the leading behavior of the entanglement entropy in the non-standard scenario of arbitrary filling and/or inhomogeneous magnetic field.Our numerical results show that any such formula must necessarily vanish when the block length is either small or close to the chain's length.Likewise, it would also be of interest to find a formula describing the entropy's complex oscillations when the Rényi parameter is greater than or equal to one, akin to the Fagotti-Calabrese formula for the homogeneous chain.As mentioned in the Introduction, in the homogeneous case the behavior of the multiblock entanglement entropy has been thoroughly analyzed (see, e.g., [8-10, 12, 38]).Again, the generalization of some of these results to the inhomogeneous case would certainly be worth exploring.Finally, another problem suggested by the present work is to understand why maps p+ (z) into a positive multiple of p− (w) with k replaced by k, and the interval (−1, 1) to (1, 1/ k), it is easier in this case to perform the change of variable in the integral for p− (z).In this way we obtain II. p has two simple real and two complex conjugate roots In this case we can take A 1 B 1 < 0 and A 2 B 2 > 0 in Eq. (A.3).We can thus write (applying, if necessary, a dilation) with ν 0 ∈ R. We can also assume w.l.o.g. that A 2 , B 2 > 0, and set III. p has four simple complex roots.
In this case we have A i B i > 0 in Eq. (A.3).Moreover, since p must be positive on some open interval we can take w.l.o.g.A i , B i > 0. We can therefore write (up to a dilation) with ν, A > 0. Applying, if needed, a projective transformation z = 1/w, we can assume w.l.o.g. that A < 1, and thus set A = k 2 = 1 − k 2 with 0 < k < 1. Hence in this case p(ξ) can be reduced to the canonical form which is positive everywhere.Performing the change of variable z = tan θ we then obtain

Figure 2 .
Figure 2. Left: constant term c α in the asymptotic approximation (4.14)-(4.15) to the Rényi entanglement entropy S α of the sextic chain with a = 10 −2 and N = 400 spins for α = i/10 and 1 i 30 (blue crosses) compared to the corresponding constant c α for the homogeneous XX chain in Eq. (4.4) (solid red line).The inset shows the difference between c α and its counterpart (4.4) for the homogeneous chain in the range 1/2 α 2 at intervals of 1/10.Right: Rényi entanglement entropy S 2 for the sextic chain with a = 10 −2 and N = 400 spins at half filling and zero magnetic field (blue crosses) compared to its asymptotic approximation (4.14) (red squares).The inset shows a blow up of the range 150 L 250 in the latter plot.

1 Figure 5 .
Figure 5. Left: difference between the von Neumann entanglement entropy S of the Krawtchouk chain with q = 1/2, N = 400, k F = π/2 and its leading asymptotic approximation S app := (1/6) log f (N, λ) as a function of the block size L. (Only the range L = 1, . . ., N/2 has been represented, since in this case S α is symmetric under L → N −L.)The height of the red (resp.green dashed) horizontal line is the value of c 1 computed from Eq. (7.8) (resp.as the average of the differences S − S app for all block sizes).Right: Difference c 1,hom − c 1 for N = 400, . . ., 600 in increments of 10.