Computation of entanglement entropy in inhomogeneous free fermions chains by algebraic Bethe ansatz

The computation of the entanglement entropy for inhomogeneous free fermions chains based on q -Racah polynomials is considered. The eigenvalues of the truncated correlation matrix are obtained from the diagonalization of the associated Heun operator via the algebraic Bethe Ansatz. In the special case of chains based on dual q -Hahn polynomials, the eigenvectors and eigenvalues are expressed in terms of symmetric polynomials evaluated on the Bethe roots.


Introduction
The characterization of entanglement in many-body systems is motivated by its numerous applications in quantum information [1,2] and its role in describing quantum critical points [3].This endeavour is usually carried out in bipartite situations, where the amount of entanglement between a region and its complement is determined.While many techniques have been developed to perform this task [4], analytical results for entanglement entropy in large systems remain rare.
For spin chains and free fermions models, this problem reduces to diagonalizing a matrix referred to as the truncated correlation matrix [5].In cases where couplings are homogeneous, for example the XX spin chain, this matrix is Toeplitz or Toeplitz+Hankel and one can use the Fisher-Hartwig conjecture to compute the bipartite entanglement in the thermodynamic limit [6,7].For more general couplings and truncated correlation matrices, applying these techniques is not possible and different approaches are required.
Inhomogeneous fermionic chains associated to hypergeometric orthogonal polynomials of the Askey-Wilson scheme [8] are solvable and describe a wide variety of models.It was observed that their truncated correlation matrix admits a commuting tridiagonal matrix, identified as a Heun-Askey-Wilson operator [9][10][11][12][13].This suggests an interesting connection with the theory of integrable systems.Indeed, these operators arise in the transfer matrices associated to solutions of the reflection equations [14].They correspond to Hamiltonians of XXZ spin chains with specific boundary fields and have been shown to be dagonalizable via the algebraic Bethe anstaz [15] (other methods have been developed in [16,17]).They are also examples of the so-called homogeneous case in the context of the modified algebraic Bethe ansatz, which has been designed to deal with generic Heun operators [14,18,19] and diagonalize integrable models with arbitrary boundary conditions (see e.g.[18,[20][21][22][23][24]).
This paper applies the algebraic Bethe ansatz framework to investigate the spectrum of truncated correlation matrices of models associated to polynomials of the Askey-Wilson scheme.In particular, the eigenvalues of the truncated correlation matrix of free fermionic chains associated to dual q-Hahn polynomials are provided in terms of solutions of a set of Bethe equations.In section 2, we recall the definition of free fermions chains associated to q-Racah polynomials and diagonalize their Hamiltonians.In section 3, we discuss the problem of computing the entanglement entropy and introduce the truncated correlation matrix.In section 4, we exhibit a commuting tridiagonal matrix referred to as the algebraic Heun operator and diagonalize it via the algebraic Bethe ansatz.This yields a set of relations known as Bethe equations.The eigenvalues of the truncated correlation matrix are then given in terms of roots of these equations.The associated TQ-relation and the thermodynamic limit are briefly discussed in section 5.

The model
Let us consider the following free fermions inhomogeneous Hamiltonian with nearest-neighbour interaction J n and magnetic field µ n , where c n and c † n are fermionic annihilation and creation operators satisfying For convenience, we enumerate the sites of the lattice from 0 to N .This model is equivalent to an inhomogeneous XX spin chain.Indeed, the Jordan-Wigner transformation allows to rewrite the canonical relations of the creation and annihilation operators (2) and the Hamiltonian (1) in terms of spin-1/2 operators, We are interested in the case where the coupling parameters J n and the local magnetic field µ n are constructed from the recurrence coefficients of the q-Racah polynomials [8]: where ε = ±1 and A n , C n are defined by The choice of such inhomogeneous interactions and magnetic fields yields analytical results for the spectrum, as shown below.It also describes a large class of models thanks to the presence of various parameters.Indeed, the constants A n , C n and µ n depend on the five real parameters q, α, β, γ and δ, restricted only by the requirement that For instance, as shown in figure 1, we can get couplings J n which are monotone in n or peaking at a certain value.Taking q < 0 also gives models with oscillating couplings, reminiscent of alternating spin chains [25].
Figure 1: Inhomogeneous free fermions chains of length N = 10, based on q-Racah polynomials, for different parameters (q, α, β, γ, δ).The vertices and edges represent respectively the sites and the couplings.The color of the edges indicates the magnitude of J n , i.e. the strength of these couplings.Darker color is associated to stronger couplings.

Diagonalization of the Hamiltonian
In order to diagonalize H, it is convenient to rewrite this operator in matrix form as where A is the hermitian (N + 1) × (N + 1) tridiagonal matrix given by with the convention J N = J −1 = 0.The vectors {|0〉, |1〉, . . ., |N 〉} are naturally associated to sites in the chain and give the canonical orthonormal basis of N +1 .They will be referred to as elements of the position basis.The spectral problem for A reads where Knowing that the entries of A are the recurrence coefficients of the q-Racah polynomials, one deduces (see eq. ( 76)) that its eigenvalues ω k are The wavefunctions φ n (ω k ) = 〈ω k |n〉 are given in terms of q-Racah polynomials R n (ω k ) [8]: The definition of R n (ω k ) and the normalisation factors W k are given in appendix A. The latter are chosen such that the wavefunctions φ n (ω k ) are orthonormal i.e.N k=0 From these wavefunctions, we can define new pairs of fermionic creation and annihilation operators in terms of which the Hamiltonian is diagonal: where Note that the operators c † k and c k are associated to the single particle excitations of the system, with energies given by the spectrum of the matrix A. One may further observe that these energies are invariant under arbitrary transformations of α, β and under This is not true of the coupling parameters J n and local magnetic field µ n which depend nontrivially on (q, α, β, γ, δ).Important properties characterizing these systems, like the entanglement entropy in the ground state, should thus depend on these parameters.

Entanglement entropy
Entanglement in a multipartite system A∪ B is measured by the entanglement entropy S A , defined as where A is a subsystem of A ∪ B with a reduced density matrix ρ A given by the trace over the degrees of freedom in B, In the following, we take A to be the first L + 1 sites of the inhomogeneous free fermionic chain introduced in the previous section.The states considered are obtained by filling up the first K + 1 single particle states, taken as the Fermi sea, where |0〉〉 is the vacuum state annihilated by all operators ck .For ω k monotone in k, |Ω〉〉 describes the ground state of Hamiltonians obtained as affine transformations of (10).In other words, it gives the state for which the single particle excitations with negative energy are filled.
As observed in [5], computing the entanglement entropy S A of free fermions can be done by diagonalizing the truncated correlation matrix.Indeed, it is known that [26] where the coefficients c are the eigenvalues of the (L + 1) × (L + 1) matrix C with entries C nm given by the 2-point correlation functions, This is a submatrix of the complete correlation matrix of the ground state C, i.e.
where π A is the projector onto the vector space associated to sites of subsystem A, The computation of the entanglement entropy is thus reduced to determining the eigenvalues c .With the help of the algebraic Heun operators, we shall see that the spectral problem for the truncated correlation matrix can be treated in the algebraic Bethe ansatz framework.

Algebraic Heun operator
In this section, we introduce a tridiagonal matrix that commutes with the truncated correlation matrix.To do so, we define an operator A * which is diagonal in the position basis Using the difference relation of the q-Racah polynomials and the expression (15), one finds the tridiagonal action of A * on the eigenbasis of A, The coefficients Jk and μk are given in appendix A. The Heun operator T is then defined as and has the property of commuting with both the projector π A and the complete correlation matrix This is shown easily by considering the commutators [T, π A ] in the position basis and [T, C] in the energy basis.Given relation (25), T also commutes with the truncated correlation C and thus share with it a common set of eigenvectors.This is a crucial observation, in particular because the Heun operator T can be identified in the transfer matrix of integrable models and can hence be diagonalized via the algebraic Bethe ansatz [14,27].

Algebraic Bethe ansatz
The matrices A and A * give a representation of the Askey-Wilson algebra [28,29]: where I is the N + 1 × N + 1 identity matrix and the constants ξ, χ, χ * , η and η * can be expressed in terms of the parameters in the Hamiltonian: η = (q − 1) 2 (q + 1)(αγ(βδ + δ + 1) + αβδ + γδ(βδ η * = (q − 1) 2 (q + 1) The so-called dynamical operators can be defined in terms of the generators of this algebra: and The functions f 1 (u, m) and f 2 (u, m) are given in the appendix.These operators verify where τ = q αβγδ .The functions f (u, v), g(u, v, m) and w(u, v, m) are given in appendix A. Relations (39)-(40) were verified using directly the Askey-Wilson relations (31)-(32).This is similar to the method used in [18,19] and distinct from the approach based on R and K matrices [14,15].The Heun operator (29) can be expressed in terms of A(u, L) and A(τu −1 , L) as where Next, let us consider the vectors |ū〉 defined as where Note that relation (39) implies that B(ū, L) does not depend on the ordering of the variables u i .
Since the vectors |ū〉 are obtained by applying L times a tridiagonal matrix on the vector |0〉, they are contained in the vector space spanned by {|0〉 , |1〉 , . . ., |L〉}.As such, they are eigenvectors of π A with eigenvalue 1, The aim is to show that for specific parameters ū, these vectors are also eigenvectors of T .This requires two results.The first is the following relation between the dynamical operator A(u, m) and product of dynamical operators B(ū, L): where This relation is obtained by computing the terms i = 1 and by using the symmetry in the indices u i induced by relation (39).The second required result is the action of A(u, 0) on the vector |0〉.
From the definition of A(u, 0) in terms of A and A * , and the action ( 11)-( 27) of these operators on the position basis, it follows that where In particular, we note that this action is diagonal.This feature shows that the modified algebraic Bethe ansatz is not necessary and that the model we deal with corresponds to the particular case developed in [15].From this observation and relation (46), it follows that where and Thus, for a set of parameters ū verifying E i (u, ū) = 0, the vector |ū〉 = B(ū, L) |0〉 is an eigenvector of T with eigenvalues Λ(ū), i.e.
Since the Heun operator T does not depend on the parameter u, the same is true of its eigenvalues Λ(ū).In particular, we find by evaluating (51) at u = 0 that the eigenvalues can be expressed as where Factorizing terms in the variable u in the conditions E i (u, ū) = 0, these reduce to the following conditions on ū, referred to as the Bethe equations, To keep the notation simple, ū = {u 1 , u 2 , . . ., u L } will refer from now on to Bethe roots, i.e. to solutions of the set of equations (55).

Diagonalization of the truncated correlation matrix
Since the Heun operator T commutes with the truncated correlation matrix and is non-degenerate, its eigenvectors |ū〉 also diagonalize the matrix C, To obtain an explicit expression for the eigenvalues c(ū) in terms of the parameters ū, we observe that the action of B(u, m) in the position basis is tridiagonal and given by and Z n,m can be computed directly from the definition (38) and the action of A and A * in the position basis (see appendix A).In the case where β = 0, i.e. the dual q-Hahn special case of the q-Racah polynomials [8], these coefficients simplify greatly, In particular, B(u, m) becomes a raising operator in the sense that 〈n − 1| B(u, m) |n〉 = 0.This allows to compute the wavefunction 〈n|ū〉 of Bethe vectors: where Ū = {U 1 , U 2 , . . ., U L } with U i = qu −2 i .The terms S r ( Ū) are symmetric polynomials of degree r in the variables U i defined by Then, one can use the representation of the truncated correlation matrix in the position basis (24) to obtain a formula for its eigenvalues in terms of Bethe roots.For any n ∈ {0, 1, . . .L}, we find where This is valid for parameters ū which are solutions of the Bethe equations (55).For β = 0, these equations reduce to

TQ-relations and thermodynamic limit
Equations ( 54) and (61) show that the spectra of the Heun operator and of the truncated correlation matrix can be obtained by solving Bethe equations.An alternative approach is given by interpreting the expression (51) for the eigenvalues of T as a q-difference equation.In the case β = 0, (51) can indeed be rewritten as where p(U) is the following polynomial in the variable U, and Q(U) is a polynomial of degree L, the zeros U i of which are expressed in terms of entries of a Bethe root ū = {u 1 , u 2 , . . .u L }: Thus, one can use the zeros of polynomial solutions of equation ( 64) to identify Bethe roots.This equation is referred to as the TQ-relation in the literature.
Let us now further fix1 δ = 0, γ ∈ [0, 1] and q < 1. Inserting the r.h.s of (66) in equation (64) yields a three term recurrence relation for the symmetric polynomials S n ( Ū): where σ n = (q + 1)q −n + (q + 1)q n − 2(q + 1) (68) ε n = −αγ(q + 1) q −K + q −L + αγ(q + 1) q K+L q n + αγ(q + 1)q −n .(70) In the thermodynamic limit N → ∞, the parameter α = q −N −1 , with 0 < q < 1, goes to infinity and (67) becomes effectively a two term recurrence with solution The condition that S −1 ( Ū) = 0 then requires the eigenvalues Λ(ū) of T to take certain values The spectrum of the Heun operator given by this approximation is compared to spectra found using other methods in Table 1.One notes that the values match up to two digits at N = 49.This suggests that exact asymptotic results may be obtainable in the thermodynamic limit.-11957.8-11583.9

Solutions of
Table 1: Eigenvalues of the Heun operator (N = 49, L = 9, K = 24, q = 0.8, α = q −N −1 , β = 0, γ = 0.5, δ = 0) obtained by three methods.The first column are zeros of S −1 ( Ū) seen as a polynomial of degree L + 1 in Λ(ū).The polynomial was obtained by solving the three term recurrence (67).The second column corresponds to the approximation (72) of the spectrum found in the thermodynamic limit.The third is the result of diagonalizing T using scipy's linear algebra package [30].

Conclusion
Computing bipartite entanglement for free fermionic chains amounts to determining the spectrum of a truncated correlation matrix.For systems associated to q-Racah polynomials, it has been shown how this matrix can be diagonalized via the algebraic Bethe ansatz.In particular, its eigenvalues and eigenvectors have been given in terms of solutions of Bethe equations.The associated Bethe roots were also found to be related to zeros of polynomial solutions of a q-difference equation, referred to as the TQ-relation.This led to an approximate expression for the eigenvalues of the commuting tridiagonal matrix in the case δ = 0 and N → ∞.While these results do not provide an explicit formula for the bipartite entanglement, it establishes a clear connection between a central problem in quantum many-body physics and a set of tools coming from the study of integrable models.Future research should thus be directed toward investigating, notably in their thermodynamic limit, the solutions of the Bethe equations and TQ-relation that were found.Derivation of asymptotic expressions for these would provide the groundwork necessary to analyse the interplay between coupling inhomogeneities in free fermions chains and the presence of entanglement in the ground state.

Submission
Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada.