# Exactly solvable inhomogeneous XY spin chain

Pierre-Antoine Bernard, Nicolas Crampé, Quentin Labriet, Lucia Morey, Luc Vinet

► To cite this version:

Pierre-Antoine Bernard, Nicolas Crampé, Quentin Labriet, Lucia Morey, Luc Vinet. Exactly solvable inhomogeneous XY spin chain. Journal of Statistical Mechanics: Theory and Experiment, 2025, 2025 (10), pp.103101. 10.1088/1742-5468/ae0429. hal-05349319

HAL Id: hal-05349319

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

Submitted on 5 Nov 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
# Exactly solvable inhomogeneous XY spin chain

Pierre-Antoine Bernard¹*, Nicolas Crampé²†, Quentin Labriet¹‡, Lucia Morey¹§, Luc Vinet¹,³¶

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

² CNRS – Université de Montréal CRM - CNRS, France.

³ IVADO, Montréal (Québec), H2S 3H1, Canada.

## Abstract

Analytical expressions for the eigenvalues of certain inhomogeneous XY spin chains are computed. These models are rewritten in terms of free-fermion models using a well-known Jordan-Wigner transformation. Finding the spectrum of such models amounts to diagonalizing a matrix whose size is equal to the number of sites in the chain. This is achieved by recognizing and exploiting contiguity relations satisfied by specific orthogonal polynomials.

## 1 Introduction

Quantum spin chains describe a wide array of phenomena. They were introduced to study magnetization in metals [1] and are, among other things, used to investigate phase transition and thermalization (see [2, 3] for general references). These endeavours are facilitated by the existence of integrable models of spin chains, which allow the exploration of these questions by exact computations. For instance, Onsager provided a rigorous demonstration of the existence of an order-disorder phase transition in the two-dimensional Ising model by exactly computing the dominant contribution to its partition function [4]. This has been generalized to more general statistical models in [5] and paved the way to the quantum inverse scattering approach, initiated by the Leningrad school [6].

Recently, significant attention has been devoted to analyzing the impact of introducing position-dependent couplings and magnetic fields in spin chain models. These inhomogeneities, which can be thought of as representing the effects of defects in materials, are known to significantly influence the properties of spin chains. They can lead, for example, to violations of the area law for entanglement entropy [7, 8], induce localization [9], affect the scaling of magnetic conductivity [10], and are essential for achieving perfect state transfer [11] and fractional revival [12]. An overview of some recent results is provided in [13].

As a means to investigate such phenomena, exactly solvable models of inhomogeneous XX spin chains based on discrete families of Askey-Wilson polynomials were introduced in [14] (see also [15, 16]). This approach relies on the three-term recurrence relation satisfied by these polynomials to diagonalize the single-particle Hamiltonian. The stationary state wavefunctions of the latter then serve to define appropriate fermionic operators, allowing the full Hamiltonian to be diagonalized, just as in the homogeneous case [17]. Closed-form expressions for the spectra and stationary states can then be obtained. Similar constructions using multivariate orthogonal polynomials have also been used to study free-fermion models, not only on chains, but on higher dimensional lattices [18] as well. Let also mention that these models are closely related to field theory and experiments using the notion of entanglement Hamiltonian [19–21].

<footnote>*E-mail: pierre-antoine.bernard@umontreal.ca</footnote>

<footnote>†E-mail: crampe1977@gmail.com</footnote>

<footnote>‡E-mail: quentin.labriet@umontreal.ca</footnote>

<footnote>§E-mail: lucia.morey@umontreal.ca</footnote>

<footnote>¶E-mail: luc.vinet@umontreal.ca</footnote>
The aim of this paper is to extend the approach developed in [14] to construct exactly solvable models of XY spin chains. This involves diagonalizing a block tridiagonal matrix whose entries encode the inhomogeneous couplings and magnetic field characterizing the model. We show that, with an appropriate choice of these parameters, the block tridiagonal matrix can be diagonalized utilizing the so-called contiguity relations of Askey-Wilson polynomials. These were recently classified in [22], where they were interpreted in terms of spectral transforms.

Section 2 contains the definition of the inhomogeneous XY model as well as its connection to the free-fermion model under the Jordan-Wigner transformation. We recall that the spectrum of this model is obtained by diagonalizing a $d \times d$ matrix, where $d$ is the number of sites. This is explored in Section 3 where it is shown that this diagonalization is closely related to the $B_2$-contiguity relations which are also defined in this section. Section 4 is devoted to the examples associated with the $q$-Racah polynomials. Two examples corresponding to two different contiguity relations of the $q$-Racah polynomials are studied in detail. Their spectrum is explicitly given. Section 5 discusses different avenues for future research that this work suggests.

## 2 XY model and free fermions

This section aims to fix the notations while recalling the definition of the XY model and the usual methods to find its spectrum.

**Inhomogeneous XY model.** The model describes a chain of $N + 1$ spins 1/2 interacting between nearest neighbors and with an external magnetic field. The strengths of the interaction and of the magnetic field depend on the site. The quantum Hamiltonian of the inhomogeneous XY model with open boundaries is given by

$$ \mathcal{H} = \sum_{j=0}^{N-1} \left( (\alpha_j + \gamma_j) \sigma_j^x \sigma_{j+1}^x + (\alpha_j - \gamma_j) \sigma_j^y \sigma_{j+1}^y \right) - \sum_{j=0}^{N} \beta_j \sigma_j^z, \quad (2.1) $$

where $\sigma_j^x, \sigma_j^y$ and $\sigma_j^z$ are the Pauli matrices:

$$ \sigma^x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma^y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma^z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \quad (2.2) $$

acting on the $j^{th}$ $\mathbb{C}^2$-space of $(\mathbb{C}^2)^{\otimes(N+1)}$ and $\alpha_j, \beta_j$ and $\gamma_j$ are real parameters. For $\gamma_j = 0$, this model reduces to the XX model. The case when $\alpha_j, \beta_j$ and $\gamma_j$ are independent of $j$ is the homogeneous model solved exactly in [17] using the Jordan-Wigner transformation [23]. The fact that the model is one-dimensional is important for this transformation to apply. However, let us remark that there exist some attempts to generalize to higher dimension [24, 25] or to Y-junction [26], but the addition of supplementary sites is necessary in these cases.

**Jordan-Wigner transformation.** This transformation is based on the introduction of the following operators

$$ c_j = \sigma_0^z \dots \sigma_{j-1}^z \sigma_j^-, \quad c_j^\dagger = \sigma_0^z \dots \sigma_{j-1}^z \sigma_j^+, \quad (2.3) $$

where

$$ \sigma^+ = \frac{1}{2}(\sigma^x + i\sigma^y) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}, \quad \sigma^- = \frac{1}{2}(\sigma^x + i\sigma^y) = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix}. \quad (2.4) $$

These operators are fermionic operators since they satisfy, for $0 \le i, j \le N$,

$$ \{c_i, c_j\} = 0, \quad \{c_i^\dagger, c_j^\dagger\} = 0, \quad \{c_i^\dagger, c_j\} = \delta_{i,j}, \quad (2.5) $$

where $\{\cdot, \cdot\}$ is the anticommutator and $\delta_{i,j}$ is the Kronecker symbol. These fermionic operators allow us to rewrite the Hamiltonian (2.1) as the following free-fermion model:

$$ \mathcal{H} = 2 \sum_{j=0}^{N-1} \left( \alpha_j (c_j c_{j+1}^\dagger + c_{j+1} c_j^\dagger) + \gamma_j (c_j c_{j+1} + c_{j+1}^\dagger c_j) \right) + \sum_{j=0}^{N} \beta_j (2c_j c_j^\dagger - 1). \quad (2.6) $$

It is called a free-fermion model because it is quadratic in the fermionic operators. This Hamiltonian $\mathcal{H}$ can be rewritten more compactly as:

$$ \mathcal{H} = c^\dagger H c, \quad (2.7) $$
where $\mathbf{c}^{\dagger} = (c_0, \dots, c_N, c_0^{\dagger}, \dots, c_N^{\dagger})$ and $H$ is the $(2N+2) \times (2N+2)$-matrix

$$ H = \begin{pmatrix} A & B \\ -B & -A \end{pmatrix}, \qquad (2.8) $$

with $A$ a symmetric $(N+1) \times (N+1)$-matrix and $B$ an anti-symmetric $(N+1) \times (N+1)$-matrix taking the following forms

$$ A = \begin{pmatrix} \beta_0 & \alpha_0 & & & \\ \alpha_0 & \beta_1 & \alpha_1 & & \\ & \ddots & \ddots & \ddots & \\ & & \alpha_{N-2} & \beta_{N-1} & \alpha_{N-1} \\ & & & \alpha_{N-1} & \beta_N \end{pmatrix}, \quad B = \begin{pmatrix} 0 & \gamma_0 & & & \\ -\gamma_0 & 0 & \gamma_1 & & \\ & \ddots & \ddots & \ddots & \\ & & -\gamma_{N-2} & 0 & \gamma_{N-1} \\ & & & -\gamma_{N-1} & 0 \end{pmatrix}. \qquad (2.9) $$

The vector $\mathbf{c}$ is a column vector with components $c_0^{\dagger}, \dots, c_N^{\dagger}, c_0, \dots, c_N$.

**Properties of H.** The matrix $H$ is real symmetric and has real eigenvalues. In addition, using linearity and invariance by permutation of the rows and columns of its characteristic polynomial

$$ \mathcal{P}(x) = \begin{vmatrix} x-A & -B \\ B & x+A \end{vmatrix}, \qquad (2.10) $$

one gets $\mathcal{P}(-x) = \mathcal{P}(x)$. Therefore its diagonalized form can be written as follows

$$ H^{\text{diag}} = T^t H T = \begin{pmatrix} \Lambda & 0 \\ 0 & -\Lambda \end{pmatrix}, \qquad (2.11) $$

with $T$ an orthogonal transition matrix (i.e. $T^t T = \mathbb{I}_{N+1}$) and $\Lambda$ a $(N+1) \times (N+1)$-diagonal matrix with diagonal entries $\Lambda_0, \Lambda_1, \dots, \Lambda_N$, which are positive real numbers.

The spectral problem associated to the eigenvalue $\Lambda_j$ can be written

$$ H \begin{pmatrix} \psi_j \\ \phi_j \end{pmatrix} = \Lambda_j \begin{pmatrix} \psi_j \\ \phi_j \end{pmatrix}, \qquad (2.12) $$

where $\psi_j$ and $\phi_j$ are two vectors of $(N+1)$ components. Using the special form of $H$, one can show that the eigenvector associated to the eigenvalue $(-\Lambda_j)$ is

$$ H \begin{pmatrix} \phi_j \\ \psi_j \end{pmatrix} = -\Lambda_j \begin{pmatrix} \phi_j \\ \psi_j \end{pmatrix}. \qquad (2.13) $$

The orthogonal transition matrix $T$ takes the following form

$$ T = \begin{pmatrix} \Psi & \Phi \\ \Phi & \Psi \end{pmatrix}, \qquad (2.14) $$

where the $j^{th}$ column of $\Psi$ (resp. $\Phi$) is $\psi_j$ (resp. $\phi_j$).

**Spectrum of H.** Using the previous transition matrix, one can rewrite the Hamiltonian $\mathcal{H}$ as follows:

$$ \mathcal{H} = \mathbf{c}^{\dagger} T T^{t} H T T^{t} \mathbf{c} = \mathbf{d}^{\dagger} H^{\text{diag}} \mathbf{d}, \qquad (2.15) $$

where $\mathbf{d}^{\dagger} = \mathbf{c}^{\dagger} T$. With $\mathbf{d}^{\dagger} = (d_0, \dots, d_N, d_0^{\dagger}, \dots, d_N^{\dagger})$, one obtains

$$ d_k = \sum_{j=0}^{N} \psi_{jk} c_j + \sum_{k=0}^{N} \phi_{jk} c_j^{\dagger}, \quad d_k^{\dagger} = \sum_{j=0}^{N} \psi_{jk} c_j^{\dagger} + \sum_{k=0}^{N} \phi_{jk} c_j. \qquad (2.16) $$

Using $TT^t = \mathbb{I}_{N+1}$ and the explicit form (2.14) of $T$, one can show that the previous transformation on fermionic operators induced by $T$ is canonical i.e. that $d_j, d_j^{\dagger}$ are fermionic operators, for $0 \le i, j \le N$,

$$ \{d_i, d_j\} = 0, \quad \{d_i^{\dagger}, d_j^{\dagger}\} = 0, \quad \{d_i^{\dagger}, d_j\} = \delta_{i,j}. \qquad (2.17) $$
Finally, the Hamiltonian can be rewritten as

$$ \mathcal{H} = \sum_{j=0}^{N} \Lambda_j (2d_j d_j^\dagger - 1). \qquad (2.18) $$

Denoting by $|\emptyset\rangle$ the vector such that, for $0 \le j \le N$,

$$ d_j |\emptyset\rangle = 0, \qquad (2.19) $$

the eigenvectors of $\mathcal{H}$ are, for $0 \le \ell \le N$ and $0 \le \epsilon_1 < \dots < \epsilon_\ell \le N$,

$$ V_{\epsilon_1, \dots, \epsilon_\ell} = d_{\epsilon_1}^\dagger \dots d_{\epsilon_\ell}^\dagger |\emptyset\rangle, \qquad (2.20) $$

with eigenvalues

$$ E_{\epsilon_1, \dots, \epsilon_\ell} = 2 \sum_{j=1}^{\ell} \Lambda_{\epsilon_j} - \sum_{j=0}^{N} \Lambda_j. \qquad (2.21) $$

The goal of this paper is to provide parameters $\alpha_j$, $\beta_j$ and $\gamma_j$ such that we can provide an explicit expression for the eigenvalues $\Lambda_j$ and for the transition matrix $T$.

## 3 Exact solution and contiguity relations

In this section, we show that the entries of the eigenvectors satisfy nice properties recognized as $B_2$-contiguity relations.

**Eigenvectors of H.** Building upon the framework established in [17], the spectral problem for $H$ (2.12) can be reformulated:

$$ A\psi_j + B\phi_j = \Lambda_j\psi_j, \qquad (3.1a) $$

$$ -B\psi_j - A\phi_j = \Lambda_j\phi_j, \qquad (3.1b) $$

which are equivalently given by

$$ (A - B)(\psi_j - \phi_j) = \Lambda_j(\psi_j + \phi_j), \qquad (3.2a) $$

$$ (A + B)(\psi_j + \phi_j) = \Lambda_j(\psi_j - \phi_j). \qquad (3.2b) $$

Introducing the following notations for the $k^{th}$ components of the vectors $\psi_j - \phi_j$ and $\psi_j + \phi_j$,

$$ (\psi_j - \phi_j)_k = P_k(j), \quad (\psi_j + \phi_j)_k = Q_k(j), \qquad (3.3) $$

equations (3.2) become, for $0 \le k \le N$,

$$ \beta_k P_k(j) + (\alpha_k - \gamma_k) P_{k+1}(j) + (\alpha_{k-1} + \gamma_{k-1}) P_{k-1}(j) = \Lambda_j Q_k(j), \qquad (3.4a) $$

$$ \beta_k Q_k(j) + (\alpha_k + \gamma_k) Q_{k+1}(j) + (\alpha_{k-1} - \gamma_{k-1}) Q_{k-1}(j) = \Lambda_j P_k(j), \qquad (3.4b) $$

with the conventions $\alpha_{-1} = \gamma_{-1} = \alpha_N = \gamma_N = 0$. An important observation at this point is that these last relations bear a similarity to the $B_2$-contiguity relations studied in the context of orthogonal polynomials.

**$B_2$-contiguity relations.** We denote by $\{R_i(x; \rho)\}_{0 \le i \le N}$ a discrete family of polynomials of the $(q)$-Askey scheme where $\rho$ is the list containing all the parameters of these polynomials (in particular, it includes the integer $N$).

The $B_2$-contiguity relations take the following form

$$ \lambda_{x;\rho}^{+} R_i(x; \rho) = \Phi_i^{+1,+} R_{i+1}(\bar{x}; \bar{\rho}) + \Phi_i^{0,+} R_i(\bar{x}; \bar{\rho}) + \Phi_i^{-1,+} R_{i-1}(\bar{x}; \bar{\rho}), \qquad (3.5a) $$

$$ \lambda_{x;\rho}^{-} R_i(\bar{x}; \bar{\rho}) = \Phi_i^{+1,-} R_{i+1}(x; \rho) + \Phi_i^{0,-} R_i(x; \rho) + \Phi_i^{-1,-} R_{i-1}(x; \rho), \qquad (3.5b) $$

where $\bar{x}$ is a linear transformation of $x$ and $\bar{\rho}$ is the list of modified parameters in $\rho$. The functions $\lambda_{x;\rho}^{\pm}$, $\Phi_i^{+1,\pm}$, $\Phi_i^{0,\pm}$ and $\Phi_i^{-1,\pm}$ that appear in the previous relations will be given explicitly in the following for the $q$-Racah polynomials.
Let normalize these polynomials as follows, for $0 \le i \le N$,

$$P_i(x) = \sqrt{\lambda_{x;\rho}^{+} \Phi_0^{0,-} \prod_{k=0}^{i-1} \frac{\Phi_k^{0,+} \Phi_k^{+1,-}}{\Phi_k^{0,-} \Phi_{k+1}^{-1,+}}} R_i(x; \rho), \quad Q_i(x) = \sqrt{\lambda_{x;\rho}^{-} \Phi_0^{0,+} \prod_{k=0}^{i-1} \frac{\Phi_k^{0,-} \Phi_k^{+1,+}}{\Phi_k^{0,+} \Phi_{k+1}^{-1,-}}} R_i(\bar{x}; \bar{\rho}). \quad (3.6)$$

The previous functions $P_i(x)$ and $Q_i(x)$ satisfy relations (3.4) if $\bar{N} = N$ and if the following constraints are satisfied, for $0 \le i \le N$,

$$\frac{\Phi_i^{0,-} \Phi_i^{+1,+} \Phi_{i+1}^{-1,+} \Phi_{i+1}^{0,-}}{\Phi_i^{0,+} \Phi_i^{+1,-} \Phi_{i+1}^{-1,-} \Phi_{i+1}^{0,+}} = 1. \quad (3.7)$$

The normalizations for $R_i(x; \rho)$ and $R_i(x; \bar{\rho})$ have been chosen such that the two different eigenvalues in (3.5) become equal. The coefficients in (3.4) read

$$\Lambda_j = \sqrt{\lambda_{j;\rho}^{+} \lambda_{j;\rho}^{-}}, \quad \beta_j = \sqrt{\Phi_j^{0,+} \Phi_j^{0,-}}, \quad 2\alpha_j = \sqrt{\Phi_{j+1}^{-1,+} \Phi_j^{+1,-}} + \sqrt{\Phi_{j+1}^{-1,-} \Phi_i^{+1,+}}, \quad (3.8) \\ 2\gamma_j = \sqrt{\Phi_{j+1}^{-1,-} \Phi_j^{+1,+}} - \sqrt{\Phi_{j+1}^{-1,+} \Phi_i^{+1,-}}.$$

In [22], the $B_2$-contiguity relations for the finite families of orthogonal polynomials in the Askey scheme have been classified. We recall and use these results in Section 4 to obtain exactly solvable XY model.

# 4 q-Racah polynomials and XY model

The q-Racah polynomials are the most general orthogonal polynomials of the finite families of the q-Askey scheme. They are defined by, for $i = 0, 1, \dots, N$,

$$R_i(x; \rho) = {}_4\phi_3 \left( \begin{array}{c} q^{-i}, \ abq^{i+1}, \ q^{-x}, \ cq^{x-N} \\ aq, \ bcq, \ q^{-N} \end{array} \middle| q; q \right), \quad (4.1)$$

where $\rho = a, b, c, N, q$ and $N$ is a non-negative integer. The standard definition of the generalized q-hypergeometric function [27, 28] is used:

$$ {}_4\phi_3 \left( \begin{array}{c} q^{-i}, \ a_1, \ a_2, \ a_3 \\ b_1, \ b_2, \ b_3 \end{array} \middle| q; z \right) = \sum_{k=0}^{i} \frac{(q^{-i}, a_1, a_2, a_3; q)_k}{(q, b_1, b_2, b_3; q)_k} z^k, \quad (4.2) $$

where

$$(b_1, \dots, b_p; q)_k = (b_1; q)_k \dots (b_p; q)_k, \quad (b_i; q)_k = \prod_{\ell=0}^{k-1} (1 - b_i q^\ell). \quad (4.3)$$

The q-Racah polynomials $R_i(x; \rho)$ are polynomials of degree $i$ with respect to the variable $\lambda_{x,\rho} = -(1-q^{-x})(1-cq^{x-N})$.

The $B_2$-contiguity relations for the q-Racah polynomials has been classified in [22] (a complete list can be found in appendix A of this paper). Focusing on those relations such that $\bar{N} = N$, we see that they all satisfy the constraints (3.7). In the following, we provide the various exactly solvable inhomogeneous XY models obtained by this procedure.

## 4.1 Model based on (qRI/III) and (qRIII/I) contiguity relations

Let relation (3.5a) be the contiguity relation (qRI/III) given in Appendix A of [22]. The transformations of the parameters are given by

$$\bar{x} = x + 1, \quad \bar{a} = a/q, \quad \bar{b} = bq, \quad \bar{c} = c/q^2, \quad \bar{N} = N, \quad (4.4)$$

and the coefficients are:

$$\Phi_i^{+1,+} = - \frac{q^i(1-q^{i-N})(1-abq^{i+1})}{(1-abq^{2i+1})(1-abq^{2i+2})}, \quad \Phi_i^{-1,+} = - \frac{q^{i-N-1}(1-q^i)(1-abq^{N+i+1})}{(1-abq^{2i})(1-abq^{2i+1})}, \\ \Phi_i^{0+} = -\Phi_i^{+1,+} - \Phi_i^{-1,+}, \quad \lambda_{x,\rho}^{+} = \frac{(1-q^{-x-1})(1-cq^{x-N-1})}{(1-bc)(1-a)}.$$
The associated relation (3.5b) is the contiguity relation (qRIII/I) given in Appendix A of [22] for which the parameters are shifted. The coefficients read

$$ \lambda_{x, \rho}^{-} = (1 - cq^{x})(1 - q^{N-x}), $$

$$ \Phi_{i}^{+1,-} = \frac{(1 - q^{N-i})(1 - aq^{i})(1 - aq^{i+1})(1 - abq^{i+1})(1 - bcq^{i})(1 - bcq^{i+1})}{(1 - a)(1 - bc)(1 - abq^{2i+1})(1 - abq^{2i+2})}, $$

$$ \Phi_{i}^{0,-} = \frac{(1 - aq^{i})(1 - bq^{i+1})(1 - bcq^{i})(aq^{i+1} - c)}{q(1 - a)(1 - bc)(1 - abq^{2i+1})} \left( \frac{q(1 - q^{N-i})(1 - abq^{i+1})}{1 - abq^{2i+2}} + \frac{(1 - q^{-i})(1 - abq^{N+i+1})}{1 - abq^{2i}} \right), $$

$$ \Phi_{i}^{-1,-} = \frac{(1 - q^{-i})(aq^{i} - c)(aq^{i+1} - c)(1 - bq^{i})(1 - bq^{i+1})(1 - abq^{N+i+1})}{q(1 - a)(1 - bc)(1 - abq^{2i})(1 - abq^{2i+1})}. $$

By direct computations, we show that the constraints (3.7) are satisfied with these previous coefficients. Therefore, the XY model with the coefficients $\alpha_j$, $\beta_j$ and $\gamma_j$ given by (3.8) is exactly solvable with the eigenvalues

$$ \Lambda_j = \sqrt{\frac{(1 - cq^j)(1 - q^{N-j})(1 - q^{-j-1})(1 - cq^{j-N-1})}{(1 - bc)(1 - a)}}. \quad (4.5) $$

## 4.2 Model based on (qRII/IV) and (qRIV/II) contiguity relations

Let relation (3.5a) be the contiguity relation (qRII/IV) given in Appendix A of [22]. The transformations of the parameters are given by

$$ \bar{x} = x, \quad \bar{a} = a/q, \quad \bar{b} = bq, \quad \bar{c} = c, \quad \bar{N} = N, \quad (4.6) $$

and the coefficients are:

$$ \Phi_{i}^{+1,+} = - \frac{(q^{N} - q^{i})(1 - abq^{i+1})(1 - bcq^{i+1})(1 - bcq^{i+2})}{(1 - bcq)(1 - abq^{2i+1})(1 - abq^{2i+2})}, $$

$$ \Phi_{i}^{-1,+} = - \frac{bq(1 - q^{i})(c - aq^{i-1})(c - aq^{i})(1 - abq^{N+i+1})}{a(1 - bcq)(1 - abq^{2i})(1 - abq^{2i+1})}, $$

$$ \Phi_{i}^{0+} = \lambda_{0, \rho}^{+} - \Phi_{i}^{+1,+} - \Phi_{i}^{-1,+}, \quad \lambda_{x, \rho}^{+} = \frac{(1 - aq^{x})(c - aq^{N-x})}{a(1 - a)}. $$

The associated relation (3.5b) is the contiguity relation (qRIV/II) given in Appendix A of [22] for which the parameters are shifted. The coefficients read

$$ \Phi_{i}^{+1,-} = \frac{(1 - aq^{i})(1 - aq^{i+1})(1 - abq^{i+1})(q^{i} - q^{N})}{(1 - a)(1 - abq^{2i+1})(1 - abq^{2i+2})}, \quad \Phi_{i}^{-1,-} = - \frac{a(1 - q^{i})(1 - bq^{i})(1 - bq^{i+1})(1 - abq^{N+i+1})}{bq(1 - a)(1 - abq^{2i})(1 - abq^{2i+1})} $$

$$ \lambda_{x, \rho}^{-} = \frac{(1 - bcq^{x+1})(1 - bq^{N-x+1})}{bq(1 - bcq)}, \quad \Phi_{i}^{0-} = \lambda_{0, \rho}^{+} - \Phi_{i}^{+1,+} - \Phi_{i}^{-1,+}. $$

By direct computations, we show that the constraints (3.7) are satisfied with these previous coefficients. Therefore, the XY model with the coefficients $\alpha_j$, $\beta_j$ and $\gamma_j$ given by (3.8) is exactly solvable with the eigenvalues

$$ \Lambda_j = \sqrt{\frac{(1 - aq^j)(c - aq^{N-j})(1 - bcq^{j+1})(1 - bq^{N-j+1})}{abq(1 - a)(1 - bcq)}}. \quad (4.7) $$

# 5 Outlooks

This paper establishes a novel connection between exactly solvable quantum models and special functions, paving the way for several interesting research avenues.

The models introduced herein depend on the parameters $a, b, c$ and $q$. It would be highly valuable to explore different limits of these parameters to connect these new models with other orthogonal polynomials within the (q-)Askey scheme. Such explorations are expected to reveal diverse behaviors of the associated eigenvalues.

We focused solely on the eigenvalues to demonstrate the feasibility of obtaining exact results for inhomogeneous XY spin chains. However, since the eigenvectors are associated with orthogonal polynomials, it should be
possible to compute other physical quantities, such as density, correlation functions, phase diagram and entan-
glement entropy. For the last quantity, we wonder if the computation involving Heun operators is also possible
as for the XX model [14].

Another natural extension of this work involves generalizing the proposed construction. In reference [22], it
was demonstrated that the $B_2$-contiguity relations correspond to spectral transforms of polynomials, specifically
Christoffel and Geronimus transforms. Exploring more general Christoffel and Geronimus transforms could
potentially yield new exactly solvable XY models.

Finally, the potential role of matrix valued orthogonal polynomials warrants depeer investigation. Notably,
the system described in equation (3.4) can be reformulated as a three term recurrence relation with matrix
coefficients, a structure closely tied to the theory of matrix or vector valued orthogonal polynomials. The
connection between integrable systems and such polynomials has been investigated, primarily in the development
of non-Abelian generalizations of classical integrable structures, such as the Painlevé equations and the Toda
lattice, see for instance [29–31]. Consequently, a systematic exploration of this matrix polynomial structure may
yield new insights for the studies of XY models.

**Acknowledgements:** The authors thank warmly G. Parez for primary discussions about this work. P.-A. Bernard holds an Alexander-Graham-Bell scholarship from the Natural Sciences and Engineering Research Council (NSERC) of Canada. N. Crampé is partially supported by the international research project AAPT of the CNRS. L. Vinet is funded in part by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada. Q. Labriet and L. Morey enjoy postdoctoral fellowships provided by this grant.

References

[1] W. Heisenberg, “Zur theorie des ferromagnetismus,” *Z. Phys.* **49**, 619–636 (1928).

[2] S. Sachdev, “Quantum phase transitions,” *Physics world* **12**, 33 (1999).

[3] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, “Colloquium: Many-body localization, thermalization, and entanglement,” *Rev. Mod. Phys.* **91**, 021001 (2019).

[4] L. Onsager, “Crystal statistics. I. A two-dimensional model with an order-disorder transition,” *Phys. Rev.* **65**, 117 (1944).

[5] R. Baxter, “Partition function of the eight-vertex lattice model,” *Ann. Phys.* **70**, 193 (1972).

[6] L. Faddeev, N. Reshetikhin, and L. Takhtajan, “Quantization of lie groups and lie algebras,” *Leningrad Math. J.* **1**, 193 (1990).

[7] G. Ramírez, J. Rodríguez-Laguna, and G. Sierra, “Entanglement over the rainbow,” *J. Stat. Mech.* **2015**, P06002 (2015).

[8] G. Vitagliano, A. Riera, and J. I. Latorre, “Volume-law scaling for the entanglement entropy in spin-1/2 chains,” *New J. Phys.* **12**, 113049 (2010).

[9] P. W. Anderson, “Absence of diffusion in certain random lattices,” *Phys. Rev.* **109**, 1492 (1958).

[10] P.-A. Bernard, I. Bussière, R. Floreanini, and L. Vinet, “Currents in nonequilibrium steady states of open inhomogeneous xx spin chains,” *Phys. Rev. A* **111**, 022208 (2025).

[11] M. Christandl, N. Datta, A. Ekert, and A. J. Landahl, “Perfect state transfer in quantum spin networks,” *Phys. Rev. Lett.* **92**, 187902 (2004).

[12] M. Christandl, L. Vinet, and A. Zhedanov, “Analytic next-to-nearest-neighbor XX models with perfect state transfer and fractional revival,” *Phys. Rev. A* **96**, 032335 (2017).

[13] P.-A. Bernard, G. Parez, and L. Vinet, “Distinctive features of inhomogeneous spin chains,”
*arXiv:2411.09487*.

[14] N. Crampé, R. I. Nepomechie, and L. Vinet, “Free-fermion entanglement and orthogonal polynomials,” *J. Stat. Mech.* **2019**, 093101 (2019).

[15] N. Crampé, R. I. Nepomechie, and L. Vinet, “Entanglement in fermionic chains and bispectrality,” *Rev. Math. Phys.* **33**, 2140001 (2021).
[16] R. Sasaki, “Exactly solvable inhomogeneous fermion systems,” *Prog. Theor. Exp. Phys.* **2024**, 123A03 (2024).

[17] E. Lieb, T. Schultz, and D. Mattis, “Two soluble models of an antiferromagnetic chain,” *Ann. Phys.* **16**, 407–466 (1961).

[18] P.-A. Bernard, N. Crampé, R. I. Nepomechie, G. Parez, L. Poulain d'Andecy, and L. Vinet, “Entanglement of inhomogeneous free fermions on hyperplane lattices,” *Nucl. Phys. B* **984**, 1159–1175 (2022).

[19] M. Dalmonte, V. Eisler, M. Falconi, and B. Vermersch, “Entanglement hamiltonians: From field theory to lattice models and experiments,” *Annalen der Physik* **534**, 2200064 (2022).

[20] R. Bonsignori and V. Eisler, “Entanglement hamiltonian for inhomogeneous free fermions,” *J. Phys. A* **57**, 275001 (2024).

[21] P.-A. Bernard, R. Bonsignori, V. Eisler, G. Parez, and L. Vinet, “Entanglement hamiltonian and orthogonal polynomials,” *arXiv:2412.12021*.

[22] N. Crampé, L. Morey, L. Vinet, and M. Zaimi, “Contiguity relations for finite families of orthogonal polynomials in the Askey scheme,” *arXiv:2504.20802*.

[23] P. Jordan and E. Wigner, “Über das Paulische Äquivalenzverbot,” *Z. Physik* **47**, 631–651 (1928).

[24] F. Verstraete and J. Cirac, “Mapping local Hamiltonians of fermions to local Hamiltonians of spins,” *J. Stat. Mech.* **P09012** (2005).

[25] V. Galitski, “Fermionization transform for certain higher-dimensional quantum spin models,” *Phys. Rev. B* **82**, 060411 (2010).

[26] N. Crampé and A. Trombettoni, “Quantum spins on star graphs and the Kondo model,” *Nucl. Phys. B* **871**, 526–538 (2013).

[27] G. Gasper and M. Rahman, *Basic Hypergeometric Series*, vol. 96 of *Encyclopedia of Mathematics and Its Applications*. Cambridge University Press, 2nd ed., 2004.

[28] R. Koekoek, P. A. Lesky, and R. F. Swarttouw, *Hypergeometric orthogonal polynomials and their q-analogues*. With a foreword by Tom H. Koornwinder. Springer Monogr. Math. Berlin: Springer, 2010.

[29] A. Deaño, L. Morey, and P. Román, “Non-Abelian Toda-type equations and matrix valued orthogonal polynomials,” *Proc. Amer. Math. Soc.* **152**, 1613–1632 (2024).

[30] M. Cafasso and M. D. de la Iglesia, “Non-commutative Painlevé equations and Hermite-type matrix orthogonal polynomials,” *Comm. Math. Phys.* **326**, 559–583 (2014).

[31] G. Ariznabarreta and M. Mañas, “Matrix orthogonal Laurent polynomials on the unit circle and Toda type integrable systems,” *Adv. Math.* **264**, 396–463 (2014).