# Electronic instabilities of kagomé metals: saddle points and Landau theory

Takamori Park  Department of Physics, University of California, Santa Barbara,
CA 93106-9530, USA  Mengxing Ye  Kavli Institute for Theoretical Physics,
University of California, Santa Barbara, CA 93106-4030, USA  Leon Balents
Kavli Institute for Theoretical Physics, University of California, Santa
Barbara, CA 93106-4030, USA  Canadian Institute for Advanced Research, 661
University Ave., Toronto, ON M5G 1M1 Canada

###### Abstract

We study electronic instabilities of a kagomé metal with a Fermi energy close
to saddle points at the hexagonal Brillouin zone face centers. Using parquet
renormalization group, we determine the leading and subleading instabilities,
finding superconducting, charge, orbital moment, and spin density waves. We
then derive and use Landau theory to discuss how different primary density
wave orders give rise to charge density wave modulations, as seen in the
AV3Sb5 family, with A=K,Rb,Cs. The results provide strong constraints on the
mechanism of charge ordering and how it can be further refined from existing
and future experiments.

##  I Introduction

Two dimensional correlated metals based on transition metal ions are a classic
subject for many body physics [1, 2], displaying diverse electronic phenomena
such as unconventional superconductivity [3, 4, 5], charge and spin order [6,
7], nematicity [8, 9], strange metallic behavior [10, 11, 12], and more. These
topics have been most heavily investigated in theory and experiment in
structures based on square lattices in the cuprates [4, 13] and Fe
superconductors [3, 14, 15]. Correlated metals with hexagonal/triangular
symmetry are much less common, with the best-known example being the
triangular lattice cobaltates [16, 17], which however are complicated by Na
vacancy ordering and water intercalation[18].

Recently, a new class of kagomé metals, with chemical formula AV3Sb5, where A
= K, Rb, or Cs, have emerged as an exciting realization of quasi-2D correlated
metals with hexagonal symmetry [19]. These materials have been shown to
display several electronic orders setting in through thermodynamic phase
transitions: multi-component (“3Q”) hexagonal charge density wave (CDW) order
below a $T_{c}\approx 90K$ [20, 21, 22, 23, 24, 25, 26, 27], and
superconductivity with critical temperature of $2.5K$ or smaller [28, 20, 21,
29, 30, 31, 32, 33, 34, 35], and some indications of nematicity and one-
dimensional charge order in the normal and superconducting states [21, 31,
36]. Other experiments show a strong anomalous Hall effect [37, 38],
suggesting possible topological physics. Furthermore, density functional
calculations identified CsV3Sb5 as a $\mathbb{Z}_{2}$ topological metal [39].
Angle resolved photoemission studies show these materials to be multi-band
systems with several Fermi surface components [39], including approximately
nested components and a Fermi energy that is close to multiple saddle points
of the dispersion, in agreement with density functional theory [39, 40].
Furthermore, strong momentum dependent charge gaps near the saddle point
momenta were observed below the CDW transition temperature [41, 42].

Many of these ingredients bring to mind a storied idea of electronic
instabilities enhanced by van Hove singularities near saddle points. This
mechanism figured heavily in early theoretical treatments of the cuprates [43,
44, 45, and references therein]. In a two-dimensional system, the divergence
of the density of states generates enhanced scattering amongst electrons near
the saddle points of the bands, which may drive not only superconducting but
also other charge and spin instabilities. The same idea emerged recently in
the context of doped graphene [46, 47, 48], and has been applied to the theory
of magic angle graphene bilayers [49, 50, 51, 52, 53]. The observations in
AV3Sb5 suggest another application. Notably, the observed period of charge
order in AV3Sb5 within the two-dimensional kagomé plane is precisely that
expected from scattering amongst the saddle points located at the hexagonal
Brillouin zone face centers (denoted as ${\bm{M}}$).

While the saddle point model on a 2D hexagonal lattice has been studied
extensively in the literature [46, 48, 47, 49], its application to a coherent
understanding of the electronic instabilities revealed in multiple experiments
remains to be understood. First, the renormalization group studies from
repulsive interactions [46, 49] suggest leading spin density wave,
superconductivity, and orbital moment instabilities, rather than the charge
density wave order which is observed. Therefore it is natural to expect that
the lattice plays some role, and thus it is important to understand the
combined effect of electron-phonon and electron correlation. Second, the
strong anomalous Hall effect observed below the charge density wave critical
temperature suggests time-reversal symmetry breaking order may develop below
$T_{c}$. It would be desirable to understand the interplay between time-
reversal symmetric and broken CDW order. Third, in addition to the $2\times 2$
charge density wave within the 2D layer, STM and x-ray studies also observed
modulation in the z-direction, i.e. $k_{z}\neq 0$, in RbV3Sb5 and CsV3Sb5. The
three-dimensional alignment of 2D charge density waves has not yet been
studied theoretically.

In this paper, we explore the electronic instabilities due to interactions
amongst electrons near the saddle points with hexagonal symmetry. We apply the
parquet renormalization group scheme [54, 55, 56, 57] to determine all the
primary and secondary instabilities, generalizing prior work [46, 49]. We
further use mean field theory to study different possible emergent density
wave states, which include not only a conventional charge density wave, but
also spin and orbital moment density waves. We investigate how they relate to
the observed charge density wave measured in AV3Sb5. Through this analysis, we
provide constraints on the interpretation of experimental observations, and
suggestions for future theoretical and experimental studies.

The rest of this paper is structured as follows: In Sec. II, we construct a
low energy continuum model by taking patches around the saddle point momenta
${\bm{M}}$ and identify its connection to the real space tight-binding model.
In Sec. III, we generalize the parquet renormalization group formulation, and
discuss all stable fixed points within the patch model. In Sec. IV, we analyze
the mean field theory for each fixed point solutions. This includes a Landau
theory analysis of the conventional charge density wave (Sec. IV.2) as well as
an analysis of other leading instabilities, i.e. the orbital moment density
wave (Sec. IV.3) and spin density wave (Sec. IV.4), and their interplay with
the conventional charge density wave. In Sec. V, we discuss the implications
of the results obtained in Sec. IV to experimental observations. The real
space pattern on the kagomé lattice are shown in Sec. V.1. To understand the
three dimensional 3Q CDW order with $k_{z}\neq 0$, we study the effects of
weak interlayer coupling and discuss the allowed $k_{z}$ for both conventional
CDW and orbital moment density wave orders. The experimentally observed
anisotropic CDW order can also be explained within this picture. In Sec. V.3,
we discuss the critical behavior for different types of charge order
instabilities. The staggered and uniform orbital moment is estimated in Sec.
V.4. The possibility of a magnetic field induced mixture of conventional CDW
and orbital moment is also briefly discussed there. Lastly, a summary of our
work and discussions on the theoretical and experimental implications are
presented in Sec. VI.

##  II Model

In this section, we introduce a simple model to study the electronic behavior
in AV3Sb5. As explained in the introduction, STM and x-ray scattering
measurements show a $2\times 2$ CDW order with wavevector
$\bm{Q}_{\alpha}=\bm{Q}_{\textrm{Bragg}}/2$ [21, 20, 24, 22]. These
wavevectors are equivalent to the momenta that connect the three
${\bm{M}}_{\alpha}$ points in the hexagonal Brillouin zone as shown in Fig.
1(b). In addition, DFT calculations show that the band structure has saddle
points at the ${\bm{M}}_{\alpha}$ points near the Fermi level as shown in Fig.
2 for CsV3Sb5. In two-dimensional systems, a saddle point is a van Hove
singularity with a logarithmically diverging density of states. Based on these
two observations, we assume that the collective electronic behavior in AV3Sb5
is determined by the saddle points located at the ${\bm{M}}_{\alpha}$ points
and the interactions between them.

Figure 1: 1Vanadium kagomé lattice structure in AV3Sb5. The actual lattice
structure has a Sb atom in the center of the hexagonal faces, but they are not
shown here since our work is only focused on the electronic behavior of
vanadium. The space group of the full lattice is $P6/mmm$. 1 Hexagonal
Brillouin zone of AV3Sb5. There are three distinct M points along the
Brillouin zone boundaries. The ${\bm{M}}_{\alpha}$ points in the Brillouin
zone are connected by three nesting vectors $\\{\bm{Q}_{\alpha}\\}$. The
nesting vectors satisfy $\bm{Q}_{\alpha}\equiv-\bm{Q}_{\alpha}$ up to a
reciprocal lattice vector. Figure 2: Band structure of CsV3Sb5 calculated
using DFT without spin-orbit coupling [58]. The Fermi energy is set to
$6.81\text{\,}\mathrm{eV}$. $\Gamma_{1}^{+},\Gamma_{3}^{+}$ are the
irreducible representations of the two saddle point bands nearest the Fermi
energy at ${\bm{M}}_{1}$. The little co-group at ${\bm{M}}_{1}$ is $D_{2h}$,
and the orientation of the symmetry axes
$C_{2},C_{2}^{\prime},C_{2}^{\prime\prime}$ are shown in Fig. 1. We use the
notation found in Ref. [59].

From these considerations, we construct a non-interacting low-energy continuum
model by taking patches around the ${\bm{M}}_{\alpha}$ points in the Brillouin
zone with cutoff radius $\Lambda$. The Hamiltonian of such a model is

$H_{0}=\sum_{\alpha}\sum_{\absolutevalue{\bm{q}}<\Lambda}c^{\dagger}_{\alpha\bm{q}}\left[\varepsilon_{\alpha}(\bm{q})-\mu\right]c_{\alpha\bm{q}},\\\
$ (1)

where $\alpha,\beta$ are patch indices, $\bm{q}$ is the momentum measured from
${\bm{M}}_{\alpha}$, $\varepsilon_{\alpha}(\bm{q})$ is the saddle point
dispersion sitting at ${\bm{M}}_{\alpha}$, and $\mu$ is the chemical potential
measured away from the saddle point. Up to quadratic order in $\bm{q}$, the
saddle point dispersions take the general form,

$\displaystyle\varepsilon_{1}(\bm{q})=$ $\displaystyle aq_{x}^{2}-bq_{y}^{2},$
$\displaystyle\varepsilon_{2}(\bm{q})=$
$\displaystyle\frac{a-3b}{4}q_{x}^{2}+\frac{\sqrt{3}(a+b)}{2}q_{x}q_{y}+\frac{3a-b}{4}q_{y}^{2},$
$\displaystyle\varepsilon_{3}(\bm{q})=$
$\displaystyle\frac{a-3b}{4}q_{x}^{2}-\frac{\sqrt{3}(a+b)}{2}q_{x}q_{y}+\frac{3a-b}{4}q_{y}^{2},$
(2)

where $\bm{q}=(q_{x},q_{y})$. The parameters $a,b$ that determine the shape of
the saddle point must have the same sign, $ab>0$. Note, that the dispersion of
the three patches are related by a three-fold rotation. The condition for
perfect nesting is given by $a/b=3$, but in our work we will not necessarily
impose this condition unless otherwise stated.

In general, there are four different cases of the form of the dispersion that
we can consider: (i) $a,b>0$, $a/b>1$, (ii) $a,b>0$, $a/b<1$, (iii) $a,b<0$,
$a/b>1$, and (iv) $a,b<0$, $a/b<1$. These cases can be all be mapped to case
(i) by using two transformations. First, rotating the coordinates by $\pi/2$
takes cases (iii), (iv) to cases (ii), (i) respectively. Second, a particle-
hole transformation combined with the same $\pi/2$ rotation maps case (ii) to
case (i) (changing the sign of $\mu$). Hence, we only need to consider case
(i). Note that the latter particle-hole-like transformation becomes a symmetry
when $a=b$ and $\mu=0$.

Next, we introduce interactions to the continuum model by listing all possible
electron-electron interactions between the fermions in the three patches.
There are four such interactions,

$\displaystyle H_{1}=$
$\displaystyle\frac{1}{2\mathcal{N}}\sideset{}{{}^{\prime}}{\sum}_{\absolutevalue{\bm{q}_{1}},\cdots,\absolutevalue{\bm{q}_{4}}<\Lambda}\bigg{[}\sum_{\alpha\neq\beta}\big{(}g_{1}c^{\dagger}_{\alpha\bm{q}_{1}\sigma}c^{\dagger}_{\beta\bm{q}_{2}\sigma^{\prime}}c_{\alpha\bm{q}_{3}\sigma^{\prime}}c_{\beta\bm{q}_{4}\sigma}$
$\displaystyle\qquad+g_{2}c^{\dagger}_{\alpha\bm{q}_{1}\sigma}c^{\dagger}_{\beta\bm{q}_{2}\sigma^{\prime}}c_{\beta\bm{q}_{3}\sigma^{\prime}}c_{\alpha\bm{q}_{4}\sigma}$
$\displaystyle\qquad+g_{3}c^{\dagger}_{\alpha\bm{q}_{1}\sigma}c^{\dagger}_{\alpha\bm{q}_{2}\sigma^{\prime}}c_{\beta\bm{q}_{3}\sigma^{\prime}}c_{\beta\bm{q}_{4}\sigma}\big{)}$
$\displaystyle\qquad+\sum_{\alpha}g_{4}c^{\dagger}_{\alpha\bm{q}_{1}\sigma}c^{\dagger}_{\alpha\bm{q}_{2}\sigma^{\prime}}c_{\alpha\bm{q}_{3}\sigma^{\prime}}c_{\alpha\bm{q}_{4}\sigma}\bigg{]},$
(3)

where
$\sideset{}{{}^{\prime}}{\sum}_{\absolutevalue{\bm{q}_{1}},\cdots,\absolutevalue{\bm{q}_{4}}<\Lambda}\equiv\sum_{\absolutevalue{\bm{q}_{1}},\cdots,\absolutevalue{\bm{q}_{4}}<\Lambda}\delta_{\bm{q}_{1}+\bm{q}_{2},\bm{q}_{3}+\bm{q}_{4}}$,
$\mathcal{N}$ is the number of unit cells in the system, and $g_{i}$ are the
interactions that are defined to be intrinsic and have units of energy. As
seen in Fig. 3, the $g_{1},g_{2},g_{3},g_{4}$ interactions represent inter-
patch exchange, inter-patch density-density, Umklapp, and intra-patch density-
density scattering processes, respectively.

The patch model introduced here has been used in previous works to study the
interaction and competition between different instabilities in the Hubbard
model and doped graphene [60, 46]. An important thing to note is that the
continuum model defined above does not include details on the underlying
structure of the lattice and the global band structure. We supplement the
results of our model with the symmetry information provided by the DFT results
for CsV3Sb5 shown in Fig. 2 [58]. This gives a one-to-one correspondence
between Brillouin zone patches, ${\bm{M}}_{\alpha}$, and the vanadium sites,
$V_{\alpha}$ (see App. A). A minimal tight binding model on the kagomé lattice
can be readily obtained from this correspondence to faithfully describe the
saddle point fermions (see Sec. V.1.2 for more discussions).

Figure 3: All possible interactions in the patch model. The cones represent
the saddle point dispersion at each ${\bm{M}}_{\alpha}$ and the arrows denote
the scattering processes described by the interactions.

##  III Renormalization group analysis

Figure 4: RG phase diagram for $g_{2}^{(0)}>0$. The leading instabilities of
each phase under RG are, Phase I: sSC, rCDW, iSDW, Phase II: iCDW and rCDW,
and Phase III: dSC, rSDW, iCDW. The phase diagram was calculated assuming
$d_{1}=1/2\,,g_{2}^{(0)}=0.1\,,g_{4}^{(0)}=0.3$.

While all four interactions $g_{1},g_{2},g_{3}$ and $g_{4}$ are marginal at
tree level, they acquire leading double logarithmic corrections from particle-
particle fluctuations at zero momentum ($\Pi_{\rm pp}(0)$) and particle-hole
fluctuations at momentum transfer ${\bm{M}}_{\alpha}$ ($\Pi_{\rm
ph}({\bm{M}}_{\alpha})$) at the one loop level and so can become marginally
relevant. To study the evolution of $g_{i}$ at energy $E$ as the fermion
fluctuations from cutoff energy $\Lambda_{\rm RG}\sim t\Lambda^{2}$ to $E$ are
integrated out, we perform a parquet renormalization group (pRG) analysis.
Following Ref. [46], the pRG equations for $g_{i}$ are

$\displaystyle\frac{\mathop{}\\!\mathrm{d}g_{1}}{\mathop{}\\!\mathrm{d}y}$
$\displaystyle=2d_{1}g_{1}(g_{2}-g_{1}),$
$\displaystyle\frac{\mathop{}\\!\mathrm{d}g_{2}}{\mathop{}\\!\mathrm{d}y}$
$\displaystyle=d_{1}(g_{2}^{2}+g_{3}^{2}),$
$\displaystyle\frac{\mathop{}\\!\mathrm{d}g_{3}}{\mathop{}\\!\mathrm{d}y}$
$\displaystyle=-g_{3}^{2}-2g_{3}g_{4}+2d_{1}g_{3}(2g_{2}-g_{1}),$
$\displaystyle\frac{\mathop{}\\!\mathrm{d}g_{4}}{\mathop{}\\!\mathrm{d}y}$
$\displaystyle=-2g_{3}^{2}-g_{4}^{2},$ (4)

where $y=\Pi_{\rm pp}(0,E)\sim\ln^{2}(\Lambda_{\rm RG}/E)$ is the RG time and
$d_{1}(y)=\mathop{}\\!\mathrm{d}\Pi_{\rm ph}(M_{a})/\mathop{}\\!\mathrm{d}y$
is the “nesting parameter”, which satisfies $0<d_{1}(y)\leqslant 1/2$. For
perfect nesting, $d_{1}(y)=1/2$ and is independent of RG time. For
concreteness, we will consider the perfect nesting case hereafter, i.e.
$d_{1}\equiv 1/2$. It has been checked that other $0<d_{1}(y_{c})<1/2$ does
not qualitatively change the fixed point solutions and leading density wave
instabilities.

The RG equations, having entirely quadratic $\beta$-functions (the right hand
sides of the pRG equations), do not have any non-trivial controlled fixed
points in the usual sense. Rather, they describe flows in the vicinity of the
free fixed point: since the RG equations are perturbative, they are strictly
valid only within some sphere of small radius of the origin in $g$-space.
Within certain domains of this space, the flows will be unstable, i.e. they
will exit the sphere of control, which indicates an instability towards a new
regime, and most likely an ordered state. Within the unstable regions of phase
space, RG flows that begin very close to the origin tend to converge toward
particular unstable trajectories which act as attractors, and are typically
straight “rays”[61, 62]. Below we follow previous works that reformulate these
rays to appear as fixed points, by projecting the trajectories to a plane of
constant value of one of the parameters. The stable “fixed rays” are expected
to describe the asymptotic limit of arbitrarily weak but non-zero bare
interaction, such that convergence to these rays is nearly perfect before the
sphere of control is exited. One should keep in mind that when the bare
interactions are small but not arbitrarily so, the deviations from these rays
become important, and the physics will be less universal and controlled more
by the actual values of the interactions, but the RG equations can still be
employed.

Before considering the stable fixed rays, we note a few features of these
equations. The $\beta$-functions for $g_{1}$ and $g_{3}$ contain an overall
factor of $g_{1}$ and $g_{3}$, respectively. This follows from symmetry: all
terms except $g_{1}$ conserve spin at each saddle point separately, and all
terms except $g_{3}$ conserve the number of electrons at each saddle point
separately. The conserving interactions cannot generate a non-conserving one.
As a result, the sign of $g_{1}$ and $g_{3}$ remain fixed through out the RG
evolution. One also notes that the $dg_{2}/dy>0$ and $dg_{4}/dy<0$ under the
RG. Thus an initially positive $g_{2}$ must remain positive, and an initially
negative $g_{4}$ remains negative.

To understand the physical meaning of the $g_{i}$, it is useful to define
interactions that parametrize particular channels of ordering, e.g. they
appear in the mean-field treatment below in Sec. IV. They were previously
defined in Ref. [46, 49]. Here, we consider the interactions in the d-wave
superconductivity ($G_{\rm dSC}=g_{3}-g_{4}$), s-wave superconductivity
($G_{\rm sSC}=-2g_{3}-g_{4}$), real charge density ($G_{\rm
rCDW}=-2g_{1}+g_{2}-g_{3}$), orbital moment density ($G_{\rm
iCDW}=-2g_{1}+g_{2}+g_{3}$) 111Here, iCDW standards for “imaginary charge
density wave”. However, note that the latter does not necessarily mean that
the charge instability must break time reversal symmetry for a generic wave
vector. But at wave vector ${\bm{M}}=-{\bm{M}}$ up to a reciprocal lattice
vector, imaginary charge density must break time-reversal symmetry, and
correspond to orbital moment density wave. For similar reasoning, we use iSDW
for spin flux order., real spin density ($G_{\rm rSDW}=g_{2}+g_{3}$), spin
flux order ($G_{\rm iSDW}=g_{2}-g_{3}$) channels. The interactions $G_{a}$ are
defined such that the $a$ instability develops only when $G_{a}>0$.

From the previous discussion, we can see a few features clearly. Real and
imaginary parts of the CDW and SDW order parameters (OPs) are degenerate if
$g_{3}=0$. This is because the umklapp interaction is the only one
transferring charge between saddle points, so that in its absence there is a
separate $U(1)$ charge rotation for each valley. Similarly, the corresponding
real and imaginary parts of the CDW and SDW orders are degenerate when
$g_{1}=0$. This is because only $g_{1}$ violates separate spin conservation at
each saddle, so that when $g_{1}=0$, independent SU(2) rotations may be made
for each flavor, which mixes CDW and SDW orders. Thus, we see that the sign of
$g_{3}$ decides between real and imaginary CDW/SDW, while the sign of $g_{1}$
decides between CDW and SDW order.

To proceed further, we determine the fixed rays and the pRG flow trajectory
near them. We rewrite interactions as $g_{i}=\gamma_{i}\mathsf{g}$, and choose
$\mathsf{g}$ as one of the interactions, which diverges as
$\mathsf{g}\sim\frac{1}{y_{c}-y}$ along the fixed trajectory (as we verify
afterwards). A proper identification of $\mathsf{g}$ ensures that $\gamma_{i}$
tends to a constant value $\gamma_{i}^{*}$ along the fixed trajectory, and
$\gamma_{i}^{*}=0$ implies that the interaction $g_{i}$ either flows to zero
or diverges slower than $\frac{1}{y_{c}-y}$. We call it a fixed point
hereafter. The solutions to $\gamma_{i}^{*}$ can be obtained by solving a set
of algebraic equations $\tilde{\beta}_{i}(\\{\gamma\\})=0$ for $i=1,2,3,4$,
where

$\displaystyle\tilde{\beta}_{i}(\\{\gamma\\})=\dot{\gamma}_{i}=\frac{1}{\mathsf{g}}\left(\dot{g}_{i}-\gamma_{i}\dot{\mathsf{g}}\right).$
(5)

Among those solutions, only the stable fixed point solutions are of physical
interest, as they do not flow away under small perturbations. To examine the
stability of a fixed point solution, we define a matrix $\mathsf{T}$ that is
determined by the flow of $\tilde{\beta}_{i}$ at the fixed point, i.e.
$\mathsf{T}_{ij}=\partial\tilde{\beta}_{i}/\partial\gamma_{j}|_{\\{\gamma^{*}\\}}$.
The RG fixed point is stable only when all the eigenvalues of $\mathsf{T}$ are
non-positive.

In addition to identifying the leading instabilities, the subleading ones are
also considered here for three reasons. _First_ , as discussed above, when the
interaction strengths are not truly infinitesimal, flow to the unstable regime
may occur before the fixed ray is reached. This may occur when bare couplings
are small but not too close to a fixed ray, and deviate from it in a
particular direction favoring a subleading instability. _Second_ , the flow
may be also cut off by imperfect nesting and/or a non-zero chemical potential
(i.e. doping from the saddle point). These effects limit the convergence to
the fixed ray, and at that point a different instability might dominate.
_Third_ , a secondary instability may develop at a lower temperature, even if
the primary one occurs first. To account for those possibilities, which depend
upon the microscopic details, we will list the leading two instabilities at
the RG fixed points.

To summarize, we list all the stable fixed point solutions for both repulsive
and attractive bare interactions. In addition, we discuss an interesting semi-
stable fixed point solution with only one weak unstable direction flowing out
of the fixed point.

When $g_{2}^{(0)}>0$, there are three (semi)stable fixed points, which we take
the liberty of denoting “phases” I, II, III – this is an abuse of terminology
since these solutions really describe unstable rays in the full phase space,
which may not correspond to a unique phase. As $g_{2}(y)$ must diverge as
$\sim\frac{1}{y_{c}-y}$, we choose
$\mathsf{g}(y)=g_{2}(y)=\frac{1}{d_{1}(1+\gamma_{3}^{2})(y_{c}-y)}$.

  * I.

When $g_{3}^{(0)}<0$, $g_{3}$ flows to negative value at the stable RG fixed
point, and we find
$\gamma_{1}=0,\gamma_{2}=1,\gamma_{3}\approx-6.1,\gamma_{4}\approx-5.5$. The
subleading divergence of $g_{1}$ goes as
$g_{1}(y)\sim(y_{c}-y)^{-\frac{2}{1+\gamma_{3}^{2}}}\sim\frac{1}{(y_{c}-y)^{0.05}}$.
Note that as the flow of $g_{1}$ is subleading here, the sign of $g_{1}^{(0)}$
is not qualitatively important to determine the fixed trajectory. The leading
instabilities are $G_{\rm sSC}=17.6\mathsf{g},G_{\rm rCDW}=G_{\rm
iSDW}=7.1\mathsf{g}$.

  * II.

When $g_{1}^{(0)}<0,g_{3}^{(0)}>0$, and for large enough $|g_{1}^{(0)}|$, the
system may flow to a semi-stable fixed point, where there is only 1 weak
unstable direction in the 4-dimensional parameter space defined by
$\\{g_{1},g_{2},g_{3},g_{4}\\}$. The fixed point solution reads
$\gamma_{1}\approx-46.7,\gamma_{2}=1,\gamma_{3}\approx
9.7,\gamma_{4}\approx-4.4$. The leading instabilities are $G_{\rm
iCDW}=104.0\mathsf{g},G_{\rm rCDW}=84.6\mathsf{g}$.

  * III.

When $g_{1}^{(0)},g_{3}^{(0)}>0$, the stable RG fixed point has been discussed
a lot in the literature [46]. The fixed point solution gives
$\gamma_{1}=0,\gamma_{2}=1,\gamma_{3}\approx 5.6,\gamma_{4}\approx-10.0$,
where the subleading divergence of $g_{1}$ can be obtained as
$g_{1}(y)\sim(y_{c}-y)^{-\frac{2}{1+\gamma_{3}^{2}}}\approx\frac{1}{(y_{c}-y)^{0.06}}$.
The leading divergent instabilities are $G_{\rm dSC}=15.6\mathsf{g},G_{\rm
rSDW}=G_{\rm iCDW}=6.6\mathsf{g}$.

In Fig. 4, the phase diagram in the space of $g_{1}^{(0)},g_{3}^{(0)}$ for
$g_{2}^{(0)}>0$ is shown.

When $g_{2}^{(0)}<0$, $g_{2}$ may instead flow to zero. We find in this way a
fourth fixed ray, which we denote phase IV. It is describe by letting
$\mathsf{g}=g_{1}(y)=-\frac{1}{2d_{1}(y_{c}-y)}$, and
$\gamma_{1}=1,\gamma_{2}=0,\gamma_{3}=0,\gamma_{4}=0$. This solution requires
$g_{1}^{(0)}<0,g_{4}^{(0)}>0$, and is only stable when $g_{3}^{(0)}>0$.
Moreover, $g_{3}(y)$ is also divergent and is only logarithmically smaller
than $g_{1}$. We find
$g_{3}(y)=\frac{1}{y_{c}-y}\left(\ln\frac{1}{y_{c}-y}\right)^{-1}$. The
leading instabilities are $G_{\rm rCDW}=G_{\rm
iCDW}=-2\mathsf{g}=\frac{1}{d_{1}(y_{c}-y)}$. The fixed point solution
indicates that the rCDW and iCDW orders are degenerate at the leading order,
but they are split by a logarithmically subdominant effect due to $g_{3}>0$
weakly in favor of the iCDW, i.e. $G_{\rm iCDW}>G_{\rm rCDW}$.

In passing, we note that purely electronic interactions generally give
repulsion for all couplings, i.e. $g_{i}^{(0)}>0$ with $i=1,2,3,4$. However,
other factors, such as orbital composition of the wave function near saddle
points, and electron-phonon coupling, may contribute to attraction for certain
$g_{i}^{(0)}$. In App. B, we consider the effect of electron-phonon coupling,
and show that the renormalization to $g_{i}^{(0)}$ can be attractive for both
$\delta g_{1}^{(0)}$ and $\delta g_{3}^{(0)}$ or only $\delta g_{1}^{(0)}$,
depending upon the strength of the coupling and and phonon modes involved.

OP | Definition |  Interaction strength $(G)$
---|---|---
rCDW | $N_{\alpha}=G_{\textrm{rCDW}}\frac{\absolutevalue{\epsilon_{\alpha\beta\gamma}}}{2\mathcal{N}}\sum_{\bm{q}}\expectationvalue{c^{\dagger}_{\beta\bm{q}}c_{\gamma\bm{q}}}$ | $-2g_{1}+g_{2}-g_{3}$
iCDW | $\phi_{\alpha}=G_{\textrm{iCDW}}\frac{\epsilon_{\alpha\beta\gamma}}{2i\mathcal{N}}\sum_{\bm{q}}\expectationvalue{c^{\dagger}_{\beta\bm{q}}c_{\gamma\bm{q}}}$ | $-2g_{1}+g_{2}+g_{3}$
rSDW | $\bm{S}_{\alpha}=G_{\textrm{rSDW}}\frac{\absolutevalue{\epsilon_{\alpha\beta\gamma}}}{2\mathcal{N}}\sum_{\bm{q}}\expectationvalue{c^{\dagger}_{\beta\bm{q}}\frac{\bm{\sigma}}{2}c_{\gamma\bm{q}}}$ | $g_{2}+g_{3}$
iSDW | $\bm{\psi}_{\alpha}=G_{\textrm{iSDW}}\frac{\epsilon_{\alpha\beta\gamma}}{2i\mathcal{N}}\sum_{\bm{q}}\expectationvalue{c^{\dagger}_{\beta\bm{q}}\frac{\bm{\sigma}}{2}c_{\gamma\bm{q}}}$ | $g_{2}-g_{3}$
sSC | $\Delta_{s}=G_{\textrm{sSC}}\frac{1}{\sqrt{3}\mathcal{N}}\sum_{\bm{q}}\expectationvalue{c_{\alpha\bm{q}\downarrow}c_{\alpha-\bm{q}\uparrow}}$ | $-2g_{3}-g_{4}$
dSC | $\begin{aligned} \Delta_{xy}=&G_{\textrm{dSC}}\frac{1}{\mathcal{N}}\\\ &\times\sum_{\bm{q}}\expectationvalue{c_{\bm{q}\downarrow}D_{xy}c_{-\bm{q}\uparrow}}\\\ \Delta_{x^{2}-y^{2}}=&G_{\textrm{dSC}}\frac{1}{\mathcal{N}}\\\ &\times\sum_{\bm{q}}\expectationvalue{c_{\bm{q}\downarrow}D_{x^{2}-y^{2}}c_{-\bm{q}\uparrow}}\\\ \end{aligned}$ | $g_{3}-g_{4}$
Table 1: List of all bilinear order parameters (OPs) in the patch model.
$\bm{\sigma}$ are the Pauli matrices in spin space, and
$\sum_{\bm{q}}\equiv\sum_{\absolutevalue{\bm{q}}<\Lambda}$.
$D_{xy},D_{x^{2}-y^{2}}$ are matrices in the patch space defined as
$D_{xy}=\sqrt{1/2}\textrm{diag}\left(0,1,-1\right)$ and
$D_{x^{2}-y^{2}}=\sqrt{2/3}\textrm{ diag}(1,-1/2,-1/2)$, and
$c_{\bm{q}\sigma}=(c_{1\bm{q}\sigma},c_{2\bm{q}\sigma},c_{3\bm{q}\sigma})$.

##  IV Mean-field theory

As shown in Sec. III, we found instabilities to superconducting, charge
density wave (rCDW), orbital moment (iCDW) and spin density wave (SDW) states.
In all three AV3Sb5 materials, experiments have observed charge density wave
order (rCDW order) setting as the first instability of the symmetric state at
a $T_{c}\sim 90K$, and superconductivity only at much lower temperatures.
Hence, in this section, we neglect superconductivity, and discuss the ways
that the remaining instabilities may lead to the specific rCDW order observed
experimentally in the AV3Sb5 materials. Notably, we find that rCDW order can
be induced even if rCDW is not the primary order parameter. Hence we study the
formation of rCDWs both when the rCDW is and is not the primary instability,
and discuss the differences in the resulting properties. We carry out the
study using mean field theory.

The renormalization group analysis tells us that in phase I the rCDW is a
leading instability. In Sec. IV.2., we will study this case by calculating a
mean-field theory in the rCDW channel. In Sec. IV.3, we consider a iCDW-rCDW
coupled mean field theory, which is relevant to phases II and IV. Here we
discuss how rCDW can be induced when the iCDW is the primary order parameter.
Finally, in Sec. IV.4, we consider the situation with a leading rSDW
instability, relevant to phase III, and show how it may induce rCDW order.

###  IV.1 Complex CDW free energy

We first obtain the Landau free energy including both rCDW and iCDW order
parameters, defined as $N_{\alpha}$, $\phi_{\alpha}$ with $\alpha=1,2,3$
labeling the interpatch momentum transfer ${\bm{Q}}_{\alpha}$ (see Fig. 1
(b)), respectively. The patch-model interaction given in Eq. (3) can be
rewritten in the form of a rCDW interaction and an iCDW interaction:

$\displaystyle H_{\textrm{rCDW}}$
$\displaystyle=-\frac{\mathcal{N}G_{\textrm{rCDW}}}{2}\sum_{\alpha}\hat{\rho}_{\textrm{rC},\alpha}\hat{\rho}_{\textrm{rC},\alpha},$
$\displaystyle H_{\textrm{iCDW}}$
$\displaystyle=-\frac{\mathcal{N}G_{\textrm{iCDW}}}{2}\sum_{\alpha}\hat{\rho}_{\textrm{iC},\alpha}\hat{\rho}_{\textrm{iC},\alpha},$
(6)

where $G_{\textrm{r/iCDW}}$ is the r/iCDW interaction strength defined in Sec.
III and
$\hat{\rho}_{\textrm{rC},\alpha}=\frac{\absolutevalue{\epsilon_{\alpha\beta\gamma}}}{2\mathcal{N}}\sum_{\absolutevalue{\bm{q}}<\Lambda}c^{\dagger}_{\beta\bm{q}}c_{\gamma\bm{q}}$,
$\hat{\rho}_{\textrm{iC},\alpha}=\frac{\epsilon_{\alpha\beta\gamma}}{2\mathrm{i}\mathcal{N}}\sum_{\absolutevalue{\bm{q}}<\Lambda}c^{\dagger}_{\beta\bm{q}}c_{\gamma\bm{q}}$
is the rCDW and iCDW operators with momentum
$\bm{Q}_{\alpha}=\bm{M}_{\beta}-\bm{M}_{\gamma}$. Using the Hubbard-
Stratonovich transformation we decouple the interactions in the two channels
and then integrate out the Fermionic degrees of freedom. This gives us the
free energy as a function of the rCDW, iCDW order parameters,
$N_{1},N_{2},N_{3},\phi_{1},\phi_{2},\phi_{3}$:

$F_{\textrm{CDW}}=\frac{\mathcal{N}}{2G_{\textrm{rCDW}}}\sum_{\alpha=1}^{3}N_{\alpha}^{2}+\frac{\mathcal{N}}{2G_{\textrm{iCDW}}}\sum_{\alpha=1}^{3}\phi_{\alpha}^{2}-(\Tr\log\mathcal{G}^{-1}),$
(7)

where $\mathcal{G}^{-1}$ is defined as

$\displaystyle\mathcal{G}^{-1}(i\omega_{n},\bm{q};\\{N_{\alpha},\phi_{\alpha}\\})$
$\displaystyle=\begin{bmatrix}-i\omega_{n}+\varepsilon_{1}(\bm{q})&-(N_{3}-\mathrm{i}\phi_{3})/2&-(N_{2}+\mathrm{i}\phi_{2})/2\\\
-(N_{3}+\mathrm{i}\phi_{3})/2&-i\omega_{n}+\varepsilon_{2}(\bm{q})&-(N_{1}-\mathrm{i}\phi_{1})/2\\\
-(N_{2}-\mathrm{i}\phi_{2})/2&-(N_{1}+\mathrm{i}\phi_{1})/2&-i\omega_{n}+\varepsilon_{3}(\bm{q})\end{bmatrix}.$
(8)

Expanding the order parameter fields in Eq. (8) perturbatively, the free
energy in terms of the complex charge density order parameter
$\Delta_{\alpha}=N_{\alpha}+\mathrm{i}\phi_{\alpha}=\absolutevalue{\Delta_{\alpha}}\mathrm{e}^{\mathrm{i}\theta_{\alpha}}$
is

$\displaystyle f_{\rm CDW}=$
$\displaystyle\frac{r_{N}+r_{\phi}}{2}\sum_{\alpha}|\Delta_{\alpha}|^{2}+\frac{r_{N}-r_{\phi}}{2}\sum_{\alpha}\absolutevalue{\Delta_{\alpha}}^{2}\cos
2\theta_{\alpha}$
$\displaystyle+K_{2}\absolutevalue{\Delta_{1}}\absolutevalue{\Delta_{2}}\absolutevalue{\Delta_{3}}\cos(\theta_{1}+\theta_{2}+\theta_{3})$
$\displaystyle+K_{4}\left(\sum_{\alpha}\absolutevalue{\Delta_{\alpha}}^{2}\right)^{2}+(K_{3}-2K_{4})\sum_{\alpha<\beta}\absolutevalue{\Delta_{\alpha}}^{2}\absolutevalue{\Delta_{\beta}}^{2}$
$\displaystyle+\mathcal{O}(\Delta^{5}).$ (9)

Here,
$r_{N}=\frac{1}{2G_{\textrm{rCDW}}}+K_{1},r_{\phi}=\frac{1}{2G_{\textrm{iCDW}}}+K_{1}$,
and $K_{1},\cdots K_{4}$ are functions of temperature and chemical potential,
that we will discuss in detail in Sec. IV.2. For simplicity, we first consider
the case $K_{3}-2K_{4}<0$, which favors 3Q CDW. The configuration of
$\theta_{\alpha}$ that minimizes the free energy depends on the sign of
$r_{N}-r_{\phi}$:

  * •

When $r_{N}<r_{\phi}$, i.e. $G_{\rm rCDW}>G_{\rm iCDW}>0$, minimizing the 2nd
term requires $\theta_{\alpha}=n_{\alpha}\pi$, with $n_{\alpha}\in\mathbb{Z}$.
By choosing the proper $n_{\alpha}$, the minimization of cubic term can be
readily achieved. For any $\theta_{\alpha}=n_{\alpha}\pi$, the iCDW order
parameter must vanish, so the ground state only contains rCDW order.

  * •

When $r_{N}=r_{\phi}$, i.e. $G_{\rm rCDW}=G_{\rm iCDW}$, the second term
vanishes, and the minimization of the cubic term gives a continuously
degenerate ground state manifold. To lift the degeneracy, other perturbations
should be considered.

  * •

When $r_{N}>r_{\phi}$, i.e. $0<G_{\rm rCDW}<G_{\rm iCDW}$, the second
quadratic and cubic terms cannot be minimized simultaneously. As a result,
$\theta_{\alpha}\neq n_{\alpha}\pi/2$, with $n_{\alpha}\in\mathbb{Z}$. This
indicates that both $N_{\alpha}$ and $\phi_{\alpha}$ are non-zero. To analyze
the Free energy with iCDW as the leading instability, the iCDW-rCDW coupling
must be considered.

Below, we discuss the above three scenarios and then the rSDW-rCDW coupling
scenario.

###  IV.2 rCDW Mean-field theory

At the RG fixed point corresponding to phase I ($g_{3}^{(0)}<0$), the rCDW is
the leading density wave instability. As argued above, it is enough to
consider only the rCDW order parameters. The mean field free energy becomes

$F_{\textrm{rCDW}}=\frac{\mathcal{N}}{2G_{\textrm{rCDW}}}\sum_{\alpha=1}^{3}N_{\alpha}^{2}-(\Tr\log\mathcal{G}^{-1}),$
(10)

where $\mathcal{G}^{-1}$ is defined as

$\displaystyle\mathcal{G}^{-1}(i\omega_{n},\bm{q};\\{N_{\alpha}\\})$
$\displaystyle=\begin{bmatrix}-i\omega_{n}+\varepsilon_{1}(\bm{q})&-N_{3}/2&-N_{2}/2\\\
-N_{3}/2&-i\omega_{n}+\varepsilon_{2}(\bm{q})&-N_{1}/2\\\
-N_{2}/2&-N_{1}/2&-i\omega_{n}+\varepsilon_{3}(\bm{q})\end{bmatrix}.$ (11)

As a function of $N_{\alpha}$, the rCDW free energy has the symmetry of the
tetrahedral point group, $T_{d}$. This can be deduced by applying the symmetry
operations of the full space group of the lattice, $P6/mmm$, to the
definitions of the rCDW order parameter. Using our knowledge of this symmetry,
we can determine what solutions of the free energy are possible. In the cases
that are physically relevant, the solutions of the mean-field theory belong in
either the 3Q$+$, 3Q$-$, or 1Q rCDW configurations which are classes of rCDW
states with directions:

3Q$+$: |  $\\{(111),(1\bar{1}\bar{1}),(\bar{1}1\bar{1}),(\bar{1}\bar{1}1)\\}$,
---|---
3Q$-$: |  $\\{(\bar{1}\bar{1}\bar{1}),(\bar{1}11),(1\bar{1}1),(11\bar{1})\\}$,
1Q: |  $\\{(100),(\bar{1}00),(010),(0\bar{1}0),(001),(00\bar{1})\\}$.

Indeed, by numerically solving the full free energy we will see that the
solutions belong in either the 3Q$\pm$ or 1Q rCDW classes. A detailed
discussion is provided in App. C. The rCDW patterns of the 3Q$\pm$ and 1Q
states are shown in Fig. 7.

By assuming we are near the rCDW transition where the rCDW order parameters
are sufficiently small, we can expand the free energy to fourth order in
$N_{\alpha}$. The resulting Landau theory is

$\displaystyle
f_{\textrm{rCDW}}=\left(\frac{1}{2G_{\textrm{rCDW}}}+K_{1}\right)\sum_{\alpha}N_{\alpha}^{2}+K_{2}N_{1}N_{2}N_{3}$
$\displaystyle\quad+K_{4}\left(\sum_{\alpha}N_{\alpha}^{2}\right)^{2}+(K_{3}-2K_{4})\sum_{\alpha<\beta}N_{\alpha}^{2}N_{\beta}^{2}+\mathcal{O}(N^{5}),$
(12)

where $f_{\textrm{rCDW}}$ is the free energy density. The definitions of the
coefficients $K_{1},\cdots,K_{4}$ which are functions of chemical potential
and temperature are given in App. D. The coefficients can be evaluated
asymptotically in the limit $\mu,k_{B}T\ll t\Lambda^{2}$ and can also be found
in App. D.

Fig. 5 shows the temperature dependence of these coefficients when
$\mu/t\Lambda^{2}=0.01$. Notice that the coefficients exhibit a change in
behavior at the cross-over temperature $k_{B}T\sim\mu$. In addition, the
$K_{2},K_{3},K_{4}$ coefficients change sign near the cross-over temperature.
When the quartic coefficients become negative, the Landau theory expression
given by Eq. (12) becomes unstable, but in those regions stability can be
restored by including sixth-order terms to the free energy.

Figure 5: Absolute value of the Landau theory coefficients,
$K_{1},\cdots,K_{4}$ evaluated at $\mu/t\Lambda^{2}=10^{-2}$ for the case of
perfect nesting ($a=9t/4\,,b=3t/4$). Solid lines indicate positive values and
dashed lines indicate negative values. $K_{2},K_{3},K_{4}$ change sign at
$\mu/k_{B}T\sim 2.14,4.05,1.91$ respectively.

Eq. (12) has a third order term, $K_{2}$, that couples all three rCDW order
parameters. This is allowed by symmetry since this term is even under time-
reversal symmetry and the sum of the three nesting wave vectors satisfies
$\sum_{\alpha}{\bm{Q}}_{\alpha}\equiv 0$. This term introduces a preference
for 3Q$+$ or 3Q$-$ rCDW states depending on the sign of $K_{2}$. On the other
hand, when $K_{3}-2K_{4}>0$, the fourth order term sets a preference for 1Q
rCDW states. When $K_{3}-2K_{4}<0$, this fourth order term prefers the 3Q$\pm$
states equally.

When $a/b=3$, Fig. 5 shows that $K_{2}<0$ when $\mu/k_{B}T\lesssim 2.14$ in
the asymptotic limit. In addition, $K_{3}-2K_{4}<0$ for $\mu\ll k_{B}T$.
Hence, in this region, all terms in the rCDW Landau theory prefer the 3Q$+$
rCDW state. The transition to the 3Q$+$ rCDW state must be a first-order
transition because the free energy can become negative before the second-order
term vanishes due to the third-order term. One thing to note is that the rCDW
order parameter is not necessarily small near the first-order rCDW transition,
so the results of the Landau theory must be treated with caution. We can avoid
this issue by numerically solving the full free energy.

The full free energy defined in Eqs. (10), (11) has four tunable parameters:
temperature $T$, chemical potential $\mu$, the rCDW interaction strength
$G_{\textrm{rCDW}}$, and the nesting ratio $a/b$. As explained in Sec. II, we
only need to consider the case $a,b>0$. Given this condition, we can introduce
a convenient reparametrization,

$\displaystyle a=$
$\displaystyle\frac{\delta+\sqrt{3+\delta^{2}}}{\sqrt{3}}t,$ $\displaystyle
b=$ $\displaystyle\frac{-\delta+\sqrt{3+\delta^{2}}}{\sqrt{3}}t,$ (13)

where $t>0\textrm{ and }\delta\in\mathbb{R}$. We see that $t$ represents the
bandwidth of the saddle point band and $\delta$ represents the degree of
nesting. Under this parametrization, $a/b\geq 1$ for $\delta\geq 0$ and
$a/b<1$ for $\delta<0$. In particular, $a/b=3$ (perfect nesting) for
$\delta=1$.

For different values of $G_{\textrm{rCDW}}$ and $\delta$, we generate a phase
diagram in the $T-\mu$ plane. The results shown in Fig. 6 are representative
samples of the the entire parameter space where the continuum model
holds—$k_{B}T\,,\mu\,,G_{\textrm{rCDW}}\ll t\Lambda^{2}$. We see that when the
system is doped to the saddle point ($\mu=0$) or below it, only the 3Q$+$ rCDW
forms. As you dope the system above the saddle point, 1Q and/or 3Q$-$ rCDW
regions can emerge. In addition, for $\delta=1.6$, the 1Q rCDW region is small
and vanishes for sufficiently small $G_{\textrm{rCDW}}$ as seen in Figs. 6, 6.
Similarly, at perfect nesting and more generally $\delta\sim 1$, the 3Q$-$
rCDW region vanishes for sufficiently small $G_{\textrm{rCDW}}$ as seen in
Fig. 6. Clearly, the 3Q$+$ rCDW region is largest in all phase diagrams which
is consistent with the Landau theory discussed above.

Figure 6: rCDW phase diagrams calculated in the $k_{B}T$–$\mu$ plane.
$\tilde{G}_{\textrm{rCDW}}=G_{\textrm{rCDW}}/t\Lambda^{2}$ is defined as the
dimensionless rCDW strength. $\delta$ parametrizes nesting as defined in Eq.
(13). The phase diagrams correspond to the case $t,\delta>0$ ($a,b>0\textrm{
and }a/b>1$). The phase diagrams for $t,\delta<0$ ($a,b>0\textrm{ and }a/b<1$)
are related to these phase diagrams under the exchange
$\textrm{3Q}+\leftrightarrow\textrm{3Q}-$ and reflection in the $\mu=0$ axis.

###  IV.3 iCDW-rCDW mean field theory

Now, we consider how a fundamental iCDW order can induce subsidiary rCDW
order. To do so, we consider a iCDW-rCDW coupled mean field theory. We are
particularly interested in the regime, motivated by the RG results for phases
II and IV, in which iCDW and rCDW are close in energy and can compete.

It is convenient to express Eq. (9) in terms $N_{\alpha}$ and $\phi_{\alpha}$:

$\displaystyle f_{\rm
CDW}=r_{\phi}\sum_{\alpha=1}^{3}\phi_{\alpha}^{2}+r_{N}\sum_{\alpha=1}^{3}N_{\alpha}^{2}$
$\displaystyle+K_{2}\left(N_{1}N_{2}N_{3}-\phi_{2}\phi_{3}N_{1}-\phi_{1}\phi_{3}N_{2}-\phi_{2}\phi_{1}N_{3}\right)$
$\displaystyle+K_{4}\left(\sum_{\alpha=1}^{3}\phi_{\alpha}^{2}+N_{\alpha}^{2}\right)^{2}$
$\displaystyle+(K_{3}-2K_{4})\sum_{\alpha<\beta}(\phi_{\alpha}^{2}+N_{\alpha}^{2})(\phi_{\beta}^{2}+N_{\beta}^{2})+\mathcal{O}(\phi^{6},N^{5}).$
(14)

Recall that $r_{\phi}=\left(\frac{1}{2G_{\textrm{iCDW}}}+K_{1}\right)$,
$r_{N}=\left(\frac{1}{2G_{\textrm{rCDW}}}+K_{1}\right)$. We assume $G_{\rm
iCDW}>G_{\rm rCDW}>0$, so that iCDW is the leading instability. Note that
time-reversal symmetry forbids the term cubic in $\phi$, i.e.
$\phi_{1}\phi_{2}\phi_{3}$. As a result, the order parameter manifold for
$\phi_{\alpha}$ alone has cubic symmetry $O_{h}$. In addition, the 3Q$\pm$
classes which were distinct for the case of rCDW are now related under time-
reversal symmetry and part of the larger 3Q class which is defined as the
union of the 3Q$\pm$ classes. Close to and below the transition temperature
for iCDW, the iCDW-rCDW coupling can be treated as a perturbation.

First, we consider the iCDW Landau theory. The order parameter manifold for
$\phi_{\alpha}$ is determined by quartic terms proportional to $K_{3}$ and
$K_{4}$. The first term $K_{4}(\sum_{\alpha=1}^{3}\phi^{2}_{\alpha})^{2}$ is
positive definite and isotropic, so the selection of the ground state
configuration is determined by the term
$(K_{3}-2K_{4})\sum_{\alpha<\beta}\phi_{\alpha}^{2}\phi_{\beta}^{2}$.
Specifically, (i) when $K_{3}>2K_{4}$, the 1Q iCDW state is favored, (ii) when
$K_{3}=2K_{4}$, the $\phi_{\alpha}$ ground state configuration is degenerate
at quartic order, and requires higher order terms to break the degeneracy,
(iii) when $K_{3}<2K_{4}$, the 3Q iCDW is favored. From the asymptotic
solution shown in Fig. 5 at $\mu/t\Lambda^{2}=0.01$, we see that when
$k_{B}T_{c}/t\Lambda^{2}\lesssim 0.01$, and the 1Q iCDW is stable. When
$k_{B}T/t\Lambda^{2}\gtrsim 0.01$, $K_{3}-2K_{4}<0$. Thus 3Q iCDW develops if
$T_{c}$ sits in this range of temperature, but it may become unstable to 1Q
iCDW as temperature lowers. On the other hand, when $\mu\ll k_{B}T$, we find
$K_{3}-2K_{4}<0$ in the low temperature regime when
$k_{B}T/t\Lambda^{2}<0.025$ (see App. D). Thus 3Q iCDW should be stable at
$\mu\sim 0$. This transition is a continuous phase transition in contrast to
the rCDW case because of the absence of a third-order term. The real space
pattern for 1Q and 3Q iCDW orders are presented in Sec. V.1.2.

Next, we substitute the iCDW order parameter in Eq. (14) with the iCDW saddle
point solution, $|\phi_{\alpha}|=\phi^{*}$, and obtain the free energy for
rCDW OPs. For 1Q iCDW, no rCDW can be induced. For 3Q iCDW, as the cubic term
$K_{2}$ contains term linear in $N_{\alpha}$ through coupling to two iCDW OPs
at another two momenta, considering up to quadratic term in $N$, we find
$\\{N_{1},N_{2},N_{3}\\}=N_{0}\\{1,\pm 1,\pm 1\\}$ up to any permutations of
the signs between $N_{\alpha}$, where
$N_{0}\sim\frac{K_{2}\phi^{*2}}{(T-T_{\rm rCDW})}$. Here, $T-T_{\rm
rCDW}\sim\left(\frac{1}{2G_{\rm rCDW}}+K_{1}\right)>0$. So the sign of $N_{0}$
is determined by the sign of $K_{2}$. Since $K_{2}<0$ for $\mu/T\lesssim
2.14$, $N_{0}<0$, the induced rCDW order must be the 3Q$-$ one. As temperature
further lowers to $T<T_{\rm rCDW}$, the quartic term in $N_{\alpha}$ should be
included so that the free energy is stable. We checked that the rCDW order
remains the 3Q$-$ one for $K_{2}<0$ and vice versa.

###  IV.4 rSDW-rCDW Mean-field theory

In phase III of the RG phase diagram, the rSDW is the leading density wave
instability. Here, we consider how rCDW order can emerge as a subsidiary order
through rSDW-rCDW coupling. By rewriting the full interaction, Eq. (3), using
the rSDW operators the rSDW interaction term is

$H_{\textrm{rSDW}}=-\frac{\mathcal{N}G_{\textrm{rSDW}}}{2}\sum_{\alpha}\hat{\bm{\rho}}_{\textrm{rS},\alpha}\cdot\hat{\bm{\rho}}_{\textrm{rS},\alpha},$
(15)

where
$\hat{\bm{\rho}}_{\textrm{rS},\alpha}=\frac{\absolutevalue{\epsilon_{\alpha\beta\gamma}}}{2\mathcal{N}}\sum_{\absolutevalue{\bm{q}}<\Lambda}c^{\dagger}_{\beta\bm{q}}\frac{\bm{\sigma}}{2}c_{\gamma\bm{q}}$
is the spin density operator. Assuming the interaction in rCDW channel is also
attractive, but much weaker than that of rSDW, we can decouple the
interactions in the rSDW and rCDW channel using a Hubbard-Stratonovich
transformation, integrate out the Fermionic degrees of freedom, and expand the
resulting mean field free energy to fourth order in the rSDW and rCDW order
parameters. This gives us the rSDW free energy (which may be added to the rCDW
free energy in Eq. (12)):

$\displaystyle f_{\textrm{rSDW}}=$
$\displaystyle\left(\frac{1}{2G_{\textrm{rSDW}}}+\frac{K_{1}}{4}\right)\sum_{\alpha}\absolutevalue{\bm{S}_{\alpha}}^{2}$
$\displaystyle+\frac{K_{2}}{4}\left(N_{1}\bm{S}_{2}\cdot\bm{S}_{3}+N_{2}\bm{S}_{3}\cdot\bm{S}_{1}+N_{3}\bm{S}_{1}\cdot\bm{S}_{2}\right)$
$\displaystyle+\frac{K_{3}}{16}\sum_{\alpha<\beta}\absolutevalue{\bm{S}_{\alpha}}^{2}\absolutevalue{\bm{S}_{\beta}}^{2}+\frac{K_{4}}{16}\sum_{\alpha}\absolutevalue{\bm{S}_{\alpha}}^{4}+$
$\displaystyle+K_{5}\sum_{\alpha<\beta}\absolutevalue{\bm{S}_{\alpha}\cdot\bm{S}_{\beta}}^{2}+K_{6}\left(\bm{S}_{1}\cdot\bm{S}_{2}\times\bm{S}_{3}\right)^{2}+\cdots.$
(16)

Here $K_{5},K_{6}$ are new symmetry-allowed coefficients. The important aspect
of this Landau theory is the third-order term which couples the rSDW and rCDW
order parameters. This term is allowed by symmetry since it is invariant under
time-reversal and lattice translation symmetry.

Close to the transition temperature $T_{c}$ for the rSDW order, the rSDW-rCDW
coupling term can be treated as a perturbation and we can first solve the rSDW
Landau theory. Like the iCDW OPs, the 3Q$+$ and 3Q$-$ classes are related
under time-reversal symmetry so we only need to consider the 3Q and 1Q rSDW
classes. This Landau theory has previously been studied in Ref. [47]. For
chemical potential sufficiently close to the saddle point, the solution to the
SDW Landau theory is a uniaxial 3Q rSDW phase where all three rSDW order
parameters have the same magnitude and are oriented along the same axis, i.e.
$\bm{S}_{\alpha}^{*}=\bm{s}n_{\alpha}$, where $n_{\alpha}=\pm 1$. If we
substitute this solution into Eq. (16) and add the terms from Eq. (12), we get
a free energy that is a function of just the rCDW OPs:

$\displaystyle
f_{\textrm{rCDW}}+f_{\textrm{rSDW}}\big{\rvert}_{\bm{S}_{\alpha}=\bm{S}_{\alpha}^{*}}=\left(\frac{1}{2G_{\textrm{rCDW}}}+K_{1}\right)\sum_{\alpha}N_{\alpha}^{2}$
$\displaystyle\quad+\frac{K_{2}}{4}\absolutevalue{\bm{s}}^{2}\left(n_{2}n_{3}N_{1}+\textrm{cyc.
perms.}\right)+K_{2}N_{1}N_{2}N_{3}$
$\displaystyle\quad+K_{4}\left(\sum_{\alpha}N_{\alpha}^{2}\right)^{2}+(K_{3}-2K_{4})\sum_{\alpha<\beta}N_{\alpha}^{2}N_{\beta}^{2}.$
(17)

For $K_{3}-2K_{4}<0$, all terms prefer the 3Q$+$ rCDW state, so the system can
develop a rCDW state.

##  V Extensions and experimental implications

In this section, we discuss various aspects of the different CDW phases, and
how they may be differentiated experimentally.

###  V.1 Real space r/iCDW patterns

It is interesting to consider the real space patterns of charges and currents
associated with the rCDW and iCDW order parameters.

Our continuum model does not carry any details of the lattice, so we rely on
symmetry information provided by DFT calculations [58]. Importantly, as the
saddle point band at the three M points are even under inversion, it can be
shown that among the d-orbitals of vanadium atoms, only the $\alpha$-th
vanadium atom in the unit cell contributes to the Bloch state at
${\bm{M}}_{\alpha}$ (see App. A). Furthermore, the DFT calculation shows that
the saddle points consist mostly of d-orbitals of vanadium atoms. Motivated by
these observations, we consider nearest neighbor hopping on the Kagomé lattice
as the minimal tight binding model that should capture the essential physics
from fermions near the saddle points. For convenience, the lattice coordinate
for the three sublattice is expressed as
${\bm{r}}_{\alpha}={\bm{R}}+{\bm{\delta}}_{\alpha}$, where ${\bm{R}}$ is the
coordinate for a unit cell, whose origin is taken at the center of the
triangular plaquette with the green sublattice facing to the left in Fig. 1.
${\bm{\delta}}_{1}=(-\frac{1}{2\sqrt{3}},0)$,
${\bm{\delta}}_{2}=(\frac{1}{4\sqrt{3}},\frac{1}{4})$ and
${\bm{\delta}}_{3}=(\frac{1}{4\sqrt{3}},-\frac{1}{4})$. The tight-binding
Hamiltonian reads

$\displaystyle\hat{H}_{0}=t_{1}\sum_{\langle{\bm{r}}_{\alpha},{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}\rangle}d^{\dagger}_{{\bm{r}}_{\alpha}}d_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}+d^{\dagger}_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}d_{{\bm{r}}_{\alpha}},$
(18)

where $\langle{\bm{r}}_{\alpha},{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}\rangle$ is
the NN bond and ${\bm{e}}$ is defined such that
${\bm{r}}_{\gamma}={\bm{r}}_{\alpha}+{\bm{e}}_{\beta}$ where
$\\{\beta,\alpha,\gamma\\}$ is a permutation of $\\{1,2,3\\}$. This gives
${\bm{e}}_{1}=\frac{1}{2}\\{-1,0\\}$,
${\bm{e}}_{2}=\frac{1}{2}\\{-\frac{\sqrt{3}}{2},\frac{1}{2}\\}$ and
${\bm{e}}_{3}=\frac{1}{2}\\{\frac{\sqrt{3}}{2},\frac{1}{2}\\}$.

Due to the one-to-one correspondence between patch label ${\bm{M}}_{\alpha}$
and sublattice label $\rm{V}_{\alpha}$ of vanadium atoms, the real space
fermions on vanadium atoms can be expressed as

$\displaystyle
d_{{\bm{r}}_{\alpha}}=\frac{1}{\sqrt{\mathcal{N}}}\sum_{\bm{k}}\mathrm{e}^{\mathrm{i}{\bm{k}}\cdot{\bm{R}}}d_{\alpha,{\bm{k}}}\approx\frac{1}{\sqrt{\mathcal{N}}}\sum_{|{\bm{q}}|<\Lambda}\mathrm{e}^{\mathrm{i}{\bm{M}}_{\alpha}\cdot{\bm{R}}}c_{\alpha{\bm{q}}},$
(19)

where $c_{\alpha{\bm{q}}}$ denotes the saddle point fermion near
${\bm{M}}_{\alpha}$ as defined in the continuous model (Eq. (1)). As both rCDW
and iCDW order parameters are condensates of inter-patch fermion density
operators, the associated real space order must be a _bond order_ on the
Kagomé lattice. For simplicity, we will consider only nearest neighbor (NN)
bond order as an example.

####  V.1.1 rCDW pattern

The rCDW order parameter is time-reversal even, and contribute to bond density
modulation as

$\displaystyle\langle\delta\hat{\rho}_{{\bm{r}}_{\alpha}{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}\rangle$
$\displaystyle=\langle\delta\hat{\rho}_{{\bm{r}}_{\alpha}{\bm{r}}^{\prime}_{\gamma}}\rangle=\operatorname{Re}\left[\langle
d^{\dagger}_{{\bm{r}}_{\alpha}}d_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}\rangle\right]$
$\displaystyle=\frac{N_{\alpha\gamma}}{G_{\rm
rCDW}}\cos({\bm{M}}_{\alpha}\cdot{\bm{R}}-{\bm{M}}_{\gamma}\cdot{\bm{R}}^{\prime}),$
(20)

where
${\bm{r}}^{\prime}_{\gamma}={\bm{r}}_{\alpha}+{\bm{e}}_{\beta}={\bm{R}}^{\prime}+{\bm{\delta}}_{\gamma}$
with $\gamma\neq\alpha\neq\beta$. In Fig. 7, the rCDW charge bond density
modulation is shown for the 3Q$\pm$ and 1Q rCDW states.

Figure 7:  Real space rCDW bond order pattern of the 3Q$\pm$ and 1Q
configurations. Configurations in the 3Q$\pm$ classes are related by
translation. Configurations in the 1Q class are related by translation and a
three-fold rotation. 7 Bond ordering corresponding to the 3Q$+$ rCDW class.
The bond ordering forms hexagonal and triangular plaquettes. 7 Bond ordering
corresponding to the 3Q$-$ rCDW class. The bond ordering forms a ‘star-of-
David’ pattern. 7 Bond ordering pattern corresponding to the 1Q rCDW class.
The pattern has alternating bond strength along one direction and uniform
strength in the other two bond directions.

####  V.1.2 iCDW pattern

Figure 8: Real space iCDW (current) bond order pattern of the (a) 1Q and (b)
3Q states. All 3Q configurations are related by translations or time-reversal
(changing the sign of the currents). Configurations in the 1Q class are
related by translation (i.e. reversing the current direction) and a three-fold
rotation. Note that the 1Q state is macroscopically time-reversal symmetric,
because it is invariant under the combination of time-reversal and
translation. (a) Bond current in the 1Q class, with
$\\{\phi_{23},\phi_{31},\phi_{12}\\}=\phi_{0}\\{1,0,0\\}$. (b) Bond current
from $\\{\phi_{23},\phi_{31},\phi_{12}\\}=\phi_{0}\\{1,1,1\\}$. The sign of
$\phi_{0}$ is taken as negative here.

The iCDW order parameter breaks time-reversal symmetry and corresponds to bond
current in real space. Note, that the current operator is well defined only
when the charge is conserved, which is indeed the case in both the high
temperature disordered phase and iCDW phase. As a result, in the equilibrium
phase, the current operator must satisfy

$\displaystyle\sum_{{\bm{r}}_{\alpha}\rightarrow{\bm{r}}^{\prime}_{\alpha^{\prime}}}\hat{j}_{{\bm{r}}_{\alpha}{\bm{r}}^{\prime}_{\alpha^{\prime}}}=0.$
(21)

By comparing the continuity equation and the equation of motion for charge
density, we find the current operator on the NN bonds as (see App. E):

$\displaystyle\hat{j}_{{\bm{r}}_{\alpha}{\bm{r}}^{\prime}_{\alpha^{\prime}}}$
$\displaystyle=\hat{j}_{{\bm{r}}_{\alpha}{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}\delta_{{\bm{r}}^{\prime}_{\alpha^{\prime}},{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}\delta_{\alpha\neq\beta}$
$\displaystyle=\frac{\mathrm{i}e}{\hbar}\,t_{1}\left(d^{\dagger}_{{\bm{r}}_{\alpha}}d_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}-d^{\dagger}_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}d_{{\bm{r}}_{\alpha}}\right),$
(22)

where $t_{1}$ is the NN hopping defined in Eq. (18). From Eq. (19), the bond
current expectation value can be expressed in terms of iCDW order parameter as

$\displaystyle\langle\hat{j}_{{\bm{r}}_{\alpha}{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}\rangle=\langle\hat{j}_{{\bm{r}}_{\alpha}{\bm{r}}^{\prime}_{\gamma}}\rangle$
$\displaystyle=$
$\displaystyle\mathrm{i}\,e\,t_{1}\,\mathrm{e}^{-\mathrm{i}({\bm{M}}_{\alpha}\cdot{\bm{R}}-{\bm{M}}_{\gamma}\cdot{\bm{R}}^{\prime})}\left(\frac{1}{\mathcal{N}}\sum_{{\bm{q}}}\left(\langle
c^{\dagger}_{\alpha,{\bm{q}}}c_{\gamma,{\bm{q}}}\rangle-\langle
c^{\dagger}_{\gamma,{\bm{q}}}c_{\alpha,{\bm{q}}}\rangle\right)\right)$
$\displaystyle=$ $\displaystyle\,\frac{-2\,e\,t_{1}}{G_{\rm
iCDW}}\mathrm{e}^{-\mathrm{i}({\bm{M}}_{\alpha}\cdot{\bm{R}}-{\bm{M}}_{\gamma}\cdot{\bm{R}}^{\prime})}\phi_{\alpha\gamma},$
(23)

where again we used
${\bm{r}}^{\prime}_{\gamma}={\bm{r}}_{\alpha}+{\bm{e}}_{\beta}={\bm{R}}^{\prime}+{\bm{\delta}}_{\gamma}$
with $\gamma\neq\alpha\neq\beta$. From Eq. (23), the real space bond current
pattern from $\phi_{23}$ are shown in Fig. 8 (a). The linear combinations of
them can form loop current, as an example, we show the bond current pattern
for $\\{\phi_{23},\phi_{31},\phi_{12}\\}=\phi_{0}\\{1,1,1\\}$ in Fig. 8 (b).

The primary order parameter for the iCDW is a loop current. However, as was
discussed in Sec. IV.3, from the form of the $K_{2}$ coupling terms in Eq.
(14), a 3Q iCDW will also induce charge order. In particular we see that a 3Q
iCDW induces either a 3Q$+$ or 3Q$-$ state depending upon the sign of $K_{2}$.
Notably, the charge order is quadratic in the iCDW order parameter,
$N_{\alpha}\sim\phi^{2}$, which could be detectable near the transition
temperature.

###  V.2 Three-dimensional coupling

Here we extend the Landau theory to consider the implications of coupling of
CDW order parameters between nearby layers, adding a layer index
$z=0,1,2\cdots$ numbered beginning from the top layer. We assume this coupling
is weak, so can be approximated by the leading terms linear in the order
parameters in each layer, and decays rapidly with the distance between layers.
Hence,

$f_{\perp}=\sum_{z=0}^{\infty}\sum_{\delta=1}^{\infty}\left(K_{\perp,\delta}N_{\alpha,z}N_{\alpha,z+\delta}+L_{\perp,\delta}\phi_{\alpha,z}\phi_{\alpha,z+\delta}\right).$
(24)

We assume that all the inter-layer interactions are weak compared to the
intra-layer terms in the free energy, so that the form of the order within
each layer is established by the latter, and the intra-layer terms serve to
select particular relative orientations of the different symmetry-breaking
states in nearby layers. Here we expect $|K_{\perp,1}|\gg|K_{\perp,\delta>1}|$
and $|L_{\perp,1}|\gg|L_{\perp,\delta>1}|$ so that terms with $\delta>1$ can
be neglected unless they are required to break degeneracies.

####  V.2.1 Three-dimensional ground states

Now we discuss the resulting three-dimensional ordered structures. First
consider the imaginary CDW, within the 3Q phase. Within a given layer,
$\phi_{\alpha}$ may take one of the values $\phi_{\alpha}=|\phi_{0}|(\pm 1,\pm
1,\pm 1)$, where all 8 signs are possible, and $|\phi_{0}|$ is fixed by
single-layer energetics. If $L_{\perp,1}<0$, the minimum energy iCDW order
parameter is identical in all layers, $\phi_{\alpha,z}=\phi_{\alpha}$. If
instead, $L_{\perp,1}>0$, the minimum energy configuration is
“antiferromagnetic”, $\phi_{\alpha,z}=(-1)^{z}\phi_{\alpha}$. The two cases
above correspond to an ordering wavevector with the z-component
$k_{z}=0,\frac{1}{2}$, in lattice units. Note that $k_{z}=1/2$ is the
wavevector for the current order, but the induced real CDW order would have
$q_{z}=0$ in both cases.

Now consider the real CDW, in either the 3Q$+$ or 3Q$-$ states. Within a
single layer, $N_{\alpha}$ may take just four values, with
$N_{\alpha}=N_{0}n_{\alpha}$, with $n_{\alpha}^{2}=1$ and $n_{1}n_{2}n_{3}=1$.
For $K_{\perp,1}<0$, we again obtain a “ferromagnetic” state, with
$N_{\alpha,z}=N_{\alpha}$. For $K_{\perp,1}>0$, however, the situation is
distinct from the iCDW case, because an overall sign change in $N_{\alpha}$ is
not permitted. Instead, for a given state $n_{\alpha,z}$, the inter-layer
coupling equally favors $n_{\alpha,z+1}$ in any of the three states not equal
to $n_{\alpha,z}$. For $N$ layers, the total degeneracy of ground states,
consider just the $K_{\perp,1}$ interaction, is $4\times 3^{N-1}$ , a
macroscopic degeneracy including states with arbitrary wavevectors $k_{z}$. We
must therefore consider further neighbor interactions. The problem can be
mapped to a 4-state Potts model, by defining the four allowed configurations
of $N_{\alpha}$ in a single layer as $\sigma=1\cdots 4$. The interaction
energy becomes

$f_{\perp}=\sum_{z=0}^{\infty}\sum_{\delta=1}^{\infty}K_{\perp,\delta}N_{0}^{2}\left(4\delta_{\sigma_{z},\sigma_{z+\delta}}-1\right),$
(25)

which illustrates the $S_{4}$ permutation symmetry of a Potts model. This is a
4-state Potts chain with competing further-neighbor interactions, which has a
rich statistical mechanics for general couplings, similar to that of the ANNNI
model, a paradigm for devil’s staircases, commensurate and incommensurate
phases, and transitions between them. Here because each Potts spin represents
an entire 2d layer, the energies involved are proportional to the area of a
layer, and hence much larger than $k_{B}T$. Therefore, we are interested only
in the ground states of the Potts chain. In this limit, the ground states are
generally commensurate, but can have very large unit cells (in the $z$
direction), and the Devil’s staircase can arise. We limit our discussion to
only the simplest cases, and assume
$|K_{\perp,\delta}|\gg|K_{\perp,\delta+1}|$ as expected on grounds of
locality.

The simplest situation is $K_{\perp,2}<0$, in which case second neighbor
layers prefer to be parallel. We have then alternating states in successive
layers, $\sigma_{z}=\sigma_{1}$ for $z$ even and $\sigma_{z}=\sigma_{2}$ for
$z$ odd, with $\sigma_{1}\neq\sigma_{2}$. In terms of the rCDW vector,,
$n_{\alpha,z}=n^{(1)}_{\alpha}$ for $z$ even and
$n_{\alpha,z}=n^{(2)}_{\alpha}$ for $z$ odd, such that
$\bm{n}^{(1)}\cdot\bm{n}^{(2)}=-1$. This corresponds to the wavevector
$k_{z}=1/2$ in lattice units. Note that the three-dimensional degeneracy of
this state is $4\times 3=12$.

If $K_{\perp,2}>0$, then we require both nearest neighbor and second neighbor
Potts spins to differ in the ground state. After choosing the first layer,
there are 3 choices for the second layer, and then 2 choices for the third
layer, which must be distinct from the first two layers. This implies the
smallest possible periodicity of the ground state is 3. It may, however, be
larger. Indeed for the fourth layer, there are still 2 choices remaining, and
the ground state is not determined. To fix this degeneracy, we may yet
consider the third neighbor coupling. If $K_{\perp,3}<0$, then we favor a
return to the original state. The entire configuration becomes determined,
with a periodicity of 3, and a representative sequence in Potts variables like
$\sigma_{z}=1,2,3,1,2,3,1,2,3\cdots$. These configurations correspond to a
“chiral” ordering of the layers, consistent with a 3-fold screw axis.
Alternatively, if $K_{\perp,3}>0$, then we obtain a four layer periodicity,
$\sigma_{z}=1,2,3,4,1,2,3,4,\cdots$. These two situations have smaller
ordering wavevectors in the $z$ direction, $k_{z}=1/3,1/4$, respectively. Note
that the relatively simple results quoted here are the result of assuming a
strict hierarchy of interactions, with couplings decaying rapidly in strength
with the separation of layers. Otherwise much more complex states may arise in
the Potts chain.

We conclude that a distinct difference between the iCDW and rCDW is the
presence of periodicities larger than 2 in the $c$ direction. Appearance of
such a periodicity (i.e. $k_{z}=1/3,1/4$) would provide clear experimental
evidence in favor of the real over the imaginary CDW.

####  V.2.2 Rotational symmetry breaking

In the previous discussion, we determined the relative ordering between layers
assuming the order parameter within each layer is rigid. Now we consider a
higher order effect: the back-influence of the inter-layer interaction on the
order parameter within a single layer. In particular, in the case of rCDW
order, this leads to a breaking of $C_{3}$ rotational symmetry within a given
layer. This can be understood as follows. All 4 of the 2d ordered states in
the 3Q$+$ or 3Q$-$ states preserve $C_{3}$ symmetry around the centers of one
of the four hexagons within the quadrupled unit cell, but not around the other
three. In the states with $k_{z}>0$, the centers of neighboring layers are not
aligned. Consequently there is no rotation axis which preserves all layers. In
the chiral $k_{z}=1/3$ states, there is instead a screw axis, which preserves
macroscopic $C_{3}$ symmetry of the crystal in the bulk. This symmetry is
however broken at the surface. In the states with $k_{z}=1/2,1/4$, there is
not even macroscopic $C_{3}$ symmetry.

In all cases, if one observes the order within a single layer, its neighbors
will influence its order and lower the symmetry. The situation is simplest for
the top layer. Assume the system develops long-range order, and therefore we
may treat the inter-layer interaction in a mean field sense. We therefore
replace the coupling to the second layer from the top (the strongest such
coupling) by a term of the form

$f_{\perp,0}=K_{\perp,1}\langle N_{\alpha,1}\rangle N_{\alpha,0}.$ (26)

This term appears as a “field” on the order parameter in the top layer
($z=0$). In for example the $k_{z}=1/2$ phase, we may take $\langle
N_{\alpha,1}\rangle=N_{0}(1,1,1)$. Then this configurations in the first layer
are “pushed” away from the $(1,1,1)$ direction. For example, if the top layer
chooses the $(1,\bar{1},\bar{1})$ state, the inter-layer coupling will shift
it to the form $N_{\alpha,0}=N_{0}(1-\delta,-1-\delta,-1-\delta)$, with
$0<\delta\ll 1$. Note that $|N_{1}|<|N_{2}|=|N_{3}|$. The consequence is that
the two dimensional Bragg peaks associated to density oscillations in the top
layer develop two unequal magnitudes. This is a sign of the rotational
symmetry breaking. Indeed, one can define a two-dimensional vector
$\bm{v}_{z,z^{\prime}}=-\bm{v}_{z^{\prime},z}$

$\bm{v}_{z,z^{\prime}}=\sum_{\alpha}\langle
n_{\alpha,z}n_{\alpha,z^{\prime}}\rangle\bm{a}_{\alpha},$ (27)

where $\bm{a}_{\alpha}$ are the 2d triangular Bravais lattice vectors with
$\bm{a}_{3}=-\bm{a}_{1}-\bm{a}_{2}$. The vector $\bm{v}_{z,z^{\prime}}$ is
oriented along one of the three principle directions and selects this axis.

This rotational symmetry breaking effect exists within each layer in all the
3Q$\pm$ rCDW phases except the uniform $k_{z}=0$ one. Rotational symmetry is
also broken in the bulk (i.e. in an infinite system in the $z$ direction)
except when it is restored macroscopically by an arrangement of layers that
constitutes a screw axis. The latter occurs only for the $k_{z}=1/3$ case
detailed above. Rotational symmetry breaking is, however, absent both in
individual layers and in bulk, in the 3Q iCDW states, providing another means
to differentiate the iCDW from the rCDW experimentally.

###  V.3 Critical behavior

We briefly discuss the expected critical behavior at the ordering temperature
for some important cases based on the symmetries and order parameters. This
question is motivated by the presence of a cubic term in the rCDW Landau
theory, which might suggest a first order transition to the CDW state, which
to our knowledge is not observed experimentally. While a full understanding is
somewhat involved, we argue below that both thermal fluctuations in two
dimensions and three-dimensional coupling stabilize a continuous transition
within most scenarios.

####  V.3.1 Two dimensions

If we neglect inter-layer coupling, the problem becomes two dimensional, and
we must be wary of applying a mean field Landau theory analysis to critical
properties, since it is well known that two dimensional systems have strong
thermal fluctuations.

For the case of the 3Q$+$ or 3Q$-$ rCDW, as discussed in Sec. V.2.1, the
ordering can be described by a four state Potts model. At mean field level,
the $q$-state Potts models have first order transitions for $q\geq 3$,
consistent with the Landau analysis. It is well-known, however, that in two
dimensions the $q=3,4$ Potts models in fact have continuous transitions
described by conformal field theory. In particular, the $q=4$ Potts model is
equivalent to the Ashkin-Teller model, and has known critical exponents with
logarithmic corrections due to marginal operators (see e.g. Ref.[63]).

In the case of the 3Q iCDW, the order parameter has a degeneracy of 8 and
transforms under full cubic $O_{h}$ symmetry, and we have seen that it can be
expressed as an O(3) vector with cubic anisotropy. In two dimensions, this
problem has been analyzed by Schick[64], who concludes that the transition to
the 3Q phase (called “corner cubic anisotropy” in this reference) may be
either continuous (and in the Ising universality class, surprisingly) or first
order.

####  V.3.2 Three dimensions

The two dimensional critical behavior is valid at best in a regime close but
not too close to the critical point, where a crossover to three-dimensional
behavior must occur. The nature of the true three-dimensional critical regime
is, however, dependent on the type of inter-layer couplings, and thereby the
three-dimensional order parameter. As discussed in Sec. V.2.1, this is rather
complicated for the rCDW problem, and we will not offer a complete analysis.
The simplest case is the $k_{z}=0$ one, which occurs when the inter-layer
coupling favors ferromagnetically aligned rCDWs. Then the symmetry remains
unchanged from the four state Potts model. However, in three dimensions, this
model has a first order transition. Thus in this case a first order transition
is predicted. If instead one of the $k_{z}>0$ orderings occur, the situation
is less clear, but it is manifestly not a Potts model transition. It is
natural to think that the Landau analysis is more correct in three dimensions.
Based on this reasoning, for the cases $k_{z}=1/2$ and $k_{z}=1/4$, a cubic
term in the order parameter is no longer allowed by momentum conservation.
Thus a continuous transition may be expected. The case $k_{z}=1/3$ seems to
allow a cubic term but we will not pursue it further here.

With three-dimensional coupling, the iCDW order develops either at $k_{z}=0$
or $k_{z}=1/2$. In both cases, the three-dimensional order parameter remains
an $O(3)$ vector with cubic anisotropy, and the ground state degeneracy is
unchanged from $8$ (this can be seen because even in the $k_{z}=1/2$ can a
translation by one layer is equivalent to time-reversal). Hence the transition
should be in the O(3) cubic universality class, which is known to be
continuous, and is discussed e.g. in Ref.[65].

###  V.4 Magnetic moment induced by iCDW order

The iCDW order can induce both staggered magnetic moment due to the loop
current and uniform magnetic moment. Here, we estimate the magnitude from each
contribution.

First, we note that the bond current obtained in Eq. (23) is linear in $\phi$,
consequently, only the Fourier component at ${\bm{M}}$ contributes to the bond
current, and there is no uniform magnetic moment induced by the bond current.
To estimate the staggered magnetic moment, we consider the 3Q iCDW order as an
example. As shown in Fig. 8 (b), the bond current forms loop current around
both the honeycomb and triangular plaquettes. Treating each plaquette as a
current loop, its magnetic moment and the magnetic field it induces can be
obtained following the standard magnetostatics [66]. Due to the $2\times 2$
sublattice structure, there are two types of honeycomb plaquettes, with
${\bf{m}}_{\rm h,1}=I_{\rm e}\hat{\mathcal{A}}_{\rm h}$ for one-quarter of the
honeycomb plaquettes, and ${\bf{m}}_{\rm h,2}=-I_{\rm e}\hat{\mathcal{A}}_{\rm
h}/3$ for the rest. Similarly, one-quarter of the triangular plaquettes have
${\bf{m}}_{\rm t,1}=-I_{\rm e}\hat{\mathcal{A}}_{\rm t}$, and the rest have
${\bf{m}}_{\rm t,2}=I_{\rm e}\hat{\mathcal{A}}_{\rm t}/3$. Here, $I_{\rm
e}=-\frac{e}{\hbar}\frac{2t_{1}}{G_{\rm iCDW}}\phi_{0}$ denotes the magnitude
of the electric loop current. $\hat{\mathcal{A}}_{\rm h,t}$ is along
$\hat{z}$, such that for $I_{\rm e}>0$, the current flows counter clockwisely.
And its magnitude is determined by the area of the honeycomb (h) and
triangular (t) plaquette. In unit of Bohr magneton
$\mu_{B}=\frac{e\hbar}{2m_{e}}$, we have ${\bf{m}}\sim I_{\rm
e}\hat{\mathcal{A}}=-\frac{2t_{1}}{G_{\rm
iCDW}}\frac{\phi_{0}}{Ry}\frac{\hat{\mathcal{A}}}{a_{b}^{2}}\mu_{B}$, where
$Ry=\frac{m_{e}e^{4}}{32\pi^{2}\epsilon_{0}^{2}\hbar^{2}}\approx 13.6eV$ is
the Rydberg energy, $r_{b}=\frac{4\pi\epsilon_{0}\hbar}{m_{e}e^{2}}\approx
5.29\times 10^{-11}m$ is the Bohr radius. In $I_{\rm e}$, the factor
$\frac{2t_{1}}{G_{\rm iCDW}}$ is dimensionless and scales with the
polarization bubble $\sim\ln^{2}(t\Lambda^{2}/T_{c})\sim\mathcal{O}(1)$, the
order parameter (in unit of energy) at low temperature can be approximated as
$\phi_{0}\sim k_{B}T_{c}$, where $T_{c}$ is the transition temperature for
iCDW order. This gives ${\bf{m}}_{\rm
h,1}\approx-0.047\operatorname{sgn}{(\phi_{0})}\mu_{B}\hat{z},{\bf{m}}_{\rm
t,1}\approx 0.008\operatorname{sgn}{(\phi_{0})}\mu_{B}\hat{z}$, where we take
$2t_{1}/G_{\rm iCDW}=1$. We also note that similar analysis has been done in
the literature for cuprates [67] and iron-pnictides [68].

Next, we discuss the uniform magnetization induced by the iCDW order. On the
symmetry ground, uniform magnetization requires iCDW order at all three
momenta ${\bm{M}}$ nonzero. Otherwise, the system is invariant under the anti-
unitary symmetry composed of time reversal and translation, which forbids any
uniform magnetization. In general, the magnetic field couples to the electrons
through both the minimal coupling and the Zeeman coupling. Here, we consider
the orbital magnetization contribution for the 3Q$\pm$ iCDW order. Following
[69], the orbital magnetization in terms of the Bloch wave function reads,

$\displaystyle\mathcal{M}_{{\rm orb}}^{z}=\frac{e}{2\hbar}\sum_{n,{\bm{q}}}$
$\displaystyle\quad\operatorname{Im}\left[\langle\partial_{{\bm{q}}}{\rm
u}_{n,{\bm{q}}}|\times(\epsilon_{n,{\bm{q}}}-\hat{H}_{\rm
MF}({\bm{q}}))|\partial_{{\bm{q}}}{\rm
u}_{n,{\bm{q}}}\rangle\right]f_{n,{\bm{q}}},$ (28)

where $\hat{H}_{\rm MF}({\bm{q}})$ is the mean field Hamiltonian one can infer
from Eq. (8), $f_{n,{\bm{q}}}=f(\epsilon_{n,{\bm{q}}}-\mu)$ is the equilibrium
Fermi distribution function.

For simplicity, we consider the perfect nesting case when the Fermi energy in
the disordered phase is at the van-hove point, this gives $1/3$ filling in the
patch model. When the 3Q iCDW order is present, the triple degenerate bands at
the ${\bm{M}}$ points are fully gapped, and only the lowest band is filled.
Noting that the unit of the summand in Eq. (28) is $\frac{\hbar^{2}}{2m_{e}}$
(coming from the kinetic energy), for a filled band, the summand must be
expressed as a function of the only dimensionless parameter in the expression,
which reads $\tilde{\epsilon}=\hbar^{2}q^{2}/(2m_{e})$, i.e. the summand
(after averaging over the angular direction in ${\bm{q}}$) must be expressed
as $\sim\frac{\hbar^{2}}{2m_{e}}\mathcal{F}(\tilde{\epsilon})$. This means
physically that for electrons within an energy of order the gap $\phi$, the
typical orbital moment is an order one fraction of a Bohr magneton.

The orbital magnetization can be expressed as

$\displaystyle\mathcal{M}_{{\rm orb}}^{z}$
$\displaystyle=\frac{e\hbar}{2m_{e}}\mathcal{N}\frac{\phi_{0}}{\hbar^{2}/(2m_{e}a_{0}^{2})}\frac{\sqrt{3}}{2\pi}\left[\int_{0}^{\Lambda}\mathop{}\\!\mathrm{d}\tilde{\epsilon}\mathcal{F}(\tilde{\epsilon})\right]$
$\displaystyle=\mathcal{N}\mu_{B}\frac{\phi_{0}}{Ry}\left(\frac{a_{0}}{a_{b}}\right)^{2}\frac{\sqrt{3}}{2\pi}\mathcal{I}_{F},$
(29)

where $\mathcal{N}$ is the number of unit cells,
$\mathcal{I}_{F}=\int_{0}^{\Lambda}\mathop{}\\!\mathrm{d}\tilde{\epsilon}\mathcal{F}(\tilde{\epsilon})$
is the value of the integral, which is bounded and can be computed numerically
as $\mathcal{I}_{F}\approx-0.87$. Again approximating $\phi_{0}\sim
k_{B}T_{c}$, the orbital magnetization per unit cell is $\mathcal{M}_{{\rm
orb}}^{z}/N\sim-0.017\operatorname{sgn}(\phi_{0})\mu_{B}$.

This orbital magnetization is small but the smallness is primarily due to the
small number of electrons within the region near the energy gap. As remarked
above, for a typical electron within this region, the orbital moment is a
substantial fraction of $\mu_{B}$. This suggests that the iCDW state may be
favored by the application of a magnetic field.

We consider therefore how an rCDW state may be converted into an iCDW state by
such a field, assuming the energy of the iCDW state is not too much higher
than that of the rCDW in zero field. For simplicity we consider low
temperature, where the above estimate is valid (near $T_{c}$, the orbital
magnetization will be proportional to $\phi^{3}$, and greatly suppressed). In
general, a magnetic field breaks time-reversal, so will always induce some
iCDW component if the rCDW is present. We thus take
$\Delta_{\alpha}=N_{\alpha}+i\phi_{\alpha}=|\Delta|e^{i\theta}$ (assuming a
3Q$+$ state), and express the energy in terms of $|\Delta|$ and $\theta$. A
simplified form for the energy density is

$\mathcal{E}=-m|\Delta|^{2}\left(1+\lambda\cos\theta\right)-\mu
hm|\Delta|\sin\theta,$ (30)

where $m$ is the mass scale for the saddle points, $h$ is the magnetic field
(in the $z$ direction), $\mu$ is a typical orbital moment for the iCDW
electrons, and is a fraction of $\mu_{B}$. The first term is the condensation
energy of the CDW, and the second term is the dipole energy associated with
the orbital magnetization. The parameter $\lambda>0$ describes the energetic
preference for rCDW over iCDW order in zero field (as well as the preference
for 3Q+ over 3Q- order). When the competition between the two is close,
$\lambda\ll 1$. Minimizing over $\theta$, we obtain

$\tan\theta=\frac{\mu h}{\lambda|\Delta|}.$ (31)

We see that the rCDW smoothly evolves into an iCDW with applied field, and
this occurs on a scale which can be a small fraction of the gap, if
$\lambda\ll 1$. We note that the energy density in Eq. (30) is simplified, and
does not capture topological physics relevant to the perfect nesting situation
when the system is doped exactly to the saddle point filling. In this case,
the system is a trivial insulator for $\theta=0$, and a Chern insulator for
$\theta=\pi/2$, and hence there must be a topological transition, associated
gap closing, and non-analyticity of the energy at some intermediate angle.
However, we expect this distinction to be washed out away from perfect nesting
and for generic filling.

##  VI Summary

In this paper, we have discussed various mechanisms to induce charge density
wave order in kagomé metals, motivated by experiments on the AV3Sb5 materials.
We began with a g-ology description based on a continuum saddle point model,
which is quite general to two dimensional materials with hexagonal symmetry,
and was originally introduced in the context of doped graphene in Ref.[46]. We
extended the renormalization group analysis of this problem to the general
case, i.e. with both repulsive and attractive bare interactions, and showed
that several distinct density wave and superconducting instabilities are
possible, depending upon parameters. We also showed that the attractive bare
interactions may come from coupling of electrons to optical phonons, which
induces attraction in certain channels ($g_{i}$) and lowers the phonon energy.
The latter effect may be relevant to the phonon softening observed in ab-
initio calculations [70]. Our focus then turned to charge density wave order
at the vectors $\bm{M}_{\alpha}$, which correspond both to the locations of
the saddle points in the Brillouin zone and the spanning vectors between them.
They are half reciprocal lattice vectors and reside at the centers of the zone
faces. We applied a mean field theory to the real charge density wave (rCDW)
states of the continuum interacting model, and found isotropic 3Q$\pm$ as well
as nematic 1Q rCDW states occur and may be tuned by chemical potential. We
also carried out a Landau theory analysis of the rCDW orders (and a detailed
derivation of the Landau theory coefficients from the continuum model) which
shows that when this is the leading instability, the above states are the only
ones that occur with a high degree of generality. Next, we considered the
alternate cases of a 3Q imaginary charge density wave (iCDW) – actually a
state of circulating currents, i.e. orbital moments – and a spin density wave
(rSDW), which are competing instabilities with the same wavevectors. We showed
that both iCDW and rSDW can induce rCDW order, though in both cases the charge
density order is then not the primary order parameter. Finally, we detailed a
number of extensions and implications of the analysis to derive explicit
charge and current patterns, three dimensional ordered states, critical
properties, and orbital magnetization of the 3Q iCDW state.

What are the specific implications of these results for the AV3Sb5 materials?
Three distinct compounds, with A=K,Cs,Rb, have been studied to our knowledge.
All are known to show rCDW order for which the projection of the ordering
wavevector to the 2d kagomé plane is of the 3Q $\bm{M}_{\alpha}$ type. A key
question is whether the rCDW is the primary order parameter, or whether it is
secondary and induced, as discussed above, by either an iCDW or rSDW order.
Notably, both the iCDW and rSDW orders break time-reversal symmetry, and hence
should be detectable via their induced local magnetic moments through e.g.
muon spin resonance or neutron scattering. To our knowledge, there is no
direct evidence of time-reversal breaking from any measurement in zero
magnetic field (and there is some counter-evidence from muon spin resonance on
KV3Sb5[26]). This does not definitively exclude states with very small
moments, as indeed may be expected in the iCDW case (c.f. Sec. V.4). The
discussion of three-dimensional order in Sec. V.2.1 provides some more clear
diagnostics. There, we showed that the iCDW is expected to display identical
induced rCDW order in all layers, i.e. the z-component of the rCDW ordering
wavevector should be zero. By contrast, a pure rCDW is consistent with
$k_{z}=0,\tfrac{1}{2},\tfrac{1}{3},\tfrac{1}{4}$. Moreover, rCDW states with
$k_{z}\neq 0$ should show “nematic” $C_{3}$ rotational symmetry breaking at
the surface, due to the three-dimensional coupling, which again is not
expected in the iCDW states.

Several experimental papers are directly related to the above aspects of the
CDW order. X-ray scattering found Bragg peaks associated with three-
dimensional CDW order with $k_{z}=1/2$ in RbV3Sb5, CsV3Sb5 in Ref.[22],
$k_{z}=1/4$ in CsV3Sb5 in Ref.[27]. Ref.[24] measured an in-plane shift of the
order in the top two layers of CsV3Sb5 across a step edge using STM,
indicating $k_{z}>0$. These observations are consistent with primary rCDW
order. Ref.[36] measured a “nematic” (i.e. $C_{2}$ symmetric) dependence of
the c-axis resistivity in CsV3Sb5 on an in-plane magnetic field in the normal
(but CDW) state up to about 60K, 2/3 of the critical temperature. Such
behavior would be expected below the CDW $T_{c}$ for any rCDW except the
“screw” state with $k_{z}=1/3$. The STM study in Ref. [20] notes differences
in the intensities of the three charge order wavevector peaks in KV3Sb5. This
reference concludes that the three intensities are all unequal, and denote
this as a “chiral” charge order. To our eyes, their Fourier transform STM data
(Fig.3abc of Ref.[20]) is better described by two strong and approximately
equal intensities and one weaker one, which is consistent with our theory of
inter-layer coupling. A very recent paper reports the latter type of
anisotropy in independent STM measurements of KV3Sb5 [71].

While the above evidence seems to favor primary rCDW order, iCDW or SDW states
might be realized in some situations, e.g. as competing orders that appear at
lower temperatures or induced through applied fields. Indeed the early
observation of large anomalous Hall effect[72] in KV3Sb5 motivated several
later works to promote the possibility of an iCDW (recently similar behavior
was observed in CsV3Sb5[38]) . As discussed in Sec. V.4, the 3Q iCDW does
indeed macroscopically break time-reversal and is accompanied by a uniform
orbital moment and a large Berry curvature in the vicinity of the saddle
points, which should induce an anomalous Hall effect. We discussed how, if
iCDW order is only slightly higher in energy than rCDW, the iCDW can be
induced by modest applied fields, which may provide a possible explanation of
Ref.[72]. Further study of the anomalous Hall effect and its correlation with
other measurements is however needed to test this possibility.

While this manuscript was in preparation, several theoretical works appeared
discussing density wave order in these materials. Ref. [73] evaluated the mean
field ground state energies for multiple density wave orders. Ref.[74]
considered a single orbital extended Hubbard model on the kagomé lattice, and
within mean field theory found dominant iCDW order. Ref.[75] discusses a
coupled mean field theory of spin and charge orders and Chern bands in a
continuum saddle point model similar to the one used here and postulates a
complex CDW state in the AV3Sb5 materials with three unequal magnitudes of CDW
at the three $\bm{M}_{\alpha}$ wavevectors.

In this paper, we have focused on density wave order, which is clearly the
dominant instability in the AV3Sb5 materials. However, the superconductivity
occurring at lower temperature is of considerable interest as well. Since it
occurs within the density wave ordered state, the understanding of the latter
should be important for developing a theory of superconductivity in this
family.

###### Acknowledgements.

We acknowledge useful discussions with Stephen Wilson, Brenden Ortiz, Sam
Teicher, and Ilija Zeljkovic. L.B. is supported by BES Award DE-SC0020305.
M.Y. is supported in part by the Gordon and Betty Moore Foundation through
Grant GBMF8690 to UCSB and by the National Science Foundation under Grant No.
NSF PHY-1748958. T.P. was supported by the National Science Foundation through
Enabling Quantum Leap: Convergent Accelerated Discovery Foundries for Quantum
Materials Science, Engineering and Information (Q-AMASE-i) award number
DMR-1906325.

## References

  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998).
  * Ngai  _et al._ [2014] J. Ngai, F. Walker, and C. Ahn, Correlated oxide physics and electronics, Annual Review of Materials Research 44, 1 (2014), https://doi.org/10.1146/annurev-matsci-070813-113248 .
  * Johnston [2010] D. C. Johnston, The puzzle of high temperature superconductivity in layered iron pnictides and chalcogenides, Advances in Physics 59, 803 (2010), https://doi.org/10.1080/00018732.2010.513480 .
  * Scalapino [2012] D. J. Scalapino, A common thread: The pairing interaction for unconventional superconductors, Rev. Mod. Phys. 84, 1383 (2012).
  * Sato and Ando [2017] M. Sato and Y. Ando, Topological superconductors: a review, Reports on Progress in Physics 80, 076501 (2017).
  * Tsen  _et al._ [2015] A. W. Tsen, R. Hovden, D. Wang, Y. D. Kim, J. Okamoto, K. A. Spoth, Y. Liu, W. Lu, Y. Sun, J. C. Hone, L. F. Kourkoutis, P. Kim, and A. N. Pasupathy, Structure and control of charge density waves in two-dimensional 1t-tas2, Proceedings of the National Academy of Sciences 112, 15054 (2015), https://www.pnas.org/content/112/49/15054.full.pdf .
  * Ye  _et al._ [2018] L. Ye, M. Kang, J. Liu, F. Von Cube, C. R. Wicker, T. Suzuki, C. Jozwiak, A. Bostwick, E. Rotenberg, D. C. Bell, _et al._ , Massive dirac fermions in a ferromagnetic kagome metal, Nature 555, 638 (2018).
  * Fisher  _et al._ [2011] I. R. Fisher, L. Degiorgi, and Z. X. Shen, In-plane electronic anisotropy of underdoped “122” fe-arsenide superconductors revealed by measurements of detwinned single crystals, Reports on Progress in Physics 74, 124506 (2011).
  * Sato  _et al._ [2017] Y. Sato, S. Kasahara, H. Murayama, Y. Kasahara, E. G. Moon, T. Nishizaki, T. Loew, J. Porras, B. Keimer, T. Shibauchi, and Y. Matsuda, Thermodynamic evidence for a nematic phase transition at the onset of the pseudogap in yba2cu3oy, Nature Physics 13, 1074 (2017).
  * Varma  _et al._ [1989] C. M. Varma, P. B. Littlewood, S. Schmitt-Rink, E. Abrahams, and A. E. Ruckenstein, Phenomenology of the normal state of cu-o high-temperature superconductors, Phys. Rev. Lett. 63, 1996 (1989).
  * Stewart [2001] G. R. Stewart, Non-fermi-liquid behavior in $d$\- and $f$-electron metals, Rev. Mod. Phys. 73, 797 (2001).
  * Gunnarsson  _et al._ [2003] O. Gunnarsson, M. Calandra, and J. E. Han, Colloquium: Saturation of electrical resistivity, Rev. Mod. Phys. 75, 1085 (2003).
  * Varma [2020] C. M. Varma, Colloquium: Linear in temperature resistivity and associated mysteries including high temperature superconductivity, Rev. Mod. Phys. 92, 031001 (2020).
  * Chubukov [2012] A. Chubukov, Pairing mechanism in fe-based superconductors, Annual Review of Condensed Matter Physics 3, 57 (2012), https://doi.org/10.1146/annurev-conmatphys-020911-125055 .
  * Si  _et al._ [2016] Q. Si, R. Yu, and E. Abrahams, High-temperature superconductivity in iron pnictides and chalcogenides, Nature Reviews Materials 1, 16017 (2016).
  * Takada  _et al._ [2003] K. Takada, H. Sakurai, E. Takayama-Muromachi, F. Izumi, R. A. Dilanian, and T. Sasaki, Superconductivity in two-dimensional coo2 layers, Nature 422, 53 (2003).
  * Foo  _et al._ [2004] M. L. Foo, Y. Wang, S. Watauchi, H. W. Zandbergen, T. He, R. J. Cava, and N. P. Ong, Charge ordering, commensurability, and metallicity in the phase diagram of the layered ${\mathrm{na}}_{x}\mathrm{Co}{\mathrm{o}}_{2}$, Phys. Rev. Lett. 92, 247001 (2004).
  * Jorgensen  _et al._ [2003] J. D. Jorgensen, M. Avdeev, D. G. Hinks, J. C. Burley, and S. Short, Crystal structure of the sodium cobaltate deuterate superconductor ${\mathrm{na}}_{x}{\mathrm{coo}}_{2}\cdot 4x{\mathrm{d}}_{2}\mathrm{O}$ $(x\approx\frac{1}{3})$, Phys. Rev. B 68, 214517 (2003).
  * Ortiz  _et al._ [2019] B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, _et al._ , New kagome prototype materials: discovery of KV3Sb5, RbV3Sb5, and CsV3Sb5, Physical Review Materials 3, 094407 (2019).
  * Jiang  _et al._ [2020] Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, J. He, X. Liu, S. S. Zhang, G. Chang, I. Belopolski, Q. Zhang, M. S. Hossain, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Guguchia, G. Xu, Z. Wang, T. Neupert, S. D. Wilson, and M. Z. Hasan, Discovery of topological charge order in kagome superconductor kv3sb5 (2020), arXiv:2012.15709 [cond-mat.supr-con] .
  * Zhao  _et al._ [2021a] H. Zhao, H. Li, B. R. Ortiz, S. M. L. Teicher, T. Park, M. Ye, Z. Wang, L. Balents, S. D. Wilson, and I. Zeljkovic, Cascade of correlated electron states in a kagome superconductor CsV3Sb5, arXiv e-prints , arXiv:2103.03118 (2021a), arXiv:2103.03118 [cond-mat.supr-con] .
  * Li  _et al._ [2021] H. X. Li, T. T. Zhang, Y. Y. Pai, C. Marvinney, A. Said, T. Yilmaz, Q. Yin, C. Gong, Z. Tu, E. Vescovo, R. G. Moore, S. Murakami, H. C. Lei, H. N. Lee, B. Lawrie, and H. Miao, Observation of Unconventional Charge Density Wave without Acoustic Phonon Anomaly in Kagome Superconductors AV3Sb5 (A=Rb,Cs), arXiv e-prints , arXiv:2103.09769 (2021), arXiv:2103.09769 [cond-mat.supr-con] .
  * Uykur  _et al._ [2021] E. Uykur, B. R. Ortiz, S. D. Wilson, M. Dressel, and A. A. Tsirlin, Optical detection of charge-density-wave instability in the non-magnetic kagome metal KV3Sb5, arXiv e-prints , arXiv:2103.07912 (2021), arXiv:2103.07912 [cond-mat.str-el] .
  * Liang  _et al._ [2021] Z. Liang, X. Hou, W. Ma, F. Zhang, P. Wu, Z. Zhang, F. Yu, J. J. Ying, K. Jiang, L. Shan, Z. Wang, and X. H. Chen, Three-dimensional charge density wave and robust zero-bias conductance peak inside the superconducting vortex core of a kagome superconductor CsV3Sb5, arXiv e-prints , arXiv:2103.04760 (2021), arXiv:2103.04760 [cond-mat.supr-con] .
  * Zhou  _et al._ [2021] X. Zhou, Y. Li, X. Fan, J. Hao, Y. Dai, Z. Wang, Y. Yao, and H.-H. Wen, Origin of the Charge Density Wave in the Kagome Metal CsV3Sb5 as Revealed by Optical Spectroscopy, arXiv e-prints , arXiv:2104.01015 (2021), arXiv:2104.01015 [cond-mat.supr-con] .
  * Kenney  _et al._ [2021] E. M. Kenney, B. R. Ortiz, C. Wang, S. D. Wilson, and M. Graf, Absence of local moments in the kagome metal kv3sb5 as determined by muon spin spectroscopy, Journal of Physics: Condensed Matter (2021).
  * Ortiz  _et al._ [2021a] B. R. Ortiz, S. M. L. Teicher, L. Kautzsch, P. M. Sarte, J. P. C. Ruff, R. Seshadri, and S. D. Wilson, Fermi surface mapping and the nature of charge density wave order in the kagome superconductor csv3sb5 (2021a), arXiv:2104.07230 [cond-mat.str-el] .
  * Ortiz  _et al._ [2021b] B. R. Ortiz, P. M. Sarte, E. M. Kenney, M. J. Graf, S. M. L. Teicher, R. Seshadri, and S. D. Wilson, Superconductivity in the ${\mathbb{z}}_{2}$ kagome metal ${\mathrm{kv}}_{3}{\mathrm{sb}}_{5}$, Phys. Rev. Materials 5, 034801 (2021b).
  * Zhao  _et al._ [2021b] C. C. Zhao, L. S. Wang, W. Xia, Q. W. Yin, J. M. Ni, Y. Y. Huang, C. P. Tu, Z. C. Tao, Z. J. Tu, C. S. Gong, H. C. Lei, Y. F. Guo, X. F. Yang, and S. Y. Li, Nodal superconductivity and superconducting domes in the topological Kagome metal CsV3Sb5, arXiv e-prints , arXiv:2102.08356 (2021b), arXiv:2102.08356 [cond-mat.supr-con] .
  * Chen  _et al._ [2021a] K. Y. Chen, N. N. Wang, Q. W. Yin, Z. J. Tu, C. S. Gong, J. P. Sun, H. C. Lei, Y. Uwatoko, and J. G. Cheng, Double superconducting dome and triple enhancement of Tc in the kagome superconductor CsV3Sb5 under high pressure, arXiv e-prints , arXiv:2102.09328 (2021a), arXiv:2102.09328 [cond-mat.supr-con] .
  * Chen  _et al._ [2021b] H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, Q. Yin, C. Gong, Z. Tu, H. Lei, S. Ma, H. Zhang, S. Ni, H. Tan, C. Shen, X. Dong, B. Yan, Z. Wang, and H.-J. Gao, Roton pair density wave and unconventional strong-coupling superconductivity in a topological kagome metal, arXiv e-prints , arXiv:2103.09188 (2021b), arXiv:2103.09188 [cond-mat.supr-con] .
  * Duan  _et al._ [2021] W. Duan, Z. Nie, S. Luo, F. Yu, B. R. Ortiz, L. Yin, H. Su, F. Du, A. Wang, Y. Chen, X. Lu, J. Ying, S. D. Wilson, X. Chen, Y. Song, and H. Yuan, Nodeless superconductivity in the kagome metal CsV3Sb5, arXiv e-prints , arXiv:2103.11796 (2021), arXiv:2103.11796 [cond-mat.supr-con] .
  * Zhang  _et al._ [2021] Z. Zhang, Z. Chen, Y. Zhou, Y. Yuan, S. Wang, L. Zhang, X. Zhu, Y. Zhou, X. Chen, J. Zhou, and Z. Yang, Pressure-induced Reemergence of Superconductivity in Topological Kagome Metal CsV3Sb5, arXiv e-prints , arXiv:2103.12507 (2021), arXiv:2103.12507 [cond-mat.supr-con] .
  * Mu  _et al._ [2021] C. Mu, Q. Yin, Z. Tu, C. Gong, H. Lei, Z. Li, and J. Luo, $s$-wave superconductivity in kagome metal CsV3Sb5 revealed by 121/123Sb NQR and 51V NMR measurements, arXiv e-prints , arXiv:2104.06698 (2021), arXiv:2104.06698 [cond-mat.supr-con] .
  * Ni  _et al._ [2021] S. Ni, S. Ma, Y. Zhang, J. Yuan, H. Yang, Z. Lu, N. Wang, J. Sun, Z. Zhao, D. Li, S. Liu, H. Zhang, H. Chen, K. Jin, J. Cheng, L. Yu, F. Zhou, X. Dong, J. Hu, H.-J. Gao, and Z. Zhao, Anisotropic superconducting properties of Kagome metal CsV3Sb5, arXiv e-prints , arXiv:2104.00374 (2021), arXiv:2104.00374 [cond-mat.supr-con] .
  * Xiang  _et al._ [2021] Y. Xiang, Q. Li, Y. Li, W. Xie, H. Yang, Z. Wang, Y. Yao, and H.-H. Wen, Nematic electronic state and twofold symmetry of superconductivity in the topological kagome metal csv3sb5 (2021), arXiv:2104.06909 [cond-mat.supr-con] .
  * Wang  _et al._ [2020] Q. Wang, Q. Yin, and H. Lei, Giant topological Hall effect of ferromagnetic kagome metal Fe3Sn2, Chinese Physics B 29, 017101 (2020).
  * Yu  _et al._ [2021] F. H. Yu, T. Wu, Z. Y. Wang, B. Lei, W. Z. Zhuo, J. J. Ying, and X. H. Chen, Concurrence of anomalous Hall effect and charge density wave in a superconducting topological kagome metal, arXiv e-prints , arXiv:2102.10987 (2021), arXiv:2102.10987 [cond-mat.str-el] .
  * Ortiz  _et al._ [2020] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, $\mathrm{Cs}{\mathrm{v}}_{3}{\mathrm{sb}}_{5}$: A ${\mathbb{z}}_{2}$ topological kagome metal with a superconducting ground state, Phys. Rev. Lett. 125, 247002 (2020).
  * Zhao  _et al._ [2021c] J. Zhao, W. Wu, Y. Wang, and S. A. Yang, Electronic correlations in the normal state of kagome superconductor KV3Sb5, arXiv e-prints , arXiv:2103.15078 (2021c), arXiv:2103.15078 [cond-mat.str-el] .
  * Liu  _et al._ [2021] Z. Liu, N. Zhao, Q. Yin, C. Gong, Z. Tu, M. Li, W. Song, Z. Liu, D. Shen, Y. Huang, K. Liu, H. Lei, and S. Wang, Temperature-induced band renormalization and Lifshitz transition in a kagome superconductor RbV3Sb5, arXiv e-prints , arXiv:2104.01125 (2021), arXiv:2104.01125 [cond-mat.supr-con] .
  * Wang  _et al._ [2021] Z. Wang, S. Ma, Y. Zhang, H. Yang, Z. Zhao, Y. Ou, Y. Zhu, S. Ni, Z. Lu, H. Chen, K. Jiang, L. Yu, Y. Zhang, X. Dong, J. Hu, H.-J. Gao, and Z. Zhao, Distinctive momentum dependent charge-density-wave gap observed in CsV3Sb5 superconductor with topological Kagome lattice, arXiv e-prints , arXiv:2104.05556 (2021), arXiv:2104.05556 [cond-mat.supr-con] .
  * Hirsch and Scalapino [1986] J. E. Hirsch and D. J. Scalapino, Enhanced superconductivity in quasi two-dimensional systems, Phys. Rev. Lett. 56, 2732 (1986).
  * Markiewicz [1997] R. Markiewicz, A survey of the Van Hove scenario for high-tc superconductivity with special emphasis on pseudogaps and striped phases, Journal of Physics and Chemistry of Solids 58, 1179 (1997), arXiv:cond-mat/9611238 [cond-mat.supr-con] .
  * Hur and Maurice Rice [2009] K. L. Hur and T. Maurice Rice, Superconductivity close to the mott state: From condensed-matter systems to superfluidity in optical lattices, Annals of Physics 324, 1452 (2009), july 2009 Special Issue.
  * Nandkishore  _et al._ [2012a] R. Nandkishore, L. Levitov, and A. Chubukov, Chiral superconductivity from repulsive interactions in doped graphene, Nature Physics 8, 158 (2012a).
  * Nandkishore  _et al._ [2012b] R. Nandkishore, G.-W. Chern, and A. V. Chubukov, Itinerant Half-Metal Spin-Density-Wave State on the Hexagonal Lattice, Physical Review Letters 108, 227204 (2012b).
  * Nandkishore and Chubukov [2012] R. Nandkishore and A. V. Chubukov, Interplay of superconductivity and spin-density-wave order in doped graphene, Physical Review B 86, 115426 (2012).
  * Lin and Nandkishore [2019] Y.-P. Lin and R. M. Nandkishore, Chiral twist on the high-${T}_{c}$ phase diagram in moiré heterostructures, Phys. Rev. B 100, 085136 (2019).
  * Chichinadze  _et al._ [2020a] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Nematic superconductivity in twisted bilayer graphene, Phys. Rev. B 101, 224513 (2020a).
  * Chichinadze  _et al._ [2020b] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Valley magnetism, nematicity, and density wave orders in twisted bilayer graphene, Phys. Rev. B 102, 125120 (2020b).
  * Classen  _et al._ [2020] L. Classen, A. V. Chubukov, C. Honerkamp, and M. M. Scherer, Competing orders at higher-order van hove points, Phys. Rev. B 102, 125141 (2020).
  * Lin and Nandkishore [2020] Y.-P. Lin and R. M. Nandkishore, Parquet renormalization group analysis of weak-coupling instabilities with multiple high-order van hove points inside the brillouin zone, Phys. Rev. B 102, 245122 (2020).
  * Zheleznyak  _et al._ [1997] A. T. Zheleznyak, V. M. Yakovenko, and I. E. Dzyaloshinskii, Parquet solution for a flat fermi surface, Phys. Rev. B 55, 3200 (1997).
  * Metzner  _et al._ [1998] W. Metzner, C. Castellani, and C. D. Castro, Fermi systems with strong forward scattering, Advances in Physics 47, 317 (1998), https://doi.org/10.1080/000187398243528 .
  * Chubukov  _et al._ [2008] A. V. Chubukov, D. V. Efremov, and I. Eremin, Magnetism, superconductivity, and pairing symmetry in iron-based superconductors, Phys. Rev. B 78, 134512 (2008).
  * Chubukov [2009] A. Chubukov, Renormalization group analysis of competing orders and the pairing symmetry in fe-based superconductors, Physica C: Superconductivity 469, 640 (2009), superconductivity in Iron-Pnictides.
  * [58] S. M. Teicher, (private communication).
  * Koster [1963] G. F. Koster,  _Properties of the thirty-two point groups_ (M.I.T. Press, Cambridge, Mass., 1963).
  * Furukawa  _et al._ [1998] N. Furukawa, T. M. Rice, and M. Salmhofer, Truncation of a Two-Dimensional Fermi Surface due to Quasiparticle Gap Formation at the Saddle Points, Physical Review Letters 81, 3195 (1998).
  * Lin  _et al._ [1997] H.-H. Lin, L. Balents, and M. P. A. Fisher, $n$-chain hubbard model in weak coupling, Phys. Rev. B 56, 6569 (1997).
  * Balents and Fisher [1996] L. Balents and M. P. Fisher, Weak-coupling phase diagram of the two-chain hubbard model, Physical Review B 53, 12133 (1996).
  * Cardy  _et al._ [1980] J. L. Cardy, M. Nauenberg, and D. Scalapino, Scaling theory of the potts-model multicritical point, Physical Review B 22, 2560 (1980).
  * Schick [1983] M. Schick, Application of the heisenberg model with cubic anisotropy to phase transitions on surfaces, Surface Science 125, 94 (1983).
  * Manuel Carmona  _et al._ [2000] J. Manuel Carmona, A. Pelissetto, and E. Vicari, $n$-component ginzburg-landau hamiltonian with cubic anisotropy: A six-loop study, Phys. Rev. B 61, 15136 (2000).
  * Jackson [1999] J. D. Jackson,  _Classical electrodynamics; 2nd ed._ (Wiley, New York, NY, 1999).
  * Lederer and Kivelson [2012] S. Lederer and S. A. Kivelson, Observable nmr signal from circulating current order in ybco, Phys. Rev. B 85, 155130 (2012).
  * Klug  _et al._ [2018] M. Klug, J. Kang, R. M. Fernandes, and J. Schmalian, Orbital loop currents in iron-based superconductors, Phys. Rev. B 97, 155130 (2018).
  * Shi  _et al._ [2007] J. Shi, G. Vignale, D. Xiao, and Q. Niu, Quantum theory of orbital magnetization and its generalization to interacting systems, Physical Review Letters 99, 10.1103/physrevlett.99.197202 (2007).
  * Tan  _et al._ [2021] H. Tan, Y. Liu, Z. Wang, and B. Yan, Charge density waves and electronic properties of superconducting kagome metals, arXiv e-prints , arXiv:2103.06325 (2021), arXiv:2103.06325 [cond-mat.supr-con] .
  * Li  _et al._ [2021] H. Li, H. Zhao, B. R. Ortiz, T. Park, M. Ye, L. Balents, Z. Wang, S. D. Wilson, and I. Zeljkovic, Rotation symmetry breaking in the normal state of a kagome superconductor KV3Sb5 (2021).
  * Yang  _et al._ [2019] S.-Y. Yang, Y. Wang, B. R. Ortiz, D. Liu, J. Gayles, E. Derunova, R. Gonzalez-Hernandez, L. Smejkal, Y. Chen, S. S. Parkin, _et al._ , Anomalous, anomalous Hall effect in the layered, Kagome, Dirac semimetal KV3Sb5, arXiv preprint arXiv:1912.12288 (2019).
  * Feng  _et al._ [2021] X. Feng, K. Jiang, Z. Wang, and J. Hu, Chiral flux phase in the kagome superconductor av3sb5 (2021), arXiv:2103.07097 [cond-mat.supr-con] .
  * Denner  _et al._ [2021] M. M. Denner, R. Thomale, and T. Neupert, Analysis of charge order in the kagome metal $a$v3sb5 ($a=$k,rb,cs) (2021), arXiv:2103.14045 [cond-mat.str-el] .
  * Lin and Nandkishore [2021] Y.-P. Lin and R. M. Nandkishore, Complex charge density waves at van hove singularity on hexagonal lattices: Haldane-model phase diagram and potential realization in kagome metals $\text{AV}_{3}\text{Sb}_{5}$ (2021), arXiv:2104.02725 [cond-mat.str-el] .
  * Dresselhaus, Mildred S. and Dresselhaus, Gene, Jorio, Ado [2008] Dresselhaus, Mildred S. and Dresselhaus, Gene, Jorio, Ado,  _Group Theory: Application to the Physics of Condensed Matter_ (Springer, 2008).
  * Milnor [1973] J. W. Milnor,  _Morse theory_ , 5th ed., Annals of mathematics studies No. 51 (Princeton Univ. Press, Princeton, NJ, 1973).

## Appendix

##  A Irreducible representations at ${\bm{M}}_{\alpha}$

As shown in Fig. 2, there are two saddle points at ${\bm{M}}_{\alpha}$. The
little co-group at ${\bm{M}}_{\alpha}$ is $D_{\textrm{2h}}$, and the irreps of
the two saddle points at ${\bm{M}}_{1}$ are $\Gamma_{1}^{+}$ and
$\Gamma_{3}^{+}$. The convention used here is the one used in Ref. [59], and
the orientation of the two symmetry axes are shown in Fig. 1.

Now, we explain the one-to-one correspondence between the Brillouin zone
patches, ${\bm{M}}_{\alpha}$, and the vanadium sites, $V_{\alpha}$. Consider
how Bloch states formed from s-orbitals placed on the vanadium sites
($V_{1},V_{2},V_{3}$) transform at ${\bm{M}}_{1}$. A s-orbital Bloch wave on
$V_{1}$ transforms as the irrep $\Gamma_{1}^{+}$. On the other hand,
s-orbitals placed on $V_{2},V_{3}$ form a reducible representation that can be
decomposed into the $\Gamma_{2}^{-}+\Gamma_{4}^{-}$ irreps. Notice that the
irrep of the s-orbital on $V_{1}$ is even under inversion, but the irreps of
the s-orbitals on the $V_{2},V_{3}$ sites are odd under inversion. Since
d-orbitals are even under inversion, this means that Bloch states formed from
d-orbitals at $V_{1}$ must be even under inversion, but Bloch states formed
from d-orbitals on $V_{2},V_{3}$ must be odd under inversion. The saddle
points at ${\bm{M}}_{1}$ are even under inversion; therefore only the
d-orbitals on the $V_{1}$ site contribute to the saddle point at
${\bm{M}}_{1}$. Using a three-fold rotation, it is easy to see that at
${\bm{M}}_{2},{\bm{M}}_{3}$, the contribution to the bands must come from
d-orbitals on the $V_{2},V_{3}$ sites respectively. Hence, there is an one-to-
one correspondence between the Brillouin zone patches, ${\bm{M}}_{\alpha}$,
and the vanadium sites, $V_{\alpha}$. Away from ${\bm{M}}_{1}$, d-orbitals on
$V_{2},V_{3}$ can start mixing with the saddle point bands, but we assume
these contributions are small sufficiently close to ${\bm{M}}_{1}$.

##  B Microscopic mechanism of attractive fermion interaction

Generally, pure electronic interactions, i.e. the screened coulomb repulsion,
give rise to repulsive interactions between patch fermions, which favor SDW
rather than CDW instabilities. Indeed, the RG phase diagram (Fig. 4) suggests
that attraction in the $g^{(0)}_{1}$ or $g^{(0)}_{3}$ channels is likely
necessary to stabilize CDW order. Here, we show that the electron-optical
phonon interaction may renormalize $g^{(0)}_{1}$ or $g^{(0)}_{3}$, and lead to
attraction in these channels. Specifically, we define
$g^{(0)}_{i}=g^{(0)}_{{\rm el},i}+\delta g^{(0)}_{i}$, where $g^{(0)}_{{\rm
el},i}$ is the bare interaction due to screened Coulomb repulsion, and compute
the renormalization $\delta g^{(0)}_{i}$.

Consider the standard coupling of electron density with an optical phonon.
Because $g_{1}$ ($g_{3}$) transfers spin (charge) between saddle points, these
two interactions cannot be generated from coupling to a zone center phonon.
Instead, both $\delta g_{1}^{(0)}$ and $\delta g_{3}^{(0)}$ require coupling
with to a phonon at momentum ${\bm{M}}$. Conversely, $\delta g_{2}^{(0)}$ and
$\delta g_{4}^{(0)}$ are generated from coupling with a zone center phonon. A
recent ab-initio calculation and experiment suggest softening of a breathing
phonon mode at ${\bm{M}}$ [70, 23]. With these considerations in mind, we
focus on the coupling with a phonon at ${\bm{M}}$.

To leading ($0^{\rm th}$) order in spatial derivatives, the electron-phonon
coupling can be expressed as

$\displaystyle\mathcal{H}_{\rm el-oph}$
$\displaystyle=\sum_{{\bm{q}},{\bm{k}},\Gamma,\alpha}g_{\Gamma}u_{\Gamma,-{\bm{k}}}({\bm{M}}_{\alpha})\Psi^{\dagger}_{{\bm{q}}+{\bm{k}}}\Lambda_{\Gamma,\alpha}\Psi_{{\bm{q}}},$
(32)

where $\alpha=1,2,3$ labels the momentum ${\bm{M}}_{\alpha}$, $\Gamma$ labels
the irreducible representation (irrep) of the little group at
${\bm{M}}_{\alpha}$. In this notation, $u_{\Gamma}({\bm{M}}_{\alpha})$ is the
optical phonon mode near ${\bm{M}}_{\alpha}$ in irrep $\Gamma$,
$\Lambda_{\Gamma,\alpha}$ is obtained such that the fermion bilinear
$\Psi^{\dagger}_{{\bm{q}}}\Lambda_{\Gamma,a}\Psi_{{\bm{q}}}$ transforms as
$\Gamma$ irrep at ${\bm{M}}_{\alpha}$. Note that in Eq. (32), the momentum
${\bm{k}}$ is small, i.e. within the continuum model, and the large momentum
transfer ${\bm{M}}_{\alpha}$ is treated via the flavor indices of the fermion
fields.

Without loss of generality, we consider ${\bm{M}}_{1}$, whose little group is
${\rm D}_{2h}=\langle C_{2}^{\prime},C_{2}^{\prime\prime},i\rangle$ (see Fig.
1), where $\langle\cdot\rangle$ indicates the group generated by “$\cdot$”
operations. Here, $C_{2}^{\prime}$ and $C_{2}^{\prime\prime}$ are two-fold
rotation along y-axis and x-axis, respectively. $i$ denotes inversion. Using
the transformation of Bloch wave functions,

$\displaystyle\mathcal{U}_{\mathcal{O}}\psi_{n,{\bm{k}}}({\bm{r}})$
$\displaystyle=\mathcal{U}_{\mathcal{O}}\mathrm{e}^{\mathrm{i}{\bm{q}}\cdot{\bm{r}}}{\rm
u}_{n,{\bm{q}}}({\bm{r}})=\mathrm{e}^{\mathrm{i}{\bm{q}}\cdot\mathcal{O}{\bm{r}}}{\rm
u}_{n,{\bm{q}}}(\mathcal{O}{\bm{r}})$
$\displaystyle=\mathrm{e}^{\mathrm{i}\mathcal{O}^{-1}{\bm{q}}\cdot{\bm{r}}}{\rm
u}_{n,\mathcal{O}^{-1}{\bm{q}}}({\bm{r}})=\psi_{n,\mathcal{O}^{-1}{\bm{q}}}({\bm{r}}),$
(33)

where $\mathcal{O}\in{\rm D}_{2h}$, we find that the fermion bilinear at
momentum ${\bm{M}}_{1}$ from the patch fermions should be in either the
$\Gamma_{1}^{+}$ or $\Gamma_{3}^{+}$ irrep of ${\rm D}_{2h}$, with,

$\displaystyle\Lambda_{\Gamma_{1}^{+}}=\begin{pmatrix}0&0&0\\\ 0&0&1\\\
0&1&0\end{pmatrix},\Lambda_{\Gamma_{3}^{+}}=\begin{pmatrix}0&0&0\\\
0&0&\mathrm{i}\\\ 0&-\mathrm{i}&0\end{pmatrix}.$ (34)

Similarly, the representation of optical phonon modes at ${\bm{M}}_{1}$ in
irrep $\Gamma_{1}^{+},\Gamma_{3}^{+}$ can be identified using [76],

$\displaystyle\Gamma_{\rm oph}=\Gamma_{\rm equiv}\times\Gamma_{\rm vec},$ (35)

where $\Gamma_{\rm equiv}$ is the equivalence representation for the kagomé
lattice sites at ${\bm{M}}$, $\Gamma_{\rm vec}$ is the representation for
lattice displacement, the same as that for a 3d vector. The $\Gamma_{\rm
oph,1}^{+}$ irrep comes from either $\Gamma_{\rm equiv,2}^{-}\times\Gamma_{\rm
vec,2}^{-}$ or $\Gamma_{\rm equiv,4}^{-}\times\Gamma_{\rm vec,4}^{-}$ optical
phonon modes. While the $\Gamma_{\rm oph,3}^{+}$ irrep comes from either
$\Gamma_{\rm equiv,2}^{-}\times\Gamma_{\rm vec,4}^{-}$ or $\Gamma_{\rm
equiv,4}^{-}\times\Gamma_{\rm vec,2}^{-}$ optical phonon modes. For example,
both the “star of David” and the “inverse star of David” breathing phonon mode
studied in Ref. [70] comes from the $\Gamma_{\rm oph,1}^{+}=\Gamma_{\rm
equiv,2}^{-}\times\Gamma_{\rm vec,2}^{-}$ optical phonon mode related by
reversing the sign of the lattice displacement field.

The effective four-fermion interaction renormalized by $\mathcal{H}_{\rm el-
oph}$ reads,

$\displaystyle\delta\mathcal{H}_{\rm
eff}=\sum_{{\bm{q}},{\bm{q}}^{\prime},\Gamma,\alpha}\mathcal{V}_{\Gamma}\left(\Psi^{\dagger}_{{\bm{q}}}\Lambda_{\Gamma,\alpha}\Psi_{{\bm{q}}}\right)\left(\Psi^{\dagger}_{{\bm{q}}^{\prime}}\Lambda_{\Gamma,\alpha}\Psi_{{\bm{q}}^{\prime}}\right),$
(36)

where
$\mathcal{V}_{\Gamma}=g_{\Gamma}^{2}\mathcal{D}_{\Gamma}({\bm{M}}_{1},\Omega_{n}=0)=-\frac{2g^{2}_{\Gamma}}{\omega_{\Gamma,{\bm{M}}_{1}}}<0$,
with $\omega_{\Gamma,{\bm{k}}}$ the phonon spectrum,
$\mathcal{D}_{\Gamma}({\bm{k}},\Omega_{n})=\frac{2\omega_{\Gamma,{\bm{k}}}}{(i\Omega_{n})^{2}-\omega_{\Gamma,{\bm{k}}}^{2}}$
the phonon propagator. This gives finally,

$\displaystyle\delta g_{1}^{(0)}$
$\displaystyle=\mathcal{V}_{\Gamma_{1}^{+}}+\mathcal{V}_{\Gamma_{3}^{+}},$
$\displaystyle\delta g_{3}^{(0)}$
$\displaystyle=\mathcal{V}_{\Gamma_{1}^{+}}-\mathcal{V}_{\Gamma_{3}^{+}}.$
(37)

It is straightforward to see that $\delta g_{1}^{(0)}<0$ while the sign of
$\delta g_{3}^{(0)}$ depends on the relative strength of
$\mathcal{V}_{\Gamma_{1}^{+}}$ and $\mathcal{V}_{\Gamma_{3}^{+}}$. As a
result, all three fixed point solutions may be approached when the electron-
phonon interaction is taken into account.

Experimentally, one may test for the significance of electron-phonon coupling
by probing the effect on the phonons. In particular, the coupling
$\mathcal{H}_{\rm el-oph}$ also renormalizes the phonon spectrum, and lowers
the optical phonon gap. Indeed, the phonon gap for
$u_{\Gamma_{1}^{+}},u_{\Gamma_{3}^{+}}$ are lowered by $\delta m_{\rm oph}\sim
g_{\Gamma}^{2}\Pi_{\rm ph}({\bm{M}})$, which may lead to a Peierls structural
transition when the renormalized phonon spectrum becomes gapless at
${\bm{M}}$.

Does the CDW order observed in the kagomé materials come from electronic or
structural instabilities? While the full answer to this question requires
knowledge of the microscopic details and self-consistent analysis of the
electron-phonon coupled system, we note that the four-fermion interaction can
be significantly enhanced in the RG flow near the fixed point. Moreover, in
the weak coupling limit, as $g^{(0)}_{{\rm el},i}$ is weak, a moderate
electron-phonon interaction could alter the sign of the bare fermion
interaction, and lead to CDW ordering due to an electronic instability.

##  C Solutions of the rCDW mean-field theory

Here, we discuss the possible solutions of the rCDW mean-field theory. The
free energy of the system as a function of the rCDW vector,
$\bm{N}=(N_{1},N_{2},N_{3})$, has a tetrahedral point group symmetry,
$T_{\textrm{d}}$. For fixed $\absolutevalue{\bm{N}}=N_{0}$, rCDW vectors can
be categorized into different classes corresponding to the Wyckoff positions
of $T_{\textrm{d}}$. In particular, there are three distinct rCDW classes that
correspond to high symmetry points where the rCDW free energy must have a
critical point on the $\absolutevalue{\bm{N}}=N_{0}$ surface. They are the
3Q$+$, 3Q$-$, and 1Q classes defined as classes of rCDW states with the
directions,

3Q$+$: |  $\\{(111),(1\bar{1}\bar{1}),(\bar{1}1\bar{1}),(\bar{1}\bar{1}1)\\}$,
---|---
3Q$-$: |  $\\{(\bar{1}\bar{1}\bar{1}),(\bar{1}11),(1\bar{1}1),(11\bar{1})\\}$,
1Q: |  $\\{(\pm 100),(0\pm 10),(00\pm 1)\\}$.

Configurations in each class are related by $T_{d}$ symmetry. The critical
points of the 3Q$\pm$ class must correspond to either a maxima or minima. This
is because these directions coincide with the three-fold symmetry axes of
$T_{d}$, so the curvature of the free energy around these points must have the
same sign in all directions. The same cannot be said for the 1Q class which
coincides with the two-fold axes of $T_{d}$, so points in the 1Q directions
can either all be saddle points or maxima/minima.

For $\absolutevalue{\bm{N}}=N_{0}$, the free energy is a smooth function
defined on a sphere. Morse theory tells us that for such a function, the
critical points must satisfy the following constraint, i.e. the number of
maxima $(n_{\textrm{max}})$, minima $(n_{\textrm{min}})$, and saddle points
$(n_{\textrm{sp}})$ must satisfy,

$n_{\textrm{max}}+n_{\textrm{min}}-n_{\textrm{sp}}=\chi=2,$ (38)

where $\chi=2$ is the Euler characteristic of the sphere [77]. If we assume
that the 3Q$\pm$ and 1Q classes are the only critical points of the free
energy, Eq. (38) is satisfied if the 1Q class corresponds to a saddle point
since there are 8 directions in the 3Q$\pm$ classes and 6 directions in the 1Q
class. Other solutions to Eq. (38) require additional critical points on the
sphere. Points on the sphere other than the high symmetry points listed above
can be categorized using Wyckoff positions. For $T_{d}$, there are two types
of low symmetry Wyckoff positions which we label $c,d$ with degeneracies
$n_{c}=12,n_{d}=24$. It is straightforward to see that the simplest way to
satisfy Eq. (38) and include a low symmetry critical point is if there is a
saddle point at the Wyckoff positions $c$ or alternatively, if there is a
minimum/maximum at Wyckoff position $c$ and a saddle point at Wyckoff position
$d$. In addition to the cases listed above, pairs of saddle points and
maxima/minima can be added to both types of low symmetry Wyckoff positions
while still satisfying Eq. (38). However, we assume that cases with a large
degeneracy of critical points are not physically relevant since they can only
come from high-order terms with large powers of $N$ that are irrelevant at low
energies, and correspond to highly oscillatory functions on the sphere.
Therefore, given this assumption, the cases where the maxima/minima are only
at the 3Q$\pm$ or 1Q positions are expected and the solution to the rCDW mean-
field theory should be in of those directions.

A similar analysis can be extended to other order parameters. The iCDW and
rSDW order parameters have a cubic point group symmetry, $O_{h}$, because they
are odd under time-reversal symmetry unlike the rCDW order parameters.
Therefore, the rCDW states in the 3Q$+$ and 3Q$-$ classes are related under
time-reversal symmetry. Hence, for iCDW and rSDW, the the 3Q$+$ and 3Q$-$
classes are subsets of the larger 3Q class which is defined as the union of
the 3Q$\pm$ classes.

##  D Landau theory coefficients

The coefficients, $K_{1},K_{2},K_{3},K_{4}$, of the Landau theory defined in
Eq. (12) are

$\displaystyle K_{1}(T,\mu)=$
$\displaystyle\frac{1}{2\Lambda^{2}}\int_{-\Lambda}^{\Lambda}\frac{\differential^{2}q}{(2\pi)^{2}}\frac{f(\varepsilon_{2}(\bm{q})-\mu)-f(\varepsilon_{3}(\bm{q})-\mu)}{\varepsilon_{2}(\bm{q})-\varepsilon_{3}(\bm{q})},$
(39) $\displaystyle K_{2}(T,\mu)=$
$\displaystyle-\frac{1}{2\Lambda^{2}}\int_{-\Lambda}^{\Lambda}\frac{\differential^{2}q}{(2\pi)^{2}}\bigg{(}\frac{f(\varepsilon_{1}(\bm{q})-\mu)}{(\varepsilon_{1}(\bm{q})-\varepsilon_{2}(\bm{q}))(\varepsilon_{1}(\bm{q})-\varepsilon_{3}(\bm{q}))}$
$\displaystyle\qquad+\frac{f(\varepsilon_{2}(\bm{q})-\mu)}{(\varepsilon_{2}(\bm{q})-\varepsilon_{1}(\bm{q}))(\varepsilon_{2}(\bm{q})-\varepsilon_{3}(\bm{q}))}+\frac{f(\varepsilon_{3}(\bm{q})-\mu)}{(\varepsilon_{3}(\bm{q})-\varepsilon_{1}(\bm{q}))(\varepsilon_{3}(\bm{q})-\varepsilon_{2}(\bm{q}))}\bigg{)},$
(40) $\displaystyle K_{3}(T,\mu)=$
$\displaystyle\frac{1}{24\Lambda^{2}}\int_{-\Lambda}^{\Lambda}\frac{\differential^{2}q}{(2\pi)^{2}}\bigg{(}\frac{f^{\prime}(\varepsilon_{1}(\bm{q})-\mu)}{(\varepsilon_{1}(\bm{q})-\varepsilon_{2}(\bm{q}))(\varepsilon_{1}(\bm{q})-\varepsilon_{3}(\bm{q}))}$
$\displaystyle\qquad+\frac{f^{\prime}(\varepsilon_{2}(\bm{q})-\mu)}{(\varepsilon_{2}(\bm{q})-\varepsilon_{1}(\bm{q}))(\varepsilon_{2}(\bm{q})-\varepsilon_{3}(\bm{q}))}+\frac{f^{\prime}(\varepsilon_{3}(\bm{q})-\mu)}{(\varepsilon_{3}(\bm{q})-\varepsilon_{1}(\bm{q}))(\varepsilon_{3}(\bm{q})-\varepsilon_{2}(\bm{q}))}\bigg{)},$
(41) $\displaystyle K_{4}(T,\mu)=$
$\displaystyle-\frac{1}{8\Lambda^{2}}\int_{-\Lambda}^{\Lambda}\frac{\differential^{2}q}{(2\pi)^{2}}\bigg{(}\frac{f(\varepsilon_{2}(\bm{q})-\mu)-f(\varepsilon_{3}(\bm{q})-\mu)}{(\varepsilon_{2}(\bm{q})-\varepsilon_{3}(\bm{q}))^{3}}-\frac{f^{\prime}(\varepsilon_{2}(\bm{q})-\mu)+f^{\prime}(\varepsilon_{3}(\bm{q})-\mu)}{2(\varepsilon_{2}(\bm{q})-\varepsilon_{3}(\bm{q}))^{2}}\bigg{)},$
(42)

where $f(\varepsilon)=1/(\exp(\varepsilon/T)+1)$ is the Fermi-Dirac function.
These integrals can be evaluated asymptotically in the limit $\mu\,,k_{B}T\ll
a\Lambda^{2}\,,b\Lambda^{2}$. If we assume perfect nesting, the patch
dispersions take the form,

$\displaystyle\varepsilon_{1}(\bm{q})=$
$\displaystyle\frac{3t}{4}\left(3q_{x}^{2}-q_{y}^{2}\right),$ (43)
$\displaystyle\varepsilon_{2}(\bm{q})=$
$\displaystyle\frac{3t}{4}2q_{y}\left(q_{y}+\sqrt{3}q_{x}\right),$ (44)
$\displaystyle\varepsilon_{3}(\bm{q})=$
$\displaystyle\frac{3t}{4}2q_{y}\left(q_{y}-\sqrt{3}q_{x}\right),$ (45)

where $t>0$. Given this definition, the asymptotic expressions of
$K_{2},K_{3},K_{4}$ in the limit $\mu,k_{B}T\ll t\Lambda^{2}$ are

$\displaystyle K_{2}\approx$
$\displaystyle-\frac{16}{\pi^{2}t\Lambda^{2}k_{B}T}H_{2}(\mu/k_{B}T),$ (46)
$\displaystyle K_{3}\approx$
$\displaystyle\frac{8}{3\pi^{2}\Lambda^{2}(k_{B}T)^{2}}H_{3}(\mu/k_{B}T),$
(47) $\displaystyle K_{4}\approx$
$\displaystyle\frac{1}{12\sqrt{3}\pi^{2}t\Lambda^{2}(k_{B}T)^{2}}H_{4}(\mu/k_{B}T)\ln(t\Lambda^{2}/k_{B}T),$
(48)

where $H_{2},H_{3},H_{4}$ are the integral functions,

$\displaystyle H_{2}(z)=$ $\displaystyle\int_{0}^{\infty}\differential
x\int_{0}^{x/\sqrt{3}}\differential
y\bigg{(}\frac{F(\tilde{\varepsilon}_{1}(x,y)-z)}{(\tilde{\varepsilon}_{1}(x,y)-\tilde{\varepsilon}_{2}(x,y))(\tilde{\varepsilon}_{1}(x,y)-\tilde{\varepsilon}_{3}(x,y))}$
$\displaystyle-\frac{F(\tilde{\varepsilon}_{2}(x,y)-z)}{(\tilde{\varepsilon}_{1}(x,y)-\tilde{\varepsilon}_{2}(x,y))(\tilde{\varepsilon}_{2}(x,y)-\tilde{\varepsilon}_{3}(x,y))}+\frac{F(\tilde{\varepsilon}_{3}(x,y)-z)}{(\tilde{\varepsilon}_{1}(x,y)-\tilde{\varepsilon}_{3}(x,y))(\tilde{\varepsilon}_{2}(x,y)-\tilde{\varepsilon}_{3}(x,y))}\bigg{)},$
(49) $\displaystyle H_{3}(z)=$ $\displaystyle\int_{0}^{\infty}\differential
x\int_{0}^{x/\sqrt{3}}\differential
y\bigg{(}\frac{F^{\prime}(\tilde{\varepsilon}_{1}(x,y)-z)}{(\tilde{\varepsilon}_{1}(x,y)-\tilde{\varepsilon}_{2}(x,y))(\tilde{\varepsilon}_{1}(x,y)-\tilde{\varepsilon}_{3}(x,y))}$
$\displaystyle-\frac{F^{\prime}(\tilde{\varepsilon}_{2}(x,y)-z)}{(\tilde{\varepsilon}_{1}(x,y)-\tilde{\varepsilon}_{2}(x,y))(\tilde{\varepsilon}_{2}(x,y)-\tilde{\varepsilon}_{3}(x,y))}+\frac{F^{\prime}(\tilde{\varepsilon}_{3}(x,y)-z)}{(\tilde{\varepsilon}_{1}(x,y)-\tilde{\varepsilon}_{3}(x,y))(\tilde{\varepsilon}_{2}(x,y)-\tilde{\varepsilon}_{3}(x,y))}\bigg{)},$
(50) $\displaystyle H_{4}(z)=$ $\displaystyle\int_{0}^{\infty}\differential
v\left(\frac{-F(v-z)+F(-v-z)}{v^{3}}+\frac{F^{\prime}(v-z)+F^{\prime}(-v-z)}{z^{2}}\right).$
(51)

In the equations above, $F(x)=1/(e^{x}+1)$, and
$\tilde{\varepsilon}_{1}(x,y)=3x^{2}-y^{2}\,,\tilde{\varepsilon}_{2}(x,y)=2y(y+\sqrt{3}x)\,,\tilde{\varepsilon}_{3}(x,y)=2y(y-\sqrt{3}x).$
The expression holds for arbitrary ratios of $\mu/k_{B}T$. The values of the
$H_{2},H_{3},H_{4}$ when $\mu\ll k_{B}T$ are

$\displaystyle H_{2}(0)\approx$ $\displaystyle 0.039571,$ (52) $\displaystyle
H_{3}(0)\approx$ $\displaystyle 0.014216,$ (53) $\displaystyle
H_{4}(0)\approx$ $\displaystyle 0.213139.$ (54)

$H_{2},H_{3},H_{4}$ each change sign once at the following values:

$\displaystyle z_{0,2}\approx$ $\displaystyle 2.14268,$ (55) $\displaystyle
z_{0,3}\approx$ $\displaystyle 4.05161,$ (56) $\displaystyle z_{0,4}\approx$
$\displaystyle 1.91067.$ (57)

The expression for $K_{1}$ can be evaluated asymptotically when $\mu\ll
k_{B}T$:

$K_{1}\approx-\frac{1}{6\sqrt{3}\pi^{2}t\Lambda^{2}}\left(\ln(t\Lambda^{2}/k_{B}T)^{2}+2\ln
3\ln(t\Lambda^{2}/k_{B}T)\right).$ (58)

For $\mu\ll T$, $K_{3}-2K_{4}$ which determines whether the quartic term of
the rCDW Landau theory prefers a 3Q$+$ or 3Q$-$ configuration is

$K_{3}-2K_{4}\approx\frac{1}{\pi^{2}t\Lambda^{2}(k_{B}T)^{2}}\left(\frac{8H_{3}(0)}{3}-\frac{H_{4}(0)}{12\sqrt{3}}\ln(t\Lambda^{2}/k_{B}T)\right).$
(59)

Therefore, $K_{3}-2K_{4}<0$ when

$\frac{\mu}{t\Lambda^{2}}\ll\frac{k_{B}T}{t\Lambda^{2}}<\exp(-32\sqrt{3}\frac{H_{3}(0)}{H_{4}(0)})\approx
0.0248.$ (60)

##  E Derivation of Eq. (22): the current operator

Here, we show more details to obtain the bond current operator, which is
defined from the continuity equation of charge density. If we define the
charge density at site ${\bm{r}}_{\alpha}$ as
$\hat{n}_{{\bm{r}}_{\alpha}}=-e\,d^{\dagger}_{{\bm{r}}_{\alpha}}d_{{\bm{r}}_{\alpha}}$,
the lattice version of the continuity equation reads,

$\displaystyle\frac{\mathop{}\\!\mathrm{d}\hat{n}_{{\bm{r}}_{\alpha}}}{\mathop{}\\!\mathrm{d}\rm
t}+\sum_{{\bm{r}}_{\alpha}\rightarrow{\bm{r}}^{\prime}_{\alpha^{\prime}}}\hat{j}_{{\bm{r}}_{\alpha}{\bm{r}}^{\prime}_{\alpha^{\prime}}}=0.$
(61)

The expression for $\hat{j}$ can be obtained from the equation of motion of
$\hat{n}_{{\bm{r}}_{\alpha}}$ which reads,

$\displaystyle\frac{\mathop{}\\!\mathrm{d}\hat{n}_{{\bm{r}}_{\alpha}}}{\mathop{}\\!\mathrm{d}\rm
t}=-\frac{\mathrm{i}}{\hbar}\left[\hat{n}_{{\bm{r}}_{\alpha}},\hat{H}\right]=-\frac{\mathrm{i}}{\hbar}\left[\hat{n}_{{\bm{r}}_{\alpha}},\hat{H}_{0}\right]+...$
(62)

where “…” denotes contribution from
$\left[\hat{n}_{{\bm{r}}_{\alpha}},\hat{H}_{1}\right]$, which only possibly
contribute at $\mathcal{O}(\phi^{2})$ to the bond-current order parameter
$\langle\hat{j}\rangle$. On the other hand, as only odd power in $\phi$ term
breaks time reversal symmetry, we conclude that $\hat{H}_{1}$ should not
contribute to the bond-current order parameter at the mean field level.
Hereafter, only the contribution from $\hat{H}_{0}$ is considered. The NN bond
order requires NN hopping term in $\hat{H}_{0}$. Noting,

$\displaystyle\left[\hat{n}_{{\bm{r}}_{\alpha}},\hat{H}_{0}\right]$
$\displaystyle=\left[\hat{n}_{{\bm{r}}_{\alpha}},t_{1}\sum_{\langle{\bm{r}}^{\prime}_{\alpha^{\prime}},{\bm{r}}^{\prime}_{\alpha^{\prime}}+{\bm{e}}_{\beta}\rangle}d^{\dagger}_{{\bm{r}}^{\prime}_{\alpha^{\prime}}}d_{{\bm{r}}^{\prime}_{\alpha^{\prime}}+{\bm{e}}_{\beta}}+d^{\dagger}_{{\bm{r}}^{\prime}_{\alpha^{\prime}}+{\bm{e}}_{\beta}}d_{{\bm{r}}^{\prime}_{\alpha^{\prime}}}\right]$
$\displaystyle=\left[\hat{n}_{{\bm{r}}_{\alpha}},t_{1}\sum_{\alpha\neq\beta}d^{\dagger}_{{\bm{r}}_{\alpha}}d_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}+d^{\dagger}_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}d_{{\bm{r}}_{\alpha}}\right]$
$\displaystyle=-e\,t_{1}\sum_{\alpha\neq\beta}\left(d^{\dagger}_{{\bm{r}}_{\alpha}}d_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}-d^{\dagger}_{{\bm{r}}_{\alpha}+{\bm{e}}_{\beta}}d_{{\bm{r}}_{\alpha}}\right),$
(63)

where the second line comes from the fact that the commutator is nonzero only
for bond satisfying ${\bm{r}}_{\alpha}={\bm{r}}^{\prime}_{\alpha^{\prime}}$ or
${\bm{r}}_{\alpha}={\bm{r}}^{\prime}_{\alpha^{\prime}}+{\bm{e}}_{\beta}$, we
find Eq. (22) in the main text.