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

Pierre-Antoine Bernard, Gauvain Carcone, Nicolas Crampé, Luc Vinet

► To cite this version:

Pierre-Antoine Bernard, Gauvain Carcone, Nicolas Crampé, Luc Vinet. Computation of entanglement entropy in inhomogeneous free fermions chains by algebraic Bethe ansatz. 34th International Colloquium on Group Theoretical Methods in Physics (ICGTMP) - Group 34, Jul 2022, Strasbourg, France. pp.018, 10.21468/SciPostPhysProc.14.018. hal-05142980

HAL Id: hal-05142980

https://hal.science/hal-05142980v1

Submitted on 3 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Computation of entanglement entropy in inhomogeneous free fermions chains by algebraic Bethe ansatz

Pierre-Antoine Bernard¹*, Gauvain Carcone¹, Nicolas Crampé² and Luc vinet¹,³

1 Centre de recherches mathématiques, Université de Montréal, P.O. Box 6128, Centre-ville
Station, Montréal (Québec), H3C 3J7, Canada

2 Institut Denis-Poisson CNRS/UMR 7013 - Université de Tours - Université d'Orléans, Parc de
Grandmont, 37200 Tours, France

3 IVADO, 6666 Rue Saint-Urbain, Montréal (Québec), H2S 3H1, Canada

* bernardpierreantoine@outlook.com

December 21, 2022

Group
ICGTMP

34th International Colloquium on Group Theoretical Methods in Physics
Strasbourg, 18-22 July 2022

doi:10.21468/SciPostPhysProc.?

Abstract

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.

Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>The model</td><td>3</td></tr><tr><td>2.1</td><td>Diagonalization of the Hamiltonian</td><td>4</td></tr><tr><td>3</td><td>Entanglement entropy</td><td>5</td></tr><tr><td>4</td><td>Algebraic Heun operator</td><td>6</td></tr><tr><td>4.1</td><td>Algebraic Bethe ansatz</td><td>7</td></tr><tr><td>4.2</td><td>Diagonalization of the truncated correlation matrix</td><td>10</td></tr><tr><td>5</td><td>TQ-relations and thermodynamic limit</td><td>11</td></tr><tr><td>6</td><td>Conclusion</td><td>12</td></tr><tr><td>A</td><td>Appendix</td><td>13</td></tr></table>
A.1 $q$-Racah polynomials 13

A.2 Functions in the algebraic Bethe ansatz 14

References 14

# 1 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–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 diagonalizable via the algebraic Bethe ansatz [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–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.
## 2 The model

Let us consider the following free fermions inhomogeneous Hamiltonian with nearest-neighbour interaction $J_n$ and magnetic field $\mu_n$,

$$ \hat{\mathcal{H}} = \sum_{n=0}^{N-1} (J_n c_n^\dagger c_{n+1} + J_n c_{n+1}^\dagger c_n) - \sum_{n=0}^{N} \mu_n c_n^\dagger c_n, \quad (1) $$

where $c_n$ and $c_n^\dagger$ are fermionic annihilation and creation operators satisfying

$$ \{c_m^\dagger, c_n^\dagger\} = \{c_m, c_n\} = 0, \quad \{c_m^\dagger, c_n\} = \delta_{m,n}. \quad (2) $$

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

$$ c_n^\dagger = \sigma_0^z \sigma_1^z \cdots \sigma_{n-1}^z \sigma_n^+, \quad c_n = \sigma_0^z \sigma_1^z \cdots \sigma_{n-1}^z \sigma_n^-, \quad (3) $$

allows to rewrite the canonical relations of the creation and annihilation operators (2) and the Hamiltonian (1) in terms of spin-1/2 operators,

$$ \hat{\mathcal{H}} = -\frac{1}{2} \sum_{n=0}^{N-1} J_n (\sigma_n^x \sigma_{n+1}^x + \sigma_n^y \sigma_{n+1}^y) - \frac{1}{2} \sum_{n=0}^{N} \mu_n (1 + \sigma_n^z). \quad (4) $$

We are interested in the case where the coupling parameters $J_n$ and the local magnetic field $\mu_n$ are constructed from the recurrence coefficients of the q-Racah polynomials [8]:

$$ J_n = \epsilon \sqrt{A_n C_{n+1}}, \quad (5) $$

$$ \mu_n = A_n + C_n - 1 - \gamma \delta q, \quad (6) $$

where $\epsilon = \pm 1$ and $A_n, C_n$ are defined by

$$ A_n = \frac{(\alpha q^{n+1} - 1)(\gamma q^{n+1} - 1)(\alpha \beta q^{n+1} - 1)(\beta \delta q^{n+1} - 1)}{(1 - \alpha \beta q^{2n+1})(1 - \alpha \beta q^{2n+2})} \quad (7) $$

$$ C_n = \frac{(\beta q^n - 1)(\alpha q^n - \delta)(\alpha \beta q^n - \gamma)(q^{n+1} - q)}{(1 - \alpha \beta q^{2n})(1 - \alpha \beta q^{2n+1})}. \quad (8) $$

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 $\mu_n$ depend on the five real parameters $q, \alpha, \beta, \gamma$ and $\delta$, restricted only by the requirement that

$$ A_n C_{n+1} > 0, \quad J_n = 0. \quad (9) $$

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><img src="image_2.png" /><figcaption>Figure 1: Inhomogeneous free fermions chains of length $N = 10$, based on q-Racah polynomials, for different parameters $(q, \alpha, \beta, \gamma, \delta)$. 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.</figcaption></figure>

## 2.1 Diagonalization of the Hamiltonian

In order to diagonalize $\hat{\mathcal{H}}$, it is convenient to rewrite this operator in matrix form as

$$ \hat{\mathcal{H}} = (c_0^\dagger, \dots, c_N^\dagger) A \begin{pmatrix} c_0 \\ \vdots \\ c_N \end{pmatrix}, \qquad (10) $$

where $A$ is the hermitian $(N + 1) \times (N + 1)$ tridiagonal matrix given by

$$ A = \sum_{n=0}^{N} (J_n |n\rangle \langle n+1| - \mu_n |n\rangle \langle n| + J_n |n+1\rangle \langle n|), \qquad (11) $$

with the convention $J_N = J_{-1} = 0$. The vectors $\{|0\rangle, |1\rangle, \dots, |N\rangle\}$ are naturally associated to sites in the chain and give the canonical orthonormal basis of $\mathbb{C}^{N+1}$. They will be referred to as elements of the position basis. The spectral problem for $A$ reads

$$ A|\omega_k\rangle = \omega_k|\omega_k\rangle, \qquad (12) $$

where

$$ |\omega_k\rangle = \sum_{n=0}^{N} \phi_n(\omega_k) |n\rangle. \qquad (13) $$

Knowing that the entries of $A$ are the recurrence coefficients of the q-Racah polynomials, one deduces (see eq. (76)) that its eigenvalues $\omega_k$ are

$$ \omega_k = q^{-k} + \gamma \delta q^{k+1}. \qquad (14) $$
The wavefunctions $\phi_n(\omega_k) = \langle\omega_k|n\rangle$ are given in terms of q-Racah polynomials $R_n(\omega_k)$ [8]:

$$ \phi_n(\omega_k) = \epsilon^n \sqrt{W_k} \prod_{j=1}^{n} \sqrt{\frac{A_{j-1}}{C_j}} R_n(\omega_k). \qquad (15) $$

The definition of $R_n(\omega_k)$ and the normalisation factors $W_k$ are given in appendix A. The latter are chosen such that the wavefunctions $\phi_n(\omega_k)$ are orthonormal i.e.

$$ \sum_{k=0}^{N} \phi_n(\omega_k) \phi_m(\omega_k) = \delta_{n,m} \quad \text{and} \quad \sum_{n=0}^{N} \phi_n(\omega_k) \phi_m(\omega_k') = \delta_{k,k'}. \qquad (16) $$

From these wavefunctions, we can define new pairs of fermionic creation and annihilation operators in terms of which the Hamiltonian is diagonal:

$$ \hat{\mathcal{H}} = \sum_{k=0}^{N} \omega_k \tilde{c}_k^{\dagger} \tilde{c}_k, \quad k \in \{0, 1, \dots, N\} \qquad (17) $$

where

$$ \tilde{c}_k = \sum_{n=0}^{N} \phi_n(\omega_k) c_n, \quad \tilde{c}_k^{\dagger} = \sum_{n=0}^{N} \phi_n(\omega_k) c_n^{\dagger}. \qquad (18) $$

Note that the operators $\tilde{c}_k^\dagger$ 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 $\alpha, \beta$ and under

$$ \delta \rightarrow \delta\kappa, \quad \gamma \rightarrow \gamma\kappa^{-1}, \quad \kappa \in \mathbb{R}. \qquad (19) $$

This is not true of the coupling parameters $J_n$ and local magnetic field $\mu_n$ which depend non-trivially on $(q, \alpha, \beta, \gamma, \delta)$. Important properties characterizing these systems, like the entanglement entropy in the ground state, should thus depend on these parameters.

### 3 Entanglement entropy

Entanglement in a multipartite system $A \cup B$ is measured by the entanglement entropy $S_A$, defined as

$$ S_A = -\mathrm{tr}_A(\rho_A \ln \rho_A), \qquad (20) $$

where A is a subsystem of $A \cup B$ with a reduced density matrix $\rho_A$ given by the trace over the degrees of freedom in B,

$$ \rho_A = \mathrm{tr}_B |\Omega\rangle\langle\Omega|. \qquad (21) $$

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,

$$ |\Omega\rangle\rangle = \prod_{k=0}^{K} \tilde{c}_{k}^{\dagger} |0\rangle\rangle, \qquad (22) $$
where $|0\rangle\rangle$ is the vacuum state annihilated by all operators $\tilde{c}_k$. For $\omega_k$ monotone in $k$, $|\Omega\rangle\rangle$ 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]

$$S_A = -\sum_{\ell} c_{\ell} \ln c_{\ell} + (1 - c_{\ell}) \ln(1 - c_{\ell}), \qquad (23)$$

where the coefficients $c_\ell$ are the eigenvalues of the $(L+1) \times (L+1)$ matrix $C$ with entries $C_{nm}$ given by the 2-point correlation functions,

$$C_{nm} = \langle\langle\Omega|c_n^\dagger c_m|\Omega\rangle\rangle = \sum_{k=0}^{K} \phi_n(\omega_k) \phi_m(\omega_k), \quad n, m \in \{0, 1, \dots, L+1\}. \qquad (24)$$

This is a submatrix of the complete correlation matrix of the ground state $\hat{C}$, i.e.

$$C = \pi_A \hat{C} \pi_A, \quad \hat{C} = \sum_{k=0}^{K} |\omega_k\rangle \langle\omega_k|, \qquad (25)$$

where $\pi_A$ is the projector onto the vector space associated to sites of subsystem A,

$$\pi_A = \sum_{n=0}^{L} |n\rangle \langle n|. \qquad (26)$$

The computation of the entanglement entropy is thus reduced to determining the eigenvalues $c_\ell$. 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.

# 4 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

$$A^* |n\rangle = \lambda_n |n\rangle, \quad \lambda_n = q^{-n} + \alpha\beta q^{n+1}. \qquad (27)$$

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$,

$$A^* |\omega_k\rangle = \bar{J}_k |\omega_{k+1}\rangle - \bar{\mu}_k |\omega_k\rangle + \bar{J}_{k-1} |\omega_{k-1}\rangle. \qquad (28)$$

The coefficients $\bar{J}_k$ and $\bar{\mu}_k$ are given in appendix A. The Heun operator $T$ is then defined as

$$T = \{A, A^*\} - (\lambda_L + \lambda_{L+1})A - (\omega_K + \omega_{K+1})A^*, \qquad (29)$$

and has the property of commuting with both the projector $\pi_A$ and the complete correlation matrix $\hat{C}$,

$$[T, \pi_A] = [T, \hat{C}] = 0. \qquad (30)$$
This is shown easily by considering the commutators $[T, \pi_A]$ in the position basis and $[T, \hat{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].

## 4.1 Algebraic Bethe ansatz

The matrices $A$ and $A^*$ give a representation of the Askey-Wilson algebra [28, 29]:

$$AAA^* - \left(q + \frac{1}{q}\right)AA^*A + A^*AA = \xi A + \chi A^* + \eta I, \quad (31)$$

$$A^*A^*A - \left(q + \frac{1}{q}\right)A^*AA^* + AA^*A^* = \chi^*A + \xi A^* + \eta^*I, \quad (32)$$

where $I$ is the $N+1 \times N+1$ identity matrix and the constants $\xi, \chi, \chi^*, \eta$ and $\eta^*$ can be expressed in terms of the parameters in the Hamiltonian:

$$\chi = -\frac{\gamma\delta(q^2-1)^2}{q}, \quad \chi^* = -\frac{\alpha\beta(q^2-1)^2}{q}, \quad (33)$$

$$\xi = -(q-1)^2(\alpha(\beta\delta + \beta + \gamma + 1) + \gamma(\beta\delta + \delta + 1) + \beta\delta), \quad (34)$$

$$\eta = (q-1)^2(q+1)(\alpha\gamma(\beta\delta + \delta + 1) + \alpha\beta\delta + \gamma\delta(\beta\delta + \beta + \gamma + 1)), \quad (35)$$

$$\eta^* = (q-1)^2(q+1)(\alpha^2\beta + \alpha(\beta^2\delta + \beta(\gamma+1)(\delta+1) + \gamma) + \beta\gamma\delta). \quad (36)$$

The so-called dynamical operators can be defined in terms of the generators of this algebra:

$$A(u, m) = \frac{q^{-2L}}{(\alpha\beta q^{2m+1} - 1)} \left( \frac{q^{m+1}\{A, A^*\}}{(q+1)} - (\alpha\beta q^{2m+2} + 1)A - \frac{(q^{2m+2} + \gamma\delta u^4)}{u^2}A^* \right) + f_1(u, m)I, \quad (37)$$

and

$$\mathcal{B}(u, m) = \frac{\alpha\beta q^{m+2} + q^{-m-1}}{2(q+1)} \{A, A^*\} - \frac{q^{-m-1} - \alpha\beta q^{m+2}}{2(1-q)} [A, A^*] \\ - \alpha\beta(q+1)A - \frac{(q + \alpha\beta\gamma\delta u^4)}{u^2}A^* + f_2(u, m)I. \quad (38)$$

The functions $f_1(u, m)$ and $f_2(u, m)$ are given in the appendix. These operators verify

$$\mathcal{B}(u, m+1)\mathcal{B}(v, m) = \mathcal{B}(v, m+1)\mathcal{B}(u, m), \quad (39)$$

$$\mathcal{A}(u, m+1)\mathcal{B}(v, m) = f(u, v)\mathcal{B}(v, m)\mathcal{A}(u, m) + g(u, v, m)\mathcal{B}(u, m)\mathcal{A}(v, m) \\ + w(u, v, m)\mathcal{B}(u, m)\mathcal{A}(\tau v^{-1}, m), \quad (40)$$
where $\tau = \sqrt{\frac{q}{\alpha\beta\gamma\delta}}$. 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 $\mathcal{A}(u, L)$ and $\mathcal{A}(\tau u^{-1}, L)$ as

$$T = r(u)\mathcal{A}(u, L) + r(\tau u^{-1})\mathcal{A}(\tau u^{-1}, L) - (r(u)f_1(u, L) + r(\tau u^{-1})f_1(\tau u^{-1}, L))\mathcal{I}, \quad (41)$$

where

$$r(u) = \frac{q^L(q+1)}{\alpha\beta\gamma\delta u^4 - q} (\alpha^2\beta^2\gamma\delta u^4 q^{2L} + 1 - (\gamma\delta q^{K+1} + q^{-K-1})\alpha\beta u^2 q^L). \quad (42)$$

Next, let us consider the vectors $|\bar{u}\rangle$ defined as

$$|\bar{u}\rangle = \mathbf{B}(\bar{u}, L)|0\rangle, \quad \bar{u} = \{u_1, u_2, \dots, u_L\}, \quad (43)$$

where

$$\mathbf{B}(\bar{u}, L) = \mathcal{B}(u_1, L-1)\mathcal{B}(u_2, L-2)\dots\mathcal{B}(u_L, 0). \quad (44)$$

Note that relation (39) implies that $\mathbf{B}(\bar{u}, L)$ does not depend on the ordering of the variables $u_i$. Since the vectors $|\bar{u}\rangle$ are obtained by applying $L$ times a tridiagonal matrix on the vector $|0\rangle$, they are contained in the vector space spanned by $\{|0\rangle, |1\rangle, \dots, |L\rangle\}$. As such, they are eigenvectors of $\pi_A$ with eigenvalue 1,

$$\pi_A |\bar{u}\rangle = \sum_{i=0}^{L} |i\rangle \langle i|\bar{u}\rangle = |\bar{u}\rangle. \quad (45)$$

The aim is to show that for specific parameters $\bar{u}$, these vectors are also eigenvectors of $T$. This requires two results. The first is the following relation between the dynamical operator $\mathcal{A}(u, m)$ and product of dynamical operators $\mathbf{B}(\bar{u}, L)$:

$$\begin{aligned} \mathcal{A}(u, m)\mathbf{B}(\bar{u}, L) = & \prod_{i=1}^{L} f(u, u_i)\mathbf{B}(\bar{u}, L)\mathcal{A}(u, m-L) \\ & + \sum_{i=1}^{L} g(u, u_i, m-1)\prod_{\substack{j=1 \\ i \neq j}}^{L} f(u_i, u_j)\mathbf{B}(\bar{u}_{\neq i}, u, L)\mathcal{A}(u_i, m-L) \\ & + \sum_{i=1}^{L} w(u, u_i, m-1)\prod_{\substack{j=1 \\ i \neq j}}^{L} f(\tau u_i^{-1}, u_j)\mathbf{B}(\bar{u}_{\neq i}, u, L)\mathcal{A}(\tau u_i^{-1}, m-L), \end{aligned} \quad (46)$$

where

$$\mathbf{B}(\bar{u}_{\neq i}, u, L) = \mathcal{B}(u_1, L-1)\mathcal{B}(u_2, L-2)\dots\mathcal{B}(u, m-i)\dots\mathcal{B}(u_L, 0). \quad (47)$$

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 $\mathcal{A}(u, 0)$ on the vector $|0\rangle$. From the definition of $\mathcal{A}(u, 0)$ in terms of $\mathbf{A}$ and $\mathbf{A}^*$, and the action (11)-(27) of these operators on the position basis, it follows that

$$\mathcal{A}(u, 0)|0\rangle = a(u)|0\rangle, \quad (48)$$

where

$$a(u) = \left( \frac{q^{-2L}}{(1 - \alpha\beta q)} \left( \frac{2\mu_0\lambda_0 q}{(q+1)} - \mu_0(\alpha\beta q^2 + 1) + \frac{q^2\lambda_0}{u^2} + \gamma\delta u^2 \lambda_0 \right) + f_1(u, 0) \right). \quad (49)$$
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

$$ T |\bar{u}\rangle = \Lambda(\bar{u}) |\bar{u}\rangle + \sum_{i=1}^{L} E_i(u, \bar{u}) \mathbf{B}(\bar{u}_{\neq i}, u, L) |0\rangle, \quad (50) $$

where

$$ \Lambda(\bar{u}) = r(u)a(u) \prod_{i=1}^{L} f(u, u_i) + r(\tau u^{-1})a(\tau u^{-1}) \prod_{i=1}^{L} f(\tau u^{-1}, u_i) \\ - (r(u)f_1(u, L) + r(\tau u^{-1})f_1(\tau u^{-1}, L)), \quad (51) $$

and

$$ E_i(u, \bar{u}) = (r(u)g(u, u_i, L-1) + r(\tau u^{-1})g(\tau u^{-1}, u_i, L-1))a(u_i) \prod_{\substack{j=1 \\ i \neq j}}^{L} f(u_i, u_j) \\ + (r(\tau u^{-1})w(\tau u^{-1}, u_i, L-1) + r(u)w(u, u_i, L-1))a(\tau u_i^{-1}) \prod_{\substack{j=1 \\ i \neq j}}^{L} f(\tau u_i^{-1}, u_j). \quad (52) $$

Thus, for a set of parameters $\bar{u}$ verifying $E_i(u, \bar{u}) = 0$, the vector $|\bar{u}\rangle = \mathbf{B}(\bar{u}, L)|0\rangle$ is an eigenvector of $T$ with eigenvalues $\Lambda(\bar{u})$, i.e.

$$ T |\bar{u}\rangle = \Lambda(\bar{u}) |\bar{u}\rangle. \quad (53) $$

Since the Heun operator $T$ does not depend on the parameter $u$, the same is true of its eigenvalues $\Lambda(\bar{u})$. In particular, we find by evaluating (51) at $u = 0$ that the eigenvalues can be expressed as

$$ \Lambda(\bar{u}) = -(\omega_K + \omega_{K+1})\lambda_L - \frac{\mu_0(q-1)}{q} \left( \alpha\beta q^2 + 1 + \frac{4\alpha\beta q^2}{\alpha\beta q - 1} \right) - \frac{(q^2-1)^2}{q(q+1)} \left( \sum_{i=1}^{L} U_i \right), \quad (54) $$

where $U_i = \frac{q}{u_i^2} + \alpha\beta\gamma\delta u_i^2$. Factorizing terms in the variable $u$ in the conditions $E_i(u, \bar{u}) = 0$, these reduce to the following conditions on $\bar{u}$, referred to as the Bethe equations,

$$ \prod_{\substack{j=1 \\ i \neq j}}^{L} \frac{q(u_i^2 - q u_j^2)(\alpha \beta \gamma \delta u_i^2 u_j^2 - 1)}{(u_i^2 - q u_j^2)(\alpha \beta \gamma \delta u_i^2 u_j^2 - q^2)} = \frac{(q^2 - \alpha \beta \gamma \delta u_i^4)(\alpha \beta \gamma \delta u_i^2 q^{K+L+1} - 1)}{q^{-2L-2}(\alpha \beta \gamma \delta u_i^4 - 1)(q^{K+L+2} - u_i^2)} \times \\ \frac{(q+1)(q^2 u_i^4 - 1)(q^{K+1} - \alpha \beta u_i^2 q^L)}{\alpha \beta (q u_i^4 - 1)(q^L - \gamma \delta u_i^2 q^K)} \frac{a(u_i)}{a(\tau u_i^{-1})}. \quad (55) $$

To keep the notation simple, $\bar{u} = \{u_1, u_2, \dots, u_L\}$ will refer from now on to Bethe roots, i.e. to solutions of the set of equations (55).
## 4.2 Diagonalization of the truncated correlation matrix

Since the Heun operator T commutes with the truncated correlation matrix and is non-degenerate, its eigenvectors $|\bar{u}\rangle$ also diagonalize the matrix C,

$$C|\bar{u}\rangle = c(\bar{u})|\bar{u}\rangle, \quad c(\bar{u}) \in \mathbb{R}. \qquad (56)$$

To obtain an explicit expression for the eigenvalues $c(\bar{u})$ in terms of the parameters $\bar{u}$, we observe that the action of $\mathcal{B}(u, m)$ in the position basis is tridiagonal and given by

$$q^{m+1}\mathcal{B}(u, m)|n\rangle = V_{n,m}|n+1\rangle + (X_{n,m} + Y_{n,m}U)|n\rangle + Z_{n,m}|n-1\rangle, \qquad (57)$$

with $U = \frac{q}{u_1^2} + \alpha\beta\gamma\delta u^2$. The coefficients $V_{n,m}$, $X_{n,m}$, $Y_{n,m}$ and $Z_{n,m}$ can be computed directly from the definition (38) and the action of $\mathbf{A}$ and $\mathbf{A}^*$ in the position basis (see appendix A). In the case where $\beta = 0$, i.e. the dual q-Hahn special case of the q-Racah polynomials [8], these coefficients simplify greatly,

$$V_{n,m} = \frac{J_n}{q^{n+1}}, \quad X_{n,m} = \alpha\gamma q(q^{m+1} - q^n), \quad Y_{n,m} = 1 - q^{m+1-n}, \quad Z_{n,m} = 0. \qquad (58)$$

In particular, $\mathcal{B}(u, m)$ becomes a raising operator in the sense that $\langle n-1|\mathcal{B}(u, m)|n\rangle = 0$. This allows to compute the wavefunction $\langle n|\bar{u}\rangle$ of Bethe vectors:

$$\langle n|\bar{u}\rangle = q^{-L(L-1)/2} \left(\prod_{\ell=1}^{L-n} 1 - q^\ell\right) \left(\prod_{i=0}^{n-1} \frac{J_i}{q^{i+1}}\right) \sum_{r=0}^{L-n} \left(\frac{\alpha\gamma q(1-q^{n+1})}{q-1}\right)^{L-n-r} S_r(\bar{U}), \qquad (59)$$

where $\bar{U} = \{U_1, U_2, \dots, U_L\}$ with $U_i = q u_i^{-2}$. The terms $S_r(\bar{U})$ are symmetric polynomials of degree $r$ in the variables $U_i$ defined by

$$S_r(\bar{U}) = \sum_{i_1 < i_2 < \dots < i_r} U_{i_1} U_{i_2} \dots U_{i_r}, \quad S_0(\bar{U}) = 1. \qquad (60)$$

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 \in \{0, 1, \dots, L\}$, we find

$$c(\bar{u}) = \frac{\langle n | C | \bar{u} \rangle}{\langle n | \bar{u} \rangle} = \frac{q^{-L(L-1)/2}}{\langle n | \bar{u} \rangle} \sum_{r=0}^{L} b_{r,n} S_r(\bar{U}), \qquad (61)$$

where

$$b_{r,n} = \sum_{k=0}^{K} \sum_{n'=0}^{L-r} \phi_n(\omega_k) \phi_{n'}(\omega_k) \left( \prod_{\ell=1}^{L-n'} (1-q^\ell) \right) \left( \prod_{i=0}^{n'-1} \frac{J_i}{q^{i+1}} \right) \left( \frac{\alpha \gamma q (1-q^{n'+1})}{q-1} \right)^{L-n'-r} . \qquad (62)$$

This is valid for parameters $\bar{u}$ which are solutions of the Bethe equations (55). For $\beta = 0$, these equations reduce to

$$\prod_{j=1}^{L} \frac{\left(u_i^2 - \frac{u_j^2}{q}\right)}{\left(u_i^2 - q u_j^2\right)} = \frac{q^K (q - \alpha u_i^2) (q - \gamma u_i^2) (q - \gamma \delta u_i^2)}{(q^{K+L+2} - u_i^2) (\alpha \gamma u_i^2 - 1) (\gamma \delta u_i^2 q^K - q^L)}. \qquad (63)$$
## 5 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 $\beta = 0$, (51) can indeed be rewritten as

$$ U^2 Q(U) \Lambda(\bar{u}) = \frac{(q+1)(U-\alpha)(U-\gamma)(U-\gamma\delta)}{q^L} Q(qU) - p(U)Q(U) + q^L(q+1)(U-q^{K-L-1})(U-\alpha\gamma q)(U-\gamma\delta q^{K-L+1})Q(U/q) \quad (64) $$

where $p(U)$ is the following polynomial in the variable $U$,

$$ p(U) = -\frac{\alpha\gamma^2\delta(q+1)^2}{q^L} + \frac{\gamma(q+1)U}{q^{K+L}}(\alpha\gamma\delta q^{2K+L+2} + q^K(\alpha\delta + \alpha + \gamma\delta + \delta) + \alpha q^L) - 2q(\alpha\gamma + \alpha + \gamma\delta + \gamma)U^2 + 2(q+1)U^3, \quad (65) $$

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 $\bar{u} = \{u_1, u_2, \dots, u_L\}$:

$$ Q(U) = \prod_{i=1}^{L} (U - U_i) = \sum_{i=1}^{L} (-1)^{L-i} S_{L-i}(\bar{U}) U^i. \quad (66) $$

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 fix $\delta = 0$, $\gamma \in [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(\bar{U})$:

$$ 0 = \sigma_{n+1}S_{n+1} + (\rho_n + \Lambda(\bar{u}))S_n + \epsilon_{n-1}S_{n-1}, \quad (67) $$

where

$$ \sigma_n = (q+1)q^{-n} + (q+1)q^n - 2(q+1) \quad (68) $$

$$ \rho_n = (q+1)q^{n-L-K} \left( \alpha\gamma q^{K+L+1} + \frac{1}{q} \right) + (q+1)(\alpha+\gamma)q^{-n} - 2q(\alpha\gamma + \alpha + \gamma) \quad (69) $$

$$ \epsilon_n = -\alpha\gamma(q+1)(q^{-K} + q^{-L}) + \frac{\alpha\gamma(q+1)}{q^{K+L}}q^n + \alpha\gamma(q+1)q^{-n}. \quad (70) $$

In the thermodynamic limit $N \to \infty$, the parameter $\alpha = q^{-N-1}$, with $0 < q < 1$, goes to infinity and (67) becomes effectively a two term recurrence with solution

$$ S_n = S_L \prod_{i=n}^{L-1} \frac{\rho_{i+1} + \Lambda(\bar{u})}{\epsilon_i} + O(\alpha^{-1}). \quad (71) $$

The condition that $S_{-1}(\bar{U}) = 0$ then requires the eigenvalues $\Lambda(\bar{u})$ of T to take certain values

$$ \Lambda(\bar{u}) \in \{-\rho_n + O(\alpha^0) \mid n \in \{0, 1, \dots, L\} \}. \quad (72) $$

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.

<footnote n="1">¹The choice of parameters $\delta = 0$, $\beta = 0$ corresponds to the affine q-Krawtchouk limit of the q-Racah polynomials [8].</footnote>
<table><caption>Table 1: Eigenvalues of the Heun operator ($N = 49, L = 9, K = 24, q = 0.8, \alpha = q^{-N-1}$, $\beta = 0, \gamma = 0.5, \delta = 0$) obtained by three methods. The first column are zeros of $S_{-1}(Ū)$ seen as a polynomial of degree $L + 1$ in $\Lambda(Ū)$. 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].</caption><thead><tr><th>Solutions of<br>S<sub>-1</sub>(Ū) = 0</th><th>-ρ<sub>n</sub> for<br>n ∈ {0, 1, ..., L}</th><th>Numerical<br>diagonalization of T</th></tr></thead><tbody><tr><td>-778916</td><td>-778741</td><td>-778916</td></tr><tr><td>-592816</td><td>-592623</td><td>-592816</td></tr><tr><td>-444746</td><td>-444544</td><td>-444746</td></tr><tr><td>-327294</td><td>-327099</td><td>-327294</td></tr><tr><td>-234579</td><td>-234418</td><td>-234579</td></tr><tr><td>-161955</td><td>-161865</td><td>-161955</td></tr><tr><td>-105783</td><td>-105813</td><td>-105783</td></tr><tr><td>-63253.2</td><td>-63460.2</td><td>-63253.2</td></tr><tr><td>-32283.3</td><td>-32687.9</td><td>-32283.6</td></tr><tr><td>-11583.9</td><td>-11957.8</td><td>-11583.9</td></tr></tbody></table>

## 6 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 $\delta = 0$ and $N \to \infty$.

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.

## Acknowledgements

We thank Pascal Baseilhac and Rodrigo A. Pimenta for stimulating discussions. PAB holds an Alexander-Graham-Bell scholarship from the Natural Sciences and Engineering Research Council of Canada (NSERC). NC is supported by the international research project AAPT of the CNRS and the ANR Project AHA ANR- 18-CE40-0001. The research of LV is founded in part by a Discovery
Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada.

# A Appendix

## A.1 q-Racah polynomials

The q-Racah polynomials are defined by [8]

$$R_n(\omega_x) = {}_4\phi_3\left(\begin{array}{c} q^{-n} \alpha\beta q^{n+1} q^{-x} \gamma\delta q^{x+1} \\ \alpha q \beta\delta q \gamma q \end{array} \middle| q; q\right) \qquad (73)$$

with

$$\omega_x = q^{-x} + \gamma\delta q^{x+1}. \qquad (74)$$

The parameters are restricted by the truncation condition $R_{N+1}(x) = 0$. For instance, one can use $\alpha$ and fix

$$\alpha = q^{-N-1}. \qquad (75)$$

These polynomials also satisfy the following recurrence relation

$$(\omega_x - 1 - \gamma \delta q) R_n(\omega_x) = A_n R_{n+1}(\omega_x) - (A_n + C_n) R_n(\omega_x) + C_n R_{n-1}(\omega_x) \qquad (76)$$

where

$$A_n = \frac{(\alpha q^{n+1} - 1)(\gamma q^{n+1} - 1)(\alpha\beta q^{n+1} - 1)(\beta\delta q^{n+1} - 1)}{(1 - \alpha\beta q^{2n+1})(1 - \alpha\beta q^{2n+2})} \qquad (77)$$

$$C_n = \frac{(\beta q^n - 1)(\alpha q^n - \delta)(\alpha\beta q^n - \gamma)(q^{n+1} - q)}{(1 - \alpha\beta q^{2n})(1 - \alpha\beta q^{2n+1})} \qquad (78)$$

The normalisation weight is

$$W_k = \frac{(\beta^{-1}\gamma q, \delta q; q)_N(\gamma\delta q, \alpha q, \beta\delta q, \gamma q; q)_k(1 - \gamma\delta q^{2k+1})}{(\gamma\delta q^2, \beta^{-1}; q)_N(q, \alpha^{-1}\gamma\delta q, \beta^{-1}\gamma q, \delta q; q)_k(\alpha\beta q)^k(1 - \gamma\delta q)}. \qquad (79)$$

These polynomials also have a difference equation of the form (28), with coefficients given by

$$\bar{J}_k = \sqrt{\frac{(1 - \alpha q^{k+1})(1 - \beta \delta q^{k+1})(1 - \gamma q^{k+1})(1 - \gamma \delta q^{k+1})}{(1 - \gamma \delta q^{2k+1})(1 - \gamma \delta q^{2k+2})}} \times \\ \sqrt{\frac{(1 - q^{k+1})(1 - \alpha^{-1} \gamma \delta q^{k+1})(1 - \beta^{-1} \gamma q^{k+1})(1 - \delta q^{k+1})(\alpha \beta q)}{(1 - \gamma \delta q^{2k+2})(1 - \gamma \delta q^{2k+3})}} \qquad (80)$$

$$\bar{\mu}_k = \frac{(1 - \alpha q^{k+1})(1 - \beta \delta q^{k+1})(1 - \gamma q^{k+1})(1 - \gamma \delta q^{k+1})}{(1 - \gamma \delta q^{2k+1})(1 - \gamma \delta q^{2k+2})} + \frac{q(1 - q^k)(1 - \delta q^k)(\beta - \gamma q^k)(\alpha - \gamma \delta q^k)}{(1 - \gamma \delta q^{2k})(1 - \gamma \delta q^{2k+1})} - 1 - \alpha \beta q \qquad (81)$$
## A.2 Functions in the algebraic Bethe ansatz

The functions in the definition of the dynamical operators are:

$$f_1(u, m) = \frac{2q^{m+1}(q + \alpha\beta\gamma\delta u^4)}{u^2(\alpha\beta q^{2L+2m+1} - q^{2L})} - \frac{u^2\eta(q+1)q^{-2L+2}}{(q^2-1)^2(q^2 - \alpha\beta\gamma\delta u^4)} + \frac{(q+1)}{(q^2-1)^2(q^{2L} - \alpha\beta q^{2L+2m+1})} \left( \frac{\eta^*(\gamma\delta u^4 q^{-2L+1} - q^{2m+4})}{(q^2 - \alpha\beta\gamma\delta u^4)} - 2\xi q^{m+2} \right), \quad (82)$$

and

$$f_2(u, m) = \frac{(\alpha\beta q^{2L+2m+3} + q^{2L})(q + \alpha\beta\gamma\delta u^4)}{u^2 q^{m+2L+1}} + \frac{(\eta^* q^{2L+m+1} + \alpha\beta\xi q^{2L+2m+3} + \xi q^{2L})}{q^{m+2L}(q-1)^2(q+1)}. \quad (83)$$

The functions in the relation between the dynamical operators are:

$$f(u, v) = \frac{(u^2 - qv^2)(\alpha\beta\gamma\delta u^2 v^2 - q^2)}{q(u^2 - v^2)(\alpha\beta\gamma\delta u^2 v^2 - q)}, \quad (84)$$

$$g(u, v, m) = \frac{(q-1)(q^2 - \alpha\beta\gamma\delta v^4)(\alpha\beta v^2 q^{2L+2m+3} - u^2 q^{2L})}{q(u^2 - v^2)(\alpha\beta q^{2L+2m+3} - q^{2L})(q - \alpha\beta\gamma\delta v^4)}, \quad (85)$$

and

$$w(u, v, m) = \frac{\alpha\beta(q-1)(\alpha\beta\gamma\delta v^4 - 1)(\gamma\delta u^2 v^2 q^{2L} - q^{2(L+m+2)})}{(q^{2L} - \alpha\beta q^{2L+2m+3})(q - \alpha\beta\gamma\delta v^4)(q - \alpha\beta\gamma\delta u^2 v^2)}. \quad (86)$$

The coefficients giving the action of $\mathcal{B}(u, m)$ on vectors in the position basis are:

$$V_{n,m} = J_n (q^{-n-1} - \alpha\beta q^{m+2} - \alpha\beta q^{m+1} + \alpha^2\beta^2 q^{2m+n+4}), \quad (87)$$

$$X_{n,m} = -\frac{\mu_n \lambda_n (\alpha\beta q^{2m+3} + 1)}{(q+1)} - \frac{\alpha\beta(q+1)}{q^{-m-1}} \mu_n + \frac{(\eta^* q^{2L+m+1} + \alpha\beta\xi q^{2L+2m+3} + \xi q^{2L})}{(q-1)^2(q+1)q^{2L-1}}, \quad (88)$$

$$Y_{n,m} = -q^{m+1}\lambda_n + \frac{(\alpha\beta q^{2L+2m+3} + q^{2L})}{q^{2L}}, \quad (89)$$

and

$$Z_{n,m} = J_{n-1} \alpha \beta (q^{2m+3-n} + q^n - q^{m+2} - q^{m+1}). \quad (90)$$

## References

[1] L. Amico, R. Fazio, A. Osterloh and V. Vedral, *Entanglement in many-body systems*, Reviews of modern physics **80**(2), 517 (2008), doi:https://doi.org/10.48550/arXiv.quant-ph/0703044.

[2] M. A. Nielsen and I. L. Chuang, *Quantum Computation and Quantum Information: 10th Anniversary Edition*, Cambridge University Press, doi:10.1017/CBO9780511976667 (2010).
[3] G. Vidal, J. I. Latorre, E. Rico and A. Kitaev, *Entanglement in quantum critical phenomena*, Physical review letters **90**(22), 227902 (2003), doi:10.1103/PhysRevLett.90.227902.

[4] J. Eisert, M. Cramer and M. B. Plenio, *Colloquium: Area laws for the entanglement entropy*, Reviews of modern physics **82**(1), 277 (2010), doi:10.1103/RevModPhys.82.277.

[5] I. Peschel, *Calculation of reduced density matrices from correlation functions*, Journal of Physics A: Mathematical and General **36**(14), L205–L208 (2003), doi:10.1088/0305-4470/36/14/101.

[6] B.-Q. Jin and V. E. Korepin, *Quantum spin chain, Toeplitz determinants and the Fisher–Hartwig conjecture*, Journal of statistical physics **116**(1), 79 (2004), doi:10.1023/B:JOSS.0000037230.37166.42.

[7] M. Fagotti and P. Calabrese, *Universal parity effects in the entanglement entropy of XX chains with open boundary conditions*, Journal of Statistical Mechanics: Theory and Experiment **2011**(01), P01017 (2011), doi:10.1088/1742-5468/2011/01/P01017.

[8] R. Koekoek, P. A. Lesky and R. F. Swarttouw, *Hypergeometric orthogonal polynomials and their q-analogues*, Springer Science & Business Media, doi:https://doi.org/10.1007/978-3-642-05014-5 (2010).

[9] D. Gioev and I. Klich, *Entanglement entropy of fermions in any dimension and the Widom conjecture*, Physical review letters **96**(10), 100503 (2006), doi:https://doi.org/10.1103/PhysRevLett.96.100503.

[10] V. Eisler and I. Peschel, *Free-fermion entanglement and spheroidal functions*, Journal of Statistical Mechanics: Theory and Experiment **2013**(04), P04028 (2013), doi:10.1088/1742-5468/2013/04/P04028.

[11] N. Crampé, R. I. Nepomechie and L. Vinet, *Free-Fermion entanglement and orthogonal polynomials*, Journal of Statistical Mechanics: Theory and Experiment **2019**(9), 093101 (2019), doi:10.1088/1742-5468/ab3787.

[12] N. Crampé, R. I. Nepomechie and L. Vinet, *Entanglement in Fermionic Chains and Bispectrality*, chap. Chapter 13, pp. 77–96, doi:10.1142/9789811210679_0013.

[13] P.-A. Bernard, N. Crampé, R. I. Nepomechie, G. Parez, L. P. d'Andecy and L. Vinet, *Entanglement of inhomogeneous free fermions on hyperplane lattices*, Nuclear Physics B **984**, 115975 (2022), doi:https://doi.org/10.1016/j.nuclphysb.2022.115975.

[14] P. Baseilhac and R. A. Pimenta, *Diagonalization of the Heun-Askey-Wilson operator, Leonard pairs and the algebraic Bethe ansatz*, Nuclear Physics B **949**, 114824 (2019), doi:https://doi.org/10.1016/j.nuclphysb.2019.114824.

[15] J. Cao, H.-Q. Lin, K.-j. Shi and Y. Wang, *Exact solution of XXZ spin chain with unparallel boundary fields*, Nuclear Physics B **663**(3), 487 (2003), doi:https://doi.org/10.1016/S0550-3213(03)00372-9.

[16] R. I. Nepomechie, *Bethe ansatz solution of the open XXZ chain with nondiagonal boundary terms*, Journal of physics A: mathematical and general **37**(2), 433 (2003), doi:10.1088/0305-4470/37/2/012.
[17] N. Crampé, E. Ragoucy and D. Simon, Eigenvectors of open XXZ and ASEP models for a class of non-diagonal boundary conditions, Journal of Statistical Mechanics: Theory and Experiment 2010(11), P11038 (2010), doi:10.1088/1742-5468/2010/11/P11038.

[18] P.-A. Bernard, G. Carcone, N. Crampe and L. Vinet, Bethe ansatz diagonalization of the Heun-Racah operator, arXiv preprint arXiv:2209.09213 (2022), doi:https://doi.org/10.48550/arXiv.2209.09213.

[19] P.-A. Bernard, N. Crampé, D. Shaaban Kabakibo and L. Vinet, Heun operator of Lie type and the modified algebraic Bethe ansatz, Journal of Mathematical Physics 62(8), 083501 (2021), doi:https://doi.org/10.1063/5.0041097.

[20] J. Avan, S. Belliard, N. Grosjean and R. Pimenta, Modified algebraic Bethe ansatz for XXZ chain on the segment-III-Proof, Nuclear Physics B 899, 229 (2015), doi:https://doi.org/10.1016/j.nuclphysb.2015.08.006.

[21] S. Belliard, N. Crampé et al., Heisenberg XXX model with general boundaries: eigenvectors from algebraic bethe ansatz, SIGMA. Symmetry, Integrability and Geometry: Methods and Applications 9, 072 (2013), doi:https://doi.org/10.3842/SIGMA.2013.072.

[22] S. Belliard, N. A. Slavnov, B. Vallet et al., Modified algebraic Bethe ansatz: twisted XXX case, SIGMA. Symmetry, Integrability and Geometry: Methods and Applications 14, 054 (2018), doi:https://doi.org/10.3842/SIGMA.2018.054.

[23] N. Crampe et al., Algebraic Bethe ansatz for the XXZ Gaudin models with generic boundary, SIGMA. Symmetry, Integrability and Geometry: Methods and Applications 13, 094 (2017), doi:https://doi.org/10.3842/SIGMA.2017.094.

[24] N. Crampé, Algebraic Bethe ansatz for the totally asymmetric simple exclusion process with boundaries, Journal of Physics A: Mathematical and Theoretical 48(8), 08FT01 (2015), doi:10.1088/1751-8113/48/8/08FT01.

[25] K. Totsuka, Magnetization processes in bond-alternating quantum spin chains, Physics Letters A 228(1-2), 103 (1997), doi:https://doi.org/10.1016/S0375-9601(97)00087-X.

[26] J. A. Carrasco, F. Finkel, A. Gonzalez-Lopez and P. Tempesta, A duality principle for the multi-block entanglement entropy of free fermion systems, Scientific Reports 7(1), 1 (2017), doi:10.1038/s41598-017-09550-1.

[27] G. Bergeron, N. Crampé, S. Tsujimoto, L. Vinet and A. Zhedanov, The Heun-Racah and Heun-Bannai-Ito algebras, Journal of Mathematical Physics 61(8), 081701 (2020), doi:https://doi.org/10.1063/5.0008372.

[28] N. Crampé, L. Frappat, J. Gaboriaud, L. P. d'Andecy, E. Ragoucy and L. Vinet, The Askey-Wilson algebra and its avatars, Journal of Physics A: Mathematical and Theoretical 54(6), 063001 (2021), doi:10.1088/1751-8121/abd783.

[29] A. S. Zhedanov, "hidden symmetry" of Askey-Wilson polynomials, Theoretical and Mathematical Physics 89(2), 1146 (1991), doi:https://doi.org/10.1007/BF01015906.
[30] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett et al., *SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python*, Nature Methods **17**, 261 (2020), doi:10.1038/s41592-019-0686-2.