The Hubbard model on a triangular lattice -- pseudogap due to spin-density-wave fluctuations

TL;DR

This paper investigates how thermal spin-density-wave fluctuations in the Hubbard model on a triangular lattice induce pseudogap features in the fermionic spectral function, revealing non-mean-field effects and conditions for pseudogap collapse.

Contribution

It introduces a non-perturbative eikonal approach to analyze thermal fluctuation effects on spectral functions in a quasi-2D Hubbard model, highlighting conditions for pseudogap formation and collapse.

Findings

Spectral function exhibits peak/dip/hump features below T_N.

Pseudogap persists up to T_p > T_N, depending on spin susceptibility ratio.

Thermal fluctuations can be logarithmically singular without pseudogap formation.

Abstract

We calculate the fermionic spectral function $A_k (\omega)$ in the spiral spin-density-wave (SDW) state of the Hubbard model on a quasi-2D triangular lattice at small but finite temperature $T$. The spiral SDW order $\Delta (T)$ develops below $T = T_N$ and has momentum ${ \bf K} = (4\pi/3,0)$. We pay special attention to fermions with momenta ${\bf k}$, for which ${\bf k}$ and ${\bf k} + {\bf K}$ are close to Fermi surface in the absence of SDW. At the mean field level, $A_k (\omega)$ for such fermions has peaks at $\omega = \pm \Delta (T)$ at $T < T_N$ and displays a conventional Fermi liquid behavior at $T > T_N$. We show that this behavior changes qualitatively beyond mean-field due to singular self-energy contributions from thermal fluctuations in a quasi-2D system. We use a non-perturbative eikonal approach and sum up infinite series of thermal self-energy terms. We show that $A_k (\omega)$ shows peak/dip/hump features at $T < T_N$, with the peak position at $\Delta (T)$ and hump position at $\Delta (T=0)$. Above $T_N$, the hump survives up to $T = T_p > T_N$, and in between $T_N$ and $T_p$ the spectral function displays the pseudogap behavior. We show that the difference between $T_p$ and $T_N$ is controlled by the ratio of in-plane and out-of-plane static spin susceptibilities, which determines the combinatoric factors in the diagrammatic series for the self-energy. For certain values of this ratio, $T_p = T_N$, i.e., the pseudogap region collapses. In this last case, thermal fluctuations are logarithmically singular, yet they do not give rise to pseudogap behavior. Our computational method can be used to study pseudogap physics due to thermal fluctuations in other systems.

Equations

$$\mathcal{H}=$$

$$+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow }$$

$$\mathcal{H}=$$

$$+U\sum _i{\tilde{c}}_{i,+}^{†}{\tilde{c}}_{i,+}{\tilde{c}}_{i,-}^{†}{\tilde{c}}_{i,-},$$

$$T_{i,j}=\left(\begin{array}{cc}\cos (\mathbf{K}\cdot {\mathbf{R}}_{ij}/2)& \sin (\mathbf{K}\cdot {\mathbf{R}}_{ij}/2)\\ -\sin (\mathbf{K}\cdot {\mathbf{R}}_{ij}/2)& \cos (\mathbf{K}\cdot {\mathbf{R}}_{ij}/2)\end{array}\right),$$

$${\mathcal{H}}_{quad}=\sum _k\left(\begin{array}{cc}{\tilde{c}}_{k,+}^{†}& {\tilde{c}}_{k,-}^{†}\end{array}\right)\left(\begin{array}{cc}({ϵ}_{k+}+{ϵ}_{k-})/2-\mu & i({ϵ}_{k+}-{ϵ}_{k-})/2\\ -i({ϵ}_{k+}-{ϵ}_{k-})/2& ({ϵ}_{k+}+{ϵ}_{k-})/2-\mu \end{array}\right)\left(\begin{array}{c}{\tilde{c}}_{k,+}\\ {\tilde{c}}_{k,-}\end{array}\right),$$

$${\mathcal{H}}_{quad}=\sum _k\left(\begin{array}{cc}{\tilde{c}}_{k,+}^{†}& {\tilde{c}}_{k,-}^{†}\end{array}\right)\left(\begin{array}{cc}{\omega }_k-\mu & i{\eta }_k\\ -i{\eta }_k& {\omega }_k-\mu \end{array}\right)\left(\begin{array}{c}{\tilde{c}}_{k,+}\\ {\tilde{c}}_{k,-}\end{array}\right).$$

$${\mathcal{H}}_{int}$$

$$\to \Delta c_{i,-}^{†}{\tilde{c}}_{i,-}-\Delta c_{i,+}^{†}{\tilde{c}}_{i,+},$$

$${\mathcal{H}}_{MF}=\left(\begin{array}{cc}{\tilde{c}}_{k,+}^{†}& {\tilde{c}}_{k,-}^{†}\end{array}\right)\left(\begin{array}{cc}{\omega }_k-\mu -\Delta & i{\eta }_k\\ -i{\eta }_k& {\omega }_k-\mu +\Delta \end{array}\right)\left(\begin{array}{c}{\tilde{c}}_{k,+}\\ {\tilde{c}}_{k,-}\end{array}\right).$$

$$V_k=\left(\begin{array}{cc}\cos {\varphi }_k& i\sin {\varphi }_k\\ i\sin {\varphi }_k& \cos {\varphi }_k\end{array}\right),$$

$${\mathcal{H}}_{MF}=\sum _k(E_k^c{\gamma }_k^{c†}{\gamma }_k^c+E_k^v{\gamma }_k^{v†}{\gamma }_k^v),$$

$${\chi }^{\tilde{x}\tilde{x}}(q,\Omega )=\frac{{\rho }_{\Vert }}{{\Omega }^2-v_{\Vert }^2{q}^2}.$$

$${\chi }^{\tilde{y}\tilde{y}}(q\pm K,\Omega )=\frac{{\rho }_{\perp }}{{\Omega }^2-v_{\perp }^2{q}^2}.$$

$${\mathcal{H}}_{xx}$$

$${\mathcal{H}}_{yy}$$

$$i\int \mathrm{d}t{e}^{i\Omega t}⟨T{e}_{-q}^{\tilde{x}}(t)e_q^{\tilde{x}}(0)⟩$$

$$=\frac{{\rho }_{\Vert }}{{\Omega }^2-v_{\Vert }^2{q}^2}$$

$$i\int \mathrm{d}t{e}^{i\Omega t}⟨T{e}_{-q\mp K}^{\tilde{y}}(t)e_{q\pm K}^{\tilde{y}}(0)⟩$$

$$=\frac{{\rho }_{\perp }}{{\Omega }^2-v_{\perp }^2{q}^2}$$

$$\displaystyle\mathcal{H}_{m-f}=-U\sqrt{\frac{2}{N}}\sum_{k,q}\Big{(}c^{\dagger}_{k+q}\sigma^{x}c_{k}\,\hat{e}^{x}_{q}+c^{\dagger}_{k+q\pm K}\sigma^{y}c_{k}\,\hat{e}^{y}_{q\pm K}\Big{)}$$

$${\mathcal{H}}_{m-f}=-U\sqrt{\frac{2}{N}}\sum _{k,q}{\hat{e}}_q^x({\gamma }_{k+q}^{c†}{\gamma }_k^v+{\gamma }_{k+q}^{v†}{\gamma }_k^c)$$

$$-iU\sqrt{\frac{2}{N}}\sum _{k,q}{\hat{e}}_{q\pm K}^y({\gamma }_{k+q\pm K}^{c†}{\gamma }_k^v-{\gamma }_{k+q\pm K}^{v†}{\gamma }_k^c).$$

$$\displaystyle G^{c,v}(k_{hs},i\omega_{n},z)=G^{c,v\,(0)}(k_{hs},i\omega_{n})\sum_{n=0}^{\infty}\mathcal{C}_{n}(z)(\beta U^{2})^{n}\big{[}G^{v,c\,(0)}(k_{hs},i\omega_{n})G^{c,v\,(0)}(k_{hs},i\omega_{n})\big{]}^{n},$$

$$G^{c,v}(k_{hs},i{\omega }_n,z)=G^{c,v(0)}(k_{hs},i{\omega }_n)+[G^{c,v(0)}(k_{hs},i{\omega }_n){]}^2\tilde{\Sigma }(k_{hs},i{\omega }_n),$$

$${\tilde{\Sigma }}^{c,v(1)}(k,i{\omega }_n)=$$

$$-U^2\frac{T}{N}\sum _{q,m,j}{\chi }^{jj}(q,i{\Omega }_m)G^{v,c(0)}(k+q,i{\Omega }_m+i{\omega }_n)$$

$${\tilde{\Sigma }}^{c,v(1a)}(k,i{\omega }_n)$$

$$=-U^2\frac{T}{N}\sum _{q,m}\frac{{\rho }_{\Vert }}{(i{\Omega }_m{)}^2-v_{\Vert }^2{q}^2}\frac{1}{i{\Omega }_m+i{\omega }_n-E_{k+q}^{v,c}}$$

$$\approx -U^2\frac{T}{N}\sum _{q,m=0}\frac{{\rho }_{\Vert }}{-v_{\Vert }^2{q}^2}\frac{1}{i{\omega }_n-E_{k+q}^{v,c}}$$

$$=U^2\frac{{\rho }_{\Vert }T}{v_{\Vert }^2}\int \frac{{\mathrm{d}}^{2}q}{{\mathcal{A}}_{BZ}}\frac{1}{q^2}\frac{1}{i{\omega }_n-E_{k+q}^{v,c}}$$

$$\approx U^2\frac{\pi {\rho }_{\Vert }T}{v_{\Vert }^2{\mathcal{A}}_{BZ}}|\log ϵ|\frac{1}{i{\omega }_n-E_k^{v,c}}=\frac{{\beta }_1{U}^2}{i{\omega }_n-E_k^{v,c}},$$

$${\tilde{\Sigma }}^{c,v(1b)}(k,i{\omega }_n)$$

$$=\frac{{\beta }_2{U}^2}{i{\omega }_n-E_{k\pm K}^{v,c}},$$

$${\tilde{\Sigma }}^{c,v(2,crossing)}(k,i{\omega }_n)=({\beta }_1-2{\beta }_2{)}^2\times$$

$$(U^2{)}^2{G}^{v,c(0)}(k_{hs},i{\omega }_n{)}^2{G}^{c,v(0)}(k_{hs},i{\omega }_n).$$

$${\tilde{\Sigma }}^{c,v(2,non-crossing)}(k,i{\omega }_n)=2({\beta }_1+2{\beta }_2{)}^2\times$$

$$(U^2{)}^2{G}^{v,c(0)}(k_{hs},i{\omega }_n{)}^2{G}^{c,v(0)}(k_{hs},i{\omega }_n).$$

Full text

I Introduction

The pseudogap behavior, observed in several classes of materials, most notably high $T_{c}$ cuprates, remains one of the mostly debated phenomenon in correlated electron systems. There are two key scenarios of the pseudogap, each supported by a set of experiments. One is that the pseudogap is a distinct state of matter with an order parameter, which is either bilinear in fermions (e.g., loop current order Varma (1997); *Varma1999), or a four-fermion composite order (e.g., a spin nematic Kivelson et al. (1998); Fradkin et al. (2015); Fernandes et al. (2012)), or a topological order that cannot be easily expressed via fermionic operators Sachdev (2018). Within this scenario, the experimentally detected onset temperature of a pseudogap, $T_{p}$, is a phase transition temperature. The other scenario is that the pseudogap is a precursor to an ordered state – SDW magnetism Schmalian et al. (1998); *Schmalian1999; Kuchinskii and Sadovskii (1999); Abanov et al. (2003); Sedrakyan and Chubukov (2010), superconductivity Norman et al. (1998); Franz and Millis (1998); Berg and Altman (2007); Wu et al. (2019), or both, with the relative strength of the two precursors set by doping (a precursor to SDW is the dominant one at smaller dopings, and a precursor to superconductivity is the dominant one at larger dopings). Within this scenario, the system retains a dynamical memory about the underlying order in some temperature range where the order is already destroyed, and this memory gradually fades and disappears at around $T_{p}$. At around this temperature the behavior of the spectral function crosses-over to that in a (bad) metal. A similar but not equivalent scenario, is for pseudogap as a precursor to Mott physics LeBlanc et al. (2015). The precursor scenario is not strictly orthogonal to the competing order scenario as, e.g., the depletion of the spectral weight at low energies in the antinodal region due to pseudogap formation does enhance the system’s tendency to develop a CDW order with axial momenta, consistent with the one observed in the cuprates Metlitski and Sachdev (2010); Wang and Chubukov (2014); Chowdhury and Sachdev (2014); Atkinson et al. (2015). The same holds for pair density-wave order Wang et al. (2015a); *WangYuxuan2015prb; Agterberg et al. (2019). Whether a pre-existing pseudogap helps the system to develop a topological order is less clear.

In this paper we analyze several aspects of the precursor scenario. There is no clear path to get a precursor behavior at $T=0$, but earlier works Schmalian et al. (1998); *Schmalian1999; Sedrakyan and Chubukov (2010); Tchernyshyov (1997); Berg and Altman (2007); Kuchinskii and Sadovskii (1999) have found that thermal (static) SDW and/or superconducting fluctuations do give rise to precursors and associated pseudogap behavior. In particular, previous studies of quasi-2D systems on a square lattice have found that the pseudogap does develop in some $T$ range above the critical $T_{N}$ towards a commensurate $(\pi,\pi)$ SDW order Schmalian et al. (1998); *Schmalian1999; Sedrakyan and Chubukov (2010). The question we address is whether this is a generic property of a system near an ordered state, or there are situations when thermal fluctuations are logarithmically singular, but do not give rise to pseudogap behavior. To analyze this, lowest-order perturbation theory is not sufficient, and one has to sum up infinite series of singular self-energy corrections due to thermal SDW fluctuations. There is a well established computational procedure for this, similar to eikonal approximation in the scattering theory Shajesh . Here we consider, within the same computational scheme, the effects of thermal SDW fluctuations for fermions on a triangular lattice. SDW order on a triangular lattice develops with the momentum ${\bf K}=(4\pi/3,0)$ (for co-planar order in e.g. $x-z$ plane). We argue that the prefactors in the diagrammatic series for the thermal self-energy depend on the ratio of in-plane and out-of-plane susceptibilities, and by changing this ratio one can control the outcome of the summation of the series. This introduces a control parameter, by which one can vary the strength of the pseudogap behavior. We show that, for a certain value of the control parameter, the system does not develop the pseudogap, despite that self-energy corrections are singular. We note in passing that there is a similarity between this last case in our model and the “supermetal” scenario for 2D fermions near a single Van-Hove point Isobe and Fu (2019). In both cases, corrections to fermionic propagator are logarithmically singular, yet the system retains a conventional Fermi liquid behavior.

In our analysis we primarily focus on the “hot spots” in momentum space, i.e., on the ${\bf k}$ points, for which ${\bf k}$ and ${\bf k}+{\bf K}$ are both on the Fermi surface without SDW. An SDW order $\Delta(T)$ opens up a spectral gap at these ${\bf k}$. In the mean-field approximation, the spectral function $A_{k}(\omega)$ at hot spots then has two peaks at $\omega=\pm\Delta(T)$ (see Fig. 1(a)). Within mean-field, the peak position scales with the magnitude of a SDW order and vanishes right at $T_{N}$, where the order disappears. At higher $T$, the spectral function of a hot fermion is peaked at $\omega=0$, as is expected for a fermion on the Fermi surface in an ordinary paramagnetic metal (see Fig. 1(b)). Thermal fluctuations can change this behavior. For a generic value of our control parameter, $A_{k}(\omega)$ in a SDW state displays a peak, a dip, and a hump. A peak is at $\Delta(T)$, a hump is near $\Delta(T=0)$, and a dip is in between these two scales (see Fig. 2(b)). The spectral function almost vanishes below the peak, i.e., a true gap is $\Delta(T)$, like in a mean-field approximation. However, the spectral weight in the peak is reduced compared to a mean-field $A_{k}(\omega)$, and the difference is transferred into a hump. Above $T_{N}$, the peak disappears, but the hump survives up to $T=T_{p}>T_{N}$. In between $T_{N}$ and $T_{p}$, the spectral function at a hot spot is non-zero at $\omega=0$, like at ${\bf k}_{F}$ in a ordinary metal, but the maximum in $A_{k}(\omega)$ remains at a finite frequency, i.e., the system displays a pseudogap behavior (see Fig. 2(d)). As $T$ increases towards $T_{p}$, the value of $A_{k}(0)$ increases, and above $T_{p}$, the maximum of $A_{k}(\omega)$ moves to $\omega=0$ (see Fig. 2(e)).

For the special value of our control parameter, the spectral function below $T_{N}$ has only a peak at $\omega=\pm\Delta(T)$, but no hump. Above $T_{N}$, the peak disappears, and the spectral function has a single maximum at $\omega=0$, like in ordinary paramagnetic metal. In this situation, precursor behavior does not develop. Still, even in this case, the spectral function dressed by thermal SDW fluctuations is quite different from the mean-field $A_{k}(\omega)$.

For our analysis, we consider the Hubbard model on a triangular lattice, in the large $U$ limit. At $T=0$, this model displays a co-planar, $120^{0}$ SDW order with ordering momentum ${\bf K}=(4\pi/3,0)$. Within mean-field, the magnitude of the SDW order is $\Delta(0)=U|\langle\vec{S}\rangle|=U/2$ at half-filling. Quantum fluctuations reduce $|\langle\vec{S}\rangle|$ by about $50\%$ Jolicoeur and Le Guillou (1989); Capriotti et al. (1999), but do not destroy the order, nor change that $\Delta(0)\sim U$. The order gaps fermions at hot spots, and the distance between the conduction and the valence band is $2U|\langle\vec{S}\rangle|$ (which in the mean-field approximation is the Hubbard $U$).

At a finite $T$, the order is strongly affected by thermal fluctuations. In a 2D system, they destroy long-range order at any finite $T$. In a quasi-2D system, which we consider, the corrections to $|\langle\vec{S}\rangle|$ scale as $(T/J)|\log{\epsilon}|$, where $J=O(t^{2}/U)$ is the exchange interaction ($t$ is the hopping), and $\epsilon$ measures the deviations from pure two dimensionality. Long-range order get destroyed at $T_{N}\sim J/|\log{\epsilon}|$.

We will use three simplifications in our analysis. First, we neglect thermal variation of the chemical potential $\mu$. In principle, $\mu(T)$ has to be computed simultaneously with the fermionic self-energy, from the condition on the total number of fermions. In the SDW state, and in the pseudogap state above $T_{N}$, $\mu(T)$ by itself evolves with temperature, and this evolution keeps the position of the hump at distance $O(J)$ from the Fermi energy. As we said, our key goal is to analyze how the pseudogap behavior varies as we change the control parameter. Because the distance between the two humps at positive and negative $\omega$ in Fig. 2 does not depend on $\mu$, we will not include the thermal evolution of $\mu$ into our analysis and just use the non-pseudogap normal state value for $\mu$. As the consequence of keeping $\mu$ fixed, the spectral function $A_{k}(\omega)$ at a hot spot is a symmetric function of frequency, and the positions of the peak and the hump in the SDW state and the hump in the pseudogap state are all set by $U$. Second, we compute the spectral function only at $T<T_{N}$. This is enough for our purpose. Indeed, it is clear from Fig. 2 that when the spectral function retains a hump at $T=T_{N}-0$, it necessary displays a pseudogap behavior at $T>T_{N}$. And, likewise, when the spectral function does not have a hump at $T=T_{N}-0$, it does not display a pseudogap behavior at $T>T_{N}$. Third, in this work we only consider the renormalizations of the fermionic propagator due to an exchange of transverse spin fluctuations. The self-energy due to an exchange of longitudinal spin fluctuations is non-singular and we neglect it. We caution that this last approximation works well at $T$ substantially smaller than $T_{N}$, when the transverse fluctuations are gapless, but the longitudinal fluctuations are gapped. At $T\approx T_{N}$, the gap for longitudinal fluctuations gets smaller, and these fluctuations may enhance the tendency towards pseudogap behavior  Schmalian et al. (1998); *Schmalian1999.

Before we proceed, we briefly outline the computational procedure. As we said, it is similar to the eikonal approximation in the scattering theory Shajesh . To our knowledge, its was first applied in the solid state context in the study of one-dimensional (1D) systems with charge density wave (CDW) fluctuations Sadovskiǐ (1974a); *Sadovskii1974b; *Sadovskii1979; *[Foradetaileddiscussionoftheformalism; see]SadovskiiBook (see also Refs. McKenzie and Scarratt (1996)). The eikonal approximation been applied to cuprates to analyze the precursors of a collinear $(\pi,\pi)$ SDW state, in the paramagnetic phase  Schmalian et al. (1998); *Schmalian1999 and in the SDW state Sedrakyan and Chubukov (2010). It has also been used in the calculations of non-analytical corrections to Fermi liquid behavior in a 2D metal Aleiner and Efetov (2006). In our analysis, we follow Ref.  Sedrakyan and Chubukov (2010), introduce valence and conduction bands in the SDW state, and derive the vertex for the interaction between fermions and magnons (Goldstone modes of the transverse fluctuations of the order parameter). The $120^{\circ}$ SDW order is co-planar, but not collinear, which implies that it fully breaks the $SU(2)$ spin rotation symmetry. As the consequence, there are three Goldstone modes. The first two are associated with transformations that rotate the plane, where the order sets in. The third one rotates the SDW order parameter within the plane. Accordingly, there are two different spin susceptibilities, $\chi_{\bot}$ and $\chi_{\|}$ Azaria et al. (1990); Chubukov et al. (1994), for out-of-plane and in-plane rotations, respectively.

At a finite temperature, the leading contributions to the renormalization of the Green’s function come from scattering of thermal bosons (in Matsubara formalism, this corresponds to scattering processes with zero transferred bosonic frequency). These thermal self-energy contributions are logarithmically singular, and scale as powers of $T/J|\log\epsilon|$, where, we remind, $\epsilon$ measures the deviations from pure two dimensionality (see Eq. (III.1) below). We use $|\log\epsilon|\gg 1$ as a parameter, which allows us to separate singular contributions due to thermal fluctuations from regular contributions of these fluctuations (the terms of order $T/J$, without $|\log\epsilon|$), and perturbative contributions of quantum fluctuations. The latter are of order one, but are not relevant from physics perspective and can be safely neglected in our analysis of the pseudogap due to precursors to SDW Abanov et al. (2003) (non-perturbative quantum corrections may be quite relevant Sachdev (2018), but this goes beyond the scope of our work). We assume that $T_{N}/J$ is small and for $n$-loop self-energy keep only the terms of order $(T/J|\log\epsilon|)^{n}$. These terms contain both self-energy and vertex corrections, which we put on equal footings.

We show that the combinatoric factor at $n$-loop order, (which we define as $\mathcal{C}_{n}$ in the text) scales factorially with $n$, and its magnitude depends on the ratio $\chi_{\|}/\chi_{\bot}$. We sum the contributions from all orders and obtain the full self-energy. We then convert from Matsubara to real axis and obtain the spectral function $A_{k}(\omega)$.

The rest of the paper is organized as follows. In Sec. II.1 we introduce the Hamiltonian, discuss mean-field solution, and obtain the dynamical magnetic susceptibility associated with the Goldstone modes, the effective 4-fermion interaction mediated by magnons, and the magnon-fermion vertex function. In Sec. III, we obtain and sum up the series of leading logarithmical diagrams for the fermionic Green’s function. In Sec. IV we obtain the spectral function for different ratios of $\chi_{\|}/\chi_{\bot}$. This is the main result of the paper. In Sec. V we discuss the results and summarize our findings.

II The model

II.1 The Hamiltonian and the mean field solution

The point of departure for our analysis is the one band Hubbard model for spin $1/2$ fermions on a triangular lattice

[TABLE]

Without loss of generality, we restrict the hopping to nearest neighbors.

We take as an input the fact that the ground state of the model for large $U$ is the $120^{\circ}$ co-planar SDW order. Without loss the generality, we set the global coordinate such that the coplanar order is in the x-z plane, and $\langle S_{i}^{z}\rangle+i\langle S_{i}^{x}\rangle=\bar{S}\exp(\pm i{\bf K}\cdot{\bf R}_{i})$, where ${\bf K}=(4\pi/3,0)$ and $\bar{S}$ is the magnitude of average magnetization ($\bar{S}\approx 1/2$ at mean field level near half-filling and in the large $U$ limit). The $\pm$ sign in the exponent determines whether the “direction” of the spiral, i.e., whether the order is $+120^{\circ}$ or $-120^{\circ}$. Without loss of generality we consider the $+$ case in the rest of the paper. We introduce the rotating reference frame Côté and Tremblay (1995), in which all spins are along the same direction $\tilde{z}_{i}$. The transformation of fermionic operators $c_{i}$ to the new basis is given by $c_{i}=T_{i}\tilde{c}_{i}$, where $T_{i}=\exp(-i{\bf K}\cdot{\bf R}_{i}\sigma_{y}/2)$. One can straightforwardly verify that in the new coordinate frame $\langle\tilde{\bf S}_{i}\rangle=(1/2)\langle{\tilde{c}}^{\dagger}_{i,\alpha}{\bf\sigma}_{\alpha\beta}{\tilde{c}}_{i,\beta}\rangle=\{0,0,\bar{S}\}$, i.e. the original $120^{\circ}$ SDW order becomes ferromagnetic. The Hubbard Hamiltonian in the rotating reference frame takes the form

[TABLE]

where $T_{i,j}=T_{i}^{\dagger}T_{j}$. In explicit form, we have

[TABLE]

where ${\bf R}_{ij}={\bf R}_{i}-{\bf R}_{j}$.

The quadratic Hamiltonian Eq. (II.1) in momentum space is:

[TABLE]

where $\epsilon_{k+}=\epsilon_{k+K/2}$, $\epsilon_{k-}=\epsilon_{k-K/2}$ and $\epsilon_{k}=-t\sum_{\langle i,j\rangle}\exp(ik\cdot r_{ij})=-2t(\cos k_{x}+4\cos k_{x}/2\cos\sqrt{3}k_{y}/2).$ To simplify notations, we define $\omega_{k}=\frac{1}{2}(\epsilon_{k+}+\epsilon_{k-})$, $\eta_{k}=\frac{1}{2}(\epsilon_{k+}-\epsilon_{k-})$. In these notations,

[TABLE]

II.1.1 The mean-field approximation

In the SDW state, $\langle{\tilde{c}}^{\dagger}_{i,+}{\tilde{c}}_{i,+}\rangle=-\langle{\tilde{c}}^{\dagger}_{i,-}{\tilde{c}}_{i,-}\rangle=\bar{S}$. The interaction term in mean field approximation reduces to

[TABLE]

where $\Delta=U\bar{S}$. The mean field Hamiltonian becomes

[TABLE]

This Hamiltonian can be diagonalized by intoducing the operators of canonical modes $\Gamma_{k}=\{\gamma^{c}_{k},\gamma^{v}_{k}\}^{T}$ via $\Gamma_{k}=V^{\dagger}_{k}\tilde{c}_{k}$, where the unitary matrix $V_{k}$ is

[TABLE]

and \cos\phi_{k}=\sqrt{\frac{1}{2}\big{(}1-\frac{\Delta}{\sqrt{\Delta^{2}+\eta_{k}^{2}}}\big{)}}, \sin\phi_{k}=-\sqrt{\frac{1}{2}\big{(}1+\frac{\Delta}{\sqrt{\Delta^{2}+\eta_{k}^{2}}}\big{)}}. The mean-field Hamiltonian in terms of canonical modes is $\mathcal{H}_{MF}=\Gamma_{k}^{\dagger}\Lambda_{k}\Gamma_{k}$. where $\Lambda_{k}=V_{k}^{\dagger}H_{MF}V_{k}$ is diagonal. In explicit form

[TABLE]

and $E^{c}_{k}=\omega_{k}+\sqrt{\Delta^{2}+\eta_{k}^{2}}$, $E^{v}_{k}=\omega_{k}-\sqrt{\Delta^{2}+\eta_{k}^{2}}$.

A comment is in order here. The expressions for $E^{c}_{k}$ and $E^{v}_{k}$ are in the rotated coordinate frame. At $\Delta=0$, $E^{c,v}_{k}=\epsilon_{\pm}=\epsilon_{{\bf k}\pm{\bf K}/2}$. A hot spot location ${\bf k}_{hs}$ in the rotated coordinate frame is defined as a point for which $\epsilon_{{\bf k}_{hs}+{\bf K}/2}=\epsilon_{{\bf k}_{hs}-{\bf K}/2}=0$, hence $E^{c}_{k_{hs}}=E^{v}_{k_{hs}}=0$. However, once we shift ${\bf k}_{hs}$ by ${\bf K}$, we find that $E^{c,v}_{k_{hs}+{\bf K}}$ is either zero or $\epsilon_{{\bf k}_{hs}+3{\bf K}/2}$. The latter is numerically small at half-filling but strictly vanishes only for a certain hole doping. In the SDW state we then have $|E^{c,v}_{k_{hs}}|=\Delta$, but $|E^{c,v}_{k_{hs}+{\bf K}}|$ is not exactly $\Delta$, with some small correction at order $t$. Below we neglect this complication and approximate $|E^{c,v}_{k_{hs}+{\bf K}}|$ by $\Delta$.

The chemical potential $\mu$ and the SDW order parameter $\Delta$ should be obtained self-consistently as a function of the interaction strength $U$. At small $U/t$, self-consistent analysis yields $\Delta=0$, i.e., a paramagnetic Fermi liquid state with large Fermi surface remains stable down to $T=0$. This is similar to the case of the Hubbard model on a square lattice with both nearest and next nearest neighbor hopping. Once the interaction exceeds a threshold, $U>U_{c}$, the SDW order develops. This changes the Fermi surface topology to a set of small electron and hole pockets. The sizes of electron and hole pockets shrinks as $U$ increases. At half-filling, both electron and hole pockets vanish the large $U$ limit, i.e., all excitations are gapped: there is a filled valence band with energy $E^{v}_{k}=\omega_{k}-\sqrt{\Delta^{2}+\eta_{k}^{2}}$ and an empty conduction band with energy $E^{c}_{k}=\omega_{k}+\sqrt{\Delta^{2}+\eta_{k}^{2}}$. Such a state can be adiabatically connected to a Mott insulator, which, strictly speaking, does not require magnetic order. Away from the half-filling, the size of the remaining electron and hole pockets is determined by the Luttinger theorem for a SDW state Altshuler et al. (1998). We show the evolution of Fermi surface geometry with increasing $U$ in Fig. 3.

Across the transition at $U_{c}$, the spectral function at the “hot spot” $A_{c,v}(k_{hs},\omega)=-\frac{1}{\pi}\text{Im}\,G_{c,v}(k_{hs},\omega)$ changes qualitatively (see Fig. 1). At $\Delta=0$, $A_{c}(k_{hs},\omega)=A_{v}(k_{hs},\omega)=1/(\omega+i\delta)$, i.e., the spectral function is strongly peaked at $\omega=0$. At a finite $\Delta$, $G^{(0)}_{c,v}(k_{hs},\omega)=\frac{1}{\omega+i\delta\mp\Delta}$, and the maximum shifts to $\omega=\pm\Delta$. In the large $U$ limit, $\Delta=U\bar{S}$, i.e., the distance between the two peaks is Hubbard $U$, like in a Mott insulator.

II.2 Magnon-fermion interaction

To obtain the Green’s function renormalized by thermal fluctuations in the SDW state, one should first find the effective magnon-fermion interaction vertex.

The magnon propagator can be obtained by either the linear spin wave analysis Chubukov et al. (1994) or by evaluating the spin susceptibility within the generalized RPA framework in the large $U$ limit Côté and Tremblay (1995). We present the details of the spin wave analysis in Appendix A and here list the results and use physical arguments to rationalize them.

We set the $120^{\circ}$ coplanar order to be in the $x-z$ plane and move to rotating coordinate frame, where the order becomes ferromagnetic.

The the local coordinates for the $A,\,B,\,C$ sublattices are shown in Fig. 4(a). A straightforward symmetry analysis shows that in the SDW state there should be three gapless Goldstone modes, one associated with in-plane fluctuations, and two associated with out-of-plane fluctuations Azaria et al. (1990); Chubukov et al. (1994). From Fig. 4(b) we see that the in-plane transverse spin wave is along $\tilde{x}$ for all sub-lattices, which means that this Goldstone mode comes from fluctuations of $S^{\tilde{x}}$ at the $\Gamma$ point. The corresponding dynamical susceptibility is

[TABLE]

The other two Goldstone modes correspond to out-of-plane spin wave fluctuations. Figs. 4(c,d) show two independent modes. In one of them spins on the $A$ sublattice are fixed, and spins on the $B$ and $C$ sublattices rotate along $y$($\tilde{y}$) in the opposite direction. In the other mode, spins on the $B$ sublattice are fixed, and spins on the $A$ and $C$ sublattices rotate along $y$ in the opposite direction. The linear combinations of the two fluctuations yield two Goldstone modes with equal velocities around $\pm{\bf K}=(\pm 4\pi/3,0)$ in momentum space. The corresponding dynamical spin susceptibility is

[TABLE]

In the rest of the paper, we drop the $\sim$ label for local coordinates unless there is ambiguity.

In Hamiltonian approach, static $\chi^{xx}(q,0)=\chi^{xx}(q)$ and $\chi^{yy}(q,0)=\chi^{yy}(q)$ determine the effective static interaction between fermions, mediated by magnons. To get this interaction, we take the Hubbard interaction $\mathcal{H}_{int}=U\sum_{i}c^{\dagger}_{i,+}c_{i,+}c^{\dagger}_{i,-}c_{i,-}=\sum_{i}\frac{U}{2}(\hat{n}_{i,+}+\hat{n}_{i,-})-\frac{U}{2}\hat{S}_{i}^{2}$, dress it by RPA renormalizations, and keep the spin part of the dressed interaction Chubukov and Wölfle (2014) in the $\sigma^{x}\sigma^{x}$ and $\sigma^{y}\sigma^{y}$ channels. This yields

[TABLE]

We now introduce the magnon operators ($e^{x}_{q},\,e^{y}_{q\pm K}$) via

[TABLE]

A little experimentation shows that Eq. (II.2) are reproduced if we set magnon-fermion interaction to be

[TABLE]

In terms of the conduction and valence fermions, the interaction near the hot spot can be approximated as

[TABLE]

We present this interaction graphically in Fig. 5. We use a double wavy line ($\aquarius$) for magnon propagator and use a filled circle ($\bullet$) for magnon-fermion vertex with outgoing valence fermion ($\gamma^{v\dagger}$) and incoming conduction fermion ($\gamma^{c}$), and an empty circle ($\circ$) for the vertex with outgoing conduction fermion ($\gamma^{c\dagger}$) and incoming valence fermion ($\gamma^{v}$). From Eq. (II.2), the magnon-fermion vertex for $\hat{e}^{y}$ is purely maginary, and is of opposite sign for $\bullet$ and $\circ$ vertices. This turns out important when we calculate the full Green’s function at the two-loop and higher orders. The interaction terms involving only conduction or only valence fermions are small in $q$ and will not be relevant to our analysis. The presence of $q$ in these terms is consistent with the Adler principle, which states that the interaction between Goldstone bosons and fermions from the same branch should be of gradient type, to preserve the Goldstone theorem (see Ref. Chubukov and Morr (1997) for more discussions). Note that the strength of magnon-fermion interaction is of order Hubbard $U$.

III The full fermionic Green’s function in the SDW state

We now use the expressions for the quadratic SDW Hamiltonian, Eq. (7), and the magnon-fermion interaction, Eq. (II.2), and obtain the expression for the full fermionic propagator $G^{c,v}(k_{hs},\omega)$ at a finite temperature $T$. We show explicitly that the leading corrections come from the exchange of thermal transverse spin wave fluctuations. These corrections are logarithmically singular in quasi-2D systems, and n-loop correction scales as $|\log\epsilon|^{n}$, where, we remind, $\epsilon$ measures the deviation from pure two-dimensionality and serves as the infrared cutoff to regularize the divergence.

The fully renormalized Green’s function at one of the hot spots is expressed as

[TABLE]

where the combinatoric factor $\mathcal{C}_{n}(z)$ increases factorially with $n$ and depends on the ratio of in-plane and out-of-plane spin-wave susceptibilities. The factor $\beta=\frac{\pi T}{\mathcal{A}_{BZ}}|(\chi_{\|}-2\chi_{\bot})||\log\epsilon|$ measures the strength of thermal fluctuations (see below).

To simplify the presentation below, we express the fully renormalized Green’s function as

[TABLE]

where ${\tilde{\Sigma}}(k_{hs},i\omega_{n})$ can be evaluated order by order in terms of $\beta$. We will use ${\tilde{\Sigma}}^{(m)}(k_{hs},i\omega_{n})$ to label the $m^{th}$-loop correction. Note that ${\tilde{\Sigma}}^{(m)}(k_{hs},i\omega_{n})$ is not equivalent to the $m^{th}$-loop self-energy as it includes both irreducible and reducible diagrams, which we will count on equal footings. Below we still use the term “self-energy”, but in reference to $\tilde{\Sigma}$.

III.1 Perturbation theory at one-loop order

The fermion self-energy in Matsubara frequency at one-loop order is

[TABLE]

The two singular contributions to $\Sigma^{c,v\,(1)}(k,i\omega_{n})$ come from $xx$ and $yy$ components of the susceptibility. The contribution from $\chi^{xx}(q,i\Omega)$ to the leading logarithmical order is

[TABLE]

where we define $\beta_{1}=\frac{\pi\rho_{\|}T}{v_{\|}^{2}\mathcal{A}_{BZ}}|\log\epsilon|=\frac{\pi\chi_{\|}T}{\mathcal{A}_{BZ}}|\log\epsilon|$. The contributions from non-zero bosonic Matsubara frequencies are finite, and we neglect them. The contribution from $\chi^{yy}(q\pm K,i\Omega)$ is, similarly,

[TABLE]

where $\beta_{2}=\frac{\pi\rho_{\bot}T}{v_{\bot}^{2}\mathcal{A}_{BZ}}|\log\epsilon|=\frac{\pi\chi_{\bot}T}{\mathcal{A}_{BZ}}|\log\epsilon|$. At a hot spot, ${\bf k}={\bf k}_{hs}$, $E^{v,c}_{{\bf k}_{hs}}=E^{v,c}_{{\bf k}_{hs}+{\bf K}}=\Delta$. The sum of the two singular contributions then gives $\tilde{\Sigma}^{c,v\,(1)}(k,i\omega_{n})=\frac{(\beta_{1}+2\beta_{2})U^{2}}{i\omega_{n}-E^{v,c}_{k_{hs}}}$. We will see later that it is convenient to define $\beta=|\beta_{1}-2\beta_{2}|$ and re-express $\tilde{\Sigma}^{c,v\,(1)}(k,i\omega_{n})$ as \tilde{\Sigma}^{c,v\,(1)}(k,i\omega_{n})=\Big{(}\frac{\beta_{1}+2\beta_{2}}{\beta}\Big{)}\beta U^{2}G^{v,c\,(0)}(k_{hs},i\omega_{n}).

III.2 Perturbation theory at two-loop order

The self-energy at the two-loop order is obtained by summing over $27=3^{2}\,3!!$ diagrams, where $(2n-1)!!\rightarrow 3!!$ counts the number of diagrams with different topology [see Fig. 7(a)]. As there are three Goldstone modes, there are $3^{n}\rightarrow 3^{2}$ diagrams in each topology. The overall factor for each diagram is $\beta_{1}$ or $\beta_{2}$, like for one-loop diagrams, but the sign is either plus or minus. To explain the origin of sign alternation, consider the crossing diagram in Fig. 7(b) as an example.

Because the magnon-fermion vertex for $\hat{e}^{y}$ is imaginary [see Eq. (II.2)], one can show that the two-loop Green’s function from magnon propagator of $\hat{e}^{y}$ (wavy line) has prefactor $1=(\pm i)(\mp i)$ when ending with vertices of the opposite type ($\circ\aquarius\bullet$ or $\bullet\aquarius\circ$) as the magnon-fermion vertices associated with $\hat{e}^{y}$ contribute to a term $(i\gamma^{c\dagger}\gamma^{v})(-i\gamma^{v\dagger}\gamma^{c})$ in the expansion. Whereas it has prefactor $-1=(\pm i)(\pm i)$ when ending with vertices of the same type ($\bullet\aquarius\bullet$ or $\circ\aquarius\circ$) as the term takes a form $(i\gamma^{c\dagger}\gamma^{v})(i\gamma^{c\dagger}\gamma^{v})$ or $(-i\gamma^{v\dagger}\gamma^{c})(-i\gamma^{v\dagger}\gamma^{c})$. On the other hand, magnon propagator of $\hat{e}^{x}$ (dashed line) has the same prefactor $U^{2}$ for both the two ways that vertices are connected. Following these rules, we find that the two-loop self-energy from crossing diagrams is

[TABLE]

Similarly, we find that the self-energy from non-crossing diagrams [the last two diagrams in Fig. 7(a)] is

[TABLE]

The total self-energy at the two-loop order is

[TABLE]

where $\beta=|\beta_{1}-2\beta_{2}|$ and \mathcal{C}_{2}=1+2\Big{(}\frac{\beta_{1}+2\beta_{2}}{\beta}\Big{)}^{2}.

III.3 perturbation theory at $n^{th}$-loop

Following the same procedure of computing the prefactors for the diagrams in which magnon propagators are connected by vertices of the same type or of opposite types, we find the self-energy at order $n$ to be

[TABLE]

where

[TABLE]

and the coefficient $\mathcal{C}_{n}(z)$ for $n=2m$ even and $n=2m+1$ odd is expressed as (see Appendix B for details):

[TABLE]

where $l$ counts the number of magnon propagators ($2l$ for $n$ even, $2l+1$ for $n$ odd) which connect vertices of opposite type $\circ\aquarius\bullet$ or $\bullet\aquarius\circ$.

Summing up contributions from all loop orders, we find that the full fermionic Green’s function can be expressed as

[TABLE]

where $\mathcal{C}_{0}=1$ (this also follows from Eq. (III.3), if we set $m=0$).

We see from Eq. (25) that the value of $z$ depends on microscopic details, which set the ratio of the susceptibilities. In the large $U$ limit, the Hubbard model is well approximated by the nearest-neighbor Heisenberg model Chubukov and Musaelian (1994); *ChubukovMusaelian1995. For the latter, spin-wave calculations done at large $S$ (the case reproduced by taking $2S$ flavors of fermions) yield $\chi_{\|}=\frac{2}{9\sqrt{3}Ja^{2}}(1-\frac{0.449}{2S})$, $\chi_{\bot}=\frac{2}{9\sqrt{3}Ja^{2}}(1-\frac{0.285}{2S})$ Chubukov et al. (1994). Using these expressions we find $z=3-\frac{0.32}{S}+\mathcal{O}(1/S^{2})$. At smaller $U$, the value of $z$ changes, because there appear additional terms in the effective spin Hamiltonian, and, in principle, can be any number. In particular, $z=1$ when $\chi_{\|}=0$ or $\chi_{\bot}=0$; $z=\infty$ when $\chi_{\bot}=\chi_{\|}/2$.

For these two limiting cases, the combinatoric factor $\mathcal{C}_{n}(z)$ can be obtained in a closed form, as a function of $n$ (as opposed to the sum for a generic $z$, as in Eq. (III.3)). For $z=1$, there is no distinction between magnon propagators which connect vertices of opposite types, $\circ\aquarius\bullet$, $\bullet\aquarius\circ$ or of the same type, $\circ\aquarius\circ$, $\bullet\aquarius\bullet$. For $n$-loop diagrams, there are $2n$ vertices, thus there are $(2n-1)!!$ topologically distinct diagrams [see Fig. 8(a)]. In this situation, $\mathcal{C}_{n}(z=1)=(2n-1)!!$. In the opposite limit $z\to\infty$, we need to introduce $\beta^{\prime}$ via $\beta=\frac{\beta^{\prime}}{z}\rightarrow 0$, and keep $\beta^{\prime}$ as a constant. The most relevant term in Eq. (III.3) at $z\to\infty$ is the one with $l=m$. The corresponding diagrams contain only magnon propagators connecting vertices of the opposite type ( $\circ\aquarius\bullet$ and $\bullet\aquarius\circ$). At $n^{th}$-loop order, there are $n!$ topologically distinct diagrams [see Fig. 8(b)], so $\beta^{n}\mathcal{C}_{n}(z=\infty)=n!\beta^{\prime n}$, i.e., $\mathcal{C}_{n}(z=\infty)\to n!$.

We note in passing that the structure of multi-loop corrections to a fermionic propagator for $z=1$ and for $z=\infty$ is the same as quasi-1D models with CDW order/fluctuations. The case $z=1$ is realized at half-filling, when the ordering wave vector is $Q=\pi$). The case $z=\infty$ is realized in generic filling Sadovskiǐ (1974a); *Sadovskii1974b; *Sadovskii1979; *[Foradetaileddiscussionoftheformalism; see]SadovskiiBook. To our knowledge, there have been no prior analysis of a generic $z$, only specific cases have been considered. In our case, the value of $z$ is determined by ratio of $\chi_{\|}/\chi_{\bot}$ and can be arbitrary in the interval $[1,\infty)$. In the next section we analyze how the spectral function behaves for different $z$.

IV The spectral function

IV.1 Evaluation

The spectral function is defined as $A^{c,v}(k_{hs},\omega)=-\frac{1}{\pi}\operatorname{Im}G^{c,v}(k_{hs},\omega+i\delta)$, where $\omega$ is a real frequency. Our goal is to evaluate $A^{c,v}(k_{hs},\omega)$ analytically, starting from Eqs. (III.3) and (27). The key technical challenge is to perform the summation over $n$ in Eq. (27) in a situation when $\mathcal{C}_{n}(z)$ is not expressed in a closed form. We note that because $\mathcal{C}_{n}(z)\sim\mathcal{O}(n!)$, a numerical computation of $A^{c,v}(k_{hs},\omega)$ is quite challenging on its own.

Our strategy is to first sum over $l$ in Eq. (III.3), and express $\mathcal{C}_{n}(z)$ in an integral form as

[TABLE]

for $n\in\mathbb{Z}$ ($\oint^{(0+)}$ means the contour integral goes around the pole at $v=0$ counter-clockwisely, see Fig. 11(a)). Eq. (28) makes the analytic summation over $n$ in Eq. (27) possible, as $n$ only appears as an overall factor $n!$ and as an exponent in the integrand. Summing over $n$ and converting from Matsubara to real frequency ($i\omega_{m}\to\omega+i\delta$), we find the spectral function in the form of a single integral (see Appendix C for details).

[TABLE]

where $u_{\omega}=\frac{\beta U^{2}}{\omega^{2}-\Delta^{2}}$ and $\Theta(u_{\omega})$ is the Heaviside step function. Eq. (29) is the main result of this paper. In the next section we analyze qualitative features of the spectral function for different $z$.

IV.2 Results

We can extract from Eq. (29) a few generic properties of the spectral function.

• •

The presence of the Heaviside step function $\Theta(\frac{\beta U^{2}}{\omega^{2}-\Delta^{2}})$ on the r.h.s. of Eq. (29) means that $A^{c,v}(k_{hs},\omega)$ vanishes for $\omega\in(-\Delta,\Delta)$, i.e., the SDW order parameter defines the real gap at a hot spot.

• •

At the point where SDW order disappers, $\Delta=0$ and $u_{\omega}=\frac{\beta U^{2}}{\omega^{2}}$. One can understand whether at this point the system displays a pseudogap behavior or a conventional Fermi liquid behavior by analyzing how the spectral function behaves at $\omega\sim 0$. When $z>1$, the integral is bounded in the ultra-violet, and a straightforward analysis shows that $A^{c,v}(k_{hs},\omega)\sim\omega$, i.e., the maximum of $A^{c,v}(k_{hs},\omega)$ is at a finite frequency. This implies that the system displays a pseudogap behavior. At $z=1$, the upper bound of the integral becomes infinite ($\frac{1}{(z-1)u_{\omega}}\rightarrow\infty$). Then $A^{c,v}(k_{hs},\omega=0)$ doesn’t vanish. One can expand in $\omega$ and check that $A^{c,v}(k_{hs},\omega)$ has a maximum at $\omega=0$. This is the expected behavior in an ordinary Fermi liquid. We verified the dichotomy between the cases $z>1$ and $z=1$ by analytical calculations for $z=1$ and $z=\infty$ and numerical calculations for an arbitrary $z\in(1,\infty)$, as we show below.

Analytical result at $z=1$.   At $z=1$, the upper bound of the integral goes to $\infty$ and the integral in Eq. (29) can be evaluated analytically. The result is

[TABLE]

The same result was obtained in Ref. Sedrakyan and Chubukov (2010). We plot $A^{c,v}(k_{hs},\omega)$ in Fig. 9(a). Taking the limit $\Delta\to 0$, we obtain

[TABLE]

We see that at $\Delta\to 0$ (i.e., at $T\to T_{N}$) the spectral function is peaked at $\omega=0$. This implies, as we anticipated, that for $z=1$ thermal fluctuations do not give rise to SDW precursors. We emphasize that this could not be anticipated from the few first terms in loop expansion of the self-energy as these terms show little difference between $z=1$ and larger values of $z$. In the SDW phase, the spectral weight is zero at $\omega\in(-\Delta,\Delta)$, has a peak at $\omega=\pm(\Delta+0)$, and gradually decays at higher frequencies, i.e., it does not show peak/dip/hump structure. This is indeed consistent with the absence of a pseudogap at $T=T_{N}$.

Analytical result at $z=\infty$.  The integral in Eq. (29) can be evaluated analytically also at $z\to\infty$. As we discussed before, in this limit one should introduce $\beta^{\prime}$ via $\beta=\frac{\beta^{\prime}}{z}$ and keep $\beta^{\prime}$ finite. Expanding the upper and lower bounds of the integral in Eq. (29) as $\frac{1}{(z\mp 1)u_{\omega}}=\frac{\omega^{2}-\Delta^{2}}{\beta^{\prime}U^{2}}(1\mp\frac{1}{z})$, and subsituting into Eq. (29), we obtain after some algebra

[TABLE]

This result is also obtained if we replace $\beta^{n}\mathcal{C}_{n}(z=\infty)$ by $n!\beta^{\prime n}$ and directly sum up over $n$. We show $A^{c,v}(k_{hs},\omega)$ in Fig. 9 (b).

At $T=T_{N}$, we have

[TABLE]

We see that now the spectral function scales as $\sim\omega$ at small frequencies and has a maximum (a hump) at a frequency $\omega\sim\sqrt{\beta^{\prime}}U$. This implies that at $T=T_{N}$ the system retains memory about an SDW state. As $T$ increases above $T_{N}$, the maximum in the spectral function remains at a finite frequency over some range of $T$. This is a canonical pseudogap (precursor to magnetism) behavior.

Numerical results at $z\in(1,\infty)$.  We found numerically that for any $z>1$ the spectral function behaves qualitatively the same as at $z=\infty$. Fig. 9 (c) shows the spectral function for $z=3$ (the value of $z$ for our system at large $U$ and large spin S, i.e., large number of fermionic flavors). As we said above, this “universality” is expected because the upper bound of the integral in Eq. (29) is finite, as long as $z>1$. Then Eq. (29) can be simplified as

[TABLE]

At $T=T_{N}$, $A^{c,v}(k_{hs},\omega)$ scales linearly with $\omega$ at small frequencies. This necessarily implies pseudogap behavior in some range of temperatures above $T_{N}$. We found numerically that a maximum (a hump) is located at $\omega\sim\sqrt{\beta}U$.

To understand the behavior near $z=1$ we calculated numerically the position of the hump ($\omega_{hump}$) at $T=T_{N}$, as a function of $z$. We show the result in Fig. 10. Observe that $\omega_{hump}$ increases quite rapidly, as $(z-1)^{0.46}$. This indicates that the pseudogap feature at $T=T_{N}$ is quite robust, as long as $z>1$, while $z=1$ should be viewed as a special case. However, we should note that varying the value of $z$ changes the temperature where the hump starts to show up in the spectral function. In this paper, though an accurate analysis that maps the temperature $T\sim\beta$ with the order parameter $\Delta$ was not done, we can see the trend from a qualitative argument. On the one hand, $\omega_{hump}\sim(z-1)^{0.46}\sqrt{\beta}U$ increases with temperature. On the other hand, $\Delta$ decreases with temperature. As the hump shows up when $\omega_{hump}\gtrsim\Delta$, a higher temperature is needed to compensate for the smallness of $(z-1)^{0.46}$ as $z\rightarrow 1$.

IV.3 Additional considerations

First, we note that the spectral function, which we obtained, doesn’t have the coherent peak at $\omega=\pm\Delta$. This is an artifact of the approximation, in which we only include logarithmically singular contributions to the self-energy from thermal spin fluctuations. A coherent peak at $\omega=\pm\Delta$ is recovered once we add contributions to the self-energy from quantum fluctuations and non-singular self-energy piece from thermal fluctuations. This issue has been addressed in Ref.  Sedrakyan and Chubukov (2010). We show the result of including these additional terms into the self-energy (and the spectral function) by dashed line in Fig. 2 (b).

Second, in this paper we considered the self-energy due to exchange of transverse spin-wave fluctuations. As we said, such an exchange gives rise to series of logarithmically singular self-energy terms. Deep inside the SDW phase, longitudinal spin fluctuations are gapped and contribute only little to the self-energy. However, as the temperature increases towards $T_{N}$, the gap in the longitudinal fluctuations gets reduced. A more careful study at $T\lesssim T_{N}$ should take into account the contribution from longitudinal channel. Such analysis has been performed for a system on a square lattice Schmalian et al. (1998); *Schmalian1999; Sedrakyan and Chubukov (2010), and the conclusion was that longitudinal fluctuations enhance the tendency towards precursor behavior. In mathematical terms, this happens because the combinatoric factor changes from $n!$ (only transverse fluctuations, $z\to\infty$), to $(2n+1)!!$ (transverse and longitudinal fluctuations). We note in this regard that our $\mathcal{C}_{n}(z)$ does not become $(2n+1)!!$ for any $z$. As a result, the behavior of spectral function near $\omega=0$ changes from $\sim\omega$ to $\sim\omega^{2}$, and the energy of the hump increases. In our case (fermions on a triangular lattice) the analysis of the self-energy from longitudinal fluctuations is more involved than in the case of square lattice, and we refrain from making a definite prediction. Still, it is possible that longitudinal fluctuations induce some precursor behavior near $T_{N}$ even for $z=1$.

Third, it is interesting to compare our non-perturbative solution for the spectral function with a perturbative solution in which one restricts with the one-loop self-energy. In our computational approach, this implies that one includes irreducible diagrams for the Green’s function at one-loop order and only reducible diagrams at higher orders. The perturbative result for the spectral function is then

[TABLE]

where, we remind, $u_{\omega}=\frac{\beta U^{2}}{\omega^{2}-\Delta^{2}}$. This result holds for any value of $z$. We see that $A^{pert}_{k}(\omega)$ has two $\delta$-functional peaks at $\omega=\pm\sqrt{\Delta^{2}+\beta U^{2}}$. The peak frequency remains finite at $\Delta=0$, which implies that some evidence for a precursor to magnetism appears already within the perturbation theory. However, the full expression (Eqs. (C) and (C) in Appendix C), is more involved:

[TABLE]

The most essential difference is that the full $A_{k}(\omega)$ has contributions not only from the pole but also from the branch cut (the $\frac{1}{\sqrt{1-4v_{0}^{2}}}$ term). The branch cut contribution gives rise to the incoherent part of the spectral function. Combining the pole and the branch cut contributions, we obtain both the gap below $\Delta(T)$ and the hump at an energy $\Delta(0)\approx U/2$ (and we recall that with respect to the renormalized $\mu(T)$, the hump is at energy of order $J$).

V Summary

In this paper we studied the effects of thermal fluctuations on the spectral function of hot fermions on a triangular lattice, in the $120^{\circ}$ SDW state (the ordering momentum is ${\bf K}=(4\pi/3,0)$).

We argued that the exchange of static Goldstone bosons between fermions in the valence and the conduction band gives rise to logarithmically singular self-energy corrections. We obtained fully renormalized Green’s function by summing up infinite series of thermal self-energy diagrams. In this sense, we went beyond perturbation theory.

The key goal of our study was to understand whether the exchange of static thermal bosonic fluctuations necessarily gives rise to pseudogap behavior, or the system may display a conventional Fermi liquid behavior despite that self-energy corrections are logarithmically singular. We argued that one can address this issue by studying fermions on a triangular lattice. Specifically, we showed that the contributions from in-plane and out-of-plane spin-wave fluctuations are not equivalent, and the strength of the self-energy renormalizations depends on the ratio of in-plane and out-of-plane spin-wave susceptibilities $\chi_{\|}/\chi_{\bot}$. This ratio is an input parameter for low-energy theory, and by varying it one can study the changes in the structure of diagrammatic series for the self-energy. When $\chi_{\|}/\chi_{\bot}\sim 1$, our calculations show that the behavior of the spectral function for a fermion at a hot spot (${\bf k}_{hs}$ and ${\bf k}_{hs}+{\bf K}$ are both on the Fermi surface) is similar to that for the case of a square lattice and collinear $(\pi,\pi)$ SDW order: there is a real gap below $\Delta(T)$ and a maximum (hump) at an energy of order Hubbard $U$ (see Fig. 2). The hump persists at $T=T_{N}$, where $\Delta$ vanishes, and survives in some range of $T$ above $T_{N}$. In this range the system displays a pseudogap behavior. On the other hand, when $\chi_{\|}/\chi_{\bot}\ll 1$ or $\chi_{\|}/\chi_{\bot}\gg 1$, the pseudogap behavior exists only near $T_{N}$. In the limiting case when $\chi_{\|}/\chi_{\bot}=0$ or $\chi_{\|}/\chi_{\bot}=\infty$, there is no pseudogap behavior at any $T>T_{N}$, despite that perturbative self-energy corrections are logarithmically singular.

The calculations, which we presented in this paper, can be readily extended to other microscopic models, and our results show that by looking at the structure of perturbation series one would be able to immediately conclude whether singular self-energy corrections lead to a pseudogap behavior, or to an ordinary Fermi liquid behavior.

VI Acknowledgements

We thank L. Fu, A. Georges, E. Gul, E. Kozik, J.P.F. LeBlank, J. Schmalian, F. Simkovic and Z. Zhang for helpful discussions. The work was supported by the NSF DMR-1523036. M.Y. also acknowledges support from Louise Dosdall Fellowship from the University of Minnesota.

Appendix A Goldstone modes

In this appendix, we give more mathematical details of how the momentum and spin structure of the low energy bosonic collective modes in the SDW state, i.e., the magnon Goldstone modes, can be obtained by studying the eigenstates of its corresponding linear spin wave Hamiltonian. As we show in the main text, these universal properties of the Goldstone modes are enough to obtain the effective magnon-fermion interaction.

To be precise, as our starting point is the SDW state out of itinerant fermions, the collective modes in the interacting fermion system, magnons in this case, should be obtained in principle by calculating the spin correlation function. On the other hand, as the low energy magnons are essentially the Goldstone modes which are uniquely determined by the pattern of the spontaneous spin rotation symmetry breaking, the universal properties of the Goldstone modes can be obtained from other models that are adiabatically connected to the Hubbard model in the large $U$ limit, e.g. the isotropic Heisenberg model. For the square lattice Hubbard model, it has been checked that the spectrum calculated from the two methods match over the whole Brillouin Schrieffer et al. (1989); Chubukov and Frenkel (1992). For the triangular lattice Hubbard model, though it turns more tricky at higher energy Côté and Tremblay (1995), the matching should work in principle for the low energy Goldstone modes for the reason we explained above. We also checked that it is indeed the case. In the following, we work with the nearest neighbor isotropic Heisenberg antiferromagnetic model using linear spin wave analysis. The large $S$ spin wave expansion can be reproduced by taking $2S$ fermion flavors in the Hubbard model.

The global and local coordinates are set up as shown in Fig. 4(a). To obtain the linear spin wave Hamiltonian, we express the spin operators in terms of Holstein-Primakoff bosons $a,\,a^{\dagger}$. In the local coordinate where the magnetic order is along $\tilde{z}$, the spin operator is:

[TABLE]

The spin operators in the global coordinate can be expressed as:

[TABLE]

where $i$ is the sublattice index $i=a,b,c$, and $\theta_{i}=0,2\pi/3,4\pi/3$. The linear spin wave Hamiltonian can be expressed as

[TABLE]

where $\Psi_{k}=\{a_{k},a^{\dagger}_{-k}\}^{T}$,

[TABLE]

$A_{k}=J(3+\frac{1}{2}\gamma_{k})$, $B_{k}=-\frac{3}{2}J\gamma_{k}$, $\gamma_{k}=\cos k_{x}+2\cos k_{x}/2\cos\sqrt{3}k_{y}/2$. The spin wave Hamiltonian can be diagonalized by the transformation $\Psi_{k}=T\Psi^{\prime}_{k}$, where $\Psi^{\prime}_{k}=\{e_{k},e^{\dagger}_{-k}\}^{T}$, such that $T^{-1}\sigma_{z}H_{k}T=\omega_{k}\sigma_{z}$. We found $\omega_{k}=\sqrt{A_{k}^{2}-B_{k}^{2}}$, and

[TABLE]

As the $SU(2)$ spin rotation symmetry of the Hamiltonian is fully broken by the magnetic order, there are three Goldstone modes, associated with the three broken symmetry generators. It is straight forward to show that there are three zero modes at momentum $\Gamma=(0,0)$ and $\pm K=\pm(\frac{4\pi}{3},0)$, respectively. To check if they are the Goldstone modes and learn the spin structure, let us analyze the eigenmodes near $\Gamma,\pm K$.

Near $\Gamma$ – The spin wave spectrum is $\omega_{q+\Gamma}=\frac{3\sqrt{3}JS}{4}q$. The transformation matrix $T_{k}$ is

[TABLE]

The part singular in $\sqrt{q}$ at order $1/\sqrt{q}$ corresponds to the Goldstone mode excitation, which contributes to the divergent static susceptibility as $q\rightarrow 0$, while the subleading term at order $\sqrt{q}$ corresponds to the soft modes, whose static susceptibility is finite at $q=0$. From the singular part, we obtain the leading order dynamical spin susceptibility (labeled by superscript “(0)”) at $S=1/2$

[TABLE]

Similarly, we obtain from the $\sqrt{q}$ order terms in Eq. (42) the next order dynamical spin susceptibility (labeled by superscript “(1)”) at $S=1/2$

[TABLE]

Note that by $\tilde{x},\tilde{y}$, we mean the spin component in the local coordinates. For our interest of obtaining the logarithmical divergent contribution to the thermal self-energy, only non-zero terms in Eq. (43) is needed, which physically means spin fluctuations along the local $\tilde{x}$ direction at the $\Gamma$ point (see Fig. 4(b)).

Near $\pm K$ – By doing the same analysis near $\pm K$, we found the spin wave spectrum is $\omega_{q\pm K}=\frac{3\sqrt{3}JS}{2\sqrt{2}}q$, and to the leading order in $q$,

[TABLE]

The leading order dynamical spin susceptibility at $S=1/2$ is

[TABLE]

Eqs. (43) and (46) are essentially what we have in Eq. (II.2) in the main text. The two Goldstone modes are shown graphically in Fig. 4(c),(d).

Appendix B Calculation of $\mathcal{C}_{n}(z)$

To obtain $\mathcal{C}_{n}(z)$ for a generic $n$, we review the rules found in Sec. III.1. These are

• •

We use $\circ$ and $\bullet$ for magnon-fermion vertex $\hat{e}\,\gamma^{v\dagger}\gamma^{c}$ and $\hat{e}\,\gamma^{c\dagger}\gamma^{v}$, respectively. As vertices like $\hat{e}\,\gamma^{v\dagger}\gamma^{v}$ are not considered to the leading logarithmical order, the renormalized Green’s function at $n$-loop order should have $n$ pairs of alternating $\circ$ and $\bullet$ vertices;

• •

Adding contributions from $\chi^{xx}$ and $\chi^{yy}$, each $\circ\aquarius\circ$ or $\bullet\aquarius\bullet$ magnon propagator contribute a factor $(\beta_{1}-2\beta_{2})U^{2}$;

• •

Similarly, each $\circ\aquarius\bullet$ or $\bullet\aquarius\circ$ magnon propagator contribute a factor $(\beta_{1}+2\beta_{2})U^{2}$.

All diagrammatic configurations at $n$-loop order can be grouped by the total number of $\bullet\aquarius\bullet$, $\circ\aquarius\circ$, $\circ\aquarius\bullet$ and $\bullet\aquarius\circ$ propagators, and the contribution from each diagram is the same within a given group. In the following, we discuss $n=2m$ even and $n=2m+1$ odd seperately.

B.1 n=2m

Each diagram is in a group (labeled by $l$) that has $m-l$ of $\bullet\aquarius\bullet$, $m-l$ of $\circ\aquarius\circ$ propagators, and a total $2l$ of $\bullet\aquarius\circ$ or $\circ\aquarius\bullet$ propagators, where $l=0,1,...,m$. For a given group labeled as $l$, the combinatoric factor contributing to the renormalized Green’s function is

[TABLE]

where the first $\{...\}$ in the first line above comes from contributions of $\bullet\aquarius\bullet$ and $\circ\aquarius\circ$ propagators, and the second $\{...\}$ comes from contributions of $\bullet\aquarius\circ$ and $\circ\aquarius\bullet$ propagators. Summing up all factors, we find \mathcal{C}_{2m}(z)=\sum_{l=0}^{m}\big{[}\frac{(2m)\,!}{(2m-2l)\,!!}\big{]}^{2}\frac{z^{2l}}{(2l)\,!}.

B.2 n=2m+1

For $n$ odd, as the total number of $\circ$ or $\bullet$ vertices used up for $\circ\aquarius\circ$ or $\bullet\aquarius\bullet$ propagators must be even, there must be a total odd number of $\bullet\aquarius\circ$ and $\circ\aquarius\bullet$ propagators. As a result, each group labeled by $l$ has $m-l$ of $\bullet\aquarius\bullet$, $m-l$ of $\circ\aquarius\circ$ propagators, and a total $2l+1$ of $\bullet\aquarius\circ$ and $\circ\aquarius\bullet$ propagators, where $l=0,1,...,m$. The combinatoric factor is

[TABLE]

Summing up all factors, we find \mathcal{C}_{2m+1}(z)=\sum_{l=0}^{m}\big{[}\frac{(2m+1)\,!}{(2m-2l)\,!!}\big{]}^{2}\frac{z^{2l+1}}{(2l+1)\,!}.

Appendix C Evaluate the spectral function

We now evaluate the spectral function defined as $A^{c,v}(k_{hs},\omega)=-\frac{1}{\pi}\operatorname{Im}G^{c,v}(k_{hs},\omega+i\delta)$ analytically starting from Eqs. (III.3) and (27). The key challenge is to perform the summation over $n$ in Eq. (27) where $\mathcal{C}_{n}(z)$ doesn’t have a simple closed form. Moreover, as $\mathcal{C}_{n}(z)\sim\mathcal{O}(n!)$, a numerical calculation of $A^{c,v}(k_{hs},\omega)$ is quite challenging on its own. Our point of departure is to sum over $l$ in Eq. (III.3) by noting that from $\frac{(p+q)\,!}{p\,!q\,!}=\frac{1}{2\pi i}\oint^{(0+)}\mathop{}\!\mathrm{d}t\,t^{-p-1}(1-t)^{-q-1}$ [$p,q\in\text{Integers}\,(\mathbb{Z})$ and $p+q\geq-2$] Xi and Ren (2000),

[TABLE]

where $\oint^{(0+)}$ means the contour integral goes around the pole at $t=0$ counter-clockwisely [see Fig. 11(a)]. For concreteness, we take $n=2m$ as an example. Plug Eq. (49) into Eq. (III.3), we have

[TABLE]

To sum over $l$ in Eq. (50), we note that from $(1+x)^{m}=\sum_{p=0}^{m}\frac{m!}{(m-p)!\,p!}x^{p}$,

[TABLE]

where we define $v=\sqrt{t(1-t)}$. Note that by summing over $l$, non-analytic branch-cuts must be introduced, which turns important to get the imaginary part of the spectral function later. Changing variable from $t$ to $v$, the integration contour changes from a circle around $t=0$ to a semi-circle around $v=0$ on the right-half-plane, and $\mathop{}\!\mathrm{d}t=\frac{2v}{\sqrt{1-4v^{2}}}\mathop{}\!\mathrm{d}v$. By adding the two terms in Eq. (C) and changing variable $v\rightarrow-v$ for the second term, Eq. (50) becomes [see Fig. 11(b)]

[TABLE]

Similarly, we found that $\mathcal{C}_{2m+1}(z)$ has the same form but changes $2m$ in the expression to $2m+1$, thus

[TABLE]

From Eq. (50) to Eq. (53), we are essentially transforming the summation over $l$ in Eq. (50) into evaluating the residue of the integrand at $v=0$ in Eq. (53). The gain we have is that the number $n$ now only appears as a simple coefficient $n!$ and as exponents in the integrand, which simplifies the summation over $n$ in Eq. (27).

To find $A^{c,v}(k_{hs},\omega)$, we analytically continue $G^{c,v}(k_{hs},\omega)$ from Eq. (27) to real frequencies by replacing $i\omega\rightarrow\omega+i\delta$. As long as the series over $n$ converge, doing the analytical continuation before or after the summation over $n$ should give the same result. For convenience, we perform the analytical continuation before the summation, and have

[TABLE]

where $u(\omega+i\delta)=\beta U^{2}G^{v,c\,(0)}(k_{hs},\omega+i\delta)G^{c,v\,(0)}(k_{hs},\omega+i\delta)$. The imaginary part inside $\{...\}$ comes from two places – from $\operatorname{Im}G^{c,v\,(0)}(k_{hs},\omega+i\delta)\times\operatorname{Re}\sum_{n=0}^{\infty}...$ and from $\operatorname{Re}G^{c,v\,(0)}(k_{hs},\omega+i\delta)\times\operatorname{Im}\sum_{n=0}^{\infty}...$. As $\operatorname{Im}G^{c,v\,(0)}(k_{hs},\omega+i\delta)=-i\pi\delta(\omega\mp\Delta)$, the first contribution should be a delta-function peak at $\omega=\pm\Delta$ if $\int_{\Delta-0}^{\Delta+0}\mathop{}\!\mathrm{d}\omega A(\omega)$ is finite. We checked that this integral actually vanishes. This implies that thermal fluctuations destroy the delta-function peak. To evaluate the second contribution, we use Eq. (53) and $n!=\int_{0}^{+\infty}\mathop{}\!\mathrm{d}t\,e^{-t}t^{n}$, express $\sum_{n=0}^{\infty}\mathcal{C}_{n}(z)[u(\omega+i\delta)]^{n}$ in Eq. (C) as

[TABLE]

where $u(\omega+i\delta)$ is replaced by $u_{\omega}$ for brevity. Importantly, by summing over $n$, the multi-pole at $v=0$ vanishes, while a single pole at $v=v_{0}=\frac{u_{\omega}\,t}{2(1-u_{\omega}\,t\,z)}$ emerges due to the non-analyticity at $v=0$. Then the contour $\oint^{(0+)}$ changes to $\oint^{(0-v_{0},+)}$, where $\oint^{(0-v_{0},+)}$ indicates a counter-clockwise contour enclosing $v=0$ and $v=v_{0}$ [see Fig. 11(c)]. Enforcing $\omega\rightarrow\omega+i\delta$, we find $u_{\omega}\rightarrow u_{\omega}-i\delta\text{sgn}(\omega)$ and $v_{0}\rightarrow v_{0}-i\delta\text{sgn}(\omega)$. As $v_{0}$ remains on the lower or upper half-plane [depending on $\text{sgn}{(\omega)}$] as $t$ varies, the integrals from [math] to $v_{0}$ and from $v_{0}$ to [math] cancel, and Eq. (C) becomes

[TABLE]

To obtain the imaginary part of Eq. (56), one can show that the only contribution comes from the residue of the integrand of $\oint^{(v_{0}+)}$ when $v_{0}$ sits at the branch cut, i.e., $|v_{0}|\geq 1/2$. By examining $v_{0}=\frac{u_{\omega}\,t}{2(1-u_{\omega}\,t\,z)}$ for $t\in(0,\infty)$, we find $|v_{0}|\geq 1/2$ only when $u_{\omega}>0$ and $t\in(\frac{1}{(z+1)u_{\omega}},\frac{1}{(z-1)u_{\omega}})$ [see Fig. 11(d)]. In particular, if $z=1$, the upper bound for $t$ is $+\infty$. Thus the imaginary part of Eq. (56) when $u_{\omega}>0$ is

[TABLE]

In the following, let us consider $\omega>0$ for concreteness, so $v_{0}\rightarrow v_{0}-i\delta$. Note that $\frac{1}{\sqrt{1-4v_{0}^{2}}}$ is pure imaginary when $t\in(\frac{1}{(z+1)u_{\omega}},\frac{1}{(z-1)u_{\omega}})$. As we explain below, it needs some care to determine the sign at different $t$. From Fig. 11(d), we see that when $t\in(\frac{1}{(z+1)u_{\omega}},\frac{1}{zu_{\omega}})$, $\operatorname{Re}v_{0}>1/2$, so $\operatorname{Im}\frac{1}{\sqrt{1-4v_{0}^{2}}}=\frac{1}{\sqrt{1+2v_{0}}\sqrt{1-2v_{0}}}=\frac{-i}{\sqrt{1+2v_{0}}\sqrt{2v_{0}-1}}=\frac{-i}{\sqrt{4v_{0}^{2}-1}}$, $1-u_{\omega}\,t\,z>0$; when $t\in(\frac{1}{z\,u_{\omega}},\frac{1}{(z-1)u_{\omega}})$, $\operatorname{Re}v_{0}<-1/2$, so $\operatorname{Im}\frac{1}{\sqrt{1-4v_{0}^{2}}}=\frac{1}{\sqrt{1+2v_{0}}\sqrt{1-2v_{0}}}=\frac{i}{\sqrt{4v_{0}^{2}-1}}$, $1-u_{\omega}\,t\,z<0$. So Eq. (C) becomes

[TABLE]

where the integral over $t$ is convergent and positive definite. Similarly, we find when $\omega<0$, $i\operatorname{Im}\sum_{n=0}^{\infty}\mathcal{C}_{n}(z)u_{\omega}^{n}=i\int_{\frac{1}{(z+1)u_{\omega}}}^{\frac{1}{(z-1)\,u_{\omega}}}\mathop{}\!\mathrm{d}t\,e^{-t}\,\frac{1}{\sqrt{(u_{\omega}\,t)^{2}-(1-u_{\omega}\,t\,z)^{2}}}$. Plug them back to Eq. (C), the spectral function is

[TABLE]

where we remind $u_{\omega}=\frac{\beta U^{2}}{\omega^{2}-\Delta^{2}}$.