# Entanglement of free fermions on Hamming graphs

Pierre-Antoine Bernard, Nicolas Crampé, Luc Vinet

► To cite this version:

Pierre-Antoine Bernard, Nicolas Crampé, Luc Vinet. Entanglement of free fermions on Hamming graphs. Nuclear Physics B, 2023, 986, pp.116061. 10.1016/j.nuclphysb.2022.116061. hal-03926291

HAL Id: hal-03926291

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

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
# Entanglement of Free Fermions on Hamming Graphs

Pierre-Antoine Bernard¹, Nicolas Crampé², and 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,

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

March 30, 2021

## Abstract

Free fermions on Hamming graphs $H(d, q)$ are considered and the entanglement entropy for two types of subsystems is computed. For subsets of vertices that form Hamming subgraphs, an analytical expression is obtained. For subsets corresponding to a neighborhood, i.e. to a set of sites at a fixed distance from a reference vertex, a decomposition in irreducible submodules of the Terwilliger algebra of $H(d, q)$ also yields a closed formula for the entanglement entropy. Finally, for subsystems made out of multiple neighborhoods, it is shown how to construct a block-tridiagonal operator which commutes with the entanglement Hamiltonian. It is identified as a BC-Gaudin magnet Hamiltonian in a magnetic field and is diagonalized by the modified algebraic Bethe ansatz.

## 1 Introduction

Entanglement entropy describes the extent to which a subsystem is correlated with its complementary part in a given state of the system. It is a key quantity in quantum information and in the investigation of quantum many-body models [21, 25]. Owning to their simplicity, free fermions offer instructive testbeds for its study.

We here examine entanglement in a system of free fermions hopping on the vertices of Hamming graphs $H(d, q)$. This family of graphs has been extensively analyzed and includes the complete graphs $H(1, q)$ as well as the $d$-cubes $H(d, 2)$.

In algebraic combinatorics, Hamming graphs provide the prototype of a $P$- and $Q$- association scheme and play a crucial role in coding theory [1, 5]. The algebra spanned by the adjacency and dual adjacency matrices of $H(d, q)$ is
referred to as its Terwilliger algebra $\mathcal{T}$ [15, 22, 29, 30, 31]. Most objects that we shall use to determine the entanglement entropy arise from this structure.

Hamming graphs are known to admit both perfect state transfer [7] and fractional revival [3]. The entanglement entropy of coupled harmonic oscillators on Hamming networks was considered in [17, 18], where it was shown that the Schmidt numbers and thus the entropy could be extracted from the potential matrix of the system expressed in the right basis. Looking at free fermions, we shall use the same basis, which is tied to the irreducible representations of the Terwilliger algebra of the Hamming scheme.

We shall examine the entanglement entropy when the system is in its ground state. For free fermions, it is known that the spectrum of the chopped correlation matrix $C$ suffices to compute the entanglement entropy [24, 26]. We shall look for the eigenvalues of $C$ for two types of subsystems.

First, we consider subsets of vertices that form by themselves a Hamming graph. For instance, it may correspond to a 2-cube in a 3-cube, i.e. to one of the faces of a three dimensional hypercube. For such cases, a general formula for the eigenvalues of $C$, their degeneracy and the entanglement entropy can be derived.

Second, we examine neighborhoods. A neighborhood is the set of vertices at a given distance from a reference site. For this type of subsystem, the restricted or chopped correlation matrix can be expressed in terms of elements in $\mathcal{T}$, the Terwilliger algebra of $H(d, q)$. The decomposition in irreducible $\mathcal{T}$-submodules of the vector space on which $C$ is acting will reduce our system of free fermions on a graph to a direct sum of fermionic chains (or paths). For a subsystem composed of a single neighborhood, this decomposition diagonalizes $C$ and allows a direct computation of the entanglement entropy. For a subsystem composed of a large number of neighborhoods, it is however not sufficient. Nevertheless, we can use a method recently introduced in [8] to construct a block-tridiagonal operator $\mathcal{T}$, referred to as a generalized Heun operator, which commutes with the chopped correlation matrix. This approach is an application to the computation of entanglement entropy of the tools associated to time and band limiting problems [9, 10, 20, 28]. For the systems considered here, the generalized Heun operator corresponds to a BC-Gaudin magnet Hamiltonian in a magnetic field and it proves possible to diagonalize this operator using the modified algebraic Bethe ansatz [4, 11]. This done, the eigenvalues of $C$ and the entropy can be extracted.

The paper is organized as follows. In section 2 and 3, we introduce the Hamiltonian associated to free fermions on a Hamming graph. We diagonalize it, define its ground state and discuss entanglement entropy. In section 4, we obtain the entropy for subsystems corresponding to Hamming subgraphs. In section 5, we consider subsystems associated to the bundles of neighborhoods. We give an overview of the Terwilliger algebra of the Hamming scheme and of its relevance to the problem at hand. We present the decomposition in irreducible $\mathcal{T}$-submodules and explain the simplification this brings. $\mathcal{T}$ will also be seen to be the quotient of the Onsager algebra [14, 23, 27] by certain Davies relations [12]. We show how the eigenvalues of the chopped correlation matrix $C$ are
obtained when the subsystem is a single neighborhood. The generalized Heun operator T is introduced in the subsection 5.3. We indicate that it commutes with C and diagonalize it using the modified algebraic Bethe ansatz.

## 2 Free fermions on Hamming graphs

The Hamming graph $H(d, q)$ can be constructed in the following way. Take the set of vertices $F_q^d$ to consist of all the d-tuples $v = (v_1, v_2, \dots, v_d)$ for which the elements are integers in $\{0, 1, \dots, q-1\}$. $v$ and $v'$ are connected by an edge when there exists a unique position $i$ in the d-tuples such that $v_i \neq v_i'$. It follows that the distance between two vertices in the resulting graph is given by the Hamming distance $\partial$ between their associated tuple:

$$ \partial(v, v') = \#\{i \in \{1, \dots, d\} \mid v_i \neq v_i'\}. \qquad (1) $$

It is worth noting that the case $d=1$ corresponds to the complete graph $K_q$ and that the case $q=2$ leads to the d-cube. We shall consider systems of free fermions living on the vertices of a Hamming graph. We shall study in particular fermionic systems for which the hopping constant between the sites $v$ and $v'$, $\alpha_{\partial(v,v')} \in \mathbb{R}$, depends only on their Hamming distance $\partial(v, v')$. More precisely, the Hamiltonian is defined as

$$ \hat{\mathcal{H}} = \sum_{v, v' \in F_q^d} \alpha_{\partial(v, v')} c_v^\dagger c_{v'}, \qquad (2) $$

where $c_v^\dagger$ and $c_v$ are fermionic creation and annihilation operators associated to the site $v$ in the graph. They satisfy the following canonical relations:

$$ \{c_v, c_{v'}\} = 0, \quad \{c_v^\dagger, c_{v'}^\dagger\} = 0 \quad \text{and} \quad \{c_v, c_{v'}^\dagger\} = \delta_{vv'}, \quad \text{for } v, v' \in F_q^d. \qquad (3) $$

Since the diameter of the graph $H(d, q)$ is $d$, the model defined by (2) contains $d+1$ independent parameters. The constant $\alpha_0$ is related to the presence of an external magnetic field.

We can give an alternative expression for $\hat{\mathcal{H}}$. For $i \in \{0, 1, \dots, d\}$, one can define the $i$th adjacency matrix $A_i$ of $H(d, q)$ as the matrix whose entry $[A_i]_{vv'}$ is

$$ [A_i]_{vv'} = \begin{cases} 1 & \text{if } \partial(v, v') = i, \\ 0 & \text{otherwise.} \end{cases} \qquad (4) $$

Furthermore, the vertex associated to the tuple $v = (v_1, \dots, v_d)$ can be represented by the following vector in $(\mathbb{C}^q)^{\otimes d}$:

$$ |v\rangle = |v_1\rangle \otimes |v_2\rangle \otimes \dots \otimes |v_d\rangle, \qquad (5) $$
where $|v_i\rangle$ is the vector in $\mathbb{C}^q$ which has its $(v_i + 1)^{\text{th}}$ entry equals to 1 as its unique non-zero entry:

$$|v_i\rangle = \underbrace{(0, 0, \dots, 0, 1, 0, \dots, 0)^T}_{v_i \text{ times}}. \qquad (6)$$

In this basis, the first adjacency matrix is expressed as

$$A_1 \equiv A = \sum_{i=1}^{d} \underbrace{\mathbb{1}_{q \times q} \otimes \dots \otimes \mathbb{1}_{q \times q}}_{i-1 \text{ times}} \otimes (J_{q \times q} - \mathbb{1}_{q \times q}) \otimes \underbrace{\mathbb{1}_{q \times q} \otimes \dots \otimes \mathbb{1}_{q \times q}}_{d-i \text{ times}}, \quad (7)$$

where $\mathbb{1}_{q \times q}$ is the identity matrix and $J_{q \times q}$ is the matrix of ones. More generally, we have

$$A_i = \sum_{\substack{b \in F_2^d \\ \text{wt}(b)=i}} \underbrace{(J_{q \times q} - \mathbb{1}_{q \times q})^{b_1} \otimes (J_{q \times q} - \mathbb{1}_{q \times q})^{b_2} \otimes \dots \otimes (J_{q \times q} - \mathbb{1}_{q \times q})^{b_d}}_{d \text{ times}}, \quad (8)$$

where the weight of a binary string $b$ is $\text{wt}(b) = \#\{n : b_n \neq 0\} = \sum_{n=1}^d b_n$. Using the vectors of operators $\hat{c}^\dagger = \sum_{v \in F_q^d} c_v^\dagger \langle v |$ and $\hat{c} = \sum_{v \in F_q^d} c_v |v\rangle$, the Hamiltonian can be rewritten as

$$\hat{\mathcal{H}} = \hat{c}^{\dagger} \left[ \sum_{i=0}^{d} \alpha_i A_i \right] \hat{c}. \qquad (9)$$

## 2.1 Diagonalization and energies

We are now interested in obtaining the single-particle excitation energies of the system. Diagonalizing $\hat{\mathcal{H}}$ amounts to diagonalizing $\sum_{i=0}^d \alpha_i A_i$, and in fact $A_1$. Indeed, since the Hamming graphs are distance-regular, it is known from the theory of association schemes that the matrix $A_i$ is a polynomial of degree $i$ of the first adjacency matrix $A$ [1]. Specifically,

$$\begin{aligned} A_i &= \binom{d}{i} (q-1)^i \sum_{j=0}^{i} \frac{(-i)_j \left( \frac{(q-1)d}{q} + \frac{A}{q} \right)_j}{(-d)_j j!} \binom{q}{q-1}^j \\ &= \binom{d}{i} (q-1)^i K_i \left( \frac{(q-1)d}{q} - \frac{A}{q}; \frac{q-1}{q}, d \right), \end{aligned} \qquad (10)$$

where $(a)_i = (a)(a+1)\dots(a+i-1)$ and $K_i$ refers to the Krawtchouk polynomial of degree $i$ [19]. Therefore, diagonalizing all the $d+1$ adjacency matrices is the same as diagonalizing $A$. This can be achieved by considering tensor products of $d$ eigenvectors of the matrix $J_{q \times q}$. This matrix has a unique eigenvector $|\theta_q\rangle$ of eigenvalue $q$:

$$|\theta_q\rangle = \frac{1}{\sqrt{q}} \sum_{i=0}^{q-1} |i\rangle, \qquad (11)$$
and its other eigenspaces have eigenvalue 0 and a degeneracy of $q-1$. We take $\{|\theta_i\rangle : i \in \{1, \dots, q\}\}$ to be an orthonormal basis of this space. Thus, we have $\langle \theta_i | \theta_j \rangle = \delta_{ij}$,

$$J_{q \times q} |\theta_q\rangle = q |\theta_q\rangle \quad \text{and} \quad J_{q \times q} |\theta_i\rangle = 0 \quad \text{if } i \neq q. \qquad (12)$$

Then, if $|\theta_{i_1} \theta_{i_2} \dots \theta_{i_d}\rangle = |\theta_{i_1}\rangle \otimes |\theta_{i_2}\rangle \otimes \dots \otimes |\theta_{i_d}\rangle$ with $i_1, \dots, i_{d-1}$ and $i_d \in \{1, \dots, q\}$, we see from (7) that

$$A |\theta_{i_1} \theta_{i_2} \dots \theta_{i_d}\rangle = (kq - d) |\theta_{i_1} \theta_{i_2} \dots \theta_{i_d}\rangle, \qquad (13)$$

where $k$ is the number of $|\theta_q\rangle$ in the tensor product. The spectrum of $A$ is thus $\omega_k = kq - d$ with $k \in \{0, 1, \dots, d\}$. As for the degeneracy of the $k^{th}$ eigenspace, it is easy to see that it is given by $D_k = \binom{d}{k}(q-1)^{d-k}$. Let us also note that since $|\theta_q\rangle \langle \theta_q| = J_{q \times q}/q$, the projection operator $E_k$ over the $(d-k)^{th}$ eigenspace of $A$ can be expressed as

$$E_k = \sum_{\substack{b \in F_2^d, \\ \text{wt}(b)=k}} (b_1 \mathbf{1} + (-1)^{b_1} \frac{J_{q \times q}}{q}) \otimes \dots \otimes (b_d \mathbf{1} + (-1)^{b_d} \frac{J_{q \times q}}{q}). \qquad (14)$$

From the spectrum of $A$ and (10), we can deduce that the eigenvalues of $\sum_{i=0}^d \alpha_i A_i$, i.e. the single-particle excitation energies $\Omega_k$ of the system, are

$$\Omega_k = \sum_{i=0}^{d} \alpha_i \binom{d}{i} (q-1)^i K_i(d-k; \frac{q-1}{q}, d). \qquad (15)$$

Indeed, if we take the set $\{|\omega_k, l\rangle\}$ ($l$ labels the degeneracy) to be an orthonormal basis of eigenvectors of $A$, we have

$$\hat{\mathcal{H}} = \sum_{k=0}^{d} \sum_{l=1}^{D_k} \Omega_k \bar{c}_{kl}^{\dagger} \bar{c}_{kl}, \qquad (16)$$

where $\bar{c}_{kl}^{\dagger} = \sum_v \langle v | \omega_k, l \rangle c_v^{\dagger}$ and $\bar{c}_{kl} = \sum_v \langle \omega_k, l | v \rangle c_v$. One can check that these new creation and annihilation operators respect the same canonical relations as the operators $c_v$ and $c_v^{\dagger}$. Formula (15) reduces to a simpler form in some useful cases. For instance, if we restrict ourselves to the nearest neighbour interactions $(\alpha_i = 0$ for $i \neq \{0, 1\})$, the energies are linear in $k$:

$$\Omega_k = \alpha_0 + \alpha_1(kq - d). \qquad (17)$$

Furthermore, in the case where the hopping constants decrease exponentially with the distance, we have $\alpha_i = e^{-ci}$ for $i > 0$, where $c \ge 0$. Thus, we find that (15) corresponds to the generating function of the Krawtchouk polynomials which grows exponentially with $k$:

$$\Omega_k = (1 - e^{-c})^{d-k} (1 + e^{-c}(q-1))^k + \alpha_0 - 1. \qquad (18)$$
## 2.2 Ground state

Let $|0\rangle\rangle$ be the vacuum state annihilated by all the operators $\bar{c}_{kl}$. The ground state $|\Psi_0\rangle\rangle$ is the state for which all the energy levels $\Omega_k < 0$ are occupied, which corresponds to filling up the Fermi sea. We denote $SE$ the set of all the integers $k \in \{0, 1, \dots, d\}$ associated to a negative single-particle excitation energy $\Omega_k$. For some fixed parameters $\alpha_i$, one can easily identify $SE$ by computing the values taken by (15). In the case where the hopping terms decrease exponentially with the distance, we know from the last section that it corresponds to a set of the form $\{0, 1, \dots, k_0\}$ for some integer $k_0$. In any case, we have:

$$|\Psi_0\rangle\rangle = \left[ \prod_{k \in SE} \prod_{l=1}^{D_k} \bar{c}_{kl}^{\dagger} \right] |0\rangle\rangle. \qquad (19)$$

As we will see in the next section, the key information we need is contained in the correlation matrix $\hat{C}$ whose components $\hat{C}_{vv'}$ are defined as

$$\hat{C}_{vv'} = \langle\langle\Psi_0|c_v^\dagger c_{v'}|\Psi_0\rangle\rangle, \quad \text{where } v, v' \in F_q^d. \qquad (20)$$

We can use the eigenbasis of $A$ to express $c_v^\dagger$ and $c_v$ in terms of $\bar{c}_{kl}^\dagger$ and $\bar{c}_{kl}$. Then, simple algebraic manipulations are sufficient to show that

$$\begin{aligned} \hat{C} &= \sum_{k \in SE} \sum_{l=1}^{D_k} |\omega_k, l\rangle \langle\omega_k, l| \\ &= \sum_{k \in SE} E_{d-k} \equiv \pi_{SE}, \end{aligned} \qquad (21)$$

where $E_{d-k}$ is the projection operator onto the $k^{th}$ eigenspace of $A$ given by (14) and $\pi_{SE}$ is the projection operator onto all the eigenspaces associated to an integer in $SE$.

## 3 Entanglement entropy

Let $SV \subset F_q^d$ be a subset of sites of the Hamming graph referred to as the subsystem 1. The projection operator over this subsystem is

$$\pi_{SV} = \sum_{v \in SV} |v\rangle \langle v|. \qquad (22)$$

We shall refer to its complement $F_q^d \setminus SV$ as the subsystem 2. In the ground state, the reduced density matrix of the subsystem 1 is defined by

$$\rho_1 = \text{tr}_2 |\Psi_0\rangle\langle\Psi_0| \qquad (23)$$

and its von Neumann entropy $S$ is

$$S = -\text{tr}(\rho_1 \ln \rho_1). \qquad (24)$$
This quantity allows one to determine to which degree the state of SV is interwined with the rest of the system. Once the eigenvalues of $\rho_1$ are determined, its computation is easy to realize. It is known that these eigenvalues can be extracted from those of the chopped correlation matrix $C$ [24, 26], which is defined as

$$C = |\hat{C}_{vv'}|_{vv' \in SV} \qquad (25)$$

and given by

$$C = \pi_{SV} \pi_{SE} \pi_{SV}, \qquad (26)$$

in terms of the projectors (21) and (22). In a nutshell, since the ground state is a Slater determinant, all the correlations can be expressed in terms of two-particle functions. This imply by Wick's theorem that the reduced density matrix $\rho_1$ is expressible as the exponential of a quadratic operator referred to as the entanglement hamiltonian $\mathcal{H}$:

$$\rho_1 = \kappa \exp(-\mathcal{H}), \quad \text{where} \quad \mathcal{H} = \sum_{v,v' \in SV} h_{vv'} c_v^\dagger c_{v'}. \qquad (27)$$

Since $C_{vv'} = \text{tr}(\rho_1 c_v^\dagger c_{v'})$, we see that $C$ and $h$ can be simultaneously diagonalized and that their spectra are related by the relation

$$h = \ln[(1-C)/C]. \qquad (28)$$

Following [6], we can give the entropy in terms of the eigenvalues $\lambda$ of $C$ and their degeneracy $D_\lambda$:

$$S = - \sum_{\lambda} D_{\lambda} [\lambda \ln(\lambda) + (1-\lambda) \ln(1-\lambda)]. \qquad (29)$$

In the following, we aim to diagonalize $C$ in two interesting cases.

# 4 Entanglement entropy for Hamming subgraphs

In this section, we take as subsystem 1 a subset of vertices which forms a Hamming graph of lower dimension. Take the full system to be the Hamming graph $H(d, q)$ associated to tuples of length $d$ constructed from elements in $\{0, 1, \dots, q-1\}$. A subset which corresponds to a Hamming subgraph is obtained by fixing some elements in the $d$-tuples. Here, without loss of generality, we fix the $d-L$ first elements to zero. The subset of $d$-tuples $SV_L$ corresponding to the subsystem 1 is thus

$$SV_L = \{v \in F_q^d : v_i = 0 \ \forall i \le d-L\}. \qquad (30)$$

For instance, in the case $d=3$ and $q=2$, we could consider the subset of binary tuples starting with one 0. This is the same as considering one of the faces of a 3-cube (see Figure 1).
<figure><img src="image_2.png" /><figcaption>Figure 1: A three dimensional cube and one of its Hamming subgraph (in black) that corresponds to $SV_2$, a two dimensional cube.</figcaption></figure>

The restriction of the action of $A$ on the vertices labeled by these $d$-tuples, i.e. $\pi_{SV_L} A \pi_{SV_L}$, corresponds to the action of the adjacency matrix of a Hamming graph associated to tuples of length $L$. We want to diagonalize the chopped correlation matrix for this choice of subsystem. The projection operator onto $SV_L$ can be expressed as

$$ \pi_{SV_L} = \underbrace{|0\rangle \langle 0| \otimes |0\rangle \langle 0| \otimes \cdots \otimes |0\rangle \langle 0|}_{d-L \text{ times}} \otimes \underbrace{\mathbb{1}_{q \times q} \otimes \cdots \otimes \mathbb{1}_{q \times q}}_{L \text{ times}}. \quad (31) $$

Using (14), it is also possible to give an explicit expression for the projection operator onto the energy eigenspaces contained in the ground state $\pi_{SE} = \sum_{k \in SE} E_{d-k}$:

$$ \pi_{SE} = \sum_{k \in SE} \sum_{\substack{b \in F_2^d, \\ \text{wt}(b)=d-k}} (b_1 \mathbb{1} + (-1)^{b_1} \frac{J_{q \times q}}{q}) \otimes \cdots \otimes (b_d \mathbb{1} + (-1)^{b_d} \frac{J_{q \times q}}{q}). \quad (32) $$

From these expressions, one can deduce the general form of the eigenvectors of $\pi_{SV_L} \pi_{SE} \pi_{SV_L}$. Indeed, let us recall that $\{|\theta_i\rangle\}$ with $i \in \{1, 2, \dots, q\}$ is a set of eigenvectors of $J_{q \times q}$ and let us note that

$$ \langle 0 | (b_i \mathbb{1} + (-1)^{b_i} \frac{J_{q \times q}}{q}) | 0 \rangle = \begin{cases} \frac{1}{q} & \text{if } b_i = 0, \\ \frac{q-1}{q} & \text{if } b_i = 1, \end{cases} \quad (33) $$

and that

$$ |\theta_q\rangle \langle\theta_q| = \frac{J_{q \times q}}{q} \quad \text{and} \quad \sum_{i=1}^{q-1} |\theta_i\rangle \langle\theta_i| = \mathbb{1} - \frac{J_{q \times q}}{q}. \quad (34) $$

Thus, if we define the following vectors in $\pi_{SV_L}(\mathbb{C}^q)^{\otimes d}$:

$$ |\theta_{i_1} \dots \theta_{i_L}\rangle = \underbrace{|0\rangle \otimes |0\rangle \otimes \cdots \otimes |0\rangle}_{d-L \text{ times}} \otimes |\theta_{i_1}\rangle \otimes \cdots \otimes |\theta_{i_L}\rangle, \quad (35) $$
we can see that

$$ \begin{aligned} \pi_{SV_L} E_{d-k} \pi_{SV_L} |\theta_{i_1} \dots \theta_{i_L} \rangle &= \pi_{SV_L} E_{d-k} |\theta_{i_1} \dots \theta_{i_L} \rangle \\ &= \begin{cases} \left(\frac{1}{q}\right)^{k-N_q} \left(\frac{q-1}{q}\right)^{d-L-k+N_q} \binom{d-L}{k-N_q} |\theta_{i_1} \dots \theta_{i_L} \rangle & \text{if } d-L \ge k-N_q \ge 0, \\ 0 & \text{otherwise,} \end{cases} \end{aligned} \quad (36) $$

where $N_q \equiv \text{number of } |\theta_q\rangle \text{ in } |\theta_{i_1} \dots \theta_{i_L}\rangle$. Therefore, the action of the chopped correlation matrix on these vectors is

$$ C |\theta_{i_1} \dots \theta_{i_L} \rangle = \sum_{\substack{k \in SE \\ d-L \ge k-N_q \ge 0}} \left(\frac{1}{q}\right)^{k-N_q} \left(\frac{q-1}{q}\right)^{d-L-k+N_q} \binom{d-L}{k-N_q} |\theta_{i_1} \dots \theta_{i_L} \rangle. \quad (37) $$

For all $N_q \in \{0, 1, \dots, L\}$, we thus have an eigenvalue $\lambda_{N_q}$ of C, given by the following expression:

$$ \lambda_{N_q} = \sum_{\substack{k \in SE \\ d-L+N_q \ge k \ge N_q}} \left(\frac{1}{q}\right)^{k-N_q} \left(\frac{q-1}{q}\right)^{d-L-k+N_q} \binom{d-L}{k-N_q}. \quad (38) $$

As for its degeneracy $\binom{L}{N_q}(q-1)^{L-N_q}$, it corresponds to the number of vectors $|\theta_{i_1} \dots \theta_{i_L}\rangle$ having the same $N_q$. This result can then be used as an input in (29) to obtain the entanglement entropy $S$.

In order to provide an explicit example, let us restrict ourselves to the case $SE = \{0, 1, \dots, k_0\}$. When $N_q > k_0$, we notice that $SE \cap \{N_q, N_q+1, \dots, d-L+N_q\}$ is empty and so $\lambda_{N_q}$ is 0. Moreover, formula (38) corresponds to a sum over all the terms in a binomial distribution when $d-L+N_q < k_0$. This implies that $\lambda_{N_q} = 1$. One can check that in these two cases, $\lambda_{N_q}$ does not contribute to (29). We are left with $d-L+1$ eigenvalues to consider and obtain:

$$ S = - \sum_{i=0}^{d-L} \binom{L}{d-k_0-i} (q-1)^{d-k_0-i} [F(i) \ln F(i) + (1-F(i)) \ln(1-F(i))], \quad (39) $$

where $F$ is the cumulative binomial distribution function

$$ F(i) = F(i; d-L, 1/q) = \sum_{i'=0}^{i} \binom{d-L}{i'} \left(\frac{1}{q}\right)^{i'} \left(\frac{q-1}{q}\right)^{d-L-i'} . \quad (40) $$

In terms of Krawtchouk polynomials, we also have

$$ F(i) = 1 - \left(\frac{1}{q}\right)^{i+1} \left(\frac{q-1}{q}\right)^{d-L-i-1} \binom{d-L}{i+1} K_{d-L-i-1}(-1; 1-q, i+2). \quad (41) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Entanglement entropy for cubes of dimension $L$ in cubes of dimension $d$. (a): entropy at $k_0 = d/2$ (half-filling) for sub-cubes of dimension $L = d/4$, $L = d/2$ and $L = 3d/4$. (b) Ratio of the entropy over the dimension $2^L$ of the sub-cubes ($L = d/4$) for different filling ratios of the Fermi sea.</figcaption></figure>

In Figure 2, the entropy of sub-cubes in hypercubes are presented for different size ratios $L/d$ and filling ratios $k_0/d$. The left figure suggests that the entropy of Hamming subgraphs grows with their volume $|SV| = 2^L$, so that their state are nearly perfectly mixed. This does not contradict the area law, since every site in a Hamming subgraph has at least one neighbor outside the subsystem and is thus on the boundary. The figure on the right shows that $S$ peaks at half-filling $k_0/d$.

# 5 Entanglement entropy for neighborhoods

Let us pick a reference tuple $v_0$ and consider all the sites labeled by tuples $v$ at a common distance $i$ from $v_0$. This set of sites is called the $i^{th}$ neighborhood of $v_0$. In this section, we want to find the entanglement entropy of subsystems which are bundles of neighborhoods of $v_0$. For a subset of distances $SD \subset \{0, 1, \dots, d\}$, we take the subsystem 1 to be the set of sites $v$ for which the distance from $v_0$ is in $SD$, i.e.

$$ SV = \{v \in F_q^d : \partial(v_0, v) \in SD\}. \qquad (42) $$

For instance, we could consider $H(3, 2)$ with $v_0 = (0, 0, 0)$ and $SD = \{1\}$. Then, the subsystem would be the sites labeled by $(0, 0, 1)$, $(0, 1, 0)$ or $(1, 0, 0)$ (see Figure 3).
<figure><img src="image_4.png" /><figcaption>Figure 3: A three dimensional cube and, in black, the 1st neighborhood of the vertex (0,0,0) (in gray).</figcaption></figure>

Now, without loss of generality, we take $v_0 = (0,0,0,\dots,0)$. The projection operator onto the space associated to vertices at a distance $i$ from $v_0$ is defined as

$$E_i^* = \sum_{\substack{v \in F_q^d \\ \partial(v_0, v) = i}} |v\rangle \langle v| . \qquad (43)$$

The projector over the subsystem is

$$\pi_{SV} = \sum_{\substack{v \in F_q^d \\ \partial(v_0, v) \in SD}} |v\rangle \langle v| = \sum_{i \in SD} E_i^* . \qquad (44)$$

Thus, the chopped correlation matrix is

$$C = \sum_{i,i' \in SD} \sum_{k \in SE} E_i^* E_{d-k} E_{i'}^* . \qquad (45)$$

To diagonalize this chopped correlation matrix, we can take advantage of the fact that $C$ is expressed only in terms of elements in the Terwilliger algebra of the Hamming scheme (denoted by $\mathcal{T}$). This means that the vector space $(\mathbb{C}^q)^{\otimes d}$ on which the operators $E_k$ and $E_i^*$ are acting corresponds to a $\mathcal{T}$-module and can be decomposed into irreducible submodules. This decomposition simplifies the problem, since it allows to work with one irreducible module at a time. In the next subsections, we give a quick overview of $\mathcal{T}$ and show how the decomposition in irreducible components of the representation on $(\mathbb{C}^q)^{\otimes d}$ is carried out.

## 5.1 The Terwilliger algebra of the Hamming scheme

Since the Hamming graphs are distance-regular, the following relations between adjacency matrices and projection operators are verified:

$$\bullet \ A_0 = 1 \text{ and } E_0 = J_{q^d \times q^d} / q^d ;$$
• $\sum_{i=0}^{d} A_i = J_{q^d \times q^d}$ and $\sum_{i=0}^{d} E_i = 1$;

• $A_i \circ A_j = \delta_{ij}A_i$ and $E_i E_j = \delta_{ij}E_i$;

• $A_i A_j = \sum_{k=0}^{d} p_{ij}^k A_k$ and $E_i \circ E_j = \frac{1}{q^d} \sum_{k=0}^{d} q_{ij}^k E_k$,

where $(A \circ B)_{ij} = A_{ij}B_{ij}$ is the entry-wise product. The commuting algebra spanned by the adjacency matrices $A_i$ is in the present case the Bose-Mesner algebra of the Hamming scheme. It is important to note that the set of primitive idempotents $E_i$ also offers a basis for this algebra. Indeed, there exists coefficients $p_i(j)$ and $q_i(j)$ such that

$$A_i = \sum_{j=0}^{d} p_i(j) E_j \quad \text{and} \quad E_i = \frac{1}{|X|} \sum_{j=0}^{d} q_i(j) A_j, \quad (46)$$

where $|X|$ is the dimension of the matrices and is equal to $q^d$ for the Hamming scheme. From there, we can choose a reference vertex $v_0$ and introduce dual matrices $A_i^*(v_0)$ and $E_i^*(v_0)$. The dual adjacency matrices and dual primitive idempotents of the Hamming scheme are hence

$$[A_i^*(v_0)]_{vv} = q^d [E_i]_{v_0 v} \quad \text{and} \quad [E_i^*(v_0)]_{vv} = [A_i]_{v_0 v}. \quad (47)$$

In the following, we use a simplified notation: $A_i^* = A_i^*(v_0)$ and $E_i^* = E_i^*(v_0)$. These matrices obey relations similar to the ones of their counterpart $A_i$ and $E_i$:

$$A_i^* A_j^* = \sum_{k=0}^{d} q_{ij}^k A_k^* \quad \text{and} \quad E_i^* E_j^* = \delta_{ij} E_i^*. \quad (48)$$

One can check that the definition of $E_i^*$ given above is equivalent to (43). In the Hamming scheme, $A^*$ is given by the following tensor products:

$$A^* = q \sum_{i=1}^{d} \left( \underbrace{1 \otimes \cdots \otimes 1}_{i-1 \text{ times}} \otimes |0\rangle \langle 0| \otimes \underbrace{1 \otimes \cdots \otimes 1}_{i-1} \right) - d. \quad (49)$$

It is easy to check that it acts diagonally on the set of vectors $|v\rangle$, $v \in F_q^d$, associated to sites in the graph.

The Terwilliger algebra is the algebra spanned by the set of adjacency matrices $\{A_0, A_1, \dots, A_d\}$ and the set of dual adjacency matrices $\{A_0^*, A_1^*, \dots, A_d^*\}$ [29, 30, 31]. In terms of the primitive idempotents, it can also be thought of as the algebra spanned by $\{E_0, E_1, \dots, E_d\}$ and $\{E_0^*, E_1^*, \dots, E_d^*\}$. Looking back at the definition of the chopped correlation matrix (45), we see that $C$ is indeed part of the Terwilliger algebra of the Hamming scheme.
## 5.2 Extracting the irreducible modules

We now show how to decompose $(\mathbb{C}^q)^{\otimes d}$ into the irreducible modules of the Terwilliger algebra of the Hamming scheme. This will allow to diagonalize $C$ one irreducible module at a time.

This decomposition has been studied in the past [15, 22]. It has been applied to the computation of entanglement entropy of coupled harmonic oscillators on Hamming networks in [17, 18] where the Schmidt numbers and thus the entropy could be extracted from the potential matrix of the system expressed in the basis associated to the decomposition.

### 5.2.1 Irreducible representations for the hypercube $H(d, 2)$

First, we consider the case $H(d, 2)$, i.e. the Terwilliger algebra of hypercubes. We see then from (7) and (51) that the adjacency and dual adjacency matrices are given in terms of Pauli matrices:

$$A = \sum_{i=1}^{d} \left( \underbrace{\mathbf{1} \otimes \cdots \otimes \mathbf{1}}_{i-1 \text{ times}} \otimes \sigma^x \otimes \mathbf{1} \otimes \cdots \otimes \mathbf{1} \right) \quad (50)$$

and

$$A^* = \sum_{i=1}^{d} \left( \underbrace{\mathbf{1} \otimes \cdots \otimes \mathbf{1}}_{i-1 \text{ times}} \otimes \sigma^z \otimes \mathbf{1} \otimes \cdots \otimes \mathbf{1} \right). \quad (51)$$

$A$ and $A^*$ can thus be interpreted in terms of elements in the representation of $\mathfrak{su}(2)$ constructed from the tensor product of $d$ irreducible representations of spin $1/2$. Indeed, if $s_i^x$, $s_i^y$ and $s_i^z$ are the following operators:

$$s_i^j = \underbrace{\mathbf{1} \otimes \cdots \otimes \mathbf{1}}_{i-1 \text{ times}} \otimes \frac{\sigma^j}{2} \otimes \mathbf{1} \otimes \cdots \otimes \mathbf{1}, \quad \text{with } j \in \{x, y, z\}, \quad (52)$$

giving for each $i$ a two dimensional representation of $\mathfrak{su}(2)$, it follows that the adjacency and dual adjacency matrices are proportional to the components $s^x$ and $s^z$ of the total spin operator, i.e.

$$A = 2 \sum_{i=1}^{d} s_i^x = 2s^x \quad \text{and} \quad A^* = 2 \sum_{i=1}^{d} s_i^z = 2s^z. \quad (53)$$

From this, one deduces that the irreducible modules of the Terwilliger algebra of $H(d, 2)$ are the vector spaces associated to the irreducible representations of $\mathfrak{su}(2)$ contained in the $d$-fold product of the fundamental representation of $\mathfrak{su}(2)$. For $d$ even (resp. odd) we know from the standard Clebsh-Gordan decomposition that for each $j$ in $\{0, 1, \dots, \frac{d}{2}\}$ (resp. in $\{1/2, 3/2, \dots, \frac{d}{2}\}$) there
exists $\frac{2j+1}{d+1} \binom{d+1}{2-j}$ subspaces of $(\mathbb{C}^2)^{\otimes d}$ spanned by vectors $\{|j, m\rangle\}_{-j \le m \le j}$ such that

$$A|j, m\rangle = \sqrt{(j+m+1)(j-m)}|j, m+1\rangle + \sqrt{(j-m+1)(j+m)}|j, m-1\rangle \quad (54)$$

and

$$A^*|j, m\rangle = 2m|j, m\rangle \quad (55)$$

The label necessary to distinguish the different spaces of dimension $2j + 1$ is kept implicit. It is also important to note that since the $|j, m\rangle$ are eigenvectors of $A^*$ with eigenvalue $2m$, they are contained in the $(\frac{d}{2} - m)$-th neighborhood of the hypercube, i.e. they are linear combinaisons of vectors $|v\rangle$ of weight $\mathrm{wt}(v) = \frac{d}{2} - m$.

### 5.2.2 Irreducible representations for the Hamming graph $H(d, q)$

Now that we have the decomposition for $q = 2$, it can be used to obtain a similar result for a general Hamming graph $H(d, q)$. More precisely, we want to identify subspaces of $(\mathbb{C}^q)^{\otimes d}$ on which the action of $A$ and $A^*$ is similar to their action in the case of a hypercube. First, let us define a set of vectors $\{|\tilde{n}\rangle\}_{0 \le n < q}$ which gives a basis of $\mathbb{C}^q$:

$$|\tilde{0}\rangle = |0\rangle, \quad |\tilde{1}\rangle = \frac{1}{\sqrt{q-1}} \sum_{i=1}^{q-1} |i\rangle; \quad (56)$$

and for $n > 1$, the degeneracy of the eigenspace of $J$ associated to the eigenvalue 0 allows to take them such that

$$\langle \tilde{n} | \tilde{m} \rangle = \delta_{mn}, \quad \langle 0 | \tilde{n} \rangle = 0 \quad \text{and} \quad J | \tilde{n} \rangle = 0. \quad (57)$$

When $q = 2$, we notice that $|\tilde{1}\rangle = |1\rangle$ and that this basis matches with the usual one. Furthermore, for any $v \in F_q^d$, we can define $|\tilde{v}\rangle = |\tilde{v}_1\rangle \otimes |\tilde{v}_2\rangle \otimes \dots \otimes |\tilde{v}_d\rangle$, which gives a new basis for $(\mathbb{C}^q)^{\otimes d}$.

Let us divide this basis of $(\mathbb{C}^q)^{\otimes d}$ in subsets. For two $d$-tuples $v$ and $v'$ in $F_q^d$, we say that $v \sim v'$ if $\forall i \in \{1, \dots, d\}$ we either have (1) $v_i = v'_i$, (2) $v_i = 0$ and $v'_i = 1$ or (3) $v'_i = 0$ and $v_i = 1$. In other words, we say that two tuples are equivalent if they differ only on their binary part. We refer to the set of $d$-tuples equivalent to $v$ as $[v]$. We can define the subspace of $(\mathbb{C}^q)^{\otimes d}$ generated by all the vectors associated to tuples in $[v]$ as $V_{[v]} = \mathrm{span}\{|\tilde{v}'\rangle\}_{v' \in [v]}$. Thus, we have the following decomposition of our vector space:

$$(\mathbb{C}^q)^{\otimes d} = \bigoplus_{v \in F_q^d / \sim} V_{[v]}. \quad (58)$$
We are interested in characterizing the action of A on a subspace $V_{[v]}$. With no loss of generality, we pick a tuple v for which the first n elements are binary. Since d-tuples in $[v]$ differ only at positions where they take the value 0 or 1, we can refer to the basis vectors $|\tilde{v}'\rangle$ of $V_{[v]}$ with a binary n-tuples $b = (b_1, \dots, b_n)$. Indeed, we have

$$ \begin{aligned} |\tilde{v}'\rangle &= \underbrace{|\tilde{v}'_1\rangle \otimes \dots \otimes |\tilde{v}'_n\rangle}_{\text{binary part}} \otimes \underbrace{|\tilde{v}'_{n+1}\rangle \otimes \dots \otimes |\tilde{v}'_d\rangle}_{\text{fixed by the choice of } [v]} \quad (59) \\ &= |\tilde{b}_1\rangle \otimes \dots \otimes |\tilde{b}_n\rangle \otimes |\psi\rangle \equiv |b\rangle, \end{aligned} $$

where $|\psi\rangle = |\tilde{v}'_{n+1}\rangle \otimes \dots \otimes |\tilde{v}'_d\rangle$ is the fixed part. To find the action of A on a vector $|b\rangle$, we notice that

$$ (J-1)|\tilde{0}\rangle = \sqrt{q-1}|\tilde{1}\rangle, \quad (J-1)|\tilde{1}\rangle = \sqrt{q-1}|\tilde{0}\rangle + (q-2)|\tilde{1}\rangle, \quad (60) $$

and that for $n > 1$

$$ (J-1)|\tilde{n}\rangle = -|\tilde{n}\rangle. \quad (61) $$

Therefore, we can express the action of A in terms of tensor products of Pauli matrices in the following way:

$$ \begin{aligned} A|b\rangle = & \left[ \sqrt{q-1} \sum_{i=1}^{n} \left( \underbrace{\mathbf{1} \otimes \dots \otimes \mathbf{1}}_{i-1 \text{ times}} \otimes \sigma^x \otimes \mathbf{1} \otimes \dots \otimes \mathbf{1} \right) |\tilde{b}_1\rangle \otimes \dots \otimes |\tilde{b}_n\rangle \right] \otimes |\psi\rangle \\ & - \left[ \frac{(q-2)}{2} \sum_{i=1}^{n} \left( \underbrace{\mathbf{1} \otimes \dots \otimes \mathbf{1}}_{i-1 \text{ times}} \otimes \sigma^z \otimes \mathbf{1} \otimes \dots \otimes \mathbf{1} \right) |\tilde{b}_1\rangle \otimes \dots \otimes |\tilde{b}_n\rangle \right] \otimes |\psi\rangle \\ & + \left( \frac{qn}{2} - d \right) |b\rangle, \end{aligned} \quad (62) $$

where

$$ \sigma^x |\tilde{0}\rangle = |\tilde{1}\rangle, \quad \sigma^x |\tilde{1}\rangle = |\tilde{0}\rangle, \quad \sigma^z |\tilde{0}\rangle = |\tilde{0}\rangle \quad \text{and} \quad \sigma^z |\tilde{1}\rangle = -|\tilde{1}\rangle. \quad (63) $$

For the action of $A^*$ on $V_{[v]}$, we see from (51) that:

$$ A^*|b\rangle = q(n - \text{wt}(b)) - d|b\rangle. \quad (64) $$

At this point, we can apply the results obtained for $q = 2$. Thus, given a subspace $V_{[v]}$ characterized by a n even (resp. odd), we know that for each j in $\{0, 1, \dots, \frac{n}{2}\}$ (resp. in $\{1/2, 3/2, \dots, \frac{n}{2}\}$) there exists $\frac{2j+1}{n+1} \binom{n+1}{\frac{n}{2}-j}$ subspaces spanned by vectors $\{|j, m\rangle\}_{-j \le m \le j}$ such that

$$ \begin{aligned} A|j, m\rangle = & \sqrt{(q-1)(j+m+1)(j-m)} |j, m+1\rangle \\ & + \sqrt{(q-1)(j+m)(j-m+1)} |j, m-1\rangle \\ & + \left( \frac{qn}{2} - d - (q-2)m \right) |j, m\rangle \end{aligned} \quad (65) $$
and

$$A^* |j, m\rangle = qm + \frac{nq}{2} - d |j, m\rangle. \qquad (66)$$

Let $\ell$ label the degenerate subspaces of dimension $2j+1$. We see that $K_{j,\ell,[v]} = \text{span}\{|j,m\rangle : m \in \{-j, \dots, j\}\}$ corresponds to the irreducible $\mathcal{T}$-submodules. Therefore, the decomposition needed is:

$$(\mathbb{C}^q)^{\otimes d} = \bigoplus_{v \in F_q^d / \sim} \bigoplus_{j, \ell} K_{j, \ell, [v]}. \qquad (67)$$

For the general Hamming case, it is worth noting that $|j, m\rangle$ is in the $(d - \frac{n}{2} - m)$-th neighborhood of the graph and that each irreducible module $K_{j,\ell,[v]}$ contains no more than one vector per neighborhood. In algebraic terms, this amounts to say that $\mathcal{T}$ is *thin* and it is well known that the Terwilliger algebra of the Hamming scheme has this property [31].

### 5.2.3 Chopped correlation matrix for neighborhoods and entanglement entropy

Now, let us see how (67) simplifies the diagonalization of the chopped correlation matrix. First, we can consider

$$\underline{c}_{j,\ell,m,v}^{\dagger} = \sum_{v' \in F_q^d} \langle v' | j, m \rangle c_{v'}^{\dagger} \quad \text{and} \quad \underline{c}_{j,\ell,m,v} = \sum_{v' \in F_q^d} \langle j, m | v' \rangle c_{v'}^{\dagger}. \qquad (68)$$

One can check that these fermionic operators respect the canonical anticommutation relations. Expressing the Hamiltonian (9) in terms of these operators, one finds that it decomposes into a sum of Hamiltonians $\hat{\mathcal{H}}_{K_{j,\ell,[v]}}$ acting independently on orthogonal subspaces:

$$\begin{aligned} \hat{\mathcal{H}} &= \sum_{v \in F_q^d / \sim} \sum_{j, \ell} \sum_{m, m' = -j}^{j} \langle j, m | [\sum_i \alpha_i A_i] | j, m' \rangle \underline{c}_{j, \ell, m, v}^{\dagger} \underline{c}_{j, \ell, m', v} \\ &= \sum_{v \in F_q^d / \sim} \sum_{j, \ell} \hat{\mathcal{H}}_{K_{j, \ell, [v]}} . \end{aligned} \qquad (69)$$

Since each module contains at most one vector per neighborhood, the operators $\hat{\mathcal{H}}_{K_{j,\ell,[v]}}$ can be interpreted as Hamiltonians of free fermions on chains. Similarly for $C$, we have

$$\begin{aligned} C &= \pi_{SV} \pi_{SE} \pi_{SV} \\ &= \sum_{v \in F_q^d / \sim} \sum_{j, \ell} [\pi_{SV} \pi_{SE} \pi_{SV}]_{K_{j, \ell, [v]}} \\ &= \sum_{v \in F_q^d / \sim} \sum_{j, \ell} [C]_{K_{j, \ell, [v]}}, \end{aligned} \qquad (70)$$
where $[C]_{K_{j,\ell,[v]}}$ is the restriction of the chopped correlation matrix to the subspace $K_{j,\ell,[v]}$. The entries of $[C]_{K_{j,\ell,[v]}}$ can be obtained in the following way. For $|j, \omega_{k'}\rangle$ an eigenvector of $[A]_{K_{j,\ell,[v]}}$, we get a three terms recurrence relations for the overlap coefficients $Q_{m,k'} = \langle j, \omega_{k'}|j, m\rangle$ by considering the matrix element $\langle j, \omega_{k'}| [A]_{K_{j,\ell,[v]}} |j, m\rangle$. We have

$$ \omega_{k'} Q_{m,k'} = \sqrt{(q-1)(j+m+1)(j-m)} Q_{m+1,k'} + \left(\frac{nq}{2} - d - (q-2)m\right) Q_{m,k'} \\ + \sqrt{(q-1)(j-m+1)(j+m)} Q_{m-1,k'} . \qquad (71) $$

This leads to the identifications $\omega_{k'} = q(j - k' + \frac{n}{2}) - d$ and to

$$ Q_{m,k'} = \sqrt{\binom{2j}{k'}} \sqrt{\binom{2j}{j-m}} q^{-j} (q-1)^{\frac{j-m+k'}{2}} K_{j-m}(k'; (q-1)/q, 2j) . \quad (72) $$

Therefore, we find for $m, m' \in SD \cap \{-j, \dots, j\}$ that

$$ \langle j, m | [C]_{K_{j,\ell,[v]}} | j, m' \rangle = \sum_{j-k'+\frac{n}{2} \in SE} Q_{m,k'} Q_{m',k'} . \qquad (73) $$

Knowing all the entries of the matrix $[C]_{K_{j,\ell,[v]}}$, we can extract its eigenvalues, i.e. find the zeros of its characteristic polynomial of degree $|SD \cap \{d - \frac{n}{2} - j, \dots, d - \frac{n}{2} + j\}|$. This degree is at most the number of neighborhoods in the subsystem. It is in general significantly less than the dimension of $C$, which is equal to the number of sites in all the neighborhoods in the subsystem. In particular, when we consider a single neighborhood, this decomposition actually yields the diagonalization. Indeed, $[C]_{K_{j,\ell,[v]}}$ is then a 1-dimensional matrix for which the only entry is its eigenvalue. If $SV$ is the $i^{th}$ neighborhood, the eigenvalues of $C$ are

$$ \lambda_{n,j} = \sum_{k \in SE} Q_{d-\frac{n}{2}-i, j-k+\frac{n}{2}}^2 . \qquad (74) $$

The degeneracy is given by the product of the number of equivalence classes $[v]$ with $d$-tuples containing $n$ binary terms and the number of representations of $\mathfrak{su}(2)$ of dimension $2j + 1$ in the hypercube $H(n, 2)$. Thus, it is given by $D_{n,j} = (q-2)^{d-n} \frac{2j+1}{d+1} \binom{d+1}{j-n} \binom{d}{n}$. In Figure 4, the entropy of neighborhoods of hypercubes are presented for different distance ratios $i/d$ and filling ratios $k_0/d$. The left figure suggests that the entropy of a neighborhood grows with its volume $|SV|$ and thus that their state are nearly perfectly mixed. The figure on the right shows that $S$ peaks at half-filling $k_0/d$. We recall that similar results were obtained for Hamming subgraphs.
<figure><img src="image_5.png" /><figcaption>Figure 4: Entanglement entropy for neighborhoods of d-cubes. (a): entropy at $k_0 = d/2$ (half-filling) for the neighborhoods $i = d/8, i = d/4$ and $i = 3d/8$. (b) Ratio of the entropy over the dimension of the neighborhood $i = d/4$ for different filling ratios of the Fermi sea.</figcaption></figure>

### 5.2.4 Interlude : $\mathcal{T}$ as a quotient of the Onsager algebra

Before considering the case where $SV$ is made of a large number of neighborhoods, let us make an additional remark about $\mathcal{T}$. It is interesting to note that the decomposition of the vector space in $\mathcal{T}$-submodules has brought up a relation between the Terwilliger algebra of $H(d, q)$ and a quotient of the Onsager algebra $O[2, 14, 23, 27]$. To see this, we first highlight the fact that the relation between $\mathcal{T}$ and $su(2)$ extends to $q \neq 2$. In terms of the matrices $s^x$ and $s^z$ from an irreducible representation of dimension $2j + 1$, we note that

$$ [A]_{K_{j,\ell,[v]}} = 2\sqrt{q-1}s^x - (q-2)s^z + \left(\frac{dq}{2} - d\right)\mathbf{1}, \quad (75) $$

and

$$ [A^*]_{K_{j,\ell,[v]}} = qs^z + \left(\frac{qn}{2} - d\right)\mathbf{1}. \quad (76) $$

Using these expression and the commutation relations of $su(2)$, one can easily check that the dual and non-dual adjacency matrices of the Hamming scheme verify the following relations:

$$ [A, [A, [A, A^*]]] = q^2[A, A^*] \quad \text{and} \quad [A^*, [A^*, [A^*, A]]] = q^2[A^*, A]. \quad (77) $$

These were first obtained by Terwilliger in [31]. Taking $A_0 = \frac{4}{q}A$ and $A_1 = \frac{4}{q}A^*$, one also sees that these operators verify the Dolan-Grady relations [13]:

$$ [A_0, [A_0, [A_0, A_1]]] = 16[A_0, A_1] \quad \text{and} \quad [A_1, [A_1, [A_1, A_0]]] = 16[A_1, A_0], \quad (78) $$
which serve as defining relations for the Onsager algebra $\mathcal{O}$. This algebra is in general of infinite dimension and the rest of its generators, $\mathcal{A}_n$ and $\mathcal{G}_n$ with $n \in \mathbb{Z}$, are defined by the equivalent presentation [23]:

$$[\mathcal{A}_n, \mathcal{A}_m] = 4\mathcal{G}_{n-m}, \quad [\mathcal{G}_n, \mathcal{A}_m] = 2\mathcal{A}_{n+m} - 2\mathcal{A}_{m-n} \quad \text{and} \quad [\mathcal{G}_n, \mathcal{G}_m] = 0. \quad (79)$$

We see that $\mathcal{T}$ gives a representation of $\mathcal{O}$. It is however of finite dimension. Indeed, one can check using (75) and (76) that the following Davies relations hold [12]:

$$\mathcal{A}_{p+2} - \mathcal{A}_{p-1} + \frac{q+4}{q}(\mathcal{A}_{p+1} - \mathcal{A}_p) = 0 \quad (80)$$

and

$$\mathcal{G}_{p+2} - \mathcal{G}_{p-1} + \frac{q+4}{q}(\mathcal{G}_{p+1} - \mathcal{G}_p) = 0, \quad (81)$$

for all $p \in \mathbb{Z}$. In other words, the generators can all be expressed as linear combinations of $\{\mathcal{A}_{-1}, \mathcal{A}_0, \mathcal{A}_1, \mathcal{G}_1\}$. Using the expression of these generators in terms of $\{s^x, s^y, s^z, (\frac{nq}{2} - d)\mathbf{1}\}$, one can also check that they are linearly independent. Thus, $\mathcal{T}$ for $H(d, q)$ yields a Lie algebra isomorphic to the quotient of $\mathcal{O}$ by (80) and (81).

## 5.3 Entanglement entropy for multiple neighborhoods

The decomposition in irreducible modules has simplified the diagonalization of the chopped correlation matrix $C$ by allowing to work with matrices of lower dimension. However, we may still have to diagonalize large matrices with few non-zero entries and many eigenvalues approaching 0 or 1 when the subsystem is composed of many neighborhoods. This happens for instance if we take as subsystem 1 all the vertices at a distance lower than some large integer $N$ from a given site. Then,

$$\pi_{SV} = \sum_{i=0}^{N} E_{i}^{*}. \quad (82)$$

In the spirit of the celebrated time and band limiting approach [16, 20, 28], we would now wish to find a simple operator $T$ with practical diagonalization properties such that

$$[C, T] = 0 \quad (83)$$

and that it thus shares with $C$ common eigenvectors.

### 5.3.1 The generalized algebraic Heun operator

Since we are working with a distance-regular graph, we can adopt the approach based on Heun operators developed in [8, 9, 10] to construct $T$. This leads us to
look at the most general symmetric block-tridiagonal operator in the Terwilliger algebra of the Hamming scheme:

$$T = \{A, A^*\} + \mu A^* + \nu A, \qquad (84)$$

These are referred to as generalized algebraic Heun operators. We want to fix $\mu$ and $\nu$ to allow $T$ to commute with both $\pi_{SV}$ and $\pi_{SE}$, assuring that $T$ also commutes with $C$. For a given irreducible module $K_{j,\ell,[\nu]}$, the action of $T$ is tridiagonal on the eigenbases of both $A$ and $A^*$:

$$T|j, m\rangle = (\nu + 2qm + nq + q - 2d)\sqrt{(q-1)(j+m+1)(j-m)}|j, m+1\rangle \\ + (\nu + 2qm + nq - q - 2d)\sqrt{(q-1)(j+m)(j-m+1)}|j, m-1\rangle \\ + \left[(\nu + 2qm + nq - 2d)\left(\frac{qn}{2} - d - (q-2)m\right) + \mu\left(qm + \frac{nq}{2} - d\right)\right]|j, m\rangle \qquad (85)$$

and

$$T|j, \omega_{k'}\rangle = (\mu + \omega_{k'} + \omega_{k'+1})\sqrt{(q-1)(2j-k')(k'+1)}|j, \omega_{k'+1}\rangle \\ + (\mu + \omega_{k'} + \omega_{k'-1})\sqrt{(q-1)(k')(2j-k'+1)}|j, \omega_{k'-1}\rangle \\ + \left[2\omega_{k'}\left(\frac{qn}{2} - d - (q-2)(j-k')\right) + \nu\omega_{k'}\right]|j, \omega_{k'}\rangle, \qquad (86)$$

where we recall that $\omega_{k'} = q(j - k' + \frac{n}{2}) - d$. As for the action of the projection operators on these vectors, they are given by

$$\pi_{SV}|j, m\rangle = \begin{cases} |j, m\rangle & \text{if } d - \frac{n}{2} - m \le N \\ 0 & \text{otherwise} \end{cases} \qquad (87)$$

and

$$\pi_{SE}|j, \omega_{k'}\rangle = \begin{cases} |j, \omega_{k'}\rangle & \text{if } j - k' + \frac{n}{2} \le k_0 \\ 0 & \text{otherwise.} \end{cases} \qquad (88)$$

From (85) and (87), we see that $[T, \pi_{SV}] = 0$ if $(\nu + 2q(d - N) - q - 2d) = 0$. Similarly, from (86) and (88), we see that $[T, \pi_{SE}] = 0$ if $(\mu + 2qk_0 - q - 2d) = 0$. Thus, if we take

$$\nu = -2q(d - N) + q + 2d \quad \text{and} \quad \mu = -2qk_0 + q + 2d, \qquad (89)$$

we see that

$$[T, C] = [T, \pi_{SV}\pi_{SE}\pi_{SV}] = 0. \qquad (90)$$

Since the choice of $\mu$ and $\nu$ does not depend on $n$ or $j$, (90) holds for all the irreducible modules. To diagonalize $T$, we recall (75) and (76) which imply

$$[T]_{K_{j,\ell,[\nu]}} = 2q\sqrt{q-1}\{s^x, s^z\} - q(q-2)\{s^z, s^z\} + 2\sqrt{q-1}\left(2\left(\frac{nq}{2} - d\right) + \nu\right)s^x \\ + (4\left(\frac{nq}{2} - d\right) + \mu q - \nu(q-2))s^z + \left(\frac{nq}{2} - d\right)\left(2\left(\frac{nq}{2} - d\right) + \mu + \nu\right). \qquad (91)$$
This operator has a well behaved spectrum, making its eigenvectors easy to obtain numerically. Computing the entropy is then a matter of acting on these vectors with the chopped correlation matrix to extract its spectrum and applying (29). Figure 5 presents results obtained with this method in the case of hypercubes containing up to $2^{120}$ sites. As expected, the entropy seems to verify an area law as we approach the thermodynamic limit.

<figure><img src="image_6.png" /><figcaption>Figure 5: Entanglement entropy of the sites at a distance $\le N$ from a fixed vertex. (a): entropy at $k_0 = d/2$ (half-filling) for $N = d/8$, $N = d/4$ and $N = d/2$. (b) Ratio of the entropy $S$ over $|\partial SV|$ (the number of sites at the boundary of the subsystem) for different filling ratios of the Fermi sea.</figcaption></figure>

The diagonalization of $T$ can also be approached analytically. Indeed, the operator (91) corresponds to a BC-Gaudin magnet Hamiltonian in a magnetic field and can be diagonalized using the modified algebraic Bethe ansatz. This is briefly presented in the following subsection (for more details, see [4, 11]).

### 5.3.2 The modified algebraic Bethe ansatz and the spectrum of $T$

In a nutshell, the modified algebraic Bethe ansatz works in the following way. From a solution $r$ of the non-standard classical Yang-Baxter equation, one constructs a $K$-matrix which satisfies the classical reflection equation and defines a transfer matrix $t(u) = \text{tr}K(u)^2$. This matrix $t(u)$ satisfies $[t(u), t(v)] = 0$, that is the matrices associated to different values of the parameter commute. It is known how to find solutions $r$ and $K$ for which $t(u)$ can be identified as a BC-Gaudin magnet Hamiltonian in a magnetic field (and thus, as $T$) when $u \to 0$. It is also known that such transfer matrix can be diagonalized by vectors
of the following form:

$$|z\rangle = B(z_1, 1)B(z_2, 2)\dots B(z_{2j+1}, 2j+1)|j, -j\rangle, \quad (92)$$

where $B(z, n)$ are operators that can be written in terms of $su(2)$ generators and are constructed to satisfy particular commutation relations with the entries of the $K$-matrix. For the ansatz (92) to be an eigenvector of $t(u)$, a set of equations depending on $z$ needs to be verified. These are the Bethe equations that are given by:

$$\frac{e^{2i\theta} j(z_k^2 - 1)}{(e^{2i\theta} z_k - 1)(e^{2i\theta} - z_k)} + \frac{(z_k^2 + 1)(1 - \frac{\nu}{2}) - \mu z_k}{(z_k^2 - 1)} + \sum_{\substack{p=1 \\ p \neq k}}^{M} \frac{z_p(z_k^2 - 1)}{(z_k - z_p)(z_k z_p - 1)} = 0, \quad (93)$$

where $\theta = \frac{1}{2} \arccot \frac{\sqrt{q-1}}{q-2}$ and the dimension of $\pi_{SV} K_{j,\ell,[v]} \pi_{SV}$ is $M = N - d + j + \frac{n}{2}$. To a given solution $\bar{z}$, there corresponds an eigenvector of $t(u)$ and thus of $T$. The associated eigenvalue of the generalized Heun operator is

$$t_{\bar{z}} = q \left( j \cos(2\theta) + \frac{\mu}{2} - \frac{1}{2} \sum_{i=1}^{M} \left( \bar{z}_i + \frac{1}{\bar{z}_i} \right) \right) + q \left( \frac{nq}{2} - d \right) (nq - 2d + \mu + \nu). \quad (94)$$

### 5.3.3 The spectrum of the chopped correlation matrix

Once the generalized Heun operator is diagonalized, we can use the results to recover the spectrum of $C$, which shares with $T$ the same eigenspaces. Usually, this is done by acting with $C$ on the eigenvectors of $T$ and by reading out the eigenvalues of the chopped correlation matrix from the outcomes. Here, we present an alternative way to recover the spectrum. We construct a polynomial $P$ of order $M-1$ such that on $SV$

$$[C]_{K_{j,\ell,[v]}} = P([T]_{K_{j,\ell,[v]}}) = \sum_{i=0}^{M-1} a_i [T]_{K_{j,\ell,[v]}}^i. \quad (95)$$

If we have $P$ and an eigenvalue $t_{\bar{z}}$ of $[T]_{K_{j,\ell,[v]}}$, it is obvious that $P(t_{\bar{z}})$ gives an eigenvalue of $C$. To prove the existence of $P$, it is sufficient to know that the Heun operator and the chopped correlation matrix commute and to notice that $[T]_{K_{j,\ell,[v]}}$ is irreducible tridiagonal on $SV$. Therefore, the first $M$ powers of $T$ are linearly independent. The fact that $T$ is tridiagonal also allows to determine $P$. Indeed, since we have that

$$\beta_{i,j+m} \equiv \langle j, -j | [T]_{K_{j,\ell,[v]}}^i | j, m \rangle = 0 \quad \text{if } i < j + m, \quad (96)$$

we are led to

$$a_{M-1} = \frac{\langle j, -j | [C]_{K_{j,\ell,[v]}} | j, -j + M - 1 \rangle}{\beta_{M-1,M-1}} \quad (97)$$
and to the following recurrence relation:

$$a_i = \frac{1}{\beta_{i,i}} \left[ \langle j, -j | [C]_{K_{j,\ell,[v]}} | j, -j + i \rangle - \sum_{i'=i+1}^{M-1} a_{i'} \beta_{i',i} \right]. \quad (98)$$

Given this relation, constructing P is straightforward. The same then goes for the eigenvalues of C and the entanglement entropy.

# 6 Concluding remarks

We have considered systems of free fermions on Hamming graphs with $d+1$ parameters. We have recalled the connection between $\mathfrak{su}(2)$ and Krawtchouk polynomials. We have provided analytical expressions for the entanglement entropy of subsystems corresponding to Hamming graphs of a lower dimension. Similar results were also obtained for subsystems corresponding to a small number of neighborhoods. This was made possible by the decomposition of Hamming graphs in chains, i.e. the identification of the irreducible modules of the Terwilliger algebra of the scheme. Finally, we have shown how to construct block-tridiagonal operators commuting with the chopped correlation matrix that proved of great assistance in the determination of the entanglement entropy for subsystems consisting in a large number of neighborhoods.

The approach we used rested on the fact that the Hamming graphs are distance-regular or equivalently, related to a P- and Q- polynomial association scheme. Exploring the entanglement entropy of free fermions on graphs of other schemes would also be warranted. Another well known family of graphs are those of the Johnson scheme whose eigenvalue matrices are given in terms of the dual Hahn polynomials. We plan on examining the entanglement entropy of free fermions on these graphs and expect to report the results soon. In the future, it may also be interesting to consider the dual polar and Grassmann graphs [5], which are also distance-regular and related to q-polynomials of the Askey-scheme [31].

# Acknowledgements

We thank Krystal Guo for discussions. PAB holds a scholarship from the Natural Sciences and Engineering Research Council of Canada (NSERC). The research of LV is supported in part by a Discovery Grant from NSERC.

# References

[1] E. Bannai and T. Ito. *Algebraic Combinatorics I: Association Schemes*. Benjamin/Cummings, Menlo Park, 1984.

[2] P. Baseilhac, S. Belliard, and N. Crampé. FRT presentation of the Onsager algebras. *Letters in Mathematical Physics*, 108(10):2189-2212, Mar 2018.
[3] P.-A. Bernard, A. Chan, E. Loranger, C. Tamon, and L. Vinet. A graph with fractional revival. *Physics Letters A*, 382(5):259-264, Feb 2018.

[4] P.-A. Bernard, N. Crampé, D. S. Kabakibo, and L. Vinet. Heun operator of Lie type and the modified algebraic Bethe ansatz. arXiv:2011.11659, 2020.

[5] A. E. Brouwer, A. M. Cohen, and A. Neumaier. *Distance-Regular Graphs*. Springer-Verlag, Berlin, 1989.

[6] J. A. Carrasco, F. Finkel, A. González-López, and P. Tempesta. A duality principle for the multi-block entanglement entropy of free fermion systems. *Scientific Reports*, 7(1), Sep 2017.

[7] M. Christandl, N. Datta, T. C. Dorlas, A. Ekert, A. Kay, and A. J. Landahl. Perfect transfer of arbitrary states in quantum spin networks. *Physical Review A*, 71(3), Mar 2005.

[8] N. Crampé, K. Guo, and L. Vinet. Entanglement of Free Fermions on Hadamard graphs. *Nuclear Physics B*, 960:115176, 2020.

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

[10] N. Crampé, R. I. Nepomechie, and L. Vinet. Entanglement in Fermionic Chains and Bispectrality. *Roman Jackiw*, page 77-96, Jul 2020.

[11] N. Crampe. Algebraic Bethe Ansatz for the XXZ Gaudin Models with Generic Boundary. *Symmetry, Integrability and Geometry: Methods and Applications*, Dec 2017.

[12] B. Davies. Onsager's algebra and superintegrability. *Journal of Physics A: Mathematical and General*, 23(12):2245-2261, jun 1990.

[13] L. Dolan and M. Grady. Conserved charges from self-duality. *Phys. Rev. D*, 25:1587-1604, Mar 1982.

[14] C. El-Chaar. The Onsager Algebra. Master's thesis, University of Ottawa, 2010.

[15] J. T. Go. The Terwilliger Algebra of the Hypercube. *European Journal of Combinatorics*, 23(4):399-429, 2002.

[16] F. A. Grünbaum, L. Vinet, and A. Zhedanov. Algebraic Heun Operator and Band-Time Limiting. *Communications in Mathematical Physics*, 364(3):1041-1068, Jul 2018.

[17] M. A. Jafarizadeh, F. Eghbalifam, and S. Nami. Entanglement entropy in the Hypercube networks. 2014.
[18] M. A. Jafarizadeh, S. Nami, and F. Eghbalifam. Entanglement entropy in the Hamming networks. 2015.

[19] R. Koekoek and R. F. Swarttuuw. The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue. Technical Report 98-17, Delft University of Technology, Faculty of Information Technology and Systems, Department of Technical Mathematics and Informatics, 1998.

[20] H. J. Landau. *An Overview of Time and Frequency Limiting*, pages 201–220. Springer US, Boston, MA, 1985.

[21] J. I. Latorre and A. Riera. A short review on entanglement in quantum spin systems. *Journal of Physics A: Mathematical and Theoretical*, 42(50):504002, Dec 2009.

[22] F. Levstein, C. Maldonado, and D. Penazzi. The Terwilliger algebra of a Hamming scheme $H(d,q)$. *European Journal of Combinatorics*, 27(1):1–10, 2006.

[23] L. Onsager. Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition. *Phys. Rev.*, 65:117-149, Feb 1944.

[24] I. Peschel. Calculation of reduced density matrices from correlation functions. *Journal of Physics A: Mathematical and General*, 36(14):L205–L208, Mar 2003.

[25] I. Peschel. Special Review: Entanglement in Solvable Many-Particle Models. *Brazilian Journal of Physics*, 42(3-4):267-291, Mar 2012.

[26] I. Peschel and V. Eisler. Reduced density matrices and entanglement entropy in free lattice models. *Journal of Physics A: Mathematical and Theoretical*, 42(50):504003, Dec 2009.

[27] S.-s. Roan. *Onsager's algebra, loop algebra and chiral Potts model*. Max-Planck-Institut für Mathematik, 1991.

[28] D. Slepian. Some Comments on Fourier Analysis, Uncertainty and Modeling. *SIAM Review*, 25(3):379–393, 1983.

[29] P. Terwilliger. The Subconstituent Algebra of an Association Scheme (Part I). *J. Algebraic Comb.*, 1(4):363–388, Dec. 1992.

[30] P. Terwilliger. The Subconstituent Algebra of an Association Scheme (Part II). *J. Algebraic Comb.*, 2(1):73–103, Mar. 1993.

[31] P. Terwilliger. The Subconstituent Algebra of an Association Scheme (Part III). *J. Algebraic Comb.*, 2(2):177–210, June 1993.