First principles theory for the magnetic and charge instabilities in AV₃Sb₅ systems

The AV₃Sb₅ (A=K, Rb, Cs) system crystallizes in a layered kagome lattice structure as shown in Fig. 1(a). In each primitive cell there are three V sublattices forming kagome plane, within which one Sb atom is intercalated into each primitive cell. There are two layers of Sb atoms forming honeycomb lattices above and below the vanadium kagome planes, forming V-Sb hexahedra frameworks. The alkali metal are further intercalated between the V-Sb layers, donating electrons and making this system highly metallic. In Fig. 1(b) we show the band structure of RbV₃Sb₅ projected onto V $3d$ orbitals and Sb $5p$ orbitals in the pristine lattice structure which are calculated by first principles density functional theory (DFT). Clearly there are Dirac-like dispersions near $K$ and $H$ points contributed by V $d$ orbitals, van Hove singularities near $M$ point contributed by both V $d$ and Sb $p$ orbitals, as well as electron pockets centered at $\Gamma$ which are mostly contributed by Sb $p$ orbitals. Such complicated Fermi surfaces and orbital characters make it difficult to realistically describe the electronic structure of this system using a simplified lattice model with few orbitals. In Fig. 1(c) we present the phonon dispersions calculated by DFT, where we see two unstable phonon modes at the $M$ and $L$ points, consistent with previous reports. A linear combination of the unstable modes at the three equivalent $M$ points result in an inverse-star-of-David structure with $2\times 2\times 1$ supercell as shown in the right panel of Fig. 1(a). The band structures of the system are strongly re-constructed in the ISD lattice structure as clearly shown by the spectral function unfolded to the original primitive Brillouin zone (Fig. 1(d)), and new van Hove singularities emerge as a result of the lattice deformation Cho et al. (2021).

The presence of van Hove singularities both in the pristine and in the ISD structures indicate that the non-interacting Fermi surfaces may be unstable to $e\textrm{-}e$ Coulomb interactions. In order to study the interaction effects on the complicated Fermi surfaces with diverse orbital characters as shown in Fig. 1(b), we construct a realistic tight-binding model in the Wannier-function basis using the Wannier90 code package Mostofi et al. (2008). To be specific, the Bloch functions generated from DFT calculations are projected onto 60 spinor Wannier functions (for each pristine primitive cell) including all the V 3$d$ orbitals and Sb $5p$ orbitals, from which we obtain the following Wannier tight-binding model denoted by $H_{0}$:

$$H_{0}=\sum_{i\alpha\sigma,j\beta\sigma^{\prime}}h_{i\alpha\sigma,j\beta\sigma^{\prime}}\,\hat{c}^{\dagger}_{i\alpha\sigma}\,\hat{c}_{j\beta\sigma^{\prime}}$$

(1)

where the sets of indices $\{i,j\}$, $\{\alpha,\beta\}$, $\{\sigma,\sigma^{\prime}\}$ denote, in turn, all the V and Sb atomic sites, orbitals and spin degrees of freedom, and $h_{i\alpha\sigma,j\beta\sigma^{\prime}}$ denotes the real-space hopping amplitude including spin-orbit coupling effects . $\hat{c}^{\dagger}_{i\alpha\sigma}$ and $\hat{c}_{j\beta\sigma^{\prime}}$ are the electron creation and annihilation operators. The Wannier interpolated band structures can exactly reproduce those obtained from DFT calculations within an energy window extending from 2 eV below to 2 eV above the Fermi energy. We include the on-site multi-orbital Coulomb interactions of the Kanamori type for the V 3$d$ orbitals Kanamori (1963),

	$\displaystyle H_{\textrm{K}}=$	$\displaystyle U\sum_{i,\alpha}\hat{n}_{i\alpha\uparrow}\hat{n}_{i\alpha\downarrow}+U^{\prime}\sum_{i,\alpha<\beta,\sigma,\sigma^{\prime}}\hat{n}_{i\alpha\sigma}\hat{n}_{i\beta\sigma^{\prime}}\;$
		$\displaystyle-J_{\textrm{H}}\sum_{i,\alpha<\beta,\sigma,\sigma^{\prime}}\hat{c}^{\dagger}_{i\alpha\sigma}\hat{c}^{\vphantom{\dagger}}_{i\alpha\sigma^{\prime}}\hat{c}^{\dagger}_{i\beta\sigma^{\prime}}\hat{c}^{\vphantom{\dagger}}_{i\beta\sigma}\;$
		$\displaystyle+J_{\textrm{P}}\sum_{i,\alpha<\beta,\sigma}\hat{c}^{\dagger}_{i\alpha\sigma}\hat{c}^{\dagger}_{i\alpha-\sigma}\hat{c}^{\vphantom{\dagger}}_{i\beta\sigma}\hat{c}^{\vphantom{\dagger}}_{i\beta-\sigma}\;,$		(2)

where $U$ and $U^{\prime}$ are the intra-orbital and inter-orbital direct Coulomb interactions. $J_{\textrm{H}}$ and $J_{\textrm{P}}$ denote the Hunds’ coupling and pair hoppings respectively, and $\hat{n}_{i\alpha\sigma}=\hat{c}^{\dagger}_{i\alpha\sigma}\hat{c}_{i\alpha\sigma}$ is the density operator. We set $U\!=\!3.6\,$eV and $J\!=\!J_{\textrm{H}}\!=\!0.72\,$eV throughout this paper, and we let $U^{\prime}\!=\!U-J$ and $\!J_{\textrm{P}}\!=0$ such that Eq. (2) has full rotational invariance Georges et al. (2013). We further consider the density-density inter-site interactions between the neighboring V-Sb sites and V-V sites, the interaction amplitudes are denoted by $V$ and $W$ respectively,

$\displaystyle H_{\textrm{V-Sb}}=V\,\sum_{\langle ij\rangle}\sum_{\alpha\beta\sigma\sigma^{\prime}}\,\hat{c}^{\dagger}_{i\alpha\sigma}\hat{c}^{\dagger}_{j\beta\sigma^{\prime}}\hat{c}_{j\beta\sigma^{\prime}}\hat{c}_{i\alpha\sigma}\;$

where “$\sum_{\langle ij\rangle}$” means a summation over all the first-neighbor V-Sb sites (including both intralayer ones and interlayer ones), and the symbol “$\sum_{\llangle ij\rrangle}$” means summing over all the in-plane nearest neighbor V-V sites. In this work $V$ and $W$ are treated as two free parameters which vary from 0 to 1 eV. The full Hamiltonian is then $H_{0}+H_{\textrm{K}}+H_{\textrm{V-Sb}}+H_{\textrm{V-V}}$.

As discussed above, the system has a quite complicated primitive cell with multiple sublattice, orbital, and spin degrees of freedom. For the 60-band tight-binding model, there are 3600 independent matrix elements of the density operator $\rho_{i\alpha\sigma,j\beta\sigma^{\prime}}=\langle\hat{c}^{\dagger}_{i\alpha\sigma}\hat{c}_{j\beta\sigma^{\prime}}\rangle$ (consider $i,j$ as sublattice sites within the same primitive cell), each of which can be considered as an “order parameter” that break certain symmetries of the system and bring about Fermi surface instability. The actual ground state may be characterized by a complicated linear combination of the density matrix elements, which is difficult to guess a priori due to such complexity. Therefore, in order to find the leading Fermi-surface instability mode, we propose to calculate the interaction-renormalized generalized susceptibility tensor defined in the sublattice-orbital-spin space, the tensor element of which is denoted as $\chi(\mathbf{q})_{\mu\nu,\mu^{\prime}\nu^{\prime}}$, with $\mu\equiv\{i\alpha\sigma\}$ being a composite index referring to all the sublattice, orbital, and spin degrees of freedom (so do $\mu^{\prime},\nu$, and $\nu^{\prime}$). Suppose there is a weak static field with wavevector $\mathbf{q}$ that is only coupled with some given element of the density operator $\rho_{\mu^{\prime}\nu^{\prime}}(\mathbf{q})$ (denoted as $A_{\mu^{\prime}\nu^{\prime}}(\mathbf{q})$) is applied to the system, a change of density matrix element would be induced by the field, i.e., $\delta\rho_{\mu\nu}(\mathbf{q})=\chi(\mathbf{q})_{\mu\nu,\mu^{\prime}\nu^{\prime}}\,A_{\mu^{\prime}\nu^{\prime}}(\mathbf{q})$. One can further consider a field $\mathbf{A}^{\kappa}$ that is only coupled to (or “parallel” to) the $\kappa$th eigenmode $u_{\mu\nu}^{\kappa}(\mathbf{q})$ of the susceptibility tensor, i.e., $\mathbf{A}^{\kappa}(\mathbf{q})\sim u_{\mu\nu}^{\kappa}(\mathbf{q})$, then the system would respond to such a field along the eigenvector direction via

$$\delta\rho_{\mu\nu}(\mathbf{q})\sim\lambda^{\kappa}(\mathbf{q})u_{\mu\nu}^{\kappa}(\mathbf{q})$$

(3)

where $\lambda^{\kappa}(\mathbf{q})$ and $u_{\mu\nu}^{\kappa}(\mathbf{q})$ are the eigenvalues and eigenvectors of the susceptibility tensor $\chi(\mathbf{q})$. If any of the eigenmodes has a diverging eigenvalue $\lambda(\mathbf{q})^{\kappa}$, it means that the system would develop a spontaneous order parameter characterized by the diverging eigenmode $u_{\mu\nu}^{\kappa}(\mathbf{q})$ upon the application of an infinitesimal field, and would enter a spontaneous symmetry-breaking state through a second order phase transition. Note that similar idea has been used in the studies of surface-state instabilities of nodal-line semietals Liu and Balents (2017) and octupolar instabilities of RuCl₃ etc. Qiu et al. (2021); Uehara et al. (2015)

Following the above argument, we calculate the generalized susceptibility tensor for the AV₃Sb₅ (A=K, Rb, Cs) system based on the Wannierized tight-binding Hamiltonian given by Eqs. (1)-(LABEL:eq:intersite), and the effects of the on-site and inter-site Coulomb interactions are treated by random phase approximation (RPA), the details of which are presented in Supplementary Information sup . For clarity’s sake, in main text we only present the results for RbV₃Sb₅, and the results of the other two compounds are given in Supplementary Information sup . In Fig. 2(a) we present the 10 eigenmodes with maximal eigenvalues (of the RPA susceptibility) as a function of the wavevector $\mathbf{q}$, with on-site interaction only and the inter-site interactions being temporarily turned off. We see that there are three diverging modes only at $\Gamma$ point. These diverging modes at $\Gamma$ point turn out to be three spin magnetic modes. The inclusion of inter-site interactions does not change the diverging nature of the three spin ferromagnetic modes.On the other hand, since the system undergoes a CDW transition forming an ISD structure around 100 K, it is worthwhile to study the microscopic origin of such CDW transition. In such a context, we also carefully study the eigenmodes of the RPA susceptibility tensor at the three $M$ points, and find that there are no notable unstable modes. We further include V-Sb and V-V inter-site interactions, and find that the eigenvalues of the susceptibility tensor at $M$ points only increase by small amounts due to the inclusion of intersite interactions, as clearly shown in Fig. 2(b), where the different markers represent data sets for different choices of $V$ and $W$ values. Therefore, the CDW phase is likely driven by phonon instability and electron-phonon couplings instead of $e$-$e$ interactions.

In order to verify the conclusions obtained from the susceptibility calculations, we further perform un-restricted self-consistent Hartree-Fock (HF) calculations including both the on-site Kanamori interaction (Eq. (2)) and the inter-site interactions (Eq. (LABEL:eq:intersite)), with the nonlocal exchangee effects being treated exactly. We first perform HF calculations for the pristine primitive cell, and we take the three leading instability modes at $\Gamma$ point (Fig. 2(a)) as initial anstza for the density matrices in the HF calculations. The HF ground states turn out to be three spin ferromagnetic states with the magnetizations primarily pointing along the crystalline $x$, $y$, and $z$ directions with slight cantings, which are denoted as FM_x, FM_y, and FM_z states respectively. With on-site interactions only ($U\!=\!3.6\,$eV, $J\!=\!0.72\,$eV, $W\!=\!V\!=\!0$), the spin magnetization is as large as $\sim\!0.6\,\mu_{\textrm{B}}$ per V atom. Real-space current loops are generated in the ferromagnetic states as a result of spin-orbit coupling, which are shown in the three panels of Fig. 3(a) for the FM_y, FM_z, and FM_x phases. The red (double) arrows indicate interlayer V-Sb currents, while the black arrows indicate V-V in-plane currents. The maximal current amplitudes $\sim 800\,$nA, generating an orbital moment $\sim\!0.2\mu_{\textrm{B}}$, e.g., for the hexagonal current loop in the middle panel of Fig. 3(a). The filling dependence of the calculated anomalous Hall conductivity (AHC) $\sigma_{xy}$ in the FM_z state is shown in Fig. 3(b), where we see a peak of $\sigma_{xy}\!\sim\!800\,$S/cm when the Fermi energy is assumed to be $\sim\!0.15\,$eV above the physical Fermi energy; otherwise $\sigma_{xy}$ is quite small $\sim\!10\textrm{-}50\,$S/cm. The HF band structures (with on-site interactions only) are shown in Fig. 3(c).

The results presented in Fig. 3 are generated for the case with on-site interactions only. The inclusion of inter-site interactions, e.g., with $V\!=\!1\,$eV and $W\!=\!0$, would significantly suppress the spin ferromagnetism, resulting in a spin moment $\sim\!0.06\,\mu_{\textrm{B}}$ per V atom, and the corresponding AHC is two orders of magnitudes smaller $\sim\!0.5-5\,$S/cm. This can be interpreted from the charge transfer induced by the inter-site Sb-V Coulomb interactions. Once the inter-site interaction is taken into account, charges tend to transfer from the V to the Sb sites to minimize the Hartree energy, resulting in a charge occupation $\sim\!0.3\,$ per V atom for $V\!=\!\,1$eV, which is reduced by a factor of 6 compared with the case when $V\!=\!0$. As a result, the spin magnetization driven by the strong on-site Coulomb interactions of the V sites is also significantly diminished. Moreover, the energies differences of the three FM states are incredibly close to each other $\sim 0.01\,$meV, which imply strong magnetic fluctuations in real space. Therefore, the giant anomalous Hall effect observed in experiments Yang et al. (2020) likely arises from skew scatterings due to strong real-space magnetic fluctuations Ishizuka and Nagaosa (2021). The detailed results for the HF calculations including inter-site interactions with pristine primitive cell are presented in Supplementary Information sup .

We have also performed un-restricted self-consistent HF calculations for the $2\!\times\!2\!\times\!1$ ISD lattice structure with a realistic tight-binding model including 240 Wannier orbitals, including both the on-site and inter-site interactions as given by Eqs. (1)-(LABEL:eq:intersite), and the non-local exchange effects are treated exactly. We take the three spin ferromagnetic modes at $\Gamma$ point and four leading modes at $M$ point (obtained from RPA susceptibility calculations) as our initial ansatz for the HF calculations. The ground states, again, involve three nearly degenerate spin ferromagnetic states with their magnetizations pointing along the crystalline $x$, $y$, and $z$ directions respectively, with the energy difference $\sim 0.01\,$meV. Similar to the case of pristine primitive cell, the magnetizations in the $2\!\times\!2\!\times\!1$ ISD structure are on the order of $1\,\mu_{\textrm{B}}$ if only on-site interactions are included, and can be suppressed by one order of magnitude due to the inclusion of inter-site interactions. Interestingly, in the $2\!\times\!2\!\times\!1$ ISD structure, the magnetization distribution is inhomogeneous within the ISD primitive cell: the spin magnetizations are mostly concentrated at the six V atoms forming two contracted kagome triangles, and the other six V atoms forming a hexagon have nearly vanishing spin magnetizations, as schematically shown in Fig. 4(b) for the FM_z phase.

As in the case of the pristine cell, current loops are generated in the FM states as a result of spin-orbit coupling. The current-loop patterns in the three FM states are shown in Fig. 4(a) for the FM_y, FM_z, and FM_x states respectively. For example, for the FM_z state, there are three current loops: two of them are counter-clockwise circulating around the two kagome triangles, and the other one is clockwise circulating around the hexagon. With the intersite interaction $V\!=\!1\,$eV, the current amplitudes are very small $\sim 80\,$nA, which can generate a small orbital moment $\sim 0.02\,\mu_{\textrm{B}}$ for each current loop. We have further calculated the anomalous Hall conductivities of the three FM states in the ISD structure. The inclusion of V-Sb inter-site Coulomb interactions significantly suppress ferromagnetism and anomalous Hall effect (AHE), and the calculated AHC $\sigma_{xy}$ in the FM_z state $\sim 1-20\,$S/cm as shown in Fig. 4(c). Again, this indicates an incredibly small intrinsic contribution to AHE, and the giant AHC measured in experiments Yang et al. (2020) are likely due to skew scatterings from the strong magnetic fluctuations in real-space by virtue of the near degeneracy of the three ferromagnetic states Ishizuka and Nagaosa (2021). In Fig. 4(d), we have shown the calculated HF band structures of the $2\!\times\!2\!\times\!1$ ISD structure (plotted in the supercell Brillouin zone) with intersite interaction V=1.0 eV, from which we see a strong renormalization effects on the electronic structures due to the inclusion of nonlocal interactions. Notably, the Fermi surfaces are mostly contributed by two bands with nearly compensated electron-type and hole-type carriers due to the onset of magnetism, which is consistent with the Hall-effect measurements below $\sim 30\,$K Yang et al. (2020). Our results indicate that there may be two transitions in the AV₃Sb₅ (A=K, Rb, Cs) system, the first one is the phonon-driven CDW transition, and the other is the magnetic transition (likely with lower $T_{c}$) marked by the emergence of giant anomalous Hall effect and the reconstruction of Fermi surfaces due to the onset of weak spin ferromagnetism.

To summarize, in this work we have comprehensively studied the magnetic and charge instabilities of the AV₃Sb₅ system based on first principles tight-binding models including effects of both on-site and inter-site Coulomb interactions. Results from both RPA susceptibility calculations and un-restricted Hartree-Fock calculations indicate that the leading Fermi-surface instabilities driven by $e$-$e$ interactions are three nearly degenerate spin ferromagnetic modes regardless whether the system is in the pristine structure or ISD structure. The competition between the nearly degenerate ferromagnetic states can lead to strong real-space magnetic fluctuations, which may be responsible for the giant anomalous Hall effect observed in experiments. Current loops can be generated in the weak spin ferromagnetic states as a result of spin-orbit coupling effects, giving rise to spin-dependent orbital magnetic fluxes. The electronic structures in the ferromagnetic states are significantly reconstructed which have nearly compensated electron and hole carriers from two bands, which may explain the recent transport data at low temperatures. On the other hand, we do not find any diverging instability mode at $M$ point driven by electron-electron Coulomb interactions. The unstable phonon modes at $M$ and $L$ points indicate that the CDW state in the system is likely driven by phonon instability and electron-phonon couplings. Our results imply separate magnetic and charge instabilities in the AV₃Sb₅ (A=K, Rb, Cs) system, and the system may undergo an additional magnetic phase transition marked by the onset of anomalous Hall effect and reconstruction of Fermi surfaces at a critical temperature lower than the CDW transition temperature. Our work sheds light on the puzzling magnetic and charge properties of the V-based kagome metal systems, and will provide useful guidelines for future studies.

We thank Yueshen Wu, Jun Li, Kun Jiang, and Zheng Han for valuable discussions. This work is supported by the National Science Foundation of China (grant no. 12174257), the National Key R & D program of China (grant no. 2020YFA0309601), and the start-up grant of ShanghaiTech University.