Introduction

Recent discovery of the kagome-lattice metal AV3Sb5 (A = K, Rb, Cs) shown in Fig. 1a has opened the way to study the unique physics of geometrically frustrated metals with strong correlation1,2,3. In CsV3Sb5, the formation of the 2 × 2 Star-of-David or Tri-Hexagonal density wave (DW) was detected by scanning tunneling microscopy (STM) at T ≈ 90 K at ambient pressure4,5. It is presumably the triple-q (3Q) bond order (BO) shown in Fig. 1b, which is the even-parity modulation in the hopping integral \(\delta {t}_{ij}^{{{{{{\rm{b}}}}}}}\) (=real)6,7,8,9,10,11,12. Below the BO transition temperature TBO, superconductivity (SC) with highly anisotropic gap emerges for A = Cs13,14, and the gap structure changes to isotropic by introducing impurities. Also, nodal to nodeless crossover is induced by the external pressure in A = Rb,K15. These results are naturally understood based on the BO fluctuation mechanism12.

Fig. 1: Lattice structure, Fermi surface, and BO form factor in kagome metal.
figure 1

a Kagome-lattice structure composed of the sublattices A, B, and C. 2alm is the minimum translation vector, and we set ∣2alm∣ = 1. The relation aAB + aBC + aCA = 0 holds. b 3Q Tri-Hexagonal bond order (BO) state. c Fermi surface (FS) at n = 0.917 with the nesting vectors q1, q2, q3. The color of the FS represents the weight of the sublattice (A = red, B = blue, C = green). The FS has single sublattice character near the van-Hove singularity (vHS) points. In kagome metals, q1 connects vHS-A and vHS-B. It is given as \({{{{{{\bf{q}}}}}}}_{1}=(2{{{{{{\bf{a}}}}}}}_{{{{{{\rm{AB}}}}}}})\times (2\pi /\sqrt{3}){{{{{{\bf{e}}}}}}}_{z}\), where ez is the unit vector perpendicular to the xy-plane.

Full size image

More recently, the non-trivial time reversal symmetry breaking (TRSB) order at TTRSB attracts considerable attention. It has been reported by μSR study15,16,17,18, Kerr rotation analysis19, field-tuned chiral transport study20, and STM measurements4,20. The transition temperature TTRSB is close to TBO in many experiments, while the TRSB order parameter is strongly magnified at T* ≈ 35 K for A = Cs16,18,20 and T* ≈ 50 K for A = Rb15. Recently, magnetic torque measurement reveals the TRSB order associated with the rotational symmetry breaking, which is called the nematic order, at T* ≈ 130 K21. In contrast, TRSB was not reported by different experimental groups using the Kerr rotation22 and STM23 measurements. Thus, the TRSB onset temperature is still under debate. The chiral cLC is driven by the additional odd-parity hopping integral \(\delta {t}_{ij}^{{{{{{\rm{c}}}}}}}\) (=imaginary), and the accompanied topological charge-current24 gives the giant anomalous Hall effect (AHE) below T ≈ 35 K25,26. The correlation-driven topological phase in kagome metals is very unique, while its mechanism is still unknown.

In addition to the cascade of quantum phase transitions, the emergent nematic order inside the BO and the cLC phases attracts great attention. The nematic transition is clearly observed by the elastoresistance27, the scanning birefringence19, and the STM5 studies. In addition, nematic SC states have been reported23,28. Thus, kagome metals provide a promising platform for exploring the interplay between electron correlations and topological nature.

To understand the rich quantum phases in kagome metals, lots of theoretical studies have been performed6,7,8,9,10,11,12,29,30,31. Each BO and cLC order is explained by introducing various off-site interactions in the mean-field approximation (MFA)9,29,32,33, while a fine-tuning of off-site interactions is necessary to explain the cascade of phase transitions. On the other hand, beyond-MFA mechanisms have been developed to explain the rich phase transitions34,35,36,37,38,39,40,41,42,43,44,45. For example, strong interplay between the magnetism, nematicity and SC in Fe-based superconductors and other strongly correlated metals were understood by beyond-MFA mechanisms36,37,38,39,40,41,42,43,44,45. It is urgent and important to elucidate why the BO and cLC orders/fluctuations coexist in the study of kagome metals. For example, these fluctuations will mediate non-BCS SC12 and exotic pair-density-wave states46,47,48.

In this paper, we reveal that the cLC order is mediated by the BO fluctuations that are abundant above TBO in kagome metals49,50. The sizable off-site Umklapp scattering by the BO fluctuations induces the odd-parity and TRSB current order (=imaginary \(\delta {t}_{ij}^{{{{{{\rm{c}}}}}}}\)). This cLC mechanism is universal because it is irrelevant to the origin of the BO. Furthermore, we discover that the coexistence of the BO and the cLC order gives rise to the novel Z3 nematicity along the three lattice directions reported in refs. 5,19,27. The present theory reveals the close relationship between the cLC, BO, nematicity, and SC state, which is significant to understand the unsolved quantum phase transitions in kagome metals.

The phase transitions in metals are described as the symmetry breaking of the normal self-energy; \(\Delta \Sigma \equiv \Sigma -{\Sigma }_{{A}_{1g}}\)45,51. ΔΣ is determined by the stationary condition of the free energy; δF[ΔΣ]/δ(ΔΣ) = 0. The DW equation enables us to derive the solution that satisfies the stationary condition, as we proved based on the Luttinger-Ward theory51. Based on the DW equation, we discover that the odd-parity and TRSB ΔΣ is driven by the BO fluctuation exchange processes. (Note that the DW equation for ΔΣ is analogous to the Eliashberg equation for the SC gap Δ).

Results

BO form factor and fluctuations

Here, we introduce the kagome-lattice tight-binding model with a single d-orbital of each vanadium site (A, B, or C) shown in Fig. 1a. (The d-orbital belongs to b3g of the D2h point group at V site, while its representation is not essential here). The kinetic term is given by \({\hat{H}}_{0}={\sum }_{{{{{{\bf{k}}}}}},l,m,\sigma }{h}_{lm}^{0}({{{{{\bf{k}}}}}}){c}_{{{{{{\bf{k}}}}}},l,\sigma }^{{{\dagger}} }{c}_{{{{{{\bf{k}}}}}},m,\sigma }\), where l, m denote the sublattices A, B, C, and \({h}_{lm}^{0}({{{{{\bf{k}}}}}})\,(={h}_{ml}^{0}{({{{{{\bf{k}}}}}})}^{*})\) is the Fourier transform of the nearest-neighbor hopping integral t in ref. 52 in addition to the inter-sublattice hopping \({t}^{{\prime} }\) shown in Fig. 1a. We set t (=−0.5eV) to fit the bandwidth, and \({t}^{{\prime} }(=-0.08\,{{{{{\rm{eV}}}}}})\) to reproduce the shape of the Fermi surface (FS). Numerical results are insensitive to the presence of \({t}^{{\prime} }\). Hereafter, the unit of energy is eV unless otherwise noted. The FS around the van-Hove singularity (vHS) point (k ≈ kA, kB, or kC) is composed of a single 3d-orbital on V ion, which is called the sublattice interference6. This simple three-site model well captures the main pure-type FS in kagome metals4,53,54,55,56,57. The FS at the vHS filling (nvHS = 0.917 per site and both spins) is shown in Fig. 1c. The wavevectors of the BO correspond to the inter-sublattice nesting vectors qn (n = 1, 2, 3) in Fig. 1c. (The equivalent square lattice kagome model is convenient for the numerical study; see Supplementary Note 1). The good inter-sublattice nesting of the FS naturally triggers the observed inter-sublattice BO at q = qn, as shown in previous theoretical studies6,7,12.

The Fourier transform of the BO modulation, \(\delta {t}_{ij}^{{{{{{\rm{b}}}}}}}\), gives the even-parity BO form factor \({g}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}})\)45,58:

$${g}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}})=\frac{1}{N}\mathop{\sum }\limits_{i}^{{{{{{\rm{sub}}}}}}-l}\mathop{\sum }\limits_{j}^{{{{{{\rm{sub}}}}}}-m}\delta {t}_{ij}^{{{{{{\rm{b}}}}}}}{e}^{i{{{{{\bf{k}}}}}}\cdot ({{{{{{\bf{r}}}}}}}_{i}-{{{{{{\bf{r}}}}}}}_{j})}{e}^{-i{{{{{\bf{q}}}}}}\cdot {{{{{{\bf{r}}}}}}}_{j}},$$
(1)

where q is the wavevector of the BO. In this study, we use the simplified BO form factor due to the nearest sites presented in Supplementary Note 2-1. The form factor at q = q1, \({g}_{{{{{{{\bf{q}}}}}}}_{1}}^{lm}\), is nonzero only when {l, m} = {A, B}, and we set \({g}_{{{{{{\bf{q}}}}}}}^{lm}={g}_{{{{{{{\bf{q}}}}}}}_{1}}^{lm}\) when q is in region I in Fig. 2b. In the same way, we set \({g}_{{{{{{\bf{q}}}}}}}^{lm}={g}_{{{{{{{\bf{q}}}}}}}_{2}}^{lm}\,({g}_{{{{{{{\bf{q}}}}}}}_{3}}^{lm})\) when q is in region II (III). \({g}_{{{{{{{\bf{q}}}}}}}_{2}}^{lm}\) (\({g}_{{{{{{{\bf{q}}}}}}}_{3}}^{lm}\)) is nonzero for {l, m} = {B, C} ({C, A}). This treatment is justified because the BO fluctuations strongly develop only for q ≈ qn in kagome metals. Furthermore, we use \({\bar{f}}_{{{{{{\bf{q}}}}}}}({{{{{\bf{k}}}}}})=({f}_{{{{{{{\bf{q}}}}}}}_{n}}({{{{{\bf{k}}}}}})+{f}_{{{{{{{\bf{q}}}}}}}_{n}}({{{{{\bf{k}}}}}}+{{{{{\bf{q}}}}}}-{{{{{{\bf{q}}}}}}}_{n}))/2\) for q ~ qn in the numerical study to improve the accuracy. Both BO and cLC form factors are Hermite \(\delta {t}_{ij}^{lm}={(\delta {t}_{ji}^{ml})}^{ * }\), which leads to the relation \({g}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}})={({g}_{-{{{{{\bf{q}}}}}}}^{ml}({{{{{\bf{k}}}}}}+{{{{{\bf{q}}}}}}))}^{ * }\)45.

Fig. 2: BO fluctuations and self-energy.
figure 2

a Expressions of the bond-order (BO) irreducible susceptibility \({\chi }_{g}^{0,lm{m}^{{\prime} }{l}^{{\prime} }}(q)\). b Obtained q-dependence for \({\chi }_{g}^{0,{{{{{\rm{A}}}}}}{{{{{\rm{B}}}}}}{{{{{\rm{A}}}}}}{{{{{\rm{B}}}}}}}(q)\), which takes the maximum at q = q1. c \({\hat{\chi }}_{g}(q)\) enlarged by the Hartree term of the electron-phonon interaction (2). [Note that \({\hat{\chi }}_{g}(q)\) is also enlarged by the Fock term of the off-site Coulomb interaction; see Supplementary Note 2-2]. d Self-energy induced by the BO fluctuations.

Full size image

To express the development of the bond order and fluctuations in kagome metals, we introduce the following effective BO interaction:

$${\hat{H}}_{{{{{{\rm{int}}}}}}}=-\frac{1}{2N}\mathop{\sum}\limits_{{{{{{\bf{q}}}}}}}\frac{v}{2}\,\,{\hat{O}}_{{{{{{\bf{q}}}}}}}^{g}\,\,{\hat{O}}_{-{{{{{\bf{q}}}}}}}^{g},$$
(2)

where \({\hat{O}}_{{{{{{\bf{q}}}}}}}^{g}\equiv {\sum }_{{{{{{\bf{k}}}}}},l,m,\sigma }{g}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}}){c}_{{{{{{\bf{k}}}}}}+{{{{{\bf{q}}}}}},l,\sigma }^{{{\dagger}} }{c}_{{{{{{\bf{k}}}}}},m,\sigma }\) is the BO operator45,58,59. and \(v\) is the effective interaction. We assume that the form factor \({g}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}})\) is normalized as \(\mathop{\max }\nolimits_{{{{{{\bf{k}}}}}},l,m}| {g}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}})|=1\) at each q, ı.e., \(| \delta {t}_{ij}^{{{{{{\rm{b}}}}}}}| \equiv 1/2\) for the nearest sites. Then, the maximum matrix element of BO interaction in Eq. (2) is v/2. The interaction (2) would originate from the combination of (i) the paramagnon-interference due to on-site U12, (ii) the bond-stretching phonon60, and (iii) the Fock term of off-site Coulomb interaction V9. In (i), Eq. (2) is induced by the spin-fluctuation-mediated beyond-RPA processes, whose diagrammatic expressions are shown in Fig. 3c, d in ref. 12. This processes give rise to the nematic BO in Fe-based SCs45. A great advantage of this theory12 is that the function of the BO form factor and the BO wavevector are automatically optimized to maximize TBO. Based on this theory, the BO at q = qn (n = 1, 2, 3) is robustly obtained based on the first principles multi-orbital model for CsV3Sb512. The effective parameter \(v\) in Eq. (2) is given as \(v\)AL ~ [gback + gum]/2, which is about 1.5 near the BO critical point (λbond ≲ 1), as we see in Fig. 3e of ref. 12. Thus, the value of \(v\) given by the AL processes is comparable to that used in the present study. In (ii), \({g}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}})\) is given by the hopping modulation due to the stretching mode and \(v\) = 2η2/ωD, where η is the electron-phonon (e-ph) coupling constant and ωD is the phonon energy at q ≈ qn. The BO interaction for the three vHS points model was derived in ref. 10. In (iii), \(v\) = 2V as we explain in the Supplementary Note 2-2. Thus, the effective interaction (2) is general. A possible driving forces of the BO have been discussed experimentally61,62.

Next, we study the susceptibility of the BO operator (per spin) defined as

$${\chi }_{g}({{{{{\bf{q}}}}}},{\omega }_{l})\equiv \frac{1}{2}\int\nolimits_{0}^{\beta }d\tau \left\langle {\hat{O}}_{{{{{{\bf{q}}}}}}}^{g}(\tau){\hat{O}}_{-{{{{{\bf{q}}}}}}}^{g}(0)\right\rangle {e}^{i{\omega }_{l}\tau }$$
(3)

where ωl is a boson Matsubara frequency. \({\hat{O}}_{{{{{{\bf{q}}}}}}}^{g}(\tau)\) is the Heisenberg representation of the BO operator. When \(v\) = 0, χg(q) is equivalent to the BO irreducible susceptibility \({\chi }_{g}^{0}(q)\)45,58:

$${\chi }_{g}^{0}(q)=\mathop{\sum}\limits_{lm{m}^{{\prime} }{l}^{{\prime} }}{\chi }_{g}^{0,lm{m}^{{\prime} }{l}^{{\prime} }}(q),$$
(4)
$${\chi }_{g}^{0,lm{m}^{{\prime} }{l}^{{\prime} }}({{{{{\bf{q}}}}}},{\omega }_{l})= \frac{T}{N}\mathop{\sum}\limits_{{{{{{\bf{k}}}}}},{\epsilon }_{n}}{g}_{{{{{{\bf{q}}}}}}}^{lm}{({{{{{\bf{k}}}}}})}^{*}{G}_{l{m}^{{\prime} }}({{{{{\bf{k}}}}}}+{{{{{\bf{q}}}}}},{\epsilon }_{n}+{\omega }_{l})\\ \times {G}_{{l}^{{\prime} }m}({{{{{\bf{k}}}}}},{\epsilon }_{n}){g}_{{{{{{\bf{q}}}}}}}^{{m}^{{\prime} }{l}^{{\prime} }}({{{{{\bf{k}}}}}}),$$
(5)

where q ≡ (q, ωl = 2πTl) and ϵn is a fermion Matsubara frequency. Equation (5) contains two form factors, so it vanishes when \(l=m\) or \({l}^{{\prime} }={m}^{{\prime} }\). Its diagrammatic expression is given in Fig. 2a. The numerical result for \({\chi }_{g}^{{{{{{\rm{0,ABAB}}}}}}}({{{{{\bf{q}}}}}},0)\) is shown in Fig. 2b, which exhibits the broad peak at the nesting vector between vHS-A and vHS-B; q = q1.

The BO susceptibility in Eq. (3) is strongly magnified by the Hartree term of Eq. (2) because of the same form factors in both equations. Its process is expressed in Fig. 2c, and its analytic expression is

$${\chi }_{g}(q)={\chi }_{g}^{0}(q)/(1-v{\chi }_{g}^{0}(q)).$$
(6)

where the notation q ≡ (q, ωl = 2πTl) is used. Here, the relation \({\chi }_{g}({{{{{{\bf{q}}}}}}}_{n},0)\propto {(1-{\alpha }_{{{{{{\rm{BO}}}}}}})}^{-1}\) holds, where \({\alpha }_{{{{{{\rm{BO}}}}}}}\equiv \mathop{\max }\nolimits_{{{{{{\bf{q}}}}}}}v{\chi }_{g}^{0}({{{{{\bf{q}}}}}})\) is the BO stoner factor. χg(qn, 0) diverges when αBO = 1. In contrast, the cLC susceptibility for the odd-parity cLC form factor, \({f}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}}-{{{{{\bf{q}}}}}}/2)=-{f}_{{{{{{\bf{q}}}}}}}^{ml}(-{{{{{\bf{k}}}}}}-{{{{{\bf{q}}}}}}/2)\), is unchanged by the Hartree term because g is even-parity.

The BO susceptibility is the largest in the Hartree-Fock (HF) approximation. As we discuss in the Supplementary Note 3. the BO and cLC susceptibilities at q = qn are \({\tilde{\chi }}_{g}\propto {(1-(v+{v}^{{\prime} }){\chi }_{g}^{0})}^{-1}\) and \({\tilde{\chi }}_{{{{{{\rm{cLC}}}}}}}\propto {(1-{v}^{{\prime\prime} }{\chi }_{g}^{0})}^{-1}\), respectively. Here, \(-{v}^{{\prime} } \sim {v}^{{\prime\prime} } \sim 0.3yv\) originates from the Fock term. (The coefficient \(y\)( ~ O(1)) depends on the origin of BO fluctuations. \(y\) = 1/2 for Himt in Eq. (2). The detailed discussion on \(y\) will be presented later). Thus, both susceptibilities are enlarged, while \({\tilde{\chi }}_{{{{{{\rm{cLC}}}}}}} \, < \, {\tilde{\chi }}_{g}\) within the HF approximation. However, we discover that \(-{v}^{{\prime} }\) and \(v^{{\prime\prime} }\) are further enlarged by the Maki-Thompson (MT) vertex corrections.

The MT term describes the scattering of electrons due to the developed bosonic fluctuations. This scattering process is important in metals near the quantum critical points. For example, in nearly antiferromagnetic metals, the d-wave SC transition is induced by the MT processes of spin fluctuations. In kagome metals, the MT term represents the strong inter vHS scattering of electrons mediated by the abundant BO fluctuations; see Fig. 3a. (The MT term also describes the s-wave SC state in kagome metals12). Here, we find that both \({\tilde{\chi }}_{g}\) and \({\tilde{\chi }}_{{{{{{\rm{cLC}}}}}}}\) are comparably enlarged due to the MT processes in the present theory.

Fig. 3: cLC order driven by BO fluctuation mechanism.
figure 3

a Density-wave (DW) equation due to the single exchange term of the bond-order (BO) fluctuations. b Eigenvalue of the DW equation λq (red-solid line) and BO Stoner factor \({\alpha }_{{{{{{\bf{q}}}}}}}^{{{{{{\rm{BO}}}}}}}\) (blue-dashed line) for \(v\) = 0.7 and T = 0.012. Both show peaks at q = qn. c Imaginary hopping modulation \({{{{{\rm{Im}}}}}}\delta {t}_{{{{{{\rm{AC}}}}}}}^{{{{{{\rm{c}}}}}}}(R)\). Its triple-q order gives the cLC pattern in (d). One can check that the clock-wise (anti-clock-wise) loop currents on hexagons (triangles) in (d) are inverted and moved by aAC under the sign change of \({f}_{{{{{{{\bf{q}}}}}}}_{3}}^{{{{{{\rm{AC}}}}}}}\).

Full size image

To understand the BO+cLC phase diagram and the energy scale of these orders accurately, we have to include the self-energy that describes the quasiparticle properties. We calculate the on-site self-energy due to BO fluctuations (see Eq. (8) in “Methods”). The fluctuation-induced self-energy is essential to reproduce the T-dependence of various physical quantities, as well-known in spin fluctuation theories63,64,65. Here, we calculate χg(q) in Eq. (6) and Σm(ϵn) in Eq. (8) self-consistently.

BO fluctuation-mediated cLC order

Next, we discuss the cLC mechanism. The HF approximation for the BO interaction (2) does not lead to the cLC order, as we explain in the Supplementary Note 3. (It is the same for off-site Coulomb interaction case; see Supplementary Note 2-1). Thus, the cLC order should be ascribed to the beyond-HF mechanism. Here, we explain that the strong electron scattering between different vHS points due to the BO fluctuations, which are described as the MT processes, causes the odd-parity cLC order \(\delta {t}_{ij}^{{{{{{\rm{c}}}}}}}=-\delta {t}_{ji}^{{{{{{\rm{c}}}}}}}\). (Note that the spin-fluctuation-exchange processes cause the cLC order in quasi-1D systems66). This process is generated by solving the following linearized DW equation41,51,66:

$${\lambda }_{{{{{{\bf{q}}}}}}}{f}_{{{{{{\bf{q}}}}}}}^{L}(k)= \frac{T}{N}\mathop{\sum}\limits_{p,{M}_{1},{M}_{2}}{I}_{{{{{{\bf{q}}}}}}}^{L,{M}_{1}}(k,p)\\ \times {\{-G(p)G(p+{{{{{\bf{q}}}}}})\}}^{{M}_{1},{M}_{2}}{f}_{{{{{{\bf{q}}}}}}}^{{M}_{2}}(p),$$
(7)

where \(L\equiv (l,{l}^{{\prime} })\) and Mi represent the pair of sublattice indices. \({I}_{{{{{{\bf{q}}}}}}}^{L,M}(k,p)\propto -{\chi }_{g}(k-p)\) is given by the BO fluctuation scattering process shown in Fig. 3a, which is called the MT process. The expression of \({I}_{{{{{{\bf{q}}}}}}}^{L,M}\) is given in Eq. (10) in “Methods” section. Note that T∑n{−G(p, ϵn)G(p + q, ϵn)} > 0.

By solving the DW equation (7), the optimized order parameter function is given as the eigenfunction \({f}_{{{{{{\bf{q}}}}}}}^{L}(k)\) for the maximum eigenvalue λq. \(\mathop{\max }\limits_{{{{{{\bf{q}}}}}}}\{{\lambda }_{{{{{{\bf{q}}}}}}}\}=1\) at the phase transition temperature. Note that \({f}_{{{{{{\bf{q}}}}}}}^{L}(k)\) represents the symmetry-breaking part in the normal self-energy \(\Delta \Sigma ({{{{{\bf{k}}}}}},{{{{{\bf{q}}}}}}) \sim \langle {c}_{{{{{{\bf{k}}}}}}+q\sigma }^{{{\dagger}} }{c}_{{{{{{\bf{k}}}}}}\sigma }\rangle \), and DW equation is directly derived from the stationary condition δF[ΔΣ]/δ(ΔΣ) = 051. We can regard the DW equation (7) as the gap equation for the optimized particle-hole (p-h) condensation45,51.

Note that the BO fluctuation-mediated interaction for the self-energy (Eq. (8)) and that for the kernel function (Eq. (10)) have the same coefficient \(y\), guaranteed by the Ward identity. The \(y\) depends on the BO fluctuation mechanism: \(y\) = 1/2 for the BO interaction \(v\) in Eq. (2) that works only in charge-channel. \(y\) ≈ 1/2 for the AL mechanism for the same reason12. \(y\) = 2 for the off-site V that induced both charge- and three spin-channel BO fluctuations as we explain in Supplementary Note 2-2. In kagome metals, both \(v\) and V coexist. In this case, BO fluctuations in charge-channel dominate over those in spin-channel, and therefore \(y\) ≳ 0.5 is expected in real kagome metals. Detailed explanation is given in Supplementary Note 2-3. Because we are interested in a general argument, we set \(y\) = 0.5 ~ 1 as a model parameter below. Note that the Aslaoazov-Larkin term is unimportant as we discuss in Supplementary Note 4.

Figure 3b shows the largest eigenvalue of the DW equation λq (red line) and BO Stoner factor \({a}_{{{{{{\bf{q}}}}}}}^{{{{{{\rm{BO}}}}}}}\equiv v{\chi }_{g}^{0}({{{{{\bf{q}}}}}})\) (blue line) as functions of q, for \(v\) = 0.7 and T = 0.012. They exhibit the maximum value at q = qn (n = 1, 2, 3). The corresponding solution of the DW equation is odd-parity: \({f}_{{{{{{\bf{q}}}}}}}^{lm}({{{{{\bf{k}}}}}}-{{{{{\bf{q}}}}}}/2)=-{f}_{{{{{{\bf{q}}}}}}}^{ml}(-{{{{{\bf{k}}}}}}-{{{{{\bf{q}}}}}}/2)\). Then, the corresponding real-space hopping modulation is odd-parity \(\delta {t}_{ij}^{{{{{{\rm{c}}}}}}}=-\delta {t}_{ji}^{{{{{{\rm{c}}}}}}}\) and pure imaginary when \(\delta {t}_{ij}^{{{{{{\rm{c}}}}}}}\) is Hermitian. The obtained \(\delta {t}_{{{{{{\rm{AC}}}}}}}^{{{{{{\rm{c}}}}}}}(R)\equiv \delta {t}_{{i}_{{{{{{\rm{A}}}}}}}{j}_{{{{{{\rm{C}}}}}}}}^{{{{{{\rm{c}}}}}}}\) for the cLC at q = q3 along the A-C direction is shown in Fig. 3c, where the odd integer R is defined as \({{{{{{\bf{r}}}}}}}_{i}^{{{{{{\rm{C}}}}}}}-{{{{{{\bf{r}}}}}}}_{j}^{{{{{{\rm{A}}}}}}}\equiv R{{{{{\bf{a}}}}}}\). In addition, the odd-parity relation \(\delta {t}_{{{{{{\rm{AC}}}}}}}^{{{{{{\rm{c}}}}}}}(R)=-\delta {t}_{{{{{{\rm{CA}}}}}}}^{{{{{{\rm{c}}}}}}}(-R)\) is verified. The obtained charge loop current pattern for the 3Q state is depicted in Fig. 3d.

Here, we discuss why the cLC order is mediated by the BO fluctuations. Let us consider the infinite series of MT terms in Fig. 4a, which is equal to \({f}_{{{{{{{\bf{q}}}}}}}_{3}}^{{{{{{\rm{AC}}}}}}}{({\lambda }_{{{{{{{\bf{q}}}}}}}_{3}}^{-1}-1)}^{-1}\) according to the DW equation (7). The first term together with other odd-order MT terms in Fig. 4a give the repulsive Umklapp interaction, \({\Gamma }_{{{{{{\rm{um}}}}}}}^{{{{{{\rm{MT}}}}}}} < 0\), which leads to the odd-parity order \({f}_{{{{{{{\bf{q}}}}}}}_{3}}^{{{{{{\rm{A}}}}}}{{{{{\rm{C}}}}}}}=-{f}_{{{{{{{\bf{q}}}}}}}_{3}}^{{{{{{\rm{C}}}}}}{{{{{\rm{A}}}}}}}\). In contrast, the second term together with other even-order MT terms give the attractive backward interaction, \({\Gamma }_{{{{{{\rm{back}}}}}}}^{{{{{{\rm{MT}}}}}}} \, > \, 0\), which gives the attraction among the same \({f}_{{{{{{{\bf{q}}}}}}}_{3}}^{lm}\). Therefore, all series of MT terms cooperatively induce the odd-parity current order form factor shown in Fig. 3c. Figure 4b exhibits the obtained \(v\)-dependence of the cLC eigenvalue as a function of T in the case of \(y\) = 1.

Fig. 4: cLC and BO transition temperatures and predicted phase diagram.
figure 4

a Series of Maki-Thompson (MT) processes produced in the density-wave (DW) equation. Yellow wavy lines represent the bond order (BO) propagators. The first-order and other odd-order terms give the repulsive Umklapp scattering \({\Gamma }_{{{{{{\rm{um}}}}}}}^{{{{{{\rm{MT}}}}}}} \, < \, 0\). The second-order and other even-order terms give the attractive backward scattering \({\Gamma }_{{{{{{\rm{back}}}}}}}^{{{{{{\rm{MT}}}}}}} \, > \, 0\). Both scatterings give the odd-parity charge loop current (cLC) order cooperatively. b Obtained \(v\)-dependence of the eigenvalue of cLC \({\lambda }_{{{{{{{\bf{q}}}}}}}_{3}}\) at \(y\) = 1. c Obtained T-dependence of \({\lambda }_{{{{{{{\bf{q}}}}}}}_{3}}\) at \(y\) = 1. αBO is shown by the color of each line, and the small black circle on each line represents TBO. The relation TcLC > TBO is satisfied in the weak-coupling region (\(v\) < \(v\)*). Obtained TcLC and TBO as functions of \(v\) for d \(y\) = 1 and e \(y\) = 0.5. f Schematic phase diagram in the present theory. The nematic 3Q BO+cLC coexisting phase appears when TcLC < TBO.

Full size image

Figure 4c exhibits the T-dependence of \({\lambda }_{{{{{{{\bf{q}}}}}}}_{3}}\) for \(v\) = 0.4−1.4 at \(y\) = 1. The cLC transition temperature TcLC is given by the relation \({\lambda }_{{{{{{{\bf{q}}}}}}}_{3}}=1\). The color on each line represents αBO: It is clearly seen that αBO at T = TcLC monotonically increases with \(v\). In 2D systems, αBO asymptotically approaches 1 with \(v\), but never exceeds 1 due to the χg-induced self-energy65,67,68. Here, TBO is defined as \({\alpha }_{{{{{{\rm{BO}}}}}}}={\alpha }_{{{{{{\rm{BO}}}}}}}^{*}\) with \({\alpha }_{{{{{{\rm{BO}}}}}}}^{*}=0.985\), which is shown as a small circle on each line in Fig. 4c, by considering the small inter-layer BO coupling ∣\(v\)⊥∣( ≪ \(v\)). (Overall results are unchanged for \({\alpha }_{{{{{{\rm{BO}}}}}}}^{*} \sim 0.99\)). The three-dimensional (3D) BO appears when \({\chi }_{g}^{{{{{{\rm{3D}}}}}}}={\chi }_{g}^{{{{{{\rm{2D}}}}}}}/(1-| {v}_{\perp }| {\chi }_{g}^{{{{{{\rm{2D}}}}}}})=\infty \), that is, ∣\(v\)⊥∣ ~ (1 − αBO)\(v\). Similar method is frequently used in deriving TSDW in spin fluctuation theories67,68. When \(v\) is small, the relation TcLC > TBO holds, which is natural because the MT term becomes large for αBO ≲ 1. With increasing \(v\), however, the opposite relation TcLC < TBO is realized due to the large self-energy effect.

The obtained TcLC and TBO as functions of \(v\) are shown in Fig. 4d \(y\) = 1 and e \(y\) = 0.5. In d, TcLC = TBO is realized at \(v\) = \(v\)* ≈ 1.03, and TcLC/TBO > 1 is realized in the weak-coupling region \(v\) < \(v\)*. The opposite relation TcLC/TBO < 1 is obtained in the strong-coupling region \(v\) > \(v\)* because the eigenvalue of DW equation (7) is suppressed by the large self-energy. In e, TcLC = TBO at \(v\) = \(v\)* ≈ 0.55.

Figure 4d, e indicate that both BO and cLC instabilities are comparable for \(v\)* ≈ \(v\). Based on the parity argument, the BO (cLC) instability is given by Γback + (−)Γum. Therefore, the relation Γback ≫ ∣Γum∣ should be satisfied for \(v\)* ≈ \(v\). In fact, the Hartree process gives positive \({\Gamma }_{{{{{{\rm{back}}}}}}}^{{{{{{\rm{H}}}}}}}={\Gamma }_{{{{{{\rm{um}}}}}}}^{{{{{{\rm{H}}}}}}} \sim v/(1-v{\chi }_{g}^{0})\), so the Hartree and MT processes strengthen each other in Γback but cancel each other in Γum. This relation is verified by the parquet RG study in Supplementary Note 5.

We discuss that Fig. 4d, e naturally explain the experimental P-T phase diagram with TBO and \({T}_{2}^{*}\,( \sim {T}_{{{{{{\rm{TRSB}}}}}}})\) given by μSR study15 for A = Rb, considering that \(v\)/Wband decreases with P. A schematic BO+cLC phase diagram derived from the present theory is depicted in Fig. 4f. (This schematic phase diagram is supported by the Ginzburg-Landau (GL) analysis in Supplemental Fig. 11a–c). The suppression of the secondary order due to the primary order is considered). The cLC phase is realized next to the BO phase because it is mediated by the BO fluctuations. This cLC+BO phase diagram is reminiscent of the SC-SDW phase diagram of spin-fluctuation-mediated superconductors, which has been reproduced by considering the self-energy67,68.

Z 3-nematic state given by the cLC-BO coexistence

To understand the cLC+BO coexisting states in Fig. 4f, the Ginzburg-Landau (GL) free energy analysis is very useful10,30,31. For example, the third-order GL term is F(3) = b1ϕ1ϕ2ϕ3 + b2(ϕ1η2η3 + η1ϕ2η3 + η1η2ϕ3), where the coefficients satisfy the relation b1 ~ −b2, and (ϕ1, ϕ2, ϕ3) [(η1, η2, η3)] is the magnitude of the BO [cLC] parameter at q = q1, q2, q3. Here, we introduce the 3Q states \({{{{{{\boldsymbol{\phi }}}}}}}_{1}=(\phi /\sqrt{3})(1,1,1)\), \({{{{{{\boldsymbol{\eta }}}}}}}_{1}=(\eta /\sqrt{3})(1,1,1)\), and \({{{{{{\boldsymbol{\eta }}}}}}}_{2}=(\eta /\sqrt{3})(1,-1,-1)\). The chiral center of η1 coincides with the center of the BO ϕ1, while the center of η2 is shifted by aBA from that of ϕ1. Thus, the coexisting state (ϕ1, η1[2]) has the C6[2]-symmetry as shown in Fig. 5a [c], and its FS in the folded Brillouin zone is in Fig. 5b [d]. As we explain in Supplementary Note 6-1, F(3) for the C2-coexisting state is lower than that for the C6-coexisting state in the case of ∣ϕ∣ ≫ ∣η∣ for fixed ∣ϕ∣, ∣η∣. (The optimized cLC order in the C2-phase is \({{{{{{\boldsymbol{\eta }}}}}}}_{2}^{{\prime} }\propto (2,-1,-1)\); see Supplementary Note 6-1). Therefore, the BO+cLC Z3 nematic state is realized when \(v\) > \(v\)*. This result is consistent with the recent observation of out-of-phase combination of bond charge order and loop currents by STM measurement69. We comment that the nematic BO+cLC phase is obtained when TBO ≳ TcLC by minimizing the GL free energy F[ϕ, η] exactly in Supplementary Note 6-2.

Fig. 5: cLC+BO coexisting states with C6 and C2 symmetries.
figure 5

a C6-symmetric bond order (BO) and charge loop current (cLC) coexisting state in real space. b Its folded Fermi surface (FS). The folded Brillouin zone is shown by dotted lines. c C2-symmetric BO+cLC coexisting state. The nematicity originates from the out-of-phase combination of bond and current orders. d Its nematic FS. Here, the director is parallel to aBA because the cLC order in (c) is shifted by aBA from the cLC order in (a). Thus, the Z3 nematic state with different three directors is realized. Here, we use large \(| \delta {t}_{ij}^{{{{{{\rm{b,c}}}}}}}| \,(=0.05)\) to exaggerate the nematicity.

Full size image

We also discuss the case of \(v\) < \(v\)*, where cLC is the primary order as shown in Fig. 4f. The C6 symmetry 3Q cLC order appears at T = TcLC when 2d2,a/d2,b > 1 as we discuss in Supplementary Note 5, where d2,a (d2,b) is the GL coefficients of the \({\eta }_{1}^{4}\) (\({\eta }_{1}^{2}{\eta }_{2}^{2}\)) term. Note that the primary 3Q cLC order induces the secondary BO parameter even above TBO through the b2-term in F(3) 10. In contrast, the 1Q cLC state is realized when 2d2,a/d2,b is smaller than unity. Thus, the electronic state becomes nematic at TcLC(>TBO). In this case, there is no secondary BO component above TBO. Recently, strong evidence of the emergence of the 1Q cLC state at ~130 K (>TBO) has been reported by the magnetic torque measurement21.

The obtained nematic BO+cLC state is TRSB and two-dimensional. Other possible nematic state is the shift-stacking of the 3Q BO layers, each of which has C6 symmetry. The shift-stacking is caused by the 3Q state composed of the 3D BO at \({{{{{{\bf{q}}}}}}}_{n}^{{{{{{\rm{3D}}}}}}}\) with \({q}_{1,z}^{{{{{{\rm{3D}}}}}}}={q}_{2,z}^{{{{{{\rm{3D}}}}}}}=\pi \) and \({q}_{3,z}^{{{{{{\rm{3D}}}}}}}=0\)10. We stress that these two different nematic states can be realized at different temperatures.

Giant AHE in cLC+BO state

Next, we discuss the transport phenomena that originate from the cLC24,70. Using the general expression of the intrinsic conductivity65,71, we calculate the Hall conductivity (σxy and σyx) due to the Fermi-surface contribution in the BO+cLC state. The expression is \({\sigma }_{\mu \nu }=\frac{1}{N}{\sum }_{{{{{{\bf{k}}}}}}}{A}_{\mu \nu }({{{{{\bf{k}}}}}})\), where \({A}_{\mu \nu }({{{{{\bf{k}}}}}})=\frac{{e}^{2}}{\hslash }\frac{1}{\pi }{{{{{\rm{Tr}}}}}}\{{\hat{v}}_{{{{{{\bf{k}}}}}},\mu }{\hat{G}}_{{{{{{\bf{k}}}}}}}(i\gamma){\hat{v}}_{{{{{{\bf{k}}}}}},\nu }{\hat{G}}_{{{{{{\bf{k}}}}}}}(-i\gamma)\}\). Here, \({\hat{G}}_{{{{{{\bf{k}}}}}}}(\epsilon)={((\epsilon+\mu)\hat{1}-{\hat{h}}_{{{{{{\bf{k}}}}}}})}^{-1}\) is the Green function matrix, where \({\hat{h}}_{{{{{{\bf{k}}}}}}}\) is the 12 × 12 tight-binding model with the 3Q BO and cLC order, and \({\hat{v}}_{{{{{{\bf{k}}}}}},\mu }=d{\hat{h}}_{{{{{{\bf{k}}}}}}}/d{k}_{\mu }\) is the velocity operator. γ( > 0) is the electron damping rate that is given by the imaginary part of the self-energy. We set n = nvHS and \(| \delta {t}_{ij}^{{{{{{\rm{b}}}}}}}|=| \delta {t}_{ij}^{{{{{{\rm{c}}}}}}}|=0.025\), where the band hybridization gap due to the BO+cLC order is about \(\Delta \, \approx \, 2\sqrt{| \delta {t}_{ij}^{{{{{{\rm{b}}}}}}}{| }^{2}+| \delta {t}_{ij}^{{{{{{\rm{c}}}}}}}{| }^{2}}=0.07\).

Figure 6a shows the obtained conductivities in the nematic BO+cLC state, in the unit of \(\frac{{e}^{2}}{\hslash }\,(=2.4\times 1{0}^{-4}{\Omega }^{-1})\). When γ ≪ Δ, the Hall conductivity \({\sigma }_{{{{{{\rm{H}}}}}}}\equiv \frac{1}{2}({\sigma }_{xy}-{\sigma }_{yx})\) is almost constant, and its magnitude is proportional to \(| \delta {t}_{ij}^{{{{{{\rm{c}}}}}}}| \). When γ ≫ Δ, in contrast, σH decreases with γ in proportion to γ−2. This crossover behavior is universal in the intrinsic Hall effect, which was first revealed in heavy fermion systems, and found to be universal in later studies65,71,72,73. Note that \(\frac{1}{2}({\sigma }_{xy}+{\sigma }_{yx})\) is nonzero in the nematic state. To understand the origin of the intrinsic Hall effect, we plot AH(k) ≡ (Axy(k) − Ayx(k))/2 at γ = 0.05 in Fig. 6b: It shows a large positive value mainly around the vHS points, due to the band-hybridization induced by the cLC order. The obtained σH ~ 1 corresponds to 4 × 103 Ω−1 cm−1 because the interlayer spacing is ~0.6 nm. Thus, giant AHE σH ~ 102 Ω−1 cm−1 reported in refs. 25,26 is understood in this theory.

Fig. 6: Giant AHE in nematic cLC+BO state.
figure 6

a Anomalous Hall conductivity in the nematic BO+cLC state (\(| \delta {t}_{ij}^{{{{{{\rm{b,c}}}}}}}|=0.025\)) as a function of the electron damping rate γ  ∝ τ−1. The full (broken) line represents the positive (negative) value. Thus, the Hall conductivity σH becomes large in the low-resistivity region (γ ≲ 0.03). b AH(k) at γ = 0.05: \(\frac{1}{N}{\sum }_{{{{{{\bf{k}}}}}}}{A}_{{{{{{\rm{H}}}}}}}({{{{{\bf{k}}}}}})={\sigma }_{{{{{{\rm{H}}}}}}}\).

Full size image

Parquet RG theory, field-induced cLC mechanism

To verify the idea of the BO fluctuation-mediated cLC, we perform the analysis of the parquet renormalization group (RG) formulation10,74 and present the results in Supplementary Note 5. A great merit of the RG method is that both particle-particle and particle-hole channels are treated on the same footing. We find that both BO and cLC fluctuations cooperatively develop in Supplementary Fig. 6. This result of the RG study strongly supports the validity of the DW equation analysis.

We comment on the complementary relationship between the present theory and the parquet RG theory. The latter theory solves a simplified 3-patch model in an unbiased way, leading to the development of both cLC order and BO, while the relationship between the two orders is not clear. On the other hand, the present theory focuses on the existence of the experimental BO phase and reveals that abundant BO quantum fluctuations lead to TRSB particle-hole condensation. Thus, the concept of the BO fluctuation-mediated cLC has been verified based on different reliable theories.

Summary

In summary, we proposed a cLC mechanism mediated by the BO fluctuations in kagome metals. This cLC mechanism is universal because it is independent of the origin of the BO. The validity of the idea of the BO fluctuation-mediated cLC has been confirmed by the parquet RG study in Supplementary Note 5. Furthermore, we revealed that novel Z3 nematicity emerges under the coexistence of the cLC and the BO reported in refs. 5,19,27 in addition to the giant AHE25,26. This theory presents a promising scenario for understanding the BO, the cLC and the nematicity in kagome metals in a unified way.

In the present study, we focus on the pure-type band composed of b3g-orbital. However, the impact of other 3d-orbitals on the cLC order has also been studied in refs. 30,75. The extension of the present theory to multi-orbital models is a very important future issue.

Here, we shortly discuss several experimental evidences of the BO+cLC coexistence. Recent transport measurement of highly symmetric fabricated CsV3Sb5 micro sample76 reveals that small magnetic field hz (<10T) or small strain gives rise to the nematic BO+cLC coexisting state below TBO. This finding is well explained by the recent GL theory under hz77: The current-bond-hz trilinear coupling caused by the orbital magnetization gives rise to the sizable hz-induced cLC order in the BO state. This theory also explains hz-induced enhancement of the cLC order observed by μSR measurements15,17,18 and field-tuned chiral transport study20. It is noteworthy that the nematic electronic state that supports the C2 BO+cLC order in Fig. 5c has been reported by recent STM measurement69.

Finally, we comment on some interesting kagome metals other than AV3Sb5. Double-layer kagome metal ScV6Sn6 shows \(\sqrt{3}\times \sqrt{3}\) charge-density wave (CDW)78. It was proposed that the CDW originates from the flat phonon modes with Sn vibrations79,80. Interestingly, ScV6Sn6 also exhibits the spontaneous TRSB state81. The mechanism of the TRSB state in ScV6Sn6 is an interesting future problem. (Note that the existence of the vHS points is not a requirement for the cLC order66). The GL free energy analysis was performed in ref. 31. Recently, very weak but definite signal of the nematic electronic order has been observed in Ti-based kagome metal CsTi3Bi582,83. To explain the observed hidden nematicity, the odd-parity BO without TRSB has been predicted theoretically84.

Methods

Self-energy due to BO fluctuations

To understand the BO+cLC phase diagram and the energy scale of these orders accurately, we have to include the self-energy that describes the quasiparticle properties. We calculate the on-site self-energy due to BO fluctuations as

$${\Sigma }_{m}({\epsilon }_{n})= \frac{T}{N}\mathop{\sum}\limits_{{{{{{\bf{k}}}}}},q,{m}^{{\prime\prime} },{m}^{{\prime\prime} {\prime} }}{G}_{{m}^{{\prime\prime} }{m}^{{\prime} }}({{{{{\bf{k}}}}}}+{{{{{\bf{q}}}}}},{\epsilon }_{n}+{\omega }_{l})\\ \times {B}_{mm,{m}^{{\prime\prime} }{m}^{{\prime} }}(k,q),$$
(8)
$${B}_{mm,{m}^{{\prime\prime} }{m}^{{\prime} }}(k,q)={g}_{{{{{{\bf{q}}}}}}}^{{m}^{{\prime} }m}({{{{{\bf{k}}}}}}){g}_{{{{{{\bf{q}}}}}}}^{{m}^{{\prime\prime} }m}{({{{{{\bf{k}}}}}})}^{ * }\cdot yv(1+v{\chi }_{g}(q)),$$
(9)

which is shown in Fig. 2d. Then, the Green function is given as \(\hat{G}(k)={(i{\epsilon }_{n}+\mu -\hat{h}({{{{{\bf{k}}}}}})-\hat{\Sigma }({\epsilon }_{n}))}^{-1}\). The effect of thermal fluctuations described by the self-energy is essential to reproduce the T-dependence of various physical quantities. Here, \(y\) = 1/2 when \({\hat{H}}_{{{{{{\rm{int}}}}}}}\) is given in Eq. (2). In the present numerical study, we calculate \({\chi }_{g}(q)={\chi }_{g}^{0}(q)/(1-v{\chi }_{g}^{0}(q))\) and Σm(ϵn) in Eq. (8) self-consistently.

Kernel function of the DW equation

The kernel function due to BO fluctuations in Eq. (7) is given as

$${I}_{{{{{{\bf{q}}}}}}}^{l{l}^{{\prime} } ,m{m}^{{\prime} }}(k,p)= -{g}_{{{{{{\bf{p}}}}}}-{{{{{\bf{k}}}}}}}^{{m}^{{\prime} }{l}^{{\prime} }}({{{{{\bf{k}}}}}})yv(1+v{\chi }_{g}(k-p)){g}_{{{{{{\bf{k}}}}}}-{{{{{\bf{p}}}}}}}^{lm}({{{{{\bf{p}}}}}}+{{{{{\bf{q}}}}}})\\ +{g}_{{{{{{\bf{q}}}}}}}^{l{l}^{{\prime} }}({{{{{\bf{k}}}}}})v{g}_{{{{{{\bf{q}}}}}}}^{m{m}^{{\prime} }}{({{{{{\bf{p}}}}}})}^{*},$$
(10)

which is expressed in Fig. 3a and Supplementary Fig. 4a. The first term, the MT term, is important when αBO ≲ 1, and its first term is the Fock term. The second term, the Hartree term, vanishes when the eigenfunction \({\hat{f}}_{{{{{{\bf{q}}}}}}}(k)\) is orthogonal to the BO form factor \({\hat{g}}_{{{{{{\bf{q}}}}}}}(k)\), like the cLC order is. Note that \(\hat{B}(k,q)=-{\hat{I}}_{{{{{{\boldsymbol{0}}}}}}}(k,k+q)\). A more detailed discussion is presented in Supplementary Note 3.

Numerical analysis

In this study, we solved the eigenvalue equation with the kernel function (7) and the integral equations (8) and (9) numerically, by dividing the Brillouin zone into 60 × 60k meshes. The number of k meshes is fine enough to achieve reliable numerical accuracy (~1%) at the calculated temperatures (T ~ 0.01).