# Sublattice Interference promotes Pair Density Wave order in Kagome Metals

Yi-Ming Wu  Stanford Institute for Theoretical Physics, Stanford University,
Stanford, California 94305, USA  Ronny Thomale  Institute for Theoretical
Physics and Astrophysics, University of Würzburg, D-97074 Würzburg, Germany
Department of Physics and Quantum Centers in Diamond and Emerging Materials
(QuCenDiEM) group, Indian Institute of Technology Madras, Chennai 600036,
India  S. Raghu  Stanford Institute for Theoretical Physics, Stanford
University, Stanford, California 94305, USA

###### Abstract

Motivated by the observation of a pair density wave (PDW) in the kagome metal
CsV3Sb5, we consider the fate of electrons near a p-type van Hove singularity
(vHS) in the presence of local repulsive interactions. We study the effect of
such interactions on Fermi surface “patches” at the vHS. We show how a feature
unique to the Kagome lattice known as sublattice interference crucially
affects the form of the interactions among the patches. The renormalization
group (RG) flow of such interactions results in a regime where the nearest
neighbor interaction $V$ exceed the onsite repulsion $U$. We identify this
condition as being favorable for the formation of charge-density-wave (CDW)
and PDW orders. In the weak coupling limit, we find a complex CDW order as the
leading instability, which breaks time reversal symmetry. Beyond RG, we
perform a Hartree-Fock study to a $V$-only model and find the pair-density-
wave order indeed sets in at some intermediate coupling.

### Introduction.$-$

Spins and electrons on the Kagome lattice [Fig.1(a)] have long been studied
due to their potential for exhibiting a panoply of exotic phases of matter.
Insulating kagome compounds, for instance, are among the most prominent
candidate spin liquid materials Han2012 , and insulating phases with non-
trivial topology have also been studied on the Kagome lattice joe-Fe3Sn2 ;
Kang2020 ; Yin2020 . With the discovery of a family of Kagome metals AV3Sb5
(A=K,Cs,Rb), a new wave of excitement has been elicited by the prospects for
intriguing density wave and superconducting ground states in these systems
PhysRevMaterials.3.094407 ; jiangii ; neupert-review .

Among the more fascinating observed phenomena in these kagome metals include
new bragg-like peaks inside the superconducting phase roton-pdw , a hallmark
of an exotic superconducting order known as a pair density wave (PDW)
doi:10.1146/annurev-conmatphys-031119-050711 ; Fradkin2015 ; Himeda2002 ;
Yang_2009 ; Raczkowski2007 ; Berg2007 ; Capello2008 ; PALee2014 ; Wang2015a ;
Wang2015b ; Edkins2019 ; Wang2018 ; Zhengzhi ; tTMD ; Shaffer1 ; Shaffer2
PDWs are superconductors with an order parameter that varies periodically in
space. Additionally, fascinating effects have been discovered in the
superconducting fluctuation spectrum, including nearly condensed excited
states with charge 4e, 6e superconducting fluctuations ge2022 . These
observations call for a greater theoretical scrutiny, and invite us to make
predictions for electronic phases of kagome metals.

In this Letter, motivated by these recent developments, we address the issue
of whether the PDW superconducting phase can in principle arise on the Kagome
lattice. Since such superconductivity requires analysis of the intermediate
coupling problem, robust pairing mechanisms for PDW formation have only
recently been uncovered WuPDWlargeN . One key requirement for PDW order is the
presence of strong repulsive interactions with somewhat suppressed onsite
interactions. While this is rather unusual in most solids, we show here that
this requirement is met when the chemical potential crosses a so-called p-type
(for “pure”) van Hove singularity (vHS) PhysRevB.86.121105 ;
PhysRevLett.127.177001 ; kang ; hu2021rich  [see Fig.1(b)].

Figure 1: Sublattice interference in the Kagome lattice. (a) Each unit cell
has $\alpha=\,$A,B,C different sublattice sites. (b) Band structure of the
tight binding model. We focus on the middle band, highlighted with red color,
which exhibits a p-type van Hove singularity. (c) For this band, the
transformation matrix $u_{\alpha}({\bm{k}})$, defined as
$c_{\alpha}({\bm{k}})=u_{\alpha}({\bm{k}})\psi_{\alpha}({\bm{k}})$ where
$c_{\alpha}$ and $\psi_{\alpha}$ are lattice and band fermions respectively,
has a modulated distribution in the Brillouin zone.

In an important theoretical study, one of us showed PhysRevB.86.121105  that
precisely at such a p-type vHS, there is the phenomena of sublattice
interference (SI) - where each of the 3 distinct sublattices has non-zero
support only on one of the 3 distinct vH points. To show this SI crucially
determines the low energy effective interactions in the system, we construct a
renormalization group (RG) theory based on the p-type vHS and show that the
onsite repulsion runs towards weak coupling, while nearest neighbor
interactions grow under the RG. This peculiarity results in a rich phase
diagram including time-reversal symmetry breaking charge-density-wave (CDW)
orders and various uniform superconductivity in the weak coupling limit; while
in the strong coupling regime it makes the Kagome system poised to exhibit the
PDW order. While more experiments are needed to characterize precisely the
normal state in this system, our analysis already establishes that the kagome
motif once again provides us with an avenue towards an exotic phase of matter,
in this case the PDW superconductor.

SI and the triviality of Hubbard $U$.$-$ To determine the fate of effective
interactions at p-type van Hove filling, we utilize parquet RG methods and
restrict attention to Fermi surface “patches” in the neighborhood of the
distinct van Hove points Dzyaloshinskii1987 ; Schulz1987  (similar RG analysis
with two-fold of van Hove singularities is considered in Ref.Scammell2023 ).
Defining the fermion destruction operators in patch $\alpha$ as
$\psi_{\alpha}$, we encode real-space 4-fermion interactions as intra- and
interpatch couplings after Fourier transformation. In a crystal with time-
reversal and/or inversion symmetry, such interactions take the form
$H_{I}=g_{1}\psi_{\alpha}^{\dagger}\psi_{\beta}^{\dagger}\psi_{\alpha}\psi_{\beta}+g_{2}\psi_{\alpha}^{\dagger}\psi_{\beta}^{\dagger}\psi_{\beta}\psi_{\alpha}+g_{3}\psi_{\alpha}^{\dagger}\psi_{\alpha}^{\dagger}\psi_{\beta}\psi_{\beta}+g_{4}\psi_{\alpha}^{\dagger}\psi_{\alpha}^{\dagger}\psi_{\alpha}\psi_{\alpha}$
where $\alpha\neq\beta$ and the momentum summation and a spin configuration
$\sigma\sigma^{\prime}\sigma^{\prime}\sigma$ is assumed. Patch models have
been applied to the square latticeFurukawa , the honeycomb lattice
Nandkishore2012 ; PhysRevB.89.144501  and moiré systemsIsobe2018 ;
PhysRevB.100.085136 ; Zhengzhi  to capture interaction-driven electronic
orders. In both cases, all the $g_{i}$’s are set by the largest Hubbard onsite
interaction $H_{U}=U\sum_{i}n_{i\uparrow}n_{i\downarrow}$.

However, adopting the same strategy on the Kagome lattice near a p-type vHS,
one finds that the Hubbard interaction U contributes only to the $g_{4}$. This
is due to SI: Near a p-type vHS, each of the 3 sublattices (A,B,C) has non-
zero support only on one of the 3 distinct van Hove points
$(\bm{M}_{\alpha},\alpha=A,B,C)$ which is visible through the sublattice
weight in Fig.1(c). The SI does not result from fine tuning, as it exists even
if longer range hoppings are considered supp . Thus, when only $H_{U}$ is
present for electrons near a p-type vHS, $g_{1}=g_{2}=g_{3}=0$, and
$g_{4}=U/t$. With this choice of bare couplings, $g_{4}$ weakens under RG
flow, eventually reaching a trivial fixed point with $g_{4}^{*}=0$ supp .
Thus, due to SI, the repulsive Hubbard U is irrelevant in the RG sense near
the p-type vHS of the Kagome lattice 111For comparison, one can also consider
the alternative kagome van Hove filling which is of m-type [see Fig.1(b)].
There, each $\psi_{\alpha}$ is a superposition of two distinct lattice
fermions and it can be directly verified that the Hubbard $U$ then has equal
contributions to all $g_{i}$..

Extended interactions and the 6-patch theory.$-$ Due to the apparent
irrelevance of $U$, it is necessary to include at least nearest neighbor
interactions. The most important such term is the nearest neighbor density-
density repulsion $H_{V}=V\sum_{\braket{ij}}n_{i}n_{j}.$ With both $H_{U}$ and
$H_{V}$, the 3-patch model described above is no longer adequate since nearest
neighbor interactions contain momentum dependence, which differentiates the
patches at $\pm\bm{M}_{\alpha}$. Therefore, we need to consider the full
6-patch theory, for which there are 16 different symmetry allowed interactions
in total, as shown in Fig.2(a). Here we adopt the convention used in
Ref.Isobe2018  and define the patch fermions as
$\psi_{\alpha+}=\psi({\bm{k}}\to\bm{M}_{\alpha})$ and
$\psi_{\alpha-}=\psi({\bm{k}}\to-\bm{M}_{\alpha})$. The interactions in this
six patch model can be written as

$H_{I}=\sum_{i,j=1}^{4}\sum_{\left\\{\alpha_{j},\tau_{j}\right\\}}g_{ij}\psi^{\dagger}_{\alpha_{1}\tau_{1}}\psi^{\dagger}_{\alpha_{2}\tau_{2}}\psi_{\alpha_{3}\tau_{3}}\psi_{\alpha_{4}\tau_{4}},$
(1)

with $g_{ij}$ being the dimensionless interaction strengths. Momentum
conservation constrains the indices as follows: the patch indices satisfy
$\alpha_{1}=\alpha_{3}\neq\alpha_{2}=\alpha_{4}$ for $i=1$,
$\alpha_{1}=\alpha_{4}\neq\alpha_{2}=\alpha_{3}$ for $i=2$,
$\alpha_{1}=\alpha_{2}\neq\alpha_{3}=\alpha_{4}$ for $i=3$, and
$\alpha_{1}=\alpha_{2}=\alpha_{3}=\alpha_{4}$ for $i=4$. The ‘valley’ indices
$\tau_{i}=\pm$, labeling whether the patch is at $\bm{M}_{\alpha}$ or
$-\bm{M}_{\alpha}$, obey the same rule associated with $j$. That is,
$\tau_{1}=\tau_{3}\neq\tau_{2}=\tau_{4}$ for $j=1$ et cetera.

Figure 2: (a) All symmetry allowed interactions in the six patch model.(b) RG
flow for the interactions and the order parameter vertex $\Gamma$ with the SI
effect. We set $U(0)=0.8t$ and $V(0)=0.3t$ and starting from some intermediate
energy scale we have $V\gg U$. The resulting weak coupling instabilities are
degenerate complex CDW- and imaginary SDW+/-. (c) The same RG flow but without
SI effect. In this case all the $g_{ij,0}$ are set by the largest Hubbard $U$,
and there does not exist a constant map $U,V$ and $g_{ij}(y)$. The leading
weak coupling instability is the $d$-wave SCNandkishore2012 ;
PhysRevB.89.144501 .

Once again, SI crucially influences the initial conditions of the RG flows. A
straightforward calculation supp  shows that $U$ contributes only to $g_{4j}$
\- a direct generalization of the 3-patch theory - while $V$ only contributes
to $g_{2j}$:

$g_{4j,0}={\textstyle{U\over
t}};g_{22,0}=g_{24,0}=-g_{21,0}=-g_{23,0}={\textstyle{2V\over t}}.$ (2)

The subscript $0$ above denotes bare interactions before running RG. From
microscopics, it is natural to expect that $U$ (and therefore $g_{4j,0}$)
should be the largest.

CDW at weak coupling.$-$ We first investigate the weak coupling limit. In this
case the fate of the fermions can be described by the one loop RG equations of
$g_{ij}$, for which we keep the most divergent “log squared” terms in
perturbation theory, i.e. the particle-hole bubble $\Pi_{ph}^{(0)}$ at momenta
$\bm{M}_{\alpha}$ and the particle-particle bubble $\Pi_{pp}^{(0)}$ at zero
momentum, and set the running parameter as the latter: $y=\Pi_{pp}^{(0)}(0)$
Isobe2018 ; supp . With this convention, the RG equations for the interactions
take the form $dg_{ij}/dy=g_{mn}R_{mn,kl}g_{kl}$. It’s easy to see that close
to the critical value $y_{c}$ all $g_{ij}$ behaves in a similar manner, namely
$g_{ij}=G_{ij}/(y_{c}-y)$. If some $G_{ij}$ is nonzero, the corresponding
$g_{ij}$ thus diverges at $y_{c}$. A typical order parameter flows as
$d\Delta_{i}/dy=d\Gamma_{i}\Delta_{i}$ under renormalization, where $d=1$ for
superconducting orders and $d=d_{1}\leq{\textstyle{1\over 2}}$ for density
wave orders, and $\Gamma_{i}$ is certain linear combination of $g_{ij}$supp .
From this we can integrate to obtain the behavior of the corresponding
susceptibility $\chi_{i}$ which scales as
$\chi_{i}\sim(y_{c}-y)^{\gamma_{i}}$. The leading instabilities will be those
with the most negative $\gamma_{i}$, and accordingly the most divergent
$\Gamma_{i}$. We mainly consider the following weak coupling instabilities,

$\displaystyle\Delta_{\text{CDW}_{\pm}}$
$\displaystyle=\braket{\psi_{\alpha+}^{\dagger}\psi_{\beta+}\pm\psi_{\alpha-}^{\dagger}\psi_{\beta-}},~{}\Delta_{\text{SC}}=\braket{f_{\alpha}\psi_{\alpha+}\psi_{\alpha-}},$
(3) $\displaystyle\Delta_{\text{SDW}_{\pm}}$
$\displaystyle=\braket{\psi_{\alpha+}^{\dagger}\bm{\sigma}\psi_{\beta+}\pm\psi_{\alpha-}^{\dagger}\bm{\sigma}\psi_{\beta-}}.$

Note that the form factor $f_{\alpha}=\pm 1$ in $\Delta_{\text{SC}}$
determines the pairing symmetry. For the density wave orders, their real and
imaginary parts flow differently, so depending on whether their imaginary
parts are zero or not, the system could either preserve or break the time
reversal symmetry.

In Fig.2(b) we show the RG results in the presence of SI, i.e. the initial
values of $g_{ij}$ are set by Eq.(2). Interestingly, we find that this map
between lattice interactions and $g_{ij}$ persists for all $y<y_{c}$, and it
is this constant map that enables us to extract the RG flows for $U$ and $V$.
We see the $U$ decays as before, while $V$ increases. The leading instability
in this case is the degenerate complex $\Delta_{\text{CDW}_{-}}$ and imaginary
$\Delta_{\text{SDW}_{\pm}}$. The presence of imaginary parts for these density
wave orders indicate the time reversal symmetry is brokenJiang2021 ;
yu2021evidence ; Mielke2022 ; FENG20211384 ; PhysRevB.104.045122 ;
PhysRevLett.129.167001 ; PhysRevB.104.035142 ,, which is due to the SI effect.
For comparison, we also show in Fig.2(c) the RG analysis for the model but
with no SI effect. In this case, all $g_{ij}$ are set by the largest
interaction $U$, and the resulting weak coupling instability is the $d$-wave
uniform superconductivity, consistent with Ref.Nandkishore2012 ;
PhysRevB.89.144501 .

PDW at intermediate coupling.–Having established the way different
interactions get renormalized in the presence of the SI effect, we now discuss
how the PDW order can emerge at some intermediate energy scale. Based on the
observation discussed above [see Fig.2(b) for example], we consider an
effective model where only a large $V$ is kept. Model similar this has also
been studied recently in Ref.Dong_2023 . We consider a sufficiently large $V$,
and perform a Hartree-Fock mean-field studyPhysRevB.104.075150 ;
PhysRevLett.124.187601 ; PhysRevB.104.214403 ; PhysRevB.102.035136 ;
PhysRevB.102.201104  for the corresponding orders222Justification of the
Hartree-Fock analysis can be achieved by generalizing the model to a large-$N$
theory, where fluctuations around the saddle point are suppressed by $1/N$ and
the Hartree-Fock mean field equations become asymptotically exact. See
doi:10.1146/annurev-conmatphys-031620-102024  and references therein for
examples..  Without loss of generality, we discuss the $AB$ bond only, as the
other bonds follow via $C_{3}$ rotation. On the AB bond, the relevant
interaction is
$2V\cos({\bm{q}}\cdot\bm{\alpha})\psi_{A\sigma}^{\dagger}({\bm{k}})\psi_{B\sigma^{\prime}}^{\dagger}({\bm{k}}^{\prime})\psi_{B\sigma^{\prime}}({\bm{k}}^{\prime}+{\bm{q}})\psi_{A\sigma}({\bm{k}}-{\bm{q}})$
where $\bm{\alpha}=\bm{a}_{1}/2$ and $\bm{a}_{1}$ is the vector connecting two
adjacent $A$ and $B$ sites [see Fig.1 (a)]. Whether this is repulsive or
attractive depends on the momentum transfer ${\bm{q}}$. In the Cooper channel,
we can write it as
$2V\cos[({\bm{k}}+{\bm{k}}^{\prime}-{\bm{q}})\cdot\bm{\alpha}]\psi_{A\sigma}^{\dagger}({\bm{k}})\psi_{B\sigma^{\prime}}^{\dagger}(-{\bm{k}}+{\bm{q}})\psi_{B\sigma^{\prime}}({\bm{k}}^{\prime})\psi_{A\sigma}(-{\bm{k}}^{\prime}+{\bm{q}})$,
which is attractive when
$\cos[({\bm{k}}+{\bm{k}}^{\prime}-{\bm{q}})\cdot\bm{\alpha}]<0$. In the patch
model this condition is met when ${\bm{k}}$ is around $\bm{M}_{A}$,
$-{\bm{k}}+{\bm{q}}$ is around $-\bm{M}_{B}$, and ${\bm{k}}^{\prime}$ is
around $-\bm{M}_{A}$. The interaction then becomes
$-2V\psi_{A}^{\dagger}\psi_{\bar{B}}^{\dagger}\psi_{B}\psi_{\bar{A}}$. When
$V$ is large, we can use

$\psi_{A}^{\dagger}\psi_{\bar{B}}^{\dagger}\psi_{B}\psi_{\bar{A}}\approx\braket{\psi_{A}^{\dagger}\psi_{\bar{B}}^{\dagger}}\psi_{B}\psi_{\bar{A}}+\psi_{A}^{\dagger}\psi_{\bar{B}}^{\dagger}\braket{\psi_{B}\psi_{\bar{A}}}$
(4)

for mean field analysis, and the gap function, defined as
$\Delta_{{\bm{Q}}}\sim\braket{\psi_{B}\psi_{\bar{A}}}$, is apparently a PDW
order with momentum $\bm{M}_{C}$. We note that, due to the presence of SI, the
uniform superconductivity is not a competing order in the large $V$ model,
since such order couples to, e.g. $\psi^{\dagger}_{A}\psi_{A}^{\dagger}$,
which can not be obtained by decomposing Eq.(4).

Similarly, the onsite CDW does not arise since this order parameter couples
to, e.g. $\psi^{\dagger}_{A}\psi_{A}$, and effectively becomes the chemical
potential. However, the bond charge density wave (bond CDW) order can indeed
arise when $V$ becomes strong PhysRevLett.110.126405 , and compete with PDW.
This occurs also when $\cos({\bm{q}}\cdot\bm{\alpha})<0$ and the density
interaction becomes
$-2V\psi_{A}^{\dagger}\psi_{\bar{B}}^{\dagger}\psi_{B}\psi_{\bar{A}}$ with
$V>0$. We can, however, consider

$-\psi_{A}^{\dagger}\psi_{\bar{B}}^{\dagger}\psi_{B}\psi_{\bar{A}}\approx\braket{\psi_{A}^{\dagger}\psi_{B}}\psi_{\bar{B}}^{\dagger}\psi_{\bar{A}}+\psi_{A}^{\dagger}\psi_{B}\braket{\psi_{\bar{B}}^{\dagger}\psi_{\bar{A}}}$
(5)

to absorb the minus sign and manifestly show it contains strong repulsion in
the particle-hole channel. Moreover, the spin indices on the left hand side of
Eq. 5 can be included explicitly, namely
$-\psi_{A\alpha}^{\dagger}\psi_{\bar{B}\beta}^{\dagger}\psi_{B\gamma}\psi_{\bar{A}\delta}\delta_{\alpha\delta}\delta_{\beta\gamma}$,
and using the $SU(2)$ identity
$\delta_{\alpha\delta}\delta_{\beta\gamma}=(\sigma_{\alpha\gamma}\delta_{\beta\delta}+\bm{\sigma}_{\alpha\gamma}\bm{\sigma}_{\beta\delta})/2$,
it is easy to see that the order parameter
$\Delta_{\bm{Q}}=\braket{\psi_{A}^{\dagger}\psi_{B}}$ represents a bond
CDW/SDW order which are directly related to those density waves that arise in
the weak coupling limit. In fact, the degeneracy is an artifact of the
oversimplification of the patch model. Considering that the SDW order has to
break a global $SU(2)$ symmetry, and upon taking the whole Fermi surface into
account, we assume for the time being that the CDW order is more likely to
occur than SDW. In the following we thus constrain ourselves to CDW order
competing with PDW at large $V$.

To inspect the competition between PDW and CDW at large $V>0$, we can
calculate the susceptibilities for both of the two orders in random-phase-
approximation (RPA). The result is
$\chi_{PDW}=\Pi_{pp}^{0}({\bm{Q}})/[1-2V\Pi_{pp}^{0}({\bm{Q}})]$ and
$\chi_{CDW}=\Pi_{ph}^{0}/[1-2V\Pi_{ph}^{0}({\bm{Q}})]$. Here the bare
susceptibilities $\Pi^{0}_{pp}>0$ and $\Pi_{ph}>0$ are obtained using free
fermion propagators. The ordering condition is given when the denominators
vanish. The problem is now reduced to comparing the strength of the bare
$\Pi_{pp}^{(0)}$ and $\Pi_{ph}^{(0)}$. Within the patch model, the particle-
particle and particle-hole bubbles can be obtained in a straightforward way,
and the results are333The apparent dependence on the artificial energy scale
$\Lambda$ is unphysical, and there should be contributions from outside the
patches which eventually add up to make the final results larger and
independent of $\Lambda$.

$\displaystyle\Pi_{pp}^{(0)}({\bm{Q}})$ $\displaystyle={\textstyle{1\over
4\sqrt{3}\pi^{2}t}}\ln{\textstyle{\Lambda\over\max\\{T,|\mu|\\}}}$ (6)
$\displaystyle\Pi_{ph}^{(0)}({\bm{Q}})$ $\displaystyle={\textstyle{1\over
8\sqrt{3}\pi^{2}t}}\ln{\textstyle{\Lambda\over\max\\{T,|\mu|\\}}}\ln{\textstyle{\Lambda\over\max\\{T,\mu,|t^{\prime}|\\}}},$

where $\mu=0$ at the van Hove filling, and $t^{\prime}$ is the next nearest
neighbor hopping which can spoil the Fermi surface nesting in the particle-
hole channel. In Fig. 3(a) we compare both $\Pi_{pp}^{(0)}$ and
$\Pi_{ph}^{(0)}$ as a function of temperature $T$ at $\mu-0$. We see that the
$\ln^{2}$ divergence of $\Pi_{ph}^{(0)}$ becomes dominant at small $T$ for the
case of perfect nesting ($t^{\prime}=0$), while at larger $T$,
$\Pi_{pp}^{(0)}$ is larger even with perfect nesting. The consequence is that,
in the weak coupling limit, one has to go to small $T$ to see the instability,
where CDW wins over PDW. For large $V$, however, a relatively smaller
$\Pi_{pp}^{(0)}$ or $\Pi_{ph}^{(0)}$ is enough to induce the instability,
which could be in the regime where $\Pi_{pp}^{(0)}>\Pi_{ph}^{(0)}$ and PDW is
the leading instability. Note that Eq.(6) is only a crude estimation about the
relative strength between $\Pi_{pp}^{(0)}({\bm{Q}})$ and
$\Pi_{ph}^{(0)}({\bm{Q}})$. Away from the van Hove filling ($\mu\neq 0$), we
expect that both bubbles are reduced, and therefore one needs to reach smaller
$T$ in order to reach the instability. This results in the schematic picture
we show as the inset of Fig.2(a).

Figure 3: (a)Comparison between $\Pi_{pp}^{(0)}({\bm{Q}})$ and
$\Pi_{ph}^{(0)}({\bm{Q}})$ as a function of temperature $T$. The fact that
$\Pi_{pp}^{(0)}({\bm{Q}})>\Pi_{ph}^{(0)}({\bm{Q}})$ at larger $T$ indicates
that PDW order may win over CDW order when interaction is sufficiently strong.
(b) The weak coupling instabilities (at perfect nesting) when other
interactions such as $J_{0}$ and $K_{0}$ are taken into account. (c) RG flows
for the lattice interactions for $U(0)=0.8t$ and $V(0)=0.3t$ and for various
$J_{0}$ and $K_{0}$. In all cases there is a wide range in energy scale where
$V$ is the largest.

Other interactions. We note that other than $U$ and $V$, further lattice
interaction might be worth analyzing. For instance, we can consider the bond
singlet pair hopping term, for which, again taking $AB$ bond as an example,
the on-bond contribution is
$J_{0}\sum_{i}P^{\dagger}_{AB}(\bm{r}_{i})P_{AB}(\bm{r}_{i})$, where the pair
operator at position $\bm{r}_{i}$ is given by
$P^{\dagger}_{AB}(\bm{r}_{i})=c_{A\uparrow}^{\dagger}(\bm{r}_{i}+\bm{\alpha})c^{\dagger}_{B\downarrow}(\bm{r}_{i})-c_{A\downarrow}^{\dagger}(\bm{r}_{i}+\bm{\alpha})c^{\dagger}_{B\uparrow}(\bm{r}_{i})$
and $\bm{\alpha}=\bm{a}_{1}/2$ is the vector connecting two adjacent $A$ and
$B$ sites. In fact, the $J_{0}$ term just reads as the exchange interaction.
The site pair hopping term such as $K_{0}\sum_{\langle
ij\rangle}P^{\dagger}_{i}P_{j}$ with
$P_{i}^{\dagger}=c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}$ is also
possible. With the SI effect, it is easy to realize that $J_{0}$ contributes
to both $g_{1j}$ and $g_{2j}$. In particular,
$g_{11}(0)=g_{14}(0)=-g_{12}(0)=-g_{13}(0)=2J_{0}$, and
$-g_{21}(0)=-g_{23}(0)=g_{22}(0)=g_{24}(0)=2V+2J_{0}$. The $K_{0}$ term only
contributes to $g_{3j}$: $g_{31}(0)=g_{32}(0)=-g_{33}(0)=-g_{34}(0)=2K_{0}$.
We expect that $J_{0}$ and $K_{0}$ are of the same order, and both are orders
of magnitudes smaller than $U$ and $V$doi:10.1098/rspa.1963.0204 ;
PhysRevLett.58.1899 ; PhysRevLett.60.71 ; PhysRevLett.60.72 . In the presence
of $J_{0}$ and $K_{0}$, both of them increase under RG, leading to different
orders in the weak coupling limit which we present in Fig.3(b). These include
$s$\- and $f$-wave uniform SC, the real CDW order with odd parity (CDW-).
Interestingly, if we keep $K_{0}=0$ there is still a constant map between
$U,V,J_{0}$ and $g_{ij}$ for all $y<y_{c}$, and the weak coupling
instabilities are degenerate imaginary $\Delta_{\text{SDW}_{\pm}}$ for
$J_{0}>0$, and a complex $\Delta_{\text{CDW}_{-}}$ for $J_{0}<0$ [see the blue
and purple line in Fig.3(b)]. A finite $K_{0}$ spoils the constant map between
lattice interactions and $g_{ij}$ at $y$ near $y_{c}$, but at some
intermediate $y^{*}<y_{c}$ we can still approximately determine the flow of
$U,V,J_{0}$ and $K_{0}$. Fig.3(c) (see also in supp ) shows that $V$ can still
be the largest at intermediate energy scale for various cases, which justifies
the applicability of our effective $V$ model discussed above. We close by
noting the possibility that when a sizable $K_{0}$ is present initially,
$K_{0}$ could become comparable to $V$ even at $y<y^{*}$. Since $K_{0}$ is
related to the formation of uniform SC, one needs to study the competition
between uniform SC and PDW in this case.

### Discussion.-

We have focused our analysis entirely on the case of electrons subject to
instantaneous repulsive interactions on the Kagome lattice near p-type van
Hove filling. At the current level of understanding, it is not yet settled to
which extent the kagome metals AV3Sb5 (A=K,Cs,Rb) are faithfully represented
by this simplified model. Previous attempts to improve on the microscopic
rigor by taking into account the multi-orbital nature at the Fermi level
PhysRevLett.127.177001 ; PhysRevLett.127.217601 , combined with the relevance
of phononic contributions PhysRevLett.127.046401 , suggest the need for
further analysis to quantitatively approach the experimental setup in AV3Sb5,
and further experiments involving Kagome metals are necessary to validate the
nature of the CDW observed at scales as high as $T\sim 90K$. In particular,
such order would of course affect the reconstructed bands out of which the
superconductivity ultimately develops. A coherent theory of such orders,
including the effects of phonons, will only be feasible once the precise
nature of the CDW order in this system is well-characterized.

Our analysis does, however, make a rare microscopically founded statement
about the realization of PDW in a physically sensible model, specifically in a
kagome metal nearby a p-type vHS with local and nearest neighbor Coulomb
repulsion. This is because SI crucially affects the renormalization of
repulsive interactions in Kagome metals, and thus the extended repulsive
forces are enhanced relative to the onsite Hubbard repulsion. As a result,
there is an increased tendency towards PDW and bond CDW order, which we
identified without any need for a potentially biased mean field analysis.
Further input from experiment is needed at this stage to bridge the ideas
described in this paper with the thus far concluced phase diagram of Kagome
metals.

### Acknowledgement.

We sincerely thank Mengxing Ye for useful discussions, and J. Beyer, M.
Dürrnagel, T. Müller, J. Potten, and T. Schwemmer for ongoing collaborations
on related topics. Y.-M. W. acknowledges the Gordon and Betty Moore
Foundation’s EPiQS Initiative through GBMF8686 for support at Stanford
University.R.T. acknowledges support from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) through QUAST FOR 5249-449872909 (Project
P3), through Project-ID 258499086-SFB 1170, and from the Würzburg-Dresden
Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat
Project-ID 390858490-EXC 2147. S.R. is supported by the Department of Energy,
Office of Basic Energy Sciences, Division of Materials Sciences and
Engineering, under Contract No. DE-AC02-76SF00515.

## References

  * (1) T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Nature 492, 406 (2012).
  * (2) L. Ye, M. Kang, J. Liu, F. von Cube, C. R. Wicker, T. Suzuki, C. Jozwiak, A. Bostwick, E. Rotenberg, D. C. Bell, L. Fu, R. Comin, and J. G. Checkelsky, Nature 555, 638 (2018).
  * (3) M. Kang et al., Nature materials 19, 163 (2020).
  * (4) J.-X. Yin et al., Nature 583, 533 (2020).
  * (5) B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. H. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, E. Ertekin, T. M. McQueen, and E. S. Toberer, Phys. Rev. Materials 3, 094407 (2019).
  * (6) Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, G. Xu, Z. Guguchia, J. He, M. S. Hossain, X. Liu, J. Ruff, L. Kautzsch, S. S. Zhang, G. Chang, I. Belopolski, Q. Zhang, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Wang, R. Thomale, T. Neupert, S. D. Wilson, and M. Z. Hasan, Nature Materials 20, 1353 (2021).
  * (7) T. Neupert, M. M. Denner, J.-X. Yin, R. Thomale, and M. Z. Hasan, Nature Physics 18, 137 (2022).
  * (8) H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, S. Ma, S. Ni, H. Zhang, Q. Yin, C. Gong, Z. Tu, H. Lei, H. Tan, S. Zhou, C. Shen, X. Dong, B. Yan, Z. Wang, and H.-J. Gao, Nature 599, 222 (2021).
  * (9) D. F. Agterberg, J. S. Davis, S. D. Edkins, E. Fradkin, D. J. Van Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky, J. M. Tranquada, and Y. Wang, Annual Review of Condensed Matter Physics 11, 231 (2020).
  * (10) E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Rev. Mod. Phys. 87, 457 (2015).
  * (11) A. Himeda, T. Kato, and M. Ogata, Phys. Rev. Lett. 88, 117001 (2002).
  * (12) K.-Y. Yang, W. Q. Chen, T. M. Rice, M. Sigrist, and F.-C. Zhang, New Journal of Physics 11, 055053 (2009).
  * (13) M. Raczkowski, M. Capello, D. Poilblanc, R. Frésard, and A. M. Oleś, Phys. Rev. B 76, 140505 (2007).
  * (14) E. Berg, E. Fradkin, E.-A. Kim, S. A. Kivelson, V. Oganesyan, J. M. Tranquada, and S. C. Zhang, Phys. Rev. Lett. 99, 127003 (2007).
  * (15) M. Capello, M. Raczkowski, and D. Poilblanc, Phys. Rev. B 77, 224502 (2008).
  * (16) P. A. Lee, Phys. Rev. X 4, 031017 (2014).
  * (17) Y. Wang, D. F. Agterberg, and A. Chubukov, Phys. Rev. Lett. 114, 197001 (2015).
  * (18) Y. Wang, D. F. Agterberg, and A. Chubukov, Phys. Rev. B 91, 115103 (2015).
  * (19) S. D. Edkins, A. Kostin, K. Fujita, A. P. Mackenzie, H. Eisaki, S. Uchida, S. Sachdev, M. J. Lawler, E.-A. Kim, J. C. S. Davis, and M. H. Hamidian, Science 364, 976 (2019).
  * (20) Y. Wang, S. D. Edkins, M. H. Hamidian, J. C. S. Davis, E. Fradkin, and S. A. Kivelson, Phys. Rev. B 97, 174510 (2018).
  * (21) Z. Wu, Y.-M. Wu, and F. Wu, Phys. Rev. B 107, 045122 (2023).
  * (22) Y.-M. Wu, Z. Wu, and H. Yao, Pair-density-wave and chiral superconductivity in twisted bilayer transition-metal-dichalcogenides. ARXIV.2203.05480, 2022.
  * (23) D. Shaffer and L. H. Santos, Triplet Pair-Density Wave Superconductivity on the $\pi$-flux Square Lattice. ARXIV.2210.16324, 2022.
  * (24) D. Shaffer, F. J. Burnell, and R. M. Fernandes, Weak-Coupling Theory of Pair Density-Wave Instabilities in Transition Metal Dichalcogenides. ARXIV.2209.14469, 2022.
  * (25) J. Ge, P. Wang, Y. Xing, Q. Yin, H. Lei, Z. Wang, and J. Wang, (2022).
  * (26) Y.-M. Wu, P. A. Nosov, A. A. Patel, and S. Raghu, Phys. Rev. Lett. 130, 026001 (2023).
  * (27) M. L. Kiesel and R. Thomale, Phys. Rev. B 86, 121105 (2012).
  * (28) X. Wu, T. Schwemmer, T. Müller, A. Consiglio, G. Sangiovanni, D. Di Sante, Y. Iqbal, W. Hanke, A. P. Schnyder, M. M. Denner, M. H. Fischer, T. Neupert, and R. Thomale, Phys. Rev. Lett. 127, 177001 (2021).
  * (29) M. Kang, S. Fang, J.-K. Kim, B. R. Ortiz, S. H. Ryu, J. Kim, J. Yoo, G. Sangiovanni, D. Di Sante, B.-G. Park, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, S. D. Wilson, J.-H. Park, and R. Comin, Nature Physics (2022).
  * (30) Y. Hu, X. Wu, B. R. Ortiz, S. Ju, X. Han, J. Ma, N. C. Plumb, M. Radovic, R. Thomale, S. D. Wilson, A. P. Schnyder, and M. Shi, Nature Communications 13, 2220 (2022).
  * (31) I. Dzyaloshinskii, Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki 93, 1487 (1987).
  * (32) H. J. Schulz, Europhysics Letters 4, 609 (1987).
  * (33) H. D. Scammell, J. Ingham, T. Li, and O. P. Sushkov, Nature Communications 14, 605 (2023).
  * (34) N. Furukawa, T. M. Rice, and M. Salmhofer, Phys. Rev. Lett. 81, 3195 (1998).
  * (35) R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Nature Physics 8, 158 (2012).
  * (36) R. Nandkishore, R. Thomale, and A. V. Chubukov, Phys. Rev. B 89, 144501 (2014).
  * (37) H. Isobe, N. F. Q. Yuan, and L. Fu, Phys. Rev. X 8, 041041 (2018).
  * (38) Y.-P. Lin and R. M. Nandkishore, Phys. Rev. B 100, 085136 (2019).
  * (39) Supplemental Information.
  * (40) For comparison, one can also consider the alternative kagome van Hove filling which is of m-type [see Fig.1(b)]. There, each $\psi_{\alpha}$ is a superposition of two distinct lattice fermions and it can be directly verified that the Hubbard $U$ then has equal contributions to all $g_{i}$.
  * (41) Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, G. Xu, Z. Guguchia, J. He, M. S. Hossain, X. Liu, J. Ruff, L. Kautzsch, S. S. Zhang, G. Chang, I. Belopolski, Q. Zhang, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Wang, R. Thomale, T. Neupert, S. D. Wilson, and M. Z. Hasan, Nature Materials 20, 1353 (2021).
  * (42) L. Yu, C. Wang, Y. Zhang, M. Sander, S. Ni, Z. Lu, S. Ma, Z. Wang, Z. Zhao, H. Chen, K. Jiang, Y. Zhang, H. Yang, F. Zhou, X. Dong, S. L. Johnson, M. J. Graf, J. Hu, H.-J. Gao, and Z. Zhao, Evidence of a hidden flux phase in the topological kagome metal CsV3Sb5, 2021.
  * (43) C. Mielke, D. Das, J.-X. Yin, H. Liu, R. Gupta, Y.-X. Jiang, M. Medarde, X. Wu, H. C. Lei, J. Chang, P. Dai, Q. Si, H. Miao, R. Thomale, T. Neupert, Y. Shi, R. Khasanov, M. Z. Hasan, H. Luetkens, and Z. Guguchia, Nature 602, 245 (2022).
  * (44) X. Feng, K. Jiang, Z. Wang, and J. Hu, Science Bulletin 66, 1384 (2021).
  * (45) Y.-P. Lin and R. M. Nandkishore, Phys. Rev. B 104, 045122 (2021).
  * (46) J.-T. Jin, K. Jiang, H. Yao, and Y. Zhou, Phys. Rev. Lett. 129, 167001 (2022).
  * (47) T. Park, M. Ye, and L. Balents, Phys. Rev. B 104, 035142 (2021).
  * (48) J.-W. Dong, Z. Wang, and S. Zhou, Physical Review B 107, (2023).
  * (49) J. Zang, J. Wang, J. Cano, and A. J. Millis, Phys. Rev. B 104, 075150 (2021).
  * (50) C. Repellin, Z. Dong, Y.-H. Zhang, and T. Senthil, Phys. Rev. Lett. 124, 187601 (2020).
  * (51) N. C. Hu and A. H. MacDonald, Phys. Rev. B 104, 214403 (2021).
  * (52) Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Phys. Rev. B 102, 035136 (2020).
  * (53) H. Pan, F. Wu, and S. Das Sarma, Phys. Rev. B 102, 201104 (2020).
  * (54) Justification of the Hartree-Fock analysis can be achieved by generalizing the model to a large-$N$ theory, where fluctuations around the saddle point are suppressed by $1/N$ and the Hartree-Fock mean field equations become asymptotically exact. See doi:10.1146/annurev-conmatphys-031620-102024  and references therein for examples.
  * (55) M. L. Kiesel, C. Platt, and R. Thomale, Phys. Rev. Lett. 110, 126405 (2013).
  * (56) The apparent dependence on the artificial energy scale $\Lambda$ is unphysical, and there should be contributions from outside the patches which eventually add up to make the final results larger and independent of $\Lambda$.
  * (57) J. Hubbard and B. H. Flowers, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 276, 238 (1963).
  * (58) S. Kivelson, W.-P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 58, 1899 (1987).
  * (59) J. T. Gammel and D. K. Campbell, Phys. Rev. Lett. 60, 71 (1988).
  * (60) S. Kivelson, W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 60, 72 (1988).
  * (61) M. M. Denner, R. Thomale, and T. Neupert, Phys. Rev. Lett. 127, 217601 (2021).
  * (62) H. Tan, Y. Liu, Z. Wang, and B. Yan, Phys. Rev. Lett. 127, 046401 (2021).
  * (63) D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, Annual Review of Condensed Matter Physics 13, 239 (2022).

ONLINE SUPPORTING MATERIAL
Sublattice Interference promotes Pair Density Wave in Kagome Metals

Yi-Ming Wu,1, Ronny Thomale2 and S. Raghu1

1 Stanford Institute for Theoretical Physics, Stanford University, Stanford,
California 94305, USA

2 Institute for Theoretical Physics and Astrophysics, University of Würzburg,
D-97074 Würzburg, Germany

In this Supplemental Material we show i) the sublattice interference is robust
against the inclusion of long range hoppings, ii) how the sublattice
interference sets the initial values for different patch interactions and iii)
the detailed calculation for the bare PDW and CDW susceptibilities.

##  I Robustness of Sublattice inteference

The tight binding model for the Kagome lattice including the nearest neighbor
(NN) hopping $t$, the next nearest neighbor (NNN) hopping $t^{\prime}$ and the
next-next nearest neighbor (NNNN) hopping $t^{\prime\prime}$ is given by

$H_{0}=\sum_{{\bm{k}}}(c_{A}^{\dagger}({\bm{k}}),c_{B}^{\dagger}({\bm{k}}),c_{C}^{\dagger}({\bm{k}}))\mathcal{H}({\bm{k}})\begin{pmatrix}c_{A}({\bm{k}})\\\
c_{B}({\bm{k}})\\\ c_{C}({\bm{k}})\\\ \end{pmatrix}$ (7)

where

$\mathcal{H}({\bm{k}})=\begin{pmatrix}2t^{\prime\prime}\sum_{i=1}^{3}\cos({\bm{k}}\cdot\bm{a}_{i})&2t\cos({\bm{k}}\cdot\bm{a}_{1}/2)+2t^{\prime}\cos({\bm{k}}\cdot\bm{A}_{1})&2t\cos({\bm{k}}\cdot\bm{a}_{3}/2)+2t^{\prime}\cos({\bm{k}}\cdot\bm{A}_{3})\\\
2t\cos({\bm{k}}\cdot\bm{a}_{1}/2)+2t^{\prime}\cos({\bm{k}}\cdot\bm{A}_{1})&2t^{\prime\prime}\sum_{i=1}^{3}\cos({\bm{k}}\cdot\bm{a}_{i})&2t\cos({\bm{k}}\cdot\bm{a}_{2}/2)+2t^{\prime}\cos({\bm{k}}\cdot\bm{A}_{2})\\\
2t\cos({\bm{k}}\cdot\bm{a}_{3}/2)+2t^{\prime}\cos({\bm{k}}\cdot\bm{A}_{3})&2t\cos({\bm{k}}\cdot\bm{a}_{2}/2)+2t^{\prime}\cos({\bm{k}}\cdot\bm{A}_{2})&2t^{\prime\prime}\sum_{i=1}^{3}\cos({\bm{k}}\cdot\bm{a}_{i})\\\
\end{pmatrix}$ (8)

and

$\displaystyle\bm{a}_{1}$ $\displaystyle={\textstyle{1\over
2}}(1,\sqrt{3}),~{}\bm{a}_{2}=(1,0),~{}\bm{a}_{3}={\textstyle{1\over
2}}(-1,\sqrt{3})$ (9) $\displaystyle\bm{A}_{1}$
$\displaystyle={\textstyle{1\over
4}}(-3,\sqrt{3}),~{}\bm{A}_{2}={\textstyle{1\over
2}}(0,\sqrt{3}),~{}\bm{A}_{3}={\textstyle{1\over 4}}(3,\sqrt{3})$

The Hamiltonian can be diagonalized as
$H_{0}=\sum_{n,{\bm{k}}}\psi_{n}^{\dagger}({\bm{k}})\psi_{n}({\bm{k}})$ by
using the following linear transformation between the band fermions and the
lattice fermions:

$c_{s}(k)=\sum_{n}u_{sn}(k)\psi_{n}(k)$ (10)

where $s\in(A,B,C)$ and $n\in(1,2,3)$ is the band index. If we focus on the
middle band of the kagome lattice where the p-type VHS locates, we can omit
the band index $n$. The weight $u_{s}({\bm{k}})$ shows strong sublattice
interference, i.e. the distribution in momentum space strongly depends on the
sublattice index $s$. For clarity, in Fig.4 we show the configuration of
$|u_{A}({\bm{k}})|^{2}$ for different situations with and without including
longer range hopping amplitudes. The results show that even in the presence of
$t^{\prime}$ and $t^{\prime\prime}$, $|u_{A}({\bm{k}})|^{2}$ changes little.

If we consider only three patches in the small vicinity of these VHS points
($\bm{M}_{A},\bm{M}_{B},\bm{M}_{C}$), we find that $u_{s}({\bm{k}})$ satisfies

$\displaystyle\text{if}~{}{\bm{k}}\in\bm{M}_{A},~{}u_{A}({\bm{k}})=1,~{}u_{B}({\bm{k}})=u_{C}({\bm{k}})=0,$
(11)
$\displaystyle\text{if}~{}{\bm{k}}\in\bm{M}_{B},~{}u_{B}({\bm{k}})=1,~{}u_{A}({\bm{k}})=u_{C}({\bm{k}})=0,$
$\displaystyle\text{if}~{}{\bm{k}}\in\bm{M}_{C},~{}u_{C}({\bm{k}})=1,~{}u_{A}({\bm{k}})=u_{B}({\bm{k}})=0.$

These results are immune to the presence of a nonzero $t^{\prime}$ or
$t^{\prime\prime}$, thus the sublattice interference is a robust phenomenon.

Figure 4: Momentum space configuration of $|u_{A}({\bm{k}})|^{2}$ for
different cases. Here we set $t=1$.

##  II Map from lattice interactions onto the patch model

For the three patch model, there are only four symmetry allowed interactions
$g_{1},g_{2},g_{3}$ and $g_{4}$, which are defined through the most general
interacting Hamiltonian

$H_{I}=g_{1}\psi^{\dagger}_{a}\psi^{\dagger}_{b}\psi_{a}\psi_{b}+g_{2}\psi^{\dagger}_{a}\psi^{\dagger}_{b}\psi_{b}\psi_{a}+g_{3}\psi^{\dagger}_{a}\psi^{\dagger}_{a}\psi_{b}\psi_{b}+g_{4}\psi^{\dagger}_{a}\psi^{\dagger}_{a}\psi_{a}\psi_{a}$
(12)

where $a,b\in(M_{1},M_{2},M_{3})$ are the patch indices. When generalized to
six patch model, we have

$H_{I}=\sum_{i,j=1}^{4}\sum_{\left\\{\alpha_{j},\tau_{j}\right\\}}g_{ij}\psi^{\dagger}_{\alpha_{1}\tau_{1}}\psi^{\dagger}_{\alpha_{2}\tau_{2}}\psi_{\alpha_{3}\tau_{3}}\psi_{\alpha_{4}\tau_{4}},$
(13)

All the interactions can be projected onto the three patch interactions
$g_{i}$ using Eq.(11). For the onsite Hubbard $U$, we have

$H_{U}=U\sum_{s}\sum_{k,k^{\prime},q}u^{*}_{s}(k)u^{*}_{s}(k^{\prime})u_{s}(k^{\prime}-q)u_{s}(k+q)\psi^{\dagger}_{\uparrow}\psi^{\dagger}_{\downarrow}\psi_{\downarrow}\psi_{\uparrow}$
(14)

where the four fermion term $\psi^{\dagger}\psi^{\dagger}\psi\psi$ has the
same momentum configurations as the $u^{*}u^{*}uu$ factors. Because all these
$u_{s}$ factors have the same sublattice index, and according to (11), all the
four fermion operators must be in the same patch, otherwise the term vanishes.
Therefore we have

$g_{4j}(0)=U$ (15)

where $g_{4j}(0)$ means the bare interaction, i.e. before RG.

For $H_{V}$ term, its corresponding form in momentum space is

$H_{V}=2V\sum_{k,k^{\prime},q}\cos(\bm{q}\cdot\bm{\alpha}_{1})u_{A}^{*}(k)u_{B}^{*}(k^{\prime})u_{B}(k^{\prime}+q)u_{A}(k-q)\psi^{\dagger}_{\sigma}\psi^{\dagger}_{\sigma^{\prime}}\psi_{\sigma^{\prime}}\psi_{\sigma}+...$
(16)

where $...$ stands for similar contributions but with $AC$ and $BC$
combinations, and $\bm{\alpha}_{1}=(1,\sqrt{3})/4$ connects two adjacent $A$
and $B$ sites. Because of the SI condition Eq.(11), both ${\bm{k}}$ and
${\bm{k}}-\bm{q}$ should be on the $\bm{M}_{1}$ patch, so $\bm{q}=0$ and
$\cos(\bm{q}\cdot\bm{\alpha}_{1})=1$. By choosing different patches, we can
also have $\cos(\bm{q}\cdot\bm{\alpha}_{1})=-1$. This interaction apparently
contributes to $g_{2j}$, thus

$g_{21}(0)=g_{23}(0)=-2V,~{}~{}g_{22}(0)=g_{24}(0)=2V$ (17)

If we include the site-to-site pair hopping term
$H_{C}=K_{0}\sum_{\braket{ij}}P_{i}^{\dagger}P_{j}$ interaction, it is
completely parallel to get something like

$K_{0}\sum_{k,k^{\prime},q}2\cos(\bm{q}\cdot\bm{\alpha}_{1})u_{A}^{*}(k)u_{A}^{*}(-k+q)u_{B}(k^{\prime})u_{B}(-k^{\prime}+q)\psi^{\dagger}_{\uparrow}\psi^{\dagger}_{\downarrow}\psi_{\downarrow}\psi_{\uparrow}+...$
(18)

Again we use the sublattice interference condition in (11) to conclude that
for this kind of interaction we must have $\bm{q}=2\bm{M}_{1}$ and it
contributes to $g_{3j}$. For the coefficients $u_{s}(\bm{q})$ and the fermions
$\psi(\bm{q})$, we must have the condition $2\bm{M}_{1}=0$. However, this is
not the case for the function $\cos(q\cdot\alpha_{1})$, because
$\bm{\alpha}_{1}$ is half of the unit cell primitive vector, and we have
$\cos(q\cdot\alpha_{1})=-1$ instead. Similarly, we can have situations when
$\cos(q\cdot\alpha_{1})=1$, depending on the combination of patches. The
resulting projections are

$g_{31}(0)=g_{32}(0)=2K_{0},~{}~{}g_{33}(0)=g_{34}(0)=-2K_{0}$ (19)

The bond pair hopping is a bit more complex. We only show the $AB$ bond for
simplicity, but other bonds are easy to obtain. First of all, the momentum
space representation of the pair operators are

$\displaystyle
P_{AB}^{\dagger}(\bm{q})=\sum_{{\bm{k}}}\left(c_{A\uparrow}^{\dagger}({\bm{k}})c_{B\downarrow}^{\dagger}(-{\bm{k}}+\bm{q})-c_{A\downarrow}^{\dagger}({\bm{k}})c_{B\uparrow}^{\dagger}(-{\bm{k}}+\bm{q})\right)e^{i(\bm{q}/2-{\bm{k}})\cdot\bm{\alpha}_{1}}$
(20) $\displaystyle
P_{AB}(\bm{q})=\sum_{{\bm{k}}}\left(c_{B\downarrow}(-{\bm{k}}+\bm{q})c_{A\uparrow}({\bm{k}})-c_{B\uparrow}(-{\bm{k}}+\bm{q})c_{A\downarrow}({\bm{k}})\right)e^{-i(\bm{q}/2-{\bm{k}})\cdot\bm{\alpha}_{1}}$
$\displaystyle
P_{BA}^{\dagger}(\bm{q})=\sum_{{\bm{k}}}\left(c_{B\uparrow}^{\dagger}({\bm{k}})c_{A\downarrow}^{\dagger}(-{\bm{k}}+\bm{q})-c_{B\downarrow}^{\dagger}({\bm{k}})c_{A\uparrow}^{\dagger}(-{\bm{k}}+\bm{q})\right)e^{i(\bm{q}/2-{\bm{k}})\cdot\bm{\alpha}_{1}}$
$\displaystyle
P_{BA}(\bm{q})=\sum_{{\bm{k}}}\left(c_{A\downarrow}(-{\bm{k}}+\bm{q})c_{B\uparrow}({\bm{k}})-c_{A\uparrow}(-{\bm{k}}+\bm{q})c_{B\downarrow}({\bm{k}})\right)e^{-i(\bm{q}/2-{\bm{k}})\cdot\bm{\alpha}_{1}}$

Using this, we obtain for the on-bond pairing hopping,

$\displaystyle
J_{0}\sum_{i}P_{AB}^{\dagger}(\bm{r}_{i})P_{AB}(\bm{r}_{i})+(A\leftrightarrow
B)=J_{0}\sum_{q}P_{AB}^{\dagger}(\bm{q})P_{AB}(\bm{q})+P_{BA}^{\dagger}(\bm{q})P_{BA}(\bm{q})$
(21)
$\displaystyle=2J_{0}\sum_{k,k^{\prime},q}\sum_{\sigma\neq\sigma^{\prime}}\cos[({\bm{k}}-{\bm{k}}^{\prime})\cdot\bm{\alpha}_{1}]u_{A}^{*}(k)u_{B}^{*}(-k+q)u_{B}(-k^{\prime}+q)u_{A}(k^{\prime})\psi^{\dagger}_{\sigma}\psi^{\dagger}_{\sigma^{\prime}}\psi_{\sigma^{\prime}}\psi_{\sigma}$
$\displaystyle+2J_{0}\sum_{k,k^{\prime},q}\sum_{\sigma\neq\sigma^{\prime}}\cos[({\bm{k}}-{\bm{k}}^{\prime})\cdot\bm{\alpha}_{1}]u_{A}^{*}(k)u_{B}^{*}(-k+q)u_{A}(k^{\prime})u_{B}(-k^{\prime}+q)\psi^{\dagger}_{\sigma}\psi^{\dagger}_{\sigma^{\prime}}\psi_{\sigma^{\prime}}\psi_{\sigma}$

In the above expression, since both ${\bm{k}}$ and ${\bm{k}}^{\prime}$ are in
the $\bm{M}_{1}$ patch, ${\bm{k}}-{\bm{k}}^{\prime}=0$ and the $\cos$ term
equals $1$. It is the apparent that the first part of the result contributes
to $g_{2j}$ while the second parts contributes to $g_{1j}$, with both positive
$2J_{0}$. The same analysis can be applied to the nearest neighbor bond pair
hopping,

$\displaystyle
J_{1}\sum_{i}P_{AB}^{\dagger}(\bm{r}_{i})P_{BA}(\bm{r}_{i}+\bm{\alpha})+P_{AB}^{\dagger}(\bm{r}_{i})P_{BA}(\bm{r}_{i}-\bm{\alpha})+(A\leftrightarrow
B)$ (22)
$\displaystyle=2J_{1}\sum_{q}\left(\cos\bm{q}\cdot\bm{\alpha}_{1}P_{AB}^{\dagger}(\bm{q})P_{BA}(\bm{q})+A\leftrightarrow
B\right)$
$\displaystyle=4J_{1}\sum_{k,k^{\prime},q}\sum_{\sigma\neq\sigma^{\prime}}\cos\bm{q}\cdot\bm{\alpha}_{1}\cos[({\bm{k}}-{\bm{k}}^{\prime})\cdot\bm{\alpha}_{1}]u^{*}_{A}(k)u^{*}_{B}(-k+q)u_{B}(k^{\prime})u_{A}(-k^{\prime}+q)\psi^{\dagger}_{\sigma}\psi^{\dagger}_{\sigma^{\prime}}\psi_{\sigma^{\prime}}\psi_{\sigma}$
$\displaystyle+4J_{1}\sum_{k,k^{\prime},q}\sum_{\sigma\neq\sigma^{\prime}}\cos\bm{q}\cdot\bm{\alpha}_{1}\cos[({\bm{k}}-{\bm{k}}^{\prime})\cdot\bm{\alpha}_{1}]u^{*}_{A}(k)u^{*}_{B}(-k+q)u_{A}(-k^{\prime}+q)u_{B}(k^{\prime})\psi^{\dagger}_{\sigma}\psi^{\dagger}_{\sigma^{\prime}}\psi_{\sigma^{\prime}}\psi_{\sigma}$

Because of the SI condition, ${\bm{k}}$ is at $\bm{M}_{1}$ while
${\bm{k}}^{\prime}$ is at $\bm{M}_{2}$, then
${\bm{k}}-{\bm{k}}^{\prime}=\bm{M}_{3}$. It then follows that
$\cos\bm{q}\cdot\bm{\alpha}_{1}\cos[({\bm{k}}-{\bm{k}}^{\prime})\cdot\bm{\alpha}_{1}]=-1$.
Similar to the $J_{0}$ term, $J_{1}$ term also has two distinct parts, of
which the first part contributes to $g_{2j}$ and the second part contributes
to $g_{1j}$.

##  III Renormalization group analysis

Having projected all the lattice interactions onto $g_{ij}$, we can
investigate how they affected one another when high energy degrees of freedom
is integrated out. This is captured by the following one-loop RG
equationsIsobe2018 ,

$\displaystyle\dot{g}_{11}$
$\displaystyle=2d_{1}(g_{11}g_{22}+g_{31}g_{32}-g_{13}^{2}-g_{33}^{2}),$ (23)
$\displaystyle\dot{g}_{12}$
$\displaystyle=2d_{1}(g_{12}g_{24}+g_{32}g_{34}-g_{12}g_{14}-g_{32}g_{34}),$
$\displaystyle\dot{g}_{13}$
$\displaystyle=2d_{1}(g_{13}g_{22}+g_{32}g_{33}-g_{11}g_{13}-g_{31}g_{33}),$
$\displaystyle\dot{g}_{14}$ $\displaystyle=2d_{1}(g_{14}g_{24}-g_{14}^{2}),$
$\displaystyle\dot{g}_{21}$ $\displaystyle=d_{1}(g_{21}g_{24}+g_{31}g_{34}),$
$\displaystyle\dot{g}_{22}$ $\displaystyle=d_{1}(g_{22}^{2}+g_{32}^{2}),$
$\displaystyle\dot{g}_{23}$ $\displaystyle=d_{1}(g_{22}g_{23}+g_{32}g_{33}),$
$\displaystyle\dot{g}_{24}$ $\displaystyle=d_{1}(g_{24}^{2}+g_{34}^{2}),$
$\displaystyle\dot{g}_{31}$
$\displaystyle=d_{1}(g_{21}g_{34}+g_{24}g_{31}+2g_{11}g_{32}+2g_{22}g_{31}-4g_{13}g_{33})-(g_{31}g_{42}+g_{32}g_{41}+g_{31}g_{32}),$
$\displaystyle\dot{g}_{32}$
$\displaystyle=d_{1}(2g_{22}g_{32}+2g_{24}g_{32}-2g_{14}g_{32})-(2g_{32}g_{42}+g_{32}^{2}),$
$\displaystyle\dot{g}_{33}$
$\displaystyle=d_{1}(3g_{22}g_{33}+g_{23}g_{32}+2g_{13}g_{32}-2g_{11}g_{33}-2g_{13}g_{31})-(g_{33}g_{44}+g_{34}g_{43}+g_{33}g_{34}),$
$\displaystyle\dot{g}_{34}$
$\displaystyle=d_{1}(4g_{24}g_{34}-2g_{14}g_{34})-(2g_{34}g_{44}+g_{34}^{2}),$
$\displaystyle\dot{g}_{41}$ $\displaystyle=-(g_{41}g_{42}+2g_{31}g_{32}),$
$\displaystyle\dot{g}_{42}$ $\displaystyle=-(g_{42}^{2}+2g_{32}^{2}),$
$\displaystyle\dot{g}_{43}$ $\displaystyle=-(g_{43}g_{44}+2g_{33}g_{34}),$
$\displaystyle\dot{g}_{44}$ $\displaystyle=-(g_{44}^{2}+2g_{34}^{2}).$

where $\dot{g_{i}}=dg_{i}/dy$ and $y=\nu_{0}\ln^{2}(\Lambda/T)$ with $\Lambda$
the UV cutoff and $\nu_{0}$ the density of states. The nesting parameters are

$d_{1}={\textstyle{d\Pi_{ph}(\bm{M}_{\alpha},T)\over
dy}}\approx{\textstyle{\Pi_{ph}(\bm{M}_{\alpha},T)\over\Pi_{pp}(0,T)}}.$ (24)

Note that $\Pi_{ph}(\bm{M}_{\alpha},T)$ contains $\ln^{2}$-divergence as well
and $d_{1}$ reaches the maximum value of ${\textstyle{1\over 2}}$ in the
perfect nesting case. From their definition we know that $d_{1}$ is related to
S/CDW.

It is interesting to note that, in the absence of SI effect, such that all
$g_{ij}$ are set by Hubbard $U$ initially, one can make take the condition
that $g_{ij}=g_{i}$. Then Eqs.(23) reduces to the well know RG equations for
the conventional 3-patch model,

$\displaystyle\dot{g}_{1}$ $\displaystyle=2d_{1}g_{1}(g_{2}-g_{1}),$ (25)
$\displaystyle\dot{g}_{2}$ $\displaystyle=d_{1}(g_{2}^{2}+g_{3}^{2}),$
$\displaystyle\dot{g}_{3}$
$\displaystyle=2d_{1}g_{3}(2g_{2}-g_{1})-(2g_{3}g_{4}+g_{3}^{2}),$
$\displaystyle\dot{g}_{4}$ $\displaystyle=-(g_{4}^{2}+2g_{3}^{2}).$

consistent with the 3-patch generalization of Ref.Furukawa , and also with
Ref.PhysRevB.100.085136 .

In the weak coupling limit, all these $g_{ij}$ behaves in a similar manner
close to the critical point $y_{c}$, namely,

$g_{ij}={\textstyle{G_{ij}\over y_{c}-y}}.$ (26)

To identify the leading order in the weak coupling limit, we first consider
the RG equations for various order parameters. We will consider charge-
density-wave, spin-density-wave, uniform superconductivity and pair-density-
wave. The renormalization for these vertices are shown in Fig. For the CDW,
the order parameter can be written as $\Delta_{\text{CDW}+}$ and
$\Delta_{\text{CDW}-}$, where we use $+$ and $-$ sign to denote the parity of
this order. The real and imaginary part of these CDW order parameters are
governed by the following RG equations,

$\displaystyle\dot{\Delta}^{R}_{\text{CDW}+}={\textstyle{d_{1}\over
2}}(g_{22}-g_{33}-2g_{11}+g_{23}+g_{32}-2g_{31}-2g_{13}){\Delta}^{R}_{\text{CDW}+}$
(27) $\displaystyle\dot{\Delta}^{R}_{\text{CDW}-}={\textstyle{d_{1}\over
2}}(g_{22}-g_{33}-2g_{11}-g_{23}-g_{32}+2g_{31}+2g_{13}){\Delta}^{R}_{\text{CDW}+}$
$\displaystyle\dot{\Delta}^{I}_{\text{CDW}+}={\textstyle{d_{1}\over
2}}(g_{22}+g_{33}-2g_{11}+g_{23}-g_{32}+2g_{31}-2g_{13}){\Delta}^{I}_{\text{CDW}+}$
$\displaystyle\dot{\Delta}^{I}_{\text{CDW}-}={\textstyle{d_{1}\over
2}}(g_{22}+g_{33}-2g_{11}-g_{23}+g_{32}-2g_{31}+2g_{13}){\Delta}^{I}_{\text{CDW}+}$

Similarly for SDW orders we have

$\displaystyle\dot{\Delta}^{R}_{\text{SDW}+}={\textstyle{d_{1}\over
2}}(g_{22}+g_{33}+g_{23}+g_{32}){\Delta}^{R}_{\text{SDW}+}$ (28)
$\displaystyle\dot{\Delta}^{R}_{\text{SDW}-}={\textstyle{d_{1}\over
2}}(g_{22}-g_{33}+g_{23}-g_{32}){\Delta}^{R}_{\text{SDW}+}$
$\displaystyle\dot{\Delta}^{I}_{\text{SDW}+}={\textstyle{d_{1}\over
2}}(g_{22}+g_{33}-g_{23}-g_{32}){\Delta}^{I}_{\text{SDW}+}$
$\displaystyle\dot{\Delta}^{I}_{\text{SDW}-}={\textstyle{d_{1}\over
2}}(g_{22}-g_{33}-g_{23}+g_{32}){\Delta}^{I}_{\text{SDW}+}$

For the uniform superconducting order, we obtain $s$\- and $f$-wave orders
from the six-patch model. In addition, we also obtain two degenerate $d$-wave
and two degenerate $p$-wave. The order parameters obey the following RG
equations,

$\displaystyle\dot{\Delta}_{s\text{SC}}$ $\displaystyle=-{\textstyle{1\over
2}}(2g_{31}+2g_{32}+g_{41}+g_{42}){\Delta}_{s\text{SC}}$ (29)
$\displaystyle\dot{\Delta}_{f\text{SC}}$ $\displaystyle=-{\textstyle{1\over
2}}(-2g_{31}+2g_{32}-g_{41}+g_{42}){\Delta}_{f\text{SC}}$
$\displaystyle\dot{\Delta}_{p\text{SC}}$ $\displaystyle=-{\textstyle{1\over
2}}(g_{31}-g_{32}-g_{41}+g_{42}){\Delta}_{p\text{SC}}$
$\displaystyle\dot{\Delta}_{d\text{SC}}$ $\displaystyle=-{\textstyle{1\over
2}}(-g_{31}-g_{32}+g_{41}+g_{42}){\Delta}_{d\text{SC}}$

Note that using Eq.(26) all the RG equations for the order parameters can be
written in the form

$\dot{\Delta}=\Gamma(g_{ij})\Delta=-\beta(G_{ij})\Delta/(y_{c}-y).$ (30)

Solving this equation, one obtains,

$\Delta(y)\propto(y_{c}-y)^{\beta(G_{ij})}$ (31)

The susceptibility for the corresponding order parameters behaves as

$\dot{\chi}=\tilde{d}|\Delta|^{2}$ (32)

where $\tilde{d}=1$ for uniform superconducting orders, $\tilde{d}=d_{s}$ for
PDW and $\tilde{d}=d$ for charge or spin density waves. It can be inferred
directly that

$\chi(y)\propto(y_{c}-y)^{\gamma(G_{ij})}$ (33)

and

$\gamma(G_{ij})=2\beta(G_{ij})+1.$ (34)

The leading order with be the one which has the most negative
$\gamma(G_{ij})$. Below are the explicit expressions for $\gamma$:

$\displaystyle\gamma_{\text{CDW}+}^{R}$
$\displaystyle=-d_{1}(G_{22}-G_{33}-2G_{11}+G_{23}+G_{32}-2G_{31}-2G_{13})+1$
(35) $\displaystyle\gamma_{\text{CDW}-}^{R}$
$\displaystyle=-d_{1}(G_{22}-G_{33}-2G_{11}-G_{23}-G_{32}+2G_{31}+2G_{13})+1$
$\displaystyle\gamma_{\text{CDW}+}^{I}$
$\displaystyle=-d_{1}(G_{22}+G_{33}-2G_{11}+G_{23}-G_{32}+2G_{31}-2G_{13})+1$
$\displaystyle\gamma_{\text{CDW}-}^{I}$
$\displaystyle=-d_{1}(G_{22}+G_{33}-2G_{11}-G_{23}+G_{32}-2G_{31}+2G_{13})+1$
$\displaystyle\gamma_{\text{SDW}+}^{R}$
$\displaystyle=-d_{1}(G_{22}+G_{33}+G_{23}+G_{32})+1,~{}~{}\gamma_{\text{SDW}-}^{R}=-d_{1}(G_{22}-G_{33}+G_{23}-G_{32})+1$
$\displaystyle\gamma_{\text{SDW}+}^{I}$
$\displaystyle=-d_{1}(G_{22}+G_{33}-G_{23}-G_{32})+1,~{}~{}\gamma_{\text{SDW}-}^{I}=-d_{1}(G_{22}-G_{33}-G_{23}+G_{32})+1$
$\displaystyle\gamma_{s\text{SC}}$
$\displaystyle=(2G_{31}+2G_{32}+G_{41}+G_{42})+1,~{}~{}\gamma_{f\text{SC}}=(-2G_{31}+2G_{32}-G_{41}+G_{42})+1$
$\displaystyle\gamma_{p\text{SC}}$
$\displaystyle=(G_{31}-G_{32}-G_{41}+G_{42})+1,~{}~{}\gamma_{d\text{SC}}=(-G_{31}-G_{32}+G_{41}+G_{42})+1$

As a result, we need to first solve the RG equations for $g_{ij}$ to obtain
$G_{ij}$, and then obtain $\gamma$’s from the above equation to identify the
leading weak coupling instability.

Figure 5: Solutions of the RG equation with finite interactions. (a), when
there is no SI effect, then all the interactions are set by the largest $U$.
In this case, all $g_{ij}$ are equal to $g_{i}$, and the RG result indicate
the most divergent order is the $d$-wave SC order, consistent with the 3-patch
model analysis for doped graphene in Ref.Nandkishore2012 . (b), with SI effect
and when keeping only the largest $U$ and $V$ as initial input, the RG flows
show that $U$ decays and $V$ grows. The resulting leading instability is
degenerate complex CDW order with odd parity ($\Delta_{\text{CDW}-}$). (c) and
(d), when the exchange interaction $J_{0}$ is further taken into account, the
RG flows show that $J_{0}$ also grows, but is subleading to $V$. The leading
instabilities in this situation is imaginary SDW with even and odd parity for
$J_{0}>0$ and the complex CDW with odd parity for $J_{0}<0$. (e) and (f), when
a finite initial $K_{0}$ is added into the RG flow, the RG results shows that
near the critical energy scale, the presence of $K_{0}$ spoils the constant
map between $g_{ij}$ and lattice interactions, but for $y<y_{c}$ the constant
map is still a good approximation. Although we cannot identify the flow of
each lattice interaction, we are still able to find the leading instabilities
in weak coupling limit, which can be $s$-wave SC or real CDW order with odd
parity.

In Fig.5 we show the weak coupling RG flows of the interactions for different
situations, together with the flows of the order parameter vertex functions
$\Gamma$ which help to illustrate the leading instability.

We first show the results without the SI effect in Fig.5(a), which can be
realized at the m-type Van Hove filling for the Kagome lattice, or other
hexagonal lattices such as triangular or honeycomb lattice. In these cases,
all the interactions are set by the largest Hubbard $U$, and therefore
$g_{ij}$ for all $j=1,..,4$ are identical. Then the resulting leading
instability is the $d$-wave uniform superconductivity, consistent with the
findings in Ref.Nandkishore2012 . The situation is quite different in the
presence of SI effect. In Fig.5(b) we show the results for $U$ and $V$ only
model. In this particular case, we are able to find a constant map between the
lattice interactions and the band fermion interactions, which clearly shows
that as the energy scale is lowered, the $U$ decays and $V$ grows and we can
have a low energy effective model where $V\gg U$. The resulting weak coupling
instability is a complex CDW order with odd parity, and imaginary SDW order
with both even and odd parity– all of these orders are degenerate.

We next add some small but finite exchange interaction $J_{0}$ at the bare
level and the results are shown in Fig.5(c) and (d). If we only include
$J_{0}$ but not $K_{0}$, we see that there is still a constant map between
$g_{ij}$ and the lattice interactions, which clearly show that $U$ decays and
both $V$ and $J_{0}$ increases, but $V$ still remains the largest. The weak
coupling instabilities are imaginary CDW with both even and odd parity for
$J_{0}>0$ and the complex CDW order with odd parity for $J_{0}<0$.

In Fig.5 (e) and (f) we show the RG results when all $U$, $V$, $J_{0}$ and
$K_{0}$ are present. We see that the presence of a finite $K_{0}$ spoils the
constant map between lattice interaction and $g_{ij}$ especially near the
critical $y_{c}$. However, this destruction is minor at some intermediate
energy scale and when the input $K_{0}$ is smaller. In this situation we can
still have some approximate constant map between lattice interaction and
$g_{ij}$, and from this we again can easily see that $V$ is the largest at
some intermediate energy scale. In the weak coupling limit when the onset of
instability is close to $y_{c}$, we can still extract the leading instability,
for which we obtain $s$-wave uniform superconductivity with the parameters in
Fig.5(e) and real CDW order with odd parity with the parameters in Fig.5(f).

We emphasize that in all these RG results, the constant map between lattice
interaction and $g_{ij}$ is present up to some intermediate energy scale, for
which we have $V$ as the largest interaction. This observation makes us
confident of using $V$-only model to study the competing orders in the strong
coupling case.

##  IV Calculation of the bare particle-hole and particle-particle bubble

Near the van Hove singularities at $\bm{M}_{\alpha}$ points, we can focus on
the small patches around them. Here without loss of generality, we study the
patches around $\bm{M}_{A}=(0,2\pi/\sqrt{3})$ and
$-\bm{M}_{B}=(\pi,\pi/\sqrt{3})$. We denote the dispersions within these two
patches as $\xi_{1}({\bm{k}})$ and $\xi_{2}({\bm{k}})$ respectively. The
particle-hole bubble and the particle-particle bubble we will evaluate are
given by

$\displaystyle\Pi_{pp}^{(0)}$
$\displaystyle=T\sum_{n}\int{\textstyle{d{\bm{k}}\over
4\pi^{2}}}{\textstyle{1\over
i\omega_{n}-\xi_{1}({\bm{k}})}}{\textstyle{1\over-i\omega_{n}-\xi_{2}({\bm{k}})}}$
(36) $\displaystyle\Pi_{ph}^{(0)}$
$\displaystyle=T\sum_{n}\int{\textstyle{d{\bm{k}}\over
4\pi^{2}}}{\textstyle{1\over i\omega_{n}-\xi_{1}({\bm{k}})}}{\textstyle{1\over
i\omega_{n}-\xi_{2}({\bm{k}})}}$

Near the vHS points, and in the model with only nearest neighbor hopping, we
can approximate the dispersion as

$\displaystyle\xi_{1}({\bm{k}})$ $\displaystyle={\textstyle{t\over
2}}(k_{x}^{2}-3k_{y}^{2})=ta_{+}a_{-},$ (37) $\displaystyle\xi_{2}({\bm{k}})$
$\displaystyle=-t(k_{x}^{2}+\sqrt{3}k_{x}k_{y})=-ta_{+}(a_{+}+a_{-}),$

where $a_{+}=(k_{x}-\sqrt{3}k_{y})/\sqrt{2}$ and
$a_{-}=(k_{x}+\sqrt{3}k_{y})/\sqrt{2}$ Using this, it is easy to see

$\displaystyle\Pi_{pp}^{(0)}$ $\displaystyle={\textstyle{-T\over
4\sqrt{3}\pi^{2}t}}\sum_{n}\int_{-\sqrt{\Lambda}}^{\sqrt{\Lambda}}da_{+}da_{-}{\textstyle{1\over
i\omega_{n}-a_{+}a_{-}}}{\textstyle{1\over i\omega_{n}-a_{+}^{2}-a_{+}a_{-}}}$
(38) $\displaystyle={\textstyle{-T\over
4\sqrt{3}\pi^{2}t}}\sum_{n}\int_{-\sqrt{\Lambda}}^{\sqrt{\Lambda}}{\textstyle{da_{+}\over
a_{+}^{3}}}\int_{-\sqrt{\Lambda}a_{+}}^{\sqrt{\Lambda}a_{+}}dx\left({\textstyle{1\over
i\omega_{n}-x-a_{+}^{2}}}-{\textstyle{1\over i\omega_{n}-x}}\right)$

here we measure both $\omega_{n}$ and $T$ in units of $t$, and we have changed
the variable as $x=a_{+}a_{-}$ in the second line. For the integration over
$x$, we can shift the integration variables to get

$\displaystyle\int_{-\sqrt{\Lambda}a_{+}}^{\sqrt{\Lambda}a_{+}}dx\left({\textstyle{1\over
i\omega_{n}-x-a_{+}^{2}}}-{\textstyle{1\over i\omega_{n}-x}}\right)$ (39)
$\displaystyle=$
$\displaystyle\int_{-\sqrt{\Lambda}a_{+}+a_{+}^{2}}^{\sqrt{\Lambda}a_{+}+a_{+}^{2}}{\textstyle{dx\over
i\omega_{n}-x}}-\int_{-\sqrt{\Lambda}a_{+}}^{\sqrt{\Lambda}a_{+}}{\textstyle{dx\over
i\omega_{n}-x}}$ $\displaystyle=$
$\displaystyle\int_{\sqrt{\Lambda}a_{+}}^{\sqrt{\Lambda}a_{+}+a_{+}^{2}}{\textstyle{dx(-x-i\omega_{n})\over
x^{2}+\omega_{n}^{2}}}-\int_{-\sqrt{\Lambda}a_{+}}^{-\sqrt{\Lambda}a_{+}+a_{+}^{2}}{\textstyle{dx(-x-i\omega_{n})\over
x^{2}+\omega_{n}^{2}}}$ $\displaystyle=$
$\displaystyle-\int_{\sqrt{\Lambda}a_{+}-a_{+}^{2}}^{\sqrt{\Lambda}a_{+}+a_{+}^{2}}{\textstyle{dxx\over
x^{2}+\omega_{n}^{2}}}$

Note we have used the fact that the imaginary part is odd in $\omega_{n}$ and
will not survive after frequency summation, and is thus dropped. In low $T$
limit, we can replace $T\sum_{n}$ with
$\int_{-\infty}^{\infty}{\textstyle{d\omega_{n}\over 2\pi}}$ and carry out the
frequency integration, this results in $\int
d\omega_{n}1/(x^{2}+\omega_{n}^{2})=\pi/|x|$. Using this, it is easy to see

$\displaystyle\Pi_{pp}^{(0)}$ $\displaystyle={\textstyle{1\over
4\sqrt{3}\pi^{2}t}}\int_{-\sqrt{\Lambda}}^{\sqrt{\Lambda}}{\textstyle{da_{+}\over
a_{+}^{3}}}a_{+}^{2}\text{sgn}(a_{+})$ (40)

Imposing an IR cutoff $\sqrt{T}$ for the integration, we obtain $\ln\Lambda/T$
for the result. Now if we assume that there is a small deviation from the van
Hove filling, characterized by a nonzero $\mu$, we end up with the following
expression,

$\Pi_{pp}^{(0)}={\textstyle{1\over
4\sqrt{3}\pi^{2}t}}\ln{\textstyle{\Lambda\over\max\\{T,|\mu|\\}}}$ (41)

For the particle-hole bubble, similar analysis shows that

$\Pi_{ph}^{(0)}={\textstyle{1\over
8\sqrt{3}\pi^{2}t}}\ln{\textstyle{\Lambda\over\max\\{T,|\mu|\\}}}\ln{\textstyle{\Lambda\over\max\\{T,\mu,|t^{\prime}|\\}}}$
(42)

Comparing to $\Pi_{pp}^{(0)}$, there is an additional $\ln$ coming from the
nesting of the Fermi surface in the particle-hole channel, which can be cut by
either a finite shift from the van Hove filling ($\mu\neq 0$), or including
the next nearest neighbor hopping term $t^{\prime}\neq 0$.


# On the solutions of second order difference equations with variable
coefficients

S. Kadyrov  Suleyman Demirel University
Mathematics and Natural Sciences
1/1 Abylay Khan str
Kaskelen
Kazakhstan<EMAIL_ADDRESS>and  A. Orynbassar  Suleyman Demirel
University
Mathematics and Science Educations
1/1 Abylay Khan str
Kaskelen
Kazakhstan<EMAIL_ADDRESS>
###### Abstract.

In this article we study solutions to second order linear difference equations
with variable coefficients. Under mild conditions we provide closed form
solutions using finite continued fraction representations. The proof of the
results are elementary and based on factoring a quadratic shift operator. As
an application, we obtain two new generalized continued fraction formulas for
the mathematical constant $\pi^{2}$.

###### Key words and phrases:

difference equations, continued fractions, non-homogenous equations, closed
form solutions, Apéry-like series, $\pi$

###### 1991 Mathematics Subject Classification:

39A06; 11J70

Both authors acknowledge the support by a grant from the Ministry of Education
and Science of the Republic of Kazakhstan within framework of the project
AP08051987.

##  1\. Introduction

Compared to differential equations, the theory of difference equations is a
less studied fields of mathematics. In this note, we consider second order
linear difference equations

(1.1) $y_{n}-b_{n}y_{n-1}-a_{n}y_{n-2}=f_{n},\,\,n\geq 1$

with variable coefficients $(b_{n})$ and $(a_{n})$ in complex numbers. When
coefficients are constants, (1.1) can be easily solved with representing the
solutions in the form $y_{n}=x^{n}$ and finding the roots of the corresponding
characteristic equation, see e.g. [4]. On the other hand, when the
coefficients are non-constant, the general solution representation becomes
very complicated [11, 10, 9, 1, 3] and far from providing any useful insights.
It is known that (1.1) has two linearly independent solutions and general
solution can be written as a linear combination of the two [6]. Besides, when
initial values $y_{0}$ and $y_{-1}$ are given, clearly the solution is unique.

In this note, we aim to provide a closed form solution for (1.1) using finite
continued fractions and use it to prove two conjectures from
www.ramanujanmachine.com. Our main result is the following.

######  Theorem 1.

Assume that there exist sequences $(c_{n})$ and $(d_{n})$ of complex numbers
satisfying

(1.2) $\displaystyle c_{n}+d_{n+1}$ $\displaystyle=b_{n},$ (1.3)
$\displaystyle c_{n}d_{n}$ $\displaystyle=-a_{n},$

for any $n\geq 1.$ Then, the solution to the second order equation (1.1) with
initial values $y_{-1}$ and $y_{0}$ is given by

(1.4)
$y_{n}=\sum_{i=1}^{n+1}\left(\sum_{j=1}^{i-1}f_{j}\prod_{k=j+1}^{i-1}c_{k}+(y_{0}-d_{1}y_{-1})\prod_{k=1}^{i-1}c_{k}\right)\prod_{j=i+1}^{n+1}d_{j}+y_{-1}\prod_{j=1}^{n+1}d_{j},\,\,\,\forall
n\geq 1,$

The idea of the proof is based on factorization of the shift operator and is
not new, see e.g. [6, 12]

As an immediate consequence we may let $f_{n}=0$ to obtain the following
solution for the corresponding homogeneous equation.

######  Theorem 2.

Assume that there exist sequences $(c_{n})$ and $(d_{n})$ of complex numbers
satisfying (1.2) and (1.3) for any $n\geq 1.$ Then, the solution to the second
order homogeneous equation

(1.5) $y_{n}-b_{n}y_{n-1}-a_{n}y_{n-2}=0,\,\,n\geq 1$

with initial values $y_{-1}$ and $y_{0}$ is given by

(1.6)
$y_{n}=(y_{0}-d_{1}y_{-1})\sum_{i=1}^{n+1}\prod_{j=1}^{i-1}c_{i}\prod_{j=i+1}^{n+1}d_{j}+y_{-1}\prod_{j=1}^{n+1}d_{j},\,\,\,\forall
n\geq 1.$

We remark that inductively it easy to see that the sequences $(c_{n})$ and
$(d_{n})$ satisfying (1.2) and (1.3) exist for most sequences $(a_{n})$ and
$(b_{n})$.

We also note that, we have flexibility in choosing $d_{1}$. This does not mean
that the difference equation (1.1) has more than one solutions. Instead, it
gives different representation of the same solution. In particular, when
$y_{-1}$ is non-zero, we take $d_{1}$ such that $y_{0}-d_{1}y_{-1}=0$ yielding
to a simpler represented solution.

######  Corollary 1.

Let $(c_{n})$ and $(d_{n})$ be given as in Theorem 1. If $y_{-1}\neq 0$, then
the difference equation (1.1) has the solution

(1.7) $y_{n}=y_{-1}\prod_{j=1}^{n+1}d_{j},\,\,\,\forall n\geq 1,$

with $d_{1}=y_{0}/y_{-1}.$

In order to provide the closed form solutions, one needs to find the closed
form solutions to (1.2) and (1.3). In the special case, when $a_{n}\neq 0$ for
any $n\geq 1,$ we clearly have $d_{n}\neq 0$ for $n\geq 0$ and can take
$d_{1}\neq 0.$ In this case, substituting into (1.3) we get
$(b_{n}-d_{n+1})d_{n}=-a_{n}.$ Solving for $d_{n+1}$ we get the recursive
relation

$d_{n+1}=b_{n}+\frac{a_{n}}{d_{n}}.$

Iterating the equation we have generalized continued fraction representations
for the sequences $(c_{n})$ and $(d_{n})$ as follows.

######  Lemma 1.

Fix any non-zero complex number $d_{1}$ and for $n\geq 1$ define

(1.8)
$d_{n+1}=b_{n}+\cfrac{a_{n}}{b_{n-1}+\cfrac{a_{n-1}}{b_{n-2}+\cdots+\cfrac{a_{2}}{b_{1}+\cfrac{a_{1}}{d_{1}}}}}\text{
and
}c_{n}=-\cfrac{a_{n}}{b_{n-1}+\cfrac{a_{n-1}}{b_{n-2}+\cdots+\cfrac{a_{2}}{b_{1}+\cfrac{a_{1}}{d_{1}}}}}.$

Then, provided the sequence $(d_{n})$ is well-defined, $(d_{n})$ and $(c_{n})$
satisfy (1.2) and (1.3).

Next, we turn to applying Theorem 1 to obtain new generalized continued
fraction representations to $\pi^{2}$. Representing various mathematical
constants with generalized continued fractions is not new. However, recently a
machine learning based computer algorithm is developed to generate lots of
seemingly new conjectures on continued fraction representations [14, 13].
These generated conjectures are posted at www.ramanujanmachine.com attracting
new investigations in the theory of continued fractions [7, 5, 2, 8, 15].
Using Theorem 2 we prove the following results stated as conjectures in the
site.

######  Theorem 3.

We have

$\frac{8}{\pi^{2}}=b_{0}+\cfrac{a_{1}}{b_{1}+\cfrac{a_{2}}{b_{2}+\ddots}},$

for the sequences given by $b_{n}=3n(n+1)+1,n\geq 0$ and
$a_{n}=−(2n−1)n^{3},n\geq 1.$

######  Theorem 4.

We have

$\frac{18}{\pi^{2}}=b_{0}+\cfrac{a_{1}}{b_{1}+\cfrac{a_{2}}{b_{2}+\ddots}},$

for the sequences given by $b_{n}=n(5n+6)+2,n\geq 0$ and
$a_{n}=−4n^{4}+2n^{3},n\geq 1.$

##  2\. Proof of Theorem 1

In this section we prove our main result. To this end we introduce the _shift
operator_ $T$ given by

$Ty_{n}=y_{n+1},$

for any sequence. We note that (1.1) can be rewritten with shift operator $T$
as follows

(2.1) $y_{n}-b_{n}y_{n-1}-a_{n}y_{n-2}=(T^{2}-b_{n}T-a_{n})y_{n-2}=0.$

######  Lemma 2.

Let $(c_{n})$ and $(d_{n})$ be sequences as in Theorem 1. Then, (1.1) is
equivalent to

(2.2) $(T-c_{n})(T-d_{n})y_{n-2}=0.$

###### Proof.

We have

$\displaystyle(T-c_{n})(T-d_{n})y_{n-2}$
$\displaystyle=(T-c_{n})(y_{n-1}-d_{n}y_{n-2})$
$\displaystyle=y_{n}-(d_{n+1}+c_{n})y_{n-1}+c_{n}d_{n}y_{n-2}.$

Comparing with (1.1) we see that $d_{n+1}+c_{n}=b_{n}$ and $c_{n}d_{n}=-a_{n}$
as required. ∎

We note non-commutative nature of $T$ that in general
$(T-c_{n})(T-d_{n})y_{n-2}=(T-d_{n})(T-c_{n})y_{n-2}$, which is the reason why
(1.2) is not symmetric. We now prove Theorem 1.

###### Proof of Theorem 1.

Let $(c_{n})$ and $(d_{n})$ be given as in the theorem statement. Using Lemma
2 we rewrite (1.1) as $(T-c_{n})(T-d_{n})y_{n-2}=f_{n}$. Let us define
$(w_{n})_{n\geq-1}$ by

(2.3) $w_{n-2}:=(T-d_{n})y_{n-2}$

so that $(T-c_{n})w_{n-2}=f_{n}.$ The latter gives
$w_{n-1}=f_{n}+c_{n}w_{n-2}$ for any $n\geq 1$. Iteration yields

$w_{n-1}=\sum_{i=0}^{n}f_{i}\prod_{j=i+1}^{n}c_{j}+w_{-1}\prod_{j=1}^{n}c_{j},$

where $w_{-1}$ is obtained from (2.3).

$w_{-1}=y_{0}-d_{1}y_{-1}.$

Now, applying (2.3) we arrive at a first order difference equation with
variable coefficients

$y_{n-1}=w_{n-2}+d_{n}y_{n-2}.$

Iterating yields the closed form solution

$y_{n-1}=\sum_{i=1}^{n}w_{i-2}\prod_{j=i+1}^{n}d_{j}+y_{-1}\prod_{j=1}^{n}d_{j}.$

This finishes the proof. ∎

##  3\. Generalized continued fraction representations of $\pi^{2}$

It is well-known that $n$’th approximants $A_{n}$ and $B_{n}$ of a generalized
continued fractions satisfy

$A_{n}=b_{n}A_{n-1}+a_{n}A_{n-2}\text{ and }B_{n}=b_{n}B_{n-1}+a_{n}B_{n-2},$
$A_{-1}=1,A_{0}=b_{0},B_{-1}=0,B_{0}=1.$

###### Proof of Theorem 3.

Since $A_{-1}\neq 0,$ we may apply Corollary 1 with
$d_{1}=A_{0}/A_{-1}=b_{0}=1$ to obtain

$A_{n}=\prod_{j=1}^{n+1}d_{j},$

with $d_{j}$ as in (1.8). We have $d_{2}=b_{1}+a_{1}=7-1=6$. Inductively, if
$d_{n}=n(2n-1)$ we see that

$d_{n+1}=b_{n}+\frac{a_{n}}{d_{n}}=3n(n+1)+1-\frac{(2n−1)n^{3}}{n(2n-1)}=(n+1)(2n+1).$

Hence, $A_{n}=\prod_{j=1}^{n+1}j(2j-1)=(n+1)!(2n+1)!!=(2(n+1))!/2^{n+1}.$

As for $B_{n}$ we have $B_{-1}=0$ and $B_{0}=1$ which gives the solution

$B_{n}=\sum_{i=1}^{n+1}\prod_{j=1}^{i-1}c_{i}\prod_{j=i+1}^{n+1}d_{j}.$

We may take $d_{1}=1$ as before and arrive at

$d_{n}=n(2n-1)=2n(2n-1)/2,$

which in turn gives

$c_{n}=b_{n}-d_{n+1}=3n^{2}+3n+1-2n^{2}-3n-1=n^{2}.$

Thus,

$\prod_{j=1}^{i-1}c_{i}\prod_{j=i+1}^{n+1}d_{j}=\frac{(i-1)!^{2}(n+1)!(2n+1)!!}{i!(2i-1)!!}=A_{n}\frac{(i-1)!^{2}}{i!(2i-1)!!}=A_{n}\frac{i!^{2}2^{i}}{i^{2}(2i)!},$

where we use the fact that $i!(2i-1)!!=(2i)!/2^{i}.$ Hence,
$B_{n}=A_{n}\sum_{i=1}^{n+1}\frac{i!^{2}2^{i}}{i^{2}(2i)!}$. Thus, using [16,
Eq.3.78] we arrive at

$\lim_{n\to\infty}\frac{A_{n}}{B_{n}}=\frac{1}{\sum_{i=1}^{\infty}\frac{i!^{2}2^{i}}{i^{2}(2i)!}}=\frac{8}{\pi^{2}}.\qed$

###### Proof of Theorem 4.

The proof of Theorem 4 is very similar to the above and we omit the details.
Taking $d_{1}=A_{0}/A_{-1}=b_{0}=2$ it is easy to see that $d_{n}=2n(2n-1)$
which leads to

$A_{n}=\prod_{j=1}^{n+1}d_{j}=(2n+2)!.$

The same $d_{n}$ can be used to find $B_{n}$. In this case, we see that
$c_{n}=b_{n}-d_{n+1}=n^{2}$ which gives

$B_{n}=\sum_{i=1}^{n+1}\prod_{j=1}^{i-1}c_{i}\prod_{j=i+1}^{n+1}d_{j}=\sum_{i=1}^{n+1}\frac{(i-1)!^{2}(2n+2)!}{(2i)!}=A_{n}\sum_{i=1}^{n+1}\frac{(i-1)!^{2}}{(2i)!}.$

Together with [16, Eq.3.75] we arrive at

$\lim_{n\to\infty}\frac{A_{n}}{B_{n}}=\frac{1}{\sum_{i=1}^{\infty}\frac{(i-1)!^{2}}{(2i)!}}=\frac{18}{\pi^{2}}.\qed$

##  4\. Discussion and Conclusion

In this note, we proposed a different approach than [9] to solve second order
homogeneous equations with variable coefficients. The main idea is to factor
quadratic shift operator into linear operators and represent the solution with
terms involving finite continued fractions. Similar methods can be studied for
higher order equations.

While the proposed approach is simple and effective in applications to
generalized continued fraction representations, it does not address all second
order linear difference equations. Thus, it is interesting to find out when
this method works, more specifically, to investigate a necessary and
sufficient condition for the sequences $c_{n}$ and $d_{n}$ satisfying (1.2)
and (1.3) to exist.

Finally, we note that while solving second order homogeneous equations with
variable coefficients helps to reduce the infinite generalized continued
fractions into series representations it does not necessarily help to find the
limit. Hence, unless the sum of the series representation is known the problem
of representing a number with continued fractions is not fully solved.

## References

  * [1] B. Buchberger, “O velichinakh $s_{i}$, opredelennykh rekursiey $s_{i+1}=u_{i}s_{i}+v_{i}s_{i-1}$,” _Tech. Rep._ , 1971.
  * [2] R. Dougherty-Bliss and D. Zeilberger, “Automatic conjecturing and proving of exact values of some infinite families of infinite continued fractions,” _The Ramanujan Journal_ , pp. 1–17, 2021.
  * [3] S. Elaydi, _An Introduction to Difference Equations_. Springer, Berlin, 2005.
  * [4] A. J. Jerri, _Linear difference equations with discrete transform methods_. Springer Science & Business Media, 2013, vol. 363.
  * [5] S. Kadyrov and F. Mashurov, “Generalized continued fraction expansions for $\pi$ and $e$,” _arXiv preprint arXiv:1912.03214_ , 2019.
  * [6] W. G. Kelley and A. C. Peterson, _Difference equations: an introduction with applications_. Academic press, 2001\.
  * [7] Z. Lu, “Elementary proofs of generalized continued fraction formulae for $e$,” _arXiv preprint arXiv:1907.05563_ , 2019.
  * [8] P. Lynch, “Derangements and continued fractions for $e$,” _arXiv preprint arXiv:2012.12692_ , 2020.
  * [9] R. K. Mallik, “On the solution of a second order linear homogeneous difference equation with variable coefficients,” _Journal of Mathematical Analysis and Applications_ , vol. 215, no. 1, pp. 32–47, 1997.
  * [10] R. K. Mallik, “Solutions of linear difference equations with variable coefficients,” _Journal of mathematical analysis and applications_ , vol. 222, no. 1, pp. 79–91, 1998.
  * [11] R. K. Mallik, “On the solution of a linear homogeneous difference equation with variable coefficients,” _SIAM Journal on Mathematical Analysis_ , vol. 31, no. 2, pp. 375–385, 2000.
  * [12] I. Parasidis and E. Providas, “Factorization method for the second order linear nonlocal difference equations,” in _International Conference Polynomial Computer Algebra, Saint-Petersburg, Russia, Euler International Mathematical Institute, April_ , no. 16-22, 2018, p. 84.
  * [13] G. Raayoni, S. Gottlieb, Y. Manor, G. Pisha, Y. Harris, U. Mendlovic, D. Haviv, Y. Hadad, and I. Kaminer, “Generating conjectures on fundamental constants with the ramanujan machine,” _Nature_ , vol. 590, no. 7844, pp. 67–73, 2021\.
  * [14] G. Raayoni, S. Gottlieb, G. Pisha, Y. Harris, Y. Manor, U. Mendlovic, D. Haviv, Y. Hadad, and I. Kaminer, “The ramanujan machine: automatically generated conjectures on fundamental constants,” _arXiv preprint arXiv:1907.00205_ , 2019.
  * [15] C. Reyes-Bustos, “A simple continued fraction expansion for $e^{n}$,” _arXiv preprint arXiv:1909.13597_ , 2019.
  * [16] T. Sherman, “Summation of Glaisher and Apéry-like series,” 2000.