# Inter-orbital nematicity and the origin of a single electron Fermi pocket in
FeSe

Daniel Steffensen1, Andreas Kreisel2, P. J. Hirschfeld1,3, Brian M. Andersen1
1Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200
Copenhagen, Denmark
2Institut für Theoretische Physik, Universität Leipzig, D-04103 Leipzig,
Germany
3Department of Physics, University of Florida, Gainesville, Florida 32611, USA

###### Abstract

The electronic structure of the enigmatic iron-based superconductor FeSe has
puzzled researchers since spectroscopic probes failed to observe the expected
electron pocket at the $Y$ point in the 1-Fe Brillouin zone. It has been
speculated that this pocket, essential for an understanding of the
superconducting state, is either absent or incoherent. Here, we perform a
theoretical study of the preferred nematic order originating from nearest-
neighbor Coulomb interactions in an electronic model relevant for FeSe. We
find that at low temperatures the dominating nematic components are of inter-
orbital $d_{xz}-d_{xy}$ and $d_{yz}-d_{xy}$ character, with spontaneously
broken amplitudes for these two components. This inter-orbital nematic order
naturally leads to distinct hybridization gaps at the $X$ and $Y$ points of
the 1-Fe Brillouin zone, and may thereby produce highly anisotropic Fermi
surfaces with only a single electron pocket at one of these momentum-space
locations. The associated superconducting gap structure obtained with the
generated low-energy electronic band structure from spin-fluctuation mediated
pairing agrees well with that measured experimentally. Finally, from a
comparison of the computed spin susceptibility to available neutron scattering
data, we discuss the necessity of additional self-energy effects, and explore
the role of orbital-dependent quasiparticle weights as a minimal means to
include them.

##  I Introduction

The research of electronic properties of iron-based superconductors has proven
to be a rich field that challenges our understanding of correlated multi-band
systems. In this respect, particularly the superconducting iron-chalcogenides,
FeSe and doped FeTe, have attracted considerable attention due to their
unusual low-energy electronic states[1, 2, 3]. FeSe enters an orthorhombic
phase near $T_{n}\sim 90$ K without concomitant static magnetic order at lower
temperatures, unlike most of the iron-pnictide superconductors. There is
strong evidence that the rotational symmetry breaking at $T_{n}$ is driven by
the electronic degrees of freedom, and thus constitutes a rare example of
electronic nematicity[4]. In particular, the spectrum of low-energy spin
excitations, thought to drive electron pairing, is extremely anisotropic in
untwinned crystals[5]. The superconducting phase sets in around $T_{c}\sim 9$
K and is therefore generated from an instability of the nematic “normal”
state. Thus, it is not surprising that the superconducting properties also
break rotational symmetry, as is indeed observed experimentally[6]. However,
the origins of the extreme normal state nematicity, as well as the reasons for
the absence of magnetism at ambient pressure, are still being debated[1, 2,
3].

The very high level of rotational anisotropy is evident, e.g., from
theoretical modelling of experimental data measured at low temperatures in
FeSe[3]. Quite generally, the anisotropy inherent to standard electronic bands
in the orthorhombic state, is not enough to explain the much more extreme
anisotropy observed experimentally[7, 8, 6, 9, 5, 10]. Thus, various
interaction-driven feedback effects have been proposed that enhance the
symmetry breaking[6, 11, 12, 13]. In this respect, several works have argued
for orbital-dependent quasiparticle weights $Z_{\alpha}$ and shown how this
effect may reconcile several experimental probes with simple models[14, 6, 9,
11, 15, 16, 17]. The existence of orbital-dependent $Z_{\alpha}$ has been
explored extensively within dynamical mean-field theory (DMFT) and related
methods, and arise naturally in multi-orbital systems with different degrees
of orbital contributions to the Fermi level electronic states[18, 14, 19, 20,
21, 22]. For the specific case of FeSe, it was shown that $Z_{xy}$ must be
strongly reduced from 1, but additionally that a large quasiparticle weight
anisotropy of $Z_{xz}/Z_{yz}\sim 0.2-0.3$ was necessary in order to obtain
agreement with experiments[6, 9, 23] within this approach. The origin of this
rather small ratio remains unsettled, but may be related to strong feedback
effects on the electronic states close to a stripe magnetic quantum critical
point. One of the important implications of these extracted values for the
orbital weights $Z_{\alpha}$ was that the electron pocket at $Y$ should be
nearly completely incoherent, and thus difficult to observe with spectroscopic
probes; indeed, no conclusive observation of this pocket has been reported by
either angular-resolved photoemission spectroscopy (ARPES)[24, 25, 26] or
scanning tunneling microscopy (STM)[9].

Subsequently, however, laser ARPES studies reported significant $d_{yz}$
content on the $\Gamma$-centered hole pocket, as well as a strong $d_{xz}$
component apparently inconsistent with the extreme “orbitally selective”
scenario[27]. Furthermore, very recent synchrotron ARPES measurements on
detwinned FeSe have proposed that the $Y$ pocket may be lifted entirely from
the Fermi surface rather than unobservable due to incoherence[28, 29]. The
explanation of the unusual spectroscopic observations of this pocket is
therefore important not only to obtain the correct low-energy electronic band
structure relevant for (detwinned) nematic FeSe, but also to ultimately
understand the origin of the strong electronic nematicity in FeSe.

Of course, any theoretical approach for explaining the large electronic
anisotropy of FeSe relies on having a reasonable starting point for inclusion
of interactions in the band structure. It was initially proposed, based on
comparison to ARPES data available at the time, that FeSe features nematicity
in the $d_{xz}$ and $d_{yz}$ intra-orbital channels with form factors
transforming as the irreducible representations (IRs) A1g and B1g of the given
point group[30, 31, 32, 33, 34, 35, 36, 37]. More recently, however, several
works have advocated for the relevance of (intra-orbital) nematicity involving
also the $d_{xy}$ orbitals[38, 39, 40, 41, 42, 43]. Until recently, first
principles approaches were unable to stabilize a nematic ground state without
concomitant stripe magnetism. However a recent density functional theory DFT+U
exploration of the energy landscape found a lowest energy nematic state
transforming according to the Eu irreducible representation of the D4h point
group, containing significant inter-orbital nematic components[44]. This form
of nematicity, which breaks also inversion symmetry, was shown to produce
Fermi surfaces containing only a single electron pocket[44]. Whether this
approach has material-specific predictive power is still an open question111
Applying a similar approach to LiFeAs, for example, produces a dramatically
distorted low-energy band structure with missing $d_{xy}$ band at the Fermi
level. R. Valenti, private communication., however, and in any case it remains
important to identify the cause of nematic ordering using more transparent
model-based methods.

These recent developments raise a number of important questions related to the
low-energy band structure and the associated superconducting gap structure of
FeSe. For example, what underlying physical interactions naturally produce
nematic order that generate single-electron-pocket bands[44, 43], and what are
the consequences of this low-energy band structure for our broader
understanding of this fascinating material?

Here, motivated by several earlier theoretical studies of the effects of
longer-range Coulomb interactions $V$ on band structure renormalizations of
iron-based superconductors[38, 46, 47, 48], we explore the nematic phase
spontaneously generated by $V$. Ref. 38 originally investigated these effects
for FeSe, and additionally found a dominant intra-orbital $d$-wave bond
nematic order arising from nearest neighbor Coulomb repulsion. These results
are, however, inconsistent with the current ARPES description of the FeSe
Fermi surface and deserve to be re-examined. To this end, we have begun with
the canonical DFT calculation of the FeSe band structure in the tetragonal
phase[49], and introduced nematic order driven by longer-range Coulomb
interactions systematically. We find that, in addition to the generation of
small electron and hole pockets by $V$, it generates inter-orbital nematicity
between $d_{yz}$ and $d_{xy}$, and $d_{xz}$ and $d_{xy}$ states in a large
region of parameter space at low temperatures. This form of nematic order is
distinct from those discussed previously in the literature[38, 44, 42, 43].
The inter-orbital nematic order components hybridize the low-energy bands near
$X$ and $Y$, respectively, and thereby allow, depending on their amplitudes,
for the lifting of one of the electron pockets. We explore this “$V$-scenario”
for nematicity and Fermi surface anisotropy, and demonstrate how it naturally
generates Fermi surfaces containing only a single electron pocket. We next
discuss the comparison of the spin excitations obtained using this
renormalized band structure with the highly anisotropic spectrum reported in
the experimental literature[5]. Finally, we compute the resulting momentum-
dependent superconducting gap structure, and discuss implications for other
experimental probes of FeSe.

##  II Model and Method

In order to explore how the electronic structure of FeSe is affected by
longer-range Coulomb interactions, specifically nearest-neighbor (NN) density
interactions, we apply the following many-body Hamiltonian

$\displaystyle\begin{split}H=&\,H_{\rm kin}+H_{\rm soc}+H_{\rm int}\\\
=&-\sum_{i\,j}\sum_{\mu\nu}\sum_{\sigma}c^{\,\dagger}_{i\mu\sigma}\left(t^{\,\mu\nu}_{ij}+\mu_{0}\delta_{ij}\delta_{\mu\nu}\right)c^{\,\phantom{\dagger}}_{j\nu\sigma}\\\
&+\lambda_{\rm
SOC}\sum_{i}\sum_{\mu\nu}\sum_{\sigma\sigma^{\prime}}c^{\dagger}_{i\mu\sigma}\,\bm{L}^{\mu\nu}\cdot\bm{S}^{\sigma\sigma^{\prime}}c^{\phantom{\dagger}}_{i\nu\sigma^{\prime}}\\\
&+\frac{V}{2}\sum_{\langle
i,j\rangle}\sum_{\mu\nu}\sum_{\sigma\sigma^{\prime}}n^{\phantom{\dagger}}_{i\mu\sigma}n^{\phantom{\dagger}}_{j\nu\sigma^{\prime}},\end{split}$
(1)

where $\langle i,j\rangle$ indicates the set of NN sites. Here, the operator
$c^{\phantom{\dagger}}_{i\mu\sigma}$ annihilates an electron in orbital
$d_{\mu}$ with spin projection $\sigma$ at lattice site $\bm{R}_{i}$, and
$n_{i\mu\sigma}$ represents the density operator. The kinetic part of the
Hamiltonian[49] $H_{\rm kin}$ contains the hopping matrix elements and the
chemical potential $\mu_{0}$, $H_{\rm soc}$ includes the atomic spin-orbit
coupling (SOC) with strength $\lambda_{\rm SOC}$, and $H_{\rm int}$ describes
the repulsive NN density interaction. For all results presented in the present
work, we adjust $\mu_{0}$ to keep a fixed electron filling of $\langle
n\rangle=6.0$. Note that for simplicity we completely discard the usual onsite
Hubbard Coulomb repulsion $H_{U}$ in the model. As is well-known such
interactions can lead to important band renormalization, magnetism, and
nematicity[4, 50, 3, 51, 52], none of which lift the $Y$ electron pocket.
Here, we assume that $H_{U}$ merely renormalizes the DFT band structure,
though we do not include such effects explicitly, and explore the nematic
instability generated solely from $H_{\rm int}$ containing only NN Coulomb
repulsion.

The symmetry properties of FeSe, and thereby $H$, are governed by the non-
symmorphic space group $P4/nmm$, which consists of eight point group elements
$\\{g\,|\,(0,0,0)\\}$ for $g\in$ D2d, and
$\\{gI\,|\,(\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}},0)\\}$ with $I$ denoting
inversion. In order for the above to constitute a closed group, the product of
space group elements is defined modulo integer lattice translations [53]. For
both analytical tractability and numerical simplicity we perform our study in
the 1-Fe unit cell, for which the space group $P4/nmm$ is lowered to the
symmorphic subgroup D2d. Although one thereby lacks certain symmetry elements
of the original system, it is still possible to capture the violation of
rotational symmetry and the emergence of nematic order. In particular, we need
to focus only on the generators of the $\rm D_{2d}$ group
$g=\\{E,\,S_{4},\,C_{2},\,C^{\prime}_{2},\,\sigma_{d}\\}$, and the associated
five IRs $\Gamma=\\{\rm A_{1},\,A_{2},\,B_{1},\,B_{2},\,E\\}$. Here, the
generator $S_{4}$ refers to the combined operation of $C_{4}\sigma_{xy}$.

We incorporate the effects of NN Coulomb interactions by performing a Hartree-
Fock mean-field (MF) decoupling of $H_{\rm int}$, and introduce the following
homogeneous bond-order fields in wave-vector space [38, 46]

$\displaystyle
N^{\mu\nu}_{\bf{k}\sigma}=-\frac{2}{\mathcal{N}}\sum_{\bf{k}^{\prime}}f_{\bf{k}-\bf{k}^{\prime}}^{\rm
A_{1}}\langle
c^{\dagger}_{\bf{k}^{\prime}\nu\sigma}c^{\phantom{\dagger}}_{\bf{k}^{\prime}\mu\sigma}\rangle,$
(2)

with $\mathcal{N}$ being the number of lattice sites, and $f_{\bf{k}}^{\rm
A_{1}}$ being the form factor of the interaction, which transforms as the IR
$\rm A_{1}$. Specifically for the nearest-neighbor interaction studied here
$f_{\bf{k}}^{\rm A_{1}}\propto\sum_{{\rm d}}\cos(k_{\rm d}\,a_{\rm d})$, where
$\rm d$ labels the set of primitive lattice vectors, and $a_{\rm d}$ is the
lattice constant in the $\rm d$-direction. Note that in the decoupling of
$H_{\rm int}$, the Hartree terms contribute only to a constant energy shift,
which can be readily absorbed in the chemical potential
$\mu_{0}\mapsto\mu_{0}-4V\sum_{\bf{k}^{\prime}\mu\sigma}\langle
n_{\bf{k}^{\prime}\mu\sigma}\rangle/\mathcal{N}$.

Figure 1: Orbitally resolved Fermi surfaces (FSs) and band structures in the
absence (a),(c) and presence (b),(d) of NN interaction effects. The latter
were obtained with $V=0.45\,$eV and $\alpha=2.8$. The band was adopted from
Ref. [49], with $k_{z}=0$ and SOC coupling in the spin $z$-direction,
$\lambda_{\rm SOC}=30\,$meV. In (d) we clearly see the band shifts and gap
openings discussed in the main text, resulting in the highly anisotropic FS
displayed in panel (b). All figures were obtained with $k_{\rm B}T=1\,$meV and
$k_{\rm B}\equiv 1$.

Upon reaching the nematic transition temperature $T_{n}$, the system will
spontaneously violate $S_{4}$ symmetry, and we therefore seek to split the
bond order fields into $S_{4}$ symmetry-preserving and symmetry-breaking
terms, in order to identify the origin of nematicity driven by NN
interactions. The symmetry-preserving terms, denoted by $N_{{\bf k}\sigma,{\rm
br}}^{\mu\nu}$, will in general lead to band renormalizing (br) terms, and is
obtained by averaging over the four cycles of $S_{4}$

$\displaystyle N_{\bf{k}\sigma,{\rm
br}}^{\mu\nu}=\frac{1}{4}\sum_{\ell=0}^{3}\sum_{\mu^{\prime}\nu^{\prime}}\big{[}D(S^{\ell}_{4})\big{]}^{\mu\mu^{\prime}}\,N^{\mu^{\prime}\nu^{\prime}}_{S_{4}^{-\ell}\bf{k}\,\sigma}\,\big{[}D^{\dagger}(S^{\ell}_{4})\big{]}^{\nu^{\prime}\nu},$
(3)

where $D(S_{4}^{\ell})$ is the representation of $S_{4}^{\ell}$ in orbital
space. As a consequence of SOC, one should in general also perform the $S_{4}$
averaging in spin space. However, by exploiting that the fields $N_{{\bf
k}\sigma}^{\mu\nu}$ are diagonal in spin space and the fact that the action of
$S_{4}$ is equivalent to a $\pi/2$-rotation about the spin $z$-axis, we find
that $S_{4}$ leaves the spin projection $\sigma$ of $N_{{\bf
k}\sigma}^{\mu\nu}$ invariant. Momentum dependent band renormalizing terms,
such as $N_{{\bf k}\sigma,{\rm br}}^{\mu\nu}$, have been previously put
forward as potential candidates for explaining the band structure
renormalizations of iron-pnictides and iron-chalcogenides[38, 46, 51].

Here, the main focus is on the $S_{4}$ symmetry breaking (sb) terms leading to
nematic order, which are simply the remainder of the total field
$N_{\bf{k}\sigma,{\rm
sb}}^{\mu\nu}=N_{\bf{k}\sigma}^{\mu\nu}-N_{\bf{k}\sigma,{\rm br}}^{\mu\nu}$.
In order to extract the form factors entering in the symmetry-preserving and
symmetry-breaking fields, we perform the following projection onto the
normalized lattice versions of the basis functions $f_{\bf{k}}^{\gamma}$
belonging to the group $\rm D_{2d}$

$\displaystyle N_{\bf{k}\sigma,{\rm
br/sb}}^{\mu\nu}=\sum_{\gamma}f^{\gamma}_{\bf{k}}N_{\gamma\sigma,{\rm
br/sb}}^{\mu\nu}.$ (4)

As a consequence of the NN interaction, the sum over $\gamma$ is restricted to
include only the lowest order lattice harmonics, i.e. only linear combinations
of $\cos(k_{\rm d}a_{\rm d})$ and $\sin(k_{\rm d}a_{\rm d})$ enter in
$f_{\bf{k}}^{\gamma}$, while a longer-ranging interaction in general allows
for higher order terms. Thus, we arrive at the following MF decoupled
interaction

$\displaystyle H^{\rm MF}_{\rm
int}=V\sum_{{\bf{k}},\gamma}\sum_{\mu\nu}\sum_{\sigma}c^{\,\dagger}_{\bf{k}\mu\sigma}f_{\bf{k}}^{\gamma}\left(N_{\gamma\sigma,{\rm
br}}^{\mu\nu}+\alpha N_{\gamma\sigma,{\rm
sb}}^{\mu\nu}\right)c^{\,\phantom{\dagger}}_{\bf{k}\nu\sigma}.$ (5)

Note that we here introduced the enhancement factor $\alpha$, since the
interaction strength $V$ in the symmetry-preserving and symmetry-breaking
channels in general will be different, since the fields belong to different
IRs of the point group, and can potentially lead to distinct couplings in a
renormalization group procedure. Thus we allow $\alpha\neq 0$ throughout,
similar to earlier works[38, 46].

Lastly, in App. A we elaborate further on the construction of
$N_{\gamma\sigma,{\rm br/sb}}^{\mu\nu}$, display their explicit matrix
structure, and list the basis functions $f_{\bf{k}}^{\gamma}$ used in the
upcoming sections.

##  III Results

Equipped with the above model Hamiltonian and symmetry considerations, we are
now able to calculate the bond order fields $N_{\gamma\sigma,{\rm
br/sb}}^{\mu\nu}$ self-consistently for a fixed electron filling $\langle
n\rangle=6.0$, (for details see App. A). In our calculations we adopt the
kinetic Hamiltonian for FeSe given in Ref. [49], after including SOC in the
spin $z$-direction of bare strength $\lambda_{\rm SOC}=30\,$meV. However, the
exact value of SOC is not qualitatively important for the emergence of inter-
orbital nematic components and the concomitant disappearance of one of the
electron pockets, as found further below. However, SOC is important, within
the current model, for the FS to contain only a single hole pocket at
$\Gamma$. Due to weak inter-layer coupling in the $c$ direction, we restrict
our study to the quasi-2D lattice of a single layer of Fe and Se atoms. The
resulting orbitally-resolved Fermi surface (FS) and band structure in the
absence of NN interactions are displayed in Fig. 1(a) and 1(c), respectively.

The effects of NN interactions, captured by the terms entering in Eq. (5),
give rise to the two aforementioned effects, namely band renormalization and
$S_{4}$ symmetry breaking. Let us momentarily discuss the former,
$N_{\gamma\sigma,{\rm br}}^{\mu\nu}$, by focusing on the low-energy $t_{2g}$
orbitals. For this case, the fields lead to the following two effects: $i)$
collective down-shift of the hole pockets at $\Gamma$ and $M$, and $ii)$ up-
shift of the Dirac points at $X$ and $Y$. The down-shift of the hole pockets
is readily attributed to the self-consistent fields in Fig. 2(a), which couple
to the form factor $f_{\bf{k}}^{\rm\,s}=\cos k_{x}+\cos k_{y}$. Specifically
the fields $N_{{\rm s}\sigma,{\rm br}}^{xz\,xz}=N_{{\rm s}\sigma,{\rm
br}}^{yz\,yz}=-0.22$ lead to a down-shift of $d_{xz}$,$d_{yz}$-dominated bands
at $\Gamma$, while $N^{xy\,xy}_{{\rm s}\sigma,{\rm br}}=0.24$ ensures the
down-shift of $d_{xy}$-dominated bands at $M$. The latter occurs since the
form factor has a relative sign on the two pockets $f_{\Gamma}^{\,\rm
s}=-f_{M}^{\,\rm s}$. Additionally we note that the effects of the $\rm
s$-wave fields on the electron pockets are minimal since $f^{\,\rm
s}_{X}=f^{\,\rm s}_{Y}=0$.

By contrast to the above, the fields presented in Fig. 2(b) are dominant at
the $X$ and $Y$ points due to the $\rm d$-wave form factor $f_{\bf{k}}^{\,\rm
d}=\cos k_{x}-\cos k_{y}$. In fact, the field $N_{{\rm d}\sigma,{\rm
br}}^{xz\,xz}=0.09$ ($N_{{\rm d}\sigma,{\rm br}}^{yz\,yz}=-0.09$) shifts
$d_{xz}$($d_{yz}$)-dominated bands at $Y$ ($X$), pushing the Dirac points up
closer to the Fermi level. Similar to the $\rm s$-wave fields, also here the
symmetry of the form factor, $f^{\,\rm d}_{X}=-f^{\,\rm d}_{Y}$, compensates
the relative sign between the two fields $N^{xz\,xz}_{\rm
d\sigma,br}/N^{yz\,yz}_{\rm d\sigma,br}=-1$. Furthermore, we point out that
the electron pockets acquire a peanut-like shape for appropriate values of
$V$, since the fields $N_{\rm d\sigma,br}^{\mu\nu}$ do not affect the $d_{xy}$
orbitals. Lastly note that the down- (up-) shift of the hole pockets (Dirac
points) leads to smaller $\Gamma$ and $M$ ($X$ and $Y$) pockets. In fact, it
was previously shown that $N_{\bf{k}\sigma,{\rm br}}^{\mu\nu}$ can completely
remove the $M$ pocket in FeSe, while for LiFeAs it was the $\Gamma$ pockets
which were pushed away from the Fermi level [51].

Figure 2: Band renormalizing (br) and $S_{4}$ symmetry breaking (sb) bond-
order fields $N^{\mu\nu}_{\rm\gamma\sigma,br/sb}$, calculated self-
consistently for the parameters and temperature used in Fig. 1(b) and (d), and
fixed electron filling $\langle n\rangle=6.0$, see App. A for further details.
The fields in (a),(c) couple to the form factor $f_{\bf{k}}^{\,\rm s}=\cos
k_{x}+\cos k_{y}$, while the ones displayed in (b),(d) couple to
$f_{\bf{k}}^{\,\rm d}=\cos k_{x}-\cos k_{y}$. The values of all matrix
elements have for clarity been rounded to the second decimal place. For
brevity we only display $\gamma=\\{\rm s,\,d\\}$, but additional plots of the
remaining fields can be found in App. A.

Turning our attention to the symmetry breaking fields, $N_{\gamma\sigma,{\rm
sb}}^{\mu\nu}$, we find that only fields that are coupled to the $\rm d$-wave
form factor play an essential role, see Fig. 2(c) and 2(d). By applying the
same logic as for the band renormalizing terms, we observe that the field
$N_{\rm d\sigma,sb}^{xz\,xz}$ ($N_{\rm d\sigma,sb}^{yz\,yz}$) leads to an up-
(down-) shift of $d_{xz}$($d_{yz}$)-dominated bands at $Y$ ($X$). While this
effect appears similar to the one encountered for $N_{\rm
d\sigma,br}^{\mu\nu}$, we stress that here the Dirac points are shifted in
opposite directions due to the violation of $S_{4}$-symmetry. This reduction
in symmetry also allows for anisotropic inter-orbital $d_{xz}-d_{xy}$ and
$d_{yz}-d_{xy}$ hybridization terms, namely $\alpha N^{xz\,xy}_{\rm
d\sigma,sb}\approx 0.10$ and a vanishing coupling $\alpha N^{yz\,xy}_{\rm
d\sigma,sb}$. It is in fact these particular couplings which lead to the
distinct hybridization gaps at $X$ and $Y$ evident from Fig. 1(d), and they
will, in synergy with all the effects discussed above, result in the highly
anisotropic FS shown in Fig. 1(b), featuring only a single electron pocket at
the $X$ point. Thus, as a function of temperature, a Lifshitz transition
necessarily takes place such that a shrinking $Y$-pocket eventually disappears
at a certain temperature below $T_{n}$. A detailed discussion of the
experimental evidence for such a temperature-induced Lifshitz transition can
be found in Ref. 43.

To gain deeper insight into the emergent nematic order, we now proceed and
study more carefully the IRs of $N_{{\bf k}\sigma,\rm br/sb}^{\mu\nu}$. By
construction, the $S_{4}$ symmetry preserving fields can only consist of terms
transforming as A1 and A2, and, in fact, we find through our self-consistent
calculations that only A1 terms are non-zero in $N_{{\bf k}\sigma,{\rm
br}}^{\mu\nu}$, see App. B for further details.

Figure 3: Leading order parameters $N_{\rm B_{1}}$ and $N_{\rm E}$ versus
temperature. The system becomes nematic at $T_{n}$, and displays a second
phase transition at $T_{n}^{\prime}$. For the latter, $N_{{\rm E}_{x}}$
becomes non-zero, allowing for hybridization terms at $Y$. In this figure, we
used the same parameters as in Figs. 1 and 2.

In contrast, $N_{{\bf k}\sigma,{\rm sb}}^{\mu\nu}$ can contain terms
transforming as the remaining three IRs of the group, i.e. $\rm B_{1}$, $\rm
B_{2}$ and $\rm E$, allowing for the nematic order to arise in any of these
three channels. As we explicitly show in App. B, we find non-zero $S_{4}$
symmetry breaking terms belonging to two distinct IRs, namely B1 and E. This
implies that the system undergoes two consecutive phase transitions upon
lowering of the temperature. In order to quantify this, we focus on the
$t_{2g}$ orbitals, and define the following leading order parameters

$\displaystyle\begin{split}N_{\rm B_{1}}&=(N_{\rm d\sigma,sb}^{xz\,xz}+N_{\rm
d\bar{\sigma},sb}^{xz\,xz}+N_{\rm d\sigma,sb}^{yz\,yz}+N_{\rm
d\bar{\sigma},sb}^{yz\,yz})/4,\\\ N_{\rm E}&=\left([N_{{\rm d}\sigma,{\rm
sb}}^{xz\,xy}+N_{{\rm d}\bar{\sigma},{\rm sb}}^{xz\,xy}]/2,[N_{{\rm
d}\sigma,{\rm sb}}^{yz\,xy}+N_{{\rm d}\bar{\sigma},{\rm
sb}}^{yz\,xy}]/2\right)\\\ &\equiv(N_{{\rm E}_{x}},N_{{\rm
E}_{y}}),\end{split}$ (6)

where $\bar{\sigma}$ is the opposite spin projection of $\sigma$. For details
see App. B. In Fig. 3 we show the values of these order parameters for various
temperatures, and indeed find that $N_{\rm B_{1}}$ becomes non-zero at
$T_{n}$, while only the single component $N_{{\rm E}_{x}}$ condenses at lower
temperatures $T^{\prime}_{n}$. An intra-orbital nematic order thus arises at
$T_{n}$, lowering the point group symmetry from $\rm D_{\rm 2d}$ to $\rm
D_{2}$, and leads to the up- (down-) shift of $d_{xz}$($d_{yz}$)-dominated
bands at $Y$ ($X$) discussed in the previous paragraph. The second transition
at $T_{n}^{\prime}$ further reduces the point group to $\rm C_{2}$ and allows
for anisotropic inter-orbital $d_{xz}-d_{xy}$ or $d_{yz}-d_{xy}$ hybridization
terms, i.e. the effects ultimately leading to the highly anisotropic FS, shown
in Fig. 1(b).

Figure 4: (a) Spin susceptibility as calculated from our model including the
nematic order induced by nearest-neighbor Coulomb interactions, but using a
fully coherent electronic structure, $Z_{\alpha}=1$, $U=0.805\,$eV, $J/U=1/6$.
(b) Spin susceptibility as obtained from a weakly correlated system by
assuming orbitally selective quasiparticle weights[15] (dark red curve)
($U=3.1\,$eV) and successively decreasing the quasiparticle weight for the
$d_{xy}$ orbital to almost negligible quasiparticle contribution (yellow,
$U=3.65\,$eV). (c) Moderate correlation with a small $Z_{xy}=0.27$, but
additionally splitting the quasiparticle weights between the $d_{xz}$ and
$d_{yz}$ orbitals yields a strongly anisotropic susceptibility with no visible
peak at $(0,\pi)$ starting from $Z_{yz}/Z_{xz}\approx 1.3$ ($U$ decreased
slightly to not cross the magnetic instability).

We stress that the effects of the $e_{g}$ orbitals and remaining fields
$N^{\mu\nu}_{\gamma\sigma,{\rm br/sb}}$ not displayed in Fig. 2, are all
incorporated in our calculations, but do not lead to qualitative changes of
the low-energy band structure and FS, and are therefore not explicitly
mentioned in the above discussion. For a complete overview of all fields
$N_{\gamma\sigma,{\rm br/sb}}^{\mu\nu}$, see App. A. Furthermore, we note that
the final nematic band structure, as seen e.g. in Fig. 1(d), obviously depends
on the starting point, i.e. the tetragonal DFT band structure. Thus, the final
quantitative energy scales, i.e. the hybridization gap at $Y$ and the required
amplitudes of $V$ and $\alpha$ for generating a FS similar to that shown in
Fig. 1(b), depend on the initial bare band structure. Lastly, we note that the
purely intra-orbital nematic order generated from NN Coulomb repulsion found
earlier [38, 46], can be reproduced here when applying the same band structure
as in Ref. 38. We have not located the exact band property that leads to the
additional inter-orbital nematic fields from the DFT bulk FeSe band model used
in the current study[49], but note that the solution presented in Fig. 1 is
quite generic at low temperatures and exists for a wide parameter range.

The low-energy electronic structure established above has important
consequences for spin excitations that should be compared with experiments. In
this respect, a series of inelastic neutron scattering experiments[54, 55, 5]
have established a remarkable set of magnetic phenomena in FeSe. At low
temperatures in the nematic phase, but high energies $\sim$ 100 meV, strong
$(\pi,\pi)$ (Néel) fluctuations dominate the spectrum, with weaker but still
prominent $(\pi,0)$ and $(0,\pi)$ (stripe-like) fluctuations. As the energy is
lowered to a few tens of meV, a spin gap develops in the $(\pi,\pi)$ spectrum,
but $(\pi,0)$ spin fluctuations strengthen. Notably, the only measurement of
detwinned  FeSe crystal finds that at low energies the intensity of the
$(0,\pi)$ fluctuations essentially vanishes[5]. This extraordinary result
should emerge from a proper theory of low energy spin and orbital degrees of
freedom. We show now that, counter-intuitively, the current proposal for low-
energy nematic electronic structure is not sufficient to explain the above-
mentioned magnetic properties, and requires the additional physics of orbital
selective correlations.

In Fig. 4(a), we first illustrate the bare magnetic susceptibility
Re$\chi_{0}({\bf q},\omega=0)$, together with the enhanced susceptibility
obtained in the random phase approximation $\chi^{RPA}$. The spectrum is
clearly dominated by intense $(\pi,\pi)$ fluctuations arising from scattering
between the $d_{xy}$ states, which are not observed in experiment.
Furthermore, $(\pi,0)$ and $(0,\pi)$ states are nearly degenerate, in
contradiction to the results of Ref. 5. Note that these issues are also common
to conventional spin-fluctuation theories of FeSe including a $Y$ pocket[23,
56], or other novel schemes to lift the $Y$ pocket[43] via nematic order. In
Ref. 23, it was proposed that they could be resolved by assuming orbitally
selective incoherence of $d_{xy}$, $d_{xz}$, and $d_{yz}$ states, such as
should take place according to previous theory[22] and as observed in some
experiments[20]. With a phenomenological insertion of orbitally dependent
quasiparticle weights $Z_{\alpha}$, with $\alpha$ an orbital index, Kreisel et
al.[23] could fit inelastic neutron data on FeSe using very strong suppression
of $d_{xy}$ weights and somewhat smaller suppression of $d_{xz}$ and $d_{yz}$
weights, assuming also a large ratio of at least 1.7 for the ratio
$Z_{yz}/Z_{xz}$.

We therefore explore how orbital incoherence could improve the agreement of
the susceptibility calculated from the current highly nematic electronic
structure with experiment. In Fig. 4(b), we present the results for weak
correlations, and the evolution of the susceptibility as correlations are
enhanced by a substantial suppression of the $d_{xy}$ quasiparticle weight
(see Appendix C for a brief discussion of the Ansatz of Kreisel et al.[23]).
As anticipated from the contribution of the orbitally resolved susceptibility
of the $d_{xy}$ orbital (see Fig. 4(a)), with reduction of $Z_{xy}$ comes the
suppression of the Néel peak and concomitant moderate enhancement of the
$(\pi,0)$ stripe peak. On the other hand the $(\pi,0)/(0,\pi)$ anisotropy is
still much weaker than that reported in Ref. 5. In Fig. 4(c), we therefore
show the effect of additionally increasing the $Z_{yz}/Z_{xz}$ anisotropy. It
is easy to see that much larger $(\pi,0)/(0,\pi)$ anisotropies are obtained in
this case, but also that substantially smaller quasiparticle weight ratios, of
order $\sim 1.3$, are required compared to Ref. 23, due to the effect of
inter-orbital nematic order introduced here.

Figure 5: (a) Superconducting gap symmetry function $g(\bf{k})$ as calculated
from a fully coherent electronic structure showing the strongly anisotropic
gap on the two pockets. The expected spectrum shows nodal features and the
eigenvalue $\lambda$ is sizeable but small. (b) $g(\bf{k})$ calculated from a
electronic structure with very incoherent $d_{xy}$ orbital; almost no effect
is visible except that the eigenvalue can be larger since the relative
contribution of the $(\pi,0)$ scattering is higher. (c) Same quantity, but
calculated including a moderately correlated $d_{xy}$ orbital and a nematic
splitting of the quasiparticle weights in the $d_{yz}$ and $d_{xz}$ orbitals,
yielding an order parameter with a tiny true gap, in contrast to nodes in (a)
and (b).

On the other hand, these various hypothetical renormalizations do not lead to
substantial changes in the gap structure obtained for FeSe within the
corresponding spin fluctuation pairing theory. This is simply because in a
multiband system even at $\omega=0$ contributions to $\chi(\bf q)$ arise from
states tens or even hundreds of meV from the Fermi level[57]. By contrast,
Fermi surface states determine the anisotropy of the effective pairing
interaction completely. This effect can be seen easily in Fig. 5, where we
plot the leading eigenvector of the linearized gap equation (see Appendix C)
for the same three cases described in Fig. 4. The overall structure of the
gaps and the density of states are seen to be virtually identical.

##  IV Discussion and conclusions

In this work we have pursued a scenario where nematicity originates entirely
from NN Coulomb repulsion, even though it is well-known that onsite
interactions alone may also drive a nematic instability as a precursor to
stripe magnetism[4]. Interestingly, orbitally resolved studies within the
spin-nematic onsite interaction-only scenario, do find sizable inter-orbital
nematic susceptibilities[58]. Such studies, however, have only been performed
in the tetragonal paramagnetic phase, and it remains to be seen whether
spontaneous breaking in the inter-orbital channel, similar to the current
proposal, can also arise from higher order processes solely from onsite
interactions. However, the absence of magnetism in FeSe and the overall
agreement of the here-presented results to the experimental facts, raises the
question of whether indeed NN-repulsion is the generator of nematicity for
FeSe. From cRPA[59] and DMFT[60] calculations it is known that interactions
are generally larger in FeSe than any of the other iron-based superconducting
systems, presumably because a lack of intervening spacer layers reduces the
screening. Onsite interactions $U$, $J$ are known to lead to band
renormalizations, in particular Fermi surface pocket shrinkage[47, 48] and
orbital-dependent band narrowing. It is tempting to speculate that this, in
turn, effectively enhances the importance of the unusually large $V$ in FeSe,
thereby boosting instabilities driven by this channel.

Recently, another theoretical study investigated the possibility of nematic
order lifting one of the electron pockets (above the Fermi level) in FeSe[43].
In agreement with the present work, a Lifshitz transition necessarily takes
place as a function of temperature, and a superconducting gap structure
consistent with experiments follows directly from the resulting FS with the
missing $Y$-pocket. A main differences between Ref. 43 and the current
approach is the nature of the nematic order; whereas Rhodes et al.[43] begin
with a phenomenological $\bf k\cdot\bf p$ expansion around $\Gamma$, $X$ and
$Y$, and explore the role of a large intra-orbital B1g $d_{xy}$-nematic order
parameter imposed by “hand” on the band structure, together with additional
onsite Hartree shifts to this orbital, we have started from a nearest-neighbor
Coulomb interaction, and shown that (self-consistently generated) inter-
orbital $d_{xy}-d_{xz/yz}$ nematic components lift the $Y$-pocket. Thus, while
a B1g intra-orbital $d_{xy}$ nematic component is also present in our approach
as seen from Fig. 2(d), the most important nematic components for generating a
FS containing no electron $Y$-pocket are the distinct inter-orbital
components.

As mentioned in the “Results” section, the nematic lowest-temperature phase
advocated in this work, resides in the $\rm C_{2}$ (monoclinic) group,
containing only a $C_{2}$ rotation around the $x$-axis as the remaining
symmetry operation. This constitutes a clear prediction within the current
scenario; the electronic sector should exhibit a double transition as the
temperature is lowered. We are not aware of such evidence, e.g. from specific
heat data[61, 62, 63], which may be because of the very small electronic
entropy change at the second transition, or simply because the two transitions
are accidentally close (in temperature) (see App. B for more details) in FeSe.
In addition, if the nematic order couples strongly enough to the atomic
lattice, such symmetry lowering could be tested by experiments sensitive to
the overall crystal point group. In this regard, however, we note that recent
experimental and theoretical studies have advocated for a more complex static
crystal microstructure in FeSe than previously thought. Locally, FeSe appears
to accommodate a myriad of inhomogeneous nanoscale lattice distortions, where
only the spatially averaged structure complies with the standard tetragonal
(orthorhombic) crystal symmetries at high (low) temperatures[64, 65, 44, 66].
This may in fact be the result of two nematic channels very close in in
energy. Lastly, it is tempting to speculate that the inter-orbital nematic
order discussed here, may be relevant for the recent experimental studies of
structural transitions in Bi1-xSrxNi2As2, which exhibit transitions from a
high-temperature tetragonal phase to a triclinic low-temperature phase[67].

Recent theoretical works explored the possibility of pinned local nematic
order[68, 69, 70]. In particular, Ref. 70 explored the local disorder-induced
nematicity from nonmagnetic impurities in the tetragonal phase at $T>T_{n}$.
Several experimental works have reported evidence for such local nematic order
above $T_{n}$[64, 65]. The results of Ref. 70 were obtained by applying a one-
band model and relied on interactions leading to a single-component $\rm
B_{\rm 1g}$ nematic order, and hence the question arises how local pinned
nematic order gets affected by the presence of the substantial inter-orbital
components found in this work? We have answered this question by performing a
real-space calculation similar to that of Ref. 70, but generalized to the
multi-orbital case. For the intra-orbital $\rm B_{1}$ order parameter $N_{\rm
B_{1}}$ studied here, we observe that it enters the Landau free energy
expansion in the exact same manner as the order parameter studied in Ref. 70.
Therefore the local impurity-induced order exhibits a ”flower-shape” pinned by
nonmagnetic impurities[70]. However, for the lower temperature phase, the
spatially dependent local nematic impurity-induced order is distinct from that
of Ref. 70, since the order parameter $N_{E}$ couples differently to the
nonmagnetic impurity. It remains an interesting future study to investigate
the detailed experimental consequences of these various local nematic orders.

In summary, we have discovered an inter-orbital nematic order generated from
nearest-neighbor Coulomb repulsion for electronic models relevant for FeSe in
particular, and perhaps the iron-based systems more broadly. A natural
property of this kind of nematic order is the generation of highly anisotropic
Fermi surfaces, featuring in some cases, only a single hole and electron
pocket. We have shown how, for FeSe, this explains photoemission data and the
experimentally extracted very anisotropic superconducting gap structure.
However, a consistent picture of the neutron scattering data, cannot be
straightforwardly obtained within this nematic scenario, without invoking
additional self-energy effects, which in the simplest case, involves sizable
corrections in the form of reduced orbitally dependent quasiparticle weights.

## Acknowledgements

We thank M. H. Christensen and I. Eremin for insightful conversations. D.S.
and B.M.A. acknowledge support from the Carlsberg Foundation. B. M. A.
acknowledges support from the Independent Research Fund Denmark grant number
8021-00047B. P. J. H. was supported by the U.S. Department of Energy under
Grant No. DE-FG02-05ER46236.

##  Appendix A Details on bond-order fields and self-consistent calculations

The explicit forms of $N_{\bf{k}\sigma,{\rm br}}^{\mu\nu}$ and
$N_{\bf{k}\sigma,{\rm sb}}^{\mu\nu}$ are found through the averaging in Eq.
(3), and the difference $N_{\bf{k}\sigma,{\rm
sb}}^{\mu\nu}=N_{\bf{k}\sigma}^{\mu\nu}-N_{\bf{k}\sigma,{\rm br}}^{\mu\nu}$,
respectively. Yet this procedure can be further simplified, by relying on the
projection of the fields onto the normalized lattice versions of the basis
functions $f_{\bf{k}}^{\gamma}$

$\displaystyle
N_{\bf{k}\sigma}^{\mu\nu}=\sum_{\gamma}f_{\bf{k}}^{\gamma}\,N_{\gamma\sigma}^{\mu\nu}.$
(7)

This projection disentangles the wave-vector and orbital dependence, and
thereby greatly simplifies Eq. (3). For NN interactions in 2D systems, which
are left invariant under the elements of $\rm D_{2d}$, the set of basis
functions entering in Eq. (7) are

$\displaystyle f_{\bf{k}}^{\,\rm s}=\cos k_{x}+\cos k_{y},$ $\displaystyle
f_{\bf{k}}^{\,\rm d}=\cos k_{x}-\cos k_{y},$ (8a) $\displaystyle
f_{\bf{k}}^{{\,\rm p}_{x}}=\sqrt{2}\,i\sin k_{x},$ $\displaystyle
f_{\bf{k}}^{\,{\rm p}_{y}}=\sqrt{2}\,i\sin k_{y},$ (8b)

with $a_{x}=a_{y}=a\equiv 1$. For NN interactions in three dimensions, the
above equations are accompanied by $f^{\,\rm s}_{k_{z}}=\sqrt{2}\cos(k_{z}c)$
and $f^{\,\rm d}_{k_{z}}=\sqrt{2}\,i\sin(k_{z}c)$, with $a_{z}\equiv c\neq a$.
Each basis function transforms according to one of the IRs $\Gamma$ of the
point group. Specifically, $f^{\,\rm s}_{\bf k}$ ($f^{\,\rm d}_{\bf k}$)
transforms as $\rm A_{1}$ ($\rm B_{1}$), while $(f_{\bf k}^{\,{\rm
p}_{x}},f_{\bf k}^{\,{\rm p}_{y}})$ transform jointly as the 2D IR E.

Next step in determining $N^{\mu\nu}_{{\bf k}\sigma,{\rm br/sb}}$, is to
represent $S_{4}$ in spin space and in the relevant orbital basis
$\\{d_{xz},\,d_{yz},\,d_{x^{2}-y^{2}},\,d_{xy},\,d_{3z^{2}-r^{2}}\\}$. The
former is needed since the SOC in Eq. (1) breaks spin-rotation symmetry.
Additionally, we also need to determine how $S_{4}^{-1}$ acts on the wave-
vectors. Straightforwardly, we find

$\displaystyle\begin{split}\mathcal{D}(S_{4})&=D(S_{4})\otimes D_{\rm
spin}(S_{4})\\\ &=\begin{pmatrix}0&1&0&0&0\\\ -1&0&0&0&0\\\ 0&0&-1&0&0\\\
0&0&0&-1&0\\\
0&0&0&0&1\end{pmatrix}\otimes\frac{\mathds{1}_{\sigma}-i\sigma_{z}}{\sqrt{2}},\\\
S^{-1}_{4}\bf{k}&=(k_{y},-k_{x},-k_{z}),\end{split}$ (9)

where $D(S_{4})$ and $D_{\rm spin}(S_{4})$ are matrix representations in
orbital and spin space, respectively, while $\mathcal{D}(S_{4})$ is the
representation in the combined space. Any symmetry operation acting in spin
space can be expressed in terms of the Pauli matrices $\sigma_{x,y,z}$
accompanied by the identity $\mathds{1}_{\sigma}$. Note, however, that the
fields $N_{{\bf k}\sigma}^{\mu\nu}$ are diagonal in spin space, and therefore
not affected by $D_{\rm spin}(S_{4})$. One can therefore solely consider the
action of $D(S_{4})$, as discussed in connection to Eq. (3). In general, when
considering a given symmetry element $g$ of the point group $\rm D_{\rm 2d}$,
one needs to apply the full representation $\mathcal{D}(g)$.

By executing the above on a 2D system, we end up with the following band
renormalizing (br) terms

$\displaystyle N_{\rm s\sigma,br}$
$\displaystyle=\begin{pmatrix}\frac{1}{2}(N^{11}_{{\rm s}\sigma}+N^{22}_{{\rm
s}\sigma})&\frac{1}{2}(N^{12}_{{\rm s}\sigma}-N^{21}_{{\rm
s}\sigma})&0&0&0\\\\[5.69046pt] \frac{1}{2}(N^{21}_{{\rm
s}\sigma}-N^{12}_{{\rm s}\sigma})&\frac{1}{2}(N^{11}_{{\rm
s}\sigma}+N^{22}_{{\rm s}\sigma})&0&0&0\\\\[5.69046pt] 0&0&N^{33}_{{\rm
s}\sigma}&N^{34}_{{\rm s}\sigma}&0\\\\[5.69046pt] 0&0&N^{43}_{{\rm
s}\sigma}&N^{44}_{{\rm s}\sigma}&0\\\\[5.69046pt] 0&0&0&0&N^{55}_{{\rm
s}\sigma}\\\\[5.69046pt] \end{pmatrix},$ (10a) $\displaystyle N_{\rm
d\sigma,{\rm br}}$ $\displaystyle=\begin{pmatrix}\frac{1}{2}(N^{11}_{{\rm
d}\sigma}-N^{22}_{{\rm d}\sigma})&\frac{1}{2}(N^{12}_{{\rm
d}\sigma}+N^{21}_{{\rm d}\sigma})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(N^{12}_{{\rm d}\sigma}+N^{21}_{{\rm
d}\sigma})&-\frac{1}{2}(N^{11}_{{\rm d}\sigma}-N^{22}_{{\rm
d}\sigma})&0&0&0\\\\[5.69046pt] 0&0&0&0&N_{{\rm d}\sigma}^{35}\\\\[5.69046pt]
0&0&0&0&N_{{\rm d}\sigma}^{45}\\\\[5.69046pt] 0&0&N^{53}_{{\rm
d}\sigma}&N_{{\rm d}\sigma}^{54}&0\end{pmatrix},$ (10b) $\displaystyle N_{{\rm
p}_{x}\sigma,{\rm br}}$ $\displaystyle=\begin{pmatrix}0&0&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{13}+N_{{\rm p}_{y}\sigma}^{23})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{14}+N_{{\rm p}_{y}\sigma}^{24})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{15}-N_{{\rm p}_{y}\sigma}^{25})\\\\[5.69046pt]
0&0&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{23}-N_{{\rm
p}_{y}\sigma}^{13})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{24}-N_{{\rm
p}_{y}\sigma}^{14})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{25}+N_{{\rm
p}_{y}\sigma}^{15})\\\\[5.69046pt] \frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{31}+N_{{\rm p}_{y}\sigma}^{32})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{32}-N_{{\rm p}_{y}\sigma}^{31})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{41}+N_{{\rm
p}_{y}\sigma}^{42})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{42}-N_{{\rm
p}_{y}\sigma}^{41})&0&0&0\\\\[5.69046pt] \frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{51}-N_{{\rm p}_{y}\sigma}^{52})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{52}+N_{{\rm p}_{y}\sigma}^{51})&0&0&0\end{pmatrix},$ (10c)
$\displaystyle N_{{\rm p}_{y}\sigma,{\rm br}}$
$\displaystyle=\begin{pmatrix}0&0&-\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{23}-N_{{\rm p}_{y}\sigma}^{13})&-\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{24}-N_{{\rm p}_{y}\sigma}^{14})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{25}+N_{{\rm p}_{y}\sigma}^{15})\\\\[5.69046pt]
0&0&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{13}+N_{{\rm
p}_{y}\sigma}^{23})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{14}+N_{{\rm
p}_{y}\sigma}^{24})&-\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{15}-N_{{\rm
p}_{y}\sigma}^{25})\\\\[5.69046pt] -\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{32}-N_{{\rm p}_{y}\sigma}^{31})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{31}+N_{{\rm p}_{y}\sigma}^{32})&0&0&0\\\\[5.69046pt]
-\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{42}-N_{{\rm
p}_{y}\sigma}^{41})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{41}+N_{{\rm
p}_{y}\sigma}^{42})&0&0&0\\\\[5.69046pt] \frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{52}+N_{{\rm p}_{y}\sigma}^{51})&-\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{51}-N_{{\rm p}_{y}\sigma}^{52})&0&0&0\end{pmatrix},$ (10d)

where we used the convenient shorthand notation
$\\{d_{xz},\,d_{yz},\,d_{x^{2}-y^{2}},\,d_{xy},\,d_{3z^{2}-r^{2}}\\}\equiv\\{1,\,2,\,3,\,4,\,5\\}$.
Similarly we find the following symmetry breaking (sb) terms

$\displaystyle N_{\rm s\sigma,sb}$
$\displaystyle=\begin{pmatrix}\frac{1}{2}(N^{11}_{{\rm s}\sigma}-N^{22}_{{\rm
s}\sigma})&\frac{1}{2}(N^{12}_{{\rm s}\sigma}+N^{21}_{{\rm s}\sigma})&N_{\rm
s\sigma}^{13}&N_{\rm s\sigma}^{14}&N_{\rm s\sigma}^{15}\\\\[5.69046pt]
\frac{1}{2}(N^{12}_{{\rm s}\sigma}+N^{21}_{{\rm
s}\sigma})&-\frac{1}{2}(N^{11}_{{\rm s}\sigma}-N^{22}_{{\rm s}\sigma})&N_{\rm
s\sigma}^{23}&N_{\rm s\sigma}^{24}&N_{\rm s\sigma}^{25}\\\\[5.69046pt] N_{\rm
s\sigma}^{31}&N_{\rm s\sigma}^{32}&0&0&N_{\rm s\sigma}^{35}\\\\[5.69046pt]
N_{\rm s\sigma}^{41}&N_{\rm s\sigma}^{42}&0&0&N_{\rm
s\sigma}^{45}\\\\[5.69046pt] N_{\rm s\sigma}^{51}&N_{\rm s\sigma}^{52}&N_{\rm
s\sigma}^{53}&N_{\rm s\sigma}^{54}&0\\\\[5.69046pt] \end{pmatrix},$ (11a)
$\displaystyle N_{\rm d\sigma,{\rm sb}}$
$\displaystyle=\begin{pmatrix}\frac{1}{2}(N^{11}_{{\rm d}\sigma}+N^{22}_{{\rm
d}\sigma})&\frac{1}{2}(N^{12}_{{\rm d}\sigma}-N^{21}_{{\rm
d}\sigma})&N^{13}_{{\rm d}\sigma}&N^{14}_{{\rm d}\sigma}&N^{15}_{{\rm
d}\sigma}\\\\[5.69046pt] \frac{1}{2}(N^{21}_{{\rm d}\sigma}-N^{12}_{{\rm
d}\sigma})&\frac{1}{2}(N^{11}_{{\rm d}\sigma}+N^{22}_{{\rm
d}\sigma})&N^{23}_{{\rm d}\sigma}&N^{24}_{{\rm d}\sigma}&N^{25}_{{\rm
d}\sigma}\\\\[5.69046pt] N^{31}_{{\rm d}\sigma}&N^{32}_{{\rm
d}\sigma}&N^{33}_{{\rm d}\sigma}&N^{34}_{{\rm d}\sigma}&0\\\\[5.69046pt]
N^{41}_{{\rm d}\sigma}&N^{42}_{{\rm d}\sigma}&N^{43}_{{\rm
d}\sigma}&N^{44}_{{\rm d}\sigma}&0\\\\[5.69046pt] N^{51}_{{\rm
d}\sigma}&N^{52}_{{\rm d}\sigma}&0&0&N^{55}_{{\rm d}\sigma}\end{pmatrix},$
(11b) $\displaystyle N_{{\rm p}_{x}\sigma,{\rm sb}}$
$\displaystyle=\begin{pmatrix}N^{11}_{{\rm p}_{x}\sigma}&N^{12}_{{\rm
p}_{x}\sigma}&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{13}-N_{{\rm
p}_{y}\sigma}^{23})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{14}-N_{{\rm
p}_{y}\sigma}^{24})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{15}+N_{{\rm
p}_{y}\sigma}^{25})\\\\[5.69046pt] N^{21}_{{\rm p}_{x}\sigma}&N^{22}_{{\rm
p}_{x}\sigma}&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{23}+N_{{\rm
p}_{y}\sigma}^{13})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{24}+N_{{\rm
p}_{y}\sigma}^{14})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{25}-N_{{\rm
p}_{y}\sigma}^{15})\\\\[5.69046pt] \frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{31}-N_{{\rm p}_{y}\sigma}^{32})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{32}+N_{{\rm p}_{y}\sigma}^{31})&N^{33}_{{\rm
p}_{x}\sigma}&N^{34}_{{\rm p}_{x}\sigma}&N^{35}_{{\rm
p}_{x}\sigma}\\\\[5.69046pt] \frac{1}{2}(N_{{\rm p}_{x}\sigma}^{41}-N_{{\rm
p}_{y}\sigma}^{42})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{42}+N_{{\rm
p}_{y}\sigma}^{41})&N^{43}_{{\rm p}_{x}\sigma}&N_{{\rm
p}_{x}\sigma}^{44}&N^{45}_{{\rm p}_{x}\sigma}\\\\[5.69046pt]
\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{51}+N_{{\rm
p}_{y}\sigma}^{52})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{52}-N_{{\rm
p}_{y}\sigma}^{51})&N_{{\rm p}_{x}\sigma}^{53}&N^{54}_{{\rm
p}_{x}\sigma}&N_{{\rm p}_{x}\sigma}^{55}\end{pmatrix},$ (11c) $\displaystyle
N_{{\rm p}_{y}\sigma,{\rm sb}}$ $\displaystyle=\begin{pmatrix}N^{11}_{{\rm
p}_{y}\sigma}&N^{12}_{{\rm p}_{y}\sigma}&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{23}+N_{{\rm p}_{y}\sigma}^{13})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{24}+N_{{\rm p}_{y}\sigma}^{14})&-\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{25}-N_{{\rm p}_{y}\sigma}^{15})\\\\[5.69046pt] N^{21}_{{\rm
p}_{y}\sigma}&N^{22}_{{\rm p}_{y}\sigma}&-\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{13}-N_{{\rm p}_{y}\sigma}^{23})&-\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{14}-N_{{\rm p}_{y}\sigma}^{24})&\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{15}+N_{{\rm p}_{y}\sigma}^{25})\\\\[5.69046pt]
\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{32}+N_{{\rm
p}_{y}\sigma}^{31})&-\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{31}-N_{{\rm
p}_{y}\sigma}^{32})&N^{33}_{{\rm p}_{y}\sigma}&N^{34}_{{\rm
p}_{y}\sigma}&N^{35}_{{\rm p}_{y}\sigma}\\\\[5.69046pt] \frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{42}+N_{{\rm p}_{y}\sigma}^{41})&-\frac{1}{2}(N_{{\rm
p}_{x}\sigma}^{41}-N_{{\rm p}_{y}\sigma}^{42})&N^{43}_{{\rm
p}_{y}\sigma}&N_{{\rm p}_{y}\sigma}^{44}&N^{45}_{{\rm
p}_{y}\sigma}\\\\[5.69046pt] -\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{52}-N_{{\rm
p}_{y}\sigma}^{51})&\frac{1}{2}(N_{{\rm p}_{x}\sigma}^{51}+N_{{\rm
p}_{y}\sigma}^{52})&N_{{\rm p}_{y}\sigma}^{53}&N^{54}_{{\rm
p}_{y}\sigma}&N_{{\rm p}_{y}\sigma}^{55}\end{pmatrix}.$ (11d)

These fields, combined with the basis functions $f_{\bf{k}}^{\gamma}$, enter
in the mean-field decoupled Hamiltonian in the following way, see also Eq. (5)

$\displaystyle\begin{split}H_{\rm int}\approx H_{\rm int}^{\rm
MF}=V\sum_{\bf{k}}\sum_{\mu\nu}\sum_{\sigma}c^{\,\dagger}_{\bf{k}\mu\sigma}\left[f_{\bf{k}}^{\rm\,s}\left(N^{\mu\nu}_{{\rm
s}\sigma,{\rm br}}+\alpha N^{\mu\nu}_{{\rm s}\sigma,{\rm
sb}}\right)+f_{\bf{k}}^{\rm\,d}\left(N^{\mu\nu}_{{\rm d}\sigma,{\rm
br}}+\alpha N^{\mu\nu}_{{\rm d}\sigma,{\rm
sb}}\right)\right.\quad\quad\quad\quad\\\ +\left.f_{\bf{k}}^{\,{\rm
p}_{x}}\left(N^{\mu\nu}_{{\rm p}_{x}\sigma,{\rm br}}+\alpha N^{\mu\nu}_{{\rm
p}_{x}\sigma,{\rm sb}}\right)+f_{\bf{k}}^{\,{\rm p}_{y}}\left(N^{\mu\nu}_{{\rm
p}_{y}\sigma,{\rm br}}+\alpha N^{\mu\nu}_{{\rm p}_{y}\sigma,{\rm
sb}}\right)\right]c^{\,\phantom{\dagger}}_{\bf{k}\nu\sigma}.\end{split}$ (12)

Figure 6: Band renormalizing bond-order fields, $N^{\mu\nu}_{\gamma\sigma,{\rm
br}}$, calculated self-consistently with the parameters used in Fig 1(b) and
1(d). A given field $N^{\mu\nu}_{\gamma\sigma,{\rm br}}$ couples to the
appropriate form factor $f_{{\bf k}}^{\gamma}$, see Eq. (8). (a) and (b) are
discussed in the main text, while (c) and (d) only lead to minimal effects on
the low-energy band structure. All numbers appearing in the matrices have for
clarity been rounded to the second decimal place. Figure 7: Symmetry breaking
fields, $N_{\gamma\sigma,{\rm sb}}^{\mu\nu}$, calculated self-consistently
with the parameters used in Fig 1(b) and 1(d). Each field is coupled to the
belonging form factor listed in Eq. (8). Only (b) and (d) display non-zero
solutions, and the former will ultimately give rise to the highly anisotropic
FS shown in Fig. 1(b). Similar to Figs. 2 and 6, also here we rounded the
values of the matrix elements to the second decimal place.

The Hamiltonian becomes bilinear in creation and annihilation operators upon
approximating $H_{\rm int}\approx H_{\rm int}^{\rm MF}$, and we can thus
easily express the fermionic operators in the diagonal basis of the
Hamiltonian $c_{{\bf k}\mu\sigma}^{\,\phantom{\dagger}}=\sum_{m}\gamma_{{\bf
k}m\sigma}^{\,\phantom{\dagger}}\langle{\bf k}\,\mu\,\sigma|{\bf
k}\,m\,\sigma\rangle\equiv\sum_{m}\gamma_{{\bf
k}m\sigma}^{\,\phantom{\dagger}}a^{\mu\sigma}_{m}({\bf k})$, where
$\gamma^{\phantom{\dagger}}_{{\bf k}m\sigma}$ are the operators related to the
eigenstates of the Hamiltonian $H\gamma^{\dagger}_{{\bf
k}m\sigma}|0\rangle=E_{{\bf k}m\sigma}|{\bf k}\,m\,\sigma\rangle$. In this
diagonal basis the electron density $\langle n\rangle$ takes the simple form

$\displaystyle\langle n\rangle$ $\displaystyle=\frac{1}{\mathcal{N}}\sum_{{\bf
k}}\sum_{\sigma}\sum_{m}n_{\rm F}(E_{{\bf k}m\sigma}),$ (13)

where $n_{\rm F}(E_{{\bf k}m\sigma})$ is the Fermi-Dirac distribution
function. Furthermore, we find that the terms entering in the matrix elements
in Eqs. (10) and (11) are expressed as

$\displaystyle\begin{split}N^{\mu\nu}_{\gamma\sigma}&=-\frac{1}{\mathcal{N}}\sum_{{\bf
k}^{\prime}}\left[f^{\gamma}_{{\bf k}^{\prime}}\right]^{*}\langle
c^{\,\dagger}_{{\bf k}^{\prime}\nu\sigma}c^{\,\phantom{\dagger}}_{{\bf
k}^{\prime}\mu\sigma}\rangle\\\ &=-\frac{1}{\mathcal{N}}\sum_{{\bf
k}^{\prime}}\sum_{m}\big{[}f^{\gamma}_{{\bf
k}^{\prime}}\big{]}^{*}\,\big{[}a^{\nu\sigma}_{m}\big{]}^{*}a^{\mu\sigma}_{m}\,n_{\rm
F}(E_{{\bf k}^{\prime}m\sigma}).\end{split}$ (14)

By calculating the above self-consistently with the parameters used in Fig.
1(b) and 1(c), we arrive at the band renormalizing and $S_{4}$ symmetry
breaking fields displayed in Fig. 6 and Fig. 7, respectively. Only fields
coupled to $f_{{\bf k}}^{\,\rm s}$ and $f_{{\bf k}}^{\,\rm d}$ are discussed
in the main text, since the remaining ones lead to minimal effects on the low-
energy $t_{2g}$ orbitals. This is obviously true for $N_{{\rm
p}_{x}\sigma,{\rm sb}}^{\mu\nu}$ since all these are zero, while for the
remaining band renormalizing terms and $N^{\mu\nu}_{{\rm p}_{y}\sigma,{\rm
sb}}$ it is the form factors which are eliminating the effects, since
$f^{\,{\rm p}_{x}}_{{\bf k}=(0,k_{y})}=f^{\,{\rm p}_{y}}_{{\bf
k}=(k_{x},0)}=0$. This can be seen by considering the field $N^{xz\,xy}_{{\rm
p}_{x}\sigma,{\rm br}}=-0.27$ ($N^{yz\,xy}_{{\rm p}_{y}\sigma,{\rm
br}}=-0.27$) which introduces inter-orbital $d_{xz}-d_{xy}$ ($d_{yz}-d_{xy}$)
hybridization terms, however, the form factor for this field goes to zero on
the orbitally-relevant electron pocket at $Y$ ($X$), thus rendering the effect
minimal. Same argument holds for $N_{{\rm p}_{y}\sigma,{\rm sb}}^{xz\,yz}$,
which should be relevant at $\Gamma$, however $f^{\,{\rm p}_{y}}_{\Gamma}=0$.

##  Appendix B Irreducible representations of bond-order fields

Although the above classification of band renormalizing and symmetry breaking
fields suffices in describing the occurrence of nematic order, it fails to
express the exact symmetry of the emergent nematic order. In other words, the
band renormalizing and symmetry breaking fields transform as a sum of IRs,
specifically

$\displaystyle N^{\mu\nu}_{{\bf k}\sigma,\rm br}$ $\displaystyle\sim{\rm
A_{1}}\oplus{\rm A_{2}},$ $\displaystyle N^{\mu\nu}_{{\bf k}\sigma,\rm sb}$
$\displaystyle\sim{\rm B_{1}}\oplus{\rm B_{2}}\oplus{\rm E},$ (15)

and it is therefore not obvious whether the nematic order parameter transforms
as $\rm B_{1}$, $\rm B_{2}$ or $\rm E$.

In order to shed light on this ambiguity, we will in the following perform a
thorough classification of the bond-order fields, and further segregate these
into IRs. In doing so, we average out the various IRs of the already
established $N_{\gamma\sigma,{\rm br/sb}}^{\mu\nu}$, similar to Eq. (3), in
the following way

$\displaystyle N_{\gamma,{\rm A_{1}}}$
$\displaystyle=\frac{1}{2}\sum_{\ell=0}^{1}\mathcal{D}(C_{2}^{\prime\ell})\,N_{\gamma,{\rm
br}}\,\mathcal{D}^{\dagger}(C_{2}^{\prime\ell}),$ (16a) $\displaystyle
N_{\gamma,{\rm A_{2}}}$ $\displaystyle=N_{\gamma,{\rm br}}-N_{\gamma,{\rm
A_{1}}},$ (16b) $\displaystyle N_{\gamma,{\rm B_{1}}}$
$\displaystyle=\frac{1}{4}\sum_{\ell=0}^{1}\sum_{\ell^{\prime}=0}^{1}\mathcal{D}(C_{2}^{\prime\ell^{\prime}}C_{2}^{\ell})\,N_{\gamma,{\rm
sb}}\,\mathcal{D}^{\dagger}(C_{2}^{\ell}C_{2}^{\prime\ell^{\prime}}),$ (16c)
$\displaystyle N_{\gamma,{\rm B_{2}}}$
$\displaystyle=\frac{1}{2}\sum_{\ell=0}^{1}\mathcal{D}(C_{2}^{\ell})\,N_{\gamma,{\rm
sb}}\,\mathcal{D}^{\dagger}(C_{2}^{\ell})-N_{\gamma,{\rm B_{1}}},$ (16d)
$\displaystyle N_{\gamma,{\rm E}}$ $\displaystyle=N_{\gamma,{\rm
sb}}-N_{\gamma,{\rm B_{1}}}-N_{\gamma,{\rm B_{2}}},$ (16e)

where $N_{\gamma,\Gamma}$ represents a matrix in combined orbital and spin
space transforming as the IR $\Gamma$ of the group $\rm D_{2d}$. The
$\gamma$-subscript indicates that the matrix couples to the form factor
$f_{\bf k}^{\gamma}$. Similarly, $N_{\gamma,{\rm br/sb}}$ is a matrix
containing the elements $N_{\gamma\sigma,{\rm br/sb}}^{\mu\nu}$. As discussed
in App. A, $\mathcal{D}(g)$ is the matrix representation of $g\in{\rm D_{2d}}$
in combined orbital and spin space.

By explicitly performing these averages for the 2D system under consideration,
we get the following matrices $N_{\gamma\sigma,\Gamma}$ in orbital space for a
given spin projection $\sigma=\\{\uparrow,\downarrow\\}\equiv\\{1,-1\\}$

$\displaystyle N_{\gamma\sigma,{\rm
A_{1}}}=\begin{pmatrix}\frac{1}{2}(n^{11}_{\gamma\sigma}+n^{22}_{\gamma\sigma})&\frac{1}{2}(m^{12}_{\gamma\sigma}-m^{21}_{\gamma\sigma})&0&0&0\\\\[5.69046pt]
-\frac{1}{2}(m^{12}_{\gamma\sigma}-m^{21}_{\gamma\sigma})&\frac{1}{2}(n^{11}_{\gamma\sigma}+n^{22}_{\gamma\sigma})&0&0&0\\\\[5.69046pt]
0&0&n^{33}_{\gamma\sigma}&m^{34}_{\gamma\sigma}&0\\\\[5.69046pt]
0&0&m^{43}_{\gamma\sigma}&n^{44}_{\gamma\sigma}&0\\\\[5.69046pt]
0&0&0&0&n^{55}_{\gamma\sigma}\end{pmatrix},$ (17a) $\displaystyle
N_{\gamma\sigma,{\rm
A_{2}}}=\begin{pmatrix}\frac{1}{2}(m^{11}_{\gamma\sigma}+m^{22}_{\gamma\sigma})&\frac{1}{2}(n^{12}_{\gamma\sigma}-n^{21}_{\gamma\sigma})&0&0&0\\\\[5.69046pt]
-\frac{1}{2}(n^{12}_{\gamma\sigma}-n^{21}_{\gamma\sigma})&\frac{1}{2}(m^{11}_{\gamma\sigma}+m^{22}_{\gamma\sigma})&0&0&0\\\\[5.69046pt]
0&0&m^{33}_{\gamma\sigma}&n^{34}_{\gamma\sigma}&0\\\\[5.69046pt]
0&0&n^{43}_{\gamma\sigma}&m^{44}_{\gamma\sigma}&0\\\\[5.69046pt]
0&0&0&0&m^{55}_{\gamma\sigma}\end{pmatrix},$ (17b) $\displaystyle
N_{\gamma\sigma,{\rm
B_{1}}}=\begin{pmatrix}\frac{1}{2}(n^{11}_{\gamma\sigma}-n^{22}_{\gamma\sigma})&\frac{1}{2}(m^{12}_{\gamma\sigma}+m^{21}_{\gamma\sigma})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(m^{12}_{\gamma\sigma}+m^{21}_{\gamma\sigma})&-\frac{1}{2}(n^{11}_{\gamma\sigma}-n^{22}_{\gamma\sigma})&0&0&0\\\\[5.69046pt]
0&0&0&0&n^{35}_{\gamma\sigma}\\\\[5.69046pt]
0&0&0&0&m^{45}_{\gamma\sigma}\\\\[5.69046pt]
0&0&n^{53}_{\gamma\sigma}&m^{54}_{\gamma\sigma}&0\end{pmatrix},$ (17c)
$\displaystyle N_{\gamma\sigma,{\rm
B_{2}}}=\begin{pmatrix}\frac{1}{2}(m^{11}_{\gamma\sigma}-m^{22}_{\gamma\sigma})&\frac{1}{2}(n^{12}_{\gamma\sigma}+n^{21}_{\gamma\sigma})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(n^{12}_{\gamma\sigma}+n^{21}_{\gamma\sigma})&-\frac{1}{2}(m^{11}_{\gamma\sigma}-m^{22}_{\gamma\sigma})&0&0&0\\\\[5.69046pt]
0&0&0&0&m^{35}_{\gamma\sigma}\\\\[5.69046pt]
0&0&0&0&n^{45}_{\gamma\sigma}\\\\[5.69046pt]
0&0&m^{53}_{\gamma\sigma}&n^{54}_{\gamma\sigma}&0\end{pmatrix},$ (17d)
$\displaystyle N_{\gamma\sigma,{\rm
E}_{x}}=\begin{pmatrix}0&0&\frac{1}{2}(m_{\gamma\sigma}^{13}+m_{\gamma\sigma}^{23})&\frac{1}{2}(n_{\gamma\sigma}^{14}+n_{\gamma\sigma}^{24})&\frac{1}{2}(m_{\gamma\sigma}^{15}-m_{\gamma\sigma}^{25})\\\\[5.69046pt]
0&0&\frac{1}{2}(-n_{\gamma\sigma}^{13}+n_{\gamma\sigma}^{23})&\frac{1}{2}(-m_{\gamma\sigma}^{14}+m_{\gamma\sigma}^{24})&\frac{1}{2}(n_{\gamma\sigma}^{15}+n_{\gamma\sigma}^{25})\\\
\frac{1}{2}(m_{\gamma\sigma}^{31}+m_{\gamma\sigma}^{32})&\frac{1}{2}(-n_{\gamma\sigma}^{31}+n_{\gamma\sigma}^{32})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(n_{\gamma\sigma}^{41}+n_{\gamma\sigma}^{42})&\frac{1}{2}(-m_{\gamma\sigma}^{41}+m_{\gamma\sigma}^{42})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(m_{\gamma\sigma}^{51}-m_{\gamma\sigma}^{52})&\frac{1}{2}(n_{\gamma\sigma}^{51}+n_{\gamma\sigma}^{52})&0&0&0\end{pmatrix},$
(17e) $\displaystyle N_{\gamma\sigma,{\rm
E}_{y}}=\begin{pmatrix}0&0&\frac{1}{2}(n_{\gamma\sigma}^{13}-n_{\gamma\sigma}^{23})&\frac{1}{2}(m_{\gamma\sigma}^{14}-m_{\gamma\sigma}^{24})&\frac{1}{2}(n_{\gamma\sigma}^{15}+n_{\gamma\sigma}^{25})\\\\[5.69046pt]
0&0&\frac{1}{2}(m_{\gamma\sigma}^{13}+m_{\gamma\sigma}^{23})&\frac{1}{2}(n_{\gamma\sigma}^{14}+n_{\gamma\sigma}^{24})&\frac{1}{2}(-m_{\gamma\sigma}^{15}+m_{\gamma\sigma}^{25})\\\\[5.69046pt]
\frac{1}{2}(n_{\gamma\sigma}^{31}-n_{\gamma\sigma}^{32})&\frac{1}{2}(m_{\gamma\sigma}^{31}+m_{\gamma\sigma}^{32})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(m_{\gamma\sigma}^{41}-m_{\gamma\sigma}^{42})&\frac{1}{2}(n_{\gamma\sigma}^{41}+n_{\gamma\sigma}^{42})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(n_{\gamma\sigma}^{51}+n_{\gamma\sigma}^{52})&\frac{1}{2}(-m_{\gamma\sigma}^{51}+m_{\gamma\sigma}^{52})&0&0&0\end{pmatrix},$
(17f) $\displaystyle N_{\gamma\sigma,{\rm
E}_{L_{x}}}=\begin{pmatrix}0&0&\frac{1}{2}(m_{\gamma\sigma}^{13}-m_{\gamma\sigma}^{23})&\frac{1}{2}(n_{\gamma\sigma}^{14}-n_{\gamma\sigma}^{24})&\frac{1}{2}(m_{\gamma\sigma}^{15}+m_{\gamma\sigma}^{25})\\\\[5.69046pt]
0&0&\frac{1}{2}(n_{\gamma\sigma}^{13}+n_{\gamma\sigma}^{23})&\frac{1}{2}(m_{\gamma\sigma}^{14}+m_{\gamma\sigma}^{24})&\frac{1}{2}(-n_{\gamma\sigma}^{15}+n_{\gamma\sigma}^{25})\\\\[5.69046pt]
\frac{1}{2}(m_{\gamma\sigma}^{31}-m_{\gamma\sigma}^{32})&\frac{1}{2}(n_{\gamma\sigma}^{31}+n_{\gamma\sigma}^{32})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(n_{\gamma\sigma}^{41}-n_{\gamma\sigma}^{42})&\frac{1}{2}(m_{\gamma\sigma}^{41}+m_{\gamma\sigma}^{42})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(m_{\gamma\sigma}^{51}+m_{\gamma\sigma}^{52})&\frac{1}{2}(-n_{\gamma\sigma}^{51}+n_{\gamma\sigma}^{52})&0&0&0\end{pmatrix},$
(17g) $\displaystyle N_{\gamma\sigma,{\rm
E}_{L_{y}}}=\begin{pmatrix}0&0&\frac{1}{2}(n_{\gamma\sigma}^{13}+n_{\gamma\sigma}^{23})&\frac{1}{2}(m_{\gamma\sigma}^{14}+m_{\gamma\sigma}^{24})&\frac{1}{2}(n_{\gamma\sigma}^{15}-n_{\gamma\sigma}^{25})\\\\[5.69046pt]
0&0&\frac{1}{2}(-m_{\gamma\sigma}^{13}+m_{\gamma\sigma}^{23})&\frac{1}{2}(-n_{\gamma\sigma}^{14}+n_{\gamma\sigma}^{24})&\frac{1}{2}(m_{\gamma\sigma}^{15}+m_{\gamma\sigma}^{25})\\\
\frac{1}{2}(n_{\gamma\sigma}^{31}+n_{\gamma\sigma}^{32})&\frac{1}{2}(-m_{\gamma\sigma}^{31}+m_{\gamma\sigma}^{32})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(m_{\gamma\sigma}^{41}+m_{\gamma\sigma}^{42})&\frac{1}{2}(-n_{\gamma\sigma}^{41}+n_{\gamma\sigma}^{42})&0&0&0\\\\[5.69046pt]
\frac{1}{2}(n_{\gamma\sigma}^{51}-n_{\gamma\sigma}^{52})&\frac{1}{2}(m_{\gamma\sigma}^{51}+m_{\gamma\sigma}^{52})&0&0&0\end{pmatrix},$
(17h)

where we introduced the following quantities for brevity

$\displaystyle n_{\gamma\sigma}^{\mu\nu}$
$\displaystyle=(N_{\gamma\sigma}^{\mu\nu}+N_{\gamma\bar{\sigma}}^{\mu\nu})/2,$
$\displaystyle m_{\gamma\sigma}^{\mu\nu}$
$\displaystyle=\sigma(N_{\gamma\sigma}^{\mu\nu}-N_{\gamma\bar{\sigma}}^{\mu\nu})/2,$
(18)

with $\sigma=-\bar{\sigma}$. Note furthermore that the matrices labeled by the
IR ${\rm E}_{x,y}$ transform as the basis functions $(x,y)$, while
$N_{\gamma\sigma,{\rm E}_{L_{x}}}$ and $N_{\gamma\sigma,{\rm E}_{L_{y}}}$
instead transform as the angular momentum pseudovector $(L_{x},L_{y})$. In
general we allow for magnetic terms, i.e. $m_{\gamma\sigma}^{\mu\nu}\neq 0$,
however, we find these fields to be orders of magnitude smaller than the
density terms $n^{\mu\nu}_{\gamma\sigma}$. Nonetheless, for completeness our
calculations and classification include all terms, so not to miss any subtle
details.

We can then, after having performed the various averages in Eq. (16), express
$N_{{\bf k}\sigma}^{\mu\nu}$ in terms of the IRs in the following way

$\displaystyle\begin{split}N_{{\bf k}\sigma}^{\mu\nu}&=\sum_{\gamma}f_{\bf
k}^{\gamma}\big{(}N_{\gamma\sigma,{\rm br}}^{\mu\nu}+N_{\gamma\sigma,{\rm
sb}}^{\mu\nu}\big{)}=\sum_{\gamma,\Gamma}f_{\bf
k}^{\gamma}N_{\gamma\sigma,\Gamma}^{\mu\nu}.\end{split}$ (19)

This informs us that a single term $f_{\bf
k}^{\gamma}N_{\gamma\sigma,\Gamma}^{\mu\nu}$ must transform as the product of
the two IRs $\gamma$ and $\Gamma$, e.g. for $\gamma={\rm d}\sim{\rm B_{1}}$
and $\Gamma={\rm A_{1}}$, the term altogether transforms as $\rm B_{1}\otimes
A_{1}=\rm B_{1}$. With this in mind, we can collect terms that transform
equivalently, and get

$\displaystyle N_{{\bf k}\sigma,\rm A_{1}}^{\mu\nu}=f_{\bf k}^{\,\rm s}N_{{\rm
s}\sigma,{\rm A_{1}}}^{\mu\nu}+f_{\bf k}^{\,\rm d}N_{{\rm d}\sigma,{\rm
B_{1}}}^{\mu\nu}$ (20) $\displaystyle+f_{\bf k}^{{\rm p}_{x}}\frac{N_{{\rm
p}_{x}\sigma,{\rm E}_{x}}+N_{{\rm p}_{x}\sigma,{\rm E}_{L_{x}}}}{2}+f_{\bf
k}^{{\rm p}_{y}}\frac{N_{{\rm p}_{x}\sigma,{\rm E}_{y}}-N_{{\rm
p}_{x}\sigma,{\rm E}_{L_{y}}}}{2}$ $\displaystyle+f_{\bf k}^{{\rm
p}_{x}}\frac{N_{{\rm p}_{y}\sigma,{\rm E}_{x}}-N_{{\rm p}_{y}\sigma,{\rm
E}_{L_{x}}}}{2}+f_{\bf k}^{{\rm p}_{y}}\frac{N_{{\rm p}_{y}\sigma,{\rm
E}_{y}}+N_{{\rm p}_{y}\sigma,{\rm E}_{L_{y}}}}{2},$ $\displaystyle
N^{\mu\nu}_{{\bf k}\sigma,\rm A_{2}}=f_{\bf k}^{\,\rm s}N_{{\rm s}\sigma,{\rm
A_{2}}}^{\mu\nu}+f_{\bf k}^{\,\rm d}N_{{\rm d}\sigma,{\rm B_{2}}}^{\mu\nu}$
(21) $\displaystyle+f_{\bf k}^{{\rm p}_{x}}\frac{N_{{\rm p}_{x}\sigma,{\rm
E}_{y}}+N_{{\rm p}_{x}\sigma,{\rm E}_{L_{y}}}}{2}-f_{\bf k}^{{\rm
p}_{y}}\frac{N_{{\rm p}_{x}\sigma,{\rm E}_{x}}-N_{{\rm p}_{x}\sigma,{\rm
E}_{L_{x}}}}{2}$ $\displaystyle-f_{\bf k}^{{\rm p}_{x}}\frac{N_{{\rm
p}_{y}\sigma,{\rm E}_{y}}-N_{{\rm p}_{y}\sigma,{\rm E}_{L_{y}}}}{2}+f_{\bf
k}^{{\rm p}_{y}}\frac{N_{{\rm p}_{y}\sigma,{\rm E}_{x}}+N_{{\rm
p}_{y}\sigma,{\rm E}_{L_{x}}}}{2},$ $\displaystyle N^{\mu\nu}_{{\bf
k}\sigma,\rm B_{1}}=f_{\bf k}^{\,\rm s}N_{{\rm s}\sigma,{\rm
B_{1}}}^{\mu\nu}+f_{\bf k}^{\,\rm d}N_{{\rm d}\sigma,{\rm A_{1}}}^{\mu\nu}$
(22) $\displaystyle+f_{\bf k}^{{\rm p}_{x}}\frac{N_{{\rm p}_{x},\sigma{\rm
E}_{x}}+N_{{\rm p}_{x}\sigma,{\rm E}_{L_{x}}}}{2}-f_{\bf k}^{{\rm
p}_{y}}\frac{N_{{\rm p}_{x}\sigma,{\rm E}_{y}}-N_{{\rm p}_{x}\sigma,{\rm
E}_{L_{y}}}}{2}$ $\displaystyle-f_{\bf k}^{{\rm p}_{x}}\frac{N_{{\rm
p}_{y}\sigma,{\rm E}_{x}}-N_{{\rm p}_{y}\sigma,{\rm E}_{L_{x}}}}{2}+f_{\bf
k}^{{\rm p}_{y}}\frac{N_{{\rm p}_{y},\sigma{\rm E}_{y}}+N_{{\rm
p}_{y},\sigma{\rm E}_{L_{y}}}}{2},$ $\displaystyle N^{\mu\nu}_{{\bf
k}\sigma,\rm B_{2}}=f_{\bf k}^{\,\rm s}N_{{\rm s}\sigma,{\rm
B_{2}}}^{\mu\nu}+f_{\bf k}^{\,\rm d}N_{{\rm d}\sigma,{\rm A_{2}}}^{\mu\nu}$
(23) $\displaystyle+f_{\bf k}^{{\rm p}_{x}}\frac{N_{{\rm p}_{x}\sigma,{\rm
E}_{y}}+N_{{\rm p}_{x}\sigma,{\rm E}_{L_{y}}}}{2}+f_{\bf k}^{{\rm
p}_{y}}\frac{N_{{\rm p}_{x}\sigma,{\rm E}_{x}}-N_{{\rm p}_{x}\sigma,{\rm
E}_{L_{x}}}}{2}$ $\displaystyle+f_{\bf k}^{{\rm p}_{x}}\frac{N_{{\rm
p}_{y}\sigma,{\rm E}_{y}}-N_{{\rm p}_{y}\sigma,{\rm E}_{L_{y}}}}{2}+f_{\bf
k}^{{\rm p}_{y}}\frac{N_{{\rm p}_{y}\sigma,{\rm E}_{x}}+N_{{\rm
p}_{y}\sigma,{\rm E}_{L_{x}}}}{2},$

for the IRs $\rm A_{1}$, $\rm A_{2}$, $\rm B_{1}$ and $\rm B_{2}$,
respectively, and finally for the 2D IR we arrive at

$\displaystyle N_{{\bf k}\sigma,\rm E}^{\mu\nu}=f_{\bf k}^{\,\rm s}(N_{{\rm
s}\sigma,{\rm E}_{x}}^{\mu\nu}+N_{{\rm s}\sigma,{\rm E}_{y}}^{\mu\nu})+f_{\bf
k}^{\,\rm s}(N_{{\rm s}\sigma,{\rm E}_{L_{x}}}^{\mu\nu}+N_{{\rm s}\sigma,{\rm
E}_{L_{y}}}^{\mu\nu})$ $\displaystyle+f_{\bf k}^{\,\rm d}(N_{{\rm
d}\sigma,{\rm E}_{x}}^{\mu\nu}+N_{{\rm d}\sigma,{\rm E}_{y}}^{\mu\nu})+f_{\bf
k}^{\,\rm d}(N_{{\rm d}\sigma,{\rm E}_{L_{x}}}^{\mu\nu}+N_{{\rm d}\sigma,{\rm
E}_{L_{y}}}^{\mu\nu})$ $\displaystyle+f_{\bf k}^{\,{\rm p}_{x}}(N_{{\rm
p}_{x}\sigma,{\rm A}_{1}}^{\mu\nu}+N_{{\rm p}_{x}\sigma,{\rm
B}_{1}}^{\mu\nu}+N_{{\rm p}_{x}\sigma,{\rm A}_{2}}^{\mu\nu}+N_{{\rm
p}_{x}\sigma,{\rm B}_{2}}^{\mu\nu})$ $\displaystyle+f_{\bf k}^{\,{\rm
p}_{y}}(N_{{\rm p}_{y}\sigma,{\rm A}_{1}}^{\mu\nu}+N_{{\rm p}_{y}\sigma,{\rm
B}_{1}}^{\mu\nu}+N_{{\rm p}_{y}\sigma,{\rm A}_{2}}^{\mu\nu}+N_{{\rm
p}_{y}\sigma,{\rm B}_{2}}^{\mu\nu}).$ (24)

From the above, we are able to pinpoint exactly what fields give rise to the
nematic order. For example, if a non-zero $\rm B_{1}$ term appears in our
self-consistent calculations, then we know that the nematic order transforms
as $\rm B_{1}$, and that the resulting crystalline point group must be $\rm
D_{2}\vartriangleleft\rm D_{2d}$.

From our self-consistent calculations, we find that all A2 and B2 terms are
identically zero, i.e. $N^{\mu\nu}_{{\bf k}\sigma,{\rm
A_{2}}}=N^{\mu\nu}_{{\bf k}\sigma,{\rm B_{2}}}=0$, which immediately implies
that $N^{\mu\nu}_{{\bf k}\sigma,{\rm br}}\equiv N^{\mu\nu}_{{\bf k}\sigma,{\rm
A_{1}}}$, thus the band renormalizing terms all transform as the IR A1. See
Fig. 6 for the values of $N_{{\bf k}\sigma,{\rm A_{1}}}^{\mu\nu}$, where we
straightforwardly conclude that

$\displaystyle N_{{\rm s}\sigma,{\rm br}}^{\mu\nu}$ $\displaystyle\equiv
N_{{\rm s}\sigma,{\rm A_{1}}}^{\mu\nu},$ $\displaystyle N_{{\rm d}\sigma,{\rm
br}}^{\mu\nu}$ $\displaystyle\equiv N_{{\rm d}\sigma,{\rm B_{1}}}^{\mu\nu}.$
(25)

Similarly we conclude

$\displaystyle\begin{split}N_{{\rm p}_{x}\sigma,{\rm
br}}^{\mu\nu}&\equiv\frac{N_{{\rm p}_{x}\sigma,{\rm E}_{x}}^{\mu\nu}+N_{{\rm
p}_{x}\sigma,{\rm E}_{L_{x}}}^{\mu\nu}+N_{{\rm p}_{y}\sigma,{\rm
E}_{x}}^{\mu\nu}-N_{{\rm p}_{x}\sigma,{\rm E}_{L_{x}}}^{\mu\nu}}{2},\\\
N_{{\rm p}_{y}\sigma,{\rm br}}^{\mu\nu}&\equiv\frac{N_{{\rm p}_{x}\sigma,{\rm
E}_{y}}^{\mu\nu}-N_{{\rm p}_{x}\sigma,{\rm E}_{L_{y}}}^{\mu\nu}+N_{{\rm
p}_{y}\sigma,{\rm E}_{y}}^{\mu\nu}+N_{{\rm p}_{x}\sigma,{\rm
E}_{L_{y}}}^{\mu\nu}}{2},\end{split}$ (26)

which can be inferred from Fig. 8(a) and (b), where we show the matrices with
the form factors $f_{\rm k}^{\,{\rm p}_{x}}$ and $f_{\rm k}^{\,{\rm p}_{y}}$
in Eq. (20).

We can furthermore extract from Fig. 8(a) and (b) that $N_{{\bf k}\sigma,{\rm
B_{1}}}^{\mu\nu}$ only involve terms coupling to $f_{\bf k}^{\rm s}$ and
$f_{\bf k}^{\rm d}$, since

$\displaystyle\begin{split}\frac{N_{{\rm p}_{x}\sigma,{\rm
E}_{x}}^{\mu\nu}+N_{{\rm p}_{x}\sigma,{\rm
E}_{L_{x}}}^{\mu\nu}}{2}-\frac{N_{{\rm p}_{y}\sigma,{\rm
E}_{x}}^{\mu\nu}-N_{{\rm p}_{x}\sigma,{\rm E}_{L_{x}}}^{\mu\nu}}{2}=0\\\
-\frac{N_{{\rm p}_{x}\sigma,{\rm E}_{y}}^{\mu\nu}-N_{{\rm p}_{x}\sigma,{\rm
E}_{L_{y}}}^{\mu\nu}}{2}+\frac{N_{{\rm p}_{y}\sigma,{\rm
E}_{y}}^{\mu\nu}+N_{{\rm p}_{x}\sigma,{\rm
E}_{L_{y}}}^{\mu\nu}}{2}=0,\end{split}$ (27)

i.e. all terms coupling to $f_{\bf k}^{\,{\rm p}_{x,y}}$ cancel in this IR
channel. From our self-consistent calculations, we also find that $N_{\rm
s\sigma,{\rm B_{1}}}^{\mu\nu}=0$, leaving us with only a single remaining term
in $N_{{\bf k}\sigma,{\rm B_{1}}}^{\mu\nu}$, namely $N^{\mu\nu}_{{\rm
d\sigma,{\rm A_{1}}}}$. See Fig. 8(c) and (d) for an illustation of the
matrices $N_{{\rm s}\sigma,{\rm B_{1}}}$ and $N_{{\rm d}\sigma,{\rm A_{1}}}$.
The non-zero matrix elements of the latter can serve as order parameters for
our nematic phase, since these break $S_{4}$ symmetry and reduce the
crystalline point group symmetry to the group $\rm D_{2}$. Specifically, we
define the $\rm B_{1}$ nematic order parameter as

$\displaystyle\begin{split}N_{\rm B_{1}}&=(N_{\rm d\sigma,A_{1}}^{11}+N_{\rm
d\sigma,A_{1}}^{22})/2\\\ &=(N_{\rm d\sigma}^{11}+N_{\rm
d\bar{\sigma}}^{11}+N_{\rm d\sigma}^{22}+N_{\rm d\bar{\sigma}}^{22})/4\\\
&=(N_{\rm d\sigma,{\rm sb}}^{11}+N_{\rm d\bar{\sigma},{\rm sb}}^{11}+N_{\rm
d\sigma,{\rm sb}}^{22}+N_{\rm d\bar{\sigma},{\rm sb}}^{22})/4.\end{split}$
(28)

Alternatively, one could also define the order parameter as $N^{44}_{{\rm
d}\sigma,{\rm A_{1}}}$, however, we choose to focus on $N_{\rm B_{1}}$ since
it acquires a slightly higher value.

For the terms transforming as the 2D IR $\rm E$, we find the non-zero matrices
displayed in Fig. 9, which ultimately result in a simplified expression for
$N_{{\bf k}\sigma,{\rm E}}$, namely

$\displaystyle N_{{\bf k}\sigma,{\rm E}}^{\mu\nu}$ $\displaystyle=$ (29)
$\displaystyle f_{\bf k}^{\,\rm d}$ $\displaystyle(N_{{\rm d}\sigma,{\rm
E}_{x}}^{\mu\nu}+N_{{\rm d}\sigma,{\rm E}_{L_{x}}}^{\mu\nu})+f_{\bf k}^{\,{\rm
p}_{y}}(N_{{\rm p}_{y}\sigma,{\rm A}_{2}}^{\mu\nu}+N_{{\rm p}_{y}\sigma,{\rm
B}_{2}}^{\mu\nu}).$

Similar to the B1 terms, also here we can adopt the non-zero elements of
$N^{\mu\nu}_{{\bf k}\sigma,{\rm E}}$ as an order parameter of the system.
Specifically we define the following two component order parameter

$\displaystyle\begin{split}N_{\rm E}&=\left(N^{14}_{{\rm d}\sigma,{\rm
E}_{x}}+N^{14}_{{\rm d}\sigma,{\rm E}_{L_{x}}},N^{24}_{{\rm d}\sigma,{\rm
E}_{y}}+N^{24}_{{\rm d}\sigma,{\rm E}_{L_{y}}}\right)\\\ &=\left([N_{{\rm
d}\sigma}^{14}+N_{{\rm d}\bar{\sigma}}^{14}]/2,[N_{{\rm d}\sigma}^{24}+N_{{\rm
d}\bar{\sigma}}^{24}]/2\right)\\\ &=\left([N_{{\rm d}\sigma,{\rm
sb}}^{14}+N_{{\rm d}\bar{\sigma},{\rm sb}}^{14}]/2,[N_{{\rm d}\sigma,{\rm
sb}}^{24}+N_{{\rm d}\bar{\sigma},{\rm sb}}^{24}]/2\right)\\\ &\equiv(N_{{\rm
E}_{x}},N_{{\rm E}_{y}}).\end{split}$ (30)

The order parameters $N_{{\rm E}_{x,y}}$ will enter on equal footing in a
Landau free energy expansion, since they are interrelated by symmetry, and the
system can thus display either $N_{{\rm E}_{x}}\neq 0$ or $N_{{\rm E}_{y}}\neq
0$. Through our self-consistent calculations we find in fact both solutions,
depending on the initialization of our computations, i.e. where in the energy
landscape the calculations start. Throughout the paper we have focused on
solutions with $N_{{\rm E}_{x}}\neq 0$.

Conclusively, we see that the presence of two non-zero order parameters,
belonging to distinct IRs, hints at the following two phase transition
scenarios: $i)$ By lowering of the temperature, $N_{\rm E}$ becomes non-zero
which imposes the symmetry point group transition $\rm D_{2d}\mapsto{\rm
C}_{2}$. For this specific point group $N_{\rm B_{1}}$ transform trivially,
and is thus allowed to enter simultaneously with $N_{\rm E}$, i.e. $N_{\rm E}$
and $N_{\rm B}$ emerges at the same transition temperature. The other scenario
is $ii)$ The system undergoes two consecutive phase transitions, with $N_{\rm
B_{1}}$ entering prior to $N_{\rm E}$ when lowering the temperature, leading
to two transition temperatures and the following point group reduction scheme
${\rm D}_{\rm 2d}\mapsto{\rm D}_{\rm 2}\mapsto{\rm C}_{\rm 2}$. For the
spicific system under consideration, we find scenario $ii)$ to be true, see
Fig. 3 where we display the order parameters $N_{\rm B}$ and $N_{\rm E}$ at
various temperatures. Nonetheless, we stress that case $i)$ potentially also
could arise in experiments.

Figure 8: (a) and (b) Matrix structure of fields coupled to the p-wave form
factors $f_{\bf k}^{{\rm p}_{x}}$ and $f_{\bf k}^{{\rm p}_{y}}$, respectively.
These fields enter both in $N_{{\bf k}\sigma,{\rm A_{1}}}$ and $N_{{\bf
k}\sigma,{\rm B_{1}}}$, but only lead to a non-zero contribution in the
former, see Eqs. (20) and (22). (c) and (d) Remaining contributions entering
in $N_{{\bf k}\sigma,{\rm B_{1}}}$. The non-zero elements in (d), can be
utilized as the nematic order parameter, since they break $S_{4}$ symmetry.
See Eq. (28) for the resulting B1 neamtic order parameter. All values in the
figure have been rounded to the second decimal place. Figure 9: Non-zero
matrices entering in the 2D IR fields in Eq. (B). The matrices in (a) and (b)
[(c) and (d)] couple to the form factor $f^{\,\rm d}_{\bf k}$ ($f^{\,{\rm
p}_{y}}_{\bf k}$). As for the B1 field, also here the non-zero elements allow
us to define an order parameter, see Eq. (30), which can be used as an
indicator for when the symmetries of the system are described by the
monoclinic point group $\rm C_{2}$. All values in the figure have been rounded
to the second decimal place.

##  Appendix C Susceptibility and spin-fluctuation pairing from a partially
incoherent electronic structure

Adopting the properties of a correlated electron gas which is characterized by
a reduced quasiparticle weight $Z({\bf k}_{\mathrm{F}})$ at the Fermi level,
it seems a good approximation at low energies to parametrize the Green
function as

$\displaystyle G_{\mu\nu}({\bf k},\omega_{n})$
$\displaystyle=\sqrt{Z_{\mu}Z_{\nu}}\;\sum_{m}a_{m}^{\mu}({\bf
k})a_{m}^{\nu*}({\bf k})G^{m}({\bf k},\omega_{n}),$ (31)

where $Z_{\mu}$ are quasiparticle weights in orbital $\mu$, $E_{{\bf k}m}$ are
the eigenenergies of band $m$ of the mean field Hamiltonian in Eq. (5),
$a_{m}^{\mu}({\bf k})$ the corresponding orbital component of the eigenstate
and $G^{m}({\bf k},\omega_{n})=[i\omega_{n}-E_{{\bf k}m}]^{-1}$ the Green
function of band $m$. Adopting the usual local interactions from the Hubbard-
Hund Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle{U}\sum_{i,\mu}n_{i\mu\uparrow}n_{i\mu\downarrow}+{U}^{\prime}\sum_{i,\nu<\mu}n_{i\mu}n_{i\nu}$
(32) $\displaystyle+$
$\displaystyle{J}\sum_{i,\nu<\mu}\sum_{\sigma,\sigma^{\prime}}c_{i\mu\sigma}^{\dagger}c_{i\nu\sigma^{\prime}}^{\dagger}c_{i\mu\sigma^{\prime}}c_{i\nu\sigma}$
$\displaystyle+$
$\displaystyle{J}^{\prime}\sum_{i,\nu\neq\mu}c_{i\mu\uparrow}^{\dagger}c_{i\mu\downarrow}^{\dagger}c_{i\nu\downarrow}c_{i\nu\uparrow},$

where the parameters ${U}$, ${U}^{\prime}$, ${J}$, ${J}^{\prime}$ are given in
the notation of Kuroki et al. [71], we stay in the regime where
$U^{\prime}=U-2J$, $J=J^{\prime}$ and use the overall interaction magnitude
$U$ as free parameter to tune close to the magnetic instability and fix
$J/U=1/6$ as it has been used earlier for FeSe[11]. Within this current
Ansatz, the orbital susceptibility in the normal state is given by

$\displaystyle\tilde{\chi}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}^{0}(q)$
$\displaystyle=-\\!\\!\\!\sum_{k,m,m^{\prime}}\\!\\!M_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}^{mm^{\prime}}({\bf
k},{\bf q})G^{m}(k+q)G^{m^{\prime}}(k),$ (33)

where we have adopted the shorthand notation $k\equiv({\bf k},\omega_{n})$ and
defined the abbreviation

$\displaystyle M_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}^{mn}({\bf k},{\bf q})$
$\displaystyle=$
$\displaystyle\sqrt{Z_{\mu_{1}}Z_{\mu_{2}}Z_{\mu_{3}}Z_{\mu_{4}}}$
$\displaystyle\\!\\!\times a_{n}^{\mu_{4}}({\bf k})a_{n}^{\mu_{2},*}({\bf
k})a_{m}^{\mu_{1}}({\bf k}\\!+\\!{\bf q})a_{m}^{\mu_{3},*}({\bf k}\\!+\\!{\bf
q}).$

To evaluate the susceptibility, we perform the frequency summation
analytically and numerically sum over the full Brillouin zone using $400\times
400$ $k$ points. Note that the susceptibility is just related by

$\tilde{\chi}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}^{0}({\bf
q})=\sqrt{Z_{\mu_{1}}Z_{\mu_{2}}Z_{\mu_{3}}Z_{\mu_{4}}}\;\chi_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}^{0}({\bf
q}),$ (35)

to the one from a fully coherent electronic structure, i.e. the quasiparticle
weights enter as prefactors and renormalize each component of the
susceptibility tensor. Finally, we treat the interactions in a random phase
approximation (RPA) to calculate the spin (1) and charge (0) susceptibilities

$\displaystyle\tilde{\chi}_{1\,\mu_{1}\mu_{2}\mu_{3}\mu_{4}}^{\rm RPA}({\bf
q})$
$\displaystyle=\left\\{\tilde{\chi}^{0}(q)\left[1-\bar{U}^{s}\tilde{\chi}^{0}(q)\right]^{-1}\right\\}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}},$
(36a) $\displaystyle\tilde{\chi}_{0\,\mu_{1}\mu_{2}\mu_{3}\mu_{4}}^{\rm
RPA}({\bf q})$
$\displaystyle=\left\\{\tilde{\chi}^{0}(q)\left[1+\bar{U}^{c}\tilde{\chi}^{0}(q)\right]^{-1}\right\\}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}.$
(36b)

The total spin susceptibility as measured experimentally is given by the sum

$\tilde{\chi}({\bf
q},\omega)=\frac{1}{2}\sum_{\mu\nu}\tilde{\chi}_{1\;\mu\mu\nu\nu}^{\rm
RPA}({\bf q},\omega)\,.$ (37)

Note that the interaction matrices $\bar{U}^{s}$ and $\bar{U}^{c}$ contain
linear combinations of the parameters $U,U^{\prime},J,J^{\prime}$, for details
see for example Ref. 72.

To calculate the superconducting instability in the spin-singlet channel (the
dominant one for the present models), we use the vertex for pair scattering
between bands $n$ and $m$,

$\displaystyle\tilde{\Gamma}_{nm}(\mathbf{k},\mathbf{k}^{\prime})=\mathrm{Re}\sum_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}a_{n}^{\mu_{1},*}(\mathbf{k})a_{n}^{\mu_{4},*}(-\mathbf{k})$
$\displaystyle\hskip
17.07182pt\times{\tilde{\Gamma}}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}({\bf k},{\bf
k}^{\prime})a_{m}^{\mu_{2}}(\mathbf{k}^{\prime})a_{m}^{\mu_{3}}(-\mathbf{k}^{\prime})\,\,,$
(38)

where ${\bf k}$ and ${\bf k}^{\prime}$ are momenta restricted to the pockets
${\bf k}\in C_{n}$ and ${\bf k}^{\prime}\in C_{m}$, and is defined in terms of
the the orbital space vertex function

$\displaystyle{\tilde{\Gamma}}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}({\bf k},{\bf
k}^{\prime})=\left[\frac{3}{2}\bar{U}^{s}\tilde{\chi}_{1}^{\rm RPA}({\bf
k}-{\bf k}^{\prime})\bar{U}^{s}\right.\,\,$ (39)
$\displaystyle\,\left.+\frac{1}{2}\bar{U}^{s}-\frac{1}{2}\bar{U}^{c}\tilde{\chi}_{0}^{\rm
RPA}({\bf k}-{\bf
k}^{\prime})\bar{U}^{c}+\frac{1}{2}\bar{U}^{c}\right]_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$

Then, the linearized gap equation

$-\frac{1}{V_{G}}\sum_{m}\int_{\text{FS}_{m}}dS^{\prime}\;\tilde{\Gamma}_{nm}({\bf
k},{\bf k}^{\prime})\frac{g_{i}({\bf k}^{\prime})}{|v_{\text{F}m}({\bf
k}^{\prime})|}=\lambda_{i}g_{i}({\bf k})$ (40)

describes the superconducting gap $\Delta({\bf k})\propto g({\bf k})$ for the
largest eigenvalue $\lambda$ at least at $T_{c}$. The integration is over the
Fermi surface $\text{FS}_{m}$, the Fermi velocity $v_{\text{F}m}({\bf
k}^{\prime})$ enters as weights in the denominator and $V_{G}$ is the volume
of the Brillouin zone. For the uncorrelated case (Figs. 4, 5 (a)) we have used
$Z_{\alpha}=1$, while for the weakly correlated case we start with
$\sqrt{Z_{\alpha}}=[0.72,0.89,0.77,0.77,0.85]$ [15], subsequently reduce
$\sqrt{Z_{xy}}$ by steps of $0.1$, Figs. 4, 5(b), and finally split
$\sqrt{Z_{xz,yz}}=0.77\mp 0.02s$, $s=1,2,3,4$, Figs. 4, 5(c).

## References

  * Böhmer and Kreisel [2018] A. E. Böhmer and A. Kreisel, Nematicity, magnetism and superconductivity in FeSe, Journal of Physics: Condensed Matter 30, 023001 (2018).
  * Coldea and Watson [2018] A. I. Coldea and M. D. Watson, The key ingredients of the electronic structure of FeSe, Annual Review of Condensed Matter Physics 9, 125 (2018).
  * Kreisel  _et al._ [2020] A. Kreisel, P. J. Hirschfeld, and B. M. Andersen, On the remarkable superconductivity of FeSe and its close cousins, Symmetry 12, 10.3390/sym12091402 (2020).
  * Fernandes  _et al._ [2014] R. M. Fernandes, A. V. Chubukov, and J. Schmalian, What drives nematic order in iron-based superconductors?, Nature Physics 10, 97 (2014).
  * Chen  _et al._ [2019] T. Chen, Y. Chen, A. Kreisel, X. Lu, A. Schneidewind, Y. Qiu, J. T. Park, T. G. Perring, J. R. Stewart, H. Cao, R. Zhang, Y. Li, Y. Rong, Y. Wei, B. M. Andersen, P. J. Hirschfeld, C. Broholm, and P. Dai, Anisotropic spin fluctuations in detwinned FeSe, Nature Materials 18, 709 (2019).
  * Sprau  _et al._ [2017] P. O. Sprau, A. Kostin, A. Kreisel, A. E. Böhmer, V. Taufour, P. C. Canfield, S. Mukherjee, P. J. Hirschfeld, B. M. Andersen, and J. C. S. Davis, Discovery of orbital-selective Cooper pairing in FeSe, Science 357, 75 (2017).
  * Böhmer  _et al._ [2013] 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 (2013).
  * Tanatar  _et al._ [2016] M. A. Tanatar, A. E. Böhmer, E. I. Timmons, M. Schütt, G. Drachuck, V. Taufour, K. Kothapalli, A. Kreyssig, S. L. Bud’ko, P. C. Canfield, R. M. Fernandes, and R. Prozorov, Origin of the resistivity anisotropy in the nematic phase of FeSe, Phys. Rev. Lett. 117, 127001 (2016).
  * Kostin  _et al._ [2018] A. Kostin, P. O. Sprau, A. Kreisel, Y. X. Chong, A. E. Böhmer, P. C. Canfield, P. J. Hirschfeld, B. M. Andersen, and J. C. S. Davis, Imaging orbital-selective quasiparticles in the Hund’s metal state of FeSe, Nature Materials 17, 869 (2018).
  * Zhou  _et al._ [2020] R. Zhou, D. D. Scherer, H. Mayaffre, P. Toulemonde, M. Ma, Y. Li, B. M. Andersen, and M. H. Julien, Singular magnetic anisotropy in the nematic phase of FeSe, arXiv e-prints , arXiv:2010.01020 (2020), arXiv:2010.01020 [cond-mat.supr-con] .
  * Kreisel  _et al._ [2017] A. Kreisel, B. M. Andersen, P. O. Sprau, A. Kostin, J. C. S. Davis, and P. J. Hirschfeld, Orbital selective pairing and gap structures of iron-based superconductors, Phys. Rev. B 95, 174504 (2017).
  * Benfatto  _et al._ [2018] L. Benfatto, B. Valenzuela, and L. Fanfarillo, Nematic pairing from orbital selective spin fluctuations in FeSe, npj Quantum Materials 3, 56 (2018).
  * Hu  _et al._ [2018] H. Hu, R. Yu, E. M. Nica, J.-X. Zhu, and Q. Si, Orbital-selective superconductivity in the nematic phase of FeSe, Phys. Rev. B 98, 220503 (2018).
  * de’ Medici  _et al._ [2014] L. de’ Medici, G. Giovannetti, and M. Capone, Selective Mott physics as a key to iron superconductors, Phys. Rev. Lett. 112, 177001 (2014).
  * Björnson  _et al._ [2020] K. Björnson, A. Kreisel, A. T. Rømer, and B. M. Andersen, Orbital-dependent self-energy effects and consequences for the superconducting gap structure in multi-orbital correlated electron systems, arXiv e-prints , arXiv:2011.04243 (2020), arXiv:2011.04243 [cond-mat.supr-con] .
  * Cercellier  _et al._ [2019] H. Cercellier, P. Rodière, P. Toulemonde, C. Marcenat, and T. Klein, Influence of the quasiparticle spectral weight in FeSe on spectroscopic, magnetic, and thermodynamic properties, Phys. Rev. B 100, 104516 (2019).
  * Biswas  _et al._ [2018] P. K. Biswas, A. Kreisel, Q. Wang, D. T. Adroja, A. D. Hillier, J. Zhao, R. Khasanov, J.-C. Orain, A. Amato, and E. Morenzoni, Evidence of nodal gap structure in the basal plane of the FeSe superconductor, Phys. Rev. B 98, 180501 (2018).
  * Georges  _et al._ [2013] A. Georges, L. d. Medici, and J. Mravlje, Strong correlations from Hund’s coupling, Annual Review of Condensed Matter Physics 4, 137 (2013).
  * van Roekeghem  _et al._ [2016] A. van Roekeghem, P. Richard, H. Ding, and S. Biermann, Spectral properties of transition metal pnictides and chalcogenides: Angle-resolved photoemission spectroscopy and dynamical mean-field theory, C. R. Phys. 17, 140 (2016), iron-based superconductors / Supraconducteurs à base de fer.
  * Yi  _et al._ [2017] M. Yi, Y. Zhang, Z.-X. Shen, and D. Lu, Role of the orbital degree of freedom in iron-based superconductors, npj Quantum Materials 2, 57 (2017).
  * Guterding  _et al._ [2017] D. Guterding, S. Backes, M. Tomić, H. O. Jeschke, and R. Valentí, Ab initio perspective on structural and electronic properties of iron-based superconductors, physica status solidi (b) 254, 1600164 (2017), 1600164.
  * [22] L. de’ Medici, The physics of correlated insulators, metals, and superconductors (modeling and simulation vol. 7) (Forschungszentrum Juelich, Juelich, 2017) Chap. Hund’s Metals Explained, pp. 377–398.
  * Kreisel  _et al._ [2018] A. Kreisel, B. M. Andersen, and P. J. Hirschfeld, Itinerant approach to magnetic neutron scattering of FeSe: Effect of orbital selectivity, Phys. Rev. B 98, 214518 (2018).
  * Watson  _et al._ [2017] M. D. Watson, A. A. Haghighirad, L. C. Rhodes, M. Hoesch, and T. K. Kim, Electronic anisotropies revealed by detwinned angle-resolved photo-emission spectroscopy measurements of FeSe, New Journal of Physics 19, 103021 (2017).
  * Pfau  _et al._ [2019] H. Pfau, S. D. Chen, M. Yi, M. Hashimoto, C. R. Rotundu, J. C. Palmstrom, T. Chen, P.-C. Dai, J. Straquadine, A. Hristov, R. J. Birgeneau, I. R. Fisher, D. Lu, and Z.-X. Shen, Momentum dependence of the nematic order parameter in iron-based superconductors, Phys. Rev. Lett. 123, 066402 (2019).
  * Rhodes  _et al._ [2020] L. C. Rhodes, M. D. Watson, A. A. Haghighirad, D. V. Evtushinsky, and T. K. Kim, Revealing the single electron pocket of FeSe in a single orthorhombic domain, Phys. Rev. B 101, 235128 (2020).
  * Liu  _et al._ [2018] D. Liu, C. Li, J. Huang, B. Lei, L. Wang, X. Wu, B. Shen, Q. Gao, Y. Zhang, X. Liu, Y. Hu, Y. Xu, A. Liang, J. Liu, P. Ai, L. Zhao, S. He, L. Yu, G. Liu, Y. Mao, X. Dong, X. Jia, F. Zhang, S. Zhang, F. Yang, Z. Wang, Q. Peng, Y. Shi, J. Hu, T. Xiang, X. Chen, Z. Xu, C. Chen, and X. J. Zhou, Orbital origin of extremely anisotropic superconducting gap in nematic phase of FeSe superconductor, Phys. Rev. X 8, 031033 (2018).
  * Yi  _et al._ [2019] 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).
  * Huh  _et al._ [2020] 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, Communications Physics 3, 52 (2020).
  * Mukherjee  _et al._ [2015] S. Mukherjee, A. Kreisel, P. J. Hirschfeld, and B. M. Andersen, Model of electronic structure and superconductivity in orbitally ordered FeSe, Phys. Rev. Lett. 115, 026402 (2015).
  * Liang Yi [2015] H. J.-P. Liang Yi, Wu Xian-Xin, Electronic structure properties in the nematic phases of FeSe, Chinese Physics Letters 32, 117402 (2015).
  * Suzuki  _et al._ [2015] 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).
  * Watson  _et al._ [2016] M. D. Watson, T. K. Kim, L. C. Rhodes, M. Eschrig, M. Hoesch, A. A. Haghighirad, and A. I. Coldea, Evidence for unidirectional nematic bond ordering in FeSe, Phys. Rev. B 94, 201107 (2016).
  * Fanfarillo  _et al._ [2016] L. Fanfarillo, J. Mansart, P. Toulemonde, H. Cercellier, P. Le Fèvre, F. Bertran, B. Valenzuela, L. Benfatto, and V. Brouet, Orbital-dependent Fermi surface shrinking as a fingerprint of nematicity in FeSe, Phys. Rev. B 94, 155138 (2016).
  * Fedorov  _et al._ [2016] A. Fedorov, A. Yaresko, T. K. Kim, Y. Kushnirenko, E. Haubold, T. Wolf, M. Hoesch, A. Grüneis, B. Büchner, and S. V. Borisenko, Effect of nematic ordering on electronic structure of FeSe, Scientific Reports 6, 36834 (2016).
  * Onari  _et al._ [2016] 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).
  * Zhang  _et al._ [2015] P. Zhang, T. Qian, P. Richard, X. P. Wang, H. Miao, B. Q. Lv, B. B. Fu, T. Wolf, C. Meingast, X. X. Wu, Z. Q. Wang, J. P. Hu, and H. Ding, Observation of two distinct ${d}_{xz}$/${d}_{yz}$ band splittings in FeSe, Phys. Rev. B 91, 214503 (2015).
  * Jiang  _et al._ [2016] K. Jiang, J. Hu, H. Ding, and Z. Wang, Interatomic Coulomb interaction and electron nematic bond order in FeSe, Phys. Rev. B 93, 115138 (2016).
  * Xing  _et al._ [2017] R.-Q. Xing, L. Classen, M. Khodas, and A. V. Chubukov, Competing instabilities, orbital ordering, and splitting of band degeneracies from a parquet renormalization group analysis of a four-pocket model for iron-based superconductors: Application to FeSe, Phys. Rev. B 95, 085108 (2017).
  * Eugenio and Vafek [2018] P. M. Eugenio and O. Vafek, Classification of symmetry derived pairing at the $M$ point in FeSe, Phys. Rev. B 98, 014503 (2018).
  * Christensen  _et al._ [2019] M. H. Christensen, J. Kang, and R. M. Fernandes, Intertwined spin-orbital coupled orders in the iron-based superconductors, Phys. Rev. B 100, 014512 (2019).
  * Christensen  _et al._ [2020] M. H. Christensen, R. M. Fernandes, and A. V. Chubukov, Orbital transmutation and the electronic spectrum of FeSe in the nematic phase, Phys. Rev. Research 2, 013015 (2020).
  * Rhodes  _et al._ [2020] L. C. Rhodes, J. Böker, M. A. Müller, M. Eschrig, and I. M. Eremin, Non-local $d_{xy}$ nematicity and the missing electron pocket in FeSe, arXiv e-prints , arXiv:2009.00507 (2020), arXiv:2009.00507 [cond-mat.supr-con] .
  * Long  _et al._ [2020] X. Long, S. Zhang, F. Wang, and Z. Liu, A first-principle perspective on electronic nematicity in FeSe, npj Quantum Materials 5, 50 (2020).
  * Note [1] Applying a similar approach to LiFeAs, for example, produces a dramatically distorted low-energy band structure with missing $d_{xy}$ band at the Fermi level. R. Valenti, private communication.
  * Scherer  _et al._ [2017] D. D. Scherer, A. C. Jacko, C. Friedrich, E. Şaşıoğlu, S. Blügel, R. Valentí, and B. M. Andersen, Interplay of nematic and magnetic orders in FeSe under pressure, Phys. Rev. B 95, 094504 (2017).
  * Zantout  _et al._ [2019] K. Zantout, S. Backes, and R. Valentí, Effect of nonlocal correlations on the electronic structure of LiFeAs, Phys. Rev. Lett. 123, 256401 (2019).
  * Bhattacharyya  _et al._ [2020a] S. Bhattacharyya, P. J. Hirschfeld, T. A. Maier, and D. J. Scalapino, Effects of momentum-dependent quasiparticle renormalization on the gap structure of iron-based superconductors, Phys. Rev. B 101, 174509 (2020a).
  * Eschrig and Koepernik [2009] H. Eschrig and K. Koepernik, Tight-binding models for the iron-based superconductors, Phys. Rev. B 80, 104503 (2009).
  * Chubukov  _et al._ [2016] 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).
  * Bhattacharyya  _et al._ [2020b] S. Bhattacharyya, K. Björnson, K. Zantout, D. Steffensen, L. Fanfarillo, A. Kreisel, R. Valentí, B. M. Andersen, and P. J. Hirschfeld, Nonlocal correlations in iron pnictides and chalcogenides, Phys. Rev. B 102, 035109 (2020b).
  * Gastiasoro and Andersen [2015] M. N. Gastiasoro and B. M. Andersen, Competing magnetic double-$q$ phases and superconductivity-induced reentrance of ${C}_{2}$ magnetic stripe order in iron pnictides, Phys. Rev. B 92, 140506 (2015).
  * Cvetkovic and Vafek [2013] V. Cvetkovic and O. Vafek, Space group symmetry, spin-orbit coupling, and the low-energy effective Hamiltonian for iron-based superconductors, Phys. Rev. B 88, 134510 (2013).
  * Rahn  _et al._ [2015] M. C. Rahn, R. A. Ewings, S. J. Sedlmaier, S. J. Clarke, and A. T. Boothroyd, Strong $(\pi,0)$ spin fluctuations in $\beta-$FeSe observed by neutron spectroscopy, Phys. Rev. B 91, 180501 (2015).
  * Wang  _et al._ [2016] Q. Wang, Y. Shen, B. Pan, Y. Hao, M. Ma, F. Zhou, P. Steffens, K. Schmalzl, T. R. Forrest, M. Abdel-Hafiez, X. Chen, D. A. Chareev, A. N. Vasiliev, P. Bourges, Y. Sidis, H. Cao, and J. Zhao, Strong interplay between stripe spin fluctuations, nematicity and superconductivity in FeSe, Nature Materials 15, 159 (2016).
  * Fanfarillo  _et al._ [2018] L. Fanfarillo, L. Benfatto, and B. Valenzuela, Orbital mismatch boosting nematic instability in iron-based superconductors, Phys. Rev. B 97, 121109 (2018).
  * Kreisel  _et al._ [2015] A. Kreisel, S. Mukherjee, P. J. Hirschfeld, and B. M. Andersen, Spin excitations in a model of FeSe with orbital ordering, Phys. Rev. B 92, 224515 (2015).
  * Christensen  _et al._ [2016] M. H. Christensen, J. Kang, B. M. Andersen, and R. M. Fernandes, Spin-driven nematic instability of the multiorbital Hubbard model: Application to iron-based superconductors, Phys. Rev. B 93, 085136 (2016).
  * Miyake  _et al._ [2010] T. Miyake, K. Nakamura, R. Arita, and M. Imada, Comparison of ab initio low-energy models for LaFePO, LaFeAsO, BaFe2As2, LiFeAs, FeSe, and FeTe: Electron correlation and covalency, Journal of the Physical Society of Japan 79, 044705 (2010).
  * Yin  _et al._ [2011] Z. P. Yin, K. Haule, and G. Kotliar, Kinetic frustration and the nature of the magnetic and paramagnetic states in iron pnictides and iron chalcogenides, Nat. Mater. 10, 932 (2011).
  * Jiao  _et al._ [2017] L. Jiao, C.-L. Huang, S. Rößler, C. Koz, U. K. Rößler, U. Schwarz, and S. Wirth, Superconducting gap structure of FeSe, Scientific Reports 7, 44024 EP (2017), article.
  * Muratov  _et al._ [2018] A. Muratov, A. Sadakov, S. Gavrilkin, A. Prishchepa, G. Epifanova, D. Chareev, and V. Pudalov, Specific heat of FeSe: Two gaps with different anisotropy in superconducting state, Physica B: Condensed Matter 536, 785 (2018).
  * Böhmer  _et al._ [2015] 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).
  * Koch  _et al._ [2019] R. J. Koch, T. Konstantinova, M. Abeykoon, A. Wang, C. Petrovic, Y. Zhu, E. S. Bozin, and S. J. L. Billinge, Room temperature local nematicity in FeSe superconductor, Phys. Rev. B 100, 020501 (2019).
  * Frandsen  _et al._ [2019] B. A. Frandsen, Q. Wang, S. Wu, J. Zhao, and R. J. Birgeneau, Quantitative characterization of short-range orthorhombic fluctuations in FeSe through pair distribution function analysis, Phys. Rev. B 100, 020504 (2019).
  * Wang  _et al._ [2019] Z. Wang, X. Zhao, R. Koch, S. J. L. Billinge, and A. Zunger, Understanding Electronic Peculiarities in Tetragonal FeSe as Local Structural Symmetry Breaking, arXiv e-prints , arXiv:1911.02670 (2019), arXiv:1911.02670 [cond-mat.mtrl-sci] .
  * Eckberg  _et al._ [2020] C. Eckberg, D. J. Campbell, T. Metz, J. Collini, H. Hodovanets, T. Drye, P. Zavalij, M. H. Christensen, R. M. Fernandes, S. Lee, P. Abbamonte, J. W. Lynn, and J. Paglione, Sixfold enhancement of superconductivity in a tunable electronic nematic system, Nature Physics 16, 346 (2020).
  * Gastiasoro  _et al._ [2014] M. N. Gastiasoro, I. Paul, Y. Wang, P. J. Hirschfeld, and B. M. Andersen, Emergent defect states as a source of resistivity anisotropy in the nematic phase of iron pnictides, Phys. Rev. Lett. 113, 127001 (2014).
  * Wang  _et al._ [2015] Y. Wang, M. N. Gastiasoro, B. M. Andersen, M. Tomić, H. O. Jeschke, R. Valentí, I. Paul, and P. J. Hirschfeld, Effects of Lifshitz transition on charge transport in magnetic phases of Fe-based superconductors, Phys. Rev. Lett. 114, 097003 (2015).
  * Steffensen  _et al._ [2019] D. Steffensen, P. Kotetes, I. Paul, and B. M. Andersen, Disorder-induced electronic nematicity, Phys. Rev. B 100, 064521 (2019).
  * Kuroki  _et al._ [2008] K. Kuroki, S. Onari, R. Arita, H. Usui, Y. Tanaka, H. Kontani, and H. Aoki, Unconventional pairing originating from the disconnected Fermi surfaces of superconducting LaFeAsO1-xFx, Phys. Rev. Lett. 101, 087004 (2008).
  * Kemper  _et al._ [2010] A. F. Kemper, T. A. Maier, S. Graser, H.-P. Cheng, P. J. Hirschfeld, and D. J. Scalapino, Sensitivity of the superconducting state and magnetic susceptibility to key aspects of electronic structure in ferropnictides, New J. Phys. 12, 073030 (2010).