A Diagrammatic Approach

# Electron spectral functions in a quantum dimer model for topological metals

## Abstract

We study single electron spectral functions in a quantum dimer model introduced by Punk, Allais and Sachdev in Ref. [M. Punk, A. Allais, and S. Sachdev, Proceedings of the National Academy of Sciences 112, 9552 (2015)]. The Hilbert space of this model is spanned by hard-core coverings of the square lattice with two types of dimers: ordinary bosonic spin-singlets, as well as fermionic dimers carrying charge +e and spin 1/2, which can be viewed as bound-states of spinons and holons in a doped resonating valence bond (RVB) liquid. This model realizes a metallic phase with topological order and captures several properties of the pseudogap phase in hole-doped cuprates, such as a reconstructed Fermi surface with small hole-pockets and a highly anisotropic quasiparticle residue in the absence of any broken symmetries. Using a combination of exact diagonalization and analytical methods we compute electron spectral functions and show that this model indeed exhibits a sizeable antinodal pseudogap, with a momentum dependence deviating from a simple d-wave form, in accordance with experiments on underdoped cuprates.

## I Introduction

Since the discovery of high-temperature superconductivity in cuprates Bednorz and Müller (1986), many efforts have been made to understand the underlying pairing mechanism. A potential solution to this problem might come from a better understanding of the so-called pseudogap state at low doping, out of which the superconducting phase develops upon further doping. By now we know from a wide range of experiments that the pseudogap state is largely dominated by antiferromagnetic fluctuations Niedermayer et al. (1998); Kastner et al. (1998); Tallon and Loram (2001); Chan et al. (2016) and has a tendency to charge density-wave ordering (CDW) Wu et al. (2011); Ghiringhelli et al. (2012); Chang et al. (2012); Blanco-Canosa et al. (2013); Comin et al. (2015).

During the last decades experiments have explored the complex physics of the pseudogap phase in much detail. Early studies of the Knight shift in NMR measurements have shown that the magnetic susceptibility decreases below a characteristic temperature scale, indicative of a spin gap Alloul et al. (1989); Curro et al. (1997). Soon measurements of the tunneling density of states Renner et al. (1998), c-axis optical conductivity Homes et al. (1993); Uchida (1997) and specific heat Loram et al. (1993, 2001) revealed that the opening of the pseudogap appears in both, charge and spin degrees of freedom.

Interestingly, some transport experiments show rather ordinary metallic properties in the pseudogap phase. For example, in-plane optical conductivity and magnetoresistance measurements show Fermi liquid-like behavior Mirzaei et al. (2013); Chan et al. (2014). However, Hall coefficient as well as the Drude weight indicate that the charge carrier density is small and proportional to the density of doped holes away from half filling Orenstein et al. (1990); Uchida et al. (1991); Cooper et al. (1993); Padilla et al. (2005); Badoux et al. (2016), rather than as expected from Luttinger’s theorem and observed in the Fermi liquid regime at large doping. On the other hand, angle resolved photoemission (ARPES) experiments show a distinct Fermi-arc spectrum and a gap opening in the vicinity of the anti-nodes in momentum space Damascelli et al. (2003); Shen et al. (2005); Hashimoto et al. (2014); Shen et al. (1993); Ding et al. (1996); Lee et al. (2007); Vishik et al. (2012, 2010). Taken together, these observations are hard to reconcile with a Fermi liquid picture, where the Fermi surface is reconstructed by some thermally fluctuating order parameter. In this case one would expect some signatures of the order parameter correlation length to be observable, e.g. in transport measurements.

A different set of theoretical ideas, based on Anderson’s resonating valence bond (RVB) picture Anderson (1973); Fazekas and Anderson (1974); Anderson et al. (1987), tries to approach the pseudogap from the Mott insulating state at low doping Lee et al. (2006). Upon doping the RVB state with holes, electrons fractionalize into neutral spin-1/2 spinon excitations, as well as spinless holons, carrying charge +e. While electron fractionalization in quasi two-dimensional systems has been a topic of considerable interest on its own, it became clear early that simple incarnations of the RVB state do not capture the sharp Fermi-arc features observed in photoemission experiments, which requires that a certain part of the low-energy spectrum behaves like electron or hole-like quasiparticles which carry both spin and charge.

One possible solution to this puzzle might be provided by the idea of fractionalized Fermi liquids (FL*), introduced originally in the context of heavy Fermion systems Senthil et al. (2003). In the context of cuprates, a fractionalized Fermi liquid can be viewed as a doped RVB liquid where spinons and holons form bound states. These hole-like bound states then form a Fermi liquid, with the size of the Fermi surface proportional to the density of doped holes. Note that the FL* inherits topological order from the RVB background, which accounts for the violation of Luttinger’s theorem Senthil et al. (2004); Kaul et al. (2008); Qi and Sachdev (2010); Mei et al. (2012); Punk and Sachdev (2012); Sachdev and Chowdhury (2016). The quantum dimer model introduced in Ref. Punk et al. (2015) provides an explicit and intuitive lattice realization of a fractionalized Fermi liquid and allows to directly compute electronic properties. Subsequent numerical work indeed showed clear signatures of a small Fermi surface, indicative of a FL* ground state Lee et al. (2016). In this model bound states between spinons and holons are represented by fermionic dimers on nearest-neighbor sites, which resonate with the background of bosonic spin-singlet dimers. At vanishing doping the model reduces to the well-known Rokhsar-Kivelson model Rokhsar and Kivelson (1988). An exact solution has been found along a special line in parameter space for an arbitrary density of fermionic dimers and the ground state can be shown to be an FL* in the vicinity of this line Feldmeier et al. (2017). We note that doped quantum dimer models have been studied previously, with doped holes as monomers occupying a single lattice site, carrying no spin Rokhsar and Kivelson (1988); Poilblanc (2008); Lamas et al. (2012).

In this work our aim is to compute electronic spectral functions for the dimer model introduced in Punk et al. (2015). Using a combination of exact diagonalization and analytical approaches we show that this model exhibits a sizable pseudogap in the antinodal region of the Brillouin zone close to momentum and symmetry related momenta, in accordance with experimental observations in the pseudogap regime of underdoped cuprates.

This paper is organized as follows: In Sec. II we give a short overview of the quantum dimer model introduced in Punk et al. (2015), define the single hole spectral function and show how it is related to dimer correlation functions. In Sec. III we discuss our numerical results for ground state properties as well as the spectral functions at zero and at finite temperature. Section IV contains an analysis of the spectral functions in terms of a two-mode approximation for the low energy spectrum of the dimer model. Finally, in Sec. V we present a diagrammatic approach to compute the electron dispersion as well as the coherent quasiparticle residuum. We conclude with a discussion in Sec. VI.

## Ii Dimer model

The Hamiltonian of the quantum dimer model introduced in Ref. [Punk et al., 2015] acts on a Hilbert space spanned by close-packed hard-core configurations of two types of dimers living on the links of a two-dimensional square lattice (see Fig. 1): the usual bosonic spin-singlet dimers, represented by bosonic operators , as well as fermionic spin-1/2 dimers carrying charge +e, represented by fermionic operators . Here denotes the lattice site, is a spin S = index and distinguishes x- and y-links on the square lattice. A fermionic dimer represents a single electron, delocalized in the bonding-orbital between two neighboring lattice sites. Alternatively one can view it as a bound state of a spinon and a holon in a doped RVB liquid.

Within the restricted Hilbert-space of the dimer model, the annihilation operator of an electron with spin on lattice site can be uniquely expressed in terms of the bosonic and fermionic dimer annihilation and creation operators as Punk et al. (2015)

 ci,α=ϵαβ2∑η(F†i−^η,η,βDi−^η,η+F†i,η,βDi,η), (1)

where is the unit antisymmetric tensor and a sum over repeated spin indices is implied.

We follow Ref. Punk et al. (2015) and consider a Hamiltonian which acts on the states by resonating dimers of both types along short, flippable loops. The corresponding quantum dimer Hamiltonian reads

 H=HRK+H1, (2)

where

 HRK= −J∑iD†i,xD†i+^y,xDi,yDi+^x,y+1 term (3) +V∑iD†i,xD†i+^y,xDi,xDi+^y,x+1 term,

is the standard Rokhsar-Kivelson Hamiltonian Rokhsar and Kivelson (1988) and

 H1= −t1∑iD†i,xF†i+^y,x,αFi,x,αDi+^y,x+3 terms (4) −t2∑iD†i+^x,yF†i,y,αFi,x,αDi+^y,x+7 terms −t3∑iD†i+^x+^y,xF†i,y,αFi+^x+^y,x,αDi,y+7 terms −t3∑iD†i+2^y,xF†i,y,αFi+2^y,x,αDi,y+7 terms.

describes dimer resonances between a bosonic and a fermionic dimer. Additional terms are related to the ones shown explicitly by lattice symmetries. We note that further terms describing resonances between two fermionic dimers can be included as well, but are not expected to play an important role at low doping, where the density of fermionic dimers is small.

The overlap between two possible dimer configurations can be caluclated using transition graphs and decreases strongly with system size Sutherland (1988); Rokhsar and Kivelson (1988). We therefore demand that two different dimer configurations are orthogonal by construction .

In the following we use periodic boundary conditions (see Fig. 1). The Hilbert space then splits into different topological sectors labeled by a set of two integer winding numbers defined by

 Wx=∑ix(−1)ix(∑αF†i,y,αFi,y,α+D†i,yDi,y), (5) Wy=∑iy(−1)iy(∑αF†i,x,αFi,x,α+D†i,xDi,x). (6)

where the sums run over a line of lattice sites in one direction, counting the staggered number of dimers that cross the reference lines in Fig. 1. Note that any local Hamiltonian like (2) has non-zero matrix elements only between states in the same topological sector. Matrix elements between states in different topological sectors vanish.

We further make use of symmetries of the quantum dimer Hamiltonian that allow to reduce the computational cost for exact diagonalization. Particle number conservation as well as spin-rotation symmetry is implicit in the dimer representation. Finally, the Hamiltonian (2) is invariant under translations as well as square lattice point group symmetries. We make use of these symmetries in our numerical implementation to reach a maximum system size of lattice sites for which we calculate the full spectrum for one fermionic dimer embedded in a background of bosonic dimers using exact diagonalization, and a maximum size of sites for the computation of ground-state wavefunctions using a Lanczos algorithm.

In this work our main quantity of interest is the single electron spectral function, which is directly measurable in ARPES experiments Damascelli et al. (2003); Shen et al. (2005); Hashimoto et al. (2014); Shen et al. (1993); Ding et al. (1996); Lee et al. (2007); Vishik et al. (2012, 2010). The mapping in Eq. (1) allows us to directly compute the hole-part of the electron spectral function

 A−,k(ω) =2πZ∑m,ne−βEm|⟨n|ck,α|m⟩|2δ(ω−Em+En+μ) (7)

where denotes a complete set of eigenstates of the dimer model with energy and

 ck,α Extra open brace or missing close brace (8)

is the electron annihilation operator in the dimer representation. For our numerical computation we fix the particle number and take to be eigenstates of the undoped Rokhsar-Kivelson model, whereas are eigenstates of the dimer model with one fermionic dimer. The full single electron spectral function follows from the hole-part via detailed balance and is usually normalized as . It is important to note here that the electron annihilation operator defined in Eq. (8) does not obey canonical anticommutation relations, because it is a composite operator. For this reason the positive definite electron spectral function computed here does not satisfy the above normalization condition. The reason for this is that the electron annihilation operator has a non-local representation in the dimer Hilbert space and is dressed by the form factor . The normalization of the spectral function at zero temperature indeed depends explicitly on momentum and is given by , which is thus an upper bound to the quasiparticle residuum shown in Fig. 2. The spectrum of the dimer model in our finite size numerics is discrete and the spectral functions are thus composed of a series of delta function peaks with different weight. For better visibility we broaden the delta functions in our figures using a Lorentzian with a width . In order to extract quantitative results for the pseudogap from our numerical spectral function data, we define the gap function as the distance between lowest energy state at fixed momentum and the Fermi energy. It is important to realize that this is a lower bound for the energy gap in the spectral function, because our numerical data shows that the lowest energy states in the spectrum carry a vanishingly small weight in the zero temperature spectral function as one approaches the antinode (see also Fig. 2 middle), shifting the apparent gap to larger energies. This is no longer true at finite temperatures, however.

## Iii Numerical Results

### iii.1 Fermi pockets and quasiparticle residue

The ground state energy of a single fermionic dimer at fixed total momentum in a background of bosonic dimers has already been computed in Ref. Punk et al. (2015), together with the corresponding quasiparticle residue . Here denotes the ground state of (2) with one fermionic dimer at fixed total momentum , and is the ground state of the undoped RK-model. In Fig. 2 we show similar data, but with strongly increased momentum resolution, obtained using a Lanczos algorithm for a lattice with twisted, rather than periodic boundary conditions, which allows us to compute and for any momentum in the Brillouin zone.

Our results for parameters , , , and , which follow from electron hopping parameters in an effective model appropriate for cuprates (see Ref. Punk et al. (2015) for details), show a pronounced dispersion minimum in the vicinity, but not directly at momentum . We note that the position of the dispersion minimum depends on microscopic details and the precise value of the amplitudes. At a finite density of fermionic dimers this would give rise to a Fermi surface comprised of small hole-pockets with an area proportional to the density of fermionic dimers, which equals the density of doped holes away from half filling, realizing a fractionalized Fermi liquid. The corresponding quasiparticle residue drops sharply for momenta larger than and the associated photoemission response, which is proportional to the hole-part of the electron spectral function, indeed shows Fermi-arc like features due to the highly anisotropic quasiparticle residue. Note that the ground state at a finite density of fermonic dimers potentially breaks symmetries in certain parameter regimes and is not necessarily a fractionalized Fermi liquid with a Fermi surface encompassing an area proportional to the density of fermionic dimers. While the precise nature of the ground state is not yet known in the full parameter space, numerical studies showed the presence of Friedel oscillations associated with a small Fermi surface in accordance with an FL* ground state Lee et al. (2016). Moreover, a recent exact solution of the dimer model also shows that the ground state is an FL* in an interesting parameter regime Feldmeier et al. (2017).

### iii.2 Spectral functions and pseudogap

In Fig. 3 we display results for the hole spectral function Eq. (7) at zero temperature along a series of momenta between the nodal point on the Fermi pocket and the antinode at for two different system sizes, and lattice sites, where we use the same parameters in the Hamiltonian as in the previous section (, , , and ). At zero temperature we only have to consider matrix elements between the RVB ground state of the Rokhsar-Kivelson model and eigenstates of (2) with one fermionic dimer that belong to the zero winding number sector . At the RVB ground state is an equal amplitude superposition of all dimer configurations and has vanishing energy. Again, we used twisted boundary conditions to access arbitrary momenta within the first Brillouin zone.

The top left panel shows at the nodal point, which has a sharp peak at the Fermi surface () independent of system size, as expected from our ground state computations. The incoherent part of the spectrum, containing of the hole spectral weight, is broadly distributed with a maximum located around . The data clearly shows that the spectral weight of the central peak vanishes and a gap opens as we go away from the nodal point towards the antinode, where the spectral weight is redistributed to lower frequencies. In the antinodal region around , the spectral function exhibits a sizable pseudogap on the order of , independent of system size (see Fig. 3 bottom right). The complete spectral weight is now distributed over a broad peak of width centered around .

In Fig. 4 we plot the pseudogap as function of the angle between the antinode () and the node () for a fixed distance from , extracted from our numerical data for a lattice size of sites. We emphasize again that corresponds to the onset of the spectral function at low energies and represents a lower bound for the pseudogap extracted from the spectral function, because the low energy states in the spectrum turn out to have a vanishingly small weight in the antinodal region, which leads to an apparently larger gap in the spectral function. This effect is only seen at zero temperature, however. Nonetheless, the dimer model features a sizeable pseudogap on the order of at the antinode. It is also important to note that the pseudogap shows a clear deviation from a simple d-wave form (see inset in Fig. 4) and is thus clearly distinguishable from the superconducting gap. This is in agreement with a wide range of experiments that found evidence for corrections to the d-wave symmetry for the pseudogap in the underdoped regime, in stark contrast to the superconducting gap Lee et al. (2007); Vishik et al. (2010); Loret et al. (2017).

Finally, in Fig. 5 we show finite temperature results for the hole spectral function in the antinodal region. At finite temperature matrix elements involving excited states of the Rokshar-Kivelson model contribute to the spectrum as well, which leads to a thermal broadening of the incoherent peak and a shift of spectral weight to lower energies. More interestingly, the finite temperature results show that the low energy states do contribute finite weight to the spectral function at the antinode, and the pseudogap is indeed given by shown in Fig. 4. This is in contrast to the zero temperature case, where low energy states have a vanishingly small weight and the pseudogap appears to be larger in the spectral function.

## Iv Two-mode approximation

In this section we present a semi-analytic description of the hole spectral function in terms of a two-mode approximation, where we assume that the low energy part of the spectrum of excited states can be captured by states of the form

 |p⟩k =~F†k+p,1~Dp,1|GS⟩ (9)

where the fermionic (bosonic) operator () is a linear combination of () and () and describes the lower of the two fermionic (bosonic) bands. denotes the ground state of (2) at a given doping. In order to determine the matrix elements relating and we use a mean-field approximation to diagonalize the Hamiltonian (2), as outlined below. A similar mean-field approach has been used previously in Ref. Goldstein et al. (2017).

It is important to note here that the two-mode approximation captures the incoherent part of the hole spectral function, but is not able to reproduce the coherent peak, as we explain in detail below. Moreover, the mean-field approach violates the hard-core constraint for the dimers, which is a very crude approximation. Nevertheless, the main features of the pseudogap are correctly reproduced when approaching the antinodal region in momentum space.

We start by outlining the mean field approximation for the description of a small, but finite density of fermionic dimers interacting with a background of bosonic dimers forming an RVB spin liquid. In total we need five real, homogeneous and isotropic fields to decouple the dimer Hamiltonian in all possible bosonic channels

 χsz=⟨D†i,ηDj,η′⟩, (10)

where denotes on-site terms (), whereas are nearest-neighbor terms (). The index denotes parallel () or perpendicular () dimer configurations.
We demand that the mean-field solutions represent a state with no broken symmetries, thus must be invariant under the symmetry group of the square lattice. From this follows that the mean-fields and have to be equal. The RK mean-field Hamiltonian thus takes the form

 Unknown environment '% (11)

Similarly, the interaction term between fermionic and bosonic dimers can be decoupled using the mean-fields (10) and reads

 HMF1=∑q(F†q,x,F†q,y)(−2t1χ1||cos(qy)−t2χ1⊥h2,q−t3χ2⊥h3,q−t2χ1⊥¯h2,q−t3χ2⊥¯h3,q−2t1χ1||cos(qx))(Fq,xFq,y), (12)

with

 h2,q =(1+eiqx+e−iqy+ei(qx−qy)), (13) h3,q =(ei(qy−2qx)+eiqx+e−i2qx+ei(qx+qy)+ei2qy+ei(2qy−qx)+e−iqy+e−i(qx+qy)). (14)

We diagonalize the mean-field Hamiltonian using a Bogoliubov transformation and , which provides the matrix elements relating and , as mentioned above. The single particle operators and describe excitations with dispersion relations , respectively , with band index .

Assuming that the electronic spectrum is well approximated by excitations of the form (9), the corresponding eigenenergies of the excitations are then given by

 HMF|p⟩k =[ξ~Fk+p,1−ξ~Dp,1+ξGS]|p⟩k,

where and .

The hole part of the spectral function at zero temperature in this mean-field approximation reads

 A−,k(ω) =∑pQk(p) δ(ω+ξ~Fk+p,1−ξ~Dp,1+ξGS), (15)

where

 Qk(p) =[fk(p)(1−nF(ξ~Fk+p,1))nB(ξ~Dp,1)]2, fk(p) =1√2NxNy[¯N11p+kM11pf^x(k)+¯N21p+kM21pf^y(k)].

Note that the spectral function in the two-mode approximation cannot have a sharp, coherent delta-function peak, because Eq. (15) always involves an integral over momenta. We use the numerical data of the dispersion in Fig. 2 to constrain the fitting of the mean field dispersions and by using the relation , which determines the onset of the spectrum. Furthermore we approximate the momentum distribution of fermionic dimers by a simple Fermi-Dirac distribution, which is appropriate within the mean-field description. For the bosonic dimers a Bose-Einstein distribution would not be appropriate, however, because it does not capture the important hard-core constraint of bosonic dimers. For this reason we take the bosonic dimer distribution to be a constant, as expected for hard core bosons in a semi-classcial limit Coletta et al. (2012).

In Fig. 3 we plot the spectral functions together with the two-mode approximation for different momenta. Even though this appraoch is not capable to reproduce the coherent peak of the spectral function, we find a finite weight at the chemical potential. Upon approaching the antinodal region the pseudogap slowly opens as function of momentum and has a similar size as in the ED results (see Fig. 3 bottom left). The semi-analytic mean field approach describes the spectral function at the antinode reasonably well, where it shows a clear siginature of the pseudogap with a dominant incoherent peak at (see Fig. 3 bottom right).

In conclusion, the two-mode approximation is able to repoduce some common freatures of the spectral function. Its drawback is that it doesn’t capture the coherent quasiparticle peak, however. Such a coherent peak would appear in the two-mode approximation only if we allow for a boson condensate. Since there is no physical basis for the appearance of a condensate in a hard-core boson system at integer filling, we refrained from such a modification. Another modification would be the inclusion of corrections to the boson distribution beyond the semiclassical limit, as discussed in Ref. Coletta et al. (2012). These corrections give rise to a term in the momentum distribution of the bosonic dimers which leads to the appearance of an additional peak in the spectral function at high energies around , but don’t change the qualitative behavior at lower energies, in particular the onset of the pseudogap.

## V Diagrammatic Results

In this section we present a systematic diagrammatic approach to compute the electron spectral function and the coherent quasiparticle residuum in particular. Since the Hamiltonian of the dimer model does not feature a quadratic part, there formally does not exist a small parameter in the system that would rigorously justify the use of such perturbative means. We pursue this approach nonetheless and for find good agreement with numerical results. The essential ingredient is the expansion around the exactly solvable RK-point Rokhsar and Kivelson (1988), where dimer-dimer correlations take a classical form, see e.g. Fisher and Stephenson (1963).

Within the domain , fermionic dimers are inserted into a RK-like bosonic background which most importantly allows to consider the ground state to be translationally invariant. We expect the coherent part of the spectral function to be induced by the part of the dimer Hamiltonian (2) in which we will conduct a perturbative expansion starting from the action of the model in momentum space

 S= ∑q1,q2,q3,q4{H[¯F(q1),F(q4),¯D(q2),D(q3)]}+ (16) +∑q,η,α¯Fq,η,α(iω)(−iω−μf)Fq,η,α(iω)+ +∑q,η¯Dq,η(iω)(−iω−μb)Dq,η(iω),

where in the notation momentum and energy conservation is understood. The fields , correspond to anticommuting Grassmann variables while , are complex fields. Since symmetry is manifest, we drop the spin index for the fermionic fields in the following. As a further approximation we do not enforce the hard-core dimer constraint exactly for the moment, but introduce chemical potentials and to fix the fermionic and bosonic dimer densities on average (we comment on how to enforce the hard-core constraint below Eq. (24)). Since the individual terms in the Hamiltonian obey the hard-core constraint locally, this turns out to be a good approximation. The bare dimer propagators thus are

 G0f/b(iω)=1iω+μf/b, (17)

with doping dependent chemical potentials fixing the average dimer densities on a given link of the lattice to and . Within the Matsubara formalism we can compute the electronic spectral function from the electronic imaginary time ordered Green’s function via analytic continuation . In the dimer Hilbert space (see Eq. (1)), can be obtained from the relation

 Gp(iωp)= 14βN∑η1,η2ϵ{x,y}(1+e−ipη1)(1+eipη2) (18) × ∑q1,q2⟨Fq2,η2Dq1+p,η1¯Dq2+p,η2¯Fq1,η1⟩,

where the average is to be evaluated in the framework of the action from Eq. (16) ( is the number of sites). The corresponding spectral function for holes is then given by .

We can write the four point dimer correlator in full generality as

 ⟨Fq2,η2Dq1+p,η1¯Dq2+p,η2¯Fq1,η1⟩= δq1,q2δη1,η2Gf(q1)Gb(q1+p)+ (19) + Gf(q2)Gb(q1+p)~Γη1,η2,η1,η2(q1,q2+p,q1+p,q2)Gb(q2+p)Gf(q1),

where we have introduced the full interaction vertex and the full (dressed) propagators and . It is important to realize that the full propagators can not develop a dispersion at any order in perturbation theory, as the interaction terms in the Hamiltonian locally respect the hard-core constraint. Any diagram with only two external lines thus necessarily has propagators on the same lattice site for the incoming and outgoing lines. Typically, we expect the full propagators to contain finite lifetimes for the two dimer species which merely lead to a Lorentzian broadening of delta-type peaks in the spectral function. We thus approximate in Eq. (19). The first term in Eq. (19) corresponds to the zeroth order contribution and upon insertion into Eq. (18) results in the electronic spectral function

 A0p(ωp)=14(cos2(px2)+cos2(py2))⋅2πδ(ωp) (20)

for the non-interaction limit . Here, , which leads to the peak position at , was used. It can be verified that this expression indeed reproduces the exact quasiparticle residuum at this specific point in parameter space. We have thus shown that our expansion around an RK-like background at reproduces the correct result right at the expansion point as required for a meaningful Ansatz.

We now seek to find a good approximation to the full vertex by considering particle-hole like ladder diagrams. These turn out to correspond to repeated exchange hoppings of a fermionic dimer in the bosonic background and are thus expected to yield good estimates for considering exchange terms only. The effective vertex is then determined by solving a Bethe-Salpter equation which is displayed diagrammatically in Fig. 6. For the vertex from Eq. (19) this equation reads

 ~Γη1η2η1η2(q1,q2+p,q1+p,q2)= Γη1η2η1η2(q1,q2+p,q1+p,q2)+ (21) +∑~ηf,~ηbβψ(iωp){∑~qΓη1~ηbη1~ηf(q1,~q+p,q1+p,~q)~Γ~ηfη2~ηbη2(~q,q2+p,~q+p,q2)},

where corresponds to the bare vertex (see Appendix) which is derived from the Hamiltonian, while the particle-hole bubble corresponds to the Matsubara sum of an antiparallel pair of bare fermionic and bosonic propagators.

To solve this integral equation for , we note that in principle only the exchange interactions and contribute to the ladder diagrams of Eq. (21). This is due to the external dimer orientations in Eq. (21) being fixed such that the orientation of the ingoing (outgoing) fermionic (bosonic) dimer () matches the orientation of the outgoing (ingoing) bosonic (fermionic) dimer. A flip term would induce a different relative orientation between ingoing and outgoing dimers and thus cannot contribute to the particle-hole ladder. However, we can include the -terms in an effective manner by substituting them with the following exchange terms,

 Ht2→~Ht2=−t2∑iD†i,yF†i+^y,x,αFi,y,αDi+^y,x+7 terms. (22)

Even though this term explicitly violates the hard-core constraint, the physical reasoning for this replacement is the following: we expect the major contributions by -processes to come via induced exchange interactions, i.e. a -flip combines with a - or - exchange and a subsequent bosonic -plaquette flip which restores the starting configuration of the bosonic background. Such a sequence effectively acts as an exchange. Since our ladder approach does not include the purely bosonic background interactions, we account for these by assigning to an effective first order exchange interaction that induces the correct exchange interactions at higher order. Note that although the terms in Eq. (22) project on non-constraint configurations, we still expect this substitution to be valid as we consider a translational invariant dimer density in our approach. Comparing our diagrammatic results to numerical data shows that this approximation indeed works very well.

For exchange terms only, the Bethe-Salpeter equation can be solved exactly as the bare vertex only depends on the electronic momentum . We can hence assume and the sum over in Eq. (21) can be evaluated trivially.

Eq. (21) can then be brought into simple matrix form by defining the vector

 →Γ(p)≡(Γxxxx(p),Γxyxy(p),Γyxyx(p),Γyyyy(p)) (23)

(and likewise for ) and writing with

 M(p)=⎛⎜ ⎜ ⎜⎝1−βNψΓxxxx0−βNψΓxyxy001−βNψΓxxxx0−βNψΓxyxy−βNψΓyxyx01−βNψΓyyyy00−βNψΓyxyx01−βNψΓyyyy⎞⎟ ⎟ ⎟⎠. (24)

Straightforward matrix inversion then yields the effective vertex.

We first shortly discuss how to implement the hard-core constraint beyond the mean field level of chemical potentials. To this end we note that for every -rung real space particle-hole ladder diagram, antiparallel bosonic and fermionic propagator lines corresponds to the same link on the lattice. This is due to ingoing bosonic links matching outgoing fermionic ones and vice versa on every interaction line in exchange processes. As the particle-hole bubble is proportional to the total average dimer density on the corresponding link, the contribution of an -rung process is proportional to the product over the total dimer densities of all links occuring in the process. Because we expand around the RK-point, we may substitute this product with the classical probability of having all the links occuring in the process occupied, denoted by . These correlations can be computed from a Grassmann field theory for the classical dimer problem Samuel (1980). Attention has to be paid only for interaction lines corresponding to , where the classical probability of having both occuring links occupied would vanish. Instead, we consider the probability of having the two dimers in a relative position which corresponds to the original -flip. By this reasoning, the classical dimer correlations can in principle be implemented exactly into our approach. We can achieve a first order approximation by the simple replacement , where correspond to the displacement vector and relative orientation change in a single process. In analogy to Punk et al. (2015) this leads to the simple replacements , , that have to be made in all our diagrammatic results presented below.

Using the effective vertex from above and introducing a finite peak width via yields the electronic spectral function

 Ap(ωp)=4γ8Z0p(16(γ2+ω2p)+K1(p))+2K2(p)(4ωp−t1cos(px)−t1cos(py))[16(γ2−ω2p)+8ωpt1(cos(px)+cos(py))+K1(p)]2+[8γ(4ωp−t1cos(px)−t1cos(py))]2, (25)

which assumes the form of a sum of two Lorentzians . The functions and are given in the Appendix while is the residuum at . The peak positions of Eq.(25) are given by the singularities in the limit and turn out to be

 ω1,2= 14t1(cos(px)+cos(py)) (26) ±14√t21(cos(px)+cos(py))2+K1(p).

Note that the spectral function (25) only has two sharp peaks and no incoherent background due to the fact that the dimer propagators cannot obtain a dispersion at any order in perturbation theory. Consequently, the particle-hole bubble of a bosonic and a fermionic dimer has a simple pole on the real frequency axis, which gives rise to the simple two-peak structure in the spectral function.

As the insertion of holes corresponds to the removal of the highest band electrons, we can relate the dispersion and residuum for the holes with the rightmost peak of Eq. (25) to obtain

 εp= −ω1(p) (27)
 Zp=Z1,p= limγ→012γAp(ω1(p)) (28)

for hole-dispersion and -residuum. Examples for different values of are compared to ED results in Fig. 7. For the chosen parameter regime the assumption of an RK-background is valid and yields good agreement with ED results. This includes parameter sets with non-vanishing -interactions, which were treated via the effective exchange interaction of Eq. (22). We note that the hole dispersion of Eq. (27) reproduces exactly the dispersion obtained in Punk et al. (2015) by a perturbative Ansatz. Our approach can hence be viewed as an extension of this Ansatz which additionally yields the quasiparticle residuum.

## Vi Discussion and Conclusions

Our numerical results show that the dimer model introduced in Ref. Punk et al. (2015) has a sizeable pseudogap in the antinodal region of the Brillouin zone close to and symmetry related momenta. Moreover, it’s momentum dependence clearly deviates from a simple d-wave form, in accordance with observations in the pseudogap regime of underdoped cuprates.

It is important to emphasize here that we always fine-tune our model to the RK-point at , where the ground-state of the undoped model is a spin liquid. Away from this point the spin liquid is confining and unstable towards symmetry broken valence bond solid states. Upon doping away from the RK-point the dimer model (2) thus realizes a -FL*, which again is expected to be confining at long length-scales. This problem can be circumvented by allowing for diagonal dimers between next-nearest neighbor sites as well. In this case the undoped RK-model has a stable spin liquid ground-state in an extended parameter regime Moessner and Sondhi (2001). Accordingly, upon doping the ground-state of the appropriately modified model in Eq. (2) is expected to be a -FL*, which is a stable phase of matter Patel et al. (2016). Nevertheless, we don’t expect a qualitatively different behavior of the single-electron spectral function when including diagonal dimers. In particular, the pseudogap in the antinodal region is expected to be a robust property of the dimer model. We leave a detailed analysis of this problem open for further study.

It’s also worthwile to contrast our results with numerical dynamical cluster approximation (DCA) studies of the Hubbard model on the square lattice, where a sizeable pseudogap in the spectral function at the antinode was found as well Gull et al. (2010, 2013). We point out here that two site cluster dynamical mean field theory (DMFT) studies of the Hubbard model indeed showed that the electron configuration on the two site cluster at low dopings is dominated by the same two states that are used here to span the Hilbert space of the dimer model on the full lattice, i.e. the bosonic and the fermionic dimer Ferrero et al. (2009). This suggests that the dimer model introduced in Ref. Punk et al. (2015) provides an effective low energy description of the square lattice Hubbard model at large and low doping.

###### Acknowledgements.
We thank D. Pimenov and F. Dorfner for valuable discussions and S. Sachdev for comments on the manuscript. This research was supported by the German Excellence Initiative via the Nanosystems Initiative Munich (NIM).

## Appendix A Diagrammatic Approach

We provide some additional information on the ladder approach of Sec.(V). The bare vertex from the Bethe-Salpeter Eq. (21) is determined from the momentum space form of the dimer Hamiltonian , including . The relevant bare vertices in Eq. (23) and (24) are

 Γxxxx(p)= 2t1βNcos(py), (29) Γyyyy(p)= 2t1βNcos(px), (30) Γxyxy(p)≡ C(p)= (31) = t3βN(ei(px+py)+eipx+ei(py−2px)+e−2ipx +e2ipy+ei(2py−px)+e−ipy+e−i(px+py))+ +t2βN(1+eipx)(1+e−ipy), Γyxyx(p)= C∗(p)=C(−p). (32)

The functions and which contribute to the spectral function from Eq. (25) are given by

 K1(p)=t23[8+4cos(px)+4cos(3px)+2cos(px−3py)+ (33) +4cos(px−2py)+4cos(2px−2py)+4cos(px−py)+ +4cos(2px−py)+2cos(3px−py)+4cos(py)+ +4cos(3py)+4cos(px+py)+4cos(2px+2py)+ +4cos(2px+py)+2cos(3px+py)+ +4cos(px+2py)+2cos(px+3py)] −t21[2cos(px−py)+2cos(px+py)]+
 +t22[4+4cos(px)+4cos(py)+ +2cos(px−py)+2cos(px+py)]+ +t2t3[4+6cos(px)+6cos(py)+4cos(2px)+ +4cos(2py)+2cos(3px)+2cos(3py)+4cos(px−3py)+ +4cos(3px−py)+2cos(2px−3py)+2cos(3px−2py)+ 4cos(px−2py)+4cos(2px−py)+4cos(2px−2py)+ Missing or unrecognized delimiter for \biggr K2(p)=t3[2+3cos(px)+2cos(2px)+cos(3</