#  A Rigorous Formalism of Unconventional Symmetry Breaking in Fermi Liquid
Theory
and Its Application to Nematicity in FeSe

Rina Tazai, Shun Matsubara, Youichi Yamakawa, Seiichiro Onari, and Hiroshi
Kontani  Department of Physics, Nagoya University, Furo-cho, Nagoya 464-8602,
Japan.

###### Abstract

Unconventional symmetry breaking due to nonlocal order parameters has
attracted considerable attention in many strongly correlated metals. Famous
examples are the nematic order in Fe-based superconductors and the star-of-
David charge density order in kagome metals. Such exotic symmetry breaking in
metals is a central issue of modern condensed matter physics, while its
theoretical foundation is still unclear in comparison with the well-
established theory of superconductivity. To overcome this difficulty, here we
introduce the “form factor” that generalizes the nonlocal order parameter into
the Luttinger-Ward (LW) Fermi liquid theory. We then construct a rigorous
formalism of the “density-wave equation” that gives the thermodynamically
stable form factor, similarly to the superconducting-gap equation. In
addition, a rigorous expression of the Ginzburg-Landau free-energy for the
unconventional order is presented to calculate various thermodynamic
properties. In the next stage, we apply the derived formalism to a typical Fe-
based superconductor FeSe, by using the one-loop LW function that represents
the free-energy gain due to the interference among paramagnons. The following
key experiments are naturally explained: (i) Lifshitz transition
(=disappearance of an electron-pocket) due to the bond+orbital order below
$T_{\rm c}$. (ii) Curie-Weiss behavior of the nematic susceptibility at higher
$T$, and the deviation from the Curie-Weiss behavior at lower $T$ near the
nematic quantum-critical-point. (iii) Scaling relation of the specific heat
jump at $T_{\rm c}$, $\Delta C/T_{\rm c}\propto T_{\rm c}^{b}$ with $b\sim 3$.
(Note that $b=0$ in the BCS theory.) These results lead to a conclusion that
the nematicity in FeSe is the bond+orbital order due to the “paramagnon
interference mechanism”. The present theory paves the way for solving various
unconventional phase transition systems.

##  I Introduction

Recently, rich symmetry breaking phenomena due to unconventional order
parameters have attracted considerable attention in various electron systems.
Famous examples are the $C_{2}$-symmetric nematic order in various Fe-based
superconductors Chubukov-nematic-rev ; Yamakawa-PRX2016 ; Onari-FeSe ;
Fernandes-Nature . It has been established that the nematic state is driven by
the electron-correlation, thanks to the “electronic nematic susceptibility”
measurements performed by the shear-modulus analysis Yoshizawa ; Fernandes1 ;
Bohmer , the Raman spectroscopy Gallais ; Gallais2 ; Raman-spectroscopy , and
the elastoresistivity measurements Fisher-review ; Hosoi ; B2g-Ishida ;
Fisher-Science2016 . Similar nematic orders are also observed in magic-angle-
twisted-bilayer-graphene Fernandes-TBG ; Onari-TBG , titanium pnictide oxide
Nakaoka-Ti , and cuprate superconductors Sato-CDW ; Murayama-CDW ; Kawaguchi-
CDW ; Tsuchiizu-CDW . Also, correlation-driven density-wave (DW) with nonzero
wavevector (${{\bm{q}}}\neq{\bm{0}}$) has attracted increasing attention
recently. Famous examples are the star-of-David DW state in kagome metal
CsV3Sb5 Tazai-kagome ; kagome-Thomale2013 ; kagome-SMFRG ; kagome-Balents2021
the CDW states in cuprate superconductors Yamakawa-CDW ; Kawaguchi-CDW ;
Tsuchiizu-CDW ; Sachdev-CDW1 ; Husemann ; Kivelson-spin-nematic ; DHLee-spin-
nematic , and the multopole DW states in heavy fermion systems Tazai-CeB6 ;
Tazai-CeCu2Si2-1 ; Tazai-CeCu2Si2-2 . Furthermore, more exotic odd-parity DW
orders that accompany the charge or spin loop current have been discovered in
kagome metals mSR_TRS ; Kerr , iridates Murayama , and cuprates Varma .

We call these DW states “unconventional” because they have non-local and
non-$s$-wave order parameters, in analogy to unconventional (non-$s$-wave)
superconductivity.  For example, the order parameter of the $d$-wave bond
order is
$O_{i,j}=\bar{O}(\delta_{|x-x^{\prime}|,1}\delta_{y-y^{\prime},0}-\delta_{x-x^{\prime},0}\delta_{|y-y^{\prime}|,1})$,
where $(x,y)$ and $(x^{\prime},y^{\prime})$ are the integer coordinates of $i$
and $j$ sites Nersesyan ; Kivelson-NJP ; Halboth:2000tt ; Davis:2013ce ;
Yamakawa-CDW ; Kawaguchi-CDW ; Tsuchiizu-CDW . (In high contrast, the
conventional magnetic order $m_{i}=n_{i\uparrow}-n_{i\downarrow}$ is local.)
However, simple local spin-density-wave (SDW) is inevitably derived within the
mean-field-level approximations for the Hubbard models with screened Coulomb
interactions Tazai-JPSJ-fRG ; Kontani-sLC . Therefore, beyond-mean-field many-
electron theories are necessary to understand the unconventional DW states.
This is a difficult but very interesting theoretical problem, and this is a
central issue of modern condensed matter physics. On the other hand, these
unconventional DW states we are interested in are metallic, so the itinerant
picture will be fruitful. In addition, the $U(1)$ gauge symmetry is preserved.
Thus, it is promising to construct the formalism of the DW states on the basis
of the microscopic Fermi liquid theory AGD ; LW ; Baym-Kadanoff ; Baym .

In general, the DW state at wavevector ${{\bm{q}}}$ originates from the
particle-hole (p-h) pairing condensation,
$D_{{\bm{k}}}^{{{\bm{q}}}{\sigma}}=(1-P_{0})\langle
c_{{{\bm{k}}}+{{\bm{q}}},{\sigma}}^{\dagger}c_{{{\bm{k}}},{\sigma}}\rangle$,
where $c_{{{\bm{k}}},{\sigma}}$ is the electron annihilation operator,
${{\bm{k}}}$ is the momentum, and ${\sigma}\ (=\pm 1)$ is the spin index
Nersesyan ; Tazai-JPSJ-fRG ; Kontani-sLC . $P_{0}$ represents the projection
onto the totally-symmetric state with respect to the space-group and the time-
reversal Tazai-JPSJ-fRG ; Kontani-sLC . Rich classes of the DW states are
determined by the symmetry of the p-h condensation
$D_{{\bm{k}}}^{{{\bm{q}}}{\sigma}}$. For example, a simple spin-DW state is
given as the ${{\bm{k}}}$-independent function
$D_{{\bm{k}}}^{{{\bm{q}}}{\sigma}}=m{\sigma}$. The realized DW state
$D_{{\bm{k}}}^{{{\bm{q}}}{\sigma}}$ should be uniquely determined as the
extremum of the free energy.

In Fermi liquids, the single-electron kinetic term between sites $i$ and $j$
is given by $t_{i,j}^{{\sigma}}=t_{i-j}^{0}+\Sigma_{i,j}^{\sigma}$, where
$t_{i-j}^{0}$ is the hopping integral of the bare Hamiltonian, and
$\Sigma_{i,j}^{\sigma}$ is the self-energy due to the correlation between
other electrons. Here, we define the symmetry breaking part of the self-energy
Tazai-JPSJ-fRG ; Kontani-sLC :

$\displaystyle\delta t_{i,j}^{\sigma}\equiv(1-P_{0})\Sigma_{i,j}^{\sigma}.$
(1)

In the absence of the DW order ($T\geq T_{\rm c}$), we obtain $\delta
t_{i,j}^{\sigma}=0$ by definition. When the DW order emerges ($T\leq T_{\rm
c}$), $\delta t_{i,j}^{\sigma}$ becomes finite due to nonzero
$D_{{\bm{k}}}^{{{\bm{q}}}{\sigma}}$. Thus, $\delta t_{i,j}^{\sigma}$ is the
energy-dimension order parameter of the DW state.

Here, we consider the DW state at a constant wavevector ${{\bm{q}}}$. For
convenience, we introduce the “form factor $\delta
t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$” that is the Fourier transform of $\delta
t_{i,j}^{\sigma}$ Tazai-JPSJ-fRG ; Kontani-sLC :

$\displaystyle\delta
t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}\equiv\frac{1}{N}\sum_{ij}\delta
t_{ij}^{\sigma}e^{-i{{\bm{k}}}\cdot(\bm{r}_{i}-\bm{r}_{j})}e^{-i{{\bm{q}}}\cdot\bm{r}_{i}},$
(2)

where $\bm{r}_{i}$ is the position of site $i$. The classification of the
symmetry of the form factor is presented in Sect. I.1.

The theoretical way to derive the form factor $\delta
t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$ has not been established yet. The aim of
the present study is to establish an exact framework to derive $\delta
t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$, based on which we can construct reliable
and useful approximate theories. In the statistical mechanics, the symmetry
breaking with $\delta t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$ occurs when the
grand potential $\Omega$ is reduced. In other words, $\delta
t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$ is uniquely determined as the stationary
state with the minimum grand potential. In strongly correlated Fermi liquids,
a rigorous formalism of the grand potential $\Omega$ is given by the
Luttinger-Ward (LW) theory LW , In the LW theory, the LW function $\Phi[G]$,
which is the functional of electron Green function $G$, plays a central role.
The self-energy $\Sigma$ and the irreducible two-particle interaction $I$ are
uniquely derived from the functional derivatives of $\Phi[G]$ LW .

In this paper, we first introduced the “form factor $\delta
t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$” into the LW theory to analyze the
symmetry breaking at wavevector ${{\bm{q}}}$. Its
$({{\bm{k}}}-{{\bm{q}}}/2)$-dependence represents the symmetry of the DW. We
next derived the “DW equation” to obtain the form factor that minimize the LW
ground potential below $T_{c}$. In this theory, the optimized $\delta
t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$ is uniquely obtained because the DW
equation is equivalent to the thermodynamic stationary condition. This
formalism enables us to study rich variety of electron-correlation-driven DW
states ($\delta t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$) without assuming any
symmetry, like the analysis of the superconducting (SC) states
($\Delta_{{\bm{k}}}^{{\sigma}{\sigma}^{\prime}}$) based on the SC gap
equation. In addition, we derive an exact expression of the Ginzburg-Landau
(GL) free energy, $F\propto a_{{\bm{q}}}\phi^{2}$, where $\phi$ is the
amplitude of the DW order. The coefficient $a_{{\bm{q}}}$ is uniquely related
to the eigenvalue of the DW equation $\lambda_{{\bm{q}}}$. Thus, we can
calculate the thermodynamic properties and the stability of the DW state.

In the next stage, we apply the derived DW equation to FeSe by using the one-
loop LW function $\Phi_{\rm FLEX}[G]$ that represents the quantum interference
among paramagnons. This theory naturally explains the following essential
experimental reports: (i) Nematic Fermi surface (FS) and the Lifshitz
transition due to bond+orbital order FeSe-ARPES-Suzuki ; FeSe-Lif1 ; FeSe-Lif2
; Eremin . (ii) Curie-Weiss (CW) behavior of the nematic susceptibility
$\chi_{\rm nem}$ at higher temperatures, and the deviation from the CW
behavior at lower temperatures near the nematic quantum-critical-point (QCP).
Hosoi ; B2g-Ishida ; Terashima ; Fisher-Science2016 . (iii) Scaling relation
of the specific heat jump at $T_{\rm c}$, $\Delta C/T_{\rm c}\propto T_{\rm
c}^{b}$ ($b\sim 3$). This relation naturally explains the smallness of (or
undetected) $\Delta C/T_{c}$ reported in several nematic systems, such as
RbFe2As2 ($T_{c}\sim 40$K) Mizukami-B2g  and cuprate superconductors
($T_{c}\sim 200$K) Sato-CDW ; Murayama-CDW . In cuprates, the nematic
transition occurs at the pseudogap temperature $T^{*}$ Sato-CDW ; Murayama-CDW
, while no anomaly in specific heat at $T=T^{*}$ has been reported previously.

Interestingly, recent experiments have revealed that the nematic QCP is
clearly separated from the magnetic QCP in Fe(Se,S) and Fe(Se,Te) Shibauchi-
nemQCP  in addition to Na(Fe,Co)As Zheng-NaFeAs . Such clear separation
between two QCPs in addition to the key points (i)-(iii) are naturally
understood by the present theory. Therefore, the nematicity in FeSe is the
bond+orbital order due to the “paramagnon interference mechanism” Onari-SCVC ;
Yamakawa-PRX2016 ; Onari-FeSe ; Tazai-JPSJ-fRG .

Here, we construct the DW equation based on the LW theory. The LW theory is
expected to be valid for various strongly-correlated metals except for the
vicinity of the localized Mott states. Thus, the present theory is expected to
pave the path to understanding the useful concept of the “vestigial nematic
order” from the itinerant picture.

###  I.1 Form Factor

Here, we discuss the rich variety of unconventional DW states Schultz ;
Affleck ; Nersesyan ; Tazai-JPSJ-fRG ; Kontani-sLC  by classifying the
symmetry of the form factor. The exotic DW states are represented by the non-
local order parameter $\delta t_{ij}^{\sigma}$, which is parameterized by
different site indices $i$ and $j$. Then, the effective hopping integral is
$t_{i,j}^{\sigma}=t_{i-j}^{0}+\delta t_{ij}^{\sigma}$, where $t_{i-j}^{0}$ is
the original hopping integral with $A_{1g}$ symmetry. In conventional charge
(spin) orders, the order parameter is expressed as $\delta
t_{ij}^{\uparrow}=+(-)\delta t_{ij}^{\downarrow}$ with $i=j$, respectively.

In recent years, in contrast, unconventional non-local orders given by $\delta
t_{ij}^{\sigma}$ with $i\neq j$ have been discovered and attracted increasing
attentions. Here, we assume the Hermitian order parameter Tazai-JPSJ-fRG ;
Kontani-sLC :

$\displaystyle\delta t_{ij}^{\sigma}=(\delta t_{ji}^{\sigma})^{*},$ (3)
$\displaystyle\delta t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}=(\delta
t^{-{{\bm{q}}}{\sigma}}_{{{\bm{k}}}+{{\bm{q}}}})^{*},$ (4)

$\delta t_{ij}^{\sigma}$ is classified into four channels in terms of parity
symmetry (${\cal P}=\pm 1$) and time-reversal symmetry (${\cal T}=\pm 1$) as
shown in Fig.1 (a). Below, we discuss the case of ${{\bm{q}}}={\bm{0}}$.
First, we consider the case of $\delta t_{ij}^{\uparrow}=\delta
t_{ij}^{\downarrow}$. When $\delta t_{ij}$ is real, the bond order with
$({\cal P},{\cal T})=(+1,+1)$ is realized. As an example, the $d$-wave bond
order in square lattice model is shown in Fig. 1 (b). When $\delta t_{ij}$ is
pure imaginary, the charge-loop current (cLC) with $({\cal P},{\cal
T})=(-1,-1)$ is realized. Its form factor in ${{\bm{k}}}$-space is $\delta
t_{{{\bm{k}}}}\propto\cos k_{x}-\cos k_{y}$. An example of the cLC order in
anisotropic triangular lattice model is shown in Fig. 1 (c).

Next, we consider the case of $\delta t_{ij}^{\uparrow}=-\delta
t_{ij}^{\downarrow}$. When $\delta t_{ij}$ is real, the spin-bond order with
$({\cal P},{\cal T})=(+1,-1)$ appears. When $\delta t_{ij}$ is pure imaginary,
the spin-loop current (sLC) order with $({\cal P},{\cal T})=(-1,+1)$ appears.
An example of the sLC order is shown in Fig. 1 (d).

In addition, the translational symmetry is violated when the DW wavevector
${{\bm{q}}}$ is nonzero. Furthermore, orbital orders Onari-SCVC , valley
orders Onari-TBG , and multipole orders Tazai-CeB6  emerge when the Wanner
functions possess multiple degrees of freedom. These rich non-local DW states
are called the quantum-liquid-crystal (QLC) order Tazai-JPSJ-fRG , We note
that, in a simple single-site model, the cLC and sLC orders in real space
($\delta t_{ij}^{\sigma}=-\delta t_{ij}^{\sigma}$) are pure imaginary
according to the Hermitian condition. However, the Fourier transformation of
such current order $\delta t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$ becomes real.

Figure 1:  (a) Classification of nonlocal order parameters (quantum liquid
crystal states). Spin-symmetric (-antisymmetric) order corresponds to the
charge (spin) channel order. The four states are labeled as $({\cal P},{\cal
T})$, where ${\cal P}=\pm 1$ (${\cal T}=\pm 1$) is the parity of the
inversion-symmetry (time-reversal-symmetry). (b) $d$-wave bond order with
$({\cal P},{\cal T})=(+1,+1)$. (c) Charge-loop-current (cLC) order with
$({\cal P},{\cal T})=(-1,-1)$ in Ref. Tazai-cLC . (d) Spin-loop-current (sLC)
order with $({\cal P},{\cal T})=(-1,+1)$ in Ref. Kontani-sLC .

###  I.2 Stoner Theory for Ferro-Magnetic Transition

Here, we briefly review the mean-field theory of the ferro-magnetic (FM)
transition and derive a simplified DW equation and GL free energy. This
explanation will be useful to understand the derivation of the DW equation
based on the LW theory in later sections. It is also understood that the mean-
field theory is insufficient to explain the nematic order. Below, we consider
the following Hubbard model,

$\displaystyle H$ $\displaystyle=$ $\displaystyle\sum_{i\neq
j,\sigma}t_{i-j}^{0}c^{\dagger}_{i\sigma}c_{j\sigma}+H_{\rm{I}},$ (5)
$\displaystyle H_{\rm{I}}$ $\displaystyle=$
$\displaystyle\sum_{i\sigma}Uc^{\dagger}_{i\sigma}c_{i\sigma}c^{\dagger}_{i\bar{\sigma}}c_{i\bar{\sigma}}=\frac{U}{4}\sum_{i}(n_{i}^{2}-m_{i}^{2}),$
(6)

where $n_{i}=n_{i\uparrow}+n_{i\downarrow}$ and
$m_{i}=n_{i\uparrow}-n_{i\downarrow}$. The conduction electron energy is
${\epsilon}_{{\bm{k}}}=\frac{1}{N}\sum_{i,j}t_{i-j}^{0}e^{i(\bm{r}_{i}-\bm{r}_{j})\cdot{{\bm{k}}}}$.
In the mean-field approximation, the magnetic order is expressed as the spin-
dependent ${\delta}t$: $\delta t_{ii}^{\uparrow}=-\delta
t_{ii}^{\downarrow}=-\frac{U}{2}m_{i}$. Here, we set $M_{i}\equiv Um_{i}/2$.
In the case of the FM ($q=0$) order, $M_{i}=M$. Thus, the grand potential is
given by

$\displaystyle\Omega_{\rm{MF}}$ $\displaystyle=$
$\displaystyle\frac{T}{N}\sum_{k,\sigma}\ln
G_{k\sigma}^{\rm{MF}}+\frac{1}{U}M^{2}$ (7) $\displaystyle=$
$\displaystyle-\frac{1}{N}\sum_{{{\bm{k}}},\sigma}\ln(1+e^{-\beta(\epsilon_{{{\bm{k}}}{\sigma}}-\mu)})+\frac{1}{U}M^{2},$

where $k\equiv({{\bm{k}}},{\epsilon}_{n}=(2n+1)\pi T)$,
$\epsilon_{{{\bm{k}}}\sigma}=\epsilon_{{{\bm{k}}}}+M\sigma$, and
$G_{k\sigma}^{\rm{MF}}=(i\epsilon_{n}-\epsilon_{{{\bm{k}}}{\sigma}}+\mu)^{-1}$
is the electron Green function. From the stationary condition of the
$\Omega_{\rm{MF}}$, which is given by
${{\partial}\Omega_{\rm{MF}}}/{{\partial}M}=0$, the mean-field equation for
$M$ is obtained as

$\displaystyle
M=-\frac{U}{2N}\sum_{{{\bm{k}}}\sigma}f(\epsilon_{{{\bm{k}}}\sigma})\sigma.$
(8)

Next, we derive the linearized mean-field equation by linearizing the right-
hand side of Eq. (8). It is given as

$\displaystyle{\alpha}_{S}M=\frac{U}{N}\sum_{{{\bm{k}}}}\left(-\frac{{\partial}n(\epsilon_{{{\bm{k}}}})}{{\partial}\epsilon_{{{\bm{k}}}}}\right)M=U\chi^{0}(q=0)M,$
(9)

where $\alpha_{S}$ is the eigenvalue, which reaches unity at the critical
temperature. $n({\epsilon})$ is the Fermi distribution function. In the mean-
field approximation, $\alpha_{S}=U\chi^{0}(0)$, where $\chi^{0}(0)$ is the
irreducible susceptibility at ${{\bm{q}}}=0$. $\alpha_{S}$ is called the spin
Stoner factor. Equation (9) is the simplest example of the spin-channel DW
equation with $I^{s}_{kk^{\prime},q}=U$.

Here, we consider the GL free energy for the magnetic transition:

$\displaystyle F_{GL}=aM^{2}+\frac{1}{2}bM^{4}.$ (10)

Here, $F_{GL}=\Omega_{\rm MF}+\mu N$, where $\mu$ is the chemical potential.
Based on Eq. (7), the coefficient $a$ is expressed as

$\displaystyle a$ $\displaystyle=$
$\displaystyle\chi^{0}(0)\left(-1+\frac{1}{\alpha_{S}}\right).$ (11)

This equation is consistent with the Stoner theory stoner .

For the magnetic state at finite wavevector ${{\bm{q}}}$, the eigenvalue of
Eq. (9) is given as ${\alpha}_{S}({{\bm{q}}})=U\chi^{0}({{\bm{q}}},0)$, and
the ${{\bm{q}}}$-dependent coefficient of the GL free energy becomes
$\displaystyle
a_{{\bm{q}}}=\chi^{0}({{\bm{q}}},0)\left(-1+\frac{1}{\alpha_{S}({{\bm{q}}})}\right)$.
Note that $\displaystyle a_{{\bm{q}}}\approx
a_{{{\bm{q}}}={\bm{0}}}+\frac{1}{2}\sum_{\mu}^{x,y,z}c_{\mu}q_{\mu}^{2}$,
where
$c_{\mu}={\partial}_{\mu}^{2}a_{{\bm{q}}}/{\partial}q_{\mu}^{2}|_{{{\bm{q}}}={\bm{0}}}$.
Here,
$\displaystyle\chi^{0}({{\bm{q}}},0)=\frac{1}{N}\sum_{{\bm{k}}}\frac{n({\epsilon}_{{{\bm{k}}}+{{\bm{q}}}})-n({\epsilon}_{{{\bm{k}}}})}{{\epsilon}_{{{\bm{k}}}}-{\epsilon}_{{{\bm{k}}}+{{\bm{q}}}}}$
is the irreducible susceptibility for general ${{\bm{q}}}$. The relation
${\alpha}_{S}>{\alpha}_{S}({{\bm{q}}})$ holds for ${{\bm{q}}}\neq{\bm{0}}$ in
ferromagnetic metals.

Considering the $T$-dependence of $\chi^{0}(0)$, the coefficient
$a_{{{\bm{q}}}=\bm{0}}$ at $T=0$ is given as

$\displaystyle a_{{{\bm{q}}}=\bm{0}}(T=0)\simeq\frac{\pi^{2}}{3}BT_{\rm
c}^{2},$ (12)

where $B\equiv D^{\prime\prime}(0)+{(D^{\prime}(0))^{2}}/{D(0)}$, and $D(0)$
is the density of states at the Fermi energy. Thus, the necessary condition
for the FM transition ($T_{\rm c}>0$) is given as $B<0$ in the mean-field
approximation.

However, recent experiments have revealed that the ferro-DW order appears even
in the case of $B>0$ in several strongly correlated electron systems such as
Fe-based superconductors. Thus, to understand the ferro-DW transition, we have
to go beyond the mean field theory. In addition, exotic nonlocal DW orders
($\delta t_{ij}^{\sigma}$ with $i\neq j$) summarized in Fig. 1 cannot be
explained within the mean-field approximation. To solve these difficulties, in
this paper, we study the DW phase transition on the bases of the LW theory.

##  II Derivation of DW-equation from Luttinger-Ward theory

The Luttinger-Ward (LW) theory in Ref. LW  provides an exact expression of the
grand potential $\Omega$, which is applicable for strongly correlated metals
unless the perturbation treatment is violated. In the first half part of this
section, we discuss the order parameter at ${{\bm{q}}}={\bm{0}}$. Hereafter,
we omit the orbital indices of the Green functions and the Coulomb
interactions to simplify the expressions, because it is straightforward to
denote them explicitly. In the LW theory, the grand potential is expressed as

$\displaystyle\Omega_{\rm{LW}}[G]$ $\displaystyle=$
$\displaystyle\Omega_{\rm{F}}[G]+\Phi[G],$ (13)
$\displaystyle\Omega_{\rm{F}}[G]$ $\displaystyle=$
$\displaystyle\frac{T}{N}\sum_{k\sigma}\left\\{\ln(G_{k\sigma})\right.$ (14)
$\displaystyle\left.-G_{k\sigma}\left((G^{\rm{\rm{free}}}_{k\sigma})^{-1}-G^{-1}_{k\sigma}\right)\right\\},$

where $k\equiv({{\bm{k}}},{\epsilon}_{n})$: ${{\bm{k}}}=(k_{x},k_{y})$ is the
wavevector and ${\epsilon}_{n}=(2n+1)\pi T$ is the fermion Matsubara
frequency. Here, $T\sum_{k}\cdots\equiv
T\sum_{{\epsilon}_{n}}\sum_{{\bm{k}}}\cdots$, and $G_{k\sigma}$ is the Green
function with the self-energy:
$G_{k\sigma}=\left(\left(G_{k\sigma}^{\rm{free}}\right)^{-1}-\Sigma_{k\sigma}\right)^{-1}$.
Now, the self-energy can be divided into

$\displaystyle\Sigma=\Sigma^{0}+\delta t,$ (15)

where $\Sigma^{0}$ is the “normal state self-energy” without any symmetry
breaking, and $\delta t$ is equal to the DW order parameter introduced in
Sec.I. Here, $\Sigma^{0}$ belongs to $A_{1g}$ symmetry, while ${\delta}\Sigma$
belongs to non-$A_{1g}$ symmetry. Thus, $\delta t=0$ for $T>T_{\rm c}$.

In Eq. (14), $\Phi[G]$ is the Luttinger-Ward function which is given by
calculating the all closed linked skeleton diagrams. Based on Eq. (14), we can
define $\Omega$ as a functional of $\Sigma$ Potthoff .

$\displaystyle\Omega[\Sigma]=-\frac{T}{N}\sum_{k\sigma}\ln(-(G_{k\sigma}^{\rm{\rm{free}}})^{-1}+\Sigma_{k\sigma})+P[\Sigma],$
(16)

where $P[\Sigma]$ is considered as the Legendre-transformation of $\Phi[G]$
introduced by Potthoff Potthoff .

$\displaystyle P[\Sigma]$ $\displaystyle\equiv$
$\displaystyle-\frac{T}{N}\sum_{k\sigma}G_{k\sigma}\Sigma_{k\sigma}+\Phi[G].$
(17)

In deriving the GL free energy, we have to analyze $\Omega[\Sigma]$ in Eq.
(16). Using the Luttinger-Ward function $\Phi$ and Potthoff function $F$
Potthoff , the self-energy and Green function are respectively expressed as

$\displaystyle\frac{{\delta}\Phi[G]}{{\delta}G_{k\sigma}}$ $\displaystyle=$
$\displaystyle\Sigma_{k\sigma}[G],$ (18)
$\displaystyle\frac{{\delta}P[\Sigma]}{{\delta}\Sigma_{k\sigma}}$
$\displaystyle=$ $\displaystyle-G_{k\sigma}[\Sigma].$ (19)

Then, the functional derivations of $\Omega[G]$ and $\Omega[\Sigma]$ are
respectively given by

$\displaystyle\frac{{\delta}\Omega[G]}{{\delta}G_{k\sigma}}$ $\displaystyle=$
$\displaystyle
G_{k\sigma}^{-1}-(G^{\rm{free}}_{k\sigma})^{-1}+\Sigma_{k\sigma}[G],$ (20)
$\displaystyle\frac{{\delta}\Omega[\Sigma]}{{\delta}\Sigma_{k\sigma}}$
$\displaystyle=$
$\displaystyle\frac{1}{(G_{k\sigma}^{\rm{free}})^{-1}-\Sigma_{k\sigma}}-G_{k\sigma}[\Sigma].$
(21)

When $\Omega$ is stationary, the following Dyson equation is satisfied,

$\displaystyle\Sigma_{k\sigma}[G]=(G_{k\sigma}^{\rm{free}})^{-1}-G_{k\sigma}^{-1}\,\,\,\,\,\,\,\,(\mbox{from
Eq.(\ref{dyson1})}=0),$ (22) $\displaystyle
G_{k\sigma}[\Sigma]=\frac{1}{(G_{k\sigma}^{\rm{free}})^{-1}-\Sigma_{k\sigma}}\,\,\,\,\,\,\,\,(\mbox{from
Eq.(\ref{dyson2})}=0).$ (23)

Based on the Luttinger-Ward theory, the ferro $(q=0)$ DW transitions are
naturally described by the self-consistent equation (We call this the DW
equation).

Here, we introduce the irreducible 4-point vertex
$I^{\sigma\sigma^{\prime}}_{kk^{\prime}}$ shown in Fig.2 (a). It is a Jacobian
connecting $\Sigma$ and $G$ as Potthoff

$\displaystyle\frac{{\delta}\Sigma_{k\sigma}[G]}{{\delta}G_{k^{\prime}\sigma^{\prime}}}$
$\displaystyle=$ $\displaystyle I^{\sigma\sigma^{\prime}}_{kk^{\prime}},$ (24)
$\displaystyle\frac{{\delta}G_{k\sigma}[\Sigma]}{{\delta}\Sigma_{k^{\prime}\sigma^{\prime}}}$
$\displaystyle=$
$\displaystyle\left\\{I^{\sigma\sigma^{\prime}}_{kk^{\prime}}\right\\}^{-1}.$
(25)

Therefore, the following analytical relation Potthoff  is obtained:

$\displaystyle\frac{T}{N}\sum_{k^{\prime\prime}\sigma^{\prime\prime}}\frac{{\delta}\Sigma_{k\sigma}[G]}{{\delta}G_{k^{\prime\prime}\sigma^{\prime\prime}}}\,\,\frac{{\delta}G_{k^{\prime\prime}\sigma^{\prime\prime}}[\Sigma]}{{\delta}\Sigma_{k^{\prime}\sigma^{\prime}}}=\delta_{kk^{\prime}}\delta_{\sigma\sigma^{\prime}},$
(26)

which is exactly satisfied when $G$ is uniquely determined from $\Sigma$ via
one-to-one correspondence.

Figure 2:  (a) Definition of the irreducible 4-point vertex function
$I^{x}_{kkq}$ ($x=s,c$). (b) Right-hand side of the linearized DW equation
composed of $I^{x}_{kkq}$. $f_{k}^{q,x}$ is the form factor at wavevector $q$.

The DW equation is derived from the following stationary conditions

$\displaystyle\left.\frac{{\delta}\Omega_{\rm{LW}}[G]}{{\delta}\Sigma_{k\sigma}}\right|_{\Sigma^{0}}$
$\displaystyle=$ $\displaystyle 0,\,\,\,\,\,\,(\mbox{at any $T$})$ (27)
$\displaystyle\left.\frac{{\delta}\Omega_{\rm{LW}}[G]}{{\delta}\Sigma_{k\sigma}}\right|_{\bar{\Sigma}}$
$\displaystyle=$ $\displaystyle 0\,\,\,\,\,\,\,\,\,(T<T_{\rm c}),$ (28)

where $\Sigma^{0}$ is the self-energy without any symmetry breaking, and
$\bar{\Sigma}$ is the stationary self-energy after the symmetry breaking. For
$T>T_{\rm c}$, the thermodynamic state is given by Eq. (27), which corresponds
to the minimum of the free energy shown in Fig. 3 (a). For $T<T_{\rm c}$, Eq.
(28) gives the symmetry breaking state shown in Fig. 3 (b). (Eq. (28)
corresponds to the unstable extremum.)

Figure 3:  Schematic picture of the GL free energy (a) above $T_{\rm c}$
($a(T)>0$) and (b) below $T_{\rm c}$ ($a(T)<0$). Stationary points are shown
by blue circles.

By using Eq.(28) and (21),

$\displaystyle\left.\frac{{\delta}\Omega_{\rm{LW}}[G]}{{\delta}\Sigma_{k\sigma}}\right|_{\bar{\Sigma}}$
$\displaystyle=$ $\displaystyle-
G_{k\sigma}^{2}\Sigma_{k\sigma}+\frac{{\delta}\Phi[G]}{{\delta}\Sigma_{k\sigma}}$
(29) $\displaystyle=$ $\displaystyle-
G_{k\sigma}^{2}\Sigma_{k\sigma}+\frac{T}{N}\sum_{k^{\prime}\sigma^{\prime}}\frac{{\delta}G_{k^{\prime}\sigma^{\prime}}}{{\delta}\Sigma_{k\sigma}}\frac{{\delta}\Phi[G]}{{\delta}G_{k^{\prime}\sigma^{\prime}}},$

where
$\frac{{\delta}G_{k^{\prime}\sigma^{\prime}}}{{\delta}\Sigma_{k\sigma}}=G^{2}_{k\sigma}\delta_{kk^{\prime}}$.
Thus, the stationary condition of Eq. (28) is rewritten as

$\displaystyle\left.\frac{{\delta}\Phi[G]}{{\delta}G_{k\sigma}}\right|_{G^{0}}$
$\displaystyle=$ $\displaystyle\Sigma_{k\sigma}^{0},$ (30)
$\displaystyle\left.\frac{{\delta}\Phi[G]}{{\delta}G_{k\sigma}}\right|_{\bar{G}}$
$\displaystyle=$ $\displaystyle\bar{\Sigma}_{k\sigma},$ (31)

where ${\hat{G}}^{0}=(\\{{\hat{G}}^{\rm
free}\\}^{-1}-{\hat{\Sigma}}^{0})^{-1}$. Equations (30) and (31) compose the
“exact DW equation” that describes the DW state below $T_{\rm c}$. Then, the
order parameter is $\delta
t_{k\sigma}=\bar{\Sigma}_{k\sigma}-\Sigma_{k\sigma}^{0}$.

It is noteworthy that superconducting gap equation is derived from the
stationary condition of the LW function in the SC state $\delta\Phi[G,F]$,
where $F$ ($G$) is the anomalous (normal) Green function Yanase . Thus, the
derived DW equation for the form factor, Eqs. (30) and (31), is well-founded
comparable with the well-known superconducting-gap equation.

Next, we derive the linearized DW equation with respect to $\delta t$. By
subtracting Eq. (30) from (31) , we obtain

$\displaystyle{\delta}t_{k\sigma}$ $\displaystyle=$
$\displaystyle\frac{T}{N}\sum_{k^{\prime}\sigma^{\prime}}\left.\frac{{\delta}^{2}\Phi[G]}{{\delta}G_{k^{\prime}\sigma^{\prime}}{\delta}G_{k\sigma}}\right|_{\bar{G}}{\delta}G_{k^{\prime}\sigma^{\prime}}$
(32) $\displaystyle=$
$\displaystyle\frac{T}{N}\sum_{k^{\prime}\sigma^{\prime}}\left.\frac{{\delta}^{2}\Phi[G]}{{\delta}G_{k^{\prime}\sigma^{\prime}}{\delta}G_{k\sigma}}\right|_{G^{0}}{\delta}G_{k^{\prime}\sigma^{\prime}}+O({\delta}G^{2}),$

where ${\delta}G\equiv\bar{G}-G^{0}$. Here, we rewrite ${\delta}t_{k\sigma}$
as

$\displaystyle{\delta}t_{k\sigma}\equiv\phi\cdot f_{k\sigma},$ (33)

where $\phi$ is a real parameter, and $f_{k}^{q}$ is the normalized order
parameter that belongs to one of the irreducible representations in
non-$A_{1g}$ symmetry. It is convenient to set $\max_{k}|f_{k\sigma}|=1$
because the relation $\phi=\max_{k}|{\delta}t_{k\sigma}|$ holds. Thus, we
derive the following “linearize DW equation for ${{\bm{q}}}={\bm{0}}$” by
introducing the eigenvalue $\lambda$ to the left-hand side of Eq. (32);

$\displaystyle\lambda
f_{k\sigma}=\frac{T}{N}\sum_{k^{\prime}\sigma^{\prime}}I^{\sigma\sigma^{\prime}}_{kk^{\prime}}(G^{0}_{k^{\prime}\sigma^{\prime}})^{2}f_{k^{\prime}\sigma^{\prime}},$
(34)

where we denote
$I^{\sigma\sigma^{\prime}}_{kk^{\prime}}\equiv\left.I^{\sigma\sigma^{\prime}}_{kk^{\prime}}\right|_{\Sigma^{0}}$
to simplify the notation. In Eq. (34), the largest eigenvalue $\lambda$
reaches unity at $T=T_{\rm c}$, and its eigenvector gives the form factor of
the DW state.

The linearized DW equation can be generalized for finite $q$ orders as
follows. First, we consider the DW order with the wavevector
$\bm{q}=\bm{g}\cdot m/n$ where $\bm{g}$ is the reciprocal lattice vector, and
$m$, $n$ are integers ($0\leq m<n$). Then, we can introduce the $n\times
n$-matrix Green function $G^{lm}_{k}=\langle k+lq|\hat{G}|k+mq\rangle$ where
$l,m=0\sim n-1$. In this case, Eq. (31) becomes

$\displaystyle\left.\frac{{\delta}\Phi[G]}{{\delta}G_{k\sigma}}\right|_{\bar{G}^{lm}}$
$\displaystyle=$ $\displaystyle\bar{\Sigma}_{k\sigma}^{lm}.$ (35)

Hereafter, we drop the overlines of $\bar{G}$ and $\bar{\Sigma}$ to simplify
the notation.

Here, we adopt the extended Brillouin zone scheme for $\vec{k}$ to simplify
the explanation. After that, the DW equation can be linearized with respect to
the ${{\bm{q}}}$-linear term in ${\delta}t$ given by
${\delta}t^{m+1,m}\equiv{\delta}t^{q}_{k}$. By introducing the $q$-dependent
eigenvalue, we obtain the following “linearize DW equation for general
${{\bm{q}}}$”,

$\displaystyle\lambda_{q}f_{k\sigma}^{q}=\frac{T}{N}\sum_{k^{\prime}\sigma^{\prime}}I^{\sigma\sigma^{\prime}}_{kk^{\prime}q}G^{0}_{k^{\prime}\sigma^{\prime}}G^{0}_{k^{\prime}+q\sigma^{\prime}}f_{k^{\prime}\sigma^{\prime}}^{q},$
(36)

where $q\equiv({{\bm{q}}},{\omega}_{l})$: ${{\bm{q}}}=(q_{x},q_{y})$ is the
wavevector and ${\omega}_{l}=2l\pi T$ is the boson Matsubara frequency. The
condition $\lambda_{q}=1$ brings the DW transition temperature $T_{\rm c}$
with wavevector $q$, which can be interpreted as the particle-hole (p-h) gap
equation. The form factor of the eigenvalue equation (36) contains the
uncertainty of the phase factor $e^{i\theta}$. The correct phase $\theta$ is
uniquely fixed by following the Hermitian condition in Eq. (4).

The DW equation (36) is further simplified by introducing the spin (s) and
charge (c) channel functions in the absence of the spin-orbit interaction
Kontani-sLC :

$\displaystyle f^{c(s)}$ $\displaystyle=$ $\displaystyle
f_{\uparrow}+(-)f_{\downarrow}.$ (37) $\displaystyle I^{c(s)}$
$\displaystyle=$ $\displaystyle
I^{\uparrow\uparrow}+(-)I^{\uparrow\downarrow}.$ (38)

Finally, we derive the simplified “linearized DW equation for $x\ (=s,c)$
channel form factor at $q$” from Eq. (36) as follows:

$\displaystyle\lambda^{x}_{q}\,\,f_{k}^{q,x}=\frac{T}{N}\sum_{k^{\prime}}I^{x}_{kk^{\prime}q}G^{0}_{k^{\prime}}G^{0}_{k^{\prime}+q}\,\,f_{k^{\prime}}^{q,x}.$
(39)

The right-hand side of Eq. (39) is shown in Fig.2 (b). Therefore, we derived
the exact expression of the DW equation in Eq. (36) or Eq. (39) composed of
the true irreducible four-point vertex $I$ and the self-energy.

We can derive the expression of $\lambda^{x}_{q}$ from Eq. (39) as

$\displaystyle\lambda^{x}_{q}=X^{x}_{q}/\chi^{0f}(q),$ (40)

where $X$ is given as

$\displaystyle
X^{x}_{q}=\frac{T^{2}}{N^{2}}\sum_{k,k^{\prime}}(f_{k}^{q,x})^{*}G^{0}_{k}G^{0}_{k+q}I^{x}_{kk^{\prime}q}G^{0}_{k^{\prime}}G^{0}_{k^{\prime}+q}f_{k^{\prime}}^{q,x},$
(41)

and $\chi^{0f}(q)$ is the irreducible susceptibility with the form factor
$f_{k\sigma}^{q}$:

$\displaystyle\chi^{0f}(q)$ $\displaystyle\equiv$
$\displaystyle-\frac{T}{N}\sum_{k\sigma}f_{k\sigma}^{q}G^{0}_{k+q\sigma}G^{0}_{k\sigma}f_{k\sigma}^{-q}.$
(42)

From Eq. (40), the relation $\lambda^{x}_{q}\propto X^{x}_{q}$ is obtained
because $\chi^{0f}(q)$ is nearly $T$-independent.

It is noteworthy that the DW equation introduced by Onari and Kontani et al.
in Refs. Onari-FeSe ; Kawaguchi-CDW ; Tsuchiizu-CDW ; Tazai-JPSJ-fRG , which
has been applied to iron-based and cuprate superconductors, is derived from
the exact DW equation given in Eq. (36). The detailed derivation is shown in
Appendix A.

Before closing this section, we reproduce the Stoner theory by applying the
mean-field approximation to Eq. (39). In the mean-field theory,

$\displaystyle
I^{s}_{kk^{\prime}q}=-U,\,\,\,\,\,\,\,\,\,\,\,\,{\delta}\Sigma^{q,s}_{k}=M\delta_{q,0}.$
(43)

Therefore, the Eq. (39) is given by

$\displaystyle\lambda^{s}M=-\frac{T}{N}\sum_{k^{\prime}}G^{0}_{k^{\prime}}G^{0}_{k^{\prime}}UM=U\chi^{0}M.$
(44)

Thus, the eigenvalue for the ferromagnetic transition corresponds to the
Stoner factor $\alpha_{S}$:

$\displaystyle\lambda^{s}=U\chi^{0}\equiv\alpha_{S}.$ (45)

##  III Derivation of GL free energy based on the linearized DW equation

###  III.1 GL free energy for DW state with form factor $f$

Here, we derive the expression of the Ginzburg-Landau (GL) free energy based
on the DW-equation. The GL free energy due to the ferro-DW transition
(${\delta}t_{k\sigma}=\phi f_{k\sigma}$) is expressed as

$\displaystyle\Omega_{\rm{DW}}(T,\mu,\phi)=a(T)\phi^{2}+\frac{1}{2}b\phi^{4}.$
(46)

Its schematic picture above $T_{\rm c}$ ($a(T)>0$) and that below $T_{\rm c}$
($a(T)<0$) are shown in Figs. 3 (a) and (b), respectively.

Note that the coefficients $a$ and $b$ are functional of the form factor
$f_{k\sigma}$. By using the $\Omega[\Sigma]$ defined in Eq.(16), the
coefficient $a$ is calculated from the second functional derivation of
$\Omega[\Sigma]$ as

$\displaystyle\Omega[\Sigma^{0}+{\delta}t]-\Omega[\Sigma^{0}]$ $\displaystyle\
\ \ \ \ \approx\frac{T}{N}\sum_{k\sigma
k^{\prime}\sigma^{\prime}}\left.\frac{{\delta}^{2}\Omega[\Sigma]}{{\delta}\Sigma_{k^{\prime}\sigma^{\prime}}{\delta}\Sigma_{k\sigma}}\right|_{\Sigma=\Sigma^{0}}{\delta}t_{k^{\prime}\sigma^{\prime}}{\delta}t_{k\sigma}.$
(47)

By using Eq. (21), we obtain

$\displaystyle\frac{{\delta}^{2}\Omega[\Sigma]}{{\delta}\Sigma_{k^{\prime}\sigma^{\prime}}{\delta}\Sigma_{k\sigma}}$
$\displaystyle=$
$\displaystyle\frac{{\delta}}{{\delta}\Sigma_{k^{\prime}\sigma^{\prime}}}\left\\{\frac{1}{(G_{k\sigma}^{\rm{free}})^{-1}-\Sigma_{k\sigma}}-G_{k\sigma}[\Sigma]\right\\}$
(48) $\displaystyle=$ $\displaystyle
G_{k\sigma}^{2}\delta_{kk^{\prime}}-\frac{{\delta}G_{k\sigma}[\Sigma]}{{\delta}\Sigma_{k^{\prime}\sigma^{\prime}}}$
$\displaystyle=$
$\displaystyle(G^{0}_{k\sigma})^{2}\delta_{kk^{\prime}}-\left\\{I_{kk^{\prime}}^{\sigma\sigma^{\prime}}\right\\}^{-1},$

where we used the relation in Eq. (25). Therefore, the Eq. (47) is rewritten
as

$\displaystyle\Omega[\Sigma^{0}+{\delta}t]-\Omega[\Sigma^{0}]\approx\frac{T}{N}\sum_{k\sigma}(G^{0}_{k\sigma})^{2}({\delta}t_{k\sigma})^{2}$
$\displaystyle\ \ \ \ \ \ -\frac{T}{N}\sum_{k\sigma
k^{\prime}\sigma^{\prime}}\left\\{I_{kk^{\prime}}^{\sigma\sigma^{\prime}}\right\\}^{-1}{\delta}t_{k\sigma}{\delta}t_{k^{\prime}\sigma^{\prime}}.$
(49)

Here, we recall that the order parameter ${\delta}t_{k\sigma}=\phi\cdot
f_{k{\sigma}}$ is determined by using the DW equation (36). By using Eq. (36)
together with Eqs. (24)-(26), Eq. (49) is rewritten as

$\displaystyle\Omega[\Sigma^{0}+{\delta}t]-\Omega[\Sigma^{0}]$
$\displaystyle\simeq$
$\displaystyle\left(1-\frac{1}{\lambda}\right)\frac{T}{N}\sum_{k\sigma}(G^{0}_{k\sigma})^{2}f_{k\sigma}^{2}\phi^{2}$
(50) $\displaystyle=$
$\displaystyle-2\chi^{0f}(0)\left\\{1-\frac{1}{\lambda}\right\\}\phi^{2},$

where the factor $2$ originates from the spin degeneracy.

We can derive the GL free energy for the order parameter at nonzero wavevector
$q$ by considering the large unit cell as we discussed in Sect. II. Thus, the
coefficient $a$ defined in the Gibbs free energy in Eq. (46) is obtained as

$\displaystyle a_{q}(T)$ $\displaystyle=$
$\displaystyle-2\chi^{0f}(q)\left(1-\frac{1}{\lambda_{q}}\right).$ (51)

As a result, we obtain the exact expression for coefficient $a_{q}$ by using
the eigenvalue $\lambda_{q}$ in the DW equation. The obtained general
expression in Eqs. (42) and (51) are meaningful to discuss the DW transition.

Finally, we stress that the Potthoff’s Legendre-transformation of the LW
formalism Potthoff  is necessary to derive the correct GL free energy
expression. In Appendix B, we explain that the expansion of (14) with respect
to $\delta t$ leads to an inaccurate expression.

###  III.2 GL free energy for BCS Superconductivity

Here, we consider the GL equation for the spin-singlet superconductivity.
Here, we express the spin-singlet SC gap function as $\Delta_{k}=\psi\cdot
f_{k}$, where $f_{{\bm{k}}}$ is the normalized form factor. Based on the LW
theory for the superconducting states, we can derive that the second order GL
parameter is given by

$\displaystyle a(T)$ $\displaystyle=$ $\displaystyle-2\chi_{\rm
pp}^{0\psi}(0)\left(1-\frac{1}{\lambda_{\rm sc}}\right),$ (52)
$\displaystyle\chi_{\rm pp}^{0\psi}(0)$ $\displaystyle=$
$\displaystyle\frac{T}{N}\sum_{k}|\psi_{k}|^{2}G^{0}_{k}G^{0}_{-k},$ (53)

where $\lambda$ is the eigenvalue of the linearized gap equation given by

$\displaystyle\lambda_{\rm
sc}f_{k}=\frac{T}{N}\sum_{k^{\prime}}V_{kk^{\prime}}G^{0}_{k^{\prime}}G^{0}_{-k^{\prime}}f_{k^{\prime}},$
(54)

where
$V_{kk^{\prime}}=\left.\frac{{\delta}^{2}\Phi}{{\delta}F{\delta}F^{\dagger}}\right|_{\Delta=0}$.
Here, $F$ and $F^{\dagger}$ are anomalous Green functions. The derivation of
Eq. (52) is essentially the same as that for the DW transition given in
previous sections. We can show that the relationship (52) is also valid for
the spin-triplet superconductivity.

##  IV numerical analysis of nematic state in FeSe

In this section, we explain the important unsolved problems in FeSe, which is
one of the most famous Fe-based superconductors. We try to understand the
following key topic on the nematicity, for both above and below the nematic
transition temperature $T_{\rm c}$, based on a unified theory. (i) Lifshitz
transition below $T_{\rm c}$. (ii) Nematic susceptibility above and below
$T_{\rm c}$. (iii) Specific heat jump at $T=T_{\rm c}$. In FeSe, $T_{\rm c}$
corresponds to the structure transition temperature $T_{S}$.

In previous sections, we derived the exact expressions of the linearized DW
equation in Eq. (39) for $T>T_{\rm c}$ and the full DW equation in Eqs.
(30)-(31) for $T<T_{\rm c}$. Here, we solve these equations for FeSe based on
the one-loop approximation for the LW function, $\Phi_{\rm FLEX}$, derived in
Appendix A. We include the normal state (=without order parameter) self-energy
$\Sigma^{0}$ into the DW equations because it is necessary to satisfy the
stationary condition (30)-(31), although it was dropped in previous studies
Onari-FeSe .

We study a realistic $d+p$ orbital Hubbard model with on-site multiorbital
Coulomb interaction $U$ for FeSe:

$\displaystyle H=H_{0}+rH_{U},$ (55)

where $H_{0}$ is the $d$+$p$ orbital tight-binding model for FeSe, and $H_{U}$
is the $d$-orbital Coulomb interaction for $d$+$p$ orbital model given by the
constrained-random-phase-approximation (cRPA) method. The matrix elements in
$H_{U}$ are composed of the intra-orbital Coulomb repulsion $U_{l,l}$, the
inter-orbital Coulomb repulsion $U_{l,m}$ $(l\neq m)$, and the exchange term
$J_{l,m}$, as we explain in Appendix C.

In Eq. (55), $r(<1)$ is the reduction factor of $H_{U}$ that represents the
screening due to $p$ orbitals. According to Ref. Miyake , the averaged intra-
orbital $U_{\rm av}\sim 7eV$ is reduced to $\sim 4$eV due to the screening
effect by $p$-orbitals. In the present study, we set $r=0.3\sim 0.4$, where
$T_{c}$ increases with $r$. In contrast, $T_{c}$ slowly decreases with $r$ for
$r>0.5$. Thus, obtained $T_{c}$ depends on $r$, while the symmetry and the
form factor of the nematic order is insensitive to the choice of $r$. We note
that the relation $\alpha_{S}<1$ is satisfied for any $r$ in the present two-
dimensional FeSe model because the FLEX approximation satisfies the Mermin-
Wagner theorem Mermin . This fact is favorable for realizing the nonmagnetic
nematic state ($\lambda>1$ and $\alpha_{S}<1$).

In the present numerical study for FeSe, we use $64\times 64$
${{\bm{k}}}$-meshes, and 4096 or 8192 Matsubara frequencies. Figure 4 (a)
represents the Fermi surface (FS) of FeSe model. We derive the normal self-
energy $\Sigma^{0}$ by applying the FLEX approximation. In the case of
$r=0.36$, the obtained orbital-dependent mass-enhancement factors at $T=10$meV
are about $z_{xy}^{-1}\sim 5$ and $z_{xz,yz}^{-1}\sim 3$, respectively. The
Stoner factor is about $0.9$ and its $T$-dependence is very weak.

###  IV.1 Above $T_{\rm c}$: linearized DW equation analysis

First, we analyze the multiorbital Hubbard model for FeSe based on the
linearized DW equation in Eq. (39), with the kernel function in Eq. (68).
Here, we incorporate the normal state self-energy $\Sigma^{0}$ given in Eq.
(67) into the DW equation to perform the conserving approximation.
$\Sigma^{0}$ is significant to derive realistic $T_{\rm c}$ and beautiful
CW/non-CW behaviors of $\chi_{\rm nem}$, although it has been dropped in our
previous analyses.

Here, we discuss the kernel function in Eq. (68). The first line in Eq. (68)
gives the Hartree term, Maki-Thompson (MT) term, and the second and the third
lines in Eq. (68) give the Aslamazov-Larkin (AL) terms. Both MT and AL terms
cause important “fluctuation-induced interaction for the DW”. In Fe-based SCs,
the nematic order mainly originates from the AL terms, which represent the
“interference between paramagnons” Onari-SCVC ; Onari-FeSe ; Yamakawa-PRX2016
; Chubukov-AL . The MT term is also important to induce the characteristic
sign-reversing in the form factor in ${{\bm{k}}}$-space Onari-FeSe ; Chubukov-
nematic-rev . On the other hand, the cLC orders in geometrically frustrated
Hubbard models mainly originate from the the MT terms Tazai-cLC . Note that
the MT terms induce striking non-Fermi-liquid transport phenomena near the
QCPs Kontani-ROP .

Figure 4 (b) shows the ${{\bm{q}}}$-dependence of the largest charge-channel
eigenvalue $\lambda_{{\bm{q}}}^{c}$ at $r=0.40$ and $T=5$meV derived from the
DW equation. (Below, we drop the superscript $c$ of $\lambda_{{\bm{q}}}^{c}$
for simplicity.) The obtained $\lambda_{{\bm{q}}}$ exhibit the maximum at
${{\bm{q}}}={\bm{0}}$ because the convolution of two $\chi_{q}^{s}$’s,
$C_{q}\equiv\sum_{p}\chi^{s}_{p}\chi^{s}_{p+q}$, included in the AL-type VCs
include is largest at $q=0$. Here, $\lambda_{{{\bm{q}}}={\bm{0}}}$ exceeds
unity, and typical transition temperature $T_{\rm c}\ (\sim 100{\rm K})$ in
Fe-based SCs is reproduced by including $\Sigma^{0}$. We note that similar
results were obtained in Fig. S4 (a) in Ref. Onari-smectic  by taking
$\Sigma^{0}$ into account.

The obtained form factors belong to $B_{1g}$ symmetry, which is shown in Fig.4
(c): The $xz,yz$-orbital form factors express the ${{\bm{k}}}$-dependent
orbital polarization, which has been reported by previous DW equation studies
without $\Sigma^{0}$ Onari-FeSe . The obtained $xz,yz$-orbital polarization
elongates the hole-pocket along the $k_{y}$-axis as experimentally reported in
Refs. FeSe-ARPES-Suzuki ; FeSe-ARPES-Zhang . In addition, the $xy$-orbital
form factor the represents the $d$-orbital bond order $f\propto\cos k_{x}-\cos
k_{y}$ emerges at the same $T_{\rm c}$. This $d$-wave order leads to the
disappearance of an electron-pocket around Y-point FeSe-Lif1 ; FeSe-Lif2 ;
Eremin . Thus, experimentally-observed ferro-nematic order in FeSe is
naturally obtained. We stress that the coexistence of the $xz,yz$-orbital
order and the $d_{x^{2}-y^{2}}$-wave bond order on $xy$-orbital was already
reported in Fig. S3 (a)-(d) in Ref. Onari-smectic . In addition, the
$d_{xy}$-wave bond order on $xy$-orbital has been studied in RbFe2As2 in Ref.
Onari-B2g .

We comment that a simple $A_{1g}$ symmetry order that accompanies the net
charge order is suppressed by the Hartree term in the kernel function. In
contrast, the $B_{1g}$ symmetry order in Fig. 4 (c) is free from the
suppression by the Hartree term due to sign reversal in the form factor.

Figure 4:  (a) FS of FeSe model. Here, the weights of the $xz$, $yz$, and $xy$
orbitals are shown by green, red, and blue colors. (b) Obtained eigenvalue
$\lambda_{{\bm{q}}}$ at $T=5$meV and $r=0.40$ derived from the linearized DW
equation. The fact that $\lambda_{{\bm{q}}}$ exhibits the maximum at
${{\bm{q}}}={\bm{0}}$ means the emergence of the ferro-DW order. (c)
Renormalized form factor at ${{\bm{q}}}={\rm 0}$, $\delta
t_{{{\bm{k}}},m,m}^{*}$ with $m=m^{\prime}=xz$, $yz$ and $xy$, derived from
the full DW equation at $T=5$meV; see Fig. 5. The deformed FS in the nematic
state is shown. (d) $T$-dependence of $\lambda_{{{\bm{q}}}={\bm{0}}}(T)$ for
$r=0.32$, $0.36$ and $0.40$. (e) Obtained $R=(-\dot{\lambda}T_{\rm c})$ as a
function of $T_{\rm c}$ for $r=0.40\sim 0.32$. The relation $R\propto T_{\rm
c}^{b}$ with $b\approx 3$ holds. The inset shows $R$ in the RbFe2As2 model.
Note that $R\approx 1$ ($b=0$) in the BCS superconductivity.

Figure 4 (d) shows the temperature-dependence of
$\lambda_{{{\bm{q}}}={\bm{0}}}(T)$ for $r=0.40\sim 0.32$. At higher
temperatures ($T\gtrsim 10$meV), $\lambda_{{{\bm{q}}}={\bm{0}}}(T)$ exhibits
the $T$-linear behavior. The nematic susceptibility above $T_{\rm c}$ is
$\chi_{\rm nem}=\chi^{0f}(0)(1-\lambda_{{{\bm{q}}}={\bm{0}}})^{-1}$, as proved
theoretically in Ref. Onari-B2g . Therefore, it is confirmed that the
experimental CW behavior of $\chi_{\rm nem}$ at higher temperature is
naturally explained in the present theory. The deviation from the CW behavior
at lower temperatures will be discussed in Sect. V.2.

Figure 4 (e) shows the obtained $R=(-\dot{\lambda}T_{\rm c})$ as a function of
$T_{\rm c}$. Here,
$\displaystyle\dot{\lambda}\equiv\left.\frac{d\lambda}{dT}\right|_{T_{\rm
c}}$, $R$ approximately corresponds to $\lambda_{{{\bm{q}}}={\bm{0}}}(T=0)-1$,
and $R$ becomes very small when $T_{\rm c}\ll 10$meV. The reason is the
recovery of the Fermi-liquid behavior $\lambda(0)-\lambda(T)\propto T^{2}$
because the system is far from the magnetic QCP. Also, the relation $R\propto
T_{\rm c}^{b}$ ($b\sim 3$) is satisfied in the present study. In Sect. V.1, we
will explain that $R$ is proportional to the jump in the heat capacity at
$T_{\rm c}$.

###  IV.2 Below $T_{\rm c}$: full DW equation analysis

Next, we analyze the full DW equation given as Eqs. (30)-(31) for $T<T_{\rm
c}$. We safely assume the uniform (${{\bm{q}}}={\bm{0}}$) order parameter
because $\lambda_{{\bm{q}}}$ takes the largest value at ${{\bm{q}}}={\bm{0}}$
as found in Fig. 4 (b). The aim of this subsection is to explain the essential
properties of the nematic state ($T<T_{\rm c}$) in FeSe based on the
paramagnon interference mechanism.

Now, we explain the procedure of the numerical study in detail: The total
self-energy is given in Eq. (15), where $\Sigma^{0}$ is the normal self-energy
without any symmetry breaking given by Eq. (67). Here, we calculate
$\Sigma^{0}$ at each $T$ by subtracting its static and Hermitian part,
$\Sigma^{0,{\rm
H}}({{\bm{k}}})\equiv(\Sigma^{0}({{\bm{k}}},+i\delta)+\Sigma^{0}({{\bm{k}}},-i\delta))/2$,
in order to fix the shape of the FS Text-SCVC . Next, we derive the symmetry
breaking part $\delta t$ self-consistently based on the following procedure:
(a) We first calculate $\displaystyle
S_{k}\equiv\frac{T}{N}\sum_{q}G_{k+q}[\Sigma]W_{q}[\Sigma]$, where
$G_{k}[\Sigma]$ and $W_{q}[\Sigma]$ are functions of the total self-energy.
(b) Next, we derive $\delta t$ as

$\displaystyle\delta t_{k}=(1-P_{0})S_{k},$ (56)

where $P_{0}$ is the projection operator for the totally-symmetric ($A_{1g}$)
channel. (c) The total self-energy is given as $\Sigma=\Sigma_{0}+\delta t$.
We repeat (a)-(c) till $\delta t$ converges.

It is easy to show that the $\delta t$-linear term of $S_{k}$ gives the right-
hand side of the linearized DW equation (39) with the kernel function in Eq.
(68). Thus, the full-DW equation is equivalent to the linearized-DW equation
when $\delta t$ is very small.

Figure 5 (a) represents the obtained renormalized order parameter $\delta
t_{m}^{*}({{\bm{k}}})=z_{m}\delta t_{m}({{\bm{k}}})$ for $m=xz$ and $m=xy$ at
Y-point. The obtained $T_{\rm c}=12$meV completely coincides with that given
by the linearized DW equation. The nematic order occurs as the second-order,
and the averaged order parameter $\delta t_{\rm av}^{*}\equiv(|\delta
t_{xz}^{*}|+|\delta t_{xy}^{*}|)/2$ at Y-point is about $2T_{\rm c}$ at
$T=5$meV. Thus, the present theory gives the ratio $\delta t_{\rm
av}^{*}/T_{\rm c}\sim 2$, which is similar to the ratio $\Delta_{0}^{*}/T_{\rm
c}^{\rm SC}\sim 2$ in the BCS theory. Thus, both the development of $\chi_{\rm
nem}$ above $T_{\rm c}$ and the nematic order parameter below $T_{\rm c}$ are
well explained by the present theory.

The relation $2|\delta t_{xz}^{*}|\approx|\delta t_{xy}^{*}|$ in Fig. 5 (a)
indicates that both the $(xz+yz)$-orbitals and $xy$-orbital equally contribute
to the nematic order. In Appendix D, we explain the relative phase between
$xz$-orbital and $xy$-orbital form factors, $\delta t_{xz,xz}(0,\pi)\delta
t_{xy,xy}(0,\pi)<0$, on the basis of the Ginzburg-Landau analysis. This
relation is significant for the Lifshitz transition below $T_{\rm c}$ as we
will explain below.

Figure 5:  (a) Obtained renormalized symmetry breaking self-energy $\delta
t_{m,m}^{*}({{\bm{k}}})$ for $m=xz,xy$ at Y-point, derived from the full DW
equation at $r=0.40$. The second-order transition occurs at $T_{\rm c}=12$meV,
which is consistent with the linearized DW equation analysis in Fig. 4 (d).
$T_{L}\ (=9{\rm meV})$ is the Lifshitz transition temperature. (b) FS in the
nematic state at $T=5{\rm meV}\ (<T_{L})$. The electron-pocket around Y-point
disappears due to the $xy$-orbital form factor. In addition, the hole-pocket
is elongated along the $k_{y}$-axis due to the $xz,yz$-orbital polarization
with the sign-reversal in ${{\bm{k}}}$-space. (c) Renormalized band structure
in the normal state. (d) Renormalized band structure in the nematic state.

Due to the nematic order parameter, the FS with $C_{4}$-symmetry in Fig. 4 (a)
is deformed to the $C_{2}$-symmetry FS depicted in Fig. 5 (b) at $r=0.40$ and
$T=5$meV. The corresponding band-dispersions in the normal state and that in
the nematic state are shown in Figs. 5 (c) and (d), respectively. They are
renormalized by the factor $z\sim 5$ for $xy$ orbital and $z\sim 3$ for
$xz$,$yz$ orbitals. The original band-dispersion is shown in Fig. 10 in
Appendix C. The hole-pocket is elongated along the $k_{y}$-axis due to the
$xz,yz$-orbital polarization Onari-FeSe ; FeSe-ARPES-Suzuki ; FeSe-ARPES-Zhang
. Interestingly, the electron-pocket around Y-point disappears in the nematic
state due to the $d_{x^{2}-y^{2}}$-wave form factor on the $xy$-orbital. This
nematic Lifshitz transition has been confirmed by many angle-resolved
photoemission spectroscopy (ARPES) studies FeSe-Lif1 ; FeSe-Lif2 . The
relative phase between two form factors, $\delta t_{xz,xz}(0,\pi)\delta
t_{xy,xy}(0,\pi)<0$, originates from the kinetic energy gain due to the
pseudogap formation by the Lifshitz transition.

In addition, in the nematic phase, the hole-pocket is elongated along the
$k_{y}$-axis due to the $xz,yz$-orbital polarization with the sign-reversal in
${{\bm{k}}}$-space. This has been also confirmed by ARPES studies FeSe-ARPES-
Suzuki ; FeSe-ARPES-Zhang . Thus, experimental key findings in the nematic
states are satisfactorily reproduced by the present “paramagnon interference
mechanism” Onari-SCVC ; Yamakawa-PRX2016 ; Onari-FeSe ; Tazai-JPSJ-fRG .

###  IV.3 Connection between above and below $T_{\rm c}$

In the previous subsection, we derived the nematic self-energy
$\Sigma\equiv\Sigma^{0}+\delta t$ based on the full DW equation. The derived
nematic state corresponds to the stationary point of the LW grand potential
$\Omega_{\rm LW}[\Sigma]$, as proved in Sect. II. Because the obtained nematic
state is thermodynamically stable, we can calculate the nematic susceptibility
below $T_{\rm c}$ on the bases of the linearized DW equation with the nematic
self-energy.

Figure 6 exhibits the eigenvalue of the linearized DW equation in the nematic
state ($\Sigma\equiv\Sigma^{0}+\delta t$), $\lambda_{\rm full-DW}$, in the
case of $r=0.40$. (We also show the DW equation eigenvalue with the normal
self-energy $\Sigma^{0}$, $\lambda^{\prime}(T)$, for reference.) We see that
$\lambda(T)$ reaches unity at $T=T_{\rm c}$, while it monotonically decreases
for $T<T_{\rm c}$. We find that $1-\lambda(T)\approx|1-\lambda^{\prime}(T)|$
for $T\leq T_{\rm c}$, as naturally expected in the GL theory. The beautiful
numerical result in Fig. 6 means that the nematic-state derived from the full-
DW method corresponds to the stationary point of $\Omega$ very accurately.
Thus, electronic states of FeSe both above and below $T_{\rm c}$ are
understood in a unified way based on the paramagnon interference mechanism.

The nematic susceptibility in the nematic state $T<T_{\rm c}$ is $\chi_{\rm
nem}(T=0)=\chi^{0f}(0)/(1-\lambda(T=0))$, where $1-\lambda(T=0)\approx R$.
Since $R\ (\propto T_{\rm c}^{3})$ is much smaller than unity as shown in Fig.
4 (e), the present study clarified that sizable nematic fluctuations remain in
the nematic phase in Fe-based superconductors. This is an important
information to understand the pairing mechanism in FeSe.

Figure 6:  Eigenvalue of the DW equation in the nematic state $\lambda(T)$
with the “ nematic self-energy $\Sigma^{0}+\delta t$”, where $\delta t$ is
derived from the full DW equation. Because $\Sigma^{0}+\delta t$ represents
the thermal equilibrium state, the nematic susceptibility $\bar{\chi}_{\rm
nem}(T)=(1-\lambda(T))^{-1}$ is positive and diverges at $T=T_{\rm c}$.

##  V Discussions

In this section, we discuss important unsolved properties in Fe-based SCs
based on the present theory. We analyze the specific heat jump at $T=T_{\rm
c}$ in Sect. V.1, and calculate the $T$-dependence of $\chi_{\rm nem}$ near
the nematic QCP in Sect. V.2.

###  V.1 Jump in the specific heat at $T_{\rm c}$ in FeSe

From the stationary point of the free energy
($\left.{\partial}\Omega_{\rm{DW}}(T,\mu,\phi)/{\partial}\phi\right|_{T,\mu}=0$),
we obtain $\phi=0$ and $\phi=\sqrt{\frac{-a(T)}{b}}$ above and below $T_{\rm
c}$, respectively. Here, we assume a simple $T$-dependence of $a$,
$a(T)=\dot{a}(T-T_{\rm c})$, where
$\dot{a}=-\left.\frac{da}{dT}\right|_{T_{\rm c}}\ (>0)$. Then, we obtain the
BCS-like order parameter for $T<T_{\rm c}$ is

$\displaystyle\phi=\sqrt{\frac{\dot{a}}{b}(T_{\rm c}-T)}.$ (57)

Then, the order parameter at $T=0$ is $\phi_{0}\approx\sqrt{(\dot{a}T_{\rm
c})/b}$.

Based on the GL free energy, we discuss the jump of the heat capacity $\Delta
C$ due to the DW transition, which is calculated by

$\displaystyle\frac{\Delta C_{DW}}{T_{\rm
c}}=-\left.\frac{d^{2}\Omega_{\rm{DW}}}{dT^{2}}\right|_{T=T_{\rm
c}}=(\dot{a}T_{\rm c})\left(\frac{\phi_{0}}{T_{\rm c}}\right)^{2}.$ (58)

Now, we calculate $\dot{a}$ based on Eq. (50). In FeSe, $\chi^{0f}(0)$ is
almost independent of $T$, while
$-\dot{\lambda}\equiv-\left.\frac{d\lambda}{dT}\right|_{T_{\rm c}}$ takes a
large positive value as shown in Fig. 4 (d), due to the AL-type VCs in the
kernel function $I$. Then, we obtain

$\displaystyle(\dot{a}T_{\rm c})$ $\displaystyle=$ $\displaystyle
2\chi^{0f}(0)(-\dot{\lambda}T_{\rm c}).$ (59)

Note that $\chi^{0f}(0)$ is equal to the DOS projected by the form factor $f$,
$D^{f}(0)\equiv\frac{1}{N}\sum_{k}\delta({\epsilon}_{{\bm{k}}}-\mu)f_{k}^{2}$,
in the absence of the self-energy.

Hereafter, we explicitly consider the mass-enhancement factor due to the self-
energy, $z^{-1}\equiv 1-{\partial}{\rm
Re}\Sigma/{\partial}{\epsilon}|_{{\epsilon}=\mu}$. The relation $z^{-1}\gg 1$
holds in general strongly correlated metals. In the Fermi liquid theory, the
Green function is given as
$G_{k}=z/(i{\epsilon}_{n}-z({\epsilon}_{{\bm{k}}}-\mu))$. Then, the DOS is
changed to $\chi^{0f}_{z}(0)=z\chi^{0f}_{z=1}(0)=zD^{f}(0)$, and the observed
renormalized order parameter is $\phi_{0}^{*}\equiv z\phi_{0}$. Thus, $\Delta
C_{\rm DW}/T_{\rm c}$ due to the nematic transition is given by

$\displaystyle\frac{\Delta C_{\rm{DW}}}{T_{\rm
c}}=2z^{-1}D^{f}(0)R\left(\frac{\phi_{0}^{*}}{T_{\rm c}}\right)^{2},$ (60)

where $R=(-\dot{\lambda}T_{\rm c})$. As we summarized in Fig. 4 (e), $R\sim
0.3$ for $r=0.40$, and $R\sim 0.1$ for $r=0.36$. In contrast, $\Delta C_{\rm
SC}/T_{\rm c}^{\rm SC}$ due to the BCS superconductivity is ${\Delta
C_{\rm{SC}}}/{T_{\rm c}^{\rm SC}}=2z^{-1}D^{f}(0)\left({\psi_{0}^{*}}/{T_{\rm
c}^{\rm SC}}\right)^{2}$ with $\Delta=\psi\cdot f$, which corresponds to $R=1$
in Eq. (60). ($\psi_{0}^{*}\equiv z\psi_{0}$ is the observed gap function.)
Because ${\psi_{0}^{*}}/{T_{\rm c}^{\rm SC}}\approx 2$, we obtain ${\Delta
C_{\rm{SC}}}/{T_{\rm c}^{\rm SC}}=8z^{-1}D^{f}(0)$, which is close to
$9.4z^{-1}D^{f}(0)$ in the BCS theory. Because ${\phi_{0}^{*}}/{T_{\rm c}^{\rm
SC}}\approx 2$ in the present numerical study, we obtain the relation

$\displaystyle\frac{\Delta C_{\rm{DW}}}{T_{\rm c}}\sim R\frac{\Delta C_{\rm
SC}}{T_{\rm c}^{\rm SC}}.$ (61)

In the present theory, $R\propto T_{\rm c}^{b}$ with $b\sim 3$ for $r=0.40\sim
0.34$ ($T_{\rm c}=12\sim 6$meV) as shown in Fig. 4 (e). Thus, the relation
$\displaystyle\frac{\Delta C_{\rm{DW}}}{T_{\rm c}}\propto T_{\rm c}^{b}$ is
predicted by the present theory.

Figure 7:  (a) Obtained eigenvalue $\lambda_{{\bm{q}}}$ for RbFe2As2 model at
$r=0.30$ Onari-B2g . The obtained form factor at ${{\bm{q}}}={\bm{0}}$ is the
$d_{xy}$-symmetry bond order on $xy$-orbital. Here, we take the
renormalization factor $z=1/2$ into account by following Ref. Yamakawa-PRX2016
. (b) $T$-dependence of $\lambda_{{{\bm{q}}}={\bm{0}}}$ in RbFe2As2 model.

Next, we discuss the nematic state in RbFe2As2, which is a heavily hole-doped
Fe-based superconductor. This system exhibits the uniform
(${{\bm{q}}}={\bm{0}}$) nematic order. Interestingly, the observed nematicity
possesses the $d_{xy}$-wave symmetry, whose director is $45$-degrees rotated
from the $d_{x^{2}-y^{2}}$-wave nematicity in FeSe. Figure 7 (a) shows the DW
equation eigenvalue in the RbFe2As2 Hubbard model, which was introduced in
Ref. Onari-B2g . The obtained form factor at ${{\bm{q}}}={\bm{0}}$ is the
$d_{xy}$-wave bond order on $xy$-orbital, as we revealed in Ref. Onari-B2g .
Here, we set $\Sigma^{0}=(1-z^{-1})(i{\epsilon}_{n}-\mu)$ by following Ref.
Yamakawa-PRX2016 , instead of calculating the FLEX self-energy. We set the
renormalization factor $z=1/2\ (=m/m^{*})$. Then, $T_{\rm c}$ is renormalized
to be $T_{\rm c}^{*}=zT_{\rm c}$, and $\lambda_{{{\bm{q}}}}^{z=1}(T)$ is equal
to $\lambda_{{{\bm{q}}}}^{z}(zT)$ Yamakawa-PRX2016 .

Figure 7 (b) exhibits the $T$-dependence of $\lambda_{{{\bm{q}}}={\bm{0}}}$.
Here, the relations $R\sim 0.1$ when $T_{\rm c}\sim 10$meV and $R\propto
T_{\rm c}^{b}$ ($b\sim 2.5$) are obtained in the RbFe2As2 model. This result
indicates that the relation $R\sim 0.01$ is satisfied at $T_{\rm c}\sim 40$K.
In fact, in FeSe model, we obtained the relation $R\propto T_{\rm c}^{b}$
($b\sim 3$) at low $T_{\rm c}$ ($=6\sim 12{\rm meV}$) by using fine
${{\bm{k}}}$-meshes ($64^{2}$) and many Matsubara frequencies ($8192$) to
obtain reliable results at low $T$; see Fig. 4 (e). Similar relation is
expected to be realized in other Fe-based superconductor models within the
same paramagnon interference mechanism. Therefore, the present theory gives a
natural explanation why $\displaystyle\frac{\Delta C_{\rm{DW}}}{T_{\rm c}}$ in
RbFe2As2 ($T_{\rm c}\approx 40$K) reported in Ref. Mizukami-B2g  is much
smaller than that in FeSe ($T_{\rm c}\approx 90$K) in Ref. FeSe-specific-heat
.

###  V.2 CW/non-CW behavior in Nematic susceptibility

Here, we discuss the nematic susceptibility $\chi_{\rm nem}$ due to the
electron correlation. According to Ref. Onari-B2g ; Onari-smectic , the
nematic susceptibility is given as

$\displaystyle\chi_{\rm nem}=zD^{f}(0)\frac{1}{1-\lambda(T)},$ (62)

where $\lambda(T)$ is the eigenvalue of the DW equation with the “normal self-
energy” $\Sigma^{0}$. Figure 8 (a) shows the normalized susceptibility
$\bar{\chi}_{\rm nem}\equiv\chi_{\rm nem}/\chi_{\rm nem}^{0}=1/(1-\lambda(T))$
for $r=0.40$ and $0.34$, which corresponds to $T_{\rm c}=12$meV and $6.2$meV
respectively. In both cases, $\bar{\chi}_{\rm nem}$ follows the CW behavior at
higher temperatures ($T>T^{*}\sim 10$meV). In contrast, at lower temperatures
($T<T^{*}$) for $r=0.34$, $\bar{\chi}_{\rm nem}$ exhibits a clear deviation
from the CW behavior. The Weiss temperature $T_{0}$ is derived from Fig. 8
(b): We see that $\lambda(T)$ changes from $T$-linear to $T^{2}$-like at
$T\sim T^{*}\sim 8$meV. Similar Fermi-liquid behavior in $\lambda(T)$ is also
recognized in Fig. 4 (d). This result is natural because the system is far
from the magnetic QCP. Thus, the present theory provides a natural explanation
for the non-CW $\chi_{\rm nem}$ near the nematic QCP with $T_{\rm c}\approx 0$
reported in Refs. Hosoi .

Figure 8:  (a) Nematic susceptibility $\bar{\chi}_{\rm nem}=1/(1-\lambda(T))$
for $r=0.34$ and $0.40$ as a function of $T$. $\bar{\chi}_{\rm nem}$ follows
the CW behavior ($C/(T-T_{0})$) at higher temperatures. However, it deviates
from the CW behavior at low temperatures. (b) Derivation of the Weiss
temperature $T_{0}$ from $\lambda(T)$. (c) The “paramagnon interference AL
term” $X_{q=0}$ that is proportional to $\lambda_{q=0}$.
$\Lambda^{f}_{{{\bm{q}}}={\bm{0}}}({\bm{Q}})$ is the three-point vertex. (d)
Obtained $\bar{X}_{q=0}$ as a function of $T$. $\bar{X}_{q=0}$ starts to
saturate at low temperatures, consistently with the deviation from the CW law
in $\bar{\chi}_{\rm nem}$.

Here, we discuss the reason why $\chi_{\rm nem}$ exhibits the CW/non-CW
behavior depending on $r$. According to Eq. (40), the eigenvalue $\lambda_{q}$
at $q=0$ is proportional to $X$ in Eq. (41) because $\chi^{0f}(q=0)$ is almost
constant. In FeSe, the nematic state is mainly caused by the “paramagnon
interference AL term” in $X_{q}$. It is approximately given as

$\displaystyle X_{q=0}$ $\displaystyle=$ $\displaystyle
T\sum_{Q,m}3|\Lambda_{Q,{q=0}}^{f,m}|^{2}(U^{2}\chi^{s,m}_{Q})^{2},$ (63)

where $\chi^{s,m}_{Q}$ is the spin susceptibility for $d$-orbital $m$, and
$\Lambda_{Q,q}^{f,m}$ is the three-point vertex composed of three Green
functions:

$\displaystyle\Lambda_{Q,{q=0}}^{f,m}$ $\displaystyle=$
$\displaystyle-T\sum_{p}({\hat{G}}^{0}_{p}{\hat{f}}_{p}^{q=0}{\hat{G}}^{0}_{p})_{m,m}({\hat{G}}^{0}_{p+Q})_{m,m}\\},$
(64)

where ${\hat{G}}^{0}_{p}$ is the matrix representation of the multiorbital
Green function with FLEX self-energy. Their expressions are shown in Fig. 8
(c). In a simple single-orbital model, the analytic expression of Eq. (64) is
given as
$\displaystyle\Lambda_{Q,q=0}^{f}=\frac{1}{N}\sum_{{\bm{k}}}\left(-\frac{{\partial}n({\epsilon}_{{\bm{k}}})}{{\partial}{\epsilon}_{{\bm{k}}}}\right)\frac{1}{{\epsilon}_{{{\bm{k}}}-{{\bm{Q}}}}-{\epsilon}_{{\bm{k}}}}f_{{\bm{k}}}^{q=0}$,
where $n({\epsilon})$ is the Fermi distribution function Kontani-Raman . When
${{\bm{Q}}}\approx{{\bm{Q}}}_{\rm nesting}$, $\Lambda_{Q,q=0}^{f}$ exhibits
strong enhancement at low temperatures due to
$(-{\partial}n({\epsilon}_{{\bm{k}}})/{\partial}{\epsilon}_{{\bm{k}}})$
Yamakawa-PRX2016 ; Kontani-Raman .

In FeSe, the spin Stoner factor ${\alpha}_{S}$ is nearly constant for
$T>T_{\rm c}$. Then, the strong $T$-dependence of $\lambda(T)$ originates from
$\Lambda_{Q,{q=0}}^{f}$, not from $V^{s}_{Q}$, as we explained in Ref.
Yamakawa-PRX2016 . The paramagnon interference magnifies the nematic
susceptibility, but its magnification is nearly constant. To confirm this
fact, we introduce a simplification of $X_{q=0}$ in Eq. (63) as follows:

$\displaystyle\bar{X}_{q=0}$ $\displaystyle=$ $\displaystyle
T\sum_{Q,m}|\Lambda_{Q,{q=0}}^{f,m}|^{2}.$ (65)

Figure 8 (d) shows the numerical result of $\bar{X}_{q=0}$ at $r=0.34$. For
$T\gtrsim 10$meV, $\bar{X}_{q=0}$ exhibit almost perfect $T$-linear behavior.
Its increment at low $T$ originates from the FS nesting. In contrast, at lower
temperatures, $\bar{X}_{q=0}$ starts to saturate when $T$ is smaller than the
nesting energy scale Yamakawa-PRX2016 . This saturation gives rise to the non-
CW behavior of $\chi_{\rm nem}$ near the nematic QCP ($T_{\rm c}\sim 0$) as
shown in Fig. 8 (a).

Thus, the paramagnon interference mechanism satisfactorily explains both the
CW behavior of $\chi_{\rm nem}$ above $T^{*}(\sim 10{\rm meV})$ and its non-CW
behavior below $T^{*}$. These behaviors are actually observed in various Fe-
based superconductors near the nematic QCP: Ba(Fe,T)2As2 with T=Co,Ni,
(Ba,A)Fe2As2 with A=K,Rb, and Fe(Se,Pn) with Pn=Te,S Fisher-Science2016 ;
Shibauchi-nemQCP ; Terashima . (In this mechanism, $\Lambda_{Q,{q=0}}^{f}$ is
the coupling constant between the nematicity and the paramagnons, and its
increment leads to large $\chi_{\rm nem}$ at low temperatures.) Once the
nematic order is established below $T_{c}$, the spin Stoner factor
${\alpha}_{S}$ increases as we discussed in Onari-FeSe ; saigo .

In the present mechanism, the deviation from the CW behavior of $\chi_{\rm
nem}$ below $T^{*}$ is equal to the Fermi-liquid behavior
$\lambda(0)-\lambda(T)\propto T^{2}$. This deviation is naturally expected
when the nematic QCP is well separated from the magnetic QCP, even in the
absence of impurities. We will discuss this point in the Summary section.

##  VI Summary

In this paper, we derived a formally exact “density-wave (DW) equation”, by
introducing the form factor of the DW state $\delta
t^{{{\bm{q}}}{\sigma}}_{{{\bm{k}}}}$ into the LW theory. Its solution
automatically satisfies the extremum condition of the grand potential. By
solving the DW equation, the optimized form factor and its wavevector are
uniquely obtained for both above and below $T_{\rm c}$. This formalism enables
us to perform the Baym-Kadanoff conserving approximation that is essential to
obtain thermodynamic stable states. In addition, we derive an exact expression
of the Ginzburg-Landau (GL) free energy, $F\propto a_{{\bm{q}}}\phi^{2}$,
where $\phi$ is the amplitude of the DW order at wavevector ${{\bm{q}}}$. The
coefficient $a_{{\bm{q}}}\ [\approx
a_{{{\bm{q}}}_{0}}+\frac{1}{2}\sum_{\mu,\nu}c_{\mu,\nu}(q^{\mu}-q_{0}^{\mu})(q^{\nu}-q_{0}^{\nu})]$
is uniquely related to the eigenvalue of the DW equation $\lambda_{{\bm{q}}}$.
This formalism enables us to calculate various thermodynamic properties of the
DW state.

In the second part, we analyzed the nematic state in FeSe based on the derived
DW equation based on a realistic multiorbital Hubbard model with one single
parameter $r$. We explained the following key experiments in Fe-based SCs: (i)
Lifshitz transition due to bond+orbital order FeSe-Lif1 ; FeSe-Lif2 . (ii-1)
The CW behavior of $\chi_{\rm nem}\propto 1/|1-\lambda(T)|$ at higher-
temperatures; $1-\lambda(T)\propto T_{0}-T$. (ii-2) Deviation from the CW
behavior of $\chi_{\rm nem}$ at low temperatures near the nematic QCP without
magnetic criticality; $\lambda(0)-\lambda(T)\propto T^{2}$ Hosoi ; B2g-Ishida
; Terashima . (iii) A scaling relation $\Delta C/T_{\rm c}\propto T_{\rm
c}^{b}$ ($b\sim 3$) that naturally explains the smallness of $\Delta C/T_{\rm
c}$ reported in several nematic systems Sato-CDW ; Murayama-CDW ; Mizukami-B2g
. This is because the gain of the free energy in the nematic transition is
much smaller than that in the SC state. In addition, we explain the (iv)
Nematic QCP away from the magnetic QCP observed in Fe(Se,S), Fe(Se,Te),
Shibauchi-nemQCP  and Na(Fe,Co)As Zheng-NaFeAs .

The present theory naturally explains the essential points (i)-(iv). Thus, it
is concluded that the nematicity in FeSe is the bond+orbital order due to the
“paramagnon interference mechanism” depicted in Fig. 8 (c) Onari-SCVC ;
Yamakawa-PRX2016 ; Onari-FeSe ; Tazai-JPSJ-fRG .

The behavior (ii-2) has been observed in various Fe-based superconductors near
the nematic QCP: Ba(Fe,T)2As2 with T=Co,Ni, (Ba,A)Fe2As2 with A=K,Rb, and
Fe(Se,Pn) with Pn=Te,S Fisher-Science2016 ; Shibauchi-nemQCP . This behavior
is frequently ascribed to the impurity-induced Griffiths phase, while it is
widely observed insensitively to the impurity potential strength. (For
example, the quantum oscillation is observed in Te- and S-doped FeSe.) In the
present theory, the behavior (ii-2) is naturally explained when the nematic
QCP is well separated from the magnetic QCP, even in the absence of
impurities. It is useful to verify the relation $\lambda(0)-\lambda(T)\propto
T^{2}$ experimentally. In the present mechanism, the increment of $\chi_{\rm
nem}$ at low $T$ originates from the $T$-dependence of
$\Lambda_{Q,{{{\bm{q}}}=\bm{0}}}$ in Eq. (64) Yamakawa-PRX2016 , and the self-
energy due to thermal spin fluctuations Kontani-ROP ; Onari-transport-FeSe  is
also important to derive a perfect CW behavior of $\chi_{\rm nem}$. The self-
energy due to nematic fluctuations will be also important ($|\Sigma_{\rm
nem}|\gtrsim|\Sigma_{\rm FLEX}|$) adjacent to the nematic QCP, unless the
dynamical nematic fluctuations are suppressed by the acoustic phonons. This is
an important future issue. It is considered that a perfect CW behavior for
$30$-$250$K observed in BaFe2(As0.3P0.7)2 Fisher-Science2016  is ascribed to
the magnetic criticality due to the coincidence of the nematic and magnetic
QCPs.

The present theory paves the way for understanding various unconventional
phase transition systems for both above and below $T_{\rm c}$. For example,
the analysis of the odd-parity DW order accompanying spontaneous current,
which has been reported in cuprates and kagome metals recently, is an
important future problem Tazai-cLC ; Kontani-sLC . The (local and/or non-
local) multipole order physics in $5d$\- and $f$-electron systems with strong
spin-orbit interaction is another important future issue Tazai-CeB6 . In
addition, it is important to develop the numerical method beyond the one-loop
approximation. the functional-renormalization-group method Metzner-RMP ;
Honerkamp ; Tsuchiizu-PRL2013 ; Tazai-cLC ; Tazai-JPSJ-fRG , which is
equivalent to the parquet equation, would be useful to obtain a reliable
kernel function of the DW equation. The exotic superconductivity mediated by
the DW fluctuations Kontani-PRL2010 ; Tazai-kagome ; Yamakawa-FeSe-underP ;
Yamakawa-FeSe-SC  would be a very interesting future problem.

###### Acknowledgements.

We are grateful to E.-G. Moon for very enlightening discussions. We also thank
Y. Matsuda, T, Shibauchi, T. Hashimoto, and Y. Mizukami for fruitful
discussions on experiments. This work is supported by Grants-in-Aid for
Scientific Research (KAKENI) Research (No. JP20K22328, No. JP20K03858, No.
JP19H05825, No. JP18H01175) from MEXT of Japan.

##  Appendix A Justification of Onari-Kontani approximation in the DW equation

Here, we derive the DW equation introduced by Onari and Kontani et al. in
Refs. Onari-FeSe ; Kawaguchi-CDW ; Tsuchiizu-CDW ; Tazai-JPSJ-fRG , which has
been applied to iron-based and cuprate superconductors, from the exact DW
equation given in Eq. (36). In the present calculation, we apply the
fluctuation-exchange (FLEX) approximation for $\Phi_{\rm{FLEX}}[G]$. It is
given as Yanase

$\displaystyle\Phi_{\rm FLEX}$ $\displaystyle=$ $\displaystyle T\sum_{q}{\rm
Tr}\left\\{\frac{3}{2}\ln(1-U^{s}\chi^{0}_{q})+\frac{1}{2}\ln(1-U^{c}\chi^{0}_{q})\right\\}$
(66) $\displaystyle+\frac{T}{4}\sum_{q}{\rm
Tr}\left\\{(U^{s}\chi^{0}_{q})^{2}+(U^{c}\chi^{0}_{q})^{2}\right\\}$
$\displaystyle+T\sum_{q}{\rm
Tr}\left\\{\frac{3}{2}U^{s}\chi^{0}_{q}+\frac{1}{2}U^{c}\chi^{0}_{q}\right\\},$

which is expressed in Fig.9 (a). Here, $U^{s(c)}$ is the spin-channel (charge-
channel) Coulomb interaction, and $U^{s}=-U^{c}=U$ in the single-orbital
Hubbard model. Their matrix expressions in multiorbital systems are introduced
in the next section.

In the framework of the conserving approximation, the first-order derivative
of $\Phi_{\rm{FLEX}}[G]$ gives the self-energy $\Sigma$. It is expressed as
Bickers ; Onari-SCVCS

$\displaystyle\Sigma^{0}_{k}=\frac{T}{N}\sum_{q}G^{0}_{k+q}W^{0}_{q},$ (67)

which is expressed in Fig.9 (b). Here, $\displaystyle
W^{0}_{q}=\left(\frac{3}{2}V^{s}_{q}+\frac{1}{2}V^{c}_{q}\right)$,
$V^{x}_{q}=U^{x}+U^{x}\chi^{x}_{q}U^{x}$ ($x=s,c$), and
$\chi^{x}_{q}=\chi^{0}(q)/(1-U^{x}\chi^{0}(q))$.

Finally, we derive the irreducible four-point vertex $I$ from the second
derive of $\Phi_{\rm FLEX}[G]$. The derived charge-channel kernel function in
the DW equation (39) for $x=c$ is given by Kawaguchi-CDW ; Kontani-sLC

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!I^{c}_{kk^{\prime}q}=-\frac{3}{2}V^{s}_{k-k^{\prime}}-\frac{1}{2}V^{c}_{k-k^{\prime}}$
$\displaystyle+\frac{T}{N}\sum_{p}[\frac{3}{2}V^{s}_{p+q}V^{s}_{p}+\frac{1}{2}V^{c}_{p+q}V^{c}_{p}]G^{0}_{k-p}G^{0}_{k^{\prime}-p}$
$\displaystyle+\frac{T}{N}\sum_{p}[\frac{3}{2}V^{s}_{p+q}V^{s}_{p}+\frac{1}{2}V^{c}_{p+q}V^{c}_{p}]G^{0}_{k-p}G^{0}_{k^{\prime}+p},$
(68)

which is depicted in Fig.9 (c). Note that the double-counting $U^{2}$-terms in
Eqs. (68) and (67) should be subtracted properly.

Figure 9:  (a) Diagrammatic expressions of $\Phi_{\rm FLEX}[G]$ in the one-
loop approximation. For simplicity, diagrams in the single-orbital Hubbard
models are shown, while we study multi-orbital Hubbard models in this paper.
(b) Self-energy $\Sigma_{k}=\delta\Phi_{\rm FLEX}/\delta G_{k}$. (c)
Linearized DW equation with the kernel function $I$ derived from the second
derivative of $\Phi_{\rm FLEX}[G]$ with respect to $G$. The Maki-Thompson (MT)
term and the Aslamazov-Larkin (AL) terms give the “fluctuation-induced
interaction for the DW”.

Although the DW equation with the kernel in Eq. (68) is an approximation, it
satisfies the Baym-Kadanoff’s conserving laws by introducing $\Sigma_{\rm
FLEX}$. That is, the solution of the DW equation is the thermal equilibrium
state derived from the stationary condition of $\Omega_{\rm{FLEX}}$. Thus, the
Onari-Kontani type DW equation Onari-FeSe ; Kawaguchi-CDW ; Tsuchiizu-CDW ;
Tazai-JPSJ-fRG  is given by dropping $\Sigma_{\rm FLEX}$ from Eq. (68).
Furthermore, the present exact DW equation is useful to go beyond the Onari-
Kontani’s approximation.

##  Appendix B GL free energy from $\Omega[G]$

In the main text, we derive the GL free energy based on the grand potential
$\Omega[\Sigma]$. Here, we show a different way to obtain GL free energy by
using $\Omega[G]$ starting from Eq. (14) in the main text. The coefficient $a$
is derived from the second functional derivation of $\Omega[G]$ as

$\displaystyle\Omega[\bar{G}]-\Omega[G^{0}]$ $\displaystyle=$
$\displaystyle\sum_{k\sigma
k^{\prime}\sigma^{\prime}}\left.\frac{{\delta}^{2}\Omega[G]}{{\delta}G_{k^{\prime}\sigma^{\prime}}{\delta}G_{k\sigma}}\right|_{G=G^{0}}$
(69)
$\displaystyle\times{\delta}G_{k^{\prime}\sigma^{\prime}}{\delta}G_{k\sigma}$
$\displaystyle\frac{{\delta}^{2}\Omega[G]}{{\delta}G_{k^{\prime}\sigma^{\prime}}{\delta}G_{k\sigma}}$
$\displaystyle=$
$\displaystyle\frac{{\delta}}{{\delta}G_{k^{\prime}\sigma^{\prime}}}\left\\{G_{k\sigma}^{-1}-(G^{\rm{\rm{free}}}_{k\sigma})^{-1}+\Sigma_{k\sigma}[G]\right\\}$
$\displaystyle=$ $\displaystyle-
G_{k\sigma}^{-2}\delta_{kk^{\prime}}-\frac{{\delta}\Sigma_{k\sigma}[G]}{{\delta}G_{k^{\prime}\sigma^{\prime}}}$
$\displaystyle\left.\frac{{\delta}^{2}\Omega[G]}{{\delta}G_{k\sigma}{\delta}G_{k^{\prime}\sigma^{\prime}}}\right|_{G=G^{0}}$
$\displaystyle=$
$\displaystyle-(G^{0}_{k\sigma})^{-2}\delta_{kk^{\prime}}-I_{kk^{\prime}}^{\sigma\sigma^{\prime}}.$
(70)

Thus, we obtain the following results

$\displaystyle\Omega[G]-\Omega[G^{0}]$ $\displaystyle=$
$\displaystyle\frac{T}{N}\sum_{k\sigma}(G^{0}_{k\sigma})^{-2}({\delta}G_{k\sigma})^{2}$
(71) $\displaystyle-\frac{T}{N}\sum_{k\sigma
k^{\prime}\sigma^{\prime}}I_{kk^{\prime}}^{\sigma\sigma^{\prime}}{\delta}G_{k\sigma}{\delta}G_{k^{\prime}\sigma^{\prime}}$
$\displaystyle=$ $\displaystyle\chi^{f}(0)\left\\{1-\lambda\right\\}\phi^{2},$

where we use the relation ${\delta}G_{k\sigma}=(G^{0}_{k\sigma})^{2}\delta
t_{k\sigma}$ and set $t_{k\sigma}=f_{k\sigma}\cdot\phi$. Therefore, the
coefficient $a$ of GL free energy from $\Omega[G]$ is given by

$\displaystyle a$ $\displaystyle=$
$\displaystyle\chi^{f}(0)\left\\{1-\lambda\right\\}\phi^{2}.$ (72)

At $T=T_{\rm c}$, the obtained $a$ by $\Omega[G]$ is the same as the one by
$\Omega[\Sigma]$ in the main text, while it becomes different at $T\neq T_{\rm
c}$. Moreover, the expression $a$ by $\Omega[G]$ does not reproduce the
results by mean-field approximation. Thus, the expression of Eq. (72) is
correct only at $T=T_{\rm c}$

##  Appendix C Eight-orbital models for FeSe

Here, we introduce the eight-orbital $d$-$p$ models for FeSe. We first derived
the first-principles tight-binding models using the WIEN2k and WANNIER90
codes. Crystal structure parameters of FeSe is given in Refs. FeSe-t . By
following Ref. Yamakawa-PRX2016 , we introduce the $k$-dependent shifts for
orbital $l$, $\delta E_{l}$, in order to obtain the experimentally observed
FSs FeSe-Lif1 ; FeSe-Lif2 . In this paper, we introduce the intra-orbital
hopping parameters into the first-principles FeSe model in order to shift the
$d_{xy}$ orbital band ($d_{xz/yz}$-orbital band) at ($\Gamma$, M, X) points by
(0,-0.35,+0.40) [(-0.22,0,+0.16)] in units of eV. Such level shifts are
introduced by the additional intra-orbital Yamakawa-PRX2016 . The band-
dispersion of the present FeSe model without self-energy is shown in Fig. 10.
Its Fermi surface is given in Fig. 4 (a) in the main text.

In this multiorbital model, the matrix expression of the non-interacting Green
function is given as

$\displaystyle{\hat{G}}_{k}^{\rm
free}=(({\epsilon}_{n}-\mu){\hat{1}}+{\hat{h}}_{{\bm{k}}}^{0})^{-1},$ (73)

where $k\equiv({{\bm{k}}},{\epsilon}_{n})$, and ${\hat{h}}_{{\bm{k}}}^{0}$ is
the matrix expression of the kinetic term, which is given by the Fourier
transformation of the tight-binding model.

Figure 10:  Band-dispersion of the present FeSe model without self-energy. The
renormalized dispersion due to the self-energy is shown in Fig. 5.

We also explain the multiorbital Coulomb interaction term $H_{U}$. The
multiorbital Coulomb interaction term is expressed as
$H_{U}=\frac{1}{4}\sum_{i,{\sigma}{\sigma}^{\prime}}\sum_{ll^{\prime}mm^{\prime}}U^{{\sigma}{\sigma}^{\prime}}_{ll^{\prime}mm^{\prime}}d_{i,l,{\sigma}}^{\dagger}d_{i,m^{\prime},{\sigma}^{\prime}}^{\dagger}d_{i,m,{\sigma}^{\prime}}d_{i,l^{\prime},{\sigma}}$,
where $l,m$ represent the orbital indices, in ${\sigma}=+1\ (-1)$ represents
the $\uparrow\ (\downarrow)$ spin, $i$ is the site index, and
${\hat{U}}^{{\sigma}{\sigma}^{\prime}}=-{\hat{U}}^{c}-{\sigma}{\sigma}^{\prime}{\hat{U}}^{s}$.
The matrix elements of ${\hat{U}}^{s}$ is given by Yamakawa-PRX2016

$({\hat{U}}^{\mathrm{s}})_{l_{1}l_{2},l_{3}l_{4}}=\begin{cases}U_{l_{1},l_{1}},&l_{1}=l_{2}=l_{3}=l_{4}\\\
U_{l_{1},l_{2}}^{\prime},&l_{1}=l_{3}\neq l_{2}=l_{4}\\\
J_{l_{1},l_{3}},&l_{1}=l_{2}\neq l_{3}=l_{4}\\\
J_{l_{1},l_{2}},&l_{1}=l_{4}\neq l_{2}=l_{3}\\\
0,&\mathrm{otherwise}.\end{cases}$ (74)

Also, the bare Coulomb interaction for the charge channel is

$({\hat{U}}^{\mathrm{c}})_{l_{1}l_{2},l_{3}l_{4}}=\begin{cases}-U_{l_{1},l_{1}},&l_{1}=l_{2}=l_{3}=l_{4}\\\
U_{l_{1},l_{2}}^{\prime}-2J_{l_{1},l_{2}},&l_{1}=l_{3}\neq l_{2}=l_{4}\\\
-2U_{l_{1},l_{3}}^{\prime}+J_{l_{1},l_{3}},&l_{1}=l_{2}\neq l_{3}=l_{4}\\\
-J_{l_{1},l_{2}},&l_{1}=l_{4}\neq l_{2}=l_{3}\\\
0.&\mathrm{otherwise}.\end{cases}$ (75)

Here, $U_{l,l}$, $U_{l,l^{\prime}}^{\prime}$ and $J_{l,l^{\prime}}$ are the
first-principles Coulomb interaction terms given in Ref. Miyake .

In the main text, we omit the orbital indices of the Green functions and the
Coulomb interactions to simplify the expressions. It is straightforward to
write the orbital indices of these expressions in multiorbital models by using
Eqs. (73)-(75).

##  Appendix D Relative phase between $xz$-orbital and $xy$-orbital form
factors

As we discussed in the main text, the relation $2|\delta
t_{xz}^{*}|\approx|\delta t_{xy}^{*}|$ in Fig. 5 (a) indicates that both the
$(xz+yz)$-orbitals and $xy$-orbital equally contribute to the nematic order in
FeSe. The $(xz+yz)$-orbital polarization [$xy$-orbital bond order] originates
from the spin fluctuations on $(xz+yz)$-orbitals [$xy$-orbital]. Here, we
discuss the relative phase between $xz$-orbital and $xy$-orbital form factors,
$\delta t_{xz,xz}(0,\pi)\times\delta t_{xy,xy}(0,\pi)<0$, based on the
Ginzburg-Landau analysis.

According to Eq. (48) or (49), the second-order inter-orbital free-energy is
$F^{(2)}_{2,4}=-2\chi^{0f}_{2,4}\phi_{2}\phi_{4}$, where $\chi_{2,4}^{0\delta
t}=-T\sum_{k}G_{2,4}(k)G_{4,2}(k)f_{2,2}(k)f_{4,4}(k)$ and $\delta
t_{m,m}(k)=f_{m,m}(k)\phi_{m}$. We verified numerically that
$\chi_{2,4}^{0f}>0$, which is consistent with the relations
$\chi^{0}_{22,44}(0)<0$ in FeSe and $\delta t_{2,2}(k)\delta t_{4,4}(k)<0$
around Y-point shown in Fig. 4 (c). Thus, the relation $\phi_{2}\phi_{4}>0$ is
realized, and therefore the Lifshitz transition occurs in FeSe. In other
words, the relation $\phi_{2}\phi_{4}>0$ (i.e., $\delta t_{2,2}(0,\pi)\delta
t_{4,4}(0,\pi)<0$) is a direct consequence of the multiorbital band structure
of Fe-based superconductors.

## References

  * (1) A. V. Chubukov, M. Khodas, and R. M. Fernandes, Magnetism, Superconductivity, and Spontaneous Orbital Order in Iron-Based Superconductors: Which Comes First and Why?, Phys. Rev. X 6, 041045 (2016).
  * (2) Y. Yamakawa, S. Onari, and H. Kontani, Nematicity and Magnetism in FeSe and Other Families of Fe-Based Superconductors, Phys. Rev. X 6, 021032 (2016).
  * (3) S. Onari, Y. Yamakawa, and H. Kontani, Sign-Reversing Orbital Polarization in the Nematic Phase of FeSe due to the ${C}_{2}$ Symmetry Breaking in the Self-Energy, Phys. Rev. Lett. 116, 227001 (2016).
  * (4) R. M. Fernandes, A. V. Chubukov, and J. Schmalian, What drives nematic order in iron-based superconductors?, Nat. Phys. 10, 97–104 (2014).
  * (5) M. Yoshizawa, D. Kimura, T. Chiba, S. Simayi, Y. Nakanishi, K. Kihou, C.-H. Lee, A. Iyo, H. Eisaki, M. Nakajima, and S.-i. Uchida, Structural Quantum Criticality and Superconductivity in Iron-Based Superconductor Ba(Fe1-xCox)2As2, J. Phys. Soc. Jpn. 81, 024604 (2012).
  * (6) R. M. Fernandes, L. H. VanBebber, S. Bhattacharya, P. Chandra, V. Keppens, D. Mandrus, M. A. McGuire, B. C. Sales, A. S. Sefat, and J. Schmalian, Effects of Nematic Fluctuations on the Elastic Properties of Iron Arsenide Superconductors, Phys. Rev. Lett. 105, 157003 (2010).
  * (7) A. E. Böhmer, P. Burger, F. Hardy, T. Wolf, P. Schweiss, R. Fromknecht, M. Reinecker, W. Schranz, and C. Meingast, Nematic Susceptibility of Hole-Doped and Electron-Doped ${\mathrm{Ba}\mathrm{F}\mathrm{e}}_{2}{\mathrm{As}}_{2}$ Iron-Based Superconductors from Shear Modulus Measurements, Phys. Rev. Lett. 112, 047001 (2014).
  * (8) Y. Gallais, R. M. Fernandes, I. Paul, L. Chauvière, Y.-X. Yang, M.-A. Méasson, M. Cazayous, A. Sacuto, D. Colson, and A. Forget, Observation of Incipient Charge Nematicity in ${\mathrm{Ba}(\mathrm{Fe}{}_{1-X}\mathrm{Co}{}_{X})}_{2}\mathrm{As}{}_{2}$, Phys. Rev. Lett. 111, 267001 (2013).
  * (9) P. Massat, D. Farina, I. Paul, S. Karlsson, P. Strobel, P. Toulemonde, M.-A. Méasson, M. Cazayous, A. Sacuto, S. Kasahara, T. Shibauchi, Y. Matsuda, and Y. Gallais, Charge-induced nematicity in FeSe, Proc. Natl. Acad. Sci. U.S.A. 113, 9177 (2016).
  * (10) T. Kissikov, R. Sarkar, M. Lawson, B. T. Bush, E. I. Timmons, M. A. Tanatar, R. Prozorov, S. L. Bud’ko, P. C. Canfield, R. M. Fernandes, and N. J. Curro, Uniaxial strain control of spin-polarization in multicomponent nematic order of BaFe2As2, Nat. Commun. 9, 1058 (2018).
  * (11) 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, Rep. Prog. Phys. 74, 124506 (2011).
  * (12) S. Hosoi, K. Matsuura, K. Ishida, H. Wang, Y. Mizukami, T. Watashige, S. Kasahara, Y. Matsuda, and T. Shibauchi, Nematic quantum critical point without magnetism in FeSe1−xSx superconductors, Proc. Natl. Acad. Sci. U.S.A. 113, 8139 (2016).
  * (13) K. Ishida, M. Tsujii, S. Hosoi, Y. Mizukami, S. Ishida, A. Iyo, H. Eisaki, T. Wolf, K. Grube, H. v. Löhneysen, R. M. Fernandes, and T. Shibauchi, Novel electronic nematicity in heavily hole-doped iron pnictide superconductors, Proc. Natl. Acad. Sci. U.S.A. 117, 6424 (2020).
  * (14) H.-H. Kuo, J.-H. Chu, J. C. Palmstrom, S. A. Kivelson, and I. R. Fisher, Ubiquitous signatures of nematic quantum criticality in optimally doped Fe-based superconductors, Science 352, 958 (2016).
  * (15) R. M. Fernandes and J. W. F. Venderbos, Nematicity with a twist: Rotational symmetry breaking in a moiré superlattice, Sci. Adv. 6, eaba8834 (2020).
  * (16) S. Onari and H. Kontani, SU(4) $\text{Valley}+\text{Spin}$ Fluctuation Interference Mechanism for Nematic Order in Magic-Angle Twisted Bilayer Graphene: The Impact of Vertex Corrections, Phys. Rev. Lett. 128, 066401 (2022).
  * (17) H. Nakaoka, Y. Yamakawa, and H. Kontani, Theoretical prediction of nematic orbital-ordered state in the Ti oxypnictide superconductor ${\mathrm{BaTi}}_{2}{(\mathrm{As},\mathrm{Sb})}_{2}\mathrm{O}$, Phys. Rev. B 93, 245122 (2016).
  * (18) 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, Nat. Phys. 13, 1074–1078 (2017).
  * (19) H. Murayama, Y. Sato, R. Kurihara, S. Kasahara, Y. Mizukami, Y. Kasahara, H. Uchiyama, A. Yamamoto, E.-G. Moon, J. Cai, J. Freyermuth, M. Greven, T. Shibauchi, and Y. Matsuda, Diagonal nematicity in the pseudogap phase of HgBa2CuO4+${\delta}$, Nat. Commun. 10, 3282 (2019).
  * (20) K. Kawaguchi, Y. Yamakawa, M. Tsuchiizu, and H. Kontani, Competing Unconventional Charge-Density-Wave States in Cuprate Superconductors: Spin-Fluctuation-Driven Mechanism, J. Phys. Soc. Jpn. 86, 063707 (2017).
  * (21) M. Tsuchiizu, K. Kawaguchi, Y. Yamakawa, and H. Kontani, Multistage electronic nematic transitions in cuprate superconductors: A functional-renormalization-group analysis, Phys. Rev. B 97, 165131 (2018).
  * (22) R. Tazai, Y. Yamakawa, S. Onari, and H. Kontani, Mechanism of exotic density-wave and beyond-Migdal unconventional superconductivity in kagome metal AV3Sb5 (A = K, Rb, Cs), Sci. Adv. 8, eabl4108 (2022).
  * (23) M. L. Kiesel, C. Platt, and R. Thomale, Unconventional Fermi Surface Instabilities in the Kagome Hubbard Model, Phys. Rev. Lett. 110, 126405 (2013).
  * (24) W.-S. Wang, Z.-Z. Li, Y.-Y. Xiang, and Q.-H. Wang, Competing electronic orders on kagome lattices at van Hove filling, Phys. Rev. B 87, 115135 (2013).
  * (25) T. Park, M. Ye, and L. Balents, Electronic instabilities of kagome metals: Saddle points and Landau theory, Phys. Rev. B 104, 035142 (2021).
  * (26) Y. Yamakawa and H. Kontani, Spin-Fluctuation-Driven Nematic Charge-Density Wave in Cuprate Superconductors: Impact of Aslamazov-Larkin Vertex Corrections, Phys. Rev. Lett. 114, 257001 (2015).
  * (27) S. Sachdev and R. La Placa, Bond Order in Two-Dimensional Metals with Antiferromagnetic Exchange Interactions, Phys. Rev. Lett. 111, 027202 (2013).
  * (28) C. Husemann and W. Metzner, Incommensurate nematic fluctuations in the two-dimensional Hubbard model, Phys. Rev. B 86, 085113 (2012).
  * (29) C. Fang, H. Yao, W.-F. Tsai, J. Hu, and S. A. Kivelson, Theory of electron nematic order in LaFeAsO, Phys. Rev. B 77, 224509 (2008).
  * (30) F. Wang, S. A. Kivelson, and D.-H. Lee, Nematicity and quantum paramagnetism in FeSe, Nat. Phys. 11, 959–963 (2015).
  * (31) R. Tazai and H. Kontani, Multipole fluctuation theory for heavy fermion systems: Application to multipole orders in ${\mathrm{CeB}}_{6}$, Phys. Rev. B 100, 241103(R) (2019).
  * (32) R. Tazai and H. Kontani, Fully gapped $s$-wave superconductivity enhanced by magnetic criticality in heavy-fermion systems, Phys. Rev. B 98, 205107 (2018).
  * (33) R. Tazai and H. Kontani, Hexadecapole Fluctuation Mechanism for $s$-wave Heavy Fermion Superconductor CeCu2Si2: Interplay between Intra- and Inter-Orbital Cooper Pairs, J. Phys. Soc. Jpn. 88, 063701 (2019).
  * (34) 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, Time-reversal symmetry-breaking charge order in a kagome superconductor, Nature 602, 245–250 (2022).
  * (35) Q. Wu, Z. X. Wang, Q. M. Liu, R. S. Li, S. X. Xu, Q. W. Yin, C. S. Gong, Z. J. Tu, H. C. Lei, T. Dong, and N. L. Wang, Revealing the immediate formation of two-fold rotation symmetry in charge-density-wave state of Kagome superconductor CsV3Sb5 by optical polarization rotation measurement, arXiv:2110.11306.
  * (36) H. Murayama, K. Ishida, R. Kurihara, T. Ono, Y. Sato, Y. Kasahara, H. Watanabe, Y. Yanase, G. Cao, Y. Mizukami, T. Shibauchi, Y. Matsuda, and S. Kasahara, Bond Directional Anapole Order in a Spin-Orbit Coupled Mott Insulator ${\mathrm{Sr}}_{2}({\mathrm{Ir}}_{1-x}{\mathrm{Rh}}_{x}){\mathrm{O}}_{4}$, Phys. Rev. X 11, 011021 (2021).
  * (37) M. E. Simon and C. M. Varma, Detection and Implications of a Time-Reversal Breaking State in Underdoped Cuprates, Phys. Rev. Lett. 89, 247003 (2002).
  * (38) A. A. Nersesyan, G. I. Japaridze, and I. G. Kimeridze, Low-temperature magnetic properties of a two-dimensional spin nematic state, J. Phys.: Condens. Matter 3, 3353 (1991).
  * (39) E. Berg, E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Striped superconductors: how spin, charge and superconducting orders intertwine in the cuprates, New J. Phys. 11, 115004 (2009).
  * (40) C. J. Halboth and W. Metzner, Renormalization-group analysis of the two-dimensional Hubbard model, Phys. Rev. B 61, 7364 (2000).
  * (41) J. C. S. Davis and D.-H. Lee, Concepts relating magnetic interactions, intertwined electronic orders, and strongly correlated superconductivity, Proc. Natl. Acad. Sci. U.S.A. 110, 17623 (2013).
  * (42) R. Tazai, Y. Yamakawa, M. Tsuchiizu, and H. Kontani, $d$\- and $p$-wave Quantum Liquid Crystal Orders in Cuprate Superconductors, ${\kappa}{}$-(BEDT-TTF)2X, and Coupled Chain Hubbard Models: Functional-renormalization-group Analysis, J. Phys. Soc. Jpn. 90, 111012 (2021).
  * (43) H. Kontani, Y. Yamakawa, R. Tazai, and S. Onari, Odd-parity spin-loop-current order mediated by transverse spin fluctuations in cuprates and related electron systems, Phys. Rev. Research 3, 013127 (2021).
  * (44) A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics, (New York: Dover Publications, Inc., 1975.)
  * (45) J. M. Luttinger and J. C. Ward, Ground-State Energy of a Many-Fermion System. II, Phys. Rev. 118, 1417 (1960).
  * (46) G. Baym and L. P. Kadanoff, Conservation Laws and Correlation Functions, Phys. Rev. 124, 287 (1961).
  * (47) G. Baym, Self-Consistent Approximations in Many-Body Systems, Phys. Rev. 127, 1391 (1962).
  * (48) Y. Suzuki, T. Shimojima, T. Sonobe, A. Nakamura, M. Sakano, H. Tsuji, J. Omachi, K. Yoshioka, M. Kuwata-Gonokami, T. Watashige, R. Kobayashi, S. Kasahara, T. Shibauchi, Y. Matsuda, Y. Yamakawa, H. Kontani, and K. Ishizaka, Momentum-dependent sign inversion of orbital order in superconducting FeSe, Phys. Rev. B 92, 205117 (2015).
  * (49) M. Yi, H. Pfau, Y. Zhang, Y. He, H. Wu, T. Chen, Z. R. Ye, M. Hashimoto, R. Yu, Q. Si, D.-H. Lee, P. Dai, Z.-X. Shen, D. H. Lu, and R. J. Birgeneau, Nematic Energy Scale and the Missing Electron Pocket in FeSe, Phys. Rev. X 9, 041049 (2019).
  * (50) S. S. Huh, J. J. Seo, B. S. Kim, S. H. Cho, J. K. Jung, S. Kim, C. I. Kwon, J. S. Kim, Y. Y. Koh, W. S. Kyung, J. D. Denlinger, Y. H. Kim, B. N. Chae, N. D. Kim, Y. K. Kim, and C. Kim, Absence of Y-pocket in 1-Fe Brillouin zone and reversed orbital occupation imbalance in FeSe, Commun. Phys. 3, 52 (2020).
  * (51) L. C. Rhodes, J. Böker, M. A. Müller, M. Eschrig, and I. M. Eremin, Non-local dxy nematicity and the missing electron pocket in FeSe, npj Quantum Mater. 6, 45 (2021).
  * (52) T. Terashima, Y. Matsushita, H. Yamase, N. Kikugawa, H. Abe, M. Imai, S. Uji, S. Ishida, H. Eisaki, A. Iyo, K. Kihou, C.-H. Lee, T. Wang, and G. Mu, Elastoresistance measurements on ${\mathrm{CaKFe}}_{4}{\mathrm{As}}_{4}$ and ${\mathrm{KCa}}_{2}{\mathrm{Fe}}_{4}{\mathrm{As}}_{4}{\mathrm{F}}_{2}$ with the Fe site of ${C}_{2v}$ symmetry, Phys. Rev. B 102, 054511 (2020).
  * (53) Y. Mizukami, O. Tanaka, K. Ishida, M. Tsujii, T. Mitsui, S. Kitao, M. Kurokuzu, M. Seto, S. Ishida, A. Iyo, H. Eisaki, K. Hashimoto, and T. Shibauchi, Thermodynamic Signatures of Diagonal Nematicity in RbFe2As2 Superconductor, arXiv:2108.13081.
  * (54) K. Ishida, Y. Onishi, M. Tsujii, K. Mukasa, M. Qiu, M. Saito, Y. Sugimura, K. Matsuura, Y. Mizukami, K. Hashimoto, and T. Shibauchi, Pure nematic quantum critical point accompanied by a superconducting dome, arXiv:2202.11674.
  * (55) C. G. Wang, Z. Li, J. Yang, L. Y. Xing, G. Y. Dai, X. C. Wang, C. Q. Jin, R. Zhou, and G.-q. Zheng, Electron Mass Enhancement near a Nematic Quantum Critical Point in ${\mathrm{NaFe}}_{1-x}{\mathrm{Co}}_{x}\mathrm{As}$, Phys. Rev. Lett. 121, 167004 (2018).
  * (56) S. Onari and H. Kontani, Self-consistent Vertex Correction Analysis for Iron-based Superconductors: Mechanism of Coulomb Interaction-Driven Orbital Fluctuations, Phys. Rev. Lett. 109, 137001 (2012).
  * (57) H. J. Schulz, Fermi-surface instabilities of a generalized two-dimensional Hubbard model, Phys. Rev. B 39, 2940(R) (1989).
  * (58) I. Affleck and J. B. Marston, Large-n limit of the Heisenberg-Hubbard model: Implications for high-${T}_{c}$ superconductors, Phys. Rev. B 37, 3774(R) (1988).
  * (59) R. Tazai, Y. Yamakawa, and H. Kontani, Emergence of charge loop current in the geometrically frustrated Hubbard model: A functional renormalization group study, Phys. Rev. B 103, L161112 (2021).
  * (60) E. P. Wohlfarth, Thermodynamic aspects of itinerant electron magnetism, Physica B+C 91, 305 (1977).
  * (61) M. Potthoff, Self-energy-functional approach to systems of correlated electrons, Eur. Phys. J. B 32, 429 (2003).
  * (62) Y. Yanase and M. Ogata, Kinetic Energy, Condensation Energy, Optical Sum Rule and Pairing Mechanism in High-$T_{c}$ Cuprates, J. Phys. Soc. Jpn. 74, 1534 (2005).
  * (63) T. Miyake, K. Nakamura, R. Arita, and M. Imada, Comparison of $Abinitio$ Low-Energy Models for LaFePO, LaFeAsO, BaFe2As2, LiFeAs, FeSe, and FeTe: Electron Correlation and Covalency, J. Phys. Soc. Jpn. 79, 044705 (2010).
  * (64) H. Kontani and M. Ohno, Effect of a nonmagnetic impurity in a nearly antiferromagnetic Fermi liquid: Magnetic correlations and transport phenomena, Phys. Rev. B 74, 014406 (2006).
  * (65) R.-Q. Xing, L. Classen, and A. V. Chubukov, Orbital order in FeSe: The case for vertex renormalization, Phys. Rev. B 98, 041108 (2018).
  * (66) H. Kontani, Anomalous transport phenomena in Fermi liquids with strong magnetic fluctuations, Rep. Prog. Phys. 71, 026501 (2008).
  * (67) S. Onari and H. Kontani, Hidden antiferronematic order in Fe-based superconductor $\mathrm{Ba}{\mathrm{Fe}}_{2}{\mathrm{As}}_{2}$ and NaFeAs above ${T}_{S}$, Phys. Rev. Research 2, 042005(R) (2020).
  * (68) Y. Zhang, M. Yi, Z.-K. Liu, W. Li, J. J. Lee, R. G. Moore, M. Hashimoto, M. Nakajima, H. Eisaki, S.-K. Mo, Z. Hussain, T. P. Devereaux, Z.-X. Shen, and D. H. Lu, Distinctive orbital anisotropy observed in the nematic state of a FeSe thin film, Phys. Rev. B 94, 115153 (2016).
  * (69) S. Onari and H. Kontani, Origin of diverse nematic orders in Fe-based superconductors: ${45}^{\circ}$ rotated nematicity in $A{\mathrm{Fe}}_{2}{\mathrm{As}}_{2}\phantom{\rule{4.0pt}{0.0pt}}(A=\text{Cs},\text{Rb})$, Phys. Rev. B 100, 020507(R) (2019).
  * (70) S. Onari and H. Kontani, Iron-Based Superconductivity, (ed. P.D. Johnson, G. Xu, and W.-G. Yin, Springer-Verlag Berlin and Heidelberg GmbH & Co. K (2015)).
  * (71) A. E. Böhmer, T. Arai, F. Hardy, T. Hattori, T. Iye, T. Wolf, H. v. Löhneysen, K. Ishida, and C. Meingast, Origin of the Tetragonal-to-Orthorhombic Phase Transition in FeSe: A Combined Thermodynamic and NMR Study of Nematicity, Phys. Rev. Lett. 114, 027001 (2015).
  * (72) H. Kontani and Y. Yamakawa, Linear Response Theory for Shear Modulus ${C}_{66}$ and Raman Quadrupole Susceptibility: Evidence for Nematic Orbital Fluctuations in Fe-based Superconductors, Phys. Rev. Lett. 113, 047001 (2014).
  * (73) H. Kontani, T. Saito, and S. Onari, Origin of orthorhombic transition, magnetic transition, and shear-modulus softening in iron pnictide superconductors: Analysis based on the orbital fluctuations theory, Phys. Rev. B 84, 024528 (2011).
  * (74) S. Onari and H. Kontani, In-plane anisotropy of transport coefficients in electronic nematic states: Universal origin of nematicity in Fe-based superconductors, Phys. Rev. B 96, 094527 (2017).
  * (75) W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schönhammer, Functional renormalization group approach to correlated fermion systems, Rev. Mod. Phys. 84, 299 (2012).
  * (76) M. Salmhofer and C. Honerkamp, Fermionic Renormalization Group Flows: Technique and Theory, Prog. Theor. Phys. 105, 1 (2001).
  * (77) M. Tsuchiizu, Y. Ohno, S. Onari, and H. Kontani, Orbital Nematic Instability in the Two-Orbital Hubbard Model: Renormalization-Group + Constrained RPA Analysis, Phys. Rev. Lett. 111, 057003 (2013).
  * (78) H. Kontani and S. Onari, Orbital-Fluctuation-Mediated Superconductivity in Iron Pnictides: Analysis of the Five-Orbital Hubbard-Holstein Model, Phys. Rev. Lett. 104, 157001 (2010).
  * (79) Y. Yamakawa and H. Kontani, Nematicity, magnetism, and superconductivity in FeSe under pressure: Unified explanation based on the self-consistent vertex correction theory, Phys. Rev. B 96, 144509 (2017).
  * (80) Y. Yamakawa and H. Kontani, Superconductivity without a hole pocket in electron-doped FeSe: Analysis beyond the Migdal-Eliashberg formalism, Phys. Rev. B 96, 045130 (2017).
  * (81) N. E. Bickers, D. J. Scalapino, and S. R. White, Conserving Approximations for Strongly Correlated Electron Systems: Bethe-Salpeter Equation and Dynamics for the Two-Dimensional Hubbard Model, Phys. Rev. Lett. 62, 961 (1989).
  * (82) S. Onari, Y. Yamakawa, and H. Kontani, High-${T}_{c}$ Superconductivity near the Anion Height Instability in Fe-Based Superconductors: Analysis of ${\mathrm{LaFeAsO}}_{1-x}{\mathrm{H}}_{x}$, Phys. Rev. Lett. 112, 187001 (2014).
  * (83) A. E. Böhmer, F. Hardy, F. Eilers, D. Ernst, P. Adelmann, P. Schweiss, T. Wolf, and C. Meingast, Lack of coupling between superconductivity and orthorhombic distortion in stoichiometric single-crystalline FeSe, Phys. Rev. B 87, 180505(R) (2013).