Entanglement of Exact Excited Eigenstates of the Hubbard Model in Arbitrary Dimension

We compute exactly the von Neumann entanglement entropy of the eta-pairing states - a large set of exact excited eigenstates of the Hubbard Hamiltonian. For the singlet eta-pairing states the entropy scales with the logarithm of the spatial dimension of the (smaller) partition. For the eta-pairing states with finite spin magnetization density, the leading term can scale as the volume or as the area-times-log, depending on the momentum space occupation of the Fermions with flipped spins. We also compute the corrections to the leading scaling. In order to study the eigenstate thermalization hypothesis (ETH), we also compute the entanglement Renyi entropies of such states and compare them with the corresponding entropies of thermal density matrix in various ensembles. Such states, which we find violate strong ETH, may provide a useful platform for a detailed study of the time-dependence of the onset of thermalization due to perturbations which violate the total pseudospin conservation.

The question of how equilibration and thermalization arise in isolated quantum (many-body) systems led to the eigenstate thermalization hypothesis (ETH) [1][2][3][4]. ETH states that in the thermodynamic limit, the eigenstate expectation value of a few-body operator in a typical eigenstate of a many-body Hamiltonian at energy E is equal to the microcanonical average at the mean energy per volume E/V . The two main interpretations of the ETH, the weak versus strong ETH, state that almost all versus all the finite energy density eigenstates of a many-body Hamiltonian appear thermal to all local measurements [5].
The ETH also has fundamental implications on quantum information-inspired quantities that characterize the excited states. More specifically, the entanglement spectrum and the resulting entanglement entropy have long become powerful diagnostics of topological order, gapless or gapped nature of ground-states, and other properties [6]. One implication of the ETH is that thermal states have volume law entanglement as opposed to area-type entanglement entropy of the ground state and low-lying excited states of the system. The volume law entanglement is then thought to return to an area law entanglement when/if the many-body localization sets in [4].
Unfortunately, the paucity of exact results makes it difficult to test or demonstrate ETH and its consequences in generic, non-integrable, many-body models in more than one space dimension with realistic electron-electron interactions. Numerical studies are limited to the very small system sizes imposed by the exact diagonalization. Motivated by the fact that a class of exact excited eigenstates of the Hubbard model is known [7,8], that the number of such states is a exponentially large in volume [9], and that their energy density differs from the ground state energy density by a finite amount, here we obtain the closed form exact expressions for the entanglement spectrum, the von Neumann entanglement and Renyi entropies of such states. The entanglement entropy for these states shows either a ln(V ) law, or a V (volume) law, or even an area-times-log law, depending on the number and the momentum space distribution of the flipped spins in the state. When their entropy is sub-extensive, such states therefore clearly violate strong ETH. Even when the entropy scales with V , the prefactor is independent of the Hubbard U , and is not expected to correspond to the entropy in the microcanonical average, which should be a non-trivial function of U . Despite being in the middle of the full Hubbard spectrum, the pure spin singlet eta-pairing states, which show ln(V ) entanglement, are simultaneously the ground-states and the most excited states in their specific quantum number sectors. Kantian et.al. proposed an interesting way to prepare the eta-pairing state with cold atoms in optical lattice [10]. If successfully implemented, our results make a concrete prediction about the reduced density matrix of a small subsystem, and the precise way that the remainder serves as a thermal bath.
The eta-pairing states with flipped spins are richer. They display either volume or area-times-log entanglement, depending on the momentum space occupation of the flipped Fermions. We find that even for the states with volume law entanglement, the entanglement Renyi entropies do not match those of the thermal density matrix in the canonical ensemble. They match the Renyi entropy of the thermal density matrix in a grand canonical ensemble, but with additional constraints on the quan-tum numbers of the states.
We consider the Fermionic Hubbard model on a hypercubic lattice in any dimension. The Hamiltonian iŝ H =T +V wherê c † rσ is the Fermionic creation operator at a site r (belonging to the hypercubic lattice) and spin projection σ =↑ or ↓. The total number of sites is M and the first sum is over the nearest neighbor links. The exact, 2N -particle, spin-singlet, normalized, etapairing eigenstate [7] ofĤ that we firstly focus on is where M !N ! , and π = (π, π, . . . , π). This follows readily from the commutator ofĤ and [7,8].
As shown by C.N. Yang [7], this state is not the ground state of the Hubbard model for either U ≶ 0. At half filing, µ = U/2 and the energy of this state vanishes. For repulsive U , the ground state at half filling is an antiferromagnetic insulator [11] with negative energy per particle (see e.g. [12]). For attractive U the ground states are an s-wave superconductor and a charge density wave [13], also with negative energy per particle. At weak coupling (U t) and near half filing, |ψ N sits near the middle of the energy spectrum. That is because the weakly perturbed filled Fermi sea with momenta centered near k = 0 is near the bottom of the many-body band and with momenta centered near k = π is near its top. In the Supplementary Information we introduce a generalization of the Hubbard model Eq [2] for which there exist similar eta-pairing eigenstates.
We partition the M sites into a group A with M A sites and a group B with M B = M − M A sites and compute the reduced density matrixρ A by tracing all the degrees of freedom in the group B. We then take the thermodynamic limit N → ∞, M → ∞ such that the boson filling N/M → ν ∼ O(1). After this limit, we then take M A 1. The system A is therefore small compared to B so that B can serve as its bath, but still large enough to allow scaling of its entanglement entropy.
We now sketch the derivation of the reduced density matrix [14]. We use integration over the contour C encircling the origin in the complex z-plane counterclockwise to re-write the eta-pairing state as: Terms in the sum r e iπ·rĉ † r↓ĉ † r↑ commute, therefore we can write the exponential of the sum as the product of the exponentials. Moreover since (ĉ † r↓ĉ † r↑ ) 2 = 0 we see that the operator part of Eq. 4 becomes r 1 + ze iπ·rĉ † r↓ĉ † r↑ |0 . We then obtain : |0 A denotes the state with all sites in the region A empty.
Only the same powers of z 1 and z 2 survive the contour integration. Expanding (1+z 1 z 2 ) M B using binomial expansion, performing the contour integration and eliminating the sum coming from the binomial expansion, gives the entanglement spectrum: We assumed M A < N . The states |k are orthonormal eta-pairing states of the A side: A Vandermonde convolution confirms that M A k=0 λ k = 1. Similar result for a ferromagnetic Heisenberg model appears in Ref. [15]. Ref. [16] also studies the η-pairing state, but uses a different normalization; an expression for λ k in which N appears only via N/M is quoted in [17]. Eq. 6 shows that for each k, the eigenvalue of the density matrix is equal to the number of ways to simultaneously place k pairs on M A sites and N − k pairs on M B = M − M A sites, divided by the number of ways to place N pairs onto M sites. The system is subject to the constraint on no double pair occupancy. In the thermodynamic limit, the largest number of configurations corresponds to the uniform particle density, i.e. λ k should be very sharply peaked about k m = M A (N/M ).
In the limit of interest, we can use the Stirling formula n! ≈ √ 2πne n(ln n−1) where n is large. Then, where κ = ν(1 − ν)M A . This form is valid as long as ν is not infinitesimally close to 0 or 1. Substituting the above Gaussian form into the von Neumann entanglement entropy S A = − M A k=0 λ k ln λ k , and replacing the discrete sum over k with an integral we obtain that S A scales as the logarithm [16,17] of the number of sites in the region A: The small value of S A seems to be in a contradiction with the ETH motivated expectation that finite energy density excited states in the middle of the many-body spectrum should thermalize with the entanglement entropy scaling as the volume ( ∼ M A ). However, it is not, due to the existence of additional pseudospin symmetry operators [8], and the eta-pairing states are the only states in their symmetry sector. The existence of a global conserved pseudospin [7,8] is special to the Hubbard model, and the corresponding operator is [7,8] . BecauseĴ 2 commutes withĴ + , the state |ψ N corresponds to the maximal eigenvalue of theĴ 2 , namely The set F consists of any of the wavevectors in the 1 st Brillouin zone. We denote the number of k's in F by N k . We normalize these states by computingĉ k N = [9]. Although this is a very large number, the total number of states in the Hilbert space is larger i.e. 4 M . Thus the relative fraction of eta-pairing states vanishes as M → ∞ [9]. The eigenenergy of |ψ where k are the energies of the kinetic term (i.e. in two dimensions k = 2t(cos k x + cos k y )). Consider first the states in Eq. 11 with N = 0; all such states can be easily constructed, as they are non-interacting. For a given N k , such fully spin polarized states are highest weight spin states S z = S = N k 2 . They also have J 0 = −J = and γ m and φ m (k) are respectively the eigenvalues and orthonormal eigenvectors of the Hermitian M r∈A e i(k−k )·r , with k and k ∈ F. The Fermion operators in Eqs. 14-15 obey â † m ,â m = b † m ,b m = δ m,m , and, because they live in different regions in real space, â † m ,b m = 0. Using Eq. 13, we can write where the sum is over all the different 2 N k ways to partition the N k "orbitals" m into those occupied by a † 's, denoted by the set {m A }, and the complementary set {m B } occupied byb † 's. If, for any given partition, there are N A "orbitals" in the set {m A }, then there are N k − N A "orbitals" in the set {m B }. Here, The states in Eq. 16 are The reduced density matrix can again be calculated with the help of the contour integral representation of |ψ where The von Neumann entanglement entropy is then Using the Vandermonde convolution, we can write the above as is the von Neumann entropy of the free Fermi gas. As computed in Fig. 1 (for the two dimensional case) and as discussed by Lai and Yang [19], it can result in either the volume or the area-times-log "law", depending on which k-points are occupied.
To analyze the second term in Eq. 21, we note that the sum over the partitions is sharply peaked around N A ≈ N k M M A , which is large in the limit of interest. Therefore, λ k j is peaked about a large value of j and we can use Stirling's approximation. Again, replacing the discrete sum by an integral, we finally find where ν k = N k /M . We see that the leading order scaling comes from the free Fermion part (i.e.S k A ), and the correction scales with the logarithm of the number of sites in the region A. Logarithmically diverging subleading contribution was also argued for in Ref. [20] , but with a different prefactor.
The Renyi entropy, S k,(n) A = 1 1−n ln (Tr Aρ n A ), can be computed for the states Eq. 11 using similar techniques. In the same limit as before, we find (24) Note that ν (1) The formulas Eq. 22 and Eq. 23 are therefore identical as n → 1. Again, the leading scaling comes from the free Fermi part and the correction scales as ∼ ln M A . In the context of the ETH, it is interesting to ask whether the entanglement entropy density -be it von Neumann or Renyi -for the above mentioned exact eigenstates of the Hubbard model match the entropy density for the thermal density matrixρ th = e −β(Ĥ−µ (1) withĤ being the Hubbard Hamiltonian. We included µ (1,2) th to separately control the average value of N k and N . If the trace ofρ th is to be performed over all the states in the Hilbert space of the Hubbard model, then they should not match, because Trρ th should depend on the interaction U while ρ {k} A is U -independent. However, if the trace is restricted to states of the type Eq. 11, and the distribution of the occupied k states results in the "volume" law (see Fig.1), then the first (leading) term in Eq. 23, indeed matches the "thermal" Renyi entropies computed in the grand canonical ensemble: provided the values of β andμ are selected so that k∈F k = k k f k and N k = k f k . For the kdistribution shown in Fig. 1, the comparison is shown in the Table I. We note in passing that if the thermal Renyi entropy density is computed in the canonical ensemble, they do not match S k,(n) A /M A for n > 1.
In conclusion, we have obtained the exact closed-form expression for the entanglement spectrum of exact manybody excited eigenstates of the Hubbard model. Despite being exact excited eigenstates with finite energy density above the ground state, these states violate strong ETH. This is either because their entanglement entropy is sub-extensive or because it is interaction independent. Nevertheless, despite an exponentially large number of these states [9], the fraction of these state in the Hilbert space vanishes in the thermodynamic limit. As such, they may provide a useful starting point for studying the onset of thermalization due to perturbations which violate the total pseudospin conservation. We

SUPPLEMENTARY INFORMATION
In this Supplementary Information, we provide detailed derivations of our results discussed in the main article. We also discusses extensions to the "generalized" Hubbard model and to a larger class of spin-flip eta-pairing states.

HUBBARD HAMILTONIAN AND SO(4) SYMMETRY
In this section, we introduce the extended Hubbard model and its symmetries. We consider a Hamiltonian where the kinetic energyT isT while the potential energy is a density-density "shifted" interaction: The case considered in the main text of the paper is β = 0, but for now we keep a generic β. For notation simplicity, bold symbols (such as r or β) represent vectors in the d-dimensional space. We also define the shifted momentum where G is a given vector on the lattice, to be determined later.

EXPLICIT EIGENFUNCTIONS OFĤ
The spectrum ofĤ can be placed in eigenvalues ofĴ 2 ,Ĵ z ,Ŝ 2 ,Ŝ z ,Ĥ,P . We will now write down a large number of exact eigenstates. Out ofĴ z ,Ŝ z we can make two (linearly dependent operators) quantum numbers,N ↑ = rĉ † r↑ĉ r↑ The numbers of ↑ and ↓ spins are independently conserved and will be used to interchangeably denote states.

Set of Eigenstates
First consider the eigenstates of the Hamiltonian for which the number N ↑ of ↑ particles is zero. For these states, the interaction do not appear and we have N ↓ = N k noninteracting fermions, each with their momenta: where F consists of any set of N k wavevectors in the 1 st Brillouin zone. The energy and momentum of these states is: We can see that these states have:η and hence these are the lowest weight states of a multiplet: with N 1 = 0, . . . , M − N k , N 2 = 0, . . . , N k . These states have the quantum numbers underĤ,P ,Ĵ z ,Ŝ z respectively: He then proceeded to building another set of states,η N βη † α |0 which he then proved were also eigenstates of the Hamiltonian. This state is, however, linearly dependent on Ψ {k} N1,N2 and hence not a new state. CN Yang's state |α =η † α |0 , ∀ α = β has the J z , S z quantum numbersĴ z |α = 2−M 2 |α ,Ŝ z |α = 0. It belongs to a multiplet |N 3 , N 4 , α = (η † β ) N3 (ζ † β ) N4 |α where and N 3 = 0, . . . , M − 2 and N 4 = −1, 0, 1 (the notation is (ζ † β ) −1 =ζ β ). By direct calculation we find that the lowest weight states |0, −1, α are (an energy −2µ) linear combination of the states Ψ {k} N1,N2 = |0, 0; k, π − k of Eq. 48:ζ This in fact had to be so because we will now prove that the states Eq. 48 are the only states (and their η andζ ) in the J + S = M/2 sector. Since |N 3 , N 4 , α also have J + S = M/2, they must hence be a linear combination of the states in Eq. 48 .

Completeness
The states in Eq. 48 have J, S quantum numbers that satisfy the relation J + S = M/2: given a configuration of momenta F, they have the same J, S quantum numbers Eq. 47 as their lowest weight counterparts Eq. 44, which immediately satisfy the aforementioned identity. We now prove that all the states with J + S = M/2 are part of multiplets where the lowest (and highest) weight states are noninteracting. In other words, the set of states in Eq. 48 saturate the Hilbert space of quantum numbers J + S = M/2 (irrespective of J z , S z , P ) Pick any states |J, S, J z , S z , P in the Hilbert space of the Hubbard model, with J + S = M/2. We can always apply theη β ,ζ β lowering operators the appropriate amount of times to bring this state to the lowest weight of both SU (2) ⊗ SU (2): |J, S, J z = −J, S z = −S, P . For this last state, we have J z + S z = −(J + S) = −M/2. Eq. 43 relates the quantum numbers to the number of ↑, ↓ particles existent in the system. It is then trivial to see that the states |J, S, J z = −J, S z = −S, P has N ↑ = 0 (without any restriction on N ↓ ). Since only N ↓ particles are present, there is no Hubbard U interaction, and all the lowest weight states |J, S, J z = −J, S z = −S, P (J + S = M/2) can be labeled by the momenta of the ↓ particle, as in Eq. 44. No other states can exist in these quantum number sectors.

States
To fully define the states, we compute their norm. We first present a method which will be used extensively in the calculations in this section. First, we note that we can write a Kronecker δ-function using contour integration as where the contour C encircles the origin in the complex z-plane counterclockwise. Using this and the commutation relations [ĉ † r ↓ĉ † r +β↑ ,ĉ r↓ĉ † r+β↑ ] = 0 we re-write the states: The product over r is taken over all lattice sites. We are now in a position to calculate the norm. Using the fact that Ψ {k} 0,0 contains only N k b-particles, and with the help of the identity we find The integrand can be massaged and provides the first Kronecker delta function of the momenta configurations F and F . Simple integration then provides for:

STATES
The entanglement spectrum of the states Ψ {k} N1,N2 can be analytically computed for β = 0. We sketch this calculation in the current section. The strategy to diagonalized the reduced density matrix will be to first perform the Schmidt decomposition of the non-interacting lowest weight states Ψ and apply techniques similar to those of Eq. 53 and Eq. 57 to diagonalize the density matrix. In all our calculations, we will assume that a partition of the space of M sites has been performed in an A and B parts.
The strategy for performing a Schmidt decomposition of a non-interacting state Ψ and B (right) regions is well known and has been first presented by Peschel in Ref. 18. We build and diagonalize the one-body density matrix of the A side: where we normalize our complete basis: k∈F φ m (k) φ m (k) = δ mm . Any eigenvalues which are 0, 1 and their respective eigenstates are discarded. Using this complete basis we want to build eigenstates with support fully in either region A or B. If for any γ m = 0, 1, we rescale and we have found normalized operators in the A and B side of the system: And similarly for the B region.
We are now ready to Schmidt decompose the state. As φ m (k) is a unitary transformation (keep all the γ m 's, even if 0, 1), we perform the canonical transformation : We now separateĉ † m into orthogonal left and right second quantized operators: Where a m = r∈A φ m (r)ĉ r , b m = r∈B φ m (r)ĉ r are canonical fermionic operators with support exclusively on the left and right hand side respectively. Hence: No ↑ fermions are present in the state. The many-body Schmidt decomposition of the state can be decomposed in sectors that contain N A particles in the A side and N k − N A particles on the B side. Each N A sector on the A side can be obtained by filling a set of {m A } single-particle eigenstates m. Written like this, the state is easily decomposed in where the sum is over all the different 2 N k ways to partition the N k ↓-"orbitals" m into those occupied by a † 's, denoted by the set {m A } -for the N A particle state |{m A } = m∈{m A } a † m |0 , and the complementary set {m B } occupied byb † 's for the N k − N A -particle state |{m B } = m∈{m B }b † m |0 . Here The entanglement entropy is then For the below, it is important to remember that N A is a good quantum number of the decomposition. Having obtained an A/B decomposition of the states Ψ {k} 0,0 , we now obtain the decomposition for the full states Ψ {k} N1,N2 . We start by building a new orthonormal basis for the A and B sides away from the lowest weight limit. We build theη andζ operators on the A and B sides respectively: Using the Schmidt decomposition of left and right parts in Ψ {k} 0,0 Eq. 65, we define the eta-pairing states of the A and B sides: Using the same steps as in Eq. 57, it is easy to prove orthonormality of |n 1 , n 2 , {m A } : We now can find, using the normalized Ψ {k} N1,N2 we find (the limits in the sum are obvious, for example in binomial coefficients, etc ): n3,n4,n5,n6 1 n 3 !n 4 !n 5 !n 6 ! z n3 1 z n4 2 z n5 3 z n6 where |N 1 − n 1 , N 2 − n 2 , {m B } is a normalized state of N k −N A −N 2 +n 2 +N 1 −n 1 ↓ particles and N 1 −n 1 +N 2 −n 2 ↑ particles on the B-side. The limits in the summations over n 1 , n 2 are implicit from the binomial formulas. The above expression gives the exact entanglement spectrum. By Vandermonde identity, one can check that the trace of the density matrix is unity. With the exact entanglement spectrum, it straightforward to obtain the expression of the Von Neumann entanglement entropy where λ k n1,n2 = Note that for N 1 = N and N 2 = 0, Eqs. 73 and 74 reduce to Eqs. 21 and 19 in the main text.

THERMODYNAMIC LIMIT AND SCALING
In this section, we give a detailed derivation of the entropy formula in the thermodynamic limit discussed in the main text. For sake of simplicity, we will focus on the case where N 1 = N and N 2 = 0.

Von Neumann Entropy in the Thermodynamic Limit
The factor {m A } α 2 {m A } in Eq. 73 is peaked about N * A ≈ M A M N k which is large. Therefore, M A − N A in λ k n,0 is large, and so is M − N k − (M A − N A ). This forces the peak of λ k n,0 to appear at large n. Thus we can use the Stirling's approximation, which gives So,