Doublon-holon origin of the subpeaks at the Hubbard band edges

TL;DR

This paper explains the origin of subpeaks at the Hubbard band edges in the spectral function of the SU(2) Fermi-Hubbard model near the Mott transition, attributing them to doublon-holon pair interactions, supported by DMFT and mean-field analysis.

Contribution

It identifies the doublon-holon pairs as the origin of subpeaks and demonstrates this through combined DMFT and mean-field approaches, extending understanding to SU(3) and SU(4) models.

Findings

Subpeaks originate from doublon-holon pair interactions.

Mean-field results agree with DMFT calculations.

Subpeaks are more pronounced in SU(3) and SU(4) models.

Abstract

Dynamical mean-field theory (DMFT) studies frequently observe a fine structure in the local spectral function of the SU(2) Fermi-Hubbard model at half filling: In the metallic phase close to the Mott transition, subpeaks emerge at the inner edges of the Hubbard bands. Here we demonstrate that these subpeaks originate from the low-energy effective interaction of doublon-holon pairs, by investigating how the correlation functions of doublon and holon operators contribute to the subpeaks. A mean-field analysis of the low-energy effective Hamiltonian provides results consistent with our DMFT calculation using the numerical renormalization group as an impurity solver. In the SU(3) and SU(4) Hubbard models, the subpeaks become more pronounced due to the increased degeneracy of doublon-holon pair excitations.

Figures

Equations

$$d_{i\nu }^{†}\equiv P_{i,\bar{n}+1}{c}_{i\nu }^{†},h_{i\nu }^{†}\equiv P_{i,\bar{n}-1}{c}_{i\nu },$$

$$H_{\mathrm{eff}}=H_{v;0}+H_{ss}+H_{dh}+H_{\text{3-site}}+O(v^3/U^2),$$

$$H_{ss}=\frac{v^2}{U}\sum _{⟨i,j⟩}4{\vec{S}}_i\cdot {\vec{S}}_j-P_{i1}{P}_{j1},$$

$$H_{dh}=\frac{2v^2}{U}\sum _{⟨i,j⟩}(c_{j1}^{†}{c}_{j2}^{†}{c}_{i2}{c}_{i1}+P_{i2}{P}_{j0})+(i\leftrightarrow j)$$

$$=\frac{v^2}{U}\sum _{⟨i,j⟩,\nu ,{\nu }^{\prime }}(h_{i\nu }^{†}{d}_{j\nu }^{†}+h_{j\nu }^{†}{d}_{i\nu }^{†})(d_{i{\nu }^{\prime }}{h}_{j{\nu }^{\prime }}+d_{j{\nu }^{\prime }}{h}_{i{\nu }^{\prime }}),$$

$$A_{\mathrm{disc}}^z(\omega )=\sum _{ij}{A}_{ij}^z\delta (\omega -{\omega }_{ij}^z),$$

$$A_{ij}^z\equiv ⟨E_i^z|O_1|E_j^z⟩⟨E_j^z|O_2|E_i^z⟩({\rho }_i^z\pm {\rho }_j^z),$$

$${\omega }_{ij}^z\equiv E_j^z-E_i^z,$$

$$\delta (\omega -{\omega }^{\prime })\Rightarrow {\tilde{\delta }}_{\sigma }(\omega ;{\omega }^{\prime })\equiv \frac{\Theta (\omega {\omega }^{\prime })}{\sqrt{\pi }\sigma |{\omega }^{\prime }|}{e}^{-{\left(\frac{\ln |\omega /{\omega }^{\prime }|}{\sigma }-\frac{\sigma }{4}\right)}^2},$$

$${\sigma }_{ij}=\frac{\alpha }{n_z}\frac{d\ln |{\omega }_{ij}^z|}{dz}=\frac{\alpha }{n_z|{\omega }_{ij}^z|}\left|\frac{d{E}_j^z}{dz}-\frac{d{E}_i^z}{dz}\right|,$$

$$A(\omega ,{ϵ}_{\mathbf{k}})=\frac{-1}{\pi }\mathrm{Im}\frac{1}{\omega +\mu -{ϵ}_{\mathbf{k}}-\Sigma (\omega )+i{0}^{+}},$$

$$A(\omega )=\int \mathrm{d}ϵ{\rho }_0(ϵ)A(\omega ,ϵ),$$

Full text

Introduction.— Dynamical mean-field theory (DMFT) Georges et al. (1996); *Kotliar2006 provides a widely successful approach in understanding strongly correlated systems. It treats a lattice problem by self-consistently solving an effective impurity model whose impurity and bath correspond to a lattice site and the rest of the lattice, respectively. Thus the performance of DMFT calculations directly depends on which particular impurity solver is chosen.

A benchmark calculation for various impurity solvers is the paramagnetic Mott transition in the half-filled $\mathrm{SU}(2)$ Hubbard model at temperature $T=0$ which is characterized by a striking change in the local spectral functions Zhang et al. (1993); Bulla (1999). In the metallic phase, the spectral function features a quasiparticle peak (QP) at the Fermi level, and two Hubbard bands (HBs) below and above the Fermi level each. In the insulating phase, the QP disappears and a gap opens between two HBs.

In the metallic phase close to the transition, many DMFT studies have observed sharp subpeaks that emerge at the inner edges of the HBs, by using different real-frequency impurity solvers: perturbative methods Zhang et al. (1993), the density-matrix renormalization group (DMRG) Karski et al. (2005, 2008); Ganahl et al. (2014); *Ganahl2015; Wolf et al. (2014), the numerical renormalization group (NRG) Žitko and Pruschke (2009), and exact diagonalization Granath and Strand (2012); Lu et al. (2014). In contrast, quantum Monte Carlo solvers, which obtain the spectral functions on the real frequency axis via (numerically ill-posed) analytic continuation, have not found these subpeaks. The subpeaks give rise to distinct features in the momentum-resolved spectral function Karski et al. (2008), measurable by photoemission spectroscopy Mo et al. (2003); Sekiyama et al. (2004). Despite these frequent consistent observations, the physical origin of the subpeaks and their relevance in more general (e.g., multi-flavor) models remained unclear.

In this Letter, we show that the subpeaks are induced by the effective doublon-holon (DH) Yokoyama et al. (2006); *Phillips2010; *Sato2014; *Zhou2014; *Prelovsek2015; *Zhou2016 pair interaction originating from a second-order virtual process, where a doublon (holon) means an excitation that one particle is added to (removed from) a lattice site with average integer filling. We compute the correlation functions of doublon and holon operators in the $\mathrm{SU}(2)$ Hubbard model, by using DMFT with NRG Wilson (1975); Bulla et al. (2008) as an impurity solver, and demonstrate that these correlation functions manifest the peak structure associated with the subpeaks. We reproduce the peak structure of doublon and holon correlators via a mean-field analysis of the low-energy effective Hamiltonian obtained by a generalized Schrieffer-Wolff transformation (SWT) Bukov et al. (2016); *Bukov2015; Lee et al. (2017). Both approaches consistently result in a linear dependence of the subpeak position vs. interaction strength. From our DMFT+NRG calculations of general $\mathrm{SU}(N)$ Hubbard models for $N=2,3,4$, we observe that the subpeaks become more pronounced with increasing $N$, since the DH pair excitations become more degenerate due to the larger $\mathrm{SU}(N)$ symmetry.

System.— The $\mathrm{SU}(N)$ Hubbard model describes $N$ flavors of fermions on a lattice with local repulsive interactions, recently realized in ultracold atom experiments with tunable $N$ Taie et al. (2012); *Hofrichter2016. The hopping amplitude $v$, the interaction strength $U$, and the chemical potential are flavor-independent, thus the system has $\mathrm{SU}(N)$ flavor symmetry. Its Hamiltonian is $H=H_{U}+H_{v}+H_{\mu}$, where $H_{U}=\frac{U}{2}\sum_{i}(\hat{n}_{i}-\bar{n})^{2}$, $H_{v}=v\sum_{\langle i,j\rangle,\nu}c_{i\nu}^{\dagger}c_{j\nu}+\text{h.c.}$, and $H_{\mu}=-\mu\sum_{i}\hat{n}_{i}$. Here $c_{i\nu}$ annihilates a particle of flavor $\nu=1,\ldots,N$ at lattice site $i$, $\hat{n}_{i}=\sum_{\nu}c_{i\nu}^{\dagger}c_{i\nu}$ is the particle number operator at site $i$, $\langle i,j\rangle$ indicates nearest neighbours, $\bar{n}$ is a parameter for the desired average occupation, and $\mu$ is a fine tuning of chemical potential to achieve $\langle\hat{n}_{i}\rangle=\bar{n}$. Throughout this paper, we focus on $T=0$ and the average occupation number as an integer closest to half filling $\bar{n}=\lfloor N/2\rfloor$, by fixing $\mu=0$ for $N=2,4$, and fine-tuning $\mu$ for $N=3$.

Doublon and holon.— For integer average occupation $\bar{n}$, we define doublon and holon creation operators as

[TABLE]

where $P_{in}$ means the projector onto the subspace in which the site $i$ has $n$ particles. For the $\mathrm{SU}(2)$ case, at half filling, $\bar{n}=1$, these operators reduce to $d_{i\nu}=c_{i\nu}\hat{n}_{i\bar{\nu}}$, $h_{i\nu}=c_{i\nu}^{\dagger}(1-\hat{n}_{i\bar{\nu}})$ with $\hat{n}_{i\nu}=c_{i\nu}^{\dagger}c_{i\nu}$ and $\bar{\nu}=3-\nu$, and they completely constitute the particle operator $c_{i\nu}=d_{i\nu}+h_{i\nu}^{\dagger}$. Then the particle correlator can be decomposed into four doublon and holon correlators, $A_{cc^{\dagger}}(\omega)=A_{dd^{\dagger}}+A_{dh}+A_{h^{\dagger}d^{\dagger}}+A_{h^{\dagger}h}$, where $A_{XY}(\omega)\equiv\tfrac{-1}{\pi}\mathrm{Im}G_{XY}$, with $G_{XY}(t)=-i\vartheta(t)\langle[X(t),Y(0)]_{\pm}\rangle_{T}$ being the retarded correlation function of the fermionic ($+$) or bosonic ($-$) local operators $X$ and $Y$ acting on the same site. In the particle-hole symmetric case, only two correlators are independent: “diagonal” correlators $A_{dd^{\dagger}}(\omega)=A_{h^{\dagger}h}(-\omega)$ which are asymmetric, and “off-diagonal” correlators $A_{h^{\dagger}d^{\dagger}}(\omega)=A_{dh}(\omega)$ which are symmetric under $\omega\leftrightarrow-\omega$. For $N>2$ flavors, the decomposition of $c_{i\nu}$ acquires more terms than $d_{i\nu}$ and $h_{i\nu}^{\dagger}$ Lee et al. (2017).

DMFT+NRG.— We use single-site DMFT which maps the Hubbard model onto the single-impurity Anderson model (SIAM) which provides paramagnetic solutions, by construction. We employ the semi-circular density of states of the Bethe lattice with half-bandwidth $D$, together with units $D=\hbar=k_{B}=1$, throughout. We solve the SIAM by the full-density-matrix NRG (fdm-NRG; Weichselbaum and von Delft (2007); *Weichselbaum2012:mps), exploiting $\mathrm{U}(1)_{\mathrm{charge}}\otimes\mathrm{SU}(N)_{\mathrm{flavor}}$ symmetry Weichselbaum (2012b). The coarse-grained discretization-averaged spectral data is broadened adaptively Sup ; Lee and Weichselbaum (2016) for best possible spectral resolution at higher energies, while preserving the intrinsic accuracy of NRG at low energies [e.g., the Luttinger pinning Müller-Hartmann (1989) $\tfrac{\pi}{2}A(\omega=0)=1$ in the metallic phase is accurately satisfied; see Fig. 1(a)-(b)].

SU(2) metallic phase.— We first consider the case $N=2$ equivalent to the spin-full one-band Hubbard model. At $T=0$ and half filling, a metallic phase exists for $U<U_{c2}=2.91(1)$, and a paramagnetic insulating phase for $U>U_{c1}=2.37(2)$. For $U_{c1}<U<U_{c2}$ the two phases coexist (e.g., see Fig. 2, or Refs. Bulla (1999); Bulla et al. (2001)).

Within the metallic phase, the local spectral function $A(\omega)$ features one QP and two HBs [cf. Fig. 1(a)-(b)]. As $U$ increases, the central QP narrows, the HBs widen, and the dips between the QP and the HBs deepen. On top of this, subpeaks are present at the inner edges of the HBs, whose position $\omega_{p}$ and width $\delta\omega$ decrease linearly with increasing $U$, as shown in Fig. 2.

Local spin (i.e., flavor) and charge susceptibilities, $\chi_{s}$ and $\chi_{c}$ Raas and Uhrig (2009), in Fig. 1 demonstrate that the QP and the HBs of $A(\omega)$ are tied to spin and charge degrees of freedom, respectively; that is, spin and charge excitations are energetically separated. The peak of $\chi_{s}$ indicates a spin-like collective mode responsible for the QP, which is analoguous to the Kondo resonance in the SIAM in that the spin susceptibility peaks at the Kondo energy scale Hanl and Weichselbaum (2014). The position $\omega_{s}$ and width of the $\chi_{s}$ peak decrease as the QP narrows with increasing $U$; especially, $\omega_{s}$ has a linear dependence vs. $1/U$, as shown in Fig. 2. In contrast, $\chi_{c}$ is suppressed within the QP region, while having long tails beyond the outer edges of the HBs.

For $T=0^{+}$, the positive and negative energy sides of a correlator $A_{XY}(\omega)$ are derived from $\langle X(t)Y(0)\rangle_{T}$ and $\langle Y(0)X(t)\rangle_{T}=\langle X^{\dagger}(t)Y^{\dagger}(0)\rangle^{*}_{T}$, respectively. Therefore the upper HB in Fig. 1, which mainly consists of $A_{dd^{\dagger}}$, originates from the dynamics of the doublon $d_{i\nu}^{\dagger}(0)$. Another significant feature of $A_{dd^{\dagger}}$ is a peak at $\omega=-\omega_{s}$. Just after the action of $d_{i\nu}(0)$ and just before $d_{i\nu}^{\dagger}(t)$, the site $i$ has only spin-$\bar{\nu}$. Its time evolution between [math] and $t$ with low frequency $|\omega|\simeq\omega_{s}$ is driven by the spin-like collective mode captured by the peak of $\chi_{s}$ at $\omega_{s}$. In contrast, the off-diagonal correlator $A_{dh}$ has a symmetric peak at $\omega=0$. This reflects the particle-hole symmetric processes of destroying at the same site first a doublon and then a holon, or vice versa. $A_{dd^{\dagger}}$ and $A_{dh}$ contribute comparably to the QP, having $A_{dd^{\dagger}}(0)=A_{dh}(0)=\tfrac{1}{2\pi}$.

In the metallic regime in Figs. 1(a)-(b) all of the doublon and holon correlators show peak-like features at $\omega=\pm\omega_{p}$. For $U\gtrsim 2.3$ Sup , their contributions to these subpeaks have relative weights $A_{dd^{\dagger}}(\omega_{p})>A_{dh}(\pm\omega_{p})>A_{dd^{\dagger}}(-\omega_{p})$. Our effective theory (described below) aims to reproduce this relative order of contributions, as well as the linear dependence of $\omega_{p}$ vs. $U$.

SU(2) insulating phase.— The QP, the subpeaks, the spin-charge separation in energy space, and the peaks of the doublon and holon correlators all disappear in the insulating phase, as depicted in Fig. 1(c). Instead, a Mott gap $\Delta$ opens, and the susceptibilities $\chi_{s}$ and $\chi_{c}$ spread over a large energy range, $|\omega|>\Delta/2$, with suppressed heights. While both $\omega_{p}$ in the metallic phase and $\Delta/2$ in the insulating phase correlate to the location of the inner HB edges, their dependences on $U$ are clearly different (see Fig. 2). Here the absence of subpeaks is consistent with previous studies Karski et al. (2005, 2008); Žitko and Pruschke (2009); Granath and Strand (2012); Lu et al. (2014); Ganahl et al. (2014); Wolf et al. (2014); Ganahl et al. (2015). Though other works Nishimoto et al. (2004); Gull et al. (2010); Granath and Schött (2014) have reported subpeaks even in the insulating phase, their observations are not numerically stable due to, e.g., ill-posed analytic continuation or underbroadening.

DH pair interaction.— We will now demonstrate that the peaks of the doublon and holon correlators at $\omega=\pm\omega_{p}$, which add up to the subpeaks of $A(\omega)$, originate from a DH pair interaction within the low-energy effective Hamiltonian of the $\mathrm{SU}(2)$ Hubbard model. Our theory is based on the separation of three energy scales, $\omega_{s}<\omega_{p}<U/2$, corresponding to the QP, the subpeaks, and the HBs, respectively. We focus on the intermediate scale $\omega_{p}$ by integrating out the larger scale $U/2$ and by approximating the physics of the smaller scale $\omega_{s}$.

We first integrate out the charge fluctuation of energy scale $U/2$, by employing a generalized SWT Bukov et al. (2016); *Bukov2015; Lee et al. (2017). We decompose the hopping term into different components $H_{v}=\sum_{m=-1}^{+1}H_{v;m}$ which cost Coulomb energy $mU$ since $mUH_{v;m}=[H_{U},H_{v;m}]$. Here $H_{v;0}\equiv v\sum_{\langle i,j\rangle,\nu}(d_{i\nu}^{\dagger}d_{j\nu}-h_{i\nu}^{\dagger}h_{j\nu})+\text{h.c.}$ describes the hopping of doublons and holons without energy cost, whereas $H_{v;1}\equiv v\sum_{\langle i,j\rangle,\nu}(d_{i\nu}^{\dagger}h_{j\nu}^{\dagger}+d_{j\nu}^{\dagger}h_{i\nu}^{\dagger})$ or ($H_{v;-1}=H_{v;1}^{\dagger}$) creates (annihilates) nearest-neighbor DH pairs by paying (gaining) energy cost $U$. Then we write the low-energy effective Hamiltonian $H_{\mathrm{eff}}$ as a power series in $v/U$,

[TABLE]

where $H_{\text{3-site}}$ is the sum of the products of operators at three nearest neighbor sites. The term $H_{ss}+H_{dh}+H_{\text{3-site}}=[H_{v;1},H_{v;-1}]/U$, of order $O(v^{2}/U)$, can be interpreted as second-order virtual processes. $H_{\mathrm{eff}}$ is similar to the $t$-$J$ model Harris and Lange (1967); *Chao1977; *MacDonald1988; *Eskes1994; *Eskes1994a; *Chernyshev2004, widely used as the effective low-energy model for a Mott insulator, but additionally contains a three-site term, $H_{\text{3-site}}$, and, importantly, the DH term $H_{dh}$. Each term in Eq. (2) respects the $\mathrm{SU}(2)_{\mathrm{charge}}\otimes\mathrm{SU}(2)_{\mathrm{spin}}$ symmetry of the system. See LABEL:Lee2017a for a detailed derivation for general $N$. Hereafter we discard the higher order $O(v^{3}/U^{2})$ terms.

The low-energy Hamiltonian $H_{\mathrm{eff}}$ in Eq. (2) describes two effective nearest-neighbor interactions whose role and relevance depend on the phase of the system: (i) $H_{ss}$ contains the Heisenberg spin-spin interaction. In our paramagnetic metallic phase, this interaction induces a spin-like collective mode of energy scale $\omega_{s}$. The interaction strength $v^{2}/U$ is consistent with the scaling of $\omega_{s}\sim 1/U$ (cf. Fig. 2). On the other hand, $H_{ss}$ becomes irrelevant in the paramagnetic insulating phase, where the spin susceptibility $\chi_{s}$ is overall suppressed. (ii) $H_{dh}$ describes a DH pair interaction which acts on the subspace with a finite number of DH pairs. Thus $H_{dh}$ is relevant (irrelevant) in the metallic (insulating) phase.

Doublon and holon peaks.— After integrating out the largest energy scale $U$, we consider the doublon and holon dynamics governed by the effective Hamiltonian $H_{\mathrm{eff}}$, aiming at the intermediate energy scale $\omega_{p}>\omega_{s}$, in the metallic phase. We simplify the physics at lower energies ($\lesssim\omega_{s}$) without exactly solving $H_{\mathrm{eff}}$, by introducing two approximations described in detail in LABEL:Lee2017a: (i) We introduce a mean field, $\Delta_{dh}\equiv\tfrac{v}{2}\sum_{\nu}\langle d_{i\nu}h_{j\nu}+d_{j\nu}h_{i\nu}\rangle$, which regards the Fermi-liquid ground state as the “condensate” of the DH pairs. Then we approximate the DH interaction term as $H_{dh}\approx\frac{v}{U}\sum_{\langle i,j\rangle,\nu}\Delta_{dh}^{*}(d_{i\nu}h_{j\nu}+d_{j\nu}h_{i\nu})+(\text{h.c.})$. The mean-field variable $\Delta_{dh}$, comprised of the expectation value of the pair annihilation operator $d_{i\nu}h_{j\nu}+d_{j\nu}h_{i\nu}$, is reminiscent of Bardeen-Cooper-Schrieffer theory. Here the situation is quite different, though, in that charge conservation is actually not broken, given that the pair annihilation operator is nothing but a summand of the decomposed hopping term $H_{v;-1}$. The DH pairs are singlets of the $\mathrm{SU}(2)_{\mathrm{charge}}\otimes\mathrm{SU}(2)_{\mathrm{spin}}$ symmetry preserved in the metallic phase, and the mean-field approximation of $H_{dh}$ also respects that symmetry Lee et al. (2017). (ii) We decouple the doublon and holon correlators from charge and spin density fluctations. This is based on the numerical results that they are characterized by different energy scales: charge fluctuations are suppressed in the regime $|\omega|\lesssim U/2$, and spin fluctuations predominantly occur at energies $|\omega|\lesssim\omega_{s}$ (see Fig. 1). As a result, the equations of motion for the correlators close.

Fig. 3(a) shows the resulting doublon and holon correlators for finite $\Delta_{dh}$ in the metallic phase. They have a pair of peaks at $\omega=\pm\omega_{dh}$, akin to their peaks at $\omega=\pm\omega_{p}$ in Fig. 1. Fig. 3(b) demonstrates that the DH peak position $\omega_{dh}$ from the effective theory and the DMFT+NRG result of the subpeak position $\omega_{p}$ agree well up to overall scaling factor $\simeq 4.7$ which may be expected to arise given the crudeness of our approximations. In contrast, in the insulating phase $H_{dh}$ is irrelevant, such that $\Delta_{dh}=0$. As a consequence, the subpeaks are absent in the insulating phase.

Predictions for photoemission spectroscopy.— The QP and the HBs of the local spectral functions have already been observed in photoemission spectroscopy Mo et al. (2003); Sekiyama et al. (2004). This technique, which probes the momentum-resolved spectral function $A(\omega,\mathbf{k})$ (whose momentum average yields the local $A(\omega)$ discussed hitherto), should also be able to reveal the DH subpeaks. We have thus computed $A(\omega,\mathbf{k})$, see Figs. S3 and S4 of LABEL:Supp. Our $T=0$ results agree with prior DMFT+DMRG results from LABEL:Karski2008, showing that the feature in $A(\omega,\mathbf{k})$, which leads to the subpeak in $A(\omega)$, has distinct dispersion, consistent with the interpretation of DH pair propagation. Going beyond LABEL:Karski2008, we also analyze finite $T$, and find that the subpeak-related features survive below the critical temperature for the Mott transition Sup . The distinct dispersion and $T$-dependence of the subpeak, correlated with those of the QP, distinguish it from other fine structure of the HBs originating from atomic levels. We suggest to search for such features in photoemission data, especially in multi-band materials where the subpeaks become more pronounced, as we discuss below.

$SU(N>2)$* models.—* We also analyze the $\mathrm{SU}(3)$ and $\mathrm{SU}(4)$ Hubbard models at integer filling $\langle\hat{n}_{i}\rangle=\lfloor N/2\rfloor$, with the results shown in Fig. 4. Similar to the case $N=2$ in Fig. 1, we again observe subpeaks on the inner edges of the HBs. While the subpeaks carry small weights compared with the rest of the HBs for $N=2$ [cf. Fig. 1(b)], the subpeaks for $N=3,4$ have significantly larger relative weights (cf. purple lines in Fig. 4). Even for $N=4$, the subpeaks are clearly higher than the rest of the HBs. Note that the QP persists more strongly at large $U\gtrsim 3$ for larger $N$, similarly to the widening of the Kondo peak in the $\mathrm{SU}(N)$ Kondo model Hewson (1993).

We interpret this enhancement of the subpeaks, as resulting from the enlarged space of DH pair excitations in the $\mathrm{SU}(N>2)$ Hubbard models. Generalizing the DH interaction $H_{dh}$ discussed above to the $\mathrm{SU}(N>2)$ cases, we find that the DH pair excitations on nearest neighbours are 3- and 15-fold degenerate in the $\mathrm{SU}(3)$ and $\mathrm{SU}(4)$ models, respectively, in contrast to the non-degeneracy in the $\mathrm{SU}(2)$ case Lee et al. (2017). A particularly promising area for studying this behaviour is ultracold atom physics, where pronounced DH correlations have been reported in the 2D Hubbard model Cheuk et al. (2016).

Conclusion.— We showed that the subpeaks at the inner HB edges can be related to the effective DH pair interaction by using a generalized SWT. By using NRG as a real-frequency impurity solver for DMFT, we uncovered detailed dynamical information on the decomposition of the local spectral function into doublon and holon correlators. By utilizing a recently developed broadening scheme Lee and Weichselbaum (2016), we efficiently resolved those spectral features at high energies which had been considered challenging for the NRG in the past due to its logarithmic coarse graining. An effective theory based on the scale separation of the characteristic energy scales $\omega_{s}$, $\omega_{p}$, and $U$ reproduces the linear $U$ dependence of $\omega_{p}$ found numerically in DMFT+NRG. Our predictions should be testable using photoemission spectroscopy of correlated materials, or in ultracold atom systems.

We thank M. Bukov, G. Kotliar, A. Mitchell, K. Penc, A. Polkovnikov, M. Punk, and R. Žitko for fruitful discussion. This work was supported by Nanosystems Initiative Munich. S.L. acknowledges support from the Alexander von Humboldt Foundation and the Carl Friedrich von Siemens Foundation, A.W. from the German Research Foundation (DFG) WE4819/2-1.

I Methods

We first summarize the procedure of the full-density-matrix NRG (fdm-NRG; Weichselbaum and von Delft (2007); *Weichselbaum2012:mps). Then we explain two techniques utilized in this work: adaptive broadening Lee and Weichselbaum (2016) for improving spectral resolution and $\Lambda$-averaging for suppressing discretization artifact. We also emphasize some technical aspects concerning the application of NRG to the DMFT.

I.1 Full-density-matrix NRG

We consider a single-impurity Anderson model (SIAM) in which the hybridization function between the impurity and the bath is given by $\Gamma^{\mathrm{in}}(\omega)$. The bath of the SIAM is discretized in energy space with logarithmic energy grid $\pm\omega_{\max}\Lambda^{-k-z}$, where $\Lambda>1$ is the NRG discretization parameter, $k\in\mathbb{N}_{0}$ the grid index, and $z\in]0,1]$ the discretization shift (“$z$-shift”) with integer $n_{z}$. The prefactor $\omega_{\max}$ specifies the maximum non-zero range of the input $\Gamma^{\mathrm{in}}(\omega)$; see Sec. I.4 for its determination. The value of $z$ is fixed for each separate NRG calculation. The representative energy of each discretization interval is determined by solving a differential equation Žitko and Pruschke (2009).

The discretized impurity-bath Hamiltonian is tridiagonalized to yield a semi-infinite tight-binding chain, so-called Wilson chain. The logarithmic discretization results in an overall exponential decay of the hopping amplitudes and the on-site energies in the chain. Note that the on-site energies vanish for particle-hole symmetric cases $\Gamma^{\mathrm{in}}(\omega)=\Gamma^{\mathrm{in}}(-\omega)$. In virtue of this exponential decay, energy scale separation ensures that the complete set of (well approximated) energy eigenvalues $\{E_{i}^{z}\}$ and eigenstates $\{|E_{i}^{z}\rangle\}$ of the Wilson chain can be constructed via the iterative diagonalization Anders and Schiller (2005); *Anders2006, where the superscript $z$ indicates the dependence on $z$-shift.

Using the complete basis of energy eigenstates, we can compute the discrete spectrum $A_{\mathrm{disc}}^{z}(\omega)$ of general correlation function $A_{O_{1},O_{2}}(\omega)$ in the Lehmann represesntation Weichselbaum and von Delft (2007); *Weichselbaum2012:mps,

[TABLE]

where $\rho_{i}^{z}=e^{-E_{i}^{z}/T}/\sum_{i^{\prime}}e^{-E_{i^{\prime}}^{z}/T}$ is the diagonal element of the density matrix at temperature $T$, $\pm$ in Eq. (S1b) takes the value $-(+)$ for (anti-)commuting operators $O_{1}$ and $O_{2}$, and we set $k_{\mathrm{B}}=\hbar=1$. The $\omega_{ij}^{z}$ values that determine the positions of the $\delta$ functions in $A_{\mathrm{disc}}^{z}$ are bunched, and the bunches occur roughly at $\omega\sim\pm\Lambda^{-k-z}$ ($k\in\mathbb{Z}$), reflecting the logarithmic discretization of the system Lee and Weichselbaum (2016). To reduce discretization artefacts, we average over different choices of $z$ and $\Lambda$, as described in more detail in Sec. I.2 below.

The computational efficiency and accuracy in the iterative diagonalization and in computing $A_{\mathrm{disc}}^{z}$ can be largely enhanced, especially for multi-band problems, by exploiting non-Abelian symmetries Weichselbaum (2012b). In this work we obtain the energy eigenstates as the multiplets of $\mathrm{U}(1)_{\mathrm{charge}}\otimes\mathrm{SU}(N)_{\mathrm{flavor}}$ symmetry. In the iterative diagonalization, we keep up to $N_{\mathrm{keep}}=3000$ multiplets at each step. After the first $\lceil 2\log_{\Lambda}100\rceil$ iterations, we discard the multiplets with rescaled energy above $E_{\mathrm{trunc}}=9$ for efficiency Weichselbaum (2012a).

I.2 Discretization averaging

A drawback of the logarithmic discretization is that the input hybridization function, $\Gamma^{\mathrm{in}}(\omega)$, is poorly resolved at high energies; for example, each Hubbard band in Fig. 1 contains only two or three discretization intervals. This is not a problem when $\Gamma^{\mathrm{in}}$ is featureless (as in usual impurity problems) or when the correlation functions are overbroadened at high energies (as in previous DMFT+NRG calculations), but it can induce some artificial features in correlation functions when $\Gamma^{\mathrm{in}}$ has structure and the correlation functions are finely resolved (as in our current DMFT+NRG calculations). Since not only physical structure but also artifacts (e.g., unphysical wriggle, overbroadening) can be self-reinforced during the self-consistency loop, it is necessary to suppress such discretization artifacts.

We address this problem by averaging over several discretization settings using a combination of two different schemes. The first, called $z$-averaging, is standard practice in NRG Žitko and Pruschke (2009); Bulla et al. (2008): one averages spectral functions computed for a fixed value of $\Lambda$ over $n_{z}$ equally distributed $z$-shifts $z\in\{1/n_{z},2/n_{z},\ldots,1\}$. The resulting discrete energy grids all have the same spacing on a logarithmic energy scale, but are shifted relative to each other. This $z$-averaging is also used to improve the spectral resolution of correlation functions (see Sec. I.3 below for detail). In addition, we use a second averaging scheme, which we refer to as $\Lambda$-averaging, which involves an average over different coarse grainings on the logarithmic energy scale.

We implemented these two schemes in the following manner. We first obtain multiple curves for the same correlation function by independent NRG runs using sets of different $z$ and $\Lambda$ values. After $z$-averaging the discrete data for the same $\Lambda$ followed by broadening, we then average these curves over $\Lambda$. The $\Lambda$-averaged curve of hybridization function is fed back into the DMFT self-consistency loop (see Sec. I.4 below for detail). With only two or three different $\Lambda$’s, the discretization artifacts can be significantly suppressed. In this work, we average over curves computed using the following combinations: $(\Lambda,n_{z},\alpha)=(1.7,6,1.47)$, $(2,8,1.5)$, $(2.3,10,1.56)$ for the $\mathrm{SU}(2)$ case in Fig. 1, and $(\Lambda,n_{z},\alpha)=(2.6,8,1.451)$, $(3.2,10,1.49)$, $(4,12,1.5)$ for the $\mathrm{SU}(N>2)$ case in Fig. 4. Here $\alpha$ is a parameter for broadening correlation functions (see Sec. I.3 below), and the tuples $(\Lambda,n_{z},\alpha)$ are chosen to have comparable ratios of $\tfrac{\alpha}{n_{z}}\ln\Lambda$ for each case. This $\Lambda$-averaging can be used to estimate error bars for discretization related artifacts: the shadings in Figs. 1, 4, and S1 depict the lower and upper bounds of the curves of different $\Lambda$’s, while the lines show the $\Lambda$-averaged curves.

I.3 Adaptive broadening

Physically, $A_{O_{1},O_{2}}$ should be a continuous function of $\omega$, since the original impurity model features a continuous bath before the discretization. For $T=0$, we broaden $A_{\mathrm{disc}}^{z}(\omega)$ by replacing $\delta$ functions with the log-Gaussian kernels,

[TABLE]

where $\sigma$ is a broadening width in log-frequency scale. For $T\neq 0$, the broadening kernels are modified from Eq. (S2), see LABEL:Lee2016 for details. Based on the observation on the bunching of $\delta$ functions roughly at $\omega\sim\Lambda^{-k-z}$, the conventional broadening scheme uses constant $\sigma$ for all spectral contributions, i.e., the broadening width for a weight $A_{ij}^{z}$ in linear-frequency scale is simply proportional to $|\omega_{ij}^{z}|$. Thus the spectral features at high energies (e.g., side peaks) are generally overbroadened by the conventional scheme; the bunches contributing to such features are distributed in an “irregular” manner that does not follow the $\Lambda^{-k-z}$ pattern.

A recipe to improve the spectral resolution within the conventional broadening scheme is $z$-averaging (see also above), where one broadens the discrete data averaged over different NRG calculations of $A_{\mathrm{disc}}^{z}$ for different $z$-shifts, $\bar{A}(\omega)=\tfrac{1}{n_{z}}\sum_{z}A_{\mathrm{disc}}^{z}(\omega)$. Since the $\delta$-function bunches at $\omega\sim\Lambda^{-k-z}$ for a given $z$-shift interlace with the bunches for the other $z$-shifts, one can use a narrower width, $\sigma\propto 1/n_{z}$, to achieve the continuity. Then the resolution improves as $n_{z}$ increases, but the improvement is limited Lee and Weichselbaum (2016).

Recently, two of the authors developed an adaptive broadening scheme Lee and Weichselbaum (2016) which further enhances the spectral resolution at high energies, in combination with $z$-averaging. The adaptive scheme broadens each discrete peak in $\bar{A}(\omega)$ with individual width: $\delta(\omega-\omega_{ij}^{z})$ is broadened to $\tilde{\delta}_{\sigma_{ij}}(\omega;\omega_{ij}^{z})$, where

[TABLE]

and $\alpha$ is overall prefactor of order $O(1)$. Here $\sigma_{ij}/\alpha$ estimates the distance on log-frequency scale, from one bunch of $\delta$ functions at $\omega_{ij}^{z}$ for a given $z$-shift to its neighbouring bunch at $\omega_{ij}^{z+(1/n_{z})}$ for the next $z$-shift; this estimate captures the irregular distribution of $\delta$-function bunches contributing to the spectral features at high energies, and assigns smaller broadening widths $\sigma_{ij}$ for these bunches.

As a result, the adaptive scheme better resolves such features, and the enhancement is more significant for larger $\Lambda$ or smaller $n_{z}$. Of course, the adaptive scheme retains the intrinsic accuracy of NRG at low energies; the Friedel sum rule and the Luttinger pinning are fulfilled with sub-1% error. Note that, in this work, we use the lower bound $\sigma_{ij}\geq(\ln\Lambda)/15$ to avoid unnecessary underbroadening.

I.4 DMFT

We use single-site DMFT based on the semi-circular density of states $\rho_{0}(\omega)=\tfrac{2}{\pi D^{2}}\sqrt{D^{2}-\omega^{2}}$ for the non-interacting lattice, corresponding to a Bethe lattice of coordination number $z\to\infty$ and $v\propto 1/\sqrt{z}$. We set the half-bandwidth $D\equiv 2v\sqrt{z}:=1$ as the unit of energy, as well as $\hbar=1$ throughout.

In the coexistence region of the metallic and insulating phases (for $U_{c1}<U<U_{c2}$ in our Hubbard models at $T=0$), the metallic (insulating) choice of $\Gamma^{\mathrm{in}}$ as the initial seed in the self-consistency loop results in a metallic (insulating) solution. As the metallic and insulating initial seeds, we chose the metallic solution $\pi D^{2}\rho_{0}(\omega)/4=\sqrt{D^{2}-\omega^{2}}/2$ (which is the exact solution for $U=0$), and the insulating solution for $U=3.2D$ (obtained by using the $U=0$ solution as the initial seed), respectively.

In each iteration (except for the first iteration), $\Gamma^{\mathrm{in}}(\omega)$ is determined by the result from the previous iteration that has exponentially decaying tails at large frequencies; we define $\omega_{\mathrm{max}}$ as the largest energy satisfying $\Gamma^{\mathrm{in}}(\omega_{\mathrm{max}})=\mathrm{max}(\Gamma^{\mathrm{in}})/100$, without loss of generality. Then we define the discretization grid $\pm\omega_{\max}\Lambda^{-k-z}$ as mentioned above. We compute the impurity self-energy as the ratio of two correlation functions, $\Sigma_{\mathrm{imp}}\equiv\langle[c_{\mathrm{imp},\nu},\tfrac{U}{2}n_{\mathrm{imp}}(n_{\mathrm{imp}}-1)]||c_{\mathrm{imp},\nu}^{\dagger}\rangle_{\omega}/\langle c_{\mathrm{imp},\nu}||c_{\mathrm{imp},\nu}^{\dagger}\rangle_{\omega}$ Bulla et al. (1998), where $\langle X||Y\rangle_{\omega}=G_{XY}(\omega)$ is the retarded correlation function, $c_{\mathrm{imp},\nu}$ annihilates a particle of flavor $\nu$ at the impurity, and $n_{\mathrm{imp}}=\sum_{\nu}c_{\mathrm{imp},\nu}^{\dagger}c_{\mathrm{imp},\nu}$.

As we consider the semi-elliptic density of states $\rho_{0}(\omega)$ of the lattice, the local spectral function $A(\omega)$ and the hybridization function $\Gamma^{\mathrm{out}}(\omega)=\pi D^{2}A(\omega)/4$ for the next iteration can be derived from $\Sigma_{\mathrm{imp}}$ via an explicit relation $A(\omega)=\tfrac{-2}{\pi D^{2}}\mathrm{Im}(\xi-\sqrt{\xi^{2}-D^{2}})$ without numerical integration, where $\xi\equiv\omega+\mu-\Sigma_{\mathrm{imp}}$. We continue the loop until the self-consistency criterion $|\Gamma^{\mathrm{in}}-\Gamma^{\mathrm{out}}|<\pi/10^{3}$ is satisfied.

The above relation between $\Sigma_{\mathrm{imp}}$ and $A(\omega)$, which originates from the semi-ellipticity of $\rho_{0}$, leads to two interesting properties at self-consistency, when $\Gamma^{\mathrm{in}}=\Gamma^{\mathrm{out}}$: (i) the impurity Green’s function $(\xi-\Gamma^{\mathrm{in}})^{-1}$, improved by using self-energy Bulla et al. (1998), is equivalent to the local lattice Green’s function $\tfrac{2}{D^{2}}(\xi-\sqrt{\xi^{2}-D^{2}})^{-1}$. (ii) Since $\Sigma_{\mathrm{imp}}=\omega+\mu-(\langle c_{\mathrm{imp},\nu}||c_{\mathrm{imp},\nu}^{\dagger}\rangle_{\omega})^{-1}-\tfrac{D^{2}}{4}\langle c_{\mathrm{imp},\nu}||c_{\mathrm{imp},\nu}^{\dagger}\rangle_{\omega}$, the kinks (i.e., big changes in the first derivatives) of $\Sigma_{\mathrm{imp}}$ Byczuk et al. (2007); *Raas2009 and $A_{cc^{\dagger}}\equiv\tfrac{-1}{\pi}\mathrm{Im}\langle c_{i\nu}||c_{i\nu}^{\dagger}\rangle_{\omega}$ are directly related. Therefore, the sharp peaks of $A_{dd^{\dagger}}(-\omega_{s})$ and $A_{h^{\dagger}h}(\omega_{s})$ [cf. Fig. 1] result in the kinks in $A_{cc^{\dagger}}$ as well as those in $\Sigma_{\mathrm{imp}}$.

II Supplementary results

Fig. S1 shows how the correlation functions in the metallic phase change for $2\leq U\leq 2.4$. This is similar to Fig. 1, except for the narrower range in $U$ and the zooms into $\omega\sim-\omega_{p}$ in the insets. The peak of the off-diagonal correlators, $A_{dh}(\pm\omega_{p})$ or $A_{h^{\dagger}d^{\dagger}}(\pm\omega_{p})$, grows with increasing $U$; while it is only a shoulder for $U\lesssim 2$, for $U\gtrsim 2.3$ it becomes a peak that actually exceeds the shallow peak of the diagonal correlators, $A_{dd^{\dagger}}(-\omega_{p})$ or $A_{h^{\dagger}h}(\omega_{p})$. Note that the off-diagonal correlators are negative for $|\omega|\gtrsim U/2$ for $U\gtrsim 2.1$ and their full integrals are zero by sum rule. In contrast, the diagonal correlators are positive, throughout, by construction.

Fig. S2 shows the data for the probability of single occupation $\langle P_{i1}\rangle$ which entered the mean-field decoupling analysis in Fig. 3 in the main text. It shows a clear linear dependence on $U$.

Up to now, we have shown the local (i.e., momentum-averaged) correlation functions at $T=0$. In photoemission spectroscopy experiments, the spectral functions are usually measured resolved in momentum $\mathbf{k}$ at finite $T$. Therefore we have also studied the momentum-resolved spectral functions,

[TABLE]

with the results presented in Figs. S3 and S4 for the half-filled $\mathrm{SU}(2)$ Hubbard model for different values of $U$ and $T$. In the DMFT, since the self-energy $\Sigma$ is approximated to be independent of $\mathbf{k}$, the $\mathbf{k}$-dependence of the spectral function appears only as the dependence on the non-interacting single-particle energy $\epsilon_{\mathbf{k}}$ with momentum $\mathbf{k}$, i.e., $A(\omega,\mathbf{k})=A(\omega,\epsilon_{\mathbf{k}})$. For comparison, we also plot the local spectral function $A(\omega)$, which is the average of $A(\omega,\epsilon_{\mathbf{k}})$ over momentum space,

[TABLE]

where $\rho_{0}(\epsilon)=\tfrac{2}{\pi D^{2}}\sqrt{D^{2}-\epsilon^{2}}$ is the density of states for non-interacting lattice considered for our DMFT calculations. Here we plot only the spectral functions for negative frequencies, since the photoemission spectroscopy mainly accesses the energy below the Fermi level. At half filling, one has for positive frequencies by particle-hole symmetry: $A(\omega)=A(-\omega)$ and $A(\omega,\epsilon_{\mathbf{k}})=A(-\omega,-\epsilon_{\mathbf{k}})$.

Fig. S3 illustrates how $A(\omega,\epsilon_{\mathbf{k}})$ evolves with increasing $U$. For the metallic results in panels (a)-(c), there are two most pronounced ridges of spectral intensity: First, one ridge extends from the origin at $(\omega,\epsilon_{\mathbf{k}})=(0,0)$. By parametrizing its ridge line as $\omega_{0}(\epsilon_{\mathbf{k}})$, we find an approximately linear dispersion, $\omega_{0}\simeq Z\epsilon_{\mathbf{k}}$. The slope $Z=m_{0}/m^{*}$ corresponds to the inverse effective mass of the quasiparticle [equivalently, $Z$ is the quasiparticle weight proportional to the peak position $\omega_{s}$ of the local spin susceptibility, i.e., $Z\approx 3\omega_{s}$]. As $U$ increases, the effective mass $m^{*}$ diverges at $U=U_{c2}$. Therefore the slope $\mathrm{d}\omega_{0}/\mathrm{d}\epsilon_{\mathbf{k}}$ goes to zero, as the QP disappears entirely. Second, a broad ridge is associated with the lower HB in $A(\omega)$ which stretches over a wide range of energies $\omega\lesssim-D$.

In Fig. S3(b)-(c), another intermediate ridge, whose ridge line is parametrized as $\omega_{1}(\epsilon_{\mathbf{k}})$, appears in the region $\omega_{1}\simeq-\omega_{p}$ and $\epsilon_{\mathbf{k}}\lesssim 1$, which contributes to a subpeak in $A(\omega)$ at $\omega=-\omega_{p}$. As $U$ increases within the metallic phase, this ridge becomes more pronounced, sharper, and better separated from the other ridges. On the other hand, in the insulating phase as in panel (d), both, the ridges for the subpeak and quasiparticle, disappear. Therefore the occurance or not of this intermediate ridge is completely tied to the behavior already seen in the spectral function $A(\omega$) itself.

Note that this result is consistent with the DMFT+DMRG calculation for $T=0$ by Karski et al. Karski et al. (2008). Compared to the DMFT+DMRG results of LABEL:Karski2008, the spectral resolution achieved here by DMFT+NRG is better at low energies and similar at high energies. While the former is to be expected, the latter is not, due to the reliance of NRG on logarithmic discretization. Here we nevertheless achieve a rather high resolution even at high energies by using the refined broadening scheme developed in LABEL:Lee2016.

The range of the dispersion of the intermediate ridge is comparable with the subpeak width $\delta\omega$ in $A(\omega)$. It is also consistent with the result in Fig. 2 where $\delta\omega\approx 3.5\omega_{s}$ for large $U>2.5$. The slope $\mathrm{d}\omega_{1}/\mathrm{d}\epsilon_{\mathbf{k}}$ for larger $\epsilon_{\mathbf{k}}$ on the order of the half-bandwidth $D=1$ is roughly half of the quasiparticle weight $Z\simeq\mathrm{d}\omega_{0}/\mathrm{d}\epsilon_{\mathbf{k}}$, which suggests that the underlying object responsible for the subpeaks has about twice the effective mass of an individual quasiparticle. In this sense our NRG-based numerical results are consistent with our interpretation of the subpeaks as arising from doublon-hole pairs.

The temperature dependence of the momentum-resolved spectral function $A(\omega,\epsilon_{\mathbf{k}})$ is analyzed in Fig. S4. As $T$ increases, the subpeak in $A(\omega)$ and the subpeak ridge in $A(\omega,\epsilon_{\mathbf{k}})$ become suppressed and blurred, in accordance with the suppression of the quasiparticle-related spectral features. This is consistent with our effective theory that the subpeaks originate from the doublon and holon excitations on top of the Fermi-liquid ground state, which serves as as the DH condensate. When the Fermi-liquid quasiparticles become ill-defined, also the DH condensate breaks down. Despite such thermal suppression, we emphasize that the subpeak-related features are still visible at temperatures as large as $T\lesssim 0.02$, i.e., close to the critical temperature is $T_{c}(U=2.4)\simeq 0.025$.