Doublon-holon binding, Mott transition, and fractionalized antiferromagnet in the Hubbard model

Doublon-holon binding, Mott transition, and fractionalized antiferromagnet in the Hubbard model

Sen Zhou, Yupeng Wang, and Ziqiang Wang State Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China Department of Physics, Boston College, Chestnut Hill, MA 02467, USA The Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
September 8, 2019
Abstract

We argue that the binding between doubly occupied (doublon) and empty (holon) sites governs the incoherent excitations and plays a key role in the Mott transition in strongly correlated Mott-Hubbard systems. We construct a new saddle point solution with doublon-holon binding in the Kotliar-Ruckenstein slave-boson functional integral formulation of the Hubbard model. On the half-filled honeycomb lattice and square lattice, the ground state is found to exhibit a continuous transition from the paramagnetic semimetal/metal to an antiferromagnetic ordered Slater insulator with coherent quasiparticles at , followed by a Mott transition into an electron-fractionalized AF phase without coherent excitations at . Such a phase structure appears generic of bipartite lattices without frustration. We show that doublon-holon binding unites the three important ideas of strong correlation: the coherent quasiparticles, the incoherent Hubbard bands, and the deconfined Mott insulator.

pacs:
71.10.-w, 71.27.+a, 71.10.Fd, 71.30.+h

I Introduction

The fundamental theoretical challenge of the strong correlation problem is the description of both the low energy coherent quasiparticles (QPs) and the higher energy incoherent excitations, and the spectral weight transfer from coherent to incoherent excitations with increasing correlation strength. Two very important ideas, the emergence of two broad incoherent features known as the Hubbard bands and the existence of renormalized QPs with a Luttinger Fermi surface (FS) were advanced by Hubbard hubbard , and Brinkman and Rice brinkman , respectively. Unfortunately, the Hubbard equation of motion scheme that produces the incoherent spectrum fails to produce QPs correctly and violates Luttinger’s theorem luttinger ; whereas the approaches based on the Brinkman-Rice-Gutzwiller wave functions gutzwiller capture a Luttinger FS of QPs, but find serious difficulties in constructing variational excited states to account for the incoherent excitations. Faced with this enigma, numerical approaches such as exact diagonalization, quantum Monte Carlo (QMC), and the dynamical mean field theory (DMFT) dmft have played a key role in recent studies of the strong correlation problem.

In this paper, we develop new analytical insights and construct a unified theory for both the coherent and incoherent excitations as well as the magnetic and the Mott transition. Our focus will be the half-filled single-band Hubbard model on bipartite lattices without frustration. As specific examples, we study the square lattice and the honeycomb lattice especially in view of the recent debate over the possible emergence of a gapped spin liquid (SL) phase on the honeycomb lattice meng ; wliu ; liebsch ; zylu ; jxli ; hur12 ; sorella ; seki ; hassan13 ; liebsch2 . With only on-site interactions, the Hilbert space of the Hubbard model is a product of the local Hilbert space on a single-site that consists of empty (holon), doubly occupied (doublon), and singly occupied states. The Brinkman-Rice-Gutzwiller approach amounts to a metallic state where the holon, denoted as a boson , and the doublon, as , condense fully with macroscopic phase coherence, as can be obtained by the Gutzwiller approximation gwa-finiteU or the slave boson mean-field theory kr . The metal-insulator transition is thus forced to follow a route where the density of doublons and holons vanish together with the condensates: such that there is exactly one electron per site. As a result, single-particle motion, coherent or incoherent, is completely prohibited. This so-called Brinkman-Rice (BR) transition is different from the Mott transition induced by the complete transfer of the coherent QP weight into the incoherent background, i.e. the depletion of the condensate while keeping the doublon/holon (D/H) density finite in the Mott insulator.

We will show that the crucial physics uniting the disparate ideas of Hubbard and BR is the binding between doublon and holon: . In the Mott insulator at large , although the D/H condensate vanishes () together with the disappearance of the coherent QP, the D/H density remains nonzero (). The motion of the QP is thus possible by breaking the doublon-holon (D-H) pairs, giving rise to the higher-energy incoherent excitations. With decreasing , the D/H density increases and the D-H binding energy decreases. At a critical , the D-H excitation gap closes and the D/H single-particle condensate develops, marking the onset of the Mott transition. On the metallic side of the Mott transition, D-H binding continues to play an important role since an added electron can propagate either as a coherent QP via the D/H condensate or incoherently via the unbinding of the D-H pairs.

The idea that D-H binding plays an important role in Mott-Hubbard systems was introduced by Kaplan, Horch, and Fulde khf82 and studied in the context of improved variational Gutzwiller wave functions yokoyama ; capello05 . The difficulty in constructing the appropriate variational wave functions for excitations has prevented further advances along these lines. More recently, field theory approaches involving the binding of charge 2 doublons with fermionic quasiparticles phillips09 as well as combining the bosons with fermions to form co-fermions imada11 have been put forth within the context of doped Mott insulators.

In this work, we will show that the physical picture presented above can be realized in the slave-boson functional integral formulation of the Hubbard model introduced by Kotliar and Ruckenstein (KR) kr by constructing new saddle point solutions that include the D-H binding. This approach also offers a treatment of the magnetism at half-filling that compares well to QMC simulations werner and has the added advantage of allowing the study of excitations and finite temperature properties otherslaveboson . As concrete examples, we studied the D-H binding in the half-filled Hubbard model on the honeycomb and the square lattice at zero temperature and obtained the phase diagram shown schematically in Fig. 1. On the honeycomb lattice, a continuous transition from the semimetal (SM) to an antiferromagnetic (AF) ordered insulator takes place at a critical Coulomb repulsion , suggesting that the gapped SL phase proposed by Meng et. al. meng may correspond to an AF ordered phase in the thermodynamic limit. Sorella et. al. sorella recently extended the QMC and the finite size scaling analysis to much larger system sizes and discovered that the signature of the gapped SL disappears and is replaced by that of a continuous SM to AF order transition at , in qualitative agreement with our results. Remarkably, we found a second quantum phase transition at a critical beyond which the D/H single-particle condensate vanishes () amid a finite density of doublons bound to holons. For , a new AF phase without coherent QP excitations, termed as AF in Fig. 1, emerges where the electrons are fractionalized and the elementary excitations do not carry the quantum numbers of an electron. We obtained a similar phase diagram on the square lattice; the transition to the Slater AF state happens for infinitesimal (i.e. ) due to the perfect nesting of the Fermi surface on the square lattice while the transition to the AF* phase takes place at .

The rest of the paper is organized as follows. Section II describes the slave-boson path integral formulation of the Hubbard model, the KR saddle point solution on the honeycomb and the square lattice, and the BR metal-insulator transition. In section III, we introduce the slave boson intersite correlations into the functional integral and construct the new saddle-point solution that includes the effects of D-H binding. The Mott transition and the spectral weight transfer between coherent and incoherent excitations will be studied to obtain the phase diagrams of the Hubbard model on the two half-filled bipartite lattices. We describe the transitions between the paramagnetic metal/semimetal, Slator AF insulator, and AF phases, elucidate the properties of the electron-fractionalized AF phase, and develop further insights into the nature of the incoherent excitations in Mott-Hubbard systems. Section IV contains the summary and conclusions.

Ii Hubbard model and slave boson functional integral representation

We start with the Hubbard model with nearest neighbor hopping and on-site Coulomb repulsion ,

 ^H=−t∑⟨i,j⟩,σ(c†iσcjσ+h.c.)+U∑i^ni↑^ni↓, (1)

where creates an electron with spin on site , and is the density operator. In the KR formulation kr , the electron operator is written as

 ciσ=^ziσfiσ,^ziσ=^Liσ(e†ipiσ+p†i¯σdi)^Riσ, (2)

where the boson operators describe the holon (), doublon (), and singly-occupied () sites, and is a fermion operator. The operators and are diagonal with unit eigenvalues in the (empty, ) and the (, doubly-occupied) subspaces, respectively lavagna ,

 ^Liσ=(1−d†idi−p†iσpiσ)α,^Riσ=(1−e†iei−p†i¯σpi¯σ)β,

where and can take any value. The Hubbard model Hamiltonian is thus given by,

 ^HKR= −t∑⟨i,j⟩,σ(^z†iσ^zjσf†iσfjσ+h.c.)+U∑id†idi. (3)

The partition function is a coherent state path integral over the quantum fields normalorder :

 Z=∫D[f,f†]D[e,e†]D[p,p†]D[d,d†]D[λ,λσ]e−∫β0Ldτ, (4)

where the Lagrangian is given by

 L= ∑i(e†i∂τei+d†i∂τdi)+∑i,σ(p†iσ∂τpiσ+f†iσ∂τfiσ) + ^HKR+i∑iλi^Qi+i∑i,σλiσ^Qiσ−μ∑iσf†iσfiσ, (5)

where is the chemical potential and and are the Lagrange multipliers introduced to enforce the local constraints for the completeness of the Hilbert space:

 ^Qi =e†iei+∑σp†iσpiσ+d†idi−1=0, (6)

and the equivalence between the fermion and boson representations of the spin-dependent density:

 ^Qiσ =f†iσfiσ−p†iσpiσ−d†idi=0. (7)

The KR saddle point corresponds to condensing all the boson fields uniformly with their values determined by minimizing the action kr . KR found that for , the saddle point solution recovers the Gutzwiller approximation kr .

The KR saddle point solutions on the half-filled honeycomb lattice GWhoneycomb are summarized in Fig. 2. Restricting to the paramagnetic (PM) phase, the doublon density decreases linearly from 1/4 at and vanishes at the BR metal-insulator transition . When magnetism is allowed, a SM to an AF insulator transition arises at . The D/H condensate remains nonzero for all finite and the AF phase is a Slater insulator with coherent quasiparticle excitations. The results on the half-filled square lattice are shown in Fig. 3. The BR metal-insulator takes place at . When magnetism is allowed, the PM metal is unstable with respect to the Slator AF insulator for any nonzero owning to the perfect nesting of the Fermi surface; the AF ordered moment develops exponentially at .

Iii Boson inter-site correlations and new saddle-point solutions with doublon-holon binding

The KR saddle-point solution, i.e. the Gutzwiller approximation, ignores all inter-site correlations and captures only the coherent QP single-particle excitations. Indeed, it has been shown GWinf-d that in the limit of infinite dimensions (infinite-), the latter becomes an exact solution of the variational Gutzwiller wave function approach. Our strategy for going beyond the Gutzwiller approximation described by the KR saddle point is to build in explicitly the inter-site correlations and boson dynamics in the functional integral and construct a new saddle point solution that includes D-H binding.

iii.1 Path integral including boson inter-site correlations

Introducing the operators for the D-H pairing , and the D/H hopping , on the nearest-neighbor bonds, as well as the density operators , on each site, we can rewrite the bosonic part in the hopping term in Eq. (3) as

 ^z†iσ^zjσ =^Y−1/2ij,σ[(^χeij)†p†iσpjσ+^χdijpi¯σp†j¯σ (8) +^Δjip†iσp†j¯σ+^Δ†ijpi¯σpjσ]^Y−1/2ji,σ,

where

 ^Yij,σ =^R−2iσ^L−2jσ=[(1−p†i¯σpi¯σ)(1−p†jσpjσ) (9) −^nei(1−p†jσpjσ)−^ndj(1−p†i¯σpi¯σ)+|^Δji|2].

Note that, due to the normal ordering of the square roots, the expression for involves explicitly the D-H pairing, but not the H/D hopping operators. The obvious challenge is how to build these correlations into the calculation of the partition function. Since they enter through the rather formidable factor , the usual procedure of introducing the corresponding correlation fields (, , , , ) via Hubbard-Stratonovich transformations in the path-integral does not work here. We found that the difficulty can be overcome by introducing in the functional integral additional Lagrange multipliers in the corresponding channel, , , , , and , such that the partition function becomes:

 Z= ∫D[f,f†]D[e,e†]D[p,p†]D[d,d†]D[λ,λσ] (10) D[Δ,χd,χe,nd,ne]D[Δv,χd,v,χe,v,ϵd,vϵe,v]e−∫β0dτL,

with the Lagrangian

 L =∑i(e†i∂τei+d†i∂τdi)+∑i,σ(p†iσ∂τpiσ+f†iσ∂τfiσ) +^HDH+i∑iλi^Qi+i∑i,σλiσ^Qiσ−μ∑iσf†iσfjσ, (11)

where is the effective D-H binding Hamiltonian

 ^HDH= −t∑⟨i,j⟩,σ(ziσzjσf†iσfjσ+h.c.)+U∑id†idi −i∑⟨i,j⟩[χd,vij(d†idj−χdij)+χe,vij(e†iej−χeij) +Δvij(diej−Δij)+Δv% ji(eidj−Δji)+h.c.] +i∑i[ϵd,vi(d†idi−ndi)+ϵe,vi(e†iei−nei)]. (12)

The factor in Eq. (12) has the same form given in Eqs. (8) and (9) but with the bosonic operators replaced by their corresponding correlation fields (). A few remarks are in order. (i) Eqs. (10-12) provide an exact representation of the Hubbard model; carrying out formally the last two functional integrals in Eq. (10) recovers the KR formulation given in Eqs. (3-5). (ii) The intersite correlations of the bosons can be included in a similar manner. For simplicity, we treat the latter as condensed fields in this paper since their densities (i.e., the density of single occupations) remain large at half-filling. (iii) From the perspective of finding a saddle point solution of the action, the effective Hamiltonian in Eq. (12) can be understood intuitively as a variational Hamiltonian describing the effects of intersite correlations of the doublons and holons, including that of D-H binding, where , , , , and are nothing but the variational parameters to be self-consistently determined.

iii.2 Saddle point solutions with D-H binding

We next discuss the D-H binding saddle point solutions of the path integral in Eqs. (10-12) which correspond to configurations of the quantum fields that minimize the action. We consider here that translation invariant PM and the two-sublattice AF saddle point solutions on the half-filled bipartite (honeycomb and square) lattices with sites. The bond variables are taken to be real and isotropic, and symmetry requires , , , and correspondingly, , , . Moreover , , and , where and denote the two sublattices on the bipartite lattice. Consequently, on the nearest neighbor bonds , the factor

 tziσzjσ=tg[2p0↑p0↓χd+(p20↑+p20↓)Δd]≡χvf, (13)

where , with

 Yσ=1−2nd−2p20σ+2p20σnd+p40σ+Δ2d. (14)

As will be shown later, this expression ensures that the new saddle point solution recovers the noninteracting limit at . Substituting these quantities into Eq. (12), we obtain the saddle point Hamiltonian,

 ^HspDH =^Hf+^Hd+4Nζ(χvdχd+ΔvdΔd) (15) −2N(ϵvd+2λ)−4Nϵvdnd+2N∑σ(λ−λσ)p20σ

where the coordination number on the honeycomb lattice, and on the square lattice. The Hamiltonian and determines the energy spectra in the fermion and boson sectors, respectively.

iii.2.1 Fermion spectrum

The fermion spectrum is given by, in terms of the wave vector defined on the reciprocal lattice,

 ^Hf=∑k,σ[fAkσfBkσ]†⎡⎢⎣λσ−μ−χvfηk−χvfη∗kλ¯σ−μ⎤⎥⎦[fAkσfBkσ], (16)

where is the dispersion due to the nearest neighbor hopping , which takes the form of

 ηk=exp(iky)+2cos(√3kx/2)exp(−iky/2)

on the honeycomb lattice and

 ηk=2(coskx+cosky)

on the square lattice. The sum over runs over the first Brillouin zone corresponding to a unit cell with two sites. The particle-hole symmetry at half-filling requires and , where becomes nonzero when the AF order develops. The fermion dispersion is thus obtained by diagonalizing Eq. (16),

 Ef±(k)=±√ε2+|χvfηk|2. (17)

A gap of would open in the fermion spectrum in the presence of AF order.

On the square lattice, the sublattices and become equivalent in the PM phase where the fermion spectrum in Eq. (16) simplifies to

 ^Hf=−χvf∑k,σηkf†kσfkσ. (18)

Here, the sum over runs over the first Brillouin zone corresponding to a unit cell containing only one site.

iii.2.2 Boson spectrum

The charged boson spectrum is governed by

 ^Hd=∑kΨ†kMkΨk,Ψk=[dAk,dBk,e†B¯k,e†A¯k]T,

where the boson Hamiltonian matrix

 Mk=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ϵvd+λ−χvdηk−Δvdηk0−χvdη∗kϵvd+λ0−Δvdη∗k−Δvdη∗k0ϵvd+λ−χvdη∗k0−Δvdηk−χvdηkϵvd+λ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (19)

Here the relations due to the particle-hole symmetry at half-filling have been applied. Note that is independent of spin. The boson dispersion is obtained by diagonalizing Eq. (19) using the standard boson Bogoliubov transformation:

 Ed±(k)=√(ϵvd+λ±|χvdηk|)2−|Δvdηk|2. (20)

Each branch of the above dispersion is doubly degenerate. The condition for a real physical dispersion requires that . When the equality is satisfied, the boson spectrum is gapless and the bosons can condense into the zero energy state. Otherwise, an energy gap

 Ξd=2√(ϵvd+λ−ζ|χvd|)2−(ζ|Δvd|)2

develops in the boson spectrum and the doublon and holon condensate would be depleted.

In the PM phase on the square lattice, the sublattices and become equivalent and the boson Hamiltonian matrix in Eq. (19) simplifies to

 ^Hd=∑k,σ[dke†¯k]†⎡⎣ϵvd+λ−χvdηk−Δ%vdηk−Δvdηkϵvd+λ−χ%vdηk⎤⎦[dke†¯k]. (21)

This results in a doubly degenerate boson dispersion relation

 Ed(k)=√(ϵvd+λ−χvdηk)2−(Δvdηk)2. (22)

iii.2.3 Self-consistency equations

The D-H binding saddle point solution can be obtained by solving the set of self-consistency equations derived from minimizing the energy with respect to the variables {, , , , , , , , }:

 χvd=2tgp0↑p0↓χf, (23) Δvd=gχf∑σ(tp20σ−gΔdχvfYσ), (24) ϵvd=−ζg2χfχvf∑σ(1−p20¯σ)Yσ, (25) p0σ=2ζtgχf(χdp0¯σ+Δdp0σ)λ−λσ−2ζg2χfχv% fY¯σ(1−nd−p20σ), (26) p20↑−p20↓=nf↑−nf↓, (27) 2nd+p20↑+p20↓=1, (28) χd=d20+12Nζ∑k′⟨ηkd†AkdBk+h.c.⟩, (29) Δd=d20+12Nζ∑k′⟨η∗kdAkeB¯k+h.c.⟩, (30) nd=d20+12N∑α={A,B}∑k′⟨dαkdαk⟩, (31)

where the fermion density and hopping per spin is readily obtained from the fermion spectrum in Eq. (16). It is instructive to examine the last three equations for the D/H hopping, the D-H binding, and the D/H density. The closing of the boson gap leads to a zero energy mode at whose occupation enables the single-boson condensate . This zero mode will be subsequently taken out of the momentum summations in Eqs. (29-31). Accordingly, the solutions to this set of self-consistency Eqs. (23-31) must be searched under two conditions: (i) assume and (ii) assume a nonzero . In the latter case, one more variable () is introduced together with one extra equation that ensures the existence of the zero energy mode,

 ϵvd+λ=ζ(|χvd|+|Δvd|). (32)

If multiple solutions are found, the one with the lowest energy should be chosen as the ground state. In practice, we found only one solution at any given .

iii.2.4 Electron spectral function and spectral density

Once the saddle point solution is obtained, the spectral function of the physical electrons can be calculated from the one-particle Green’s function . The detailed derivation of the spectral function and the integrated spectral function (ISF), i.e. the tunneling density of states,

 Nα(ω)=−Im∫β0dτeiωτ∑k,σGασ(k,τ), (33)

are given in the Appendix. Note that since the spectral function involves convolutions of the boson normal and the anomalous (due to pairing) Green’s functions with that of the -fermion, the single-particle energy gap for the physical electron is the sum of the fermion and boson gaps . More importantly, the coherent QP excitations would only emerge with the D/H condensate that recombines the charge and spin degrees of freedom, and can be detected by the QP coherent peaks in .

iii.3 Ground state wavefunctions

Before presenting the results on the honeycomb and the square lattice, it is instructive to discuss the possible phase structure in terms of the general form of the ground state wave function of the D-H binding saddle point. Since the Hilbert space is represented by those of the fermion and the slave bosons, the electron ground state wavefunction is a product of the ground state wave functions for the bosons and fermions,

 Ψ(→r1σ1,…,→rNσN)= ΨB(→r1,…,→rNd;→r1,…,→rNe) ⊗ΨF(→r1σ1,…,→rNσN). (34)

Here labels the spins of electrons, while and are the number of doublons and holons, respectively. From Eq. (16), it is clear that the fermion wavefunction is given by a Slater determinant, i.e. in both the PM and the AF ordered phases. Comparing to the conventional wavefunction form for an interacting many-body electron system, , the variational Jastrow factor has been promoted to full-fledged boson wavefunctions, thus allowing possible new electronic phases. The key physics in our theory is the boson inter-site correlations. The corresponding boson ground state wavefunction in second quantized form is thus a direct product of single-boson condensates and the pairing of uncondensed doublons and holons rokhsar94 .

On the square lattice, the boson wavefunction from the Hamiltonian in Eq. (21) thus takes the form

 |ΨB⟩=(e†0)N0e(d†0)N0d|0⟩B⊗∏kexp(−gkd†ke†¯k)|0⟩B, (35)

where is the vacuum of the boson sector, and the pairing function . The boson wavefunction on the honeycomb lattice described by the Hamiltonian in Eq. (19) has a similar, slightly more complicated form due to the doubling of the unit cell. In the first part of the boson wavefunction, and are the numbers of condensed holons and doublons that are determined by the condensate density and . This is the only part retained in the KR saddle point or the Gutzwiller approximation. Together with , they describe the coherent QP part in the excitation spectrum of a correlated Fermi fluid. The second part is due to D-H binding which, together with , describes the incoherent excitations and plays an integral part in the Mott transition. So long as the condensate part is present, the charge and spin degrees of freedom combine such that the elementary excitations carry the quantum numbers of an electron and appear as QP poles in the single-particle Green’s function. Thus, the coexistence of the condensate and the binding parts heralds the coherent QP and the incoherent Hubbard excitations in Mott-Hubbard systems before the Mott transition. This is the case in both the PM () and the Slater AF phase () where in the condensed part of the boson wavefunction while changes from a PM to AF Slater determinant. However, as we will show, when , the quantum fluctuations due to D-H binding destroy the single-particle condensate, i.e. , as all doublons and holons are bound together, giving rise to a charge gap. The boson wave function is given entirely by the second part in Eq. (35). Interestingly, this boson wave function is just the wave function of a resonating valence bond (liquid) state, which in the present context can also be understood as that of an excitonic insulator, since the doublon and holon carry opposite charges. Because all the doublons are bound to the holons, the elementary excitations do not carry the quantum numbers of an electron and the entire single-particle excitations are incoherent as the charge and spin cannot recombine to form a coherent quasiparticle. Had this Mott transition taken place before the AF order, this insulating phase would be a spin liquid. However, as we will see that on the bipartite lattices without frustration, AF order happens before the Mott transition, i.e. . We thus term the phase for the AF* phase, which is indeed an example of spin-charge separation above one-dimention, albeit taking place inside the AF ordered phase. Note that what distinguishes the AF* phase from the Slater AF insulator is the complete depletion of the single-particle condensates of the holons and doublons above such that all doublons are bound to holons.

iii.4 Mott, Slater AF, and AF Mott transitions

iii.4.1 Results on the honeycomb lattice

The D-H binding saddle point solutions on the half-filled honeycomb lattice are summarized in Figs. 4 and 5. The variational parameters and , the order parameter for the D/H hopping and D-H pairing are plotted in Figs. 4a and 4b as function of the Hubbard . At , , thus all doublons and holons are single-particle condensed with such that , recovering the noninteracting limit. The SM phase remains stable at small . With increasing , the doublon density decreases as shown in Fig. 4c. Due to the increase in D-H binding, the D/H condensate decays faster than in the KR saddle point solution shown in Fig. 2a. To study the Mott transition, we first restrict the solution to be in the PM phase by enforcing , which amounts to suppressing possible magnetically ordered states. As shown in Fig. 4c, the Mott transition takes place at , which is considerably smaller than that of for the BR transition (Fig. 2a). The condensate vanishes and all doublons are bound with the holons in the Mott insulating phase for , accompanied by the opening of a charge gap that is linear in (Fig. 4d). The ISF of the physical electrons is shown in Fig. 5a. Notice the transfer of the coherent QP weight to the incoherent part with increasing and the complete suppression of the coherent QPs in favor of two broad incoherent spectral features beyond the Mott transition that originate from the bosonic excitations in Eq. (19) to be discussed later. Since the -fermion spinon remains gapless, the insulating phase is a gapless SL. Thus, we find no evidence on the honeycomb lattice for the proposed gapped SL phase meng .

Next, we allow magnetism and study the interplay between AF order, D-H binding, and the Mott transition in the ground state. Fig. 4c shows that the SM phase on the honeycomb lattice remains stable until a critical , where the staggered magnetization () onsets. We find that for , where , although a single-particle gap opens in the fermion sector (Fig. 4d), the zero energy mode remains stable in the - sector and continues to support the D/H condensate. Thus, the spin and charge continues to recombine in this regime and there are coherent excitations corresponding to the sharp QP peaks in the ISF shown in Fig. 5b at and . This phase is thus an Slater AF insulator whose wavefunction would overlap well with an AF Slater determinant.

Remarkably, a Mott transition in the presence of AF order takes place at . For , an AF Mott phase (i.e., the AF phase) emerges with the opening of the boson gap in the - sector (Fig. 4d) as the D/H condensate vanishes. Since all doublons are bound to holons, the charge and spin cannot recombine and the electrons are thus fractionalized in the AF phase. A direct consequence for the lack of elemental excitations carrying the electron quantum number is that the lack of coherent QP peaks in an entirely incoherent excitation spectrum, as can be seen from the broad ISF at in Fig. 5b at and . Unlike in the Slater AF phase, the vanishing of the D/H condensate in the AF phase, enables the deconfinement the spin and charge degrees of freedom, such that the ground wavefunction has no overlap with Slater determinant-like states. The excitation energy gap for the physical electron, , exhibits a derivative discontinuity at (Fig. 4d) due to the opening of the Mott gap in the AF phase. However, the magnetization remains analytic across the AFAF transition in Fig. 4c, which is a topological confinement-deconfinement transition associated with the Ising-like global symmetry (, ) that is broken in the Slater AF phase by the D/H condensate and restored in the AF phase.

The continuous SM to AF transition at compare well to the most recent QMC calculations on large system sizes by Sorella et. al. sorella that finds the onset of AF order and a single-particle excitation gap at . Since we have not included the inter-site spin fluctuations described by the dynamics of the -boson, our magnetic gap is somewhat larger than the QMC values, and we will not attempt quantitative comparisons to results obtained by other numerical methods such as cluster dynamical mean-field theory (CDMFT) calculations with continuous-time QMC (CTQMC) or exact diagonalization (ED) impurity solvers. While the CTQMC-CDMFT wliu ; hur12 is performed at relative high temperatures and not very suitable for extracting small energy gaps in the quantum states, the ED-CDMFT liebsch ; zylu as well as the ED-VCA (variational cluster approximation) jxli ; seki revealed spurious excitation gaps at very small , before the emergence of AF order. This was viewed as supporting evidence for the proposed gapped SL phase meng . Recently, Hassan and Sénéchal hassan13 noticed that these ED-CDMFT and ED-VCA calculations use only a single bath orbital per cluster site which they argued is insufficient and leads to the artificial excitation energy gaps for all nonzero values of . Their calculations with two bath orbitals connecting each cluster site show that the PM Mott transition and thus the SL phase is indeed preempted by a magnetic transition occurring at a lower value of . Liebsch and Wuliebsch2 further pointed out that the spurious excitation gap at very small originates from the breaking of the translation symmetry in these cluster calculations.

iii.4.2 Results on the square lattice

Figs. 6 and 7 show the results obtained on the half-filled square lattice, which are qualitatively the same as those obtained on the honeycomb lattice. In the PM phase, the Mott transition is at in contrast to the Brinkman-Rice transition at without taking into account of D-H binding. Because the band structure leads to perfectly nested Fermi surface at half-filling, the PM metallic phase is unstable towards AF order for infinitesimal on the square lattice. The AF order therefore emerges at with an exponentially small staggered magnetization, as shown in Fig. 6c. This Slater AF insulator with coherent QP excitations is stable until where a transition into the AF* phase with the vanishing of holon/doublon condensate and the opening of the charge gap (Fig. 6c,d) and the disappearance of coherent QP peaks in favor of two broad incoherent features in the integrated spectral density Fig. (7). It is important to note that the Fermi level density of states, whose two limiting behaviors, vanishing or divergent, are presented by the unfrustrated honeycomb and square lattices respectively, while affecting the PM to Slater AF transition, does not play an essential role in determining the Mott transition from the Slater AF into the AF* phase. This is because the latter is tied to the opening of the charge gap in the boson sector near the doublon and holon band bottom above a finite , as can be seen from the boson dispersions shown in Fig. 9.

To further explore the generality of these predictions, we have studied the case where the noninteracting band has a semicircular density of state where is the bandwidth. We found an identical phase structure with a PM metal to Slater AF insulator at , followed by the Mott transition into that AF* phase at . Note that although the semicircle density of states can be realized in the infinite-d unfrustrated Bethe lattice, these results should not be considered as obtained for the infinite-d Hubbard model, since taking the infinite-d limit would suppress all inter-site correlations, including the inter-site D-H binding considered here. Thus, in the infinite-d limit, we would only recover the KR saddle point solution, i.e. the Gutzwiller approximation which is exact for the Gutzwiller wave function approach in infinite dimensions GWinf-d . In this sense, our approach can be viewed as going beyond the Gutzwiller approximation by including the inter-site correlations in physical dimensions.

With this difference in mind, we proceed to compare in Fig. 8 the local spectral function in the PM phase on the square lattice of our D-H binding theory with the results obtained from the single-site DMFT on the square lattice zitko09 at different Hubbard . Overall, we find remarkable agreement in the incoherent part of the spectrum for all values of on both sides of the Mott transition. The most significant deviations of the results are in the low energy QP peaks around the Mott transition. The main cause of the latter can be traced to a particular property of the single-site DMFT formulation that the spectral density at is independent of in the infinite-d limit of the Hubbard model muller-hartmann , which holds the QP peak at constant height until its width goes to zero at the Mott transition, seen from the DMFT data in Fig. 8. The latter is no longer true in the presence of inter-site correlations beyond the infinite-d limit, as shown in the D-H binding results in Fig. 8, where the QP peaks, agreeing with the DFMT result for moderate , are suppressed in both height and width and disappears completely at the Mott transition whose critical value is significantly reduced by intersite correlations park08 . We also find broad qualitative agreement with the local spectral density of cluster DMFT calculations that captures certain short-range correlations, although finding consistency in CDMFT results is difficult due to the different cluster embedding procedure (including cluster shapes and sizes) and the choice of the impurity solver. We find that it is particularly intriguing that in the CDMFT study of the PM phase by Park, et. al. park08 , the metallic phase has an ISF consistent with our result in the PM phase; whereas on the insulating side near the metal-insulator transition, the local spectral function displayed a small gap with very pronounced peaks at the gap edge that closely resemble our findings in the AF Slater insulator. Indeed, these peaks at the edge of the magnetic gap are clear hallmark of the coherence QP peaks characteristic of a Slater spin density wave insulator. We thus conjecture that the CDMFT findings of a small gap PM insulating state with pronounced gap-edge coherence peaks are due to the fluctuating or short-range Slater AF order, and with increasing , a true Mott transition would emerge with the suppression of the QP peaks and the opening of the charge gap.

iii.4.3 Stability of the D-H binding saddle point

Next we comment on the D-H binding saddle point stability with respect to gauge field fluctuations. It is known that KR formulation introduces three gauge fieldsrasul88 since the action is invariant under: , , , , and ,. The condensate breaks two of the symmetries and turns the gauge fields associated with massive by the Anderson-Higgs mechanism. The remaining symmetry is also broken in the SM and the AF phase by the D/H condensate, making the -gauge field massive. In the AF phase, it is the D-H pairing that breaks the symmetry and the -gauge field remains massive, as does its staggered component due to the D/H hopping fields . The absence of gapless gauge field fluctuations supports the stability of the obtained phases.

iii.5 Nature of incoherent Mott-Hubbard excitations

It is enlightening to discuss the energy spectrum and the spectral function of the doublons and holons in connection to the nature of the incoherent Mott-Hubbard excitations in the local spectral function. The dispersion of the holons and doublons in Eqs. (20) and (22) and the corresponding density of states are shown in Fig. 9 in the PM phases of the honeycomb and square lattice Hubbard models, respectively, at several values of across the Mott transition. Their behaviors are similar in the AF and AF phases.

Honeycomb Lattice: At a fixed , the boson spectrum shows two doubly-degenerate dispersive branches given in Eq. (20). There are several noteworthy features. (i) Both dispersive branches are flat near the point of the hexagonal Brillouin zone, leading to the two Van Hove singularities in the boson ISF plotted in the right panel of Fig. 9a. (ii) The two branches cross and produce the Dirac cone at the high symmetry point at a finite energy that increases with increasing , leading to the V-shaped density of states. Remarkably, (i) and (ii) combine to form a Dirac-cone like dispersion that is similar and can be regarded as a “ghost” band of the bare electron dispersion carried by the excitations of the D/H complex. This property was pointed out in the systematic large- expansion study of the - model for doped Mott insulators largeN . It is remarkable that the ghost Dirac-cone feature manifests itself in the broad peak-dip-peak structure in incoherent part of the ISF for the physical electrons shown in Fig. 5a, which can now be identified as the D-H excitations. (iii) The low energy properties of the boson dispersion near the -point is also intriguing. For , i.e. on the metallic side of the Mott transition, the lower energy branch is gapless, i.e. and disperses linearly away from the point. The existence of the zero-energy mode, together with the vanishing of the ISF , enables the finite-temperature D/H condensation in a two-dimensional system such that at zero temperature. In contrast, for , or when magnetism is allowed, a charge gap opens up at the point, indicating the emergence of the Mott insulating or the AF phase with complete suppression of the D/H condensate. Note also that the gapped is parabolic near , giving rise to a finite ISF at the band bottom.

Square Lattice: The above findings on the connection between the incoherent Mott-Hubbard excitations and the doublon-holon spectrum applies in a straightforward manner to the square lattice case as well. In contrast to the intrinsic two-sublattice structure of the honeycomb lattice, the boson spectrum on the square lattice has only one branch given in Eq. (22), which is shown in Fig. 9b. The corresponding density of states has a single Van Hove peak tied to the dispersion near the X point. Similar to the honeycomb lattice case, this branch of D-H excitations manifest itself in the single broad peak on the particle and the hole side of the electron local spectral function shown in Fig. 7.

The change in the bosonic dispersion across the AF Mott transition, i.e. from the Slater AF insulator to the AF* phase, is qualitatively the same as those displayed in Fig. 9 for the PM Mott transition. Namely, the boson spectrum is gapless with linear dispersion supported by the D/H single-particle condensate on the Slater AF side and develops an energy gap, when all holons are bound to doublons, above which a quadratic dispersion is found for the Bogoliubov quasiparticles. The linear dispersion in the PM phase is a bosonic representation of the collective zero-sound excitations in the Landau Fermi liquid. In the Slater AF phase, despite the opening of a single-particle magnetic gap, the absence of a charge gap is reflected in the existence of gapless collective excitations represented by the bosons. It is only after entering the AF* phase, the charge excitations are gapped out, leaving only the spin-waves as the low energy excitations inside the single-particle energy gap.

Iv Conclusions

In summary, we have shown that the binding between doublons and holons plays an essential role in describing the incoherent excitations and the Mott transition in strongly correlated Mott-Hubbard systems. For the honeycomb lattice Hubbard model, we showed that the SM to AF Slater insulator transition is followed by a Mott transition into a fractionalized AF phase with increasing . Interestingly, a different AF phase of a fractionalized antiferromagnet was proposed in the effective Z gauge theory description of doped Mott insulators in the projected () Hilbert space where spinons are paired into a Néel state and doublons are absent fishersenthil . In contrast, the incoherent charge excitations through D-H binding is essential in the AF phase proposed here, which is more inline with the importance of doublons in describing Mottness emphasized recently phillips10 .

The most practical way to test our predictions is to measure the energy gap for single-particle excitations using spectroscopic probes such as ARPES and STM. Our theory shows that with increasing (which can be varied experimentally by applying pressure or isoelectronic chemical substitution), the system goes from a gapless PM state to an AF insulator where the single-particle gap is controlled by the magnetic gap, followed by a transition to the AF* phase where a charge gap opens and adds on top of the magnetic gap. Thus, there is a singularity (kink) in the evolution of the gap as a function of . Perhaps even more directly, the single-particle spectral function as measure by ARPES shows well defined QP peaks above the magnetic gap in the Slater AF phase, but would exhibit no coherent excitations in the AF* phase. Such an AF phase on the square lattice may have been observed in the parent AF insulating compound of the high- cuprates by ARPES experiments arpesundoped , which find no coherent QP excitations at all energies.

As a concrete example, we propose to revisit the time-honored Mott-Hubbard system, i.e. the transition metal oxide VO, under chemical pressure achieved by Cr or Ti substitutions. In this case, a finite temperature Mott metal-insulator transition above the low-temperature AF insulating ground state has been well established as a function of chemical substitution kuwamoto ; mcwhan71 ; mcwhan73 . Our theory predicts that hidden inside the AF insulating ground state is a transition from a Slater AF to the AF* phase. Moreover, melting the AF order in the Slater AF insulator would result in the metallic state, whereas melting the AF* phase at higher Cr substitutions gives rise to the Mott insulator at finite temperatures. Performing the experiments described above in these materials would either provide support or disprove our theory.

We thank F. Wang and Y. Yu for useful discussions. This work is supported by DOE grants DE-FG02-99ER45747 and DE-SC0002554, and the Thousand Youth Talents Plan of China (SZ). ZW thanks Aspen Center for Physics for hospitality.

Appendix A ISF of physical electrons

The ISF, which equals the tunneling density of states (DOS), for the physical electrons is given by

 Nα(ω)=−∑k,σIm∫β0dτeiωτGασ(k,τ),

where retarded single-particle Green’s function manhan

 Gασ(k,τ)=−⟨Tτcαkσ(τ)c†αkσ(0)⟩,

with the sublattice index. In the KR slave boson formulation kr , the electron operator is composed of . Within our saddle point solution, and are approximated by their saddle point average for the local Green’s functions. The electron operator in momentum space is thus given by,

 cαkσ=rασ∑q,q′(e†α¯qpαq′σ+p†α¯q′¯σdq)fα,k−q−q′,σ,

with the normalization factor

 rασ=⟨^Lασ^Rασ⟩=[(1−nd−p2α0σ)(1−ne−p2α0¯σ)]−1/2.

Therefore, the electron Green’s function

 Gασ(k,τ)=r2ασ∑q