The Hubbard model provides a framework for investigating strongly correlated electron systems [1], particularly in two dimensions where it gives valuable insights into the mechanisms underlying high-temperature superconductivity [2]. Its Hamiltonian has two terms: the first term describes the hopping of electrons between different sites, and the second term describes the interaction when two electrons are located at the same site. The Hubbard model is difficult to solve, a challenge stemming from the competition between the two terms in its Hamiltonian. A prominent approach to studying it is quantum simulation using ultracold fermions confined in optical lattices [3–8]. Recent advances have realized such a simulator featuring approximately 800 000 sites [9].

Although the Hubbard model is very difficult to solve, there has been work aimed at obtaining its exact solutions, in addition to numerical and simulation approaches. Yang makes a pioneering work to the Hubbard model on a hypercubic lattice by constructing some exact eigenstates of the model [10]. His approach relies on the introduction of the η-pairing operator, which serves as an eigenoperator of the system’s Hamiltonian. There are also some efforts to find new exact eigenstates of the Hubbard model using other methods [11]. The η-pairing method has been successfully extended to construct exact eigenstates for certain Hubbard models on triangular lattices [12]. More recently, the conditions under which the η-pairing method can effectively identify exact eigenstates for Hubbard models on general graphs have been explored [13], and there are discussions about the η-pairing method on non-Hermitian Hubbard models [14, 15] and multicomponent Hubbard models [16].

To further enrich the Hubbard model and explore related physics, additional interactions can be incorporated. One such interaction is the bond–charge interaction, which modifies the hopping strengths depending on the occupation across the bond [17]. The inclusion of bond–charge interactions in the Hubbard model has been used as a model of hole superconductivity [17], and some exact eigenstates for this extended model have also been constructed [11]. Quantum many-body scar states have received many attentions in recent years, as these states do not obey the eigenstate thermalization hypothesis [18, 19]. It has been demonstrated that the Hubbard models with bond–charge interactions on hypercubic lattices can support η-pairing states as their eigenstates [20]. These states are true scar states, as the bond–charge interaction breaks the η-pairing SU(2) symmetry inherent in the original Hubbard models [20].

In this paper, we give the condition for the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of a Hubbard model with bond–charge interaction on a general graph, which is based on a finding that the condition for the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of a Hubbard model with bond–charge interaction is identical to the condition for it to be an eigenstate of a Hubbard model without bond–charge interactions. We also identify several eigenstates of the Hubbard model that are distinct from the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $. Our results, derived for general graphs, are applicable not only to periodic lattices in condensed matter physics but also to single molecules described by a Hubbard model.

Consider a Hubbard model with bond–charge interaction defined on a general graph, the Hamiltonian is where m and n denote sites on the graph, a_n and b_n are annihilation operators for spin-up and spin-down electrons on site n, and ${a}_{n}^{\dagger }$ and ${b}_{n}^{\dagger }$ are the corresponding creation operators. The ordinary hopping T_mn is The bond–charge interaction B_mn is The notation ⟨m, n⟩ means the sum is over the pair of sites connected on the graph.

Divide the bond–charge interaction B_mn as B_mn = B_mn1 + B_mn2 with and We can write where T_mn + B_mn1 has the same form as the hopping T_mn but with the strength t_mn + χ_mn. So adding the bond charge interaction B_mn to the Hubbard model can be understood as changing the hopping strength from t_mn to t_mn + χ_mn, and adding the leftover interaction B_mn2.

Now we discuss the condition for the η-pairing states to be eigenstates of the Hamiltonian H. The η-pairing operator is defined as There is a good relation between the operator B_mn2 and η^†, i.e.which is proved as follows. In the following n ≠ m, simple calculations show that there are The sum of the above four equations leads to In equation (13), interchanging m and n leads to In equation (13), interchanging a and b leads to In equation (15), interchanging m and n leads to Notice the expression of B_mn2 in equation (5), a linear combination of the above four equations leads to where As there is the relation in equation (8) is proved.

The relation in equation (8) is important as it can lead to the result that the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of B_mn2 with zero eigenvalue. This is because that from equation (8) we can get where N is an integer and ∣0⟩ is the vacuum state. The above equation means that there is which itself can be used N times to get Therefore the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of B_mn2 with zero eigenvalue.

To see the consequence that the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of B_mn2 with zero eigenvalue, define which is just the Hamiltonian of a Hubbard model with the hopping strength t_mn + χ_mn and without bond–charge interactions. As there is if the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of H it will also be an eigenstate of H_m and vice versa. We can draw the conclusion that the condition for the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of H is the same as that for it to be an eigenstate of H_m.

To find the condition for the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of H_m, it is needed to calculate the communicator [H_m, η^†] as Yang’s method [10]. The on-site interactions of H_m satisfy which can be easily verified. To calculate the communicator [T_mn + B_mn1, η^†], we note that there is Direct calculation shows that there is In equation (28), interchanging m and n leads to A linear combination of equations (28) and (29) gives with Based on the results in equations (26, 30), we can make the summary that When there is we have [H_m, η^†] = Uη^†, which can lead to so that ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of H_m and H with the eigenvalue NU. Therefore equation (33) is the condition for the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of H_m and H.

The condition in equation (33) gives a constraint on the hopping strength t_mn, the bond–charge interaction strength χ_mn and the phases ${{\rm{e}}}^{{\rm{i}}{\phi }_{n}}$ in η^†. In the special case that t_mn + χ_mn = 0 for all pairs ⟨m, n⟩, equation (33) gives no requirement on the phases in η^†, so η-pairing states with any phases in η^† are eigenstates of H. Such a special case is studied in detail in one dimension [21, 22]. In the following we assume that t_mn + χ_mn ≠ 0, and define where ${\tilde{t}}_{mn}$ and θ_mn are real. Now the condition in equation (33) for the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of H_m and H is equivalent to which will be discussed in the following from two different views.

Based on the condition in equation (36), we ask the following question: assume that the phases in η^† for all sites are specified, can we find a Hubbard model with bond–charge interaction such that the state ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of H? The answer is always YES. Such a model can be achieved when the phase of t_mn + χ_mn is where φ_m and φ_n are phases in η^†, as it satisfies the condition in equation (36).

Two interesting examples can be used to demonstrate this result. In the first example, define which takes φ_n = 0 for all n. If we want ${({\eta }_{1}^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of H, according to equation (37) we can set ${{\rm{e}}}^{{\rm{i}}{\theta }_{mn}}=\pm {\rm{i}}$, which means t_mn + χ_mn are all pure imaginary numbers. In the second example, define where the sets of sites A and B form a bipartition of the graph. In this example there is φ_n = 0 for all n ∈ A and φ_n = π for all n ∈ B. If we want ${({\eta }_{2}^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of H, according to equation (37) we can set ${{\rm{e}}}^{{\rm{i}}{\theta }_{mn}}=\pm {\rm{i}}$ when both m, n are in A or B, and ${{\rm{e}}}^{{\rm{i}}{\theta }_{mn}}=\pm 1$ when n ∈ A, m ∈ B, or m ∈ A, n ∈ B. This requirement means that the hopping strength t_mn + χ_mn within A or B are all pure imaginary numbers, and between A and B are real numbers. Therefore, for a Hubbard model defined on square lattices with nearest-neighbor hoppings and next-nearest-neighbor hoppings, if the nearest-neighbor hopping coefficients are real numbers and the next-nearest-neighbor hopping coefficients are pure imaginary numbers, the η-pairing state ${({\eta }_{2}^{\dagger })}^{N}| 0\rangle $ will be an eigenstate of the model Hamiltonian.

Based on the condition in equation (36), we can also ask the following question: assume that the parameters in H are specified, can we find an eta-pairing operator η^† such that the state ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of H? The answer depends on the phase θ_mn of t_mn + χ_mn.

Suppose sites (n₁, n₂, ⋯, n_k, n₁) form a loop on the graph, when the condition in equation (36) is satisfied, there is Conversely if the condition in equation (40) is satisfied for any loop (n₁, n₂, ⋯, n_k, n₁), we can always find proper ${{\rm{e}}}^{{\rm{i}}{\phi }_{n}}$ in η^† for any site n so that the condition in equation (36) can be satisfied and thus there is [H_m, η^†] = Uη^†. The method to find such ${{\rm{e}}}^{{\rm{i}}{\phi }_{n}}$ in η^† is as follows. First set ${\phi }_{{n}_{1}}=0$ for an arbitrarily chosen site n₁. Then according to equation (36) we set ${{\rm{e}}}^{{\rm{i}}{\phi }_{m}}=-{{\rm{e}}}^{{\rm{i}}{\phi }_{{n}_{1}}}{{\rm{e}}}^{2{\rm{i}}{\theta }_{m{n}_{1}}}$ for any site m connected to site n₁. And then according to equation (36) we set ${{\rm{e}}}^{{\rm{i}}{\phi }_{k}}=-{{\rm{e}}}^{{\rm{i}}{\phi }_{m}}{{\rm{e}}}^{2{\rm{i}}{\theta }_{km}}$ for any site k connected to site m. Because the condition in equation (40) is assumed to be satisfied for any loop, we can get a unique ${{\rm{e}}}^{{\rm{i}}{\phi }_{n}}$ for any site n by repeating this procedure.

Now we can answer the question raised at the beginning of this section. For a given H with hopping strength t_mn and bond–charge interaction strength χ_mn, if the phase θ_mn of t_mn + χ_mn satisfies the condition in equation (40) for any loop (n₁, n₂, ⋯, n_k, n₁) on the graph, we can always find an eta-pairing operator η^† such that the state ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of H_m and H.

In this section, we show that when the condition in equation (40) is satisfied for any loop (n₁, n₂, ⋯ , n_k, n₁) on the graph, we can find new exact eigenstates of the Hubbard-type Hamiltonian H_m in equation (24) that are different from the eta-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $.

When the condition in equation (40) is satisfied for any loop, as discussed in the above section we can find an η^† so that there is [H_m, η^†] = Uη^†. It will lead to the fact that ${({\eta }^{\dagger })}^{N}| 0\rangle $ is an eigenstate of H_m. Define T = T_a + T_b with two operators where T_a describes the hopping of spin-up electrons and T_b describes the hopping of spin-down electrons. We can write the Hubbard-type Hamiltonian H_m in equation (24) as We will show that there is for odd integer M, i.e. ${({T}_{{\rm{b}}})}^{M}{\eta }^{\dagger }| 0\rangle $ is an eigenstate of H_m with zero eigenvalue when M is an odd integer. This result together with [H_m, η^†] = Uη^† can lead to the conclusion that is an eigenstate of H_m with eigenvalue (N − 1)U, where M is an odd integer. It is obvious that this new eigenstate is different from ${({\eta }^{\dagger })}^{N}| 0\rangle $.

Now we prove the result that there is ${H}_{{\rm{m}}}{({T}_{{\rm{b}}})}^{M}{\eta }^{\dagger }| 0\rangle =0$ for odd integer M. As the condition in equation (40) is satisfied for any loop, there is [T, η^†] = 0, which can lead to Tη^†∣0⟩ = 0. This result together with the fact TT_b = T_bT can lead to for any positive integer M. Notice the expression of H_m in equation (43) and the result in equation (46), to show that ${H}_{{\rm{m}}}{({T}_{{\rm{b}}})}^{M}{\eta }^{\dagger }| 0\rangle =0$ for odd M we only need to show that there is no double occupation in the state ${({T}_{{\rm{b}}})}^{M}{\eta }^{\dagger }| 0\rangle $. Recall that ${\eta }^{\dagger }={\sum }_{n}{{\rm{e}}}^{{\rm{i}}{\phi }_{n}}{a}_{n}^{\dagger }{b}_{n}^{\dagger }$, there is which means that no double occupation in the state ${({T}_{{\rm{b}}})}^{M}{\eta }^{\dagger }| 0\rangle $ is equivalent to that ${({T}_{{\rm{b}}})}^{M}{b}_{n}^{\dagger }| 0\rangle $ is orthogonal to ${b}_{n}^{\dagger }| 0\rangle $. When M = 1 it is obvious that ${({T}_{{\rm{b}}})}^{M}{b}_{n}^{\dagger }| 0\rangle $ is orthogonal to ${b}_{n}^{\dagger }| 0\rangle $. For a general integer M bigger than 1, the state ${({T}_{{\rm{b}}})}^{M}{b}_{n}^{\dagger }| 0\rangle $ can be regarded as a spin-down electron initially at site n and then hopping M times. What is the amplitude for the electron going back to its initial site n after hopping M times? When M is an odd integer the way back to its starting site must contain a loop that has an odd number of sites, the electron can go back to its initial site along two opposite directions of the loop, the sum of the amplitude from these two opposite directions will be zero according to equation (40). Therefore when M is an odd integer, the state ${({T}_{{\rm{b}}})}^{M}{b}_{n}^{\dagger }| 0\rangle $ is orthogonal to ${b}_{n}^{\dagger }| 0\rangle $.

Finally, we note that if T_b is replaced by T_a in equation (45), the state ∣ψ_NM⟩ is still an eigenstate of H_m due to the symmetry between spin-up and spin-down electrons.

We have demonstrated that the condition for the η-pairing state ${({\eta }^{\dagger })}^{N}| 0\rangle $ to be an eigenstate of H in equation (1), which describes a Hubbard model with bond–charge interaction, is identical to the condition for it to be an eigenstate of the Hubbard-type Hamiltonian H_m in equation (24) without bond–charge interaction. For the Hamiltonian H given in equation (1), we have provided the conditions on its parameters under which the η-pairing method can be employed to construct its exact eigenstates. Specifically, if the parameters of H in equation (1) satisfy the condition in equation (40) for any loop, then the η-pairing method becomes applicable for constructing the exact eigenstates ${({\eta }^{\dagger })}^{N}| 0\rangle $, which are quantum many-body scar states [20]. When the condition in equation (40) is satisfied for any loop, we have shown that, in addition to ${({\eta }^{\dagger })}^{N}| 0\rangle $, there are new eigenstates ∣ψ_NM⟩ given in equation (45) for the Hubbard-type Hamiltonian H_m. These findings improve the knowledge about Hubbard models and the η-pairing method.