Atomistic Origin of Diverse Charge Density Wave States in CsV₃Sb₅

The kagome lattice provides a fertile playground to explore intriguing correlated phenomena [1, 2, 3, 4]. Recently, the vanadium-based kagome metals $A$V₃Sb₅ ($A$ = K, Rb, or Cs) [5] was reported to feature charge density wave (CDW), as well as superconductivity and Z₂ topology [6, 7]. The interplay between CDW and other properties triggered enormous research interests [8, 9, 10].

As an example of $A$V₃Sb₅, CsV₃Sb₅ exhibits diverse CDW instabilities, enriching their couplings with different electronic orders [6, 7, 9, 8, 11, 12, 10, 13, 14, 15]. The high temperature pristine phase of CsV₃Sb₅ crystallizes in the layered structure of V-Sb sheets intercalated by Cs atoms with the space group of P6/mmm [5], as depicted in Fig.1(a, b). In the kagome layer, the V atoms form a kagome net interspersed by Sb1 atoms and sandwiched by two honeycomb Sb2 layers. CsV₃Sb₅ exhibits a CDW transition at T_CO $\approx$ 94 K, corresponding to a 2$\times$2$\times$1 inverse star of David (ISD) pattern with $D_{6h}$ point group [6, 16]. At an intermediate temperature of T_SO $\approx$ 50-60 K, a unidirectional 4a₀ or 5a₀ stripe order has been observed on Sb-terminated surface, the atomic pattern of which is still elusive [8, 17, 11, 18, 19, 20]. With a slightly lower temperature of T_N $\approx$ 40 K, a nematic state occurs, reducing the rotation symmetry from C₆ to C₂ [11]. Also, the 2$\times$2$\times$2 and 2$\times$2$\times$4 CDW phases have been observed, associated with the stacking of the adjacent kagome layers [22, 21]. Moreover, such CDW instabilities compete with each other and the ground state varies with doping, surfacing and pressure [23, 24, 25, 26]. Among different CDW states, the origin of ISD state was successfully explained with soft phonons in a collective manner [16]. However, it still lacks an overall understanding of the known and possible hidden CDW instabilities in CsV₃Sb₅, especially in an atomistic perspective.

In this work, we apply the first-principles-based effective Hamiltonian approach [27, 28], which is widely used in ferroelectric systems [29, 30], to investigate the atomic CDW instabilities in CsV₃Sb₅. Significantly, a series of breathing-kagome-like D_3h phases are predicted with mild tensile epitaxial strain. The atomic origins for the new D_3h phases and the known ISD state are comprehensively understood in terms of V-V and V-Sb pairs. Particularly, the atomic Hamiltonian reveals the existence of ionic Dzyaloshinskii–Moriya interaction (DMI), which is found to be the driving force of the formation of these nonlinear CDW patterns.

CDW transition from DFT. Phonon dispersion of the pristine phase at zero strain is displayed in Fig. 1(f). The imaginary phonon mode at M point [e.g., $2\pi$(1/2, 0, 0)] indicates a 2$\times$1$\times$1 modulation of the structure, which mainly involves the V-V in-plane movements as shown in Fig. 1(d). A linear combination of the soft modes at the three symmetry-equivalent M points give rise to the star of David (SD) or ISD [Fig. 2(a)] pattern within a 2$\times$2$\times$1 supercell. Note that the ISD structure exhibits inverse atomic displacements compared to the SD structure. The ISD phase is energetically favored (see Fig. S1 in Supplemental materials (SM) [27]), which is consistent with previous results [16] and will be explained in terms of the microscopic interactions discussed below. Similarly, the three L-point soft modes indicate a 2$\times$2$\times$2 structural reconstruction, which is manifested by the ISD + ISD stacking pattern with $\pi$-phase shift (see Fig. S2 [27]), agreeing well with previous studies [16, 17]. For simplicity, we will focus on the 2$\times$2$\times$1 CDW orders in the following discussion unless stated otherwise.

Let us then look at the strain effects. Recent experiments report different critical temperatures between bulk and nano flakes of CsV₃Sb₅, indicating surface strain releasing plays an important role [31, 32, 33]. We thus investigate the phonons of the CsV₃Sb₅ pristine phase with in-plane biaxial strain (see Method in Part I of SM [27]). As shown in Fig. 1(g, h), the phonon instability changes rapidly with increasing tensile strain. An additional imaginary phonon mode appears at $\Gamma$ point at 2% strain and is found to correspond to a breathing kagome mode, in which the V atoms move inward in one triangle and outward in the adjacent triangles [Fig. 1(c)]. Superposing this mode onto the ISD structure leads to a new 2$\times$2$\times$1 structural reconstruction, characterized by shrinked triangle of V atoms with surrounding V-V dimers (D_3h-2 phase in Fig. 2(b), here, the $n\times n$ supercell with D_3h symmetry is denoted by D_3h-n, [see Part VIII-X, XIV of SM for band structures, scanning tunneling microscopy (STM), structural information and XRD patterns of these states [27]]. In addition, the lowest imaginary frequency moves from M point under zero strain to approximately the middle of $\Gamma$-M for 2.8$\%$ tensile strain, corresponding to the mode characterized by the V-V in-plane dimerization within a 4$\times$1$\times$1 supercell [Fig. 1(e)]. Such dimerization with a phase shift along ${\bm{a}}$ axis strongly implies a new possible 4a₀ stripe order [8, 17, 11, 20, 19], consistent with the observation of Ref. [8] (see Part IX for STM images [27]). Considering the multiplicity of above imaginary mode, the new D_3h-4 state [Fig. 2(c)] can be constructed. Interestingly, as strain goes more tensile, it is found that the lowest frequency moves more and more toward $\Gamma$ point, indicating larger and larger supercells. For example, the D_3h-6 state [Fig. 2(d)] can emerge at the strain of 3.5% (Fig. S4 [27]). Such results thus imply a continuum evolution of CDW states with increasing strain.

In order to determine the ground states at different strains, we perform density functional theory (DFT) calculations on various CDW configurations, which are condensated from the soft phonons (see Part V for all considered structures [27]). Figure 2(e) shows the energy-vs-strain curves of selected low-energy states. It indicates that, for the compressive strains exceeding -3.4$\%$, CsV₃Sb₅ prefers the pristine phase [Fig.1(a, b)], which is in line with the observation that pristine phase becomes ground state under hydrostatic pressure [20]. From -3.4$\%$ compressive strain to 0.7$\%$ tensile strain, the ISD state shown in Fig. 2(a) is energetically favorable, which is also consistent with observations [6]. Particularly, between tensile strain of 0.7$\%$ and 2.7%, the aforementioned new D_3h-2 CDW state is stabilized (see Fig. S8 for phonon dispersion and XIII for the prediction of CsV₃Bi₅ which exhibits D_3h-2 phase in the absence of any strain [27]). As the tensile strain increases, the V-V dimers form larger and larger triangle, which encircles more and more shrinked triangles of V atoms, corresponding to the D_3h-4 phase [Fig. 2(c)], D_3h-6 phase [Fig. 2(d)], etc. It is thus reasonable to speculate that, for even stronger tensile strain, the ground state will transform to a phase with only shrinked triangles, which is actually the breathing kagome state. Note that the D_3h-n phases also prefer that adjacent kagome layers have a $\pi$-shift. Such $\pi$-shift further lowers the rotation symmetry from C₃ to C₂ (Fig. S3 [27]).

Atomic driving forces from Effective Hamiltonian. In order to get more insights into the microscopic origins of the CDW states, we develop first-principles-based effective Hamiltonians for CsV₃Sb₅ (see Table. S5 for KV₃Sb₅ and RbV₃Sb₅). Here we adopt the symmetry-adapted cluster expansion method, as implemented in the PASP software [28]. Such method enables defining local modes ${\bm{u}}$ (atomic displacements) and generating all possible forms of interactions, i.e., the invariants. With input energies of random structures from DFT, the coefficients of such invariants can be fitted using machine learning method for constructing Hamiltonian [34]. Note that the effect of electron-phonon coupling is implicitly included in these coefficients. For the present case of CsV₃Sb₅, local modes are defined as (i) in-plane displacements of V atoms and (ii) out-of-plane displacements of Sb2 atoms, while other minor movements are neglected. The initial model contains enough invariants, including the self-energy of local mode up to fourth-order and the cluster interactions up to the four-body and fourth-order. After repeated fitting and refining, it arrives at below Hamiltonian,

where $\varepsilon^{\rm self}$ is the single site energy and $\varepsilon^{\rm pair}$ is the energy of pair interaction, which indicates that three-body and four-body interactions are not important here. Herein the two effective Hamiltonians are constructed separately for strains of 0% and 2% (see more strains in Table. S3), with the primary coefficients listed in Table 1, while a full model and details can be found in Part VII in SM [27]. The energies from DFT and effective Hamiltonians are in good consistency (see Fig. S11). Note that the forms of interactions are expressed in local basis [see Fig. 3(a-c)] for simplicity. As shown in Fig. 3(e), Monte Carlo simulations with the obtained effective Hamiltonians indicate that (i) at zero strain, ISD state forms below 100 K, which is rather close to the measured T${}_{\textup{CO}}$ = 94 K [6]; and (ii) at 2% strain, the D_3h-2 state emerges below T${}_{\textup{D}}$ = 20 K, indicating that the tensile strain largely lowers the CDW critical temperature.

We now work on understanding the formation of ISD and D_3h-2 states by analysing the invariants and their coefficients. According to the experience of Hamiltonians for ferroelectrics, the second-order interactions usually act as the driving force of the phase transition, while that the fourth-order interactions prevent too large atomic displacements [29]. Such mechanisms lead to the well-known double-well energy landscape. Here, as shown in Table 1, single site terms, either second or fourth orders, all have positive coefficients, indicating that neither ISD nor D_3h-2 states are initiated by on-site instability. On the other hand, three second-order pair interactions are found to drive the deformations to CDW states (note that the local basis {XYZ} is adopted): (i) the first-nearest neighbor (1NN) V-Sb2 interaction, with the form of $u_{1Y}u_{2Z}$ and a positive coefficient, drives V atoms toward in-plane trimerization (shrinked triangle) and Sb2 atoms moving toward out-of-plane direction [Fig. 3(a)]; (ii) the 1NN V-V interaction, with the form of $u_{1Y}u_{2Y}$ and also a positive coefficient, also tends to form trimers of V atoms [Fig. 3(b)]; and (iii) the second-nearest neighbor (2NN) V-V interaction, with the form of $u_{1Y}u_{2Y}$ and a negative coefficient, drives trimerization of V atoms at a large distance [Fig. 3(c)] (the sign of the coefficient of these pair interactions can be related to the band splitting near Fermi level, which may lower the energy relative to the band crossing at the inverse sign, see Figs. S13-14 [27]). Note that (i) and (ii) lead to trimers with the same orientation, while (iii) results in trimers with opposite orientations. Such analyses indicate that (i) and (ii) tend to stabilize breathing kagome state, while (iii) favors ISD state. Interestingly, the 2NN V-V pair has a third-order interaction of u${}_{1Y}^{2}$u_2Y+ u_1Yu${}_{2Y}^{2}$, which naturally lifts the energy degeneracy between SD and ISD, or between D_3h-n and inverse D_3h-n. We then look at the strain effects. At zero strain, the 2NN V-V pairs play the dominant role, which results in the ISD deformation. When a 2$\%$ tensile strain is applied, values of $J_{\textup{V-Sb2,2}}$ and 1NN $J_{\textup{V-V,2}}$ are significantly enhanced (Table 1), which tends to stabilize the breathing kagome pattern. This is understandable since the larger lengths of 1NN V-V pairs under tensile strain enhance the tendency close to each other. Meanwhile, the 2NN $J^{\prime}_{\textup{V-V,2}}$ barely changes with strain and thus still favors the formation of ISD phase. The interplay of such two strain effects is shown in Fig. 3(d), which clearly demonstrate that the D_3h-2 reconstruction (yellow bonds) originates from superposing breathing kagome (green bonds) and the ISD (red bonds) deformations (see Part IV for the superposition of D_3h-4 and D_3h-6 states).

Let us now further understand the nature of V-V interactions in the global basis. As shown in Fig. 4, between 1NN V-V pairs, e.g., the pair of V1-V2, our Hamiltonians indicate an interaction in the form of U_1xU_2y-U_1yU_2x, which can be written as ${\textbf{D}}_{12}$$\cdot$(U₁$\times$U₂), where ${\textbf{D}}_{12}$ has only $z$ component and the displacements of ${\textbf{U}}_{1}$ and ${\textbf{U}}_{2}$ possess only in-plane components. Such interaction is analogous to the magnetic DMI and is thus named as ionic DMI (i-DMI) here. As U₁$\times$U₂ has the same symmetry transformation properties as S₁$\times$S₂ (where S represents spins), the direction of the DM vector in the i-DMI case is determined by the crystal structure in the same way as the magnetic DMI case. According to Moriya’s rules [35], DM vector of V1-V2 pair is along $\pm z$ direction, as a result of an in-plane mirror including V1-V2 pair and a mirror perpendicular to V1-V2 pair. Such i-DMI favors related displacements being perpendicular to each other with the $xy$ plane, which is similar to the effects of magnetic DMI. However, it finally forms the 120^∘ patterns, as shown in Fig. 4, due to the competition associated with kagome lattice (see Part XI of SM [27]). Note that the i-DMI of 2NN V-V pairs is not discussed here due to their negligible values.

From the aspect of i-DMI, we now further understand the CDW transition from ISD to D_3h state. At zero strain, the i-DMIs exhibit a positive value of $D^{z}_{12}=158$ meV, which stabilizes a ISD-like state, as shown in Fig. 4a, with an energy gain of $E_{a}=-\sqrt{3}D^{z}_{12}|{\textbf{U}}_{1}|^{2}$ per V atom. For the exact ISD pattern that CsV₃Sb₅ displays, such i-DMI leads to an energy gain of $-\frac{1}{2}D^{z}_{12}|{\textbf{U}}_{1}|^{2}$, while the pattern shown in Fig. 4b results in an energy cost of $-E_{a}$ per V atom. Note that it is assumed that all V atoms share the same amount of displacements. In contrast, with a tensile strain of 2%, the i-DMIs change their signs and render a negative value of $D^{z}_{12}=-39$ meV, which favors the formation of breathing kagome pattern, as shown in Fig. 4b. In such case of negative $D^{z}_{12}$, the breathing kagome pattern yields an energy gain of $E_{b}=\sqrt{3}D^{z}_{12}|{\textbf{U}}_{1}|^{2}$ per V atom, while the ISD-like pattern cost energy of $-E_{b}$ per V atom. It is thus clear that, above a critical point of tensile strain, the i-DM vectors flip their directions, which deforms the ISD phase and stabilizes the $D_{3h}$ breathing kagome state.

The present discovery of i-DMI as a driving force to induce CDW states provides a new perspective to understand structural phase transitions. We note that the existence of i-DMI was first pointed out in the low energy distorted perovskites [36]. Here, we find that CsV₃Sb₅ is the first system where the i-DMI exists in the high-symmetry high temperature structure, and is responsible for the structure phase transition to the low energy structure. Considering that the kagome lattice intrinsically lacks inversion center between first nearest neighbors, the i-DMI is likely to be a common mechanism to induce structural distortions there.

In summary, combining DFT calculations, we apply the effective Hamiltonian approach to CsV₃Sb₅ and predict a series of D_3h-n CDW phases and a possible 4a₀ pattern, besides the known ISD structure. The nature of such CDW states is revealed as that the 2NN V-V pairs dominate the formation of ISD state, while that the 1NN V-V and V-Sb2 interactions tends to stabilize the D_3h-n breathing kagome state. The ionic DMI is identified as a new driving force to induce structural distortions, which is expected to commonly exist in different kagome CDW systems. Experimental verification of the new atomics mechanisms and new CDW states revealed in this work is called for.