Entanglement of inhomogeneous free fermions on hyperplane lattices

We introduce an inhomogeneous model of free fermions on a $(D-1)$-dimensional lattice with $D(D-1)/2$ continuous parameters that control the hopping strength between adjacent sites. We solve this model exactly, and find that the eigenfunctions are given by multidimensional generalizations of Krawtchouk polynomials. We construct a Heun operator that commutes with the chopped correlation matrix, and compute the entanglement entropy numerically for $D=2,3,4$, for a wide range of parameters. For $D=2$, we observe oscillations in the sub-leading contribution to the entanglement entropy, for which we conjecture an exact expression. For $D>2$, we find logarithmic violations of the area law for the entanglement entropy with nontrivial dependence on the parameters.


Introduction
This paper studies the entanglement entropy of free fermions on inhomogeneous lattices in various dimensions.
In quantum many-body physics, bipartite entanglement relates to how much a part denoted A of a system is correlated with its complementĀ. Quantifying this property is naturally of interest and warrants systematic investigation [1][2][3] as it plays in particular a key role in the characterization of critical points [4][5][6][7] and is of relevance in quantum information [8]. One such way of quantifying entanglement is through the entanglement entropy, provided by the von Neumann entropy S vN of the reduced density matrix ρ A of A. For a system in the pure state |Ψ , the entanglement entropy is defined by S vN = − Tr(ρ A log ρ A ), ρ A = TrĀ |Ψ Ψ|. (1.1) The entanglement entropy computation for an N -body system in the large-N limit is generally a challenging problem. Considerable attention has been focused on free-fermion models in one dimension, where this computation becomes tractable. Indeed, owing to the quadratic nature of such models, the spectrum of the reduced density matrix can be obtained solely from the truncated two-point correlation matrix C [9], which readily yields the entanglement entropy. For one-dimensional homogeneous chains, analytic results are available, see e.g. [10][11][12][13][14][15][16]. In particular, for the spin-1/2 XX spin chain with open boundary conditions, exact calculations based on the generalized Fisher-Hartwig conjecture [13] show that the entanglement entropy of a block of contiguous sites adjacent to a boundary scales as S vN = 1 6 log + . . . (1.2) in the limit of a large system embedded in an infinite chain. This is in agreement with results of conformal field theory (CFT) [7,17,18] with central charge c = 1.
Over recent years, a strong theoretical effort has also been aimed at better understanding entanglement properties of spatially-inhomogeneous systems [19][20][21][22][23][24][25][26]. While analytical methods used for homogeneous chains do not generalize to the inhomogeneous case, results have been obtained using CFT in curved background [20,21,26]. In Ref. [26], the authors use this formalism to argue that the entanglement entropy of a large class of inhomogeneous chains scales as in Eq. (1.2), where the inhomogeneous nature of the chains affects the sub-leading terms.
These one-dimensional models, both homogeneous and inhomogeneous, have a deep connection to orthogonal polynomials (more precisely, families of uni-variate orthogonal polynomials that belong to the Askey scheme [27]). An additional interesting feature of these models, made possible by the underlying bispectral context, is the existence of a tridiagonal Heun operator that commutes with C. The identification of this commuting operator exploits the parallel between the diagonalization of C (which is constructed out of operators projecting onto (i) the set of sites forming part A and (ii) states corresponding to energies that are present in the Fermi sea) with discrete-discrete problems in signal processing (that look at optimizing the concentration in time of a band-limited function), see [28].
Motivated by these results, we introduce here an inhomogeneous model of free fermions that hop on a (D − 1)-dimensional hyperplane sublattice of a D-dimensional hypercubic lattice. This model has D(D − 1)/2 continuous parameters that control the strength of hopping between nearest-neighbor sites, as well as D discrete parameters. For D = 2, this model reduces to the Krawtchouk chain [22]. We solve the general model exactly, and find that the eigenfunctions are given by multidimensional generalizations of Krawtchouk polynomials. These polynomials were introduced by Griffiths more than 50 years ago [29], rediscovered at the turn of the millennium and much studied thereafter [30][31][32][33][34][35][36][37][38]. For simplicity, we focus here on special cases (with one or more parameters set to zero) that were studied by Tratnik [39] 1 . Moreover, we construct a multidimensional generalization of the D = 2 Heun operator found in [22] which commutes with the chopped correlation matrix C.
We compute the entanglement entropy numerically for D = 2, 3, 4, for a wide range of different parameters. For D = 2, where the model has just one continuous parameter p ∈ (0, 1), we confirm recent results for p = 1/2 [26], and obtain new results for p = 1/2. In particular, we observe oscillations in the sub-leading contribution to the entanglement entropy, for which we conjecture an exact expression. For D > 2, we find logarithmic violations of the area law for the entanglement entropy with nontrivial dependence on the parameters. The paper will unfold as follows. The multidimensional inhomogeneous free-fermion model is introduced in Sec. 2, and various special cases are spelled out. Their analytic solutions are provided. A quick review of the characterization of the Griffiths polynomials is also given. The framework for the entanglement study is established in Sec. 3 by specifying the eigenstate in which the system will be taken, and by formulating the correlation matrix in that state. Upon fixing what has been called the part A of the system (the part whose entanglement we wish to concentrate on), the appropriate multidimensional generalization of the Heun operator is constructed and shown to commute with the chopped correlation matrix. This sets the stage for the determination of the entanglement entropy of various cases which is carried out in Sec. 4. Conclusions and perspectives are offered in Sec. 5.

The inhomogeneous free-fermion model
In this section we define an inhomogeneous free-fermion model on a hyperplane lattice, and derive its exact solution. We introduce the model in Sec. 2.1, and point out special cases for D > 2, namely the Tratnik and the one-parameter cases. We diagonalize the general Hamiltonian in Sec. 2.2 with the help of eigenfunctions that are investigated in more detail in Sec. 2.3. In particular, we obtain closed-form formulas for the various special cases under consideration in terms of Krawtchouk polynomials, which will be used for the computations in Sec. 4.

Definition of the model
We consider free fermions defined on the set of points which defines a (D − 1)-dimensional hyperplane. The number of such vertices is The parameter N gives the diameter of the hyperplane, and two sites x and y are nearest neighbors if there exist distinct indices i, j ∈ {1, . . . , D} such that As examples of the hyperplane lattices we consider, the case D = 2 corresponds to a linear lattice, and the case D = 3 corresponds to a planar triangular lattice, as illustrated in Fig. 1. In these figures, the vertices are denoted by dots, and nearest neighbors are connected by lines. The maximum possible number of nearest neighbors that a lattice site can have in these examples is 2 and 6, respectively, and in general is given by D(D − 1). In contrast, on a D-dimensional hypercubic lattice without the constraint (2.1), the maximum number of nearest neighbors is only 2D. For D = 4, the corresponding lattice is a tetrahedron obtained as the superposition of planar triangular lattices, but we do not represent it here.
The Hamiltonian is given by where the fermion annihilation and creation operators obey the usual anticommutation relations and R is an SO(D) matrix, namely In the Hamiltonian (2.5), for each value of i, one has the same types of hopping between nearest neighbors but with different coupling parameters. However, these parameters are related by virtue of being elements of a rotation matrix R. This matrix can be expressed as R = e B , where B is a real antisymmetric D × D matrix which has D(D − 1)/2 independent parameters. Finally, we restrict the real parameters α = (α 1 , . . . , α D ) to discrete values α i ∈ {0, 1}, which is required for the construction of the Heun operator in Sec. 3.2. The first term in the Hamiltonian (2.5) corresponds to nonuniform nearest-neighbor hopping, while the second term corresponds to a nonuniform chemical potential. In the following, we choose R such that (2.5) reduces to known inhomogeneous models in one dimension (for D = 2) and gives solvable generalizations thereof in higher dimensions.

The case D = 2
For D = 2, with x = (x, N − x), α = (1, 0), p ∈ R such that 0 < p < 1 and the Hamiltonian (2.5) becomes We recognize this model as one associated with the Krawtchouk polynomials, see e.g. [22]. For general values of D, we shall see that the Hamiltonian (2.5) is similarly associated with multivariate Krawtchouk polynomials.

Tratnik cases
The case D = 3 with R 12 = 0 is associated with the so-called bivariate Krawtchouk polynomials of Tratnik type [37][38][39]. For this case, the matrix R depends on 2 (rather than 3) independent parameters, which we denote by p 1 = R 2 13 and p 2 = R 2 23 : (2.10) Similarly, for D = 4, the Tratnik-type rotation matrix has 3 zero matrix elements, and 3 indepen-dent parameters that we denote by p 1 = R 2 14 , p 2 = R 2 24 and p 3 = R 2 34 : (2.11) In Fig. 2 we illustrate the lattice corresponding to the Tratnik case for D = 3 for different values of p 1 , p 2 and α. The color on each edge represents the magnitude of the hopping interaction between the two respective adjacent sites in the Hamiltonian (2.5). The actual values of the couplings are irrelevant, but darker colors represent stronger hopping interactions. The inhomogeneous nature of the model is clear. The color magnitude represents the strength of the couplings between neighboring sites, and darker colors correspond to stronger interactions.

One-parameter cases
The D = 3 rotation matrix with a single parameter p, is associated with a single uni-variate Krawtchouk polynomial. This matrix can be obtained from (2.10) by setting p 1 = 0 and p 2 = p. Similarly, the D = 4 rotation matrix (2.11) reduces for p 1 = p 2 = 0 and p 3 = p to 13) and is also associated with a single uni-variate Krawtchouk polynomial.

Solution
In order to solve this model, we begin by rewriting the Hamiltonian (2.5) in terms of a set of n(N, D) × n(N, D) matrices (H i ) x , x , where We now proceed to diagonalize the matrices (H i ) x , x . To this end, let us introduce, following [38], a set of D harmonic oscillator operators with commutation relations The algebra (2.16) has a representation on the orthonormal states defined by We define the operators X i and H i in terms of the harmonic oscillator operators by It is important to note that these operators are related by a unitary transformation, where U (R) is defined by [38] U (R) = exp and whose action on the harmonic oscillator operators is Moreover, for generic parameters R, the two sets of operators Xi and Hi generate together the whole su(D)-algebra realized by the operators a † j a k . Indeed, we have see also Lemma 2.2 in [33].
In view of (2.18), we see that and For a given N , we restrict ourselves to values of x that satisfy (2.1). The result (2.25) implies that the matrices (2.15) can be expressed as matrix elements of H i , From (2.21) and (2.24), we see that where the eigenvectors | k are given by 28) and the corresponding eigenvalues are given by We can also obtain the result , obtaining a result similar to (2.25), and then using H i = U (R) X i U (R) † and (2.28). We see from Eqs. (2.24), (2.25), (2.27) and (2.30) that the overlaps x| k are solutions to a bispectral problem defined from the set of equations obtained by equating in x|X i | k and x|H i | k the actions of X i and H i respectively on the bra and on the ket 3 .
In terms of the overlaps, the Hamiltonian (2.5) reads where we used (2.14), (2.26) and (2.27). We introduce Fourier-transformed fermionic operatorsc k , which satisfy the anticommutation relations Performing the sums over x and x in (2.31), we finally arrive at the diagonal form of the Hamiltonian, The eigenvectors of the Hamiltonian are given by where the vectors in F ⊆ V are pairwise distinct, and |0 is the vacuum state, The number of fermions in the state |Ψ (2.35) is given by |F|, the cardinality of the set F. We observe that the model's spectrum is integer valued, and does not depend on the matrix R. This can be understood from the fact that the Hamiltonian is constructed from operators H i that are related to X i by a unitary transformation (2.21). For the trivial rotation R = I (i.e., R ij = δ i,j ), we see that H i = X i , and the Hamiltonian (2.5) becomes diagonal in the | x basis, which leads to the spectrum in (2.37). Although the spectrum does not depend on R, the eigenfunctions x| k do depend on R, and are quite nontrivial, as described below. We conclude with two remarks regarding the spectrum. First, because the spectrum is integer valued, it appears to be gapped. However, because of the inhomogeneous nature of the model, one needs to investigate the gap in the scaling limit with respect to a rescaled system size, and for D = 2 the system is in fact gapless [25,26]. Second, we note that the spectrum is highly degenerated. For a given α, the components of k that correspond to vanishing components of α do not contribute to the energies and hence create the degeneracy. A refined understanding of the origin of the degeneracies and their physical implications is left for forthcoming investigations.

Eigenfunctions
We discuss here some further properties of the eigenfunctions k| x , corresponding to the eigenvectors of the H i operators (2.27). We have where the first equality follows from (2.28), and the second equality follows from results in [38]. The functions Q x ( k) are multivariate Krawtchouk polynomials, which obey the recurrence relation [38] (2.40) (While the notation of reference [38] will be adopted in the following, the reader might wish to consult the following papers [29][30][31][32][33][34] for more on these polynomials.) The polynomials Q x ( k) can be obtained from a generating function [38] where z D ≡ 1. Moreover, the function W ( x, k) in (2.39), which is given by is the discrete measure appearing in the orthogonality relation [38], Let us perform a consistency check by using (2.39) to rewrite the recurrence relation (2.40) in terms of k| x , Using (2.42), we observe the following simplification, Hence, the recurrence relation (2.44) takes the form which coincides with the result obtained by applying k| to (2.25) and then using (2.27). The functions Q x ( k) also satisfy the difference equation [38] x Similarly as before, we can use (2.39) to rewrite the difference equation in terms of x| k , and obtain which coincides with the result obtained by applying x| to (2.30) and then using (2.24).

The case D = 2
For the case D = 2 with the matrix R given by (2.8), the functions Q x ( k) can be expressed in terms of Krawtchouk polynomials, and the Pochhammer (or shifted factorial) symbol (a) k is defined by

Tratnik cases
For the case D = 3 with R 12 = 0 in Eq. (2.10), Q x ( k) can be expressed as a product of two uni-variate Krawtchouk polynomials [38] (see also [35] and the appendix of [37]), where k x (k; p; N ) is given by (2.50).
For the case D = 4 with the rotation matrix given by (2.11), the functions Q x ( k) can be expressed as a product of three uni-variate Krawtchouk polynomials [35]

One-parameter cases
For the D = 3 case with one parameter (2.12), the normalized eigenfunctions are given in terms of a single uni-variate Krawtchouk polynomial [38], Similarly, for the D = 4 case with one parameter (2.13), the normalized eigenfunctions are given by To avoid indeterminate values when evaluating these functions numerically, it is convenient to define the uni-variate Krawtchouk polynomials kx(k; p; N ) (2.50) to be 0 if N < min(x, k).

A derivation of the overlap coefficients
The explicit results for overlap coefficients given above, which will be needed for the entanglement entropy computations in Sec. 4, were borrowed from the literature. For completeness, we sketch here a way of deriving such results, see also [30,33].
Denote by v 1 , . . . , v D the standard orthonormal basis vectors of C D and (E ij ) i,j=1,...,D the elemen- (2.56) We will use the following notation for the natural orthonormal basis of (C D ) ⊗N : The basis vectors | x in bijection with the lattice sites are realized as an orthonormal basis of the symmetric part in the tensor product (C D ) ⊗N . More precisely, consider the N -th symmetrizer, that is the normalized sum over all permutations in S N of the factors in the tensor product: Note that P commutes with E ij . Then we set The prefactor ensures that the vectors | x are orthonormal. Note that | x is also obtained by applying P , with the same prefactor, on v (i 1 ,...,i N ) whenever (i 1 , . . . , i N ) contains x 1 times 1, x 2 times 2, and so on. It is easy to check that the action of E ij on these vectors is where a † i , a j are as in (2.18). Now recall that R ∈ SO(D) and define U (R) on (C D ) ⊗N by R ⊗ · · · ⊗ R, that is N v (a 1 ,...,a N ) . (2.61) The operators P and U (R) commute, so that U (R) restricts to the symmetrized product. Therefore, U (R) acts on | x . The operators X i and H i acting on | x are defined by in agreement with (2.19). To obtain the eigenvectors | k of H i , we need to apply U (R) † to the eigenbasis | x of X i , and for this we are going to use (2.59), the fact that U (R) commutes with P and (2.61). D = 2. Let us detail the calculation for D = 2. We have, for k 1 + k 2 = N , For the last equality, we collect all terms in the sum with x 1 occurrences of 1 in (a 1 , . . . , a k 1 , b 1 , . . . , b k 2 ). We can choose i of them among (a 1 , . . . , a k 1 ), resulting in R i 11 , and ( We recover (2.39) for D = 2 with (2.49), see e.g. [43].
Arbitrary D. Define the binomial, and more generally multinomial coefficients, as We can write the formulas for the overlap coefficients for D = 2 in (2.63) as √ where the sum is over nonnegative integers a 11 , a 12 , a 21 , a 22 satisfying: (the sums of the columns are respectively k 1 and k 2 ).
Of course, here there is only one independent index a 11 , which can go from 0 to x 1 . Reproducing straightforwardly the reasoning made for D = 2, one obtains the following result for general D, where the sum is now over nonnegative integers (a ij ) i,j=1,...,D such that the sum of the r-th line gives x r and the sum of the r-th column is k r .

Correlation matrix and algebraic Heun operator
In this section, we define the correlation matrix and express its entries in terms of overlaps, which are given in Sec. 2.3. We then construct, for any D, an operator that commutes with the chopped correlation matrix, thereby generalizing the D = 2 tridiagonal Heun operator found in [22,23].

Correlation matrix
Let us consider a fermionic eigenstate |Ψ as in (2.35), where the set F ⊆ V is given by The correlation matrix in this state, which will be needed later for the entanglement entropy computation, is given by the matrix elementŝ We find using (2.33) and (2.35) that The correlation matrix elements (3.2) are therefore given bŷ which implies that the correlation matrix iŝ

Algebraic Heun operator
Recall that the set F is defined in (3.1), and let us now define a similar set A ⊆ V by where β = (β 1 , . . . , β D ) with β i ∈ {0, 1}. If D = 2 and β is a unit vector, A = {0, . . . , L} is a segment of length L + 1. For D = 3 and if β is a unit vector, subsystem A is the triangular lattice from which a triangular part has been removed, see Fig. 3. The projection operators on subsets A and F are respectively. We consider a general symmetric bilinear expression of the bispectral operators X i and H i introduced in Sec. 2.2, This provides an extension to higher dimensional bispectral problems of the algebraic Heun operator introduced in [28], see also [14,16,22,23]. We solve for the coefficients µ i , ν i , ρ ij such that T commutes with both π 1 and π 2 . Let us first consider the commutativity of T with π 1 . We observe that x ∈ A iff i∈b x i ≤ L, for some set b ⊆ {1, . . . , D} corresponding to the choice of β. Furthermore, [T , π 1 ] = 0 iff x |T | x = 0 for all x and x such that (3.9) The fact that T is linear in H i implies that T is linear in E jk := a † j a k with j = k, which has the property (3.10) Let us consider x and x as in (3.10). Then We choose It then follows from (3.11) that which vanishes for ν i = −(2L + 1). For the commutativity of T with π 2 , we repeat the calculation in the | k basis. In particular, k ∈ F iff i∈f k i ≤ K, where f ⊆ {1, . . . , D} corresponds to the choice of α. The fact that T is linear in X i implies that k |T | k = 0 unless i∈f k i = K , i∈f k i = K + 1 . (3.13) Setting we obtain µ j = −(2K + 1).
In conclusion, the operator T in (3.8) commutes with both projectors, for the coefficient values where C is the chopped correlation matrix C = π 1 π 2 π 1 . (3.18) We have checked numerically for various examples that the eigenvalues (apart from 0) of the chopped Heun operator π 1 T π 1 = π 1 T are not degenerate. For D = 2, the result in [22,23] for the Heun operator is recovered. While the eigenvalues of C are expected to agglomerate near 0 and 1, the spectrum of T is usually well spaced and thus free of such problem. Since both matrices share a common eigenbasis, the Heun operator T offers an interesting tool to diagonalize the chopped correlation matrix C. It also opens the door to the application of Bethe ansatz methods, which was done for the case D = 2 [44]. For these reasons, we highlight Eqs. (3.8) and (3.16) as a main result of this paper.

Entanglement entropy
In this section, we investigate the entanglement entropy of certain eigenstates |Ψ of the Hamiltonian (2.5). Because the Hamiltonian is quadratic in terms of the fermionic operators, the reduced density matrix associated to a region A is a Gaussian operator whose eigenvalues can be obtained [9,45,46] from the chopped correlation matrix C in (3.18). Accordingly, the entanglement entropy (1.1) reads S vN = − Tr(C log C + (I − C) log(I − C)). (4.1) In studying the entanglement entropy as a function of N for given values of D, α, and β, it is also necessary to specify K and L which define the sets F in (3.1) and A in (3.6), respectively.
Let us note that if α = (1, 1, . . . , 1), then the entanglement entropy trivially vanishes. Indeed, in that case we have α · k = D i=1 k i = N for all k ∈ V . From the definition (3.1) of F, it follows that for this case |F| = 0 if K < N and |F| = n(N, D) if K = N , where n(N, D) is given by (2.2). These situations correspond to a filling fraction |F|/n(N, D) equal to 0 and 1, respectively, and in both cases the system is in a product state with zero entanglement entropy.
A natural choice is to consider K such that the system is at half filling, namely that the number of elements in F is |F| = 1 2 n (N, D). Moreover, for simplicity, we always choose α such that α = i for some value of i (see Eq. (2.4)). In that case, the choice of K does not depend on i.
For D = 2 (the one-dimensional chain), the notion of half filling is straightforward, and we choose K = N/2. Looking at the definition of F in (3.1), this implies that |F| = N/2 + 1 (remember that the chain has N + 1 sites for D = 2). Hence, we have a perfect half filling for odd values of N , whereas for even N the filling is N +2 2(N +1) and tends to half filling in the limit N → ∞. For D > 2 however, it is not always possible to choose a value of K such that we are perfectly at half filling. In these cases, we choose the largest K such that |F| is less than or equal to 1 2 n(N, D). For D = 3, 4, we find that the value of K that satisfies this property is K = N/D − 1. The formula becomes more involved for larger values of D, but we do not consider these cases.
Similarly, we take β = j for some value of j, and systematically set L = K, for simplicity.

The case D = 2
For D = 2, the system is a one-dimensional chain of length N + 1, and we fix α = β = 1 . We note that the entropy does not change if we choose α = β = 2 instead. The entropy also depends on the parameter p (see Eq. (2.9)), and if instead of α = β we had α = β, this would be equivalent to substituting p with 1 − p in our results. Since we investigate the dependence of the entanglement entropy on p, our choice for α and β is thus generic. We compute the entanglement entropy for L = K = N/2 as a function of N for various values of p. Numerically, we use the eigenfunctions given in (2.49) and (2.50) to compute the chopped correlation matrix and use Eq. (4.1). We conjecture the following result for the entanglement entropy, a(p) is a non-universal constant with respect to N , and the ellipsis indicate terms of order smaller than N −1 in the large-N limit. We compare our conjecture with numerical results in the left panel of Fig. 4 and find excellent agreement, already for moderate values of N . We computed the entanglement entropy for numerous values of p, and systematically found a similar match between the numerical data and the conjecture of Eq. (4.2). However, we did not include all the corresponding graphs in Fig. 4 for clarity. The leading term 1 6 log N +1 2 in Eq. (4.2) is in agreement with previous results from Ref. [26]. In that paper, the authors used methods from CFT in a curved background [20] to extract the leading contribution for the special case p = 1/2 at half filling. In particular, the system is described by a CFT of massless Dirac fermions in a curved (1 + 1)-dimensional background. We find that this leading contribution holds for any values of p at half filling.
The constant term with respect to N in Eq. (4.2), that we denote a(p), depends non-trivially on the parameter p. For p = 1/2, we find a(1/2) = 0.4786, which corresponds to the exact results for the homogeneous XX chain [13] at half filling. This was already observed in Ref. [26]. We report our results for a(p) in the right panel of Fig. 4. The dots are obtained from numerical fits of our results for the entropies, and the solid line is a guide to the eye. The latter does not reflect an analytical prediction or a conjecture, and therefore we do not report its expression here.
Finally, we observe oscillations in the sub-leading contribution of order N −1 in Eq. (4.2). Similar oscillations have already been studied extensively for homogeneous chains [12,13,[47][48][49][50] and in the context of the so-called Rainbow chain [21]. For homogeneous chains, the frequency of the oscillations depends only on the filling fraction. Here, for a fixed filling fraction (we are always at half filling), we observe that the frequency of the oscillations depends on the model through the parameter p and the function m(p) in Eq. (4.3). As a particular case, for p = 1/2 we have m(1/2) = 1 and we recover the results of Ref. [26]. The non-trivial dependence of the frequencies of sub-leading oscillations on the parameters of the model is one of our main results, and it calls for further investigations at different filling fractions and for different inhomogeneous chains. These issues will be addressed in a forthcoming publication.

The case D = 3
In this subsection, we investigate the entanglement entropy for D = 3, where the whole system is a planar triangle, as depicted on the right of Fig. 1. As discussed as the beginning of Sec. 4, we systematically take L = K = N/3 − 1.

Numerical results for Tratnik with D = 3
We start with the Tratnik case D = 3 and R 12 = 0, where the eigenfunctions x| k have a simple closed form, see Eqs. (2.39), (2.42) and (2.52). The entropy depends on p 1 , p 2 , as well as the system size, the eigenstate under consideration and the subsystem, through the respective choices of N , K and α, and L and β. We compute the entanglement entropy numerically for various values of p 1 and p 2 with N = 3n + 1, L = K = N/3 − 1, α = 2 , and β = 1 . We display our results on the top panels of Fig. 5. The dots are the numerical data obtained by direct diagonalization of the chopped correlation matrix, and the solid lines correspond to fits of the form where the ellipsis indicates terms that are sub-leading compared to N log N in the large-N limit. For critical free fermions in d dimensions, the entanglement entropy exhibits a logarithmic violation of the area law [2] and scales as S vN ∼ N d−1 log N , where N is a characteristic length scale of the subsystem along one spatial dimension [51][52][53]. In our case, the fact that the leading contribution to the entanglement entropy in Eq. (4.4) scales as N log N thus indicates that, similarly to the onedimensional case discussed in Sec. 4.1, our inhomogeneous system has a similar behavior as critical free fermions in two dimensions. However, here the leading coefficient depends on the parameters p 1 and p 2 in a non-trivial way, and we restrict ourselves to numerical evaluation of γ. In the bottom left panel of Fig. 5, we show the dependence of this coefficient on p 2 for fixed values of p 1 . The dots are obtained from our numerical data, and the solid lines are fits that serve as guides to the eyes. We do not report their expressions.
Another difference with the one-dimensional chain is that here we do not have parameter-dependent oscillations. In the bottom right panel of Fig. 5, we plot the entanglement entropy with the same parameters as in the top panels, but for all values of N (instead of focusing on N = 3n + 1). We observe plateaus where the entropy has almost the same value for consecutive values of N . This phenomenon is very different in nature from the oscillations we observe in the one-dimensional chain, since here the plateaus do not depend on p 1 , p 2 , but instead appear because of our choices for L and K. In fact, we verified that for L = K = N/M − 1 (hence away from half filling for M = 3), we have plateaus of length M . Here they have length 3 because L = K = N/3 − 1. These plateaus are thus artifacts of the geometry of the model, and we remove them by considering N = 3n + 1.  For the D = 3 Tratnik case with parameter values p 1 = 1/2, p 2 = 1/4, the nonzero matrix elements of the rotation matrix (2.10) are ±1/2 , ±1/ √ 2, and both Krawtchouk polynomials in (2.52) have the same parameter value (namely, 1/2). For this case, we find numerically that the entanglement entropies corresponding to various pairs ( α, β) are equal. For example, all the pairs {( 1 , 1 ) , ( 2 , 2 ) , ( 3 , 2 ) , ( 1 , 3 )} correspond to the same entanglement entropy, and similarly for {( 2 , 1 ) , Many of the above equivalences in entanglement entropy can be explained by the existence of simple changes of basis for this special case. Indeed, consider the transformation P (1) defined by whose action on the k-basis is given by Conjugation by P (1) therefore gives Similarly, conjugation by P (1) on the projectors gives where π X i and π H i denote the projectors π 1 with β = i , and π 2 with α = i , respectively. Consequently, the chopped correlation matrices π X 1 π H 1 π X 1 and π X 3 π H 1 π X 3 are related by conjugation, which explains the fact that ( 1 , 1 ) and ( 1 , 3 ) have the same entanglement entropy. With this argument, we find the following equivalences: Similarly, the transformation P (2) defined by has the following action on the k-basis Therefore P (2) : 12) and similarly for the corresponding projectors. It follows that the last two lines of (4.9) also have the same entanglement entropy. In this case, hopping is not allowed between the subsystem and its complement.

Decoupled Tratnik case
We observe that for the D = 3 Tratnik case (2.10) with R 12 = 0, the entanglement entropy exactly vanishes for a specific set of vectors ( α , β), namely for all allowed values of K, L, N . This is consistent with the fact that, for this case, the subsystem does not interact with its complement, as illustrated in Fig. 6. More generally, it follows from the definitions (2.19) that (4.14) This gives 0 if R il = 0. We consider the following particular situation. Given a matrix R, we choose α and β such that the only pairs (α i , β l ) equal to (1, 1) correspond to coefficients R il = 0, that is, Then we deduce that in this case: α · H , β · X = 0 . (4.16) Now recall that π 1 is the projector on the sum of certain eigenspaces of β · X, and as such it is a polynomial in β · X. Similarly, π 2 is a polynomial in α · H. From (4.16), we have that π 1 and π 2 commute. Therefore they can be simultaneously diagonalized (with eigenvalues either 0 or 1) and thus C = π 1 π 2 π 1 = π 1 π 2 has eigenvalues which are either 0 or 1. Therefore, the entanglement entropy vanishes in the particular situation (4.15).

One-parameter case
Let us now consider the one-parameter case for D = 3, with Hamiltonian (2.12) whose eigenfunctions are given by (2.54), and let us set α = 2 . (The case α = 3 is similar, while here α = 1 is trivial.) For β = 1 , the system becomes decoupled and has 0 entanglement entropy, as in Sec. 4.2.3. However, for β = α = 2 , the entanglement entropy is nonzero. As for the Tratnik case discussed in Sec. 4.2.1, we observe plateaus if we consider all values of N . We impose N = 3n + 1, and again observe that the entropy exhibits a logarithmic violation of the area law, see the left panel of Fig. 7. The leading coefficient γ from the expansion (4.4) is displayed on the right panel of Fig. 7. Similarly to the Tratnik case, we restrict ourselves to a numerical evaluation of this coefficient and the solid lines are simple guides to the eye. A notable difference with the Tratnik case discussed of Sec. 4.2.1 is that γ appears to depend only weakly on p. Most of the curves in the left panel of Fig. 7 are close to each other, as compared to the ones in the top panels of Fig. 5. In particular, the curves for p = 1/3, 1/4 and 1/5 are almost indistinguishable. A potential explanation for this behavior is the following. In the case of homogeneous free fermions on the square lattice, the scaling S vN ∼ N log N stems from the fact that the Hamiltonian in two dimensions can be expressed as a sum of Hamiltonians for one-dimensional chains via dimensional reduction [49,54]. The entanglement entropy of a two-dimensional strip is exactly given as a sum of N entropies of one-dimensional chains at different filling fraction, and each of them scales as log N . The dimensional reduction thus allows one to not only understand the logarithmic violation of the area law, but also to compute the coefficient of the leading term, as well as the sub-leading ones. In our model, considering the one-parameter case with α = 2 implies that the Hamiltonian (2.5) only contains hopping interactions along lines parallel to 2 − 3 , and hence can be seen as a sum of one-dimensional chains. The fact that the leading term of the entanglement entropy in one dimension is independent of p may help explain the weak dependence on p here. However, because of the geometry of the problem, the computation of the leading coefficient γ from one-dimensional chains appears more involved than on the square lattice, and we have not investigated that issue further.

The case D = 4
In this section, we investigate the entanglement entropy for D = 4. Based on our results for D = 2, 3, we expect the entropy to exhibits a logarithmic violation of the area law, and to scale as We consider two particular cases, (i) the Tratnik case, see Eq. (2.11), for which the eigenfunctions are given by (2.53), and (ii) the one-parameter case discussed in Sec. 2.3.3. Similarly as for the case D = 3, we choose L = K = N/4 − 1. Without restrictions on the allowed values for N , we also observe plateau-like effects, that are only related to the geometry of the system and our choices of L and K. Hence, we restrict ourselves to N = 4n + 1. We display our numerical results in the left and right panels of Fig. 8 for cases (i) and (ii), respectively. The dots are obtained by direct diagonalization of the chopped correlation matrix, whereas the solid lines correspond to fits of the form of Eq. (4.17).
For the one-parameter case, the curves for p = 1/3, 1/4, 1/5 and 1/6 are almost indistinguishable, and correspond to a leading term of γ 0.06. Our heuristic explanation for this phenomenon is similar as the one we discuss in Sec. 4.2.4 for the one-parameter case for D = 3.

Conclusion and outlook
In this paper, we introduced and solved a new inhomogeneous model of free fermions in arbitrary spatial dimensions. The model coincides with the Krawtchouk chain [22] in the one-dimensional case, and its solution in arbitrary dimensions relies on multivariate Krawtchouk polynomials. Moreover, using bispectral properties of the problem, we constructed an operator that commutes with the chopped correlation matrix and provides a natural extension of the Heun operator obtained in [22] to higher dimensions.
We investigated the entanglement properties of half-filled eigenstates of the Hamiltonian (2.5) for D = 2, 3, 4, which correspond to systems in one, two and three spatial dimensions, respectively. We found that the entanglement entropy scales as S vN = γ N D−2 log N + . . . , (5.1) and thus exhibits a logarithmic violation of the area law. For D = 2, the leading coefficient is γ = 1/6, in agreement with previous results obtained in the framework of curved-space CFT [26]. We also observed and conjectured the exact form of sub-leading oscillations, see Eq. (4.3), but an analytical proof is beyond the scope of this paper. Crucially, the frequencies of those oscillations depend on the model through the parameter p. This property is in stark contrast with known sub-leading oscillations in homogeneous systems, which depend only on the filling fraction, and represents one of the main results of this paper. For D = 3, 4, the leading coefficient γ depends on the parameters of the model, and this dependence is much weaker for the one-parameter cases. We restricted ourselves to numerical evaluations of this coefficient.
There are many open problems that are worth investigating in the future. A first one would be to provide an analytical proof for the conjecture of Eqs. (4.2), (4.3), and hence extend the results of Ref. [13] to inhomogeneous systems. However, we expect this problem to be very challenging, since the inhomogeneity of the model breaks the Toeplitz and Hankel structure of the correlation matrix. Another complementary idea is to use curved-space CFT methods to understand the sub-leading oscillations, similarly to the Rainbow chain [21], as well as the leading terms for p = 1/2. In parallel to analytical results, we believe that the model-dependent sub-leading oscillations deserve further investigations in the context of other inhomogeneous chains. For instance, preliminary results suggest that a similar phenomenon occurs for the anti-Krawtchouk chain [55]. It would also be interesting to better understand the dependence of γ on the parameters of the model for D = 3, 4, either via analytical computations or with field-theoretical arguments in higher dimensions. Another intriguing question is whether one could use the Heun operator to compute entanglement entropies faster and with higher precision than via direct diagonalization of the chopped correlation matrix. Preliminary results seem to suggest so. Moreover, the Heun operator could also be used to derive analytical results for the entanglement entropy via Bethe ansatz techniques [44]. Other natural extensions of our work include the investigation of inhomogeneous models based on other representations of su(D) and on graphs of P-and Q-polynomial association schemes [56,57]. The entanglement entropy of free fermions on the celebrated Hamming scheme was found to be related to the entanglement in Krawtchouk chains [58], and similarly the entanglement of free fermions on ordered Hamming graphs [41] is expected to be related to the entanglement in the higher-dimensional models introduced in this paper.