Magnetic and charge susceptibilities in the half-filled triangular lattice Hubbard model

We study magnetic and charge susceptibilities in the half-filled two-dimensional triangular Hubbard model within the dual fermion approximation in the metallic, Mott insulating, and crossover regions of parameter space. In the \textcolor{black}{insulating state}, we find strong spin fluctuations at the K point at low energy corresponding to the \textcolor{black}{120$^{\circ}$} antiferromagnetic order. These spin fluctuations persist into the metallic phase and move to higher energy. We also present data for simulated neutron spectroscopy and \textcolor{black}{spin-lattice} relaxation times, and perform direct comparisons to inelastic neutron spectroscopy experiments on the triangular material Ba$_8$CoNb$_6$O$_{24}$ and to the relaxation times on $\kappa$-(ET)$_2$Cu$_2$(CN)$_3$. Finally, we present charge susceptibilities in different areas of parameter space, which should correspond to momentum-resolved electron-loss spectroscopy measurements on triangular compounds.

Experimentally, much of our knowledge about correlated triangular systems is obtained from single-and two-particle scattering experiments such as photoemission [30], Raman spectroscopy [31], nuclear magnetic resonance (NMR) [2,32], or inelastic neutron scattering [9,33]. To understand these experimental results, it is necessary to calculate the corresponding response functions as a function of energy and momentum. For neutron spectroscopy and angular-resolved photoemission spectroscopy, in particular, both fine momentum and energy resolutions are desired. Such results are difficult to obtain, as computational methods formulated on finite lattices (such as ED, DMRG, and cluster DMFT) pro-vide limited momentum resolution. In addition, quantum Monte Carlo approaches are impeded by a sign problem in frustrated systems [34]. Results for these quantities are therefore often obtained from fits to quantum spin models, which are only justified in the large Coulomb interaction limit.
In this paper, we provide results for the momentum and energy dependence of the spin and charge spectra of the half-filled triangular lattice Hubbard model. We use the dual fermion (DF) approximation, which is a diagrammatic extension of the DMFT and recovers continuous momentum dependence without suffering from a sign problem. We perform simulations from weak to strong interactions and systematically study spin and charge spectra in different areas of parameter space. We then relate our results back to experiments on triangular lattice compounds and calculate nuclear magnetic resonance relaxation times.
Model. The Hubbard model is defined as · · · denotes a summation over nearest neighbors; c † iσ (c iσ ) creates (annihilates) an electron with spin σ on site i;n iσ = c † iσ c iσ is the particle number operator; U is the on-site Coulomb interaction; and t is the nearestneighbor hopping integral. We set t = 1 throughout this paper and restrict ourselves to half filling.
Method. We study the model in the ladder dual fermion approximation [35][36][37][38][39][40] using the open source code of Ref. [41]. The DF method is a diagrammatic extension of the DMFT [42] which treats all local correlations in a non-perturbative manner and perturbatively adds nonlocal charge and spin correlations [40,41]. DMFT calculations are performed with the continuous time auxiliary field quantum Monte Carlo method [43,44] with submatrix updates [45]. DF is accurate at high temperature [46] but uncontrolled in practice in the sense that adding systematic corrections, while possible in theory [47,48], is not feasible for the parameters studied here. A detailed assessment of the approximation errors of the susceptibility and the single-particle properties on the square lattice [49] showed that while dopingand interaction dependent scaling effects were present, the overall momentum dependence was accurate.
DF calculations are performed on a momentum space grid -here we choose a square 24 × 24 cluster, resulting in 288 points in the triangular lattice Brillouin zone. Both the single-particle Green's function and two-particle susceptibilities are defined on that grid. To examine the spectral properties, we use the ALPS implementation [50,51] of the maximum-entropy method [52] to perform the analytic continuation of Matsubara data to the real frequency space.
Results. The half-filled Hubbard model exhibits a metal-insulator phase transition on the triangular lattice at low temperature, which has been widely studied by DF [28,53,54] and other methods [8,20]. In this work, we set T = t/6, which is above the critical temperature of this transition, such that the system exhibits metallic behavior for U ≤ 8t, crossover behavior for 8t < U < 9.5t, and insulating behavior for U ≥ 9.5t [55]. To illustrate these three behaviors, we choose one point in each region: U = 6t in the metallic phase, U = 8.2t in the crossover region, and U = 12t in the insulating phase.
We first focus on the magnetic properties. Figure 1 shows the momentum-resolved static and dynamical magnetic susceptibility for U = 6t, U = 8.2t, and U = (U = 6t) and crossover (U = 8.2t) regions. The strong peak at U = 12t indicates the formation of 120 • antiferromagnetic (AFM) spin fluctuations [8], which will magnetically order at low temperature [8,25,56]. Figure 1(d), (e), and (f) show the imaginary part of the dynamical magnetic susceptibility Imχ m (q, ω) as a function of energy ω and momentum q along a high symmetry path in the Brillouin zone (see labels in panel (a)). White dots show the energy ω m (q) of the maximum intensity at each momentum, referred to as the spin-wave dispersion. It is clear that, except for momenta near the Γ point, the intensity of spin excitations is enhanced and the spin excitation energy decreases as U increases. At U = 6t there is no dominant spin excitation and spin fluctuations occur in a large part of the Brillouin zone. At U = 8.2t and U = 12t the spin excitation energy at the K point is smaller compared to other momenta, and most spin fluctuations occur at the K point. For U = 8.2t and U = 12t, the spin excitation energy at the Γ point is nonzero, violating the total spin conservation. This is an artifact of the DF approximation [49,57]. Figure 2(a) shows the spin excitation energy ω m (q) at the K point as a function of U . ω m (K) approaches to two different values at small and large U , and a sharp decrease occurs as U increases from 6t to 7t. Figure 2(b) shows Imχ m (q, ω) at q = K for various different U values. The maximum value of Imχ m (K, ω) increases very little as U increases from 6t to 7t, while it increases rapidly as U continues to increase. These results suggest that spin fluctuations start to condense at the K point around U = 7t. The density of states (DOS) at the Fermi surface, plotted in Fig. 2(a), shows that the system is still metallic at U = 7t. Our results therefore suggest that the strong spin fluctuations not only exist in the insulating phase but also extend into the metallic phase.
In Fig. 3 we compare our numerical data with the experimental magnetic susceptibilities obtained from Ba 8 CoNb 6 O 24 [9]. To the left of the green line, no experimental data is available. Our simulations were obtained for U = 12t and T = t/6, and we set t = 3 meV to fit the experimental data. Both our numerical data and the experimental data show that the intensity of the spin excitation around the K point is strong. We also note that the spin gap at the K point and the spin excitation energy at the M point are similar. In Ref. [9], a Heisenberg model was used to fit Fig. 3(a), giving an estimated spin-spin interaction J = 0.144 meV. In our case J = 1 meV if J = 4t 2 /U is used. Two sources contribute to this discrepancy. First, an energy scaling factor inherent to the DF approximation [49] may change the overall energy scale of the DF results, leading to a larger fit value. More importantly, the spin spectra of the Hubbard and Heisenberg models, when compared at zero temperature on a 3 × 3 lattice using ED, show agreement only for U > 20t and differ markedly at U = 12t [55]. Imχm(q,ω) ω [59].
It is shown that (T 1 T ) −1 is enhanced as U increases. At U = 12t, (T 1 T ) −1 increases rapidly as temperature decreases, indicating the formation of a magnetic order at low temperature with the transition temperature T ≈ 0.125t. For U = 8.2t, (T 1 T ) −1 increases very slowly at low temperature, consistent with the previous result that the magnetic order is absent at low temperature [53]. The increase of (T 1 T ) −1 also implies that spin fluctuations are not negligible. At U = 6t, (T 1 T ) −1 is almost independent of temperature, consistent with the weak spin fluctuations shown in Fig. 1(d).  shows the spin-lattice relaxation rate of κ-(ET) 2 Cu 2 (CN) 3 measured under different pressures, extracted from Ref. [2]. The x-axis has been rescaled by t, which is about 0.055 eV for κ-(ET) 2 Cu 2 (CN) 3 [58]. At 0 GPa, corresponding to the spin liquid region, (T 1 T ) −1 monotonously increases as temperature decreases. At 0.4 GPa, corresponding to the metallic phase near the phase boundary, (T 1 T ) −1 increases as T decreases and reaches the maximum value at about 0.03t. Continuing to decrease temperature, (T 1 T ) −1 decreases due to the appearance of a superconducting state, which is not plotted in Fig. 4(b) and absent in our calculations. At 0.8 GPa, (T 1 T ) −1 is independent of temperature. We notice that the temperature dependent behavior of (T 1 T ) −1 for these three pressures is similar to that for the three values of U we calculated. The main difference is the behavior of (T 1 T ) −1 at high temperature. In experiment, (T 1 T ) −1 has a smaller value at 0.8 GPa than that for 0 GPa and 0.4 GPa at high temperature. In our calculations, (T 1 T ) −1 approaches the same value at high temperature. This difference may be due to pressure changes of the lattice geometry in κ-(ET) 2 Cu 2 (CN) 3 [60].
We next examine the charge properties. Figure 5(a), (b), and (c) show the static charge susceptibility χ c (q, iν 0 ) for the same three values of U . It is clearly seen that χ c (q, iν 0 ) is suppressed as U increases and is invisible in the insulator (U = 12t). At U = 6t the maximum value of χ c (q, iν 0 ) is located at the Γ point, indicating a uniform charge distribution. χ c (q, iν 0 ) along the Γ → Γ direction (Γ is the Γ point in the second Brillouin zone) is larger compared to the other momenta. These features are weaker at U = 8.2 and invisible at U = 12t.  [49,57]. We also note that the maximum energy of the charge excitation does not change much as U increases before the Mott gap is opened, while it increases rapidly as the gap is opened. Our predicted charge spectra may be observed in momentum-resolved electron-loss spectroscopy measurements on triangular compounds.
Finally, we compare our magnetic and charge susceptibilities to the bare susceptibility Imχ 0 (q, ω) [55], which is evaluated by a multiplication of two Green's functions [55]. The low energy spectra of Imχ 0 (q, ω) and Imχ m (q, ω) are consistent at U = 6t but inconsistent at U = 8.2 and U = 12t. The high energy spectra of Imχ 0 (q, ω) are consistent with Imχ c (q, ω) only near the Brillouin zone boundary for these three values of U . These discrepancies suggest that the many-body effects or vertex corrections are essential.
Summary. We have studied the Hubbard model on a triangular lattice and presented the momentum and energy dependence of the spin and charge spectra in the metallic, Mott insulating, and crossover regimes. We find that the strong spin fluctuations at the K point at low energy exist in not only the insulator but also the metallic phase. We also compared our simulated data of neutron spectroscopy and relaxation times to inelastic neutron spectroscopy experiments on the triangular material Ba 8 CoNb 6 O 24 and to the relaxation times on κ-(ET) 2 Cu 2 (CN) 3 .
Our work employed the fermion Hubbard model instead of spin models which are typically used to study spin excitation spectra in frustrated systems. Unlike spin models, which are a low energy limit of the Hubbard model at large U (U > 20t), our results are valid both in the metallic and the insulating regime. Fig. S1 shows the quasiparticle weight Z, the double occupancy D = n ↑n↓ , and the local density of states as a function of interaction strength at temperature T = t/3 and T = t/6. Both of these temperatures are in a crossover regime above the metal-insulator phase transition. Z is approximately determined as Z(k F , iω 0 ) = 1/(1 − ImΣ(k F , iω 0 )/ω 0 ), where Σ(k F , iω 0 ) is the self-energy at the lowest Matsubara frequency ω 0 , and k F = (5π/6a, 0) is a momentum on the Fermi surface of the non-interacting system. The double occupancy is obtained via D = 1 2 χ loc c (τ = 0) − 2χ loc s (τ = 0) + 2 n ↑ n ↓ , where χ loc c(s) (τ ) is the local charge (spin) susceptibility. Figure S1(c) plots the local density of states (DOS) at the Fermi surface, obtained via analytic continuation of the local electron Green's function. All three quantities are large at small U , consistent with metallic behavior, and approach a small value at a large U , consistent with insulating behavior. Figure S2(a)-(c) show the momentum resolved single-particle spectral functions A(k, ω) for U = 6t, 8.2t, and 12t at T = t/6. The white dashed line indicates the Fermi surface. At U = 6t, when it is a metal, there is a strong quasiparticle peak at the Fermi surface. At U = 8.2t, the intensity of this quasiparticle peak is suppressed and the lower (upper) Hubbard band forms around K (Γ) point. At U = 12t, a Mott gap is fully opened. Our calculations show that the Mott gap is opened as U > 10t [see Fig. S1]. We notice that there is no superstructure in Fig. S2 (c) as long-ranged magnetic order is absent. All these features we observed in the spectral function are consistent with previous cluster perturbation theory results [S20].

Comparison between the Hubbard and Heisenberg models
In the strong Coulomb interaction limit (U t), the Hubbard Hamiltonian simplifies to the Heisenberg Hamiltonian with a spin-spin interaction J = 4t 2 /U . To valid this simplification, we calculate the imaginary part of spin susceptibilities Imχ m (q, ω) on a 3 × 3 lattice using exact diagonalization at zero temperature for the Hubbard and Heisenberg models, respectively. Figure S3 plot Imχ m (K, ω) as a function of frequency ω at U = 12t, 20t, and 30t. At U = 12t (J = 0.333t) Imχ m (K, ω) of the Hubbard model is located at lower energy compared to the Heisenberg model, and can be fitted by the Heisenberg model with J = 0.186t, implying that a spin-spin interaction of 4t 2 /U is overestimated by a factor of two. Imχ m (K, ω) of the Heisenberg model moves toward to Imχ m (K, ω) of the Hubbard model as U increases, and they overlap as U > 30t. These results indicate that the low energy physics of the Hubbard model can be described by the Heisenberg model only as U > 20t. The bare susceptibility is defined as where ω m (ν n ) is the fermion (bosonic) Matsubara frequency. We use the maximum-entropy method to perform the analytical continuation. Figure S4 shows the imaginary part of the bubble susceptibility χ 0 (q, ω). At U = 6t the spectrum is split into two branches. The low-energy spectrum corresponds to the spin excitation as shown in Fig. 1(b) in the main text; the high-energy spectrum corresponds to the charge excitation as shown in Fig. 4(b) in the main text. The intensity of the low-energy spectrum is suppressed at U = 8.2t. At U = 12t the low-energy spectrum completely disappears. The high-energy spectrum has a weak momentum dependence and moves to higher energy as U increases. All these results are different from magnetic and charge susceptibilities shown in the main text, implying that vertex corrections are important.