Condensation of preformed charge density waves in kagome metals

Charge density wave (CDW) is a spontaneous spatial modulation of charges in solids whose general microscopic descriptions are yet to be completed. Kagome metals of AV3Sb5 (A = K, Rb, Cs) provide a chance to realize CDW intertwined with dimensional effects as well as their special lattice. Here, based on a state-of-the-art molecular dynamics simulation, we propose that their phase transition to CDW is a condensation process of incoherently preformed charge orders. Owing to unavoidable degeneracy in stacking charge orders, phases of preformed orders on each layer are shown to fluctuate between a limited number of states with quite slower frequencies than typical phonon vibrations until reaching their freezing temperature. As the size of interfacial alkali atom increases, the fluctuations are shown to counterbalance the condensation of orderings, resulting in a maximized transition temperature for RbV3Sb5. Our results resolve controversial observations on their CDWs, highlighting a crucial role of their interlayer interactions.


INTRODUCTION
The kagome lattice has been considered as a fertile ground to realize exotic quantum phases because of its high degree of geometric frustration and nontrivial band structures [1][2][3].Particularly, when the number of electrons per site is 2/3±1/6, the Fermi energy (E F ) is at the van Hove singularity (vHS) such that the density of states diverges to induce various broken-symmetry phases [1][2][3].The recent discovery of kagome metals of KV 3 Sb 5 , RbV 3 Sb 5 and CsV 3 Sb 5 (abbreviated as AV 3 Sb 5 where A denotes alkali atoms of K, Rb and Cs) [4] bring a chance to realize those phases because their E F 's are close to vHS [4,5].Indeed, several studies so far have demonstrated noteworthy states such as charge density waves (CDWs) [4,[6][7][8][9], superconductivity [10], loop currents states [11][12][13] and electronic nematicity [14], to name a few.
On the other hand, a recent progress in manipulating stacked two-dimensional (2D) crystals opens a new window to study phase transitions [15][16][17].Owing to their anisotropic lattice structures and interactions, distinct electronic phases are realized in the same material by pressure [18], a number of layers [17,19] and dopings [20] etc.Unlike simple lattice structures of hexagonal or rectangular shapes in typical layered materials [16], newly discovered kagome metals of AV 3 Sb 5 add another complexity of frustrating intralayer orders in studying phase transitions, thus providing a unique chance to study their interplay with extreme anisotropic interactions.
Owing to the vHS in kagome metals, CDWs occur with decreasing temperature (T ) from normal metallic phase [21,22].In spite of their well-established presence, a few enigmatic observations raise questions on their formation mechanism.One of anomalous behaviors is the absence of phonon softening in inelastic X-ray [8], neutron [21] and Raman scattering experiments [23] while these results are incompatible with unstable phonon modes at M -and L-point of the Brillouin zone (BZ) in recent ab initio calculations [5,24,25].Moreover, from the fact that the formation energy of CDW (E CDW ) [5] monotonically decreases with the size of the alkali atoms, one could naively expect that the transition temperature for CDW (T CDW ) increases accordingly.However, as the size of alkali atom increases from K to Rb and to Cs, the observed T CDW first increases from 78 K to 102 K and then decreases to 94 K [2], respectively, contradicting the conventional behaviors [26,27].
In this work, we theoretically showed that kagome metals first develop lattice distortions within each layer at much higher temperature of T * than T CDW .Then, their phases are not fixed simultaneously with their amplitudes but fluctuate, which is shown to be inevitable owing to multiple ways of stacking CDWs between adjacent layers with the same energetic cost.As T decreases further, the phases of preformed orders eventually stop varying at T CDW .The characteristic timescale for the fluctuation is at least 10 5 slower than that a typical vibrational frequency of vanadium lattice.For CsV 3 Sb 5 , it is shown that a phonon softening signal exists only around L-point in a narrow temperature range around T * .Between T * and T CDW , the thermodynamic behaviors of preformed CDWs can be described by 4-states Potts model which undergoes a first-order phase transition at T CDW [28].Our quantitative analysis demonstrates that the fluctuations increase as the size of alkali metal increases, thus competing with E CDW for preformed orders.So, resulting T CDW maximizes not with Cs but with Rb atom and agrees well with the experimental trend [2], highlighting roles of weak interlayer interactions in determining phase transitions of layered kagome metals.
Our computations for two markedly disparate orders are made possible by developing a new large-scale molecular dynamics (MD) simulations method.A new method is specifically designed to describe lattice distortions using the relative displacements between atoms as basic variables (See Methods).Polynomial interatomic potentials are constructed using the linear regression for fitting about 2,500 basis functions to our first-principles calculation results.So, our method is accurate enough to reproduce our ab initio results very well, thus enabling us to perform large scale accurate lattice dynamics in a very long timescale.

2D charge orders and their interlayer interactions
Figure 1a shows the crystal structure of AV 3 Sb 5 .Vanadium atoms form 2D kagome lattice at vertices of dashed lines in Figs.1b and c.One thirds of Sb atoms occupy the centers of dashed hexagons and the rest two thirds are at the upper and lower plane of the kagome lattice, respectively.Alkali atoms play as spacer atoms between the layers and determine the interlayer interactions.As shown in Fig. 1c, the CDW phase has been known to form the inverse star-of-David (iSOD) structure where V atoms in 2 × 2 units are modulated into one hexagon and two triangles instead of the star-of-David structure constructed by inverting the atomic displacements of the iSOD [5,25].Our ab initio calculations also prefer the former over the latter, agreeing with previous studies [24,29].We note that, without interlayer interactions, the 2 × 2 CDW in each layer can take any of the four translationally equivalent structures or phases indicated by four rhombi in Fig. 1c.The interlayer interaction, however, lifts the degeneracy between those random stacking phases.In Fig. 1d, we display possible stacking orders where the different phases are denoted by different colors corresponding to those in Fig. 1c.
We find that the same phases for neighboring layers are hardly realizable owing to the large energetic cost while the different phases are allowed and their energies for couplings are equivalent (See Supplementary Information (SI) Section 1).This implies that second-nearestneighboring interlayer interactions should play an important role for stacking order.Without it, long-range stacking orders are absent and random stacking of CDWs would be dominant as long as the neighboring layers avoid to be the same phase (SI.Sec. 2), which contradicts experiments [8,10,30].Our first-principles calculation estimates that its magnitude is order of 1 meV for 2 × 2 units in Fig. 1c.Such a small interaction may be related with possible CDW stacking faults or variations in interlayer ordering periods [8,10,23,30,31].Notwithstanding its small magnitude, as we will demonstrate hereafter, it is an important interaction in freezing fluctuating charge orders of kagome metals.
From our MD simulations, we find that a interlayer correlation at temperature far above T CDW is not negligible in spite of the small negative long-rage interaction.Specifically, for CsV 3 Sb 5 at T = 130 K, the probability that the pairs of second-nearest-neighboring layers have the same phases in the thermal ensemble turns out to be 0.66, revealing its crucial role (without it, the probability should be 0.25).From these considerations, the interlayer interactions are shown to have ambivalent characteristics depending on the interaction ranges such that nearest-neighboring layer should avoid the same phase while second-nearest-neighboring layer favors the same one.
Aforementioned formation of intralayer CDWs are captured explicitly by collecting the MD trajectories of V atom density of ρ V (r; T ) for a given T .In Fig. 2a, it is shown that ρ V of CsV 3 Sb 5 at T = 200 K with a simulation time of 0.4 nanoseconds perfectly match with ideal kagome lattice points.In enlarged views around the lattice point in Fig. 2b, it starts to deviate from the lattice points at T ≃ 160 K.As the temperature decreases further, it continuously deviates from the lattice points and the deviation reaches ∼0.1 Å at T = 20 K.

Phonon instability
The structural instability is also reflected in the temperature-dependent phonon dispersions.At finite T , the phonon spectrum can be extracted from the density-density correlation function of S ρρ (k, ω; T ) where k = (k x , k y , k z ) is a crystal momentum and ω is a frequency of phonon (See Methods for details).In Fig. 2c, we plotted S ρρ (k, ω; T ) of CsV 3 Sb 5 along the high symmetric lines on (k x , k y ) plane with k z = b z /2 where b z is the out-of-plane reciprocal lattice vector.At T = 250 K, a soft mode denoted by an arrow is clearly shown at L-point.As T decreases, the frequency of the soft mode approaches zero.Eventually, at T * ≃ 160 K, the structural instability develops, that may indicate a phase transitions with lattice distortions [32], compatible with our simulation of ρ V (r; T ) in Fig. 2b.However, the structural instability here does not guarantee phase transition.Even though the local 2×2 CDW orders within each layer and their π-phase shift between adjacent layers can fully develops, the local CDW orders change their phases dynamically and do not spontaneously break any symmetry of the crystal.As we will see later, a global broken symmetry state occurs at much lower temperature of T CDW so that we could interpret T * as the preformation temperature of CDW.
We note that right below T * , all the softening signatures near L-point disappear as shown in Fig. 2c.At T CDW < T < T * , the preformed CDW domains have orders-of-magnitude longer lifetimes than vibrational periods as will be shown later.So, the effective phonon frequencies is essentially from the curvature of the potential around one of the four degenerate CDW phases.The phonon mode at L-point corresponding to this new metastable atomic configurations becomes harder as T decreases such that near T CDW , no softening signature can be found.
We also find that no soft modes exist along ΓM (Fig. S2a in SI.Sec.3), which is consistent with recent experiments reporting unexpected stability of phonon modes along ΓM at the transition [8,21,23].Even at T > T * , the lowest M -phonon remains hard unlike L-phonon, ex- cluding a possible preformation of CDW.At this high temperature, anharmonicity-induced phonon hardening occurs for both phonons but the effect is much stronger for M -phonon due to the enhanced interlayer force constants.If not considering the anharmonic effect, both M -and L-phonon are expected to be unstable for the structure as were demonstrated by recent first principles calculations [5,24,25].We can expect that the hardening effect for phonon at U -point of BZ (2 × 2 × 4 CDW) is in-between M -and L-phonon.Indeed, the softening of dispersions at U -point around T * is less conspicuous than one at L-point (Fig. S2b in SI.Sec.3), not fully developing into instability.

Fluctuation of preformed orders
From the results so far, it seems that CsV 3 Sb 5 undergoes a typical phase transition at T * .However, owing to large degeneracy for pairing charge orders between layers, a slow fluctuation between them proliferates across the whole layers.We note that our simulation time of 0.4 nanoseconds for the apparent ordering is long enough for simulating usual structural transitions [33,34].For kagome metals, it turns out to be not.With ordersof-magnitude longer simulations time of 12 nanoseconds, we indeed found a clear signature for the fluctuation.We also confirm that the fluctuation is not an artifact from the size effect of simulation supercell (See detailed discussions in SI.Sec. 4).
The phase fluctuations are unambiguously identified by eigenmode decomposition analysis.As shown in Fig. 1b, the in-plane CDW decomposes into three symmetrically equivalent eigenmodes of − → M i (i = 1, 2, 3).It can be easily checked that their linear combinations of m 1 (t) reproduce the four phases of 2 × 2 iSOD structures as shown in Fig. 1c.Here m i (t) (i = 1, 2, 3) is a time-dependent coefficient for eigenmodes of − → M i where t denotes time.By inspecting temperature-dependent time evolution of m i (t) for a specific 2 × 2 supercell, we can identify slow CDW fluctuations in kagome metals (see further details in SI.Sec. 5).
Figure 3a shows a temporal evolution of m 1 (t) for T = 140 K.In addition to picoseconds phonon vibrations, there are clear sign changes of m 1 (t) in timescale of a few nanoseconds without altering its averaged absolute magnitude.This indeed indicates phase flips of CDWs.The rate of phase flips can be quantified by computing a Pearson correlation coefficient (PCC) [35] for time interval of ∆t, r 11 (∆t , where E[m 1 ] and σ[m 1 ] are time average and standard deviation of m 1 , respectively.As shown in Fig. 3b, the fast phonon vibrations and the slow phase flip are manifested as picoseconds oscillation and the subsequent slow decay, respectively (For comprehensive analysis of PCCs, See SI.Sec. 5).The decay rates of r 11 in Fig. 3b imply that the phase-coherence time does not exceed a few nanoseconds at 140 K and becomes longer as T decreases further.With inclusion of the phase flips, for T CDW < T < T * , the averaged local density profile of vanadium atoms of ρ V should change from a Gaussian thermal ellipsoid shown in Fig. 2b into a four-peaked or non-gaussian distribution as shown in Fig. 3c.As T approaches T CDW , each peak becomes to be more distanced and eventually, the fluctuations are quenched to one of the four peaks at T = T CDW .This implies that in experiments incapable of resolving nanosecond dynamics, only the averaged properties will be observed without any symmetry-breaking feature between T CDW and T * .However, those phase fluctuations are reflected as zero frequency peak at S ρρ (k, ω; T ) (see Fig. S3 in SI.Sec.3), that may be accessible by inelastic neutron scattering.

Critical temperatures of charge orders
Our simulations hitherto reveal that in kagome metals, T CDW is nothing but a critical temperature for a globally ordered phase by condensing the CDWs that are preformed at T * (> T CDW ).Although the physical process of the phase transition is now understood, a reliable quantitative MD calculation of T CDW is undoable because the time for quenching preformed orders well exceeds our attainable simulation time.This motivates us to map our systems into an anisotropic 4-states Potts model [28] to compute T CDW quantitatively using a well-established statistical method [36].Due to the two well-separated time scales of dynamics (thermal vibration and phase fluctuation in Fig. 3a), if we average the atomic trajectories of preformed CDW's over 100 ps, the averaged snapshot will look like one of four 2×2 CDW phase in Fig. 1c, but with temperature-dependent amplitudes.The four phases become 4-states 'spin' variables on lattice points of layered triangular lattices as shown in Fig. 4a.So, an effective Hamiltonian for the interacting spins therein can be written as where s i,α is an effective four-states spin at i-th site of αth layer, J ∥ their effective intralayer nearest-neighboring (n.n.) interaction, J ⊥ interlayer n.n. and J ⊥2 every other interlayer n.n.interactions, respectively as described in Fig. 4a.Here, δ(s i,α , s j,β ) is zero if states of two effective spins of s i,α and s j,β are different, otherwise, it is one.We note from considerations above that signs and magnitudes of the interactions should be Unlike typical Potts models [28], interaction parameters in Eq. 1 are explicit functions of temperature because amplitudes of CDW vary as the temperature changes.For T = 0, J ∥ (0) and J ⊥ (0) can be readily calculated from the energy cost forming domain wall and energy differences between different stacking, respectively (See detailed procedures in SI.Sec. 6).For T > 0, they are similarly obtained from thermally averaged potential energy of domain structure.Here, the thermal average is approximately calculated using harmonic and rigid approximations for density correlation functions (See Methods).Then, we have temperature-dependent exchange parameters of J ∥ (T ), J ⊥ (T ), J ⊥2 (T ) in our 4-states Potts model.So, to compute T CDW , we need to solve the model self-consistently because temperature and the exchange parameters depend on each other.For kagome metals with alkali atoms of K and Rb, we can bypass timeconsuming reference calculations by rescaling polynomial interatomic potential obtained for CsV 3 Sb 5 without degrading accuracy (See SI.Sec. 7).We also note that the accuracy of J ⊥2 from our ab initio calculation of CsV 3 Sb 5 seems to be marginal.While the absolute magnitude of J ⊥2 determines the ground state stacking structure, it turns out to play a minor role in determining T CDW .So, we set it as a parameter of J ⊥2 = −0.5J⊥ and checked that T CDW is lowered by at most 7 K when we decrease J ⊥2 to be 0.
In Figs.4b and c, the self-consistently obtained J ∥ (T ) and J ⊥ (T ) are shown for three kagome metals.The critical temperature for preformed orders of T * can be determined from a condition that the amplitude of intralayer CDW vanishes or that J ∥ (T ) and J ⊥ (T ) become zero simultaneously.As expected, the intralayer coupling of J ∥ (T ) is largest (smallest) for CsV 3 Sb 5 (KV 3 Sb 5 ), being consistent with the lowest (largest) E CDW [5].Accordingly, T * monotonically increases as the size of alkali atom increases as shown in Fig. 4d.
The J ⊥ displays different temperature dependence compared with J ∥ (T ).Since all the compounds here share the same kagome lattice layer formed by V and Sb atoms and differ by their alkali atoms, J ⊥ will be dominantly set by size of the alkali atoms.Indeed, the interlayer distance for Cs compound is longest among them (SI.Sec. 1) so that it has the smallest J ⊥ (0) as shown in Fig. 4c.This implies that the weakest J ⊥2 of Cs compound will invoke the most severe fluctuation of CDW phase among three kagome metals.We also note that as temperature approaches T * , the magnitude of J ∥ is comparable to those of J ⊥ as well as J ⊥2 .Therefore, we identify two competing features to determine thermodynamic states, i.e., the phase fluctuation of CDW counterbalances the E CDW .So, even though CsV 3 Sb 5 has the highest T * , its T CDW can be lower than the other.Estimated T * and T CDW from Eq. 1 indeed confirm such a competitive interplay between inter-and intralayer orders, agreeing well with experimental trend as shown in Fig. 4d.

DISCUSSION
The preformed charge order discussed so far is not measurable as discontinuities in thermodynamics variables such as a specific heat or susceptibility.Also, we expect that Raman scattering [23,24,37,38] or nuclear magnetic resonance measurements [39,40] may not be so easy to capture its explicit signals owing to its fluctuating phases and dynamical nature.But there might be at least three experimental signatures already indicating its existence.First, a recent X-ray diffraction experiment [22] reports that the integrated peak intensity of 2 × 2 × 2 CDW order in CsV 3 Sb 5 survives up to 160 K, well above T CDW = 94 K.In addition to that, the thermal expansion of the in-plane lattice constant shows a slope change twice at T CDW and at T = 160 K, respectively [22], implying the qualitative structural changes at higher T than T CDW .Second, the coexistence of the ordered and disordered phases near T CDW is measured through nuclear magnetic resonance measurements [39,40].This is consistent with a first order phase transition from our anisotropic 4-states Potts model in Eq. 1. Third one is the absence of phonon softening around T CDW [8,21,23].As we already discussed, this absence is not anomaly but a consequence of preformation of CDW at much higher T * .
In addition to these indirect evidences, we may have a way to detect fluctuating charge orders directly.Because the preformed CDW not only changes the vibrational properties but also introduces slow cluster dynamics, we expect that it can be also detected as a central peak in the dynamic form factor of density fluctuation, or apparent symmetry-forbidden signals in Raman scattering reminiscing the precursor formation in ferroelectric materials [41,42].
We note that our current MD methods from firstprinciples approaches do not show any state related with broken time-reversal-symmetry states.Within our methods, the vanadium d-orbitals have considerable out-ofplane band dispersion across the Fermi energy [43] so that spin-spilt states related with them hardly develop.From these, our MD could not touch upon the various experimental signatures [1][2][3][11][12][13] on time-reversalsymmetry broken phases.However, if exotic electronic and magnetic orders can couple to bond orders below T CDW , we expect that our method could reflect the corresponding structural evolutions at lower temperature.
Lastly, we point out that the tiny magnitude of J ⊥2 is essential in controlling the phase fluctuation with fixed CDW amplitudes across the kagome layers.This strongly suggests that a facile engineering of thermodynamic states of kagome metals could be possible through external controls of interlayer interactions [44,45] and that our current understanding of asynchronism in thermodynamics transitions of kagome metals can be applied to understand various subtle phase transitions in layered 2D crystals [15][16][17][18][19][20]46].

Construction of interatomic potential
Our interatomic potential is based on the linear regression with a set of polynomial basis functions which is adequate for small displacement phonon calculations [47].We note that a similar method has been used for dynamical properties of kagome metal recently [48].For the interatomic potential to meet a few physical conditions such as (1) continuous translational symmetry, (2) point and space group symmetry of the crystal, (3) permutation of variables, and (4) asymptotic stability for a large displacement, basis functions are constructed as follows.
First, to meet the condition (1), we replaced the variable from the atomic displacement u αi to the relative displacement u αi − u α ′ i where α i is a condensed index for Cartesian components, basis atom in a unitcell and the Bravais lattice point.The index of α ′ i shares the same Cartesian component with α i but has different basis atom index and Bravais lattice point, respectively.Then, the interatomic potential of V can be written as (2) where the brace indicates that the summation runs for all indices of α i and α ′ i in the summands.The condition (2) can be straightforwardly incorporated by applying all symmetry operations R of a given crystal to V n in Eq. 2, that can be written as where N R is a number of symmetry operations.The condition (3) will be automatically met from linear dependence check of basis function.The condition (4) is usually ignored in phonon calculations but becomes to be critical for potentials with a huge number of variables.Since our interatomic potential has thousands of variables, without explicit consideration of the condition (4), optimization processes do not end but diverge.In general, it is hard to prove that some multivariable polynomial is bounded below but we avoid the negative divergence by forcing the highest-order basis function to be multiples of squares and by constraining their coefficients to be nonnegative during the linear regression process.Rigorously, though this approach does not guarantee the condition (4), we find that the V in Eq. 2 constructed in this way did not cause any practical problem.The resulting basis functions are sorted in ascending order by the largest distance between the two atoms in the basis function and linearlydependent basis functions are removed by Gram-Schmidt orthogonalization.
For CsV 3 Sb 5 , we truncated the polynomial out at the 4 th order and cutoff the relative distance at 10 Å for the 2 nd -order basis functions and 6 Å for 3 rd -and 4 th -order ones.Importantly, only V displacements are used for the variables in the 3 rd -and 4 th -order basis function because we found that the anharmonic effect is prominent only between the V atoms due to their large displacements in CDW phases.The number of basis functions made in this way is 434, 521, and 1453 for 2 nd -, 3 rd -and 4 th -order, respectively.
The reference configurations are collected from multiple sources to ensure reliable sampling of potential energy surfaces.We performed first-principles MD simulations of 4 × 4 × 2 supercell at 100 and 200 K for 1600 steps with a time step of 10 fs.We then collected 53 configurations from each temperature that were equally spaced in time of 300 fs.Additional 30 configurations are collected from randomly displacing (< 0.04 Å) atomic positions from the optimized 2 × 2 × 2 CDW structure.Lastly, the linearly interpolated configurations between high-symmetric structures without CDW and the optimized 2 × 2 × 2 CDW structures are also used to collect 9 more configurations.These diverse sets of 145 configurations are found to be enough to reproduce energy landscape of low-energy stacking structures fruitfully.
The coefficients of basis functions are fitted to the calculated atomic forces.We first perform the linear regression procedure to minimize the following loss function, where C is the m-dimensional coefficient vector of basis functions of c α1α ′ n in Eq. 2, F i is the 3Ndimensional vector of reference atomic forces in the i th configuration including N atoms, and A i is the 3N × m matrix whose j th column is atomic forces in the configuration i calculated from the interatomic potential by fixing the coefficient of j th basis function as c j = 1 and by forcing all the other coefficients to be zeros.Here, m = m 2 + m 3 + m 4 where m k (k = 2, 3, 4) is the number of k th -order basis functions.
For the coefficients of 4 th -order basis function to be nonnegative, the coefficient vector of C obtained from Eq. 3 is again optimized using a nonlinear conjugate gradient method for the modified loss function, where D is a m × m diagonal matrix with conditions of D ii = −1 only when i > m 2 + m 3 and c i < 0 and otherwise, D ii = 1.This enforces nonnegative condition for the coefficients of 4 th -order basis function.

Molecular dynamics simulation
Using the obtained interatomic potential, our molecular dynamics (MD) simulations are performed with the velocity Verlet algorithm [49] for the integration of Newton's equation of motions and simple velocity rescalings are applied for every time steps of 5 fs to incorporate the temperature effect.We have used 60×60×12 supercell for calculations in the Fig. 2 and 12×12×12 supercell for Fig. 3.All thermodynamic ensembles are collected after 100 ps of thermalization steps.The density-density correlation function is given by S ρρ (k, ω; T ) = |ρ(k, ω; T )| 2 and ρ(k, ω; T ) is defined as where ρ(k, t; T ) = 1 √ N l,κ e −ik•R l,κ (t;T ) and R l,κ (t; T ) is the position vector of basis atom κ in a unitcell R l at the time t and temperature T .

First-principles calculations
Ab initio calculations based on density functional theory (DFT) are performed using the Vienna ab initio simulation package (VASP) [50] using a planewave basis set with a kinetic energy cutoff of 300 eV.The generalized gradient approximation with Perdew-Burke-Ernzerhof scheme [51] and dispersion correction using DFT-D3 [52] are adopted to approximate exchange-correlation functional and the projector augmented wave method [53] is used for the ionic potentials.For structural optimizations and first-principles MD, 12 × 12 × 12 and 3 × 3 × 2 k-points are sampled in the Brillouin zone of CsV 3 Sb 5 , respectively and the internal atomic coordinates of CDW states were optimized until the Hellmann-Feynman forces exerting on each atom becomes less than 0.01 eV/ Å.

Thermal average of potential energy
For anharmonic systems without light elements, ab initio MD simulation is a quite accurate method to investigate its thermal properties.Nevertheless, it is very timeconsuming so that less demanding computational methods have been developed.The lists of those methods can be found in the recent literatures [47,54].The main idea of those methods is that the original system can be approximated by the harmonic systems that minimizes the free energy of the system, and equilibrium structures and effective phonon frequencies (or effective potential) can be obtained from them.Our interest is, however, the dynamics between metastable configurations (i.e.CDW phase domains) so that the effective potential should be obtained by thermally averaging fast motions of atoms not only in the ground state configurations but also in the metastable ones.
To obtain such an average, we first assume that Nbody density correlation function of ρ does not vary when a ground state thermal ensemble experience a static spatial displacement, i.e., a rigid density approximation.Within this approximation, the static displacement corresponds to the metastable configuration and the thermal average can be performed using the exact ρ.Our second approximation is to replace the ρ with ρ H of an approximate harmonic system, which can be written as a closed form of normal modes [54].
Procedures to obtain the effective potential from the zero-temperature interatomic potential V can be formulated as follows.For periodic systems, V can be expanded with Taylor series of atomic displacements as where α i is a condensed index for Cartesian components, basis atom in a unitcell and the Bravais lattice point.
Here, ϕ α1•••αn is a n th -order force constant and u αi is a Cartesian component of a displacement vector from a reference position chosen to be a local minimum or saddle point of V.The brace indicates that the summation runs for all indices of α i in the summands.At finite T , the n th -order thermally averaged potential of V (n) T can be written as where the bracket denotes the ensemble average with ρ.Without varying ρ, the average potential energy for adding a static displacement of ⃗ U to the thermal ensemble of u αi can be written as where k th -order effective force constant is defined as 6, we replace ρ with ρ H for which the approximate harmonic system is chosen to have the similar density correlations with the original anharmonic system.Then, using a thermally averaged potential in Eq. 6, we compute a set of harmonic frequencies from which we reconstruct a new ρH .So, if a self-consistency for ρ H is fulfilled, we can obtain the desired thermal averaged potential of a rigidly shifted thermal ensemble.We use this method to compute the temperature dependent interaction parameters of J ∥ (T ) and J ⊥ (T ) as discussed in Sec.6

Size effect
To confirm that the sign changing events in Fig. 3a are not just a local fluctuation but a phase flip of whole layer, the events between 7 and 8 nanosecond are presented in Fig. S4.Each panel is obtained from the snapshot of molecular dynamics (MD) simulations taken every 50 femtoseconds.Here, 2 × 2 unit cells and the amplitudes m 1 for the cells in the panel are presented with the circles and their graded colors, respectively.The phase flip in a  M (i = 1, 2, 3) at M points.Blues dots are kagome lattice points and arrows denotes atomic displacements.Linear combinations of three q i M or m1q 1 M + m2q 2 M + m3q 3  M denoted as [m1, m2, m3] form (b) inverse star of David (iSOD) pattern or (c) star of David (SOD) pattern.We note that when a product of three coefficients, m1m2m3, is positive (negative), the resulting configuration becomes iSOD (SOD).
layer follows a typical nucleation and growth mechanism; phase flips at 2 × 2 cells happen to coalesce to form a domain (black circle) and the growth/merge (ellipse) of the domains results in the phase flip of the whole layer.This implies that the finite size effect makes the rate of phase flips to be somewhat overestimated because domains can have more chances of merging while its melting probability is ignored when the size is larger than the simulation cell.Nevertheless, the thermodynamic stability of domain formation is still valid because the initial nucleation process is insensitive to the boundary condition.

Eigenmode decomposition and correlation analysis of CDW
To investigate the fluctuation dynamics of preformed CDW, the instantaneous local CDW order should be distinguished from trivial atomic vibrations.We found the 2 × 2 × 2 CDW in CsV 3 Sb 5 can be accurately described with a few low-energy eigenmodes and 94% of its atomic displacement are projected to three symmetrically equivalent lowest eigenmodes at M points (q 1 M , q 2 M , q 3 M in Fig. S5).For CsV 3 Sb 5 , all q i M (i = 1, 2, 3) are real such that their phases can be absorbed to amplitude m i (i = 1, 2, 3) as a sign, and we will use three signed amplitudes of m i as a descriptor for a local CDW order in 2×2×1 supercell as shown in Fig. 1c.Note we choose only one of the time-reversal pair of eigenmodes because their phases are complex conjugate of each other from the constraint that atomic displacements are real.Three   Typical time evolutions of m 1 are demonstrated in Fig. S6(a) for systems with preformed CDW at 120 K (blue), without CDW at 250 K(red) and with both preformed CDW and slow domain fluctuations at 140 K (green), where the preformation of CDWs can be identified as nonzero average value of m 1 .
The correlation and dynamics of local 2 × 2 CDW are quantified using Pearson correlation coefficients (PCC) r xy defines as where x denotes the average of x and x = m i (R, t), y = m j (R + ∆R, t + ∆t).By sufficiently sampling R and t, the correlation function r ij (∆R, ∆t) can be defined as PCC of x and y.Figures S6(b) and (c) shows the schematic scatter plots of m 1 (R, t) and m 1 (R, t+∆t) and corresponding temporal correlation function r 11 (∆R = 0, ∆t) ≡ r 11 (∆t), respectively.In common, r 11 (∆t) has oscillations with typical phonon frequency.For a system with CDW (blue), the asymptotic value of r 11 (∆t) remains finite while it decays fast to zero (red) if no static CDW exists.If there is slower dynamics such as sign flips of m i (green), it is manifested as the slow decay of r 11 (∆t), which is a good measure of phase flip rate.

Formation energy of CDW domains
We found that the considered CDW domain here is a local minimum in neither our interatomic potential nor thermally averaged potential.So, there seems to be no energy barrier to stabilize the finite-sized CDW domains.Therefore, during the optimization step, we fixed vanadium atoms in a CDW domain and were able to stabilize various shapes of CDW domains.We also have checked that fixing Sb atoms gives the similar domain energy within 5%.
From the formation energy of CDW domain, the intraand inter-layer interaction parameters J ∥ and J ⊥ are calculated as follows.We assume the formation energy E can be accurately approximated as E = AJ ⊥ +LJ ∥ where A and L are the area and circumference (more accurately, the number of neighboring 2×2 cells) of the domain.The fitted J ∥ and J ⊥ are plotted in Fig. S7(c) and (d), respectively, where the shape dependences are negligible.
As T increases, the ratio between inter-and intralayer interaction parameter J ⊥ /J ∥ is also sharply increases as shown in Fig. S8(a).One of the reason for the different temperature dependence of J ∥ and J ⊥ is the different  behaviors of domain wall.Figure S8(b) shows a rhombic CDW domain of CsV 3 Sb 5 at 0 K drawn with the same scheme as Fig. S4.The line profiles along the arrow are compared for the case of T = 0 K and T = 120 K in Fig. S8(c) and (d), respectively, and we can see the width of domain wall is extended to reduce the formation energy when the temperature is high (120 K).On the contrary, the width of domain wall along the out-of-plane direction cannot exceed one unitcell so that no degree of freedom is available to reduce the formation energy for J ⊥ .

Scaled interatomic potentials
In Fig. S9, we display total energy changes of 2 × 2 × 2 CDWs of AV 3 Sb 5 when ground state CDW distortions are linearly scaled.The computed DFT total energies are fit with E = −bρ 2 + cρ 3 + dρ 4 .Fitted vales of c are lesser than 0.003 for all cases and are not shown here.As shown in Fig. S9, the interatomic potential of CsV 3 Sb 5 well reproduces the profile so that we constructed the interatomic potential of KV 3 Sb 5 and RbV 3 Sb 5 by scaling the force constants of CsV 3 Sb 5 .
FIG. 1. 2 × 2 CDW of kagome metals and its low-energy stacking sequences.a, Crystal structure of AV3Sb5 where A denotes alkali atoms of K, Rb and Cs.b, Three eigenmodes of the CDW at M -points of Brillouin zone ( − → M 1, − → M 2, and − → M 3) for kagome shaped lattice of vanadium atoms.Blues dots are vanadium atoms and arrows denotes their displacements corresponding to each − → M i (i = 1, 2, 3).c, Schematic 2 × 2 CDW structures constructed from linear combinations of − → M i's where [m1, m2, m3] ≡ m1 − → M 1 + m2 − → M 2 + m3 − → M 3. Kagome lattices (dashed lines) are modulated into hexagons (blue lines) and triangles (red) forming inverse star of David structure.There are four equivalent phases of 2 × 2 CDW (rhombi with four different colors).d, Low-energy stacking sequences constructed by four phases in c.The phases of CDW are distinguished by colors matching those of rhombi in c.

FIG. 2 .
FIG. 2. Preformation of CDW in CsV3Sb5.a, Thermal distribution of vanadium atom density ρV (r; T ) at T = 200 K. High (low) density regions are colored in yellow (black).b, ρV (r; T ) within the dotted circle in a are enlarged at the temperature range from 20 to 160 K. Green open dots denote kagome lattice points and the space of the grid (dotted lines) is 0.05 Å. c, Temperature-dependent phonon spectra obtained from Sρρ(k, ω; T ).In top left panel, Brillouin zone and high symmetric points of CsV3Sb5 are drawn.The spectrum at T = 20 K is multiplied by 10 for a visual clarity.The arrow in spectrum at T = 250 K indicates a soft mode.Acoustic phonon branches are invisible due to their weak intensities.

FIG. 3 .
FIG. 3. Fluctuation of preformed CDW.a, Temporal fluctuation of the amplitude m1 at a specific 2 × 2 supercell for T = 140 K.The sign changes correspond to phase flips of preformed CDW.b, Temporal Pearson correlation coefficient of r11(∆t) at T = 130 and 140 K. c, Thermal evolution of ρV (r; T ) incorporating slow phase flips.As the temperature decreases, the shape of ρV (r; T ) changes from a Gaussian thermal ellipsoid (right panel, T > T * ) to four-peaked distribution (two middles, TCDW < T < T * ), and finally be collapsed into one of the four peaks (left) at T < TCDW.

FIG. 4 .
FIG. 4. Critical temperatures for CDW.a, Effective lattice structure for 4-state Potts model.Effective 4-state spins locate at vertices of triangles (black circles) corresponding to 2 × 2 supercells shown in Fig. 1c.Effective intralayer (red arrow), first-(blue arrow) and second-nearest-neighboring interlayer (green arrow) interactions between spins are denoted with J ∥ , J ⊥ and J ⊥2 , respectively.b-c, Temperature dependent J ∥ (T ) and J ⊥ (T ) of AV3Sb5 (A = K, Rb, Cs).d, Comparison of the experimental CDW temperatures (Texp) [2] and calculated critical temperatures of T * (open squares) and TCDW (open circles) for different alkali atoms of A on abscissa.

1 FIG
FIG. S4.Typical nucleation and growth mechanism leadingto CDW phase flips.The snapshots of 12 × 12 supercells are taken every 50 femtoseconds and presented from the left top to the right bottom as increasing the time step.The circles and their colors represent 2 × 2 supercells and the amplitudes m1, respectively.The domain starts from cluster (large circle) and grows/merges into one-dimensional strip (ellipse) and finally, the whole layer flips its phase.
FIG. S5. (a)From left to right panels, eigenvectors of three eigenmodes of q i M (i = 1, 2, 3) at M points.Blues dots are kagome lattice points and arrows denotes atomic displacements.Linear combinations of three q i M or m1q 1 M + m2q 2 M + m3q3  M denoted as [m1, m2, m3] form (b) inverse star of David (iSOD) pattern or (c) star of David (SOD) pattern.We note that when a product of three coefficients, m1m2m3, is positive (negative), the resulting configuration becomes iSOD (SOD).
FIG. S6.(a) Time evolutions of the amplitudes m1 of eigenmode qM at T = 120 K (blue), 250 K (red) and 140 K with domain flipping dynamics (green).For former two cases, domain flip dynamics cannot be captured due to short simulation time.(b) schematic scatter plots between m1(R, t) and m1(R, t + ∆t).Large sampling density regions are colored.(c) Pearson correlation coefficients r11(∆R = 0, ∆t).
q i M (i =1, 2, 3) are shown in Fig. S5(a) and 2 × 2 CDW displacement patterns composed of the linear combinations of them are in Fig. S5(b) and (c).When the three amplitudes are the same in the magnitudes, depending on the sign of m 1 m 2 m 3 , they have either star-of-David or inverse star-of-David pattern.
Figure S7(a) is the three types of CDW domains we have considered and Fig. S7(b) summarizes the geometry and the energy of the domains.For each shape of domain, we varied the side length n from 5 to 12 and the calculated domain energies are least-square-fitted with E(n).
FIG. S7.From the left, triangular, rhombic and hexagonal CDW domain on triangular lattices.Each lattice point corresponds 2 × 2 supercell of AV3Sb5 and CDW phases are distinguished by red and black colors.(b) Areas A(n), number of neighboring cells L(n) and energy E(n) of the three domains.n denotes the length of side.(c) intralayer and (d) interlayer interaction parameters calculated from the least square fitting of domain energies.
FIG. S8.(a) Temperature-dependent ratio between inter-and intralayer interaction parameters of KV3Sb5 (red), RbV3Sb5 (green) and CsV3Sb5 (blue).(b) Rhombic CDW domain of CsV3Sb5 at 0 K.Each circle is a 2 × 2 supercell and m1 of the cell is represented with graded colors as shown in the color bar.(c) Line profile indicated by the arrow in (b).Gradual CDW phase changes at the domain wall are reflected in the sign changes of m1 and m2.The unit of horizontal axis is 2a (a is in-plane lattice constant of CsV3Sb5).(d) The same line profile at T = 120 K.The width of domain wall corresponding the slope of m1 (or m2) becomes increased.
FIG. S9.Total energy changes of CDWs in AV3Sb5.From left to right panels, A = K, Rb and Cs, respectivley.Horizontal axes of ρRMS are root mean squares of atomic displacements scaled by maximum values.Blue dots are from DFT and red lines are least square fits on E = −bρ 2 + cρ 3 + dρ 4 .
of SI.