Attraction Versus Repulsion Between Doublons or Holons in Mott-Hubbard Systems

For the Mott insulator state of the Fermi-Hubbard model in the strong-coupling limit, we study the interaction between quasi-particles in the form of doublons and holons. Comparing different methods – the hierarchy of correlations, strong-coupling perturbation theory, and exact analytic solutions for the Hubbard tetramer – we find an effective interaction between doublons and/or holons to linear order in the hopping strength which can display attractive as well as repulsive contributions, depending on the involved momenta. Finally, we speculate about the implications of our findings for high-temperature superconductivity.


Introduction
Understanding strongly interacting quantum many-body systems is one of the major challenges of contemporary physics.In order to achieve progress in that direction, it is often useful to apply the same principles as for weakly interacting systems.Following this strategy, the first step is to find or characterize the ground or thermal equilibrium state [1][2][3].As the second step, one should identify the relevant quasi-particle excitations describing linearized perturbations around this equilibrium state and determine their properties, such as dispersion relations [4][5][6][7][8][9].Going beyond this linearized level, one can then study the interactions of these quasi-particle excitations among each other and with other degrees of freedom [10][11][12][13][14][15][16][17][18][19][20].
Since exact solutions are typically limited to special cases or small systems [21][22][23], approximations are necessary in most cases.Ideally, these approximation schemes should be based, at least in principle, on a systematic expansion into powers of some small control parameter.In contrast to weakly interacting systems, the large coupling strength prohibits its use as perturbation parameter, but one could use its inverse [24][25][26] (strong-coupling perturbation theory) or the inverse of some other large number, such as spin S ≫ 1 [27][28][29] or coordination number Z ≫ 1 [30][31][32][33][34][35], which typically leads to some sort of mean-field theory.
In the following, we consider the Fermi-Hubbard model as the drosophila of strongly interacting quantum many-body systems [38][39][40] where ĉ † µs and ĉνs denote the fermionic creation and annihilation operators at the lattice sites µ and ν with spin s ∈ {↑, ↓} while ns ν are the associated number operators.The lattice structure is encoded in the hopping matrix T µν which equals the tunneling strength T for nearest neighbors µ and ν and is zero otherwise.The coordination number Z counts the number of nearest neighbors µ for a given lattice site ν and is assumed to be large Z ≫ 1.Finally, U denotes the on-site repulsion and we focus on the strong-coupling limit U ≫ T in the following.
Let us briefly recapitulate the relevant symmetries of the Fermi-Hubbard Hamiltonian (1).In addition to the total particle numbers N s = µ ns µ , the total spin is also conserved, where σ ss ′ are the elements of the Pauli spin matrices reflecting the global SU (2)-invariance [41].Here, bold-face symbols such as Ŝ = ( Ŝx , Ŝy , Ŝz ) represent vectors.Another interesting symmetry is the particle-hole duality: If we exchange all creation and annihilation operators ĉ † µs ↔ ĉµs which implies ns µ ↔ 1− ns µ , we find that the Hamiltonian (1) is mapped to the same form with a negative hopping strength T ↔ −T up to an irrelevant shift containing the total particle number N = N ↑ + N ↓ .In order to avoid this shift, one could consider the grand-canonical Hamiltonian Ĥgc = Ĥ − µ N with the chemical potential µ = U/2 which is then mapped onto itself with T ↔ −T .
A prominent example for the crucial differences between weakly and strongly interacting systems is the Mott insulator [38,42].For weak interactions U ≪ T , the state at half filling for both spin species would be metallic, only the Fermi surface would be deformed a bit by the coupling U .The Mott insulator [43] is realized in the other limit U ≫ T however, where the ground state is insulating and basically one particle occupies each lattice site -up to small virtual hopping corrections with probabilities of order T 2 /U 2 .In order to facilitate transport, one has to excite a doublon-holon pair which requires a minimum energy given by the Mott gap ∆E Mott ≈ U .Note that, in contrast to these real and long-lived doublon-holon pairs (whose creation requires a minimum energy given by the Mott gap ∆E Mott ≈ U ), the hopping corrections mentioned above are sometimes pictured as virtual and short-lived doublon-holon pairs (which do not require such an excitation energy and are present in the ground state).
In the following, we shall study the properties of these quasi-particle excitations on top of the Mott insulating state [44].As explained above, this includes the single-particle characteristics such as their dispersion relation -but also two-particle properties describing their interaction among each other.

Hierarchy of Correlations
In order to pursue the strategy described in the Introduction, let us first employ the method of the hierarchy of correlations, see also [30,33,45].To this end, we consider the reduced density matrices of one ρµ , two ρµν , and three ρµνλ lattice sites, etc., and split up the correlated parts via ρcorr µν = ρµν − ρµ ρν , and so on.Now, based on the assumption Z ≫ 1, we may employ an expansion into powers of 1/Z where we find that higher-order correlators are successively suppressed as , and so on.This hierarchy facilitates an iterative approximation scheme, where we may start from the exact evolution equations and so on for even higher orders, where the functions F n are determined by the Hamiltonian (1).
To lowest order O(Z 0 ), we may approximate the first equation by i∂ t ρµ = F 1 (ρ µ , 0) + O(1/Z).The solution to this equation obeying the required boundary conditions then yields the mean-field ansatz ρ0 µ as the starting point for calculating the higher orders in 1/Z.
In order to describe the Mott insulator state at half filling in the strong-coupling limit U ≫ T , we use the simple mean-field ansatz at zero temperature) In principle, the aforementioned virtual hopping corrections with small probabilities ∼ T 2 /U 2 for an empty |0 µ or full lattice site |↑↓ µ could be included as well, but we neglect them here.
Note that the above mean-field ansatz (4) is invariant under the particle-hole duality transformation mentioned after Eq. ( 2) and does not include any spin ordering to lowest order.A staggered mean-field ansatz which does display spin ordering could be introduced for bi-partite lattices via where A and B denote the two sub-lattices.This state describes an Ising type anti-ferromagnet where Ŝz µ Ŝz ν is minimized and the Z 2 symmetry ĉµ↑ ↔ ĉµ↓ is spontaneously broken.
Note, however, that the ground state of the Fermi-Hubbard model (1) does not display this Ising type but rather Heisenberg type anti-ferromagnetic order where Ŝµ • Ŝν is minimized instead of Ŝz µ

Ŝz
ν as in the Ising case.This ground state is invariant under the SU (2)-invariance generated by the total spin (2) instead of the broken Z 2 symmetry of the Ising case.Formally, the reduced density matrix of a single lattice site µ is given by the ansatz (4).The correlations between neighboring lattice sites µ and ν can be taken into account via ρcorr µν .As a more intuitive picture, the Heisenberg type anti-ferromagnet can be visualized as lying somewhere in between the fully ordered Ising-type state (5) and the state (4) without any spin order, see also Eq. ( 30) below.
The correlations ρcorr µν can be further suppressed for finite temperatures.For example, if the temperature is much larger than the effective anti-ferromagnetic interaction O(T 2 /U ) but still way below the Mott gap ∆E = O(U ), the ansatz (4) would basically reproduce the exact thermal density matrix.As another possibility, the coupling to an environment can effectively steer the system towards the state (4), see, e.g., [32,46].

Doublons and holons
To next order in 1/Z, we may derive the quasi-particle excitations by approximating the second equation (3) via i∂ t ρcorr µν = F 2 (ρ 0 µ , ρcorr µν , 0) + O(1/Z 2 ) which yields a linear equation for ρcorr µν .To solve this linear equation, it is useful to split the original annihilation operator into the annihilation operator fµ↑ of a full lattice site |↑↓ µ and the creation operator ê † µ↑ of an empty lattice site |0 µ .After a spatial Fourier transform (assuming infinitesize lattices), the linear equation for ρcorr µν can be mapped onto a set of linear equations for the operators fks and ê † For convenience, we have re-scaled all length scales with respect to the lattice spacing ℓ and thus the wave-numbers used here are dimensionless.Note that the wave-numbers k are also vectors, but -depending on the dimensionality of the lattice -possibly in a vector space of different dimension than the three-dimensional vector Ŝ in Eq. ( 2).Diagonalizing the above 2 × 2-matrix, the eigenvalues λ ± k yield the quasi-particle energies E ± k via where T k denotes the Fourier transform of the hopping matrix T µν .In the strongcoupling limit U ≫ T , these quasi-particle energies simplify to Starting from the grand-canonical Hamiltonian including the chemical potential, the matrix in Eq. ( 7) would contain ±U/2 on the diagonal and thus its eigenvalues would be lowered by U/2 such that the quasi-particle energies in Eq. ( 8) assume a more symmetric form In the following, we shall use the convention (8).Note that, starting from the mean-field ansatz (5) reflecting the perfect Ising type anti-ferromagnetic order, the quasi-particle energies would not contain such a contribution linear in T k but scale quadratically O(T 2 k /U ).As an intuitive picture, since neighboring sites always have opposite spins, the propagation of doublons or holons on such a perfectly spin-ordered background can only occur via second-order tunneling processes.In contrast, for Heisenberg type spin order (or an unordered state), there is a finite probability that neighboring lattice sites are occupied by particles with the same spin, such that doublons or holons can propagate via first-order tunneling processes.
The eigenvectors of the matrix in Eq. (7) determine the Bogoliubov transformation to the quasi-particle operators dks and ĥ † ks dks ĥ † with the rotation angle where dks is the annihilation operator for a doublon with energy E + k while ĥ † ks is the creation operator of a holon with energy E − k .Of course, the above derivation of the quasi-particle picture is not unique, one can also derive it via other means, e.g., the Hubbard approximation [4,38].However, the 1/Z-expansion provides a clear and controlled path to incorporate higher orders consistently.

Boltzmann equation
To first order in 1/Z, the time evolution of the operators dks and ĥks is simply governed by the trivial phase factors exp{−iE ± k t} such that their populations d † ks dks and ĥ † ks ĥks remain constant.Interactions such as collisions between these quasi-particles leading to a finite energy and momentum transfer would induce a redistribution of these populations and are thus not described within this first-order approach.To incorporate such interactions, one has to include higher orders in 1/Z.
To second order 1/Z 2 , one should take the three-point correlator ρcorr µνλ in the second equation (3) into account.Its time-derivative does also contain the four-point correlator ρcorr µνλκ , which is of order 1/Z 3 .Truncating the set of evolution equations (3) at this order, i.e., neglecting all terms scaling with 1/Z 4 or higher, we may apply basically the same steps (Markov approximation etc.) as for weakly interacting systems and arrive at a Boltzmann equation describing the redistribution of the quasi-particle populations d s k = d † ks dks and h s k = ĥ † ks ĥks .Focusing on the holon sector for simplicity, we find in the strong-coupling limit U ≫ T (where E − k ≈ T k /2) for the mean-field background (4), see App.A and B [45,47] Thus, even in the strongly interacting limit, the quasi-particle distributions obey a Boltzmann equation which has the usual interpretation: Two holons with opposite spins and initial momenta k and p collide with each other and are scattered to the final momenta k + q and p − q where q is the momentum transfer.Note that the scattering cross section ∝ (T k +T p ) 2 is actually independent of the momentum transfer q.For two holons with the same spin, we found a vanishing scattering cross section, i.e., they do not interact at this order.The term in the third line of Eq. ( 11) represents the inverse process and ensures the conservation of probability or total holon number.Energy conservation is implied by the Dirac delta distribution in the first line of Eq. (11).Since the above Boltzmann equation (11) assumes the standard form, it entails the usual consequences, such as the H-theorem describing thermalization etc.
Focusing on the doublon sector instead, one obtains precisely the same form of the Boltzmann equation (11) for d s k = d † ks dks instead of h s k = ĥ † ks ĥks , as expected from the particle-hole duality mentioned in the Introduction.Taking both sectors into account simultaneously also accounts for collisions between doublons and holons, see Appendices A and B.
Note that initial states which are spin polarized in σ x direction, for example, would also induce off-diagonal terms such as ĥ † k↑ ĥk↓ , see also [48].In the absence of such a spin polarization, however, these terms vanish initially and thus stay zero throughout the evolution because our equations of motion do not contain symmetry-breaking contributions such as magnetic fields.Thus, we omit these off-diagonal terms here.

Effective Hamiltonian
In order to compare the Boltzmann equation ( 11) obtained via the 1/Z-expansion with the standard derivation of Boltzmann equations for weakly interacting systems, let us construct an effective Hamiltonian which would reproduce Eq. ( 11) in this way.To this end, let us start with the usual fermionic creation and annihilation operators â † ks and âks and the standard ansatz for such an effective Hamiltonian If we now set V ↑↓ kpq = −(T k + T p + T k+q + T p−q )/2 as well as E k = T k /2, we would indeed recover Eq. ( 11) via the usual Born-Markov approximation.
However, a few cautionary remarks are in order.First, the standard derivation of Eq. ( 11) from Eq. ( 12) is based on the usual fermionic commutation relations between the operators â † ks and âks .In contrast, neither the doublon d † ks and dks nor the holon operators ĥ † ks and ĥks satisfy these commutation relations, see Eqs. ( 6) and ( 9).Second, in contrast to the weakly interacting case, both V ↑↓ kpq and E k scale with T k and are thus not really independent, which requires special care when justifying the Born-Markov approximation.It should also be noted here that the insertion of the simple replacement âµ↑ → âµ↑ (1 − â † µ↓ âµ↓ ) into the free Hamiltonian µνs T µν â † µs âνs does not yield the correct effective Hamiltonian (12).
As another point, the scattering cross section in the Boltzmann equation ( 11) is given by the square of the interaction matrix element |V ↑↓ kpq | 2 and thus does not uniquely determine the sign (or phase) of V ↑↓ kpq , e.g., whether the interaction is attractive or repulsive.For example, for doublons one should insert 2 into the effective Hamiltonian (12), which does, however, yield the same Boltzmann equation (11).

Perturbation Theory in T /U
In order to settle the sign ambiguity mentioned above, let us compare our results to strong-coupling perturbation theory, i.e., a power expansion in the small control parameter ǫ = T /U ≪ 1.To this end, we split the Hamiltonian (1) via Ĥ = ĤU + ĤT = Ĥ0 + Ĥ1 into an undisturbed part Ĥ0 = ĤU = O(ǫ 0 ) plus a perturbation Ĥ1 = ĤT = O(ǫ 1 ).For general matrix elements we employ the same power expansion of the states and analogously for |Ψ out .
Because all the states considered in this section satisfy Ĥ0 |Ψ in 0 = 0 and Ĥ0 |Ψ out 0 = 0, the first-order matrix elements simplify to Apart from the power expansion in ǫ, we have not made any assumptions regarding the states |Ψ in and |Ψ out , e.g., regarding their degeneracy.They could be the same states, where M would yield the energy expectation value, or they could be different states, where M would describe a transition matrix element.

Mott state
Let us start with the Mott state |Mott , which we take to be the ground state of the Fermi-Hubbard Hamiltonian (1) at half filling (but other choices would also be possible).In the quasi-particle picture, it describes the state without doublons dks |Mott = 0 and holons ĥks |Mott = 0.After a power expansion in ǫ the zeroth order |Mott 0 has exactly one particle per site, i.e., êks |Mott 0 = 0 and fks |Mott 0 = 0.The virtual hopping corrections mentioned in the Introduction are included in the first-order correction consistent with Bogoliubov transformation ( 9) between dks and ĥks on the one hand and fks and êks on the other hand.Obviously, the first-order energy shift vanishes such that the ground-state energy is of order T 2 /U .

One-holon state
The quasi-particle picture described above motivates the ansatz ĥ † ks |Mott for the state containing one holon.However, one should be a bit careful because the operators ĥks and ĥ † ks do not obey the usual commutation relations.Fortunately, the calculation of the first-order matrix elements (15) only requires the zeroth-order states where most of these difficulties are absent because the operators ĉks and ĉ † ks do satisfy the standard commutation relations.The normalization N k↑ can be derived from Mott| ĉ † α↑ ĉβ↑ |Mott 0 = δ αβ Mott| nα↑ |Mott 0 and is -independently of k -just determined by the total number of particles with spin ↑.
In analogy, we use the same ansatz for |Ψ out 0 with k ′ and same spin ↑ (all other matrix elements vanish) Since the hopping matrix T µν is only non-zero for µ = ν and the state |Mott 0 has exactly one particle per site, we may set α = µ and β = ν or vice versa in the sum In addition to the number correlator in the second line, we obtain the spin-flip term Ŝ− µ Ŝ+ ν 0 in the third line.If the lattice and the state |Mott 0 obey translational invariance, the expectation values only depend on the relative coordinate r µ −r ν and thus the sum over the centerof-mass coordinate r µ + r ν corresponds to momentum conservation δ kk ′ .In case of rotational invariance, the expectation values yield the same result for all pairs of neighbors µ and ν and thus the remaining sum over r µ −r ν just yields the Fourier transform T k of the hopping matrix, i.e., M ∝ δ kk ′ T k + O(ǫ 2 ).For the mean-field ansatz (4), we find M = δ kk ′ T k /2 + O(ǫ 2 ) which reproduces the holon energy (8) to lowest order for k = k ′ and vanishes for k = k ′ , reflecting momentum conservation.As an outlook, one could study the scattering of holons (i.e., k = k ′ ) by spin inhomogeneities via inserting a mean-field ansatz which breaks translational invariance.
If we replace the mean-field ansatz (4) by the Ising type anti-ferromagnet (5), we find that the first-order matrix elements vanish M = O(ǫ 2 ).Again, this is consistent with the quasi-particle picture because the quasi-particle energies do not contain a linear contribution in this case, as discussed after Eq. (8).

Two-holon state
Now let us consider initial |Ψ in 0 and final |Ψ out 0 states containing two holons, where we start with the case of opposite spins, as motivated by the Boltzmann equation (11).As usual in scattering theory, we envisage initial and final holon wave-packets which do not overlap but interact in an intermediate space-time region.Then, in straightforward generalization of the one-holon case, we use the following ansatz for their Fourier components and analogously for |Ψ out 0 with k 3 and k 4 .The resulting matrix elements read The expectation values in the second line are only non-zero if α, β, and ν are mutually different, and the same for γ, δ, and µ.We only get non-vanishing contributions if the triple {α, β, ν} is a permutation of the triple {γ, δ, µ}.In view of T µν = 0 for µ = ν, we are left with four permutations (and the sum over spin s).
This result is consistent with the effective Hamiltonian (12) where the first term on the right-hand side of Eq. ( 24) corresponds to the free propagation of the two holons with their quasi-particle energies E k while the second term describes their scattering with the effective interaction potential V ↑↓ kpq .The origin of this effective interaction potential is the fact that the sums are not independent of each other, e.g., the sum over α = γ is not independent of the remaining sums over µ = β and ν = δ because α, β, and ν must be mutually different to yield a non-zero expectation value (as explained above).As an intuitive picture, the presence of the ↑-holon at site α may effectively inhibit the hopping of the ↓-holon from site ν to µ and thus changes its energy -which implies an effective interaction.

Two-holon triplet state
For comparison, let us consider the state of two holons with the same spin.In complete analogy to Eq. ( 22), we use the ansatz Following the same steps as in the previous subsection, including the insertion of the mean-field ansatz (4), we find that only the matrix elements corresponding to the free propagation survive Again, this is consistent with the Boltzmann equation ( 11) which also did not contain scattering between holons of equal spin.

Spin correlations
So far, our results were based on the zeroth-order mean-field ansatz (4) which neglects all correlations between the lattice sites.Including such correlations leads to corrections to these results.For the two-holon triplet state, for example, an effective interaction V ↑↑ kpq can be obtained if we include correlations between lattice sites, i.e., go beyond the mean-field ansatz (4).Taking into account these correlations between two lattice sites µ and ν as encoded in ρcorr µν = O(1/Z), but neglecting all three-point correlators ρcorr µνλ = O(1/Z 2 ), we find where C ↑↑ q denotes the Fourier transform of the number correlations Mott| n↑ The sign of the correlations depends on the spin order of the background state.For anti-ferromagnetic order, n↑ is negative for nearest neighbors α and β but positive for next-to-nearest neighbors -while for (locally) ferromagnetic order, it would also be positive for nearest neighbors.
In analogy, we may derive the correlation corrections to the interaction between two holons of opposite spin where we have used the SU (2)-symmetry [41] of the Mott state Ŝx

Hubbard tetramer
Let us exemplify the above results for an analytically solvable example, the Hubbard tetramer consisting of four lattice sites in the form of a square (Z = 2) [46,49].Already in this simple case, the total Hilbert space has 4 4 = 256 dimensions and thus the Hamiltonian (1) can be represented by a 256 × 256-matrix.However, by using the conserved quantities such as the particle numbers N ↑ and N ↓ , the total spin (2) and pseudo-spin (see Appendix C), as well as the spatial symmetries, one may cast this Hamiltonian into a block-diagonal form consisting of matrices with maximum rank four -admitting analytic solutions.
In contrast to the previous sections, which were devoted to the case of half filling (only marginally disturbed by one or two holons), we shall now also consider filling factors of 3/8 (one holon) and 1/4 (two holons).

Mott state
In the sector N ↑ = N ↓ = 2, the ground state has an energy of −3T 2 /U + O(ǫ 3 ) and vanishing total spin S = 0 as well as pseudo-spin η = 0, see also [41].We identify this state with the Mott state |Mott , which then displays the lowest-order structure The first-order hopping corrections can be obtained from Eq. ( 17).Note that one should be careful with the above representation because the sign of the basis vectors such as depends on the chosen order of the fermionic operators.Note that there is also another state in this singlet sector with S = η = 0 and N ↑ = N ↓ = 2, which has a slightly higher energy of −T 2 /U + O(ǫ 3 ).

One-holon state
The one-holon states -as single quasi-particle excitations around the Mott state -are then identified with the eigenstates in the doublet sector with S = 1/2 and N ↑ = 2 and N ↓ = 1 (or N ↑ = 1 and N ↓ = 2).They have eigen-energies of ±T /2 + O(ǫ 2 ) and ± √ 3T /2 + O(ǫ 2 ).In the following, we shall omit the symbols O(ǫ 2 ) for brevity and just state the energies to first order in T .

Two-holon state
To study the two-holon states, let us first consider the case N ↑ = N ↓ = 1.These twoholon states lie in the singlet (S = 0) or triplet (S = 1) sector and have eigen-energies ± √ 2T , ±T , and zero.The fact that all eigen-energies obey the reflection symmetry T → −T is a consequence of the staggered gauge transformation mentioned in the Introduction.
Already on the level of the eigen-energies, we find that not all two-holon energies can be written as a sum of two one-holon energies -which can be interpreted as a signature of their interactions.For example, adding one ↑ holon with energy − √ 3T /2 to another ↓ holon with the same energy − √ 3T /2, one would expect a total energy of − √ 3T in the non-interacting case.However, such an energy is not contained in the spectrum.Instead, the lowest two-holon energy is − √ 2T .As an intuitive picture, the presence of the ↓ holon reduces the options for the ↑ holon to lower the energy via tunneling and vice versa.As a result, these two quasi-particles effectively repel each other in this case.
In addition to the eigen-energies, we may also consider the eigen-states.To this end, let us introduce the operators Acting on the lowest-order Mott state (30), these operators generate the lowest-order one-holon states with energies ∓T /2.However, if we generate a two-holon state by applying these operators twice ĉ↑+ ĉ↓+ |Mott 0 , we obtain an eigen-state with zero energy.This again supports the interpretation that one holon disturbs the hopping options for the other holon such that they repel each other.Note that the same zero-energy state can be obtained via ĉ↑− ĉ↓− |Mott 0 which is consistent with the staggered gauge transformation mentioned in the Introduction and would then lead to the interpretation that these two holons attract each other.Another zero-energy eigenstate can be obtained by ĉ↑+ ĉ↓− |Mott 0 which would fit to the non-interacting case.

Two-holon triplet state
To complete the picture, let us discuss the two-holon states for the case N ↑ = 2 and N ↓ = 0. Obviously, they are in the triplet sector S = 1 and the repulsion U does not play any role in this case.Thus the eigen-energies are the same as in the noninteracting case, i.e., ±T and zero.If we try to write these two-holon eigen-energies as the sum of two one-holon eigen-energies, we see that this works for some of the oneholon states (with energies ±T /2), but not for the others (with energies ± √ 3T /2), which can again be interpreted as a signature of their interactions.Even though the repulsion U does not play any role for the two-holon states with N ↑ = 2 and N ↓ = 0, it is important for the one-holon states.
As an example for the states, we can obtain a two-holon eigen-state via ĉ↓+ ĉ↓− |Mott 0 which has (exactly) zero eigen-energy.Consistent with the above considerations, this would correspond to a case where two holons do not interact.

Bardeen-Cooper-Schrieffer (BCS) Theory
Having obtained repulsive as well as attractive contributions to the interaction between the quasi-particles such as holons, let us now investigate possible implications for Bardeen-Cooper-Schrieffer (BCS) like pairing, which might be relevant for our understanding of high-temperature superconductivity [50][51][52].To this end, we assume a small but finite density of holons -corresponding to a filling factor slightly below half filling.
As one possible approach, one could start from an effective Hamiltonian such as in Eq. ( 12) and then follow a procedure very analogous to the standard BCS theory of superconductivity [53], see Section 6.2 below.However, as already explained in Section 3, one might object that the effective quasi-particle operators do not obey the standard commutation relations.

Variational ansatz
Thus, we shall first pursue a more conservative approach and employ a variational ansatz in order to see whether and when BCS like pairing could lead to a reduction of the energy.To this end, we use the following ansatz for the zeroth-order BCS state with the pairing (squeezing) operator ξ µν and a normalization N which is required because the above exponential is not unitary.At a first glance, this ansatz may appear a bit unusual, but using translational invariance of the ξ µν , we may cast it into a more familiar form where we have used the fact that the exponential factorizes and its Taylor expansion terminates after the first order due to the Pauli principle.
Unfortunately, the Mott state does not factorize in the k basis, rendering the calculation of expectation values difficult.Thus, we use a Taylor expansion for small ξ k in order to study in which direction the energy could be reduced.In addition, we employ strong-coupling perturbation theory as in Sec. 4 which yields expectation values that we have already calculated there In the first term, |ξ k | 2 just gives the number of holon pairs with the holon eigenenergies E − k which are, up to small correlation induced corrections, given by T k /2.The second term corresponds to their interaction.
In order to avoid disturbing the Mott background too much, we consider a small number of holons, which is consistent with the assumption of small |ξ k | 2 .Then, only states close the minimum energies E − k ≈ T k /2 should be occupied by holons.Assuming a square lattice where T k behaves as cos k x + cos k y , states around (k x , k y ) = (π, π) are filled up first in order to minimize the energy.For these states, the lowest-order interaction term (T k + T p ) is repulsive, such that the usual s-wave pairing mechanism would not lead to a reduced energy.
However, for d-wave order parameters ξ k , which behave as cos k x − cos k y , both k ξ k and k ξ k T k vanish due to the angular average and thus the lowest-order repulsion term (T k + T p ) cancels.The remaining correlations C ↑↓ k+p can then indeed favor dwave pairing of holons since it corresponds to an effectively attractive contribution.In order to see how such a d-wave pairing could lower the energy, let us Taylor expand C ↑↓ q for small momenta After insertion into the variational ansatz (34), the constant c 0 and quadratic c 2 contributions vanish after their convolution with the d-wave order parameters ξ k and ξ * p , but the quartic term c 4 does indeed generate a reduction of the energy Ĥ provided that it is positive c 4 > 0. This condition c 4 > 0 is satisfied for anti-ferromagnetic nearest-neighbor correlations, for which C ↑↓ q behaves as cos q x +cos q y such that c 4 > 0 and c4 = 0. On the other hand, a non-zero c4 could also support tilted d-wave pairing where ξ k behaves as sin k x sin k y .
In summary, starting with the Mott state and adding a small amount of holon pairs suggests an instability towards d-wave pairing (but not s-wave pairing).This is generated by the effectively attractive contribution to the interaction between holons stemming from the correlation C ↑↓ q .

Effective Hamiltonian
Of course, it would be desirable to go beyond lowest order in ξ k and to include the chemical potential µ etc.Note that µ is now meant to describe the effective chemical potential associated to the finite density of holons, not the chemical potential µ = U/2 in the grand-canonical Hamiltonian for the original fermions as discussed after Eq. ( 2).
Ignoring the problems associated with the effective Hamiltonian ( 12) for a moment, let us treat the holons as fundamental particles as described by the creation and annihilation operators â † k,s and âk,s which obey the usual fermionic commutation relations.Then we may start from the effective Hamiltonian ( 12) together with the effective interaction (29) and perform the same steps as in standard BCS theory, including the derivation of a gap equation.To this end, we make an ansatz for the BCS state which is of the standard form [53] Here |0 denotes the vacuum state âk,s |0 = 0 and the variational coefficients fulfil As usual, non-trivial solutions △ p are obtained if the interaction V ↑↓ k,−k,p−k contains attractive contributions, which are the correlation terms C ↑↓ k in Eq. (29).In order to estimate these correlations, we exploit the SU (2)-symmetry of the Mott state and the Lieb theorem [41] which states that Ŝ2 |Mott = 0.Then, neglecting correlations beyond neighboring sites implies C ↑↓ k ≈ T k /(16T ) in two dimensions.Again assuming a square lattice where T k behaves as cos k x + cos k y , holon states around the minimum at (k x , k y ) = (π, π) are filled up first.Shifting the origin to the minimum k x,y → k x,y + π, the energies E k scale quadratically for small k.Then we may seek for d-wave solutions of the gap equation ( 37) which do also scale quadratically k 2 y − k 2 x for small k.As usual, pairing is expected to be most pronounced in the vicinity around the Fermi momentum k F such that we restrict the integral (37) Linearizing the energies E k in this interval, we find Since we have re-scaled all length scales with respect to the lattice spacing ℓ, the above Fermi momentum k F is dimensionless.Restoring physical units would correspond to the replacement k F → ℓk F .As in the usual BCS theory, we obtain an exponential suppression of the gap, but now the exponent is not inversely proportional to the coupling strength (because the kinetic and the interaction energy both scale linearly in T ) but to the fourth power of the Fermi momentum, i.e., the holon number density squared.Two powers of k F stem from the volume element, the other two from the quadratic scaling of the dwave order parameter.The strong exponential suppression for small k F might indicate that a certain holon number density is required to observe d-wave pairing [54][55][56][57][58], but further investigations are needed to settle this issue.
Nonetheless, comparing our findings with the well-known phase diagram of cuprates, for example, we find qualitative consistency as superconductivity is usually associated with a region of finite holon doping at low temperatures.Of course, the range of applicability of the simple single-band Fermi-Hubbard model (1) must be taken into account in this regard.This becomes even more important for the opposite case of electron doping.The particle-hole duality discussed in the Introduction implies that BCS states for doublons should exist in the same way as for holons.However, the asymmetry of the phase diagram of cuprates with respect to electron doping versus hole doping already shows that these systems do not display this particle-hole duality (as is also well known) and thus requires a description beyond the single-band Fermi-Hubbard model (1).

Conclusions
Via a combination of approaches, we studied the interaction between doublons or holons as quasi-particle excitations (i.e., charge modes) of the Mott insulator state in the strongly interacting Fermi-Hubbard model.Using the hierarchy of correlations and the simple mean-field ansatz (4), we derived a Boltzmann equation ( 11) with a scattering cross section which is quadratic in the hopping strength T for doublons or holons of opposite spin (and zero otherwise).
This motivates an effective interaction V ↑↓ kpq whose strength is linear in T and which can be represented by an effective Hamiltonian of the form (12).Note that this effective Hamiltonian should be treated with special care: First, the doublon and holon quasi-particle operators dks and ĥks do not satisfy the standard commutation relations.Second, the Boltzmann equation does only contain the absolute value squared of the interaction strength |V ↑↓ kpq | 2 and thus does not determine its sign (attractive or repulsive) uniquely.
Although one might use continuity arguments to demonstrate that the effective interaction V ↑↓ kpq can be attractive as well as repulsive (depending on the momenta), we employed strong-coupling perturbation theory to infer V ↑↓ kpq including its sign.Inserting the simple mean-field ansatz (4), which neglects the correlations between lattice sites, we indeed recover the interaction V ↑↓ kpq in the effective Hamiltonian ( 12) and thus the Boltzmann equation (11) to lowest order in T /U .These calculations motivate the following intuitive picture: In the Mott insulator state, hopping is suppressed due to the Mott gap, such that the tunneling probabilities scale with T 2 /U 2 .Inserting a holon, however, the system can lower its energy by tunneling -which gives rise to the single-holon quasi-particle energies of order T .Two holons far away from each other lower the energy according to the sum of their quasi-particle energies.However, if they come too close, the presence of one holon can influence (suppress) the tunneling of the other holon and vice versa, such that the energy reduction changes -giving rise to an effective interaction.Obviously, starting from the Mott state (containing one particle per site), two holons cannot occupy the same lattice site.
Consistent with this picture, two holons with momenta k and p which both lower the energy separately T k < 0 and T p < 0 would repel each other while two holons which increase the energy separately T k > 0 and T p > 0 would attract each other.The above line of argument specifically applies to holons, but the particle-hole duality mentioned in the Introduction implies the analogous behavior for doublons after the substitution T → −T .(The interaction between a doublon and a holon is discussed in Appendix B.) As a result, if two holons with momenta k and p attract each other, two doublons with the same momenta would repel each other and vice versa.
For the simple example of the Fermi-Hubbard model on a square (Hubbard tetramer) admitting an analytic solution, we could confirm the above picture -at least qualitatively.While the energy of the Mott state (30) scales quadratically O(T 2 /U ), the eigen-energies of the states corresponding to one and two holons are linear in T to lowest order.Furthermore, the lowest (highest) two-holon eigen-energies cannot be written as a sum of two one-holon eigen-energies, indicating an effective repulsion (attraction).
Going beyond the simple mean-field ansatz (4) and taking spin correlations between the lattice sites into account, we obtain corrections to the quasi-particle energies E k as well as to their effective interaction.For example, these correlations do also lead to an interaction (27) between holons of the same spin.Furthermore, for two holons of opposite spin, the effective interaction V ↑↓ kpq , which is repulsive for low-energy holons, does also acquire attractive corrections (29) due to the spin correlations.As an intuitive picture, the fact that two holons cannot occupy the same lattice site leads to an effective on-site repulsion whereas the spin correlations can induce a finite range attraction: If a holon with spin ↑ occupies the lattice site µ, there must have been an electron with that spin ↑ in the Mott state at that lattice site µ.Then, in the presence of (even short-ranged) anti-ferromagnetic order, the probability for having an electron with the other spin ↓ in the Mott state at a neighboring lattice site ν is larger than average.Thus, this neighboring lattice site ν can support a holon with the other spin ↓.In addition to the on-site repulsion explained above, one can visualize this as a nearest-neighbor attraction.
Note that the observed scale separation between the fast frequency scale T of the propagation and interaction of the doublons and holons on the one hand and the slow frequency scale T 2 /U of the spin fluctuations on the other hand allows us to approximately treat the (fast) evolution of the doublons and holons as taking place on a background with a fixed spin structure.
Finally, we discussed the implications of our results for high-temperature superconductivity.Using a BCS-like variational ansatz, we found that the usual s-wave pairing would not lower the energy (due to the effective on-site repulsion) but d-wave pairing could actually reduce the energy as a result of the nearest-neighbor attraction.Within the effective Hamiltonian approach, we deduced a gap equation.Its solution for the d-wave gap displays the usual non-perturbative structure, but in terms of the Fermi momentum of the holons instead of a coupling strength.
Of course, the effective interaction between doublons and/or holons has already been discussed in many publications, see, e.g., [39] and references therein.By now it is commonly expected that the spin degrees of freedom play an important role in that respect.The major points specific to the present work are: First, the derivation of the Boltzmann equation (based on the 1/Z expansion) displaying scattering cross sections which scale quadratically in T and thus point to an effective interaction linear in T .Second, the derivation of this effective interaction (based on the 1/U expansion) which is indeed linear in T and contains attractive as well as repulsive contributions.Third, the resulting gap equation whose solution is also exponentially suppressed, but the exponent merely contains the holon density.

Outlook
There are many ways to generalize our results.As one example, we focused on the leading order (in 1/Z or T /U ).Including higher orders would lead to modifications in several places.For instance, the lowest-order mean-field ansatz (5) could be modified by including small probabilities for an empty or doubly occupied lattice site or that this lattice site is occupied by the "wrong" spin.In this way, the back-reaction of the quantum or thermal fluctuations onto the mean field can be taken into account.This, in turn, would change the lowest-order quasi-particle energies a bit, which corresponds to a renormalization of the involved quantities, quite analogous to the case of weakly interacting systems, see, e.g., [59][60][61][62].In a similar manner, one could include higher orders in T /U in Sec. 4.
As a somewhat related point, we considered the zero-temperature limit here.Finite temperatures can also be taken into account in the approach based on the hierarchy of correlations (as it deals with density matrices), for example via the double-time correlator, see, e.g., [63].The expected impact of finite temperatures can be discussed in terms of general arguments.If the temperature is well below the typical spin energy of order T 2 /U , one would expect that our results are basically unaffected.Once the temperature is above this energy scale T 2 /U , it is expected to wash out the antiferromagnetic correlations and thus the finite-range attraction (responsible for d-wave pairing) is suppressed while the on-site repulsion remains.The next characteristic scale is reached when the temperature approaches the hopping rate T leading to thermal broadening of the holon distribution functions.Finally, once the temperature reaches or even exceeds the Mott gap of order U , thermal excitations in the form of real doublon-holon pairs change the background (4) considerably such that the insulating behavior of the Mott phase disappears.
In this context, one should also remember that we took the magnetic order of the Mott background as given, i.e., fixed.As explained above, the rationale behind that is the separation of scales between the scale T of propagation and interaction of the holons and the characteristic scale T 2 /U of the spin fluctuations.However, a complete picture would also require a more detailed treatment of the spin fluctuations and the origin of the magnetic order.For example, even though T is much larger than T 2 /U in the strong-coupling limit considered here, the superconducting gap (39) scales as T exp{−32π/(3k 4 F )} and thus it could be smaller than T 2 /U for a very low density of holons.In this case, the spin fluctuations might even destroy superconductivity.
Closely related to the magnetic order is the lattice structure.Our approach can basically be applied to quite general lattices, as long as they obey the usual (discrete) translational symmetries.The pseudo-spin η and the anti-ferromagnetic order such as in Eq. ( 5) require bi-particle lattices.For lattices which are not bi-particle (e.g., a triangular lattice), the anti-ferromagnetic order would be suppressed due to frustration, but short-range anti-ferromagnetic correlations should still persist (although on a weaker level) and thus the finite-range attraction (responsible for d-wave pairing) may survive.Apart from the discussion of the Hubbard tetramer in Sec. 5 which is obviously devoted to this specific example, we assumed a square lattice in Sec.(6).For other lattice structures (e.g., hexagonal), one should adapt the Fourier components T k accordingly, which might then alter the rotational symmetries of the superconducting gap.
It should also be illuminating to compare the results of our approach with other methods based on a large-Z expansion such as dynamical mean-field theory (DMFT) [64][65][66][67] or its time-dependent version (t-DMFT) [68][69][70][71].As a first difference, this method usually considers a different scaling limit, i.e., a factor of 1/ √ Z instead of 1/Z in front of the hopping term in the Hamiltonian (1).As a consequence, already the limit Z → ∞ becomes non-trivial, while we are mostly interested in the corrections of order 1/Z or higher.Furthermore, such methods which are based on the mapping to an effectively single lattice site or a finite cluster of sites are quite suitable for deriving frequency-dependent quantities such as the self-energy -but are less adapted to the problem considered here, where the spatial structures and the momentum dependence play an important role.
Note that our considerations are solely based on the Fermi-Hubbard model without invoking any effective descriptions (such as the t-J model).However, it would be interesting to generalize our findings to other model Hamiltonians (see, e.g., [72,73]) and to compare the results.
Along these lines the Markov solutions of (A4) and (A6) can be obtained.Finally we employ (A2) which leads in the long-time limit to ∂ t G

Appendix B Doublons and holons
The explicit form of the scattering kernel S ABCD is a rather complicated expression.However, in the limit of strong interactions the Boltzmann dynamics simplifies to The simple Boltzmann form allows to construct an effective Hamiltonian from which (B10) and (B11) can be recovered using leading order perturbation theory together with the usual Markov approximation, Ĥeff = Here the operators âk,s and bk,s are the annihilation operators for holons and doublons, respectively.As shown in section (4.3), the effective interaction potential V ↑↓,hh kpq for holon-holon scattering can also be justified from leading order perturbation theory involving two-holon states.Similarly, taking as initial and final state two-doublon wave packets, one finds the effective interaction potential V ↑↓,pp kpq = −V ↑↓,hh kpq where the minus sign originates from the particle-hole symmetry.The holon-doublon scattering can be justified from holon-doublon wavepackets of the form which is consistent with the fourth line of the effective Hamiltonian (B12).The exchange interaction in the last line of equation (B12) was choosen in order to achieve full consistency with the Boltzmann equations (B10) and (B11).