Charge dynamics of the antiferromagnetically ordered Mott insulator

We introduce a slave-fermion formulation in which to study the charge dynamics of the half-filled Hubbard model on the square lattice. In this description, the charge degrees of freedom are represented by fermionic holons and doublons and the Mott-insulating characteristics of the ground state are the consequence of holon–doublon bound-state formation. The bosonic spin degrees of freedom are described by the antiferromagnetic Heisenberg model, yielding long-ranged (Néel) magnetic order at zero temperature. Within this framework and in the self-consistent Born approximation, we perform systematic calculations of the average double occupancy, the electronic density of states, the spectral function and the optical conductivity. Qualitatively, our method reproduces the lower and upper Hubbard bands, the spectral-weight transfer into a coherent quasiparticle band at their lower edges and the renormalisation of the Mott gap, which is associated with holon–doublon binding, due to the interactions of both quasiparticle species with the magnons. The zeros of the Green function at the chemical potential give the Luttinger volume, the poles of the self-energy reflect the underlying quasiparticle dispersion with a spin-renormalised hopping parameter and the optical gap is directly related to the Mott gap. Quantitatively, the square-lattice Hubbard model is one of the best-characterised problems in correlated condensed matter and many numerical calculations, all with different strengths and weaknesses, exist with which to benchmark our approach. From the semi-quantitative accuracy of our results for all but the weakest interaction strengths, we conclude that a self-consistent treatment of the spin-fluctuation effects on the charge degrees of freedom captures all the essential physics of the antiferromagnetic Mott–Hubbard insulator. We remark in addition that an analytical approximation with these properties serves a vital function in developing a full understanding of the fundamental physics of the Mott state, both in the antiferromagnetic insulator and at finite temperatures and dopings.


Introduction
The mechanism underlying the Mott metal-insulator transition [1] stands as a fundamental theoretical challenge in condensed matter physics.In 1937, de Boer and Verwey [2] reported that a class of transition-metal oxides with partially filled bands, specifically NiO and MnO, are semiconductors or insulators in direct contradiction to predictions by conventional band theory.This motivated Mott and Peierls [3] to point out the importance of the electrostatic interaction between the electrons, and Mott later introduced the concept of the metal-insulator transition that bears his name [4,5] to describe insulating behaviour arising as a result of strong electron-electron correlations.The discovery of high-T c superconductivity in a class of doped antiferromagnetic Mott insulators [6] revived an enormous and lasting interest in understanding the Mott phase and the associated metal-insulator transition.
The Hubbard model [7] is the minimal model describing the competition between the kinetic energy of the electrons and their on-site Coulomb interaction.It captures many characteristic features of strongly correlated systems and thus serves as a paradigm for numerous phenomena in condensed matter physics.It is believed that the Hubbard model contains all the basic physics of the Mott metal-insulator transition and, in some quarters, that it may reveal the mechanism of high-T c superconductivity.However, despite the simplicity of the Hubbard model, exact results can be obtained only from the Bethe Ansatz [8] in one dimension and from Dynamical Mean-Field Theory (DMFT) [9,10] in infinite dimensions.
There exist many proposals for the primary mechanism driving the Mott transition.Hubbard's equation-of-motion methods [7,[11][12][13][14] attribute the charge gap to the formation of the incoherent lower and upper Hubbard bands.These provided the first example for a metal-insulator transition in which the insulating behaviour is not accompanied by the onset of magnetic order.Brinkman and Rice [15] applied the Gutzwiller variational method [16] to treat the metal-insulator transition out of the Fermi-liquid metallic phase, and ascribed the transition to the vanishing of the quasiparticle residue, Z, and the divergence of the quasiparticle effective mass, m * .The Hubbard approximation captures the incoherent part of the physics while the Brinkman-Rice approximation captures the coherent part.However, neither approximation takes the effect of spin fluctuations into account.
For the half-filled single-band Hubbard model on the square lattice, quantum Monte Carlo simulations [17,18] have shown that the ground state is an antiferromagnetic insulator, although by the Mermin-Wagner theorem its Néel temperature is zero.In the weak-coupling limit, Fermi-surface nesting and the proximity to a van Hove singularity in the density of states act to induce a spin-density-wave state and thus to produce a gap [17].An asymptotically exact weak-coupling solution for the Hubbard model was given in reference [19].In the strong-coupling regime, it is the large on-site Coulomb repulsion energy, U, for double site occupancy that suppresses electron mobility and determines the Mott gap.
For the intermediate-coupling regime, where no well-controlled theoretical solution exists, many numerical methods have been applied to the two-dimensional (2D) Hubbard model, including exact diagonalisation [20][21][22][23][24], quantum Monte Carlo [25][26][27][28][29][30][31][32], cluster perturbation theory [33][34][35][36][37][38], the variational cluster approximation [39][40][41][42][43] and cluster DMFT [44][45][46][47][48][49].A detailed review, including further results from density-matrix renormalisation-group (DMRG) calculations, may be found in reference [50].However, all of these methods suffer from different intrinsic limitations.Cluster perturbation theory provides an approximate lattice Green function for a continuous wave-vector space but is not self-consistent and cannot describe broken-symmetry states, which are known to be present for the half-filled square lattice.The variational cluster approximation can be viewed as an extension of cluster perturbation theory, which allows for broken symmetries by introducing Weiss fields, but remains limited by the cluster size.In cluster DMFT, the quantum impurity model can be solved by quantum Monte Carlo or exact diagonalisation.The former operates at finite temperature and imaginary time, requiring extrapolation to recover zero-temperature information and analytic continuation methods to obtain real-frequency results, neither of which is well controlled; further, the ubiquitous fermion sign problem affecting quantum Monte Carlo methods becomes severe when the system is doped.The latter is implemented at zero temperature and gives direct real-frequency dynamical information, but can access only small cluster sizes.DMRG is inherently 1D in nature and can be applied only on a narrow cylinder; the ongoing development of higher-dimensional analogues based on tensor-network states has progressed to the point where an infinite projected entangledpair state (iPEPS) method has been used very recently to obtain very competitive ground-state energies [51].Anderson [52] has argued that the half-filled 2D Hubbard model is fundamentally nonperturbative in nature, in the same way as the 1D case, with a Mott gap present for all U > 0 and robust against temperature.Thus despite all of the theoretical and numerical progress made to date, the nature of the Mott gap at half-filling and the properties of the 2D Hubbard model remain as challenging open questions.
In the strong-coupling limit, below half-filling the dominant on-site Coulomb repulsion implies the absolute exclusion of doubly occupied sites.In this case, the Hubbard model can be mapped to the t-J model at the level of second-order perturbation theory [53].At half-filling, no empty sites remain and this model reduces to the antiferromagnetic Heisenberg model with only spin-fluctuation degrees of freedom.As U decreases, charge fluctuations play an increasingly important role in the Hubbard model [54,55], and in the half-filled case the elementary charge excitations are holons (empty sites) and doublons (doubly occupied sites), in equal numbers.However, even at weak coupling, insulating behavior remains guaranteed if the holons and doublons have a tendency to form bound states, which results in the presence of a charge gap.Variational Monte Carlo results [56][57][58][59][60] have shown that a variational wave function including holon-doublon binding effects can lower the ground-state energy and that the Mott transition can be characterised as an unbinding transition of holons and doublons.Several theoretical proposals for the mechanism of Mott physics contain holon-doublon binding as an important element, including the "hidden charge-2e boson" mechanism [61][62][63], the reconstruction of poles and zeros of the Green function [64][65][66], composite fermion theory [67] and the Kotliar-Ruckenstein slave-boson theory [68].The zeros of the Green function at the chemical potential in momentum space can be taken to define the "Luttinger surface," which is closely connected to the non-interacting Fermi surface [69].
The motion of a single hole in an antiferromagnetic background has been studied extensively within the self-consistent Born approximation (SCBA) [70][71][72][73][74][75], where the neglect of Feynman diagrams with crossing propagators is equivalent to neglecting the distortion of the spin background caused by the presence of the hole.This formulation is similar to the retraceable path (Brinkman-Rice) approximation [76] and the resulting single-hole spectral function is composed of two components, a sharp peak corresponding to coherent quasiparticle motion and an incoherent background.The coherent peak arises from the coupling between the hole and the spin excitations.Recent experiments have shown that the SCBA yields excellent agreement with resonant inelastic X-ray scattering measurements performed on the quasi-2D spin-1/2 antiferromagnet Sr 2 IrO 4 [77].
To make a meaningful contribution to such a complex and deeply studied problem, here we aim to provide an analytical framework which captures the essential features of the charge dynamics of the single-band Hubbard model.To introduce this framework, we restrict our considerations to the square-lattice model with only nearest-neighbour coupling, to half-filling, and to zero temperature.While it is unrealistic to expect to find much new physics in this most generic situation, our goals are to demonstrate the qualitative power of a suitably chosen mean-field description, to establish the semiquantitative accuracy of our results by benchmarking against the plethora of available numerical studies, and to lay a foundation for development in the more experimentally relevant directions of finite temperatures, extended bandstructures and finite dopings.
An accurate description of the Mott insulator is based on the strong-coupling limit, where its properties are robust.As we discuss in more detail below, this leads to a formalism where the charge degrees of freedom are represented by fermionic holons and doublons and the Mott-insulating state involves the formation of holon-doublon bound states.The spin degrees of freedom, represented by bosonic magnons, order magnetically at zero temperature for any finite interaction strength, but their quantum fluctuations act to renormalise the charge sector.Thus the task at hand is to consider the antiferromagnetically ordered Mott insulator and to describe accurately both the holon-doublon binding process and the spin renormalisation of the charge dynamics.
For specificity, we declare here that we take the on-site interaction to be the origin of Mott physics for all dimensionalities, temperatures or dopings, and antiferromagnetic fluctuations to be one consequence.It is true that this assumption remains unproven for the square lattice (2D) with nearest-neighbour hopping at half-filling: in this somewhat pathological case, the perfect Fermi-surface nesting means that a gap is opened in the charge sector at any finite interaction strength, and in a weak-coupling picture this can be interpreted as a process driven by the onset of antiferromagnetic order, i.e. not by the charge sector but by the spin sector.This result has led to significant confusion over whether an antiferromagnetic insulator can exist independently of a Mott insulator.We use the fact that there is no transition at any finite interaction strength in the model at hand to deduce that the two possible states are different manifestations of the same physics and are connected by a crossover.For practical purposes, here we take the extensive numerical calculations on the half-filled Hubbard model to indicate that the ground state is a magnetically ordered Mott insulator for all intermediate (and experimentally relevant) interaction strengths.The Mott (charge excitation) gap in our framework is a consequence of holon-doublon binding and its presence ensures a finite spectral (single-particle excitation) gap, while the spin excitations are gapless.The same holon-doublon binding mechanism is equally applicable to the Hubbard model at finite temperature or doping, where the spin sector is present only as short-range fluctuations.
Our approach is based on a slave-particle formalism [78][79][80][81][82] in which electron operators are expressed as a combination of "slave" fermionic and bosonic operators preserving the net fermionic statistics, and spin-charge separation is assumed.A degree of arbitrariness exists in ascribing the fermionic statistics to the spin (known as the slaveboson approach) or to the charge degrees of freedom (slave-fermion decomposition).Keeping the importance of spin fluctuations at the forefront of our considerations, we assume that the ground state has Néel order, which is known in the large-U limit, and that the spin degrees of freedom are described by the antiferromagnetic Heisenberg model.The ground state of this model on the square lattice for S = 1/2 is better described by the Schwinger boson (slave-fermion) formulation [81,84].Of equal importance, for an investigation of the holon-doublon binding mechanism it is physically much more intuitive to ascribe the fermionic statistics to the charge degrees of freedom, in analogy with the electron binding mechanism of Bardeen-Cooper-Schrieffer (BCS) superconductivity.
In this slave-fermion framework, antiferromagnetic long-range order corresponds to the condensation of one of the slave bosons on each sublattice in the ground state.On this basis, we treat the spin dynamics within the linear spin-wave approximation [53], where the elementary excitations are magnons.Because the motion of holons and doublons distorts the antiferromagnetic background even at half-filling, a consistent account of spin-fluctuation effects is of key importance in describing the charge dynamics.We treat the interactions among holons, doublons and magnons within the SCBA to calculate important physical quantities including the double occupancy, the spectral function, the electronic density of states, the quasiparticle Green function and the optical conductivity.Our results show a non-zero double occupancy for any finite U, that the Mott-insulating state results from holon-doublon binding, that spin fluctuations modify the size of the Mott gap and that this gap can be probed accurately by measurements of the optical conductivity.
This paper is organised as follows.In section 2 we introduce formally the model and the methods we use to perform our calculations.In section 3 we compute the doublon density for all intermediate values of U and in section 4 we present the spectral function to discuss the coherent and incoherent components of the charge response, the density of states and the Mott gap.In section 5 we calculate the electron Green function and associated Luttinger surface, and deduce the effective quasiparticle bandstructure.Section 6 contains our results for and conclusions from the optical conductivity and a summary is provided in section 7.

Model and Method
The single-band Hubbard model is defined by the Hamiltonian where c iσ (c † iσ ) denotes the annihilation (creation) operator of an electron with spin σ on site i of the square lattice and i, j indicates that we restrict the hopping terms to nearest-neighbour sites i and j only.We set the hopping parameter as t = 1 to establish the energy units of our calculations and U represents the on-site Coulomb repulsion.
In the slave-fermion formalism, the electron operator is written as where d i and e i are fermionic operators denoting the charge degrees of freedom and s iσ are bosonic operators describing the spin degrees of freedom, with σ = 1 for spin ↑ and −1 for spin ↓.The operator e † i creates an empty (unoccupied) site, a holon, at lattice point i, d † i creates a doubly occupied site, a doublon, and s † iσ (s † iσ ) is the creation operator for a singly occupied site i with spin σ (σ).In this formulation, the local Hilbert space is enlarged and the constraint should be satisfied to eliminate unphysical states.Substituting equation (2) into equation (1) gives the form where δ denotes the lattice vectors (a, 0) and (0, a).Here we do not include a Lagrangemultiplier term to enforce the local constraint at a global level (introducing a chemical potential), but instead implement equation (3) as a self-consistent condition in our treatment.This procedure has the same effect, in that the constraint is satisfied only on average, which is a primary shortcoming of all analytical approaches to locally constrained problems unless gauge fluctuations can be included.Otherwise, the efficacy of this type of approximation is difficult to assess by any means other than comparing the results it yields with numerical calculations where the constraint can be enforced exactly.
From the first line of equation ( 4), holons and doublons can hop between nearestneighbour sites with accompanying creation and annihilation of singly-occupied states.This term serves as the starting point for applying the SCBA, which provides a proper treatment of the coupling between the charge degrees of freedom and the spin fluctuations [70][71][72][73][74][75].This procedure requires the assumption of a Néel-ordered ground state, which is equivalent to a Bose condensation of the s † iσ operators and is discussed in detail below.The second line of equation (4) shows that the kinetic term of the Hubbard model contains a pairing interaction between holons and doublons in the slave-fermion representation.As noted in section 1, the analogy with BCS theory motivates both the attribution of fermionic statistics to the charge degrees of freedom and the formation of bosonic bound pairs of fermionic holons and doublons as the process underlying the opening of a charge gap and contributing to the charge dynamics of the Mott insulator.In the last line of equation ( 4), the Coulomb repulsion, U, appears as a mass term for holons and doublons, which gain a dispersive nature both through their pairing and through their interactions with the magnetic background, as described by the SCBA.
As appropriate for a method based on the large-U limit, we assume that the ground state for the spin degrees of freedom is the Néel antiferromagnet, also for finite U, and that their fluctuations are described by the antiferromagnetic Heisenberg Hamiltonian, with the coupling constant taken for simplicity as J = 4t 2 /U.In the Schwinger boson representation, the spin operators are given by where σ = (σ x , σ y , σ z ) denotes the Pauli matrices.The Néel antiferromagnetic state corresponds to a Bose-Einstein condensation of one of the two types of bosonic operator on each sublattice [83], which we describe by the uniform mean-field assumption for the two sublattices A and B. The constraint (3) is then recast as Because there is only one Schwinger boson on each sublattice, henceforth we simplify the notation s iσ to s i .The revised form (8) of the constraint is the self-consistency condition in our calculation.We include the spin excitations of the antiferromagnetic Heisenberg model using the straightforward linear spin-wave approximation; although this method is only valid strictly for large spin and high-dimensional lattice geometries, it has been shown [84] that physical quantities, including the ground-state energy and sublattice magnetisation, obtained for the S = 1/2 square-lattice model within this approximation are qualitatively correct and readily renormalised to the values given by numerical calculations.A more detailed treatment of the Heisenberg model is provided in Appendix A.1.Within this approximation, equation ( 4) can be expressed as where ) is the Nambu spinor, which contains the charge degrees of freedom.The explicit forms of εk and M(k, q) may be found in Appendix A.2 along with full details of the SCBA for the Hubbard model in this form.The first term of equation ( 9) describes the unperturbed charge dynamics, with holon-doublon binding appearing in the off-diagonal part of the matrix.The second term describes the interaction between the charge and spin degrees of freedom.We define the full charge Matsubara Green function as and calculate this at the level of the SCBA.This process includes the coupling between the charge and spin dynamics and is equivalent to the series of Feynman diagrams shown in figure 1, where F (0) and F are respectively the bare and interacting charge Green functions and D is the magnon Green function, whose explicit form is given in Appendix A.1.The SCBA has been used to calculate the motion of a single hole in an antiferromagnetic background in the t-J model, where a consistent account of the mutual effects of charge motion and spin fluctuations is similarly essential.The results of these studies show that the coupling between the holon and the spin waves induces a quasiparticle-type response often labelled a spin polaron [70][71][72][73][74][75].Although a proper treatment of charge and spin fluctuations can be obtained in this way, we comment that the approximation does not include vertex corrections.The self-consistent Dyson equation for the charge Green function is calculated by standard techniques, which yield The solid, red line is our result, calculated for a system size of 48×48 and with broadening parameter η = 0.08.For comparison, the dashed lines show analogous results obtained from cluster perturbation theory (CPT, blue, up-pointing triangles) [86], quantum Monte Carlo (QMC) for a 4 × 4 lattice at inverse temperature β = 16 (green, down-pointing triangles) [18], variational Monte Carlo with different trial wave functions [|ψ pow (VMCPow, maroon diamonds) and |ψ B (VMCB, magenta squares) in the notation of reference [60]] and DMRG (solid, black line) with extrapolation to the thermodynamic limit [50].
where the self-energy Σ(k, iω n ) is given in Appendix A.2.The retarded charge Green function can be obtained from the Matsubara Green function (10) by analytic continuation through iω n → ω + iη.The effect of the parameter η is to provide a finite width to peaks in the density of states.

Double occupancy
We begin the discussion of our results for the Mott-insulating state of the square-lattice Hubbard model by considering the averaged double-occupancy parameter, D. For the half-filled system, the average number of doubly occupied sites is equal to the number of empty ones, and is given by D is exactly zero only when the on-site Coulomb repulsion, U (1), is infinite, but charge fluctuations are intrinsic to the Hubbard model for all finite U values, making D finite.It therefore reflects average information about the effects of charge fluctuations and as such can be used to characterise the Mott transition [15,54,55,[58][59][60]85]. The dependence of D on U at zero temperature is shown in figure 2, where our results are compared with calculations by quantum Monte Carlo [18], variational Monte Carlo [60], cluster perturbation theory [86] and DMRG [50].At the qualitative level, all results lie in the same general range of values and, with the exception of one VMC approach, show a broadly similar functional form.Quantitatively, it is immediately clear that the different numerical results differ from each other quite significantly throughout the regions of weak and intermediate coupling, which may be taken as a signal of how difficult the Hubbard problem is in this regime.We stress that finite-size extrapolation of our results (not shown) confirms that our 48 × 48 calculations are fully representative of the infinite system, whereas none of the numerical data shown in figure 2 have been extrapolated to the thermodynamic limit other than the DMRG results.However, extrapolated results obtained by a range of methods may be found in a very recent review [50].
Only in the strong-coupling regime do the differences in D values obtained from all of these methods become small.Not only is our result in full agreement here, but because it is constructed in the strong-coupling limit, it can be expected to give the correct asymptotic form of D as U becomes large.We have calculated D in the SCBA for a number of large U values up to U = 64, which we show in figure 3.By considering the doublon Green function, −k+Q (0) , in the limit of large U and extracting the zero-temperature doublon occupation from the spectral function Further details may be found in Appendix B.2.This asymptotic form is shown in figure 3, allowing us to conclude that the double occupation function is suppressed at large U according to D ∝ 1/U 2 .The results from a number of independent numerical studies, also shown in figure 3 are fully consistent with 1/U 2 scaling at large U.At half-filling, the value of D should change monotonically from 0 to 1/4, which correspond respectively to the fully localised (U = ∞) and the completely delocalised cases (U = 0).As U decreases, charge fluctuations are enhanced and the on-site Coulomb interaction is screened, making localisation effects smaller and reducing the Mott gap.The result is larger D values for smaller U, indicating an increasing mixing of the lower and upper Hubbard bands.Still, it is only at very small values, U < 2, of the on-site repulsion that our results begin to indicate the breakdown of the approximation, which benchmarks the limits of its applicability.They do not approach D = 1/4 at U = 0, which is no surprise because neither the strong-coupling slave-fermion formalism nor the single-mode spin-wave approximation to a Néel antiferromagnetic ground state as a description of the spin dynamics is appropriate in this limit.
In general, it is not true that D should increase continuously with decreasing U.In some variational Monte Carlo calculations [58,60], the behaviour of D appears to show an abrupt change at a critical interaction strength, U c , which has been interpreted as a first-order Mott transition.However, neither our result nor those obtained from cluster perturbation theory [86], quantum Monte Carlo simulations [18], or DMRG [50] contain any evidence for a Mott transition at finite U.As a consequence of the perfect nesting, the nearest-neighbour square-lattice Hubbard model is a special case where even a weak-coupling treatment gives an insulating state for all U, i.e.U c = 0.As noted in section 1, this small-U result has led to intense debate over the question of whether insulating behaviour could be driven by antiferromagnetism rather than by the Coulomb interaction, and whether there could be a transition between the two regimes at finite U.However, the result that U c = 0 in this model means that all values of U are continuously connected to the strong-coupling limit, where the answers are clear.Indeed, detailed numerical calculations have recently been used to argue [43] that the model also has no Mott-Hubbard transition at any finite U in the paramagnetic phase at low but finite temperatures, meaning that no other effective terms are generated.We conclude from our calculations of D that the holon-doublon description yields the correct functional form and semi-quantitative accuracy throughout the regime of intermediate and strong coupling (specifically, U > 2).

Derivation and Calculation
We calculate the spectral function of the original electron operators, c kσ , which in the slave-fermion framework are decomposed into convolutions of the holon operator, e, the doublon operator, d, and the Schwinger boson operator, s.The electron Green function is defined by and the spectral function the imaginary part of the retarded Green function, contains implicitly all information necessary to describe single-particle excitations.The electron density of states is obtained from the sum over all wavevectors, Figure 4 shows the single-particle spectral function, A(k, ω), for some specific highsymmetry directions in the Brillouin zone of the square lattice.These calculations were performed with U = 8 on a 48 × 48 lattice and we used a broadening parameter η = 0.08; we stress again that this system size is effectively in the thermodynamic limit.The particle-hole symmetry of the spectral function, It is clear from the spectral-intensity contour map of figure 4(a) that the spectral weight lies in two separate bands, the lower and upper Hubbard bands, which are separated by the Mott gap; this weight appears predominantly in the lower Hubbard band in regions of the Brillouin zone around (0, 0), and is transferred to the upper Hubbard band as k moves towards (π, π).
We begin our comments on the nature of these results by noting that they are in quantitative agreement with calculations performed by the variational cluster approximation [40], in which the authors included long-range antiferromagnetic correlations by adding some Weiss fields.That we obtain the self-consistent interaction effects of the spin fluctuations on the charge degrees of freedom is one of the key qualitative features of our approach and will be discussed further below.The shifting of spectral weight between Hubbard bands signals that the real part of the electron Green function changes sign between (0, 0) and (π, π), which implies that the selfenergy diverges and the Green function has a zero surface falling within this region [64].The importance of zeros and poles of the Green function has been stressed by many authors [24,[64][65][66][67]87] in the discussion of Mott physics, particularly in the context of pseudogap phenomena in the doped system.We defer a discussion of the zero surface of the electron Green function to section 5.
Results for the spectral function, A(ω), are shown in figure 4(b) for a sequence of wavevectors k along the high-symmetry directions of the Brillouin zone.The spectral function is always a multi-peak structure, and the four-peak form reported by references [28] and [40] is not very distinctive in our results.Most of the spectra consist of a peak at low energies, meaning closer to the chemical potential, accompanied by a broad, weak high-energy part.The low-energy peak can be interpreted as the coherent motion of the quasiparticles, which are formed from both holons and doublons and are renormalise by spin fluctuations.The high-energy part originates in incoherent charge excitations and forms the remainder of the Hubbard bands.
Returning to the details of figure 4(b), we observe that for k = (π, 0) (the X point), the spectral weight is symmetrical about ω = 0.As k is moved along X → Γ,  weight is transferred from the positive-to the negative-ω regime, i.e. to lower energies, while along X → M the opposite happens; this behavior is a reflection of particle-hole symmetry.Because (π, 0) lies on the noninteracting Fermi surface, when k lies inside this surface, most of the spectral weight has the properties of a particle, while outside it has the properties of a hole.An analogous spectral-weight redistribution is evident as k is moved from M to Γ, where the transfer is from positive to negative ω.At k = (π/2, π/2) ≡ S, the distribution is again symmetrical in ω and the gap takes the minimum value obtained in our approximation.We note that all of these features are very similar to the results of exact diagonalisation [20][21][22][23] and quantum Monte Carlo calculations [25][26][27][28] of the Hubbard model at half-filling on the square lattice.
Clearly the leading momentum-dependence of the spectral function is the consequence of the underlying non-interacting electronic bands.Around X ≡ (π, 0), where this dispersion is rather flat [figure 4(a)], the high density of states and strong interactions are believed to be the origin of the pseudogap and Fermi-arc phenomena found in doped Hubbard models [36,37].To gauge how much of the momentumdependence may be caused by the coupling between spin and charge degrees of freedom, we return to the question of the coherent and incoherent response.For the lower Hubbard band, there is an apparent suppression of spectral weight around ω ≈ −3.5t, visible near Γ, which separates the band into two parts.The "low-energy" part (meaning closer to ω = 0) is a quasiparticle-type band, which disperses from Γ to S.
The emergence of a coherent quasiparticle band in the dynamical response is a highly nontrivial consequence of the interactions between the charge and spin sectors in the Hubbard model.In general, hole motion in a Néel ordered state is energetically unfavourable because the movement distorts the antiferromagnetic background.However, many calculations within the SCBA [70][71][72][73][74][75] have shown that a quasiparticle, named the "spin polaron" [73,74], forms as a result of holon-magnon coupling.The bandwidth of the spin polaron is governed by the spin exchange, J, and this feature forms the low-energy part of the lower Hubbard band.For the high-energy part, hole motion is thought to originate from effective three-site hopping processes [38], which allow hole propagation on the same sublattice that does not distort the antiferromagnetic background and therefore is less affected by spin fluctuations.We postpone a more detailed discussion of the Hubbard bands to the next section.

Density of States
The density of states, ρ(ω), obtained by integrating the spectral function is shown in figure 5, where we compare the results calculated in the SCBA with those of a simple holon-doublon mean-field approximation.The dominant feature in ρ(ω) is the strong Mott gap for the values of U illustrated.In the mean-field results, shown in figure 5(a), the effects of spin fluctuations are neglected and the Mott gap is U.However, in our calculations [figure 5(b)] this gap is renormalised downwards, quite significantly at intermediate values of U.This renormalisation is related directly to the other significant feature of ρ(ω), the width of the lower and upper Hubbard bands, which clearly broadens (relative to U) as the on-site interaction decreases and spin fluctuations strengthen.As we will discuss in section 5, the bandwidth reflects the spinrelated renormalisation of the underlying quasiparticle bands and may be treated as an effective hopping parameter, t eff , which is a fraction of the electron hopping parameter, t; by contrast, the Hubbard bands given by the mean-field approximation have a width given only by the spin interaction, J = 4t 2 /U, and therefore appear much narrower in figure 5(a).
In the SCBA, only the quasiparticle band remains of width J.In general, the spin polaron is a feature of the retraceable path approximation [76] and is well-defined in a mean-field treatment with high lattice coordination.The energetics of the Mott gap and the Hubbard bands have been discussed in the DMFT, which is infinite-dimensional, and the same physics of a narrow quasiparticle band, a renormalised Mott gap and broader Hubbard bands is found [88].That our 2D holon-doublon description reproduces these features with high accuracy serves both as an indication of the degree to which it captures the key ingredients of Mott physics and as a means to verify the relevance of features observed in numerical calculations with different limitations.In the context of DMFT comparisons, we caution that the small peaks visible around ω/U ≈ 0.5 in figure 5(b) are size effects, which we can demonstrate by smoothing them away if the broadening factor η is set to a larger value.
Although a number of experiments have demonstrated the insulating behaviour of the parent compounds of the high-T c superconductors, it has proven to be very difficult to observe both the lower and upper Hubbard bands.Angle-resolved photoemission spectroscopy [89] probes the occupied states while optical spectroscopy [90,91] gives information concerning two-particle correlations.Measurements by scanning tunnelling spectroscopy (STS) probe the local single-particle density of states for both occupied and unoccupied bands.Only recently was the full electronic spectrum for the undoped Mott insulator Ca 2 CuO 2 Cl 2 measured by STS [92], showing the Mott-Hubbard gap, and in particular the upper Hubbard band, for the first time.A comparison with the robust features of our results, namely the spin polaron, the renormalised Mott gap (section 4.3), and the width of the Hubbard bands scaling with t, implies that the Hubbard model does in fact capture the essential properties of the undoped Mott insulator measured by STS.

Mott Gap in the Large-U Limit
We reiterate that the SCBA includes the interaction effects of the holons and doublons with the magnons, causing significant renormalisation effects at intermediate values of U.Only in the limit of large U do these effects vanish, causing the Mott gap to tend towards U and the bandwidths to vanish.We define the Mott gap, ∆, from the separation of the peaks in the density of states shown in figure 5(b) and present the results in figure 6.We comment that, although this definition is completely valid in the large-U regime, it may be called into question for values of U/t at the far right-hand side of figure 6: for U/t = 4 it is clear in figure 5(b) that there is a significant density of states "inside the Mott gap" and a detailed account of its effects could lower the estimate of the effective ∆ [32].Spectral information is significantly more difficult to obtain from numerical calculations than are ground-state properties, and in figure 6 we show only Mott-gap estimates obtained from exact-diagonalisation studies [21].These lie very close to our SCBA results at large U and then fall with a similar form, but to values slightly smaller than SCBA, as U/t is decreased through the intermediate-coupling regime.
At the left-hand side of figure 6, similar to our extraction of the limiting behaviour of D in section 3, we may also investigate analytically the approach of the Mott gap to U as U becomes large.As shown in Appendix B.3, by neglecting off-diagonal terms (F 12 and F 21 ) in the Green function of equation ( 11) and using the bare holon and doublon Green functions to obtain the spectral function, and hence the self-energies, one may construct the full Green function at the level only of the first step in the Born approximation.The Mott gap is simply the minimum of the holon and doublon gaps, which are identical and, on retaining only terms of order unity, are given by where the gap minimum occurs at k = (π/2, π/2) = S.We show this function in figure 6 in the form ∆/U = 1 − 3.77/U.The agreement is satisfactory, confirming that the limiting form of the approach is indeed 1 − ∆/U ∝ 1/U, although the constant of proportionality could be computed more accurately by retaining more orders in the Born approximation.We comment that this type of limiting behavior has been suggested in some numerical studies [88].It is also worth noting [32] that the Mott gap obtained in 1D from the Bethe Ansatz has the large-U limit ∆/U = 1 − 4/U to first order in 1/U.In the high-U (atomic) limit, the spectral function is entirely incoherent.The mean-field approximation gives the "unperturbed" charge Green function, meaning in the absence of spin fluctuations, and the two Hubbard bands in figure 5(a) are also incoherent.In a self-consistent treatment, the spin fluctuations cause a spectral-weight transfer from the incoherent background into a coherent quasiparticle band, allowing both the (incoherent) lower and upper Hubbard bands and the (coherent) low-energy quasiparticle band to be reconstructed, as shown in figure 4(a).These results emphasise again that a consistent inclusion of the spin degrees of freedom is intrinsic to a full understanding of the Mott state.

Spectral-Weight Integral
We comment here that, from the properties of the Green function, the density of states of electrons with both spin orientations should obey the sum rule In our results, because one of the two bosonic operators on each sublattice undergoes Bose-Einstein condensation [equation ( 7)], the sum rule is violated and instead takes the form where D is the average double-occupancy parameter calculated in section 3.In figure 7 we show both the spectral-weight integral and its value corrected by the average double occupancy, for a range of values of U.Because our calculation is performed using a finite energy interval, the corrected spectral-weight integral is less than 1, but it is clear that the deviations are very small for all intermediate values of U. In the atomic limit, D = 0, the sum rule is satisfied by the spectral weight alone, while for finite U the spectralweight integral violates the sum rule by an amount equivalent to D and approaches 1 as U increases.For an intermediate interaction strength such as U = 8, we observe that the spectral-weight integral is 0.96, indicating that our strong-coupling approximation remains robust, and capable of capturing the primary physics of the Hubbard model, in this regime.Only at weak coupling, where the spin-wave treatment of the Heisenberg model is not suitable, higher-order perturbations are important and charge fluctuations become strong, does the SCBA fail to reproduce the complex charge dynamics.

Luttinger Volume
Landau's Fermi-liquid theory is one of the fundamental theories in quantum many-body physics.Its most important single qualitative result is that the Fermi surface still exists, even in the presence of interactions, in all dimensions higher than 1.For the normal metals that can be described as Fermi liquids, the shape of the Fermi surface is essential in determining the thermodynamic and transport properties of the system, because only electrons close to the Fermi surface can be excited at low energies.For a given electron density, regardless of the shape of the Fermi surface, the volume it encloses is invariant even when interactions between particles are taken into consideration.This result, the Luttinger theorem [93][94][95], relates the total electron density to the volume in momentum space where the Green function is positive, with d the spatial dimension of the system.Although the original proof by Luttinger was based on perturbation theory, and thus the theorem may in principle be violated by nonperturbative effects, Oshikawa [95] has provided a nonperturbative proof based on topological arguments.In essence, the insensitivity of the result to any interactions should be considered as a quantisation phenomenon [95], and therefore the Luttinger theorem may be the first example of topological quantisation discovered for the quantum many-body problem.
In general, the poles of the electron Green function, G(k, ω), determine the band dispersion and define the Fermi surface at the chemical potential, where ω = 0. From the asymptotic behaviour G(k, ω → ∞) ∝ 1/ω, it is obvious that the limiting Green function is positive for positive ω and conversely.For a Fermi liquid, as the frequency approaches the Fermi surface (ω = 0), by definition G(k, ω = 0 + ) → +∞ while G(k, ω = 0 − ) → −∞, implying that the Green function must change sign at the Fermi surface through an infinity in G(k, 0) if there are no singularities in any other regions of k and ω.The Luttinger theorem then concerns precisely the volume enclosed by the surfaces where G(k, 0) changes sign.
However, this change of sign is not always connected with an infinity.It may also occur through zeros of the Green function [24,[64][65][66][67]87], and in fact these are the relevant quantities for insulating systems, where there is no Fermi surface.By contrast, all systems have a Luttinger surface, and this may be defined from the zeros.Examples of this type of physics include the BCS superconductor, where G(k, 0) changes sign through the chemical potential in the absence of a Fermi surface, and also the Mott insulator.
For systems with a full gap spanning some energy interval (figures 4 and 5), the self-energy inside the gap must be infinite as otherwise spectral features would appear.For the Green function, the divergence of the self-energy implies (11)  which in terms of Re[G −1 (k, ω)] corresponds to its points of divergence, whereas the band dispersion is given by its zeros.We note that the self-energy can have only one pole, ξ k , in an energy regime where a gap is present; a situation with two successive poles, ] in the gap region, contradicting the definition.
In figure 8 we show the real part of the inverse Green function in the form of colour contours.The finite broadening factor, η, in our calculation converts the divergence into a sharp minimum followed by an equally sharp maximum as a function of |k|.It is clear that Re[G −1 (k, ω = 0)] changes sign from negative to positive through a diamondshaped boundary, which therefore marks the zeros of the Green function at ω = 0, i.e. the Luttinger surface.This surface is specified precisely by the condition cos k x + cos k y = 0, and therefore we have found that the Luttinger surface is exactly the Fermi surface of non-interacting electrons.This result is a consequence of particle-hole symmetry at half-filling.Although the system is completely gapped in the Mott-insulating phase, the Luttinger theorem remains applicable and for the half-filled Hubbard model can be interpreted as the integral of the region within the Luttinger surface [69].
To gain more insight into the behaviour of the self-energy, in figure 9 we show Re[G −1 (k, ω)] calculated with U = 8 for the primary high-symmetry directions in momentum space.The divergence of the inverse Green function can be regarded as a dispersion relation defined by the self-energy, and it forms a rather clear quasiparticle band.The black, dashed curve in figure 9 shows the function 2t eff (cos k x + cos k y )  with t eff = 0.8t, which obviously provides a reasonable, if not perfect, fit to the zeros of the calculated Green function.This form, which resembles an inverted freeelectron dispersion with a renormalised hopping parameter, is quite similar to the results obtained by exact diagonalisation [24], except in that the effective hopping is smaller.
To interpret this result, we note first that the inverted nature of the quasiparticle band signifies its holonic origin.From equation (A.11), we expect that the doublon band will have the same form as the holon band but with a momentum shift of Q.Thus the poles of the electron self-energy are given directly by the dispersion relation of the quasiparticles, which is in turn connected with the non-interacting dispersion, 2t(cos k x + cos k y ), by a renormalisation factor.Second, from equation (A.16) one observes that this renormalisation of the effective quasiparticle dispersion, 4tb 2 0 γ k = 2tb 2 0 (cos k x + cos k y ), is related directly to the holon-doublon pairing interaction term, which in turn is determined by the magnetic ordering (condensation) parameter, b 0 .In the limit of large U, the Mott gap saturates at ∆ = U (section 4.3) and, as discussed in Appendix B.1, magnetic order becomes robust, with an ordered moment of m s = 0.321 and, from equation ( 8), b 2 0 = 0.821 in our calculations.The effective quasiparticle hopping parameter, t eff , is shown in figure 10, which compares the results extracted from the Green function in a simple mean-field approximation, with only holon-doublon interactions but without spin-fluctuation effects, to the results of the SCBA.Also shown for comparison is the result of a direct determination of b 2 0 , which allows us to deduce that the effective hopping parameter saturates at t eff = b 2 0 t = 0.821t as U → ∞ in the mean-field approximation.However, the effect of the additional spin fluctuations contained in the SCBA is to allow the effective hopping to be stronger than the value constrained by the condensation parameter.
Away from strong coupling, the departure of the quasiparticle mass from the high-U limit may be taken as a further indication of the effects of charge and spin fluctuations at intermediate U, which includes the value U = 8 shown in figure 9.The reduction of t eff computed in the SCBA with decreasing U is clear in figure 10, where again the comparison with the mean-field form and b 2 0 allows us to benchmark the effects of spin fluctuations and of self-consistency at intermediate U within our calculation.The Mott gap is effectively the mass term, U, in equation (A.16) and its downward renormalisation with decreasing U is also the consequence of holon-doublon interactions, renormalised at intermediate U by spin fluctuations, as discussed in section 4 and shown in figures 4 and 5.However, this information is not visible in the inverse Green function, which therefore can be taken as a measure only of the spin-fluctuation renormalisation effect within the SCBA, independently of the Mott-gap scale.
We conclude that the holon-doublon description contains valuable insight into the static and dynamic properties of the charge degrees of freedom.The characteristic feature of the Mott insulator is not the Fermi surface but the Luttinger surface.The Luttinger theorem for the Fermi liquid holds for the half-filled Hubbard model, despite its insulating nature, with the sign-change of the Green function at ω = 0 caused by its zeros instead of its poles.This Luttinger surface appears inside the Mott gap, which separates the lower and upper Hubbard bands, and it defines a quasiparticle band whose renormalisation characterises the effects of spin interactions on the charge sector.

Optical Conductivity
The final quantity we calculate is the optical conductivity.Optical conductivity measurements [90,91], which reveal the carrier number, the size of the energy gap, the dynamics of quasiparticle excitations and their scattering processes, have played an important role in the study of many classes of correlated electronic materials.For the insulating parent compounds of the high-T c superconductors, optical conductivity measurements performed on La 2 CuO 4 at a temperature of 300 K show an excitonic absorption peak at 2 eV [96].This behaviour was ascribed to the strongly correlated but charge-transfer-dominated character of the undoped CuO 2 plane [97].When holes are doped into the system, optical conductivity measurements show a number of features that deviate strongly from conventional band theory.The most striking phenomenon is the reconstruction of the electronic spectral weight at low doping, which involves the transfer of weight from the charge-transfer excitation regime to a midinfrared band, centred at approximately 0.5 eV.This mid-infrared band is consistent with the appearance of in-gap states [86,92] but to date there is no theory for its microscopic origin.As the doping is increased, a Drude-type response develops at farinfrared frequencies and decays much more slowly than band theory would predict.This characteristic spectral-weight transfer is evidently intrinsic to the Mott insulator and cannot be described by a theory ignoring electronic correlations.It has been argued [63,98] that this property of the Mott phase can be explained within the Hubbard model, but there is as yet no consensus on the complete dynamics of the weight-transfer phenomenon.
The optical conductivity in the x direction can be calculated from the imaginary part of the current-current correlation function, as discussed in detail in Appendix C. The results we obtain within the slave-particle framework and the SCBA are shown in figure 11 for a number of U values spanning the intermediate-coupling regime.The optical conductivity has a clear charge-excitation gap, as expected for an insulating system (but in contrast to the Drude peak predicted by conventional band theory).The first peak is expected at the "optical gap," ∆ Optical , which is determined by a two-particle process; as U increases, its location clearly moves to higher energy, even when normalised to U, and its peak height drops.This is consistent with experimental results for the optical reflectivity spectra [96], and with the expectation of a larger quasiparticle gap with increasing U contained in figure 5 Exact diagonalisation results [23,99,100] and quantum Monte Carlo calculations [26] also confirm the insulating character of the undoped Hubbard model.The optical conductivity given by exact diagonalisation [100] shows a sharp peak at the Mott-gap edge, which is attributed to spin-polaron formation in the photoexcited state, as also revealed by DMFT [101].Although these results suggest a separation of the spin-polaron band from the high-energy feature, we caution that the system sizes used are very small; in our calculations, the spectral function has a continuous and multipeak structure for most wavevectors k, with no clear band separation and thus no discrete peaks in the optical conductivity.We do not attempt a quantitative comparison of our results with other numerical methods because the effect of vertex corrections [102,104] cannot be neglected in the calculation of optical conductivity.
As a consequence of fundamental conservation laws, the optical conductivity obeys the sum rule [90,91,104 where T x is the average kinetic energy in the x direction.In our calculations for the half-filled Hubbard model, the results do not satisfy this sum rule exactly, although the deviation is within 10% for all intermediate U values.This discrepancy is presumably caused by the spectral-weight loss discussed in section 4, as well as by our neglect of vertex corrections.The inclusion of the latter, and of higher-order physical processes, may be expected to produce superior results, and to be increasingly important for calculations at lower U or at any finite temperatures.

Summary
We have used a holon-doublon slave-fermion representation in the self-consistent Born approximation to investigate the charge dynamics of the single-band Hubbard model at half-filling on the square lattice.We show that holon-doublon binding is intrinsic to the Mott-insulating state.As expected in a nearest-neighbour hopping model, the system is always insulating and the monotonically decreasing doublon density shows no critical value, U c , where a metal-insulator transition occurs.The electronic density of states we compute has a coherent quasiparticle band accompanied by incoherent lower and upper Hubbard bands.The Mott gap is the excitation gap of the holons and doublons, and so is determined by their binding energy.We calculate the Luttinger surface of the Green function at the chemical potential and demonstrate that, for the Mott insulator, the Luttinger theorem can be interpreted as the integral of the region enclosed by this Luttinger surface rather than a Fermi surface.We show further that the poles of the self-energy define an effective quasiparticle dispersion similar to that of free electrons, but with a renormalised hopping parameter.The optical conductivity displays clearly both the Mott gap and the incoherent Hubbard bands.
In this slave-fermion description, the Mott-gap state is accompanied by the formation of holon-doublon bound states.Holon-doublon binding provides the energy scale of the Mott gap in the half-filled Hubbard model.In greater detail, the pairing interaction of holons and doublons is a k-space phenomenon, occurring across the Fermi surface.We comment that both holons and doublons are components of the quasiparticles making up the lower Hubbard band, as they are of the upper, and thus the paired state should not be considered as an exciton; indeed, holon-doublon pairs would remain present at finite doping, when the particle-hole symmetry of the halffilled system is lost.The presence of these bound states, together with the interactions between the charge carriers and the magnons, reflects the fact that the high-and lowenergy degrees of freedom in the Hubbard model are intrinsically mixed, with the spin fluctuations playing an essential role in the construction of coherent "spin-polaron" features in the lower and upper Hubbard bands.
Our technique is a strong-coupling approach, valid at large values of U, and can be used to obtain the analytical limiting forms of all physical quantities in this regime.At intermediate U values, it provides accurate results for the evolution of the charge dynamics as a consequence of the renormalising effects of spin fluctuations.Only on the approach to the weak-coupling regime do we find systematic discrepancies in the framework, which we may characterise by violations of the spectral sum rule.In this regime, the single-mode approximation to the spin dynamics of the Néel antiferromagnetic ground state begins to break down, the self-consistent Born approximation requires the inclusion of vertex corrections and even the decoupling of electron operators into separate spin and charge parts may no longer be justified.As noted in section 1, a degree of confusion exists in the literature over whether the insulating properties of the weakly coupled system can be ascribed to magnetic rather than to electrostatic interactions, and we conclude from the absence of a transition that the two pictures are different sides of the same coin, driven by the same fundamental processes.However, the slave-fermion framework adopted here is clearly not the appropriate method for investigating the crossover from the robustly Mott-insulating regime to the noninteracting limit.
Despite this deficiency far from its regime of validity, we have demonstrated that the holon-doublon representation provides both a full qualitative understanding of the underlying physics of the Mott insulator and a semi-quantitative account of its physical properties across the full range of intermediate and strong interactions.Our results for the Mott state have a transparent origin, not obscured by the complexities of a many-body numerical calculation, and indeed can be used to verify the features and parameter-dependences found in numerical studies.Although one may argue that the half-filled system at zero temperature is already fully characterised, the holon-doublon description we have introduced here provides a valuable foundation for understanding both the charge dynamics of the half-filled Hubbard model at finite temperatures and the evolution of the Mott gap as the system is doped away from half-filling.In the strong-coupling limit, to lowest order one may consider only the charge part of the Hamiltonian, H C in equation (A.12), and neglect the terms H CS describing the interaction of the charge degrees of freedom with the spin fluctuations.In this approximation, the doublon Green function may be expressed as one obtains the spectral function and thus the zero-temperature doublon occupation function By Taylor expansion of this expression, and the net doublon occupation is showing asymptotic 1/U 2 scaling in the large-U limit.

Appendix B.3. Mott Gap
In the limit of large U, one may assume that the strongest contributions to the renormalisation of the Mott gap are obtained at first order in the SCBA, with higherorder contributions being heavily suppressed.In this case it is safe to assume that the off-diagonal Green-function components, F 12 and F 21 in equation (A.17), may be neglected in computing the self-energy corrections determining the Mott gap.From the doublon spectral function of equation (B.8) and its holon counterpart the self-energy components [equation (A.18)] calculated in the SCBA are and where a k = q g 2 1 (k, q).(B.15) The SCBA Green function [equations (A.16) and (A.17)] is then These functions both have four poles, which are the same in both cases, where we have defined b k = (4tb 2 0 γ k ) 2 .Thus the holon and doublon gaps are both given by in the limit of large U.
The electron Green function is the convolution of the holon and doublon Green functions, which include the magnon renormalisation effects, and hence the Mott gap is given by the holon and doublon gaps.From figure 4(a), it is clear that the minimum value of this gap corresponds to the high density of states giving the clear peaks in figure 5 in the strong-coupling limit.We conclude that the approach of the Mott gap, ∆/U = 1 − 3.77/U, to unity due to the suppression of spin fluctuations in this regime is proportional to 1/U.Because this derivation considers only the first order in the self-consistent renormalisation scheme, one may expect the coefficient to be larger than 3.77 if higher-order processes are included.

Appendix C. Optical conductivity
The conductivity tensor σ µν (q, ω) is defined as the linear response of the charge current density in a solid to the total electric field, ∧ J µ (q, ω) = ν σ µν (q, ω)E ν (q, ω), (C.1) where µ and ν denote Cartesian coordinates and . . .denotes the equilibrium average; we use units where e = 1, = 1 and c = 1.The long-wavelength (q = 0) limit of the real part, Re[σ µν (ω)], is referred to as the optical conductivity, because the wavelength of electromagnetic waves far exceeds the characteristic length scales of condensed matter systems, so here we derive this limit.For a lattice system, the hopping parameter in the presence of an electromagnetic field is expressed by the Peierls substitution [104] as where j i dl • A(r, t) is the integral of the vector potential along the hopping path.For a spatially uniform electric field, the vector potential can be taken as independent of position, A(r, t) → A(t).By restricting to the linear-response regime, expanding the phase factors to second order in A(t) and taking the Fourier transform, one obtains the current operator in the form ∧ J µ (k, t) = ∧ J (1)  µ (k, t) + ∧ J (2)  µ (k, t), (C.3)where ∧ J (1)  µ (q, t) = − µ (k, t) is known as the paramagnetic current and ∧ J µ (k, t) as the diamagnetic one.The latter is already linear with respect to the applied field but for the former the Kubo formula is required.To linear order in the applied field, the interaction term in the Hamiltonian is given by µ (−q)E µ q e −iωt . (C.6) The imaginary-time current-current correlation function is defined as Π αβ (q, τ ) = − 1 N T τ ∧ J (1)  α (q, τ ) ∧ J both of which are obtained from the charge Green function of equation (A.17).Finally, the conductivity in the x direction is given by whence the real part of the optical conductivity in the limit q = 0 is

Figure 1 .
Figure 1.Feynman diagrams for the SCBA; fermion and magnon propagators are represented respectively by straight and wavy lines.

Figure 3 .
Figure 3. Double occupancy parameter, D, shown as a function of 1/U 2 on logarithmic axes.Red points show our large-U results, calculated for a system size of 48 × 48 and with broadening parameter η = 0.08.Results from CPT, QMC and DMRG, shown in the same colour and symbol scheme as in figure 2, confirm the trend towards the limiting large-U functional form, D = 2.70/U 2 (dashed purple line).

Figure 4 .
Figure 4. Electron spectral function, A(k, ω), of the Hubbard model with U = 8, obtained for lattice dimensions 48 × 48 and with broadening parameter η = 0.08.(a) Spectral intensity along the symmetry directions Γ → M → X → Γ in the first quadrant of the first Brillouin zone.(b) Spectral function A(ω) for wavevectors k along the high-symmetry directions X → Γ → M → X → S.

Figure 5 .
Figure 5. Electronic density of states, ρ(ω), shown for different values of U ; data are normalised to U on the frequency axis.A holon-doublon mean-field approximation (a) is compared with the results of the SCBA (b).

Figure 6 .
Figure 6.Mott gap, ∆, normalised to U and shown as a function of 1/U to illustrate its asymptotic large-U limit.The purple stars indicate results for U = 4, 8 and 40 extracted from the exact diagonalisation (ED) calculations of reference [21].

Figure 8 .
Figure 8. Real part of the inverse Green function, Re[G −1 (k, ω = 0)], for the half-filled Hubbard model on the square lattice, calculated for U = 8.Re[G −1 (k, ω = 0)] changes sign from a negative to a positive maximum through a diamond-shaped boundary.The sign-change gives the Luttinger surface at ω = 0, and coincides with the locus of solutions of cos k x + cos k y = 0, marked by the black, dashed line.

Figure 9 .
Figure 9. Real part of the inverse Green function, Re[G −1 (k, ω)], calculated using U = 8 for values of k chosen along the high-symmetry directions.The black, dashed curve shows the function 1.6t(cos k x + cos k y ) (see text).

Figure 10 .
Figure 10.Effective hopping parameter, t eff , defined by the poles of the electron self-energy, 2t eff (cos k x + cos k y ), calculated in a mean-field approximation with no spin-fluctuation effects (blue) and in the SCBA (red); for comparison we show also the condensation parameter, b 2 0 (green).

Figure 11 .
Figure 11.Real part of the optical conductivity as a function of frequency, shown for a number of U values.The results were calculated on a 48 × 48 lattice at zero temperature.The red, dashed line connects the initial peaks on each curve, which determine the optical gap, ∆ Optical .The inset compares ∆ Optical with the Mott gap, ∆, determined from the peak separation in figure 5.
. The Mott gap is the one-particle charge gap and is determined from the separation of the leading peaks in the lower and upper Hubbard bands, as shown in the spectral function [figure 4(a)] and the density of states [figure 5(b)].The inset of figure 11 compares the Mott gap with the optical gap obtained in the main panel, and it is clear that the two are almost identical for the Mott insulator.