Charge order in the kagome lattice Holstein model: A hybrid Monte Carlo study

The Holstein model describes electrons coupled to local dispersionless phonon modes in a lattice through an on-site electron-phonon interaction [7]. Its Hamiltonian is

	$\displaystyle\hat{H}=$	$\displaystyle-t\sum_{\langle i,j\rangle,\sigma}\left(\hat{c}^{\dagger}_{i\sigma}\hat{c}^{\phantom{\dagger}}_{j\sigma}+h.c.\right)-\mu\sum_{i\sigma}(\hat{n}_{i\sigma}-\tfrac{1}{2})$
		$\displaystyle+\frac{1}{2}\sum_{i}\hat{P}_{i}^{2}+\frac{\omega_{0}^{2}}{2}\sum_{i}\hat{X}_{i}^{2}+\lambda\sum_{i\sigma}\hat{n}_{i\sigma}\hat{X}_{i}\,\,,$		(1)

where $\hat{c}^{\dagger}_{i\sigma}$ $(\hat{c}^{\phantom{\dagger}}_{i\sigma})$ are creation (destruction) operators for an electron at site $i$ with spin $\sigma=\{\uparrow\downarrow\}$, $\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{i\sigma}\hat{c}^{\phantom{\dagger}}_{i\sigma}$ is the electron number operator, and $\mu$ is the chemical potential, which controls the overall filling fraction. The first term describes itinerant electrons hopping between nearest-neighbor sites of the lattice, with a fixed hopping parameter $t=1$ setting the energy scale. In the non-interacting limit, the electronic bandwidth is $W=6$ for the kagome lattice. On each site $i$ are local oscillators of fixed frequency $\omega_{0}$, with $\hat{X}_{i}$ and $\hat{P}_{i}$ the corresponding phonon position and momentum operators, respectively, with the phonon mass normalized to $M=1$. The local electron density $\hat{n}_{i\sigma}$ is coupled to the displacement $\hat{X}_{i}$ through an on-site electron-phonon interaction $\lambda$, which we report here in terms of a dimensionless parameter $\lambda_{\mathrm{D}}=\lambda^{2}/\omega_{0}^{2}\,W$.

The kagome lattice vectors $\mathbf{a}_{1}=(1,0)$ and $\mathbf{a}_{2}=(\frac{1}{2},\frac{\sqrt{3}}{2})$ are shown in Fig. 1(a), with corresponding reciprocal lattice vectors $\mathbf{b}_{1}=(2\pi,-\frac{2\pi}{\sqrt{3}})$ and $\mathbf{b}_{2}=(0,\frac{4\pi}{\sqrt{3}})$, where we have set the lattice constant $a=1$. There are three sites per unit cell with basis vectors $\mathbf{u}_{\mathrm{A}}=(0,0)$, $\mathbf{u}_{\mathrm{B}}=(\frac{1}{2},0)$, and $\mathbf{u}_{\mathrm{C}}=(\frac{1}{4},\frac{\sqrt{3}}{4})$, forming a network of corner sharing triangles with three sublattices, as shown in Fig. 1(a). Each site $i$ may instead be indexed by unit cell and the sublattice $\{A,B,C\}$, such that e.g. $n_{\mathbf{i},\alpha}$ denotes the electron density at the site belonging to sublattice $\alpha$ within the unit cell at position $\mathbf{i}$. In this work, we study finite size lattices with periodic boundary conditions, with linear dimension $L$ (up to $L=15$), $N=L^{2}$ unit cells, and $N_{\mathrm{s}}=3N$ total sites. Note that discrete momentum values are given by $\mathbf{k}=\frac{m_{1}}{L}\mathbf{b}_{1}+\frac{m_{2}}{L}\mathbf{b}_{2}$ where $m_{i}$ is an integer and $0\leq m_{i}<L$.

There are multiple ways to break the sublattice symmetry of the kagome lattice. It is, therefore, important to construct an order parameter that will detect charge ordering independent of the charge distribution within the unit cell. For example, at a filling fraction of $1/3$, electrons may localize by doubly occupying only one site per unit cell, breaking a $Z_{3}$ symmetry. We therefore define an order parameter $\rho_{\mathrm{cdw}}$ that with perfect CDW order takes on one of three values $\mathrm{e}^{\mathrm{i}2\pi\left(\frac{s}{3}\right)}$, where $s=\{0,1,2\}$ corresponds to which way this symmetry is broken. The order parameter $\rho_{\mathrm{cdw}}$ should also be zero in the completely disordered state, where for any unit cell $\bf{i}$ we have $\langle\hat{n}_{\textbf{i},\mathrm{A}}\rangle=\langle\hat{n}_{\textbf{i},\mathrm{B}}\rangle=\langle\hat{n}_{\textbf{i},\mathrm{C}}\rangle$. Hence we define

$$\rho_{\mathrm{cdw}}=\frac{n_{\mathrm{c}}}{2N}\sum_{\textbf{i}}\mathrm{e}^{-\mathrm{i}(\textbf{q}\cdot\textbf{i})}\left(\langle\hat{n}_{\textbf{i},\mathrm{A}}\rangle+\mathrm{e}^{\mathrm{i}\frac{2\pi}{3}}\langle\hat{n}_{\textbf{i},\mathrm{B}}\rangle+\mathrm{e}^{\mathrm{i}\frac{4\pi}{3}}\langle\hat{n}_{\textbf{i},\mathrm{C}}\rangle\right)$$

(2)

where i is a unit cell index, $N$ is the total number of unit cells, q is the ordering wavevector, and $n_{\mathrm{c}}$ is a normalization constant included to fix $|\rho_{\mathrm{cdw}}|=1$ in the case of perfect CDW order. A structure factor that scales with system size can then be defined as $S_{\mathrm{cdw}}(\textbf{q})\propto N\langle|\hat{\rho}_{\mathrm{cdw}}|^{2}\rangle$, where again a proportionality constant can be included to fix $S_{\mathrm{cdw}}(\textbf{q})=N$ for the case of perfect CDW order.

For any pair of sites in the kagome lattice, we denote their density-density correlation in position space by

$$c_{\alpha,\nu}(\textbf{r})=\frac{1}{N}\sum_{\textbf{i}}\langle\hat{n}_{\textbf{i}+\textbf{r},\alpha}\hat{n}_{\textbf{i},\nu}\rangle,$$

(3)

where $\alpha$ and $\nu$ label the sublattice $\{A,B,C\}$ of the two sites, and r is the displacement vector between their unit cells. The Fourier transform of $c_{\alpha,\nu}(\textbf{r})$ gives a generic charge structure factor

$$S_{\alpha,\nu}(\textbf{q})=\sum_{\textbf{r}}e^{\mathrm{i}\textbf{q}\cdot\textbf{r}}c_{\alpha,\nu}(\textbf{r}),$$

(4)

which provides information about the nature of an emergent CDW phase, where $\mathbf{q}$ is a discrete momentum value within the first Brillouin zone. For an ideal CDW pattern with ordering wavevector q, $S_{\alpha,\alpha}(\textbf{q})$ will reach a maximal value proportional to the number of sites, while for $\alpha\neq\nu$ the structure factor will vanish.

In the following section, we show evidence of CDW ordering on the kagome lattice where electrons localize on only one site per unit cell but alternates cyclically between the $\{A,B,C\}$ sublattices from one unit cell to the next. To study the onset of this phase, we set $n_{\mathrm{c}}=1$ in Eq. (2) and define a charge structure factor

	$\displaystyle S_{\mathrm{cdw}}(\textbf{q})$	$\displaystyle=3N\langle|\hat{\rho}_{\mathrm{cdw}}|^{2}\rangle$
		$\displaystyle=\frac{3}{4}\sum_{\alpha}\left(S_{\alpha,\alpha}(\textbf{q})-\frac{1}{2}\sum_{\nu\neq\alpha}S_{\alpha,\nu}(\textbf{q})\right).$		(5)

Additional details are given in Supplementary Discussion. Note that we employ a $\mu$-tuning algorithm [69] to determine the chemical potential for any desired target density.

For the kagome lattice, the noninteracting tight-binding electronic structure with $t>0$ has three separate bands, including one flat band at the highest energy ($E=2t$). The lower bands touch at two inequivalent Dirac points in the Brillouin zone, which we denote $K=(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}})$ and $K^{\prime}=(\frac{4\pi}{3},0)$. The lower band is completely occupied at an average electron density per site of $\langle n\rangle=2/3$ (i.e. an overall filling fraction of $f=1/3$), while the upper band is fully occupied at $\langle n\rangle=4/3$ ($f=2/3$). There are also saddle points in the band structure at the point $M=(\pi,\frac{\pi}{\sqrt{3}})$, which produce singularities in the density of states and sit at the Fermi level for average electron densities of $\langle n\rangle=1/2$ ($f=1/4$) and $\langle n\rangle=5/6$ ($f=5/12$). Fig. 1(b) plots the non-interacting band structure and density of states for the kagome lattice, illustrating these features.

To begin, we study the variation of local quantities as a function of electron density, at fixed $\omega_{0}$ and $\lambda_{\mathrm{D}}$. We set $\omega_{0}/t=0.1$ to facilitate CDW ordering in the Holstein model, as bipolarons should localize more readily in the limit $\omega_{0}/t\rightarrow 0$ due to reduced quantum fluctuations. We also fix a moderate value of the electron-phonon coupling $\lambda_{\mathrm{D}}=0.4$. We will discuss the rationale for this choice of parameters shortly.

In Fig. 2(a) we show the average electron density per site $\langle n\rangle$ as a function of chemical potential $\mu$ for an $L=12$ lattice, as the inverse temperature is varied from $\beta=2\mbox{--}14$. We observe the formation of a plateau at $\langle n\rangle=2/3$ as the temperature is lowered, signalling the opening of a gap. No signatures of CDW ordering is observed at fillings away from $\langle n\rangle=2/3$ for these parameters. We also calculate the average electron kinetic energy as a function of electron density as shown in Fig. 2(b). We observe a sharp change at $\langle n\rangle=2/3$, where the magnitude of the electron kinetic energy becomes maximal. This is a signature of a CDW phase transition, since a configuration of doubly occupied sites surrounded by empty nearest-neighbor sites maximizes the number of bonds along which electron hopping is permitted (and corresponds to an average electron density per site $\langle n\rangle=2/3$ on the kagome lattice). Note that since the kagome lattice is not bipartite, particle-hole symmetry is not present and thus both the kinetic energy and average filling are not symmetric about half-filling.

To further study the opening of a CDW gap as the temperature is lowered, we calculate the momentum integrated spectral function $A(\omega)$, which is related to the imaginary time dependent Green’s function through the integral equation

$$G(\mathbf{k},\tau)=\langle\hat{c}(\mathbf{k},\tau)\hat{c}^{\dagger}(\mathbf{k},\tau)\rangle=\int d\omega A(\mathbf{k},\omega)\frac{\mathrm{e}^{-\omega\tau}}{1+\mathrm{e}^{-\beta\omega}},$$

(6)

which we invert using the maximum entropy method to obtain $A(\omega)$ [70]. In Fig. 3, we show the momentum integrated spectral function for an $L=15$ lattice ($N_{\mathrm{s}}=775$) for a range of temperatures down to $\beta=24$, again fixing $\omega_{0}=0.1$, $\lambda_{\mathrm{D}}=0.4$, and an average electron density per site of $\langle n\rangle=2/3$. We observe three peaks in the spectral function corresponding to the three-band structure. As the temperature is lowered, $A(\omega)$ reaches zero and a finite gap begins to open at $\beta\gtrsim 18$, as shown in the bottom panel, indicating a transition to an insulating CDW phase.

At an average electron density per site of $\langle n\rangle=2/3$, the lower energy band is completely filled and touches the upper band at the Dirac points $K$ and $K^{\prime}$. To study the onset of CDW order at this filling, we therefore calculate the charge structure factor $S_{\mathrm{cdw}}$ [Eq. (Kagome lattice Holstein model)] evaluated at $\mathbf{q}=\mathbf{K}$, as a function of phonon frequency, electron-phonon coupling, and temperature.

In Fig. 4 we show the variation of $S_{\mathrm{cdw}}(\mathbf{K})$ as the phonon frequency $\omega_{0}$ is increased from $0.1$ to $1.0$. In the antiadiabatic limit ($\omega_{0}\rightarrow\infty$), deformation of the lattice is weakened as sites respond more quickly to electron hopping and bipolarons do not readily localize, inhibiting the formation of a stable CDW pattern [22]. In addition, quantum fluctuations are enhanced at large $\omega_{0}$, further suppressing CDW order [23]. For $\omega_{0}\gtrsim 0.4$ we observe no significant growth in $S_{\mathrm{cdw}}(\mathbf{K})$ as the temperature is lowered from $\beta=2$ to $\beta=20$. However, for $\omega_{0}\lesssim 0.3$, the structure factor begins to increase in magnitude as the temperature is reduced, growing more rapidly with $\beta$ as $\omega_{0}\rightarrow 0$. We therefore fix $\omega_{0}=0.1$, and vary the dimensionless electron-phonon coupling $\lambda_{\mathrm{D}}$, in order to determine the region in which CDW order at $\langle n\rangle=2/3$ is most enhanced and subsequently estimate $T_{\mathrm{c}}$ for these parameters.

At small values of $\lambda_{\mathrm{D}}$, we find no enhancement in $S_{\mathrm{cdw}}(\mathbf{K})$ from $\beta=2$ to $\beta=20$ i.e. for $\lambda_{\mathrm{D}}\lesssim 0.3$ there is no sign of CDW order in this temperature range, as shown in Fig. 5. This may be due to the critical temperature becoming exponentially suppressed as $\lambda_{\mathrm{D}}\rightarrow 0$. However, another possibility is a finite $\lambda_{\mathrm{D}}$ is necessary for CDW formation, as is the case in the honeycomb lattice Holstein model at half-filling [25], which similarly has Dirac cones and a vanishing density of states at the Fermi surface. As $\lambda_{\mathrm{D}}$ increases, the effective electron-electron attraction is enhanced, and we observe an increase in the charge structure factor as pairs of electrons arrange themselves into a periodic CDW. As the temperature is reduced, we find that there is a maximum in $S_{\mathrm{cdw}}(\mathbf{K})$ at approximately $\lambda_{\mathrm{D}}\approx 0.4$. At larger $\lambda_{\mathrm{D}}$, the CDW structure factor is smaller, and eventually no significant growth is observed as the temperature is lowered from $\beta=2$ to $\beta=20$. This behavior might originate from the higher effective bipolaron mass at large $\lambda_{\mathrm{D}}$, which will hinder their arrangement into an ordered CDW phase, as the energy barrier associated with moving from site to site is proportional to $\lambda_{\mathrm{D}}$, thus promoting self-trapping. Consequently, $T_{\mathrm{c}}$ rapidly decreases as $\lambda_{\mathrm{D}}$ becomes much larger than its optimal value. We note that similar behavior has been observed in the honeycomb, square, and Lieb lattice Holstein models [25, 27].

The momentum dependence of $S_{\mathrm{cdw}}(\mathbf{q})$ is shown in Fig. 6, where the charge structure factor at $\langle n\rangle=2/3$ is evaluated over the first Brillouin zone for an $L=12$ lattice. An enhancement in the structure factor is observed at the Dirac points as the temperature is lowered, corresponding to the onset of an ordered CDW phase, with the magnitude of $S_{\mathrm{cdw}}$ increasing rapidly around $\beta\gtrsim 17$. For all other momentum values, including at the $M$ and $\Gamma$-points, we find no enhancement in charge correlations with inverse temperature $\beta$, at this filling.

A real-space depiction of the CDW correlations at $\langle n\rangle=2/3$ is shown in Fig. 7, which plots density-density correlations $\langle\hat{n}(\mathbf{r})\hat{n}(\mathbf{0})\rangle$ over an $L=12$ lattice with periodic boundary conditions. Here $\mathbf{r}=0$ is the position of a fixed reference site belonging to the $A$ sublattice. Hence Fig. 7 depicts $c_{\alpha,\nu}(\mathbf{r})$ with the origin fixed at this reference site. The CDW pattern is characterized by the localization of electron pairs on only one site per unit cell, which belongs to either the $A$, $B$, or $C$ sublattice, alternating cyclically between these from one unit cell to the next (in both the $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ directions). The fact that $\mathbf{K}$ and $\mathbf{K}^{\prime}$ are the ordering wavevectors for this pattern can be understood as follows. In terms of the reciprocal lattice vectors, we have $\mathbf{K}=\frac{1}{3}(\mathbf{b}_{1}-\mathbf{b}_{2})$ and $\mathbf{K^{\prime}}=\frac{1}{3}(2\mathbf{b}_{1}+\mathbf{b}_{2})$. If the doubly-occupied sites are separated by a displacement $\mathbf{r}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}$, then the Fourier transform of the density-density correlation function will have peaks at $\mathbf{K}$ or $\mathbf{K^{\prime}}$ if $\mathbf{K}\cdot\mathbf{r}=2m\pi$ or $\mathbf{K^{\prime}}\cdot\mathbf{r}=2m\pi$, where $m\in\mathbb{Z}$. This is satisfied if $(n_{1}-n_{2})\mod 3=0$ (for $\mathbf{K}$) or $(2n_{1}+n_{2})\mod 3=0$ (for $\mathbf{K^{\prime}})$, which are equivalent conditions. In other words, moving along either the $\mathbf{a}_{1}$ or $\mathbf{a}_{2}$ directions, density-density correlations will repeat with a periodicity of three unit cells, i.e. for each unit cell, the site on which the electron pairs localize will alternate cyclically between the $\{A,B,C\}$ sublattices. For any given unit cell, the onset of this type of CDW order therefore breaks a $Z_{3}$ symmetry, and the phase transition should belong to the 3-state Potts model universality class.

For an approximate estimate of the critical temperature we can examine when the correlations become long-ranged on a finite-size lattice. As shown in Fig. 7, at $\beta=16$ the charge order is emerging, but it is not quite long-ranged. At $\beta=20$ however, we see the periodic density correlations persist over the whole lattice. This suggests $T_{\mathrm{c}}$ should lie at an intermediate temperature between these two values; however, for a more accurate estimate, we must study the onset of charge order for several different lattice sizes.

In Fig. 8 we show the variation of the charge structure factor $S_{\mathrm{cdw}}(\mathbf{K})$ with inverse temperature $\beta$, for lattices with linear dimension $L=6,9,12$ and $15$, for a range of temperatures down to $\beta=24$. At high temperatures, $S_{\mathrm{cdw}}(\mathbf{K})$ is relatively small and independent of lattice size. However, as the temperature is reduced, $S_{\mathrm{cdw}}(\mathbf{K})$ grows and becomes dependent on the lattice size for $\beta\gtrsim 18$. This signals that correlations are becoming long-ranged and thus sensitive to system size on a finite lattice, and suggests a critical temperature of $\beta_{\mathrm{c}}\approx 18$. A more accurate determination of $T_{\mathrm{c}}$ can be made by studying the correlation ratio

$$R_{\mathrm{c}}=1-\frac{S_{\mathrm{cdw}}(\mathbf{q}+d\mathbf{q})}{S_{\mathrm{cdw}}(\mathbf{q})},$$

(7)

where the ordering wavevector $\mathbf{q}=\mathbf{K}$ here, and $|d\mathbf{q}|$ is the spacing between discrete momentum values for a lattice of linear dimension $L$. For the kagome lattice we average over the six nearest neighbors of the $K$-point in momentum space to obtain $S(\mathbf{K}+d\mathbf{q})$. The correlation ratio $R_{\mathrm{c}}$ is defined such that in the CDW phase, $R_{\mathrm{c}}\rightarrow 1$ as $L\rightarrow\infty$, (since $S_{\mathrm{cdw}}(\mathbf{q})$ will diverge with $L$ if there is long-range order), while $R_{\mathrm{c}}\rightarrow 0$ if there is no long-range order. When plotted for different lattice sizes, the crossing of $R_{\mathrm{c}}$ curves gives an estimate of the critical point. In Fig. 9 we plot $R_{\mathrm{c}}$ for lattices with $L=6,9,12$, and $15$, for the same parameters as in Fig. 8 ($\langle n\rangle=2/3,\omega_{0}=0.1,\lambda_{\mathrm{D}}=0.4$). There is a crossing at $\beta_{\mathrm{c}}\approx 18$, which is consistent with our previous estimates of $\beta_{\mathrm{c}}$ obtained from observing the opening of a finite gap in $A(\omega)$, the onset of long-ranged density-density correlations, and the temperature at which $S_{\mathrm{cdw}}$ becomes dependent on lattice size.

Thus far we have studied the emergence of CDW order on the kagome lattice at a fixed electron density of $\langle n\rangle=2/3$ per site. This choice was motivated by the observation of a CDW gap at $\langle n\rangle=2/3$ and a sharp change in electron kinetic energy during sweeps of $\mu$ and $\langle n\rangle$, and the fact that this filling corresponds to a completely filled lower band, which meets the middle band at the Dirac points $K$ and $K^{\prime}$. However, we also considered fillings of $\langle n\rangle=1/2$ and $\langle n\rangle=5/6$, i.e. densities at which the saddle points in the non-interacting band structure (at the $M$-points) and their van Hove singularities are at the Fermi energy. We also consider $\langle n\rangle=4/3$, which corresponds to completely filled lower and middle bands, with a quadratic touching at the $\Gamma$-point between the flat and middle bands (see Fig. 1). In all of these cases, we find no evidence for the formation of a CDW. For example, there are no anomalous features in components of the total energy, or any indications of a plateau in the $\langle n\rangle$ vs. $\mu$ plots near these fillings, as shown in Fig. 2. Moreover, as the temperature is lowered ($\beta$ increases) the charge structure factor $S_{\mathrm{cdw}}(\mathbf{q})$ does not grow significantly and remains relatively small in magnitude, as shown in Fig. 10 for several high-symmetry points $\mathbf{q}$ in the Brillouin zone [$\mathbf{\Gamma}=(0,0)$, $\mathbf{K}=(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}})$, and $\mathbf{M}=(\pi,\frac{\pi}{\sqrt{3}})$]. We fix $\omega_{0}=0.1$ here to avoid suppression of any potential CDW order, which occurs in the antiadiabatic limit. These results thus suggest an absence of any charge ordering at these fillings, at least for inverse temperatures $\beta<20$. In other words, our results show no evidence for other varieties of CDW order in the kagome lattice Holstein model other than at the $K$-points at $\langle n\rangle=2/3$.

Previous finite temperature studies of the Holstein model have typically employed DQMC [71, 72]. In this method, the inverse temperature $\beta=L_{\mathrm{t}}\Delta\tau$ is discretized along an imaginary time axis with $L_{\mathrm{t}}$ intervals of length $\Delta\tau$, and the partition function is expressed as $Z=\operatorname{Tr}\mathrm{e}^{-\beta\hat{H}}=\operatorname{Tr}\mathrm{e}^{-\Delta\tau\hat{H}}\mathrm{e}^{-\Delta\tau\hat{H}}\ldots\mathrm{e}^{-\Delta\tau\hat{H}}$. Since Eq. (Kagome lattice Holstein model) is quadratic in fermionic operators, these can be traced out, giving an expression for $Z$ in terms of the product of two identical matrix determinants $\det M(x_{i,\tau})$, which are functions of the space and time-dependent phonon displacement field only. Monte Carlo sampling using local updates to the phonon field $\{x_{i,\tau}\}$ is performed and physical quantities can be measured through the fermion Green’s function $G_{ij}=\langle c_{i}^{\dagger}c_{j}\rangle=\left[M^{-1}\right]_{ij}$. Although there is no sign problem [73] for the Holstein model, these studies have been limited for two main reasons. First, the computational cost of DQMC scales as $N_{\mathrm{s}}^{3}L^{\phantom{3}}_{\mathrm{t}}$, where $N_{\mathrm{s}}$ is the total number of lattice sites, prohibiting the study of large system sizes. Secondly, the restriction to local updates results in long autocorrelation times at small phonon frequencies. This aspect has limited simulations to phonon frequencies of $\omega_{0}\gtrsim t$, which is unrealistic for most real materials, and is far from the regime where CDW order in the Holstein model is typically the strongest ($\omega_{0}\ll t$).

Significant efficiency gains are possible by using a dynamical sampling procedure that updates the entire phonon field at each time-step [74, 75]. In this work, we use a recently developed collection of techniques to perform finite temperature simulations on extremely large clusters [68]. Our HMC-based approach achieves a near-linear scaling with system size [74, 76, 77], allowing us to study lattices of up to $N_{\mathrm{s}}=775$ sites at temperatures as low as $T=t/24$. Our algorithm efficiently updates the phonon field simultaneously, allowing study of a realistic phonon frequency $\omega_{0}/t=0.1$.

Near-linear scaling is achieved by rewriting each matrix determinant $\det M$ as a multi-dimensional Gaussian integral involving auxiliary fields $\Phi_{\sigma}$ that will also be sampled. Here, the partition function becomes

$$\mathcal{Z}\approx\left(2\pi\right)^{N_{\mathrm{s}}L_{\mathrm{t}}}\int\mathcal{D}\Phi_{\uparrow}\mathcal{D}\Phi_{\downarrow}\mathcal{D}x\,\mathrm{e}^{-S(x,\Phi_{\sigma})},$$

(8)

where the total action is

$\displaystyle S(x,\Phi_{\sigma})=$

$\displaystyle S_{\mathrm{B}}(x)+S_{\mathrm{F}}(x,\Phi_{\sigma})$

(9)

with the fermionic (F) and bosonic (B) contributions

	$\displaystyle S_{\mathrm{F}}\left(x,\Phi_{\sigma}\right)=$	$\displaystyle\frac{1}{2}\sum_{\sigma}\Phi_{\sigma}^{T}\left(M^{T}M\right)^{-1}\Phi_{\sigma}$		(10)
	$\displaystyle S_{\mathrm{B}}(x)=$	$\displaystyle\frac{\Delta\tau}{2}\sum_{i,\tau}\left[\omega_{0}^{2}x_{i,\tau}^{2}+\left(\frac{x_{i,\tau+1}-x_{i,\tau}}{\Delta\tau}\right)^{2}\right].$		(11)

A Gibbs sampling procedure is then adopted where $\Phi_{\sigma}$ and $x$ are alternately updated. The auxiliary field $\Phi_{\sigma}$ may be directly sampled. Using HMC, global updates to the phonon fields $x$ can be performed by introducing a conjugate momentum $p$ and evolving a fictitious Hamiltonian dynamics using a symplectic integrator [68].

For bipartite geometries such as the square lattice, checkerboard CDW order can occur in the Holstein model at half-filling. In these cases, electron pairs localize on one of the two sublattices ($A$ or $B$), breaking a $Z_{2}$ symmetry, with ordering wavevector $\mathbf{q}=(\pi,\pi)$. A charge structure factor that scales with system size can be defined as $S_{\mathrm{cdw}}=\sum_{\mathbf{r}}\mathrm{e}^{-\mathrm{i}(\pi,\pi)\cdot\mathbf{r}}c(\mathbf{r})$, where the sum is over all unit cells, and $c(\mathbf{r})=\frac{1}{N}\sum_{\mathbf{i}}\langle\hat{n}_{\mathbf{i}+\mathbf{r}}\hat{n}_{\mathbf{i}}\rangle$ is the real space density-density correlation function. However, one can also express $S_{\mathrm{cdw}}$ in terms of an order parameter $\rho_{\mathrm{cdw}}$ i.e. $S_{\mathrm{cdw}}=N\langle|\rho_{\mathrm{cdw}}|^{2}\rangle$, such that with perfect CDW order $\rho_{\mathrm{cdw}}=\pm 1$ depending on which sublattice the electrons localize on, and $\rho_{\mathrm{cdw}}=0$ in the disordered phase. To detect checkerboard order on the square lattice no matter how the $Z_{2}$ symmetry is broken, we should consider the difference between $\langle\hat{n}_{\mathrm{A}}\rangle$ and $\langle\hat{n}_{\mathrm{B}}\rangle$, i.e. define

	$\displaystyle\rho_{\mathrm{cdw}}$	$\displaystyle=\langle\hat{\rho}_{\mathrm{cdw}}\rangle=\frac{1}{2}\left(\langle\hat{n}_{\mathrm{A}}\rangle-\langle\hat{n}_{\mathrm{B}}\rangle\right)$
		$\displaystyle=\frac{1}{2}\left(\frac{2}{N}\sum_{\mathbf{i}\in A}\langle\hat{n}_{\mathbf{i}}\rangle-\frac{2}{N}\sum_{\mathbf{i}\in B}\langle\hat{n}_{\mathbf{i}}\rangle\right)$
		$\displaystyle=\frac{1}{N}\left(\sum_{\mathbf{i}\in A}\langle\hat{n}_{\mathbf{i}}\rangle-\sum_{\mathbf{i}\in B}\langle\hat{n}_{\mathbf{i}}\rangle\right)$
		$\displaystyle=\frac{1}{N}\sum_{\mathbf{i}}(-1)^{i_{x}+i_{y}}\langle\hat{n}_{i}\rangle$
		$\displaystyle=\frac{1}{N}\sum_{\mathbf{i}}\mathrm{e}^{\mathrm{i}(\pi,\pi)\cdot\mathbf{i}}\langle\hat{n}_{\mathbf{i}}\rangle,$		(S1)

where $\mathbf{i}=i_{x}\hat{\mathbf{x}}+i_{y}\hat{\mathbf{y}}$ are the locations of sites in a square lattice, with the lattice constant normalized to $a=1$. Taking the squared magnitude of $\rho_{\mathrm{cdw}}$, a structure factor that scales with system size can then be expressed as

	$\displaystyle S_{\mathrm{cdw}}$	$\displaystyle=N\langle|\hat{\rho}_{\mathrm{cdw}}|^{2}\rangle=N\langle\hat{\rho}_{\mathrm{cdw}}^{\dagger}\hat{\rho}_{\mathrm{cdw}}\rangle$
		$\displaystyle=N\Bigg{\langle}\bigg{(}\frac{1}{N}\sum_{\mathbf{j}}\mathrm{e}^{-\mathrm{i}(\pi,\pi)\cdot\mathbf{j}}\hat{n}_{\mathbf{j}}\bigg{)}\bigg{(}\frac{1}{N}\sum_{\mathbf{i}}\mathrm{e}^{\mathrm{i}(\pi,\pi)\cdot\mathbf{i}}\hat{n}_{\mathbf{i}}\bigg{)}\Bigg{\rangle}$
		$\displaystyle=\frac{1}{N}\sum_{\mathbf{i},\mathbf{r}}\mathrm{e}^{-\mathrm{i}(\pi,\pi)\cdot\mathbf{r}}\langle\hat{n}_{\mathbf{i}+\mathbf{r}}\hat{n}_{\mathbf{i}}\rangle$
		$\displaystyle=\sum_{\mathbf{r}}\mathrm{e}^{-\mathrm{i}(\pi,\pi)\cdot\mathbf{r}}c(\mathbf{r}),$		(S2)

where $\mathbf{r}=\mathbf{j}-\mathbf{i}$, i.e. the expression in Eq. (A. Derivation of CDW order parameter) is equivalent to the Fourier transform of the real space density-density correlation function $c(\mathbf{r})$ for the square lattice.

For the kagome lattice, since each unit cell consists of three sites, we introduced a generic density-density correlation function $c_{\alpha,\nu}(\mathbf{r})$ which has Fourier transform $S_{\alpha,\nu}(\mathbf{q})$, where each lower index denotes a sublattice $A$, $B$, or $C$. In our paper we discuss a CDW pattern in which electrons localize on one site per unit cell, alternating cyclically between the $A$, $B$, and $C$ sites, as shown in Fig. 7. Unlike checkerboard order on the square lattice, CDW order of this type can occur in three ways, breaking a $Z_{3}$ symmetry. As in the square lattice case, we can define an order parameter $\rho_{\mathrm{cdw}}$ which takes on a different value depending on how the symmetry is broken, but with $|\rho_{\mathrm{cdw}}|=1$ in the case of perfect order and $|\rho_{\mathrm{cdw}}|=0$ in the disordered phase. Since a $Z_{3}$ symmetry is broken, we should have $\rho_{\mathrm{cdw}}=\mathrm{e}^{\mathrm{i}2\pi\left(\frac{s}{3}\right)}$ (where $s=\{0,1,2\}$) in the ordered phase, and $\rho_{\mathrm{cdw}}=0$ in the disordered phase. However, unlike the simpler checkerboard order, the electron densities on each sublattice $\langle\hat{n}_{\mathrm{A}}\rangle$, $\langle\hat{n}_{\mathrm{B}}\rangle$, and $\langle\hat{n}_{\mathrm{C}}\rangle$ will vary from unit cell to unit cell, with a periodicity set by the ordering wavevector $\mathbf{q}$. We can therefore define an order parameter

$$\rho_{\mathrm{cdw}}=\frac{1}{2N}\sum_{\textbf{i}}\mathrm{e}^{-\mathrm{i}(\textbf{q}\cdot\textbf{i})}\left(\langle\hat{n}_{\textbf{i},\mathrm{A}}\rangle+\mathrm{e}^{\mathrm{i}\frac{2\pi}{3}}\langle\hat{n}_{\textbf{i},\mathrm{B}}\rangle+\mathrm{e}^{\mathrm{i}\frac{4\pi}{3}}\langle\hat{n}_{\textbf{i},\mathrm{C}}\rangle\right),$$

(S3)

which satisfies these properties, where the sum is over all unit cells. As before, we can now write a structure factor that scales with system size,

$\displaystyle S_{\mathrm{cdw}}(\mathbf{q})=3N\langle|\hat{\rho}_{\mathrm{cdw}}|^{2}\rangle=3N\langle\hat{\rho}_{\mathrm{cdw}}^{\dagger}\hat{\rho}_{\mathrm{cdw}}^{\phantom{\dagger}}\rangle,$

(S4)

where a constant factor has been inserted to ensure $S_{\mathrm{cdw}}=N$ in the case of perfect CDW order. We can now expand the above equation to obtain an expression for $S_{\mathrm{cdw}}$ in terms of the generic structure factors $S_{\alpha,\nu}$, which yields

	$\displaystyle S_{\mathrm{cdw}}(\mathbf{q})$	$\displaystyle=\frac{3}{4N}\sum_{\mathbf{r}}\bigg{[}\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}}\sum_{\alpha}\bigg{(}c_{\alpha,\alpha}(\mathbf{r})-\frac{1}{2}\sum_{\alpha\neq\nu}c_{\alpha,\nu}(\mathbf{r})\bigg{)}\bigg{]}$
		$\displaystyle=\frac{3}{4}\sum_{\alpha}\bigg{(}S_{\alpha,\alpha}(\mathbf{q})-\frac{1}{2}\sum_{\alpha\neq\nu}S_{\alpha,\nu}(\mathbf{q})\bigg{)},$		(S5)

where $c_{\alpha,\nu}(\mathbf{r})$ is the generic real space density-density correlation function. This expression for $S_{\mathrm{cdw}}$ is the measure we use to detect CDW ordering in the kagome lattice, and is the quantity we show in Fig. 8 at $\langle n\rangle=2/3$ as a function of $\beta$ for different lattice sizes. For a particular CDW pattern on the kagome lattice under study, Eq. (S5) can be multiplied by a constant factor (if necessary) to fix $S_{\mathrm{cdw}}=N$ in the case of perfect CDW order.

An additional signature of the CDW transition can be observed in the renormalized phonon energy $\Omega(\mathbf{q},\mathrm{i}\nu_{n}=0)$, where a softening of the phonon dispersion is expected to occur at the ordering wavevector $\mathbf{q}_{\mathrm{cdw}}$ as the temperature is lowered. The renormalized phonon energy is given by $\Omega(\mathbf{q},0)=[\omega_{0}^{2}+\Pi(\mathbf{q},0)]^{1/2}$, where $\Pi(\mathbf{q},0)$ is a function defined in terms of the momentum-space phonon Green’s function $D(\mathbf{q},\nu_{n})$, through the relation

$$D(\mathbf{q},\nu_{n})=\frac{2\omega_{0}}{(\mathrm{i}\nu_{n})^{2}-\omega_{0}^{2}-\Pi(\mathbf{q},\nu_{n})},$$

(S6)

where $\nu_{n}=2\pi nT$ (and we set $\hbar=1$). In Fig. S1 we show $\Omega(\mathbf{q},0)/\omega_{0}$ along a closed triangular path $\Gamma$–$K$–$K^{\prime}$–$\Gamma$ within the Brillouin zone. Results are shown for an $L=12$ lattice with electron-phonon coupling $\lambda_{\mathrm{D}}=0.4$ and bare phonon frequency $\omega_{0}=0.1$, at a filling of $\langle n\rangle=2/3$, where we previously observed evidence of CDW ordering.

At high temperature ($\beta\ll\beta_{\mathrm{c}}$) we find that the renormalized phonon dispersion is relatively flat, e.g. for $\beta=2$, we have $\Omega(\mathbf{q},0)/\omega_{0}\approx 0.6$ for all momenta $\mathbf{q}$. This is indicative of strong electron-phonon coupling, where the renormalized phonon frequency is uniformly suppressed relative to the bare phonon frequency even at temperatures well above $T_{\mathrm{c}}$. This signature of strong coupling is expected behavior in the Holstein model, and has been observed in the square lattice at similarly large values of $\lambda_{\mathrm{D}}$ [19]. As the temperature is lowered, we observe a softening of the phonon dispersion at the expected ordering wavevectors $K$ and $K^{\prime}$. We observe sharp dips in the phonon dispersion as the temperature is reduced to $\beta\approx 18$, consistent with our estimate of $\beta_{\mathrm{c}}$ obtained from the crossing of $R_{\mathrm{c}}$ curves shown in Fig. 9. Due to an expected finite-size effect, $\Omega(\mathbf{q},0)$ does not reach exactly zero below $T_{\mathrm{c}}$ in our simulations [77].

Here we provide additional details of the parameter values used in our HMC simulations, which are explained further in Ref. [68]. We fixed the imaginary-time discretization at $\Delta\tau=0.05$, and performed updates to the phonon field using HMC trajectories typically of $N_{\mathrm{t}}=50$ time-steps of size $\Delta t=0.2$. We initially thermalize our systems by performing $N_{\mathrm{therm}}=2000\mbox{--}3000$ trial updates. This is followed by $N_{\mathrm{sim}}=2500\mbox{--}4000$ updates, where measurements are taken at the end of each trajectory. Each time-step we numerically solve for the fictitious force $\dot{p}=-\frac{\partial S}{\partial x}$ using the conjugate gradient method with a relative residual error threshold $\epsilon_{\mathrm{max}}=10^{-5}$.

The forces in our HMC dynamics can be separated into fermionic and bosonic components $-\frac{\partial S_{\mathrm{F}}}{\partial x}$ and $-\frac{\partial S_{\mathrm{B}}}{\partial x}$ [see Eqs. (10) and (11)], where the bosonic part is much less expensive to calculate. We thus employ time-step splitting where trajectories evolve with time-step $\Delta t^{\prime}=\Delta t/n_{\mathrm{t}}$ with $n_{\mathrm{t}}=10$ using the bosonic force alone, followed by a single step of size $\Delta t$ using the fermionic force. We also use Fourier acceleration via a dynamical mass matrix with regularization parameter $m_{\mathrm{reg}}=\omega_{0}$ to further reduce autocorrelation times.