# Classification of electronic nematicity in three-dimensional crystals and
quasicrystals

Matthias Hecker  School of Physics and Astronomy, University of Minnesota,
Minneapolis 55455 MN, USA  Anant Rastogi  School of Physics and Astronomy,
University of Minnesota, Minneapolis 55455 MN, USA  Daniel F. Agterberg
Department of Physics, University of Wisconsin-Milwaukee, Milwaukee 53201 WI,
USA  Rafael M. Fernandes  School of Physics and Astronomy, University of
Minnesota, Minneapolis 55455 MN, USA

(August 28, 2024)

###### Abstract

Electronic nematic order has been reported in a rich landscape of materials,
encompassing not only a range of intertwined correlated and topological
phenomena, but also different underlying lattice symmetries. Motivated by
these findings, we investigate the behavior of electronic nematicity as the
spherical symmetry of three-dimensional (3D) space is systematically lowered
by the lattice environment. We consider all $32$ crystallographic point groups
as well as $4$ major classes of quasicrystalline point groups, given the
recent observations of electronic phases of interest in quasicrystalline
materials and artificial twisted quasicrystals. Valuable insights are gained
by establishing a mapping between the five-component charge-quadrupolar
nematic order parameter of the electronic fluid and the 3D tensorial order
parameter of nematic liquid crystals. We find that a uniaxial nematic state is
only generically realized in polyhedral point groups (icosahedral and cubic),
with the nematic director pointing along different sets of rotational symmetry
axes. Interestingly, icosahedral point groups are the only ones in which the
five nematic order parameter components transform as the same irreducible
representation, making them the closest analog of 3D isotropic nematics. In
axial point groups, one of the nematic components is always condensed, whereas
the other four components decompose into an in-plane and an out-of-plane
nematic doublet, resulting in biaxial nematic ground states. Because these two
nematic doublets behave as $Z_{q}$-clock order parameters, this allows us to
identify the types of crystals and quasicrystals that can host interesting
electronic nematic phenomena enabled by the critical properties of the $q\geq
4$ clock model, such as emergent continuous nematic fluctuations in 3D,
critical phases with quasi-long-range nematic order in 2D, and Ashkin-Teller
nematicity in 2D.

##  I Introduction

Fingerprints typical of electronic nematic behavior have been experimentally
seen in a wide range of settings, such as unconventional cuprates [1, 2, 3, 4,
5] and iron-based superconductors [6, 7, 8, 9, 10, 11], quantum Hall systems
[12, 13], correlated oxides [14], doped topological insulators [15, 16],
twisted moiré systems [17, 18, 19, 20], $f$-electron materials [21, 22, 23,
24], kagome metals [25], colossal magnetoresistance compounds [26], triangular
antiferromagnets [27, 28, 29, 30, 31], topological semimetals [32], and
optical lattices [33]. As its name suggests, this state of matter is the
quantum analog of classical nematic liquid crystals [34, 35], as it causes the
spontaneous breaking of a discrete rotational symmetry of the system while
preserving its properties under translations [36]. At first, the term
electronic nematicity was coined to describe the ordered state that emerges
upon the partial melting of an underlying charge or spin stripe state that
restores the translational symmetry of a correlated metal or insulator, in
close analogy to the smectic-to-nematic transition of liquid crystals [36, 37,
38, 39, 40]. More recently, this concept of a so-called vestigial nematic
phase was extended to the case of multi-component superconductors, which were
shown in certain circumstances to support a partially-melted state that
restores the U(1) gauge symmetry but keeps the rotational symmetry of the
pairing state broken [41, 42]. Since the seminal work of Ref. [36], the term
electronic nematicity has been broadly employed to refer to any spontaneous
rotational-symmetry breaking phase that is driven by electronic interactions
(as opposed to elastic interactions) and that does not break additional
symmetries – for recent reviews, see [43, 44, 45, 46]. This includes the case
of weakly-interacting Fermi liquids that undergo a charge $l=2$ Pomeranchuk
instability via a Stoner-like mechanism [47, 48, 49, 50, 51], which is
particularly favored when the electronic band structure is close to a van Hove
singularity and thus has an enhanced density of states [52, 53]. Another
example of electronic nematicity under this definition are insulating and
metallic multi-orbital systems that display ferroquadrupolar order [54].
Extensions of this concept to states that break rotational and time-reversal
or inversion symmetries have also been studied [55, 56, 57, 58, 59, 60].

Figure 1: Schematics of the systematic classification performed in this paper
of electronic nematicity in crystalline and quasicrystalline systems. The
analysis starts from the ideal isotropic three-dimensional system, described
by the orthogonal group $\mathsf{O(3)}$ (Sec. II). The spherical symmetry is
systematically lowered, leading to two different sets of groups: the seven
polyhedral groups (Sec. III) and the seven classes of axial groups. The latter
are subgroups of $\mathsf{D_{\infty h}}$, which describes a two-dimensional
isotropic system. Axial groups include 27 crystallographic point groups (Sec.
IV) as well as the non-crystallographic dodecagonal, decagonal, and octagonal
point groups that describe quasicrystals (Sec. V).

The main difference between electronic nematic materials and nematic liquid
crystals is that, in the former, rotational symmetry is already explicitly
broken by the underlying lattice, whereas the latter has full rotational
symmetry. As a result, the electronic nematic order parameter always breaks a
discrete symmetry that lowers the point group of the crystal, and all nematic
collective modes in the ordered state are gapped. Historically, much of the
theoretical investigations on electronic nematicity have focused on strongly
anisotropic layered systems and “planar” nematic order parameters [37, 48, 61,
49, 50, 51, 62, 63]. Perhaps the most well recognized example is the
electronically-driven breaking of $C_{4z}$ symmetry (i.e. fourfold rotations
about the $z$-axis), in which case the electronic nematic order parameter is
Ising-like [64, 65, 66, 67, 68, 69, 70, 71, 39, 72, 73]. We emphasize that
although such a nematic transition is, on symmetry grounds, equivalent to a
tetragonal-to-orthorhombic structural transition, the microscopic mechanisms
responsible for the spontaneously broken symmetry are completely different
[9]. The fact that electronic nematic phenomena have been heavily studied in
tetragonal lattices is not surprising, since the initial experimental
observations were on tetragonal (or nearly tetragonal) compounds such as
bilayer ruthenates [14], cuprates [1], and iron pnictides [6, 7]. More
recently, experimental reports of the spontaneous breaking of $C_{3z}$
symmetry (i.e. threefold rotations about the $z$-axis) by electronic degrees
of freedom in hexagonal and trigonal systems such as doped
$\mathrm{Bi_{2}Se_{3}}$ [15, 16], twisted multi-layer graphene [17, 18, 19,
20], and triangular antiferromagnets [27, 28, 29, 30, 31] have motivated a
deeper theoretical investigation of 3-state Potts electronic nematic order
parameters [41, 74, 57, 75, 76, 77, 78, 79, 80, 81].

The diversity of systems in which electronic nematicity has been observed
indicates that a broader description and systematic classification of this
phenomenon is timely. Indeed, there is a rich landscape of crystalline
symmetries whose effects on electronic nematic order remain to be explored in
depth. This includes cubic systems such as $\mathrm{EuB_{6}}$ [26] and
$\mathrm{CaSn_{3}}$ [32], which have been reported to display anisotropic
properties consistent with nematicity. Besides crystals, quasicrystals [82,
83] have been recently found to display electronically ordered states of
interest such as superconductivity [84, 85], magnetism [86], and quantum
criticality [87], raising interesting questions about how electronic
nematicity would be manifested in systems that lack periodicity but possess
orientational order. This includes not only quasicrystal materials with
icosahedral, dodecagonal, decagonal, and octagonal symmetries [88, 89, 90,
91], but also artificial quasicrystals that can be assembled in the laboratory
by twisting two periodic crystals by specific angles [92, 93, 94, 95, 96, 85].
Moreover, besides the planar nematic order parameters studied in layered
systems, there are additional nematic channels available. Indeed, in terms of
the usual electronic operators $\hat{\psi}_{\boldsymbol{k}\sigma}$, the
nematic order parameters can be expressed as the expectation values of the
quadrupolar charge operator [48],
$d_{i}\equiv\sum\limits_{\boldsymbol{k},\sigma}\left\langle
f_{i}\left(\hat{\boldsymbol{k}}\right)\hat{\psi}_{\boldsymbol{k}\sigma}^{\dagger}\hat{\psi}_{\boldsymbol{k}\sigma}\right\rangle$.
Since there are five $d$-wave form factors
$f_{i}\left(\hat{\boldsymbol{k}}\right)$, one must consider, in the most
general case, a nematic order parameter $\boldsymbol{d}$ with five components:

$\boldsymbol{d}=\left(\begin{array}[]{c}d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}\\\
d_{x^{2}-y^{2}}\\\ d_{2yz}\\\ d_{2xz}\\\ d_{2xy}\end{array}\right),$ (1)

where the subscript denotes the corresponding $d$-wave form factor in
Cartesian coordinates. This opens the possibility of realizing “out-of-plane”
nematic order in materials that are not strongly anisotropic, such as
$d_{2xz}$ and $d_{2yz}$, beyond the widely investigated “in-plane” nematic
states corresponding to $d_{x^{2}-y^{2}}$ and $d_{2xy}$.

In this paper, we perform a systematic and thorough classification of
electronic nematic order in all crystalline and quasicrystalline point groups.
Specifically, we derive the Landau expansions and determine the universal
properties of the nematic transitions associated with the different allowed
types of in-plane and out-of-plane nematic order parameters $d_{i}$ for a
given point group. To gain further insight into how the reduced rotational
symmetry of the underlying crystalline or quasicrystalline environment
constrains nematicity, we establish a direct relationship between the five-
component charge-quadrupolar nematic order parameter $\boldsymbol{d}$ in Eq.
(1), which describes the electronic fluid, and the full rank-2 traceless
symmetric tensor $Q_{\mu\mu^{\prime}}$ that describes nematic liquid crystals
in general [34, 35]. We find it illuminating to write the latter in the
eigenstate basis:

$Q_{\mu\mu^{\prime}}=q_{1}\Big{(}n_{\mu}n_{\mu^{\prime}}-\frac{1}{3}\delta_{\mu\mu^{\prime}}\Big{)}-q_{2}\left(m_{\mu}m_{\mu^{\prime}}-l_{\mu}l_{\mu^{\prime}}\right),$
(2)

where $q_{1}$ and $q_{2}$ are scalars and $\boldsymbol{n}$, $\boldsymbol{m}$,
$\boldsymbol{l}$ form a complete set of orthonormal vectors in three-
dimensional space. When $q_{2}=0$ or $q_{2}=q_{1}$, $Q_{\mu\mu^{\prime}}$
describes a uniaxial nematic state with a uniquely defined nematic director
$\boldsymbol{n}$ or $\boldsymbol{m}$, respectively; in all other cases,
$Q_{\mu\mu^{\prime}}$ describes a biaxial nematic state [97]. While in liquid
crystals these axes determine the orientation of the constituent molecules, in
metallic systems they are manifested in the Fermi surface of the nematic
state. Thus, by establishing the relationship between the five free parameters
encompassed by $q_{1}$, $q_{2}$, $\boldsymbol{n}$, $\boldsymbol{m}$,
$\boldsymbol{l}$ in Eq. (2) and the five charge-quadrupolar order parameters
encoded in $\boldsymbol{d}$ in Eq. (1), we determine the characteristic Fermi
surface distortion patterns in the nematic states of every point group
studied.

To perform the analysis in a transparent and insightful way, we start from the
isotropic three-dimensional system, described by the orthogonal group
$\mathsf{O(3)}=\mathsf{SO(3)}\times\left\\{E,I\right\\}$, from which all
crystallographic and non-crystallographic point groups can be obtained by
systematically reducing the symmetry according to group theory. Here,
$\mathsf{SO(3)}$ contains arbitrary rotations about any three-dimensional
axis, $I$ denotes inversion, and $E$ denotes the identity operation. Fig. 1
illustrates the well-known path by which the $\mathsf{O(3)}$ symmetry is
systematically lowered to yield all point groups. First, one needs to
distinguish between two types of point groups: the polyhedral groups and the
axial groups. There are only $7$ polyhedral groups: $3$ tetrahedral and $2$
octahedral groups, which form the cubic crystal system, and $2$ icosahedral
groups, which describe a large number of quasicrystals.

In contrast, there is an infinite number of axial groups, all of which are
subgroups of $\mathsf{D_{\infty h}=SO(2)}\times\\{E,I\\}\times\\{E,C_{2x}\\}$,
which describes a two-dimensional isotropic system and contains arbitrary
rotations about the $z$-axis, twofold rotations about any in-plane axis,
inversion, and reflection with respect to the plane. The axial groups are
further subdivided into the seven different groups outlined in Fig. 1:
$\mathsf{D_{nh}}$, $\mathsf{D_{nd}}$, $\mathsf{C_{nh}}$, $\mathsf{D_{n}}$,
$\mathsf{S_{2n}}$, $\mathsf{C_{nv}}$, and $\mathsf{C_{n}}$, which contain
discrete $n$-fold rotations around the axis perpendicular to the plane and, in
some cases, twofold rotations about an axis parallel to the plane. Enforcing
the crystallographic restriction theorem leads to $27$ crystallographic axial
point groups, which are then divided in six crystal systems: hexagonal ($7$
groups), trigonal ($5$ groups), tetragonal ($7$ groups), orthorhombic ($3$
groups), monoclinic ($3$ groups), and triclinic ($2$ groups). Together with
the $5$ polyhedral cubic groups, they constitute the $32$ crystallographic
point groups and $7$ crystal systems. Importantly, the non-crystallographic
point groups can still describe materials that are non-periodic and that can
host non-trivial electronic states [88, 89]. Thus, in this paper, we also
focus on the non-crystallographic point groups that describe quasicrystals,
namely, the polyhedral icosahedral groups as well as the axial dodecagonal,
decagonal, and octagonal point groups.

Of course, several of our results for the character of the nematic transition
in certain point groups recover well-established results obtained in previous
symmetry analyses of two-dimensional nematics [43, 44, 46], structural and
ferroelastic transitions [98, 99, 100], multipolar order [56, 101], and
intertwined orders [102]. We emphasize that our main goal here is to provide a
complete, self-contained classification of electronic nematicity in crystals
_and_ quasicrystals by systematically determining how the removal of
symmetries impacts the general structure of the tensorial nematic order
parameter (2) via its relationship with the five-component charge-quadrupolar
order parameter (1). Our results are summarized in Tables 1 and 2 in Section
VI, and their derivations are shown in Sections II-V and the Appendices D-G.
For the convenience of the reader not interested in the technical details of
our analysis, we outline here the main results.

  * •

The icosahedral quasicrystalline point groups are the only ones for which the
five components of the nematic order parameter (1) transform together as a
single five-dimensional irreducible representation. The nematic transition is
first-order within mean-field due to a cubic invariant in the Landau
expansion. Depending on the signs of the Landau coefficients, the nematic
director that characterizes the uniaxial nematic ground state aligns itself
with either the $10$ axes of threefold rotational symmetry or the $6$ axes of
fivefold rotational symmetry of the quasicrystal. In the former case, the
residual point group is crystalline (trigonal) whereas in the latter, it is
quasicrystalline (decagonal).

  * •

The components of the nematic order parameter $\boldsymbol{d}$ in the cubic
crystal groups split into separate triplet and doublet order parameters. The
triplet behaves as a $Z_{4}$-Potts order parameter, giving rise to a uniaxial
nematic state in which the director is parallel to one of four axes of
threefold rotational symmetry of the crystal (i.e. the space diagonals). The
doublet behaves as a $Z_{3}$-clock/Potts order parameter, and the three
nematic axes point along the coordinate axes. In the cubic groups for which
the coordinate axes are also axes of fourfold rotational symmetry, the nematic
ground state is uniaxial, otherwise, it is biaxial. The nematic transitions
are, again, first-order within mean-field.

  * •

In any axial group, the nematic component
$d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}$ is necessarily non-zero, which
changes the shape of the Fermi surface in the symmetry-unbroken phase from
spherical to cylindrical. The other four components are decomposed into at
least two independent nematic doublets, the in-plane
$\boldsymbol{d}^{\mathrm{ip}}=\big{(}d_{x^{2}-y^{2}},\,d_{2xy}\big{)}^{T}$ and
the out-of-plane
$\boldsymbol{d}^{\mathrm{op}}=\left(d_{2yz},\,d_{2xz}\right)^{T}$. The nematic
ground state is always biaxial for the axial groups, unless the Landau
parameters are fine tuned. In most cases, at least one of the nematic axes
aligns with a high-symmetry in-plane direction.

  * •

Each nematic doublet in the axial groups behaves as a $Z_{q}$-clock order
parameter. If the group has $\left(2n\right)$-fold rotational symmetry, $q=n$
for $\boldsymbol{d}^{\mathrm{ip}}$ and $q=2n$ for
$\boldsymbol{d}^{\mathrm{op}}$. If the group has $\left(2n+1\right)$-fold
rotational symmetry, $q=2n+1$ for both $\boldsymbol{d}^{\mathrm{ip}}$ and
$\boldsymbol{d}^{\mathrm{op}}$. In three dimensions, the $q$-state clock model
is known to undergo an XY transition for $q\geq 4$, similarly to the isotropic
2D nematic case. In two dimensions, the $Z_{q}$-clock model displays an
intermediate critical phase with quasi-long-range order for $q\geq 5$ and an
Ashkin-Teller phase transition with non-universal critical exponents for
$q=4$. A first-order transition can only occur for $q=3$ and above the upper
critical dimension $d_{u}\apprge 2$. If the axial group lacks in-plane
rotational symmetry axes, or if $\boldsymbol{d}^{\mathrm{ip}}$ and
$\boldsymbol{d}^{\mathrm{op}}$ have the same transformation properties, the
clock term in the Landau expansion acquires a non-universal offset, which we
denote as the $Z_{q}^{*}$-clock model.

  * •

The main difference between hexagonal and trigonal crystals is that, in the
latter, $\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ are
not independent, which has a significant impact on the character of the
nematic transition. In tetragonal crystals, $\boldsymbol{d}^{\mathrm{ip}}$ is
further decomposed into two one-component Ising-like order parameters.
Orthorhombic, monoclinic, and triclinic crystals are described by Abelian
point groups, which only admit one-dimensional irreducible representations. As
a result, any nematic transition must be Ising-like.

  * •

In dodecagonal, decagonal, and octagonal quasicrystals, both nematic doublets
behave as $Z_{q}$-clock order parameters with $q\geq 4$. This makes axial
quasicrystals interesting platforms to realize exotic nematicity enabled by
the critical properties of the $q\geq 4$ clock model, such as emergent XY
nematic fluctuations and critical nematic phases displaying quasi-long-range
order.

The organization of the paper, schematically shown in Fig. 1, is as follows:
Sec. II introduces the formalism and presents the results for nematicity in
the three-dimensional isotropic system. The case of polyhedral groups, i.e.,
cubic crystal system and icosahedral quasicrystals, is discussed in Sec. III.
Sec. IV presents the properties of electronic nematic order in the isotropic
two-dimensional system and in the axial crystallographic groups (hexagonal,
trigonal, tetragonal, orthorhombic, monoclinic, and triclinic crystal
systems). The cases of dodecagonal, decagonal, and octagonal quasicrystalline
axial groups are presented in Sec. V. Sec. VI contains our concluding remarks
as well as a summary of our results for the crystalline point groups in Table
1 and for the quasicrystalline point groups in Table 2. Appendix A contains
additional details about isotropic nematicity. Appendix B explains the concept
of symmetrized decomposition of products of irreducible representations,
whereas Appendix C shows details of the minimization of the nematic Landau
expansion in icosahedral systems. Appendix D presents details of electronic
nematic order in cubic crystals without fourfold rotational symmetry axes.
Appendices E, F, and G derive the properties of electronic nematicity,
respectively, in hexagonal, trigonal, and tetragonal point groups that lack
in-plane rotational symmetry axes.

##  II Electronic nematicity in a 3D isotropic system

###  II.1 Representations of the nematic order parameter

To set the stage for the remainder of the paper, we first review the
properties of electronic nematicity in 3D isotropic systems, i.e. systems that
are fully rotational invariant. In classical nematics, the order parameter is
a “headless vector” associated with the orientation of the elongated molecules
that form the liquid crystal. It is conveniently expressed in terms of the
symmetric traceless tensor $Q_{\mu\mu^{\prime}}=\left\langle
a_{\mu}a_{\mu^{\prime}}-\frac{1}{3}\delta_{\mu\mu^{\prime}}\boldsymbol{a}^{2}\right\rangle$,
with $\mu,\mu^{\prime}=1,2,3$ and the director
$\boldsymbol{a}=(a_{1},a_{2},a_{3})$ [34, 35]. The generalization to the
quantum (i.e. electronic) case has been widely discussed in the literature
[48, 97]. In terms of the electronic annihilation (creation) operators
$\hat{\psi}_{\boldsymbol{k}}$ ($\hat{\psi}_{\boldsymbol{k}}^{\dagger}$), the
tensorial order parameter is given by

$\displaystyle Q_{\mu\mu^{\prime}}$
$\displaystyle=\frac{1}{k_{F}^{2}}\sum_{\boldsymbol{k},\sigma}\left\langle\mathcal{F}_{\mu\mu^{\prime}}\left(\boldsymbol{k}\right)\hat{\psi}_{\boldsymbol{k}\sigma}^{\dagger}\hat{\psi}_{\boldsymbol{k}\sigma}\right\rangle,$
(3)

with momentum $\boldsymbol{k}$, Fermi momentum $k_{F}$, spin
$\sigma=\uparrow,\downarrow$, and the $d$-wave form factor
$\mathcal{F}_{\mu\mu^{\prime}}\left(\boldsymbol{k}\right)=2\big{(}k_{\mu}k_{\mu^{\prime}}-\frac{1}{3}\delta_{\mu\mu^{\prime}}\boldsymbol{k}^{2}\big{)}$.
Generalizations to multi-orbital systems are straightforward, but will not be
covered here [71, 101]. Thus, the components of the electronic nematic order
parameter correspond to quadrupolar charge order which, in the theory of
interacting Fermi liquids, corresponds to an angular momentum $l=2$
Pomeranchuk instability of the Fermi liquid in the singlet channel [47]. The
condensation of the tensor order parameter $Q$ leads to a distortion of the
otherwise spherically symmetric Fermi surface, thus breaking the rotational
invariance of the isotropic system:

$\displaystyle\epsilon_{\boldsymbol{k}}$
$\displaystyle=\frac{\boldsymbol{k}^{2}}{2m}-\mu_{0}+\frac{1}{2}\,\mathrm{tr}\left[Q\mathcal{F}\left(\boldsymbol{k}\right)\right].$
(4)

Here, $m$ is the effective electron mass and $\mu_{0}$, the chemical
potential. For concreteness, and to better visualize the fingerprints of
nematicity on the electronic degrees of freedom, our analysis will focus on
metals. Of course, the symmetry properties of the nematic state would be the
same in insulators.

Being a traceless symmetric tensor, $Q$ has five independent components. Thus,
it is also convenient to express nematic phenomena in terms of a five-
component “vector”
$\boldsymbol{d}=\left(d_{1},d_{2},d_{3},d_{4},d_{5}\right)$. To derive the
symmetry properties of this vector, we note that a 3D isotropic system is
described by the continuous (orthogonal) group
$\mathrm{O(3)}=\mathrm{SO(3)}\times\left\\{E,I\right\\}$, which combines the
full rotation group $\mathrm{SO(3)}$ with the inversion operation $I$ (here,
$E$ denotes the identity operation). Thus, in group-theory notation, the five-
component nematic vector $\boldsymbol{d}$ transforms according to the five-
dimensional irreducible representation (IR) $\Gamma_{j=2}^{+}$ of the
orthogonal group $\mathrm{O(3)}$. Note that the superscript $+$ ($-$)
indicates an inversion-even (inversion-odd) IR.

Such a five-component order parameter can be obtained in a straightforward way
as a bilinear $\boldsymbol{M}$ constructed from the coordinate vector
$\boldsymbol{r}=(x,y,z)$, which in turn transforms according to the three-
dimensional vector IR $\Gamma_{j=1}^{-}$. That such a bilinear
$\boldsymbol{M}$ in the $\Gamma_{j=2}^{+}$-channel exists follows from the
product decomposition
$\Gamma_{j=1}^{-}\otimes\Gamma_{j=1}^{-}=\Gamma_{j=0}^{+}\oplus\Gamma_{j=2}^{+}\oplus\Gamma_{j=1}^{+}$.
Writing the bilinear components as
$M_{i}=r_{\mu}\lambda_{\mu\mu^{\prime}}^{j=2,i}r_{\mu^{\prime}}$ with
$i\in\\{1,\dots,5\\}$ and $\mu\in\\{1,2,3\\}$, we obtain the five matrices
$\boldsymbol{\lambda}^{j=2}$ from the transformation condition

$\displaystyle\mathcal{R}_{-,j=1}^{T}(g)\lambda^{j=2,i}\mathcal{R}_{-,j=1}(g)$
$\displaystyle=\mathcal{R}_{+,j=2}(g)_{ii^{\prime}}\lambda^{j=2,i^{\prime}},$
(5)

where $\mathcal{R}_{\pm,j}(g)$ denotes the $5\times 5$ transformation matrix
of a symmetry element $g$ associated with the IR $\Gamma_{j}^{\pm}$. The
symmetry elements
$g=\big{(}\vartheta,\hat{\boldsymbol{\ell}},\mathcal{I}\big{)}$ are
parametrized in terms of the rotation angle $\vartheta$ around the unit
rotation axis $\hat{\boldsymbol{\ell}}$ and the index $\mathcal{I}=\pm 1$ for
inversion being applied ($\mathcal{I}=-1$) or not ($\mathcal{I}=+1$).
Therefore, the transformation matrices are given by
$\mathcal{R}_{+,j}(g)=\mathcal{R}_{j}(\vartheta,\hat{\boldsymbol{\ell}})$ and
$\mathcal{R}_{-,j}(g)=\mathcal{I}\mathcal{R}_{j}(\vartheta,\hat{\boldsymbol{\ell}})$
where
$\mathcal{R}_{j}(\vartheta,\hat{\boldsymbol{\ell}})=\exp\left(-\mathsf{i}\vartheta\,\boldsymbol{J}^{(j)}\cdot\hat{\boldsymbol{\ell}}\right)$
are the well-known $(2j+1)\times(2j+1)$-dimensional rotation matrices. Using
the $\boldsymbol{J}^{(j)}$ matrices outlined in Appendix A, one finds the five
matrices $\boldsymbol{\lambda}^{j=2}$ in (5) to be identical to the five
symmetric Gell-Mann matrices,

$\displaystyle\boldsymbol{\lambda}^{j=2}$
$\displaystyle=\left\\{\frac{1}{\sqrt{3}}\left(\begin{smallmatrix}\text{-}1&0&0\\\
0&\text{-}1&0\\\
0&0&2\end{smallmatrix}\right),\left(\begin{smallmatrix}1&0&0\\\
0&\text{-}1&0\\\
0&0&0\end{smallmatrix}\right),\left(\begin{smallmatrix}0&0&0\\\ 0&0&1\\\
0&1&0\end{smallmatrix}\right),\left(\begin{smallmatrix}0&0&1\\\ 0&0&0\\\
1&0&0\end{smallmatrix}\right),\left(\begin{smallmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{smallmatrix}\right)\right\\}.$ (6)

We can now label the five nematic components of $\boldsymbol{d}$ as
polynomials of the coordinate vector $\boldsymbol{r}=(x,y,z)$, according to
their symmetry properties encoded in the bilinear $\boldsymbol{M}$:

$\displaystyle\boldsymbol{d}$ $\displaystyle=\left(\begin{array}[]{c}d_{1}\\\
d_{2}\\\ d_{3}\\\ d_{4}\\\
d_{5}\end{array}\right)=\left(\begin{array}[]{c}d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}\\\
d_{x^{2}-y^{2}}\\\ d_{2yz}\\\ d_{2xz}\\\ d_{2xy}\end{array}\right).$ (17)

Unsurprisingly, the five components of $\boldsymbol{d}$ correspond to the five
$d$-wave form factors written as tesseral harmonics. The symmetric Gell-Mann
matrices (6) also conveniently establish the relationship between the tensor
notation $Q$ and the vector notation $\boldsymbol{d}$:

$\displaystyle Q$
$\displaystyle=\boldsymbol{d}\cdot\boldsymbol{\lambda}^{j=2},$
$\displaystyle\boldsymbol{d}$
$\displaystyle=\mathrm{tr}\left[Q\boldsymbol{\lambda}^{j=2}\right]\,\big{/}2,$
(18)

which, in explicit form, gives:

$\displaystyle Q$
$\displaystyle=\left(\begin{smallmatrix}d_{2}-\frac{1}{\sqrt{3}}d_{1}&d_{5}&d_{4}\\\
d_{5}&-d_{2}-\frac{1}{\sqrt{3}}d_{1}&d_{3}\\\
d_{4}&d_{3}&\frac{2}{\sqrt{3}}d_{1}\end{smallmatrix}\right).$ (19)

In this work, we interchangeably use the tensor notation $Q$ (19) and the
five-component vector notation $\boldsymbol{d}$ (17). Moreover, throughout
this work, we use the following transformation relation for the nematic order
parameter

$\displaystyle\boldsymbol{d}$
$\displaystyle\overset{g}{\longrightarrow}\mathcal{R}_{+,j=2}(g)\,\boldsymbol{d},$
(20)

such that the transformation matrices are fixed by the set of given symmetry
elements $g$. This is particularly useful as it uniquely defines the multi-
component nematic channels within reduced symmetry systems where
$\mathcal{R}_{+,j=2}(g)$ is block-diagonal, such as crystal lattices.

A convenient representation of the nematic tensor $Q$ (19) is in terms of its
eigenbasis. In this representation, which we denote the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation, the nematic
tensor becomes

$\displaystyle Q_{\mu\mu^{\prime}}$
$\displaystyle=\sqrt{3}|\boldsymbol{d}|\cos\alpha\Big{(}n_{\mu}n_{\mu^{\prime}}-\frac{1}{3}\delta_{\mu\mu^{\prime}}\Big{)}$
$\displaystyle-|\boldsymbol{d}|\sin\alpha\left(m_{\mu}m_{\mu^{\prime}}-l_{\mu}l_{\mu^{\prime}}\right),$
(21)

where the angle $\alpha$ determines the eigenvalues while the unit vectors
$\boldsymbol{n}$, $\boldsymbol{m}$, $\boldsymbol{l}$ span the eigenspace of
the matrix (19). As a set of orthonormal eigenvectors, $\boldsymbol{n}$,
$\boldsymbol{m}$, $\boldsymbol{l}$ obey the orthogonality relations
$\boldsymbol{n}\cdot\boldsymbol{m}=\boldsymbol{n}\cdot\boldsymbol{l}=\boldsymbol{m}\cdot\boldsymbol{l}=0$,
as well as the completeness relation

$\displaystyle
n_{\mu}n_{\mu^{\prime}}+m_{\mu}m_{\mu^{\prime}}+l_{\mu}l_{\mu^{\prime}}$
$\displaystyle=\delta_{\mu\mu^{\prime}},$ (22)

valid for any $\mu,\mu^{\prime}$. One reason why the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation is particularly
useful is because the tensor (21) is readily diagonalized using the orthogonal
matrix $U=\left(\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right)$,

$\displaystyle Q^{d}$
$\displaystyle=U^{T}QU=\frac{2|\boldsymbol{d}|}{\sqrt{3}}\left(\begin{smallmatrix}\cos\left(\alpha\right)&0&0\\\
0&\cos\left(\alpha+\frac{2\pi}{3}\right)&0\\\
0&0&\cos\left(\alpha+\frac{4\pi}{3}\right)\end{smallmatrix}\right).$ (23)

Moreover, in this representation, we can use Eq. (18) to write the five-
component vector $\boldsymbol{d}$ (17) as

$\displaystyle\boldsymbol{d}$
$\displaystyle=\frac{|\boldsymbol{d}|}{\sqrt{3}}\Big{\\{}\cos\left(\alpha\right)\left(n_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}n_{\mu^{\prime}}\right)+\cos\left(\alpha+\frac{2\pi}{3}\right)\left(m_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}m_{\mu^{\prime}}\right)$
$\displaystyle\quad+\cos\left(\alpha+\frac{4\pi}{3}\right)\left(l_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}l_{\mu^{\prime}}\right)\Big{\\}},$
(24)

where summation over $\mu$,$\mu^{\prime}$ is implied. The parameterization of
the orthonormal eigenvectors
$\boldsymbol{n}$,$\boldsymbol{m}$,$\boldsymbol{l}$ involves three angles and
it is not unique, see Appendix A for details. Importantly, to ensure a one-to-
one mapping between Eq. (19) and Eq. (21), the eigenvalue-angle $\alpha$ must
be restricted to the range $[0,\pi/3]$. Plotting the three eigenvalues as a
function of $\alpha$ in Fig. 2(a) makes it clear that the smallest eigenvalue
is always related to $\boldsymbol{m}$ whereas the largest one is associated
with $\boldsymbol{n}$, see also Eq. (24). Moreover, the figure also reveals
two special points, $\alpha=0,\,\pi/3$, for which two eigenvalues are
degenerate – either the smallest eigenvalue (for $\alpha=0$) or the largest
eigenvalue (for $\alpha=\pi/3$), corresponding to an uniaxial nematic state.
For later convenience, we note that in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation of Eq. (24), a
sign change in $\boldsymbol{d}$ corresponds to changing
$\alpha\rightarrow\frac{\pi}{3}-\alpha$ and swapping
$\boldsymbol{n}\leftrightarrow\boldsymbol{m}$:

$\displaystyle-\boldsymbol{d}\Big{[}\\!\left|\boldsymbol{d}\right|,\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\Big{]}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}\right|,\frac{\pi}{3}-\alpha,\boldsymbol{m},\boldsymbol{n},\boldsymbol{l}\right].$
(25)

###  II.2 Minimization of the nematic free energy

The value of the angle $\alpha$ that determines the eigenvalues of the nematic
order parameter $Q$ or $\boldsymbol{d}$, see Eqs. (23) and (24), can be
obtained by minimizing the corresponding Landau expansion. The symmetry-
allowed terms in the Landau expansion can be obtained by decomposing the
products of the non-trivial IR that defines the order parameter and then
picking the terms in the decomposition that transform trivially under the
group operations. One way to unambiguously determine these terms, which will
be very useful once we consider point groups, is through the so-called
decomposition of the symmetrized product, see Appendix B for details. This
special decomposition removes any redundancy related to the anti-symmetric
channels and avoids the double-counting that one would encounter by
considering the non-symmetrized product. With the order parameter
$\boldsymbol{d}$ transforming according to the IR $\Gamma_{j=2}^{+}$, the
decomposition of the symmetrized products for each Landau expansion order
becomes

$\displaystyle\big{[}\otimes_{l=1}^{2}\Gamma_{2}^{+}\big{]}_{s}$
$\displaystyle=\Gamma_{0}^{+}\oplus\Gamma_{2}^{+}\oplus\Gamma_{4}^{+},$ (26)
$\displaystyle\big{[}\otimes_{l=1}^{3}\Gamma_{2}^{+}\big{]}_{s}$
$\displaystyle=\Gamma_{0}^{+}\oplus\Gamma_{2}^{+}\oplus\Gamma_{3}^{+}\oplus\Gamma_{4}^{+}\oplus\Gamma_{6}^{+},$
(27) $\displaystyle\big{[}\otimes_{l=1}^{4}\Gamma_{2}^{+}\big{]}_{s}$
$\displaystyle=\Gamma_{0}^{+}\oplus 2\Gamma_{2}^{+}\oplus
2\Gamma_{4}^{+}\oplus\Gamma_{5}^{+}\oplus\Gamma_{6}^{+}\oplus\Gamma_{8}^{+}.$
(28)

Here, we use the abbreviated tensor product notation, e.g.
$\otimes_{l=1}^{3}\Gamma_{2}^{+}=\Gamma_{2}^{+}\otimes\Gamma_{2}^{+}\otimes\Gamma_{2}^{+}$,
and the subscript $s$ to indicate that only the symmetrized product is
considered. The symmetrized decompositions (26)-(28) imply the existence of
one Landau invariant per expansion order, since there is exactly one trivial
channel ($\Gamma_{0}^{+}$) per expansion order. To systematically identify
these invariants, it is useful to first determine the bilinears associated
with each IR in Eq. (26), which we denote by $D^{j=0}=|\boldsymbol{d}|^{2}$
($\Gamma_{0}^{+}$), $\boldsymbol{D}^{j=2}$ ($\Gamma_{2}^{+}$) and
$\boldsymbol{D}^{j=4}$ ($\Gamma_{4}^{+}$). Solving the respective
transformation conditions, similar to Eq. (5), leads to the five-component
vector

$\displaystyle\boldsymbol{D}^{j=2}$
$\displaystyle=\left(\begin{smallmatrix}d_{1}^{2}-d_{2}^{2}+\frac{1}{2}\left(d_{3}^{2}+d_{4}^{2}-2d_{5}^{2}\right)\\\
-2d_{1}d_{2}-\frac{\sqrt{3}}{2}\left(d_{3}^{2}-d_{4}^{2}\right)\\\
\left(d_{1}-\sqrt{3}d_{2}\right)d_{3}+\sqrt{3}d_{4}d_{5}\\\
\left(d_{1}+\sqrt{3}d_{2}\right)d_{4}+\sqrt{3}d_{3}d_{5}\\\
-2d_{1}d_{5}+\sqrt{3}d_{3}d_{4}\end{smallmatrix}\right),$ (29)

and the nine-component vector $\boldsymbol{D}^{j=4}$ given in Eq. (267) in
Appendix A. Since $\boldsymbol{D}^{j=2}$ is quadratic in $d_{i}$, it is
straightforward to construct the cubic invariant that appears in the
symmetrized decomposition in Eq. (27). Specifically, since $\boldsymbol{d}$
and $\boldsymbol{D}^{j=2}$ transform as the same IR $\Gamma_{j=2}^{+}$, their
scalar product must transform trivially and thus appear in the free-energy
expansion:

$\displaystyle\boldsymbol{d}\cdot\boldsymbol{D}^{j=2}$
$\displaystyle=\frac{\sqrt{3}}{2}\mathrm{tr}\left[Q^{3}\right]=\left|\boldsymbol{d}\right|^{3}\cos\left(3\alpha\right).$
(30)

In the last step, we employed the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation (21)-(24) to
further simplify the expression. As for the quartic decomposition in Eq. (28),
since there is only one term in the symmetrized product that transforms
trivially, we can readily identify it as
$D^{j=0}D^{j=0}=|\boldsymbol{d}|^{4}$. Therefore, the resulting nematic Landau
expansion is given by the action

$\displaystyle\mathcal{S}\left[\left|\boldsymbol{d}\right|,\alpha\right]$
$\displaystyle=\int_{\mathsf{x}}\left\\{r_{0}\left|\boldsymbol{d}\right|^{2}+g\left|\boldsymbol{d}\right|^{3}\cos\left(3\alpha\right)+u\left|\boldsymbol{d}\right|^{4}\right\\},$
(31)

where $\mathsf{x}=(\boldsymbol{r},\tau)$ comprises position and imaginary
time, $\int_{\mathsf{x}}\equiv\int d^{3}r\int_{0}^{1/T}d\tau$, $r_{0}$ is the
control parameter that tunes the system across a nematic transition, and $g$,
$u$ are cubic and quartic Landau parameters, respectively. We opted to
represent the Landau expansion in terms of an action rather than a free energy
in order to explicitly account for the temporal dependence of the nematic
order parameter, which is necessary in the case of a quantum phase transition.
For a thermal transition, we can write $r_{0}=a_{0}(T-T_{0})$ with $a_{0}>0$
and $T_{0}$ a reference temperature, such that the free energy is given by
$F=ST_{0}$.

We note that the action (31) does not depend on the orientation of the nematic
axes $\boldsymbol{n}$, $\boldsymbol{m}$, and $\boldsymbol{l}$, which is a
manifestation of the full rotational invariance of the 3D isotropic system.
Although the expansion (31) resembles that of a 3-state Potts/clock model
($Z_{3}$ model), it is important to emphasize that
$\alpha\in[0,\frac{\pi}{3}]$ and that $\boldsymbol{d}$ is a five-component
vector. Minimization with respect to $\alpha$ gives the mean-field ground-
state angle:

$\displaystyle\alpha_{0}$
$\displaystyle=\frac{\pi}{3}\left(\frac{1+\mathrm{sign}\left(g\right)}{2}\right),$
(32)

i.e. the angle is either $\alpha_{0}=0$ (for $g<0$) or $\alpha_{0}=\pi/3$ (for
$g>0$). In either case, the effective action in terms of
$\left|\boldsymbol{d}\right|$ alone assumes the form:

$\mathcal{S}_{\mathrm{eff}}\left[\left|\boldsymbol{d}\right|\right]=\int_{\mathsf{x}}\left\\{r_{0}\left|\boldsymbol{d}\right|^{2}-\left|g\right|\left|\boldsymbol{d}\right|^{3}+u\left|\boldsymbol{d}\right|^{4}\right\\}.$
(33)

The existence of a negative cubic term implies that, within a mean-field
solution, the isotropic nematic transition is first order. Another important
property of the $\alpha_{0}=0,\,\pi/3$ solutions is that they correspond to
uniaxial nematic states, i.e., states that only depend on one eigenvector and
for which two eigenvalues are degenerate [97]. This can be seen directly from
Eq. (21): when $\alpha=0$, the second term vanishes and we obtain

$Q_{\mu\mu^{\prime}}\left[\alpha=0\right]=\sqrt{3}\left|\boldsymbol{d}\right|\Big{(}n_{\mu}n_{\mu^{\prime}}-\frac{1}{3}\delta_{\mu\mu^{\prime}}\Big{)},$
(34)

which corresponds to a nematic director along $\boldsymbol{n}$.
Correspondingly, for $\alpha=\pi/3$, we have
$\sqrt{3}\cos\alpha=\sin\alpha=\sqrt{3}/2$; using the completeness relation
(22) we find:

$Q_{\mu\mu^{\prime}}\left[\alpha=\pi/3\right]=-\sqrt{3}\left|\boldsymbol{d}\right|\Big{(}m_{\mu}m_{\mu^{\prime}}-\frac{1}{3}\delta_{\mu\mu^{\prime}}\Big{)}$
(35)

corresponding to a nematic director along $\boldsymbol{m}$. Similarly,
inserting the two $\alpha_{0}$ values into the nematic order parameter (24)
gives:

$\displaystyle\boldsymbol{d}\left[\alpha=0\right]$
$\displaystyle=\big{(}\sqrt{3}\big{/}2\big{)}\,|\boldsymbol{d}|\;n_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}n_{\mu^{\prime}},$
(36) $\displaystyle\boldsymbol{d}\left[\alpha=\pi/3\right]$
$\displaystyle=-\big{(}\sqrt{3}\big{/}2\big{)}\,|\boldsymbol{d}|\;m_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}m_{\mu^{\prime}},$
(37)

where we used, once again, the completeness relation (22) as well as the fact
that $\mathrm{tr}\left(\lambda^{j=2,i}\right)=0$. Note that (37) is just the
negative of (36), in accordance with (25). Clearly, the continuous rotational
symmetry of the isotropic system is spontaneously broken in the nematic ground
state, as the nematic director $\boldsymbol{n}$ (or $\boldsymbol{m}$) can
point in any direction.

In the case of a metallic system, these distortion patterns are manifested in
the Fermi surface of the nematic state, and the nematic instability is nothing
but an $l=2$ Pomeranchuk instability in the charge channel. The corresponding
electronic dispersion (4) can be conveniently rewritten in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation. Using Eq. (23),
we find

$\epsilon_{\boldsymbol{k}}=\frac{\boldsymbol{k}^{2}}{2m}-\mu_{0}+\frac{1}{2}\,\mathrm{tr}\left[Q_{d}U^{T}\mathcal{F}\left(\boldsymbol{k}\right)U\right].$
(38)

The transformed matrix $U^{T}\mathcal{F}\left(\boldsymbol{k}\right)U$ is given
by:

$\displaystyle U^{T}\mathcal{F}\left(\boldsymbol{k}\right)U$
$\displaystyle=\left(\begin{smallmatrix}\boldsymbol{n}^{T}\\\
\boldsymbol{m}^{T}\\\
\boldsymbol{l}^{T}\end{smallmatrix}\right)\mathcal{F}\left(\boldsymbol{k}\right)\left(\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right)$
$\displaystyle=\left(\begin{smallmatrix}\boldsymbol{n}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{n}&\boldsymbol{n}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{m}&\boldsymbol{n}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{l}\\\
\boldsymbol{m}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{n}&\boldsymbol{m}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{m}&\boldsymbol{m}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{l}\\\
\boldsymbol{l}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{n}&\boldsymbol{l}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{m}&\boldsymbol{l}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{l}\end{smallmatrix}\right).$
(39)

Using the fact that

$\boldsymbol{v}^{T}\mathcal{F}\left(\boldsymbol{k}\right)\boldsymbol{w}=2\big{(}\boldsymbol{k}\cdot\boldsymbol{v}\big{)}\big{(}\boldsymbol{k}\cdot\boldsymbol{w}\big{)}-\frac{2}{3}\boldsymbol{k}^{2}\big{(}\boldsymbol{v}\cdot\boldsymbol{w}\big{)},$
(40)

and performing the matrix product, we find:

$\displaystyle\epsilon_{\boldsymbol{k}}$
$\displaystyle=\frac{\boldsymbol{k}^{2}}{2m}-\mu_{0}+\frac{2|\boldsymbol{d}|\boldsymbol{k}^{2}}{\sqrt{3}}\Big{[}\big{(}\hat{\boldsymbol{k}}\cdot\boldsymbol{n}\big{)}^{2}\cos\left(\alpha\right)$
$\displaystyle+\big{(}\hat{\boldsymbol{k}}\cdot\boldsymbol{m}\big{)}^{2}\cos\Big{(}\alpha+\frac{2\pi}{3}\Big{)}+\big{(}\hat{\boldsymbol{k}}\cdot\boldsymbol{l}\big{)}^{2}\cos\Big{(}\alpha+\frac{4\pi}{3}\Big{)}\Big{]},$
(41)

with $\hat{\boldsymbol{k}}=\boldsymbol{k}/\left|\boldsymbol{k}\right|$. Here,
we used $\sum_{\nu=0}^{2}\cos\left(\alpha+\nu\frac{2\pi}{3}\right)=0$.

The Fermi wave-vector of the dispersion (41) can be readily obtained as a
function of $\hat{\boldsymbol{k}}$:

$\displaystyle k_{F}\left(\hat{\boldsymbol{k}}\right)$
$\displaystyle=k_{F,0}\,\Big{\\{}1+\frac{4m|\boldsymbol{d}|}{\sqrt{3}}\Big{[}\big{(}\hat{\boldsymbol{k}}\cdot\boldsymbol{m}\big{)}^{2}\cos\Big{(}\alpha+\frac{2\pi}{3}\Big{)}$
$\displaystyle+\big{(}\hat{\boldsymbol{k}}\cdot\boldsymbol{l}\big{)}^{2}\cos\Big{(}\alpha+\frac{4\pi}{3}\Big{)}+\big{(}\hat{\boldsymbol{k}}\cdot\boldsymbol{n}\big{)}^{2}\cos\left(\alpha\right)\Big{]}\Big{\\}}^{-\frac{1}{2}},$
(42)

where we defined the isotropic Fermi momentum $k_{F,0}=\sqrt{2m\mu_{0}}$.
Using the results of Fig. 2(a), we conclude that the Fermi wave-vector is
longest along the $\boldsymbol{m}$ axis,
$k_{F,\mathrm{max}}=k_{F}\big{(}\hat{\boldsymbol{k}}=\pm\boldsymbol{m}\big{)}$,
and shortest along the $\boldsymbol{n}$ directions,
$k_{F,\mathrm{min}}=k_{F}\big{(}\hat{\boldsymbol{k}}=\pm\boldsymbol{n}\big{)}$.
Therefore, we identify $\boldsymbol{m}$ and $\boldsymbol{n}$ as the long and
short nematic axes of the distorted Fermi surface, respectively. Moreover,
when $\alpha=0$ ($g<0$), the Fermi momentum along the $\boldsymbol{l}$
direction is as large as the Fermi momentum along $\boldsymbol{m}$, resulting
in a Fermi surface with the shape of an oblate spheroid; we dub this a
“compressive” nematic deformation, and associate the nematic director to
$\boldsymbol{n}$. In contrast, when $\alpha=\pi/3$ ($g>0$), the Fermi momenta
along $\boldsymbol{l}$ and along the short axes $\boldsymbol{n}$ are
equivalent, resulting in a Fermi surface with the shape of a prolate spheroid,
which we associate with a “tensile” nematic deformation. The nematic director
in this case is parallel to $\boldsymbol{m}$.

We plot the Fermi surface (42) associated with the two mean-field nematic
ground-states $\alpha=\pi/3,\,0$ in Figs. 2(c)-(d) with nematic magnitude
$|\boldsymbol{d}|=1/(3m)$, together with the undistorted Fermi surface in Fig.
2(b). The same nematic magnitude is employed in all figures in this work.
Additionally, we also plot the corresponding eigenvectors $\boldsymbol{n}$,
$\boldsymbol{m}$, $\boldsymbol{l}$ rescaled by $2.25\,k_{F}(\boldsymbol{n})$,
$2.25\,k_{F}(\boldsymbol{m})$, $2.25\,k_{F}(\boldsymbol{l})$, respectively, to
better visualize the long and short axes in each case. Regardless of the value
of $g$, the resulting Fermi surface is always a uniaxial ellipsoid, which has
the shape of either an oblate spheroid ($\alpha=0$ , $g<0$, resulting in the
compressive distortion of panel (d)) or a prolate spheroid ($\alpha=\pi/3$ ,
$g>0$, resulting in the tensile distortion of panel (c)). Importantly, the
uniaxial vector $\boldsymbol{n}$ or $\boldsymbol{m}$ can point in any
direction, reflecting the spontaneous breaking of the continuous rotational
symmetry below the nematic transition.

Figure 2: (a) Eigenvalues of the tensorial nematic order parameter $Q$ (23),
in units of $2|\boldsymbol{d}|/\sqrt{3}$, as a function of the parameter
$\alpha$. The eigenvalues are colored according to their corresponding
eigenvector $\boldsymbol{n}$ (blue), $\boldsymbol{m}$ (green), and
$\boldsymbol{l}$ (red), see Eq. (24). Note that for $\alpha=0,\,\pi/3$, two
eigenvalues are degenerate and the corresponding nematic state is uniaxial.
(b)-(d) Fermi surfaces (42) of the 3D isotropic system in the absence (panel
(b)) or presence of nematic order in the ground state $\alpha=\pi/3$ (panel
(c), corresponding to a tensile nematic distortion) and $\alpha=0$ (panel (d),
corresponding to a compressive nematic distortion). The nematic axes
$\boldsymbol{n}$,$\boldsymbol{m}$,$\boldsymbol{l}$ are also shown; for visual
purposes, they are rescaled by the corresponding value of the Fermi momentum
along that direction,
$2.25\,k_{F}\left(\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}\right)$.

##  III Electronic nematicity in polyhedral point groups

Once full rotational symmetry is explicitly broken, the system is described in
terms of point groups, which, in contrast to the continuous orthogonal group
$\mathsf{O(3)}$, are finite groups. When investigating electronically ordered
states, it is customary to focus on the $32$ crystallographic point groups,
which have at most twofold, threefold, fourfold, or sixfold rotation symmetry.
However, quasicrystals [82, 83] have been recently shown to realize various
electronically-driven phenomena observed in periodic crystals [90, 91], such
as superconductivity [84, 85], magnetism [86], and quantum criticality [87].
While nematicity has not yet been observed in quasicrystalline environments,
it is interesting to analyze this possibility not only to establish valuable
predictions for future experiments, but also to gain deep insights about the
structure of electronic nematicity as the symmetries of the non-isotropic
system are systematically reduced. Therefore, in this paper, we will not
restrict our analysis to crystallographic point groups only, but will also
consider electronic nematicity in point groups that describe icosahedral,
dodecagonal, decagonal, and octagonal quasicrystals, as well as twisted
quasicrystals [92, 93, 94, 95, 96, 85].

There are two different classes of point groups, namely, polyhedral groups
(which do not have axial symmetry) and axial groups (which, as the name
implies, have cylindrical symmetry). As we show below, in what concerns
nematicity, the polyhedral groups are special, since none of the five elements
of the vector $\boldsymbol{d}$ transform as a trivial IR. In fact, they always
transform as multi-dimensional IRs, which can be either 2-dimensional,
3-dimensional, or 5-dimensional. This makes a description of electronic
nematicity in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation particularly
insightful and convenient. In contrast, in the axial groups, at least one of
the elements of $\boldsymbol{d}$ transforms trivially, and the other
components transform either as 2-dimensional IRs or as 1-dimensional IRs. As
we will see, this does not preclude a description in terms of the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation.

The seven existing polyhedral groups, which are the focus of this section, are
depicted in Fig. 1 and include 2 octahedral, 3 tetrahedral, and 2 icosahedral
point groups. The octahedral and tetrahedral groups define the cubic crystal
system, whereas the icosahedral groups describe a large class of quasicrystals
[89].

The two icosahedral groups $\mathsf{I}$ and
$\mathsf{I_{h}}=\mathsf{I}\times\left\\{E,I\right\\}$ have a total of $60$ and
$120$ symmetry elements. Those are constructed from twofold ($15$ axes),
threefold ($10$ axes), and fivefold ($6$ axes) rotations. For later
convenience, we explicitly list the sets of $6$ fivefold and the $10$
threefold rotation axes:

$\displaystyle\mathcal{V}_{5}^{\mathrm{ico}}$
$\displaystyle=\Bigg{\\{}\left(\begin{smallmatrix}\pm\alpha_{+}^{(5)}\\\
\alpha_{-}^{(5)}\\\ 0\end{smallmatrix}\right),\left(\begin{smallmatrix}0\\\
\pm\alpha_{+}^{(5)}\\\
\alpha_{-}^{(5)}\end{smallmatrix}\right),\left(\begin{smallmatrix}\alpha_{-}^{(5)}\\\
0\\\ \pm\alpha_{+}^{(5)}\end{smallmatrix}\right)\Bigg{\\}},$ (43)
$\displaystyle\mathcal{V}_{3}^{\mathrm{ico}}$
$\displaystyle=\Bigg{\\{}\left(\begin{smallmatrix}\pm\alpha_{-}^{(3)}\\\
\alpha_{+}^{(3)}\\\ 0\end{smallmatrix}\right),\left(\begin{smallmatrix}0\\\
\pm\alpha_{-}^{(3)}\\\
\alpha_{+}^{(3)}\end{smallmatrix}\right),\left(\begin{smallmatrix}\alpha_{+}^{(3)}\\\
0\\\ \pm\alpha_{-}^{(3)}\end{smallmatrix}\right),\mathcal{V}_{111}\Bigg{\\}}.$
(44)

Here, we defined $\alpha_{\pm}^{(3)}=\frac{1}{\sqrt{6}}\sqrt{3\pm\sqrt{5}}$,
$\alpha_{\pm}^{(5)}=\frac{1}{\sqrt{10}}\sqrt{5\pm\sqrt{5}}$, as well as the
set

$\displaystyle\mathcal{V}^{111}$
$\displaystyle=\frac{1}{\sqrt{3}}\left\\{\left(\begin{smallmatrix}1\\\ 1\\\
1\end{smallmatrix}\right),\left(\begin{smallmatrix}\text{-}1\\\ \text{-}1\\\
1\end{smallmatrix}\right),\left(\begin{smallmatrix}1\\\ \text{-}1\\\
\text{-}1\end{smallmatrix}\right),\left(\begin{smallmatrix}\text{-}1\\\ 1\\\
\text{-}1\end{smallmatrix}\right)\right\\},$ (45)

containing the four corners of a tetrahedron. The five octahedral and
tetrahedral point groups, which form the cubic crystal system, have much fewer
elements, and can be conveniently expressed as

$\displaystyle\mathsf{T}$ $\displaystyle=\\{E,C_{3\alpha}^{\pm
1}\\}\times\\{E,C_{2z}\\}\times\\{E,C_{2x}\\},$ (46)
$\displaystyle\mathsf{T_{h}}$ $\displaystyle=\mathsf{T}\times\\{E,I\\},$ (47)
$\displaystyle\mathsf{T_{d}}$ $\displaystyle=\mathsf{T}\times\\{E,IC_{4z}\\},$
(48) $\displaystyle\mathsf{O}$ $\displaystyle=\mathsf{T}\times\\{E,C_{4z}\\},$
(49) $\displaystyle\mathsf{O_{h}}$
$\displaystyle=\mathsf{T}\times\\{E,C_{4z}\\}\times\\{E,I\\}.$ (50)

Within the cubic crystal system (46)-(49), $\mathsf{O_{h}}$ is the supergroup,
since all the other ones are subgroups to $\mathsf{O_{h}}$. The characteristic
elements of the $\mathsf{T}$-group comprise $4$ threefold rotation axes
$\hat{\boldsymbol{\ell}}\in\mathcal{V}^{111}$, see Eq. (45), as well as $3$
twofold rotation axes $\hat{\boldsymbol{\ell}}\in\mathcal{V}^{100}$, where

$\displaystyle\mathcal{V}^{100}$
$\displaystyle=\left\\{\hat{\boldsymbol{e}}_{x},\hat{\boldsymbol{e}}_{y},\hat{\boldsymbol{e}}_{z}\right\\}.$
(51)

For instance, the element $C_{3\alpha}^{\pm 1}$ denotes a rotation by an angle
$\vartheta=\pm\frac{2\pi}{3}$ about the axis
$\hat{\boldsymbol{\ell}}=(1,1,1)/\sqrt{3}$. Similarly, $C_{2z}$ corresponds to
a $\vartheta=\frac{2\pi}{2}$ rotation about
$\hat{\boldsymbol{\ell}}=\hat{\boldsymbol{e}}_{z}$.

###  III.1 Icosahedral quasicrystals

Within the polyhedral point groups, the icosahedral class contains the largest
number of symmetry elements. In fact, the large number of symmetry elements
forces the five-component order parameter $\boldsymbol{d}$ to still transform
as a five-dimensional IR. In other words, in an icosahedral environment the
nematic order parameter is not symmetry-decomposed but remains the five-
dimensional expression (17), i.e.

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{icosahedral}}{\longrightarrow}\quad\boldsymbol{d}.$
(52)

This property is unique to icosahedral quasicrystals as no other
crystallographic or non-crystallographic point group allows for a five-
component IR. For concreteness, in this section we focus on the larger group
$\mathsf{I_{h}}$, but our results also apply to the group $\mathrm{I}$. Within
$\mathsf{I_{h}}$, the order parameter transforms according to the five-
dimensional IR $H_{g}$. While both the isotropic $\mathsf{O(3)}$ and the
icosahedral $\mathsf{I_{h}}$ environments host five-component nematic order
parameters, the preferred directions of the nematic director in the latter
should be severely restricted due to the finite number of symmetry axes.
Therefore, because the system lacks a continuous rotational symmetry, the
nematic ground state of an icosahedral quasicrystal should be qualitatively
different from the isotropic case. Mathematically, this difference should be
reflected in the corresponding nematic Landau expansion, which we now derive.

Here, and in all subsequent sections of this paper, we follow the same
strategy to derive the Landau expansion. First, we compute the decomposition,
for each expansion order, of the symmetrized products of the IR according to
which the nematic order parameter transforms (see Appendix B). We find:

$\displaystyle\big{[}\\!\otimes_{j=1}^{2}H_{g}\big{]}_{s}$
$\displaystyle=A_{g}\oplus G_{g}\oplus 2H_{g},$ (53)
$\displaystyle\big{[}\\!\otimes_{j=1}^{3}H_{g}\big{]}_{s}$
$\displaystyle=2A_{g}\oplus T_{1g}\oplus T_{2g}\oplus 3G_{g}\oplus 3H_{g},$
(54) $\displaystyle\big{[}\\!\otimes_{j=1}^{4}H_{g}\big{]}_{s}$
$\displaystyle=2A_{g}\\!\oplus\\!2T_{1g}\\!\oplus\\!2T_{2g}\\!\oplus
4G_{g}\\!\oplus 8H_{g}.$ (55)

Therefore, the Landau expansion up to fourth order contains a total of $5$
invariants (i.e. $5$ terms that transform as the trivial IR $A_{g}$), in
contrast to the $3$ invariants of the isotropic case, see Eqs. (26)-(28). To
construct these five invariants it is convenient to determine the bilinear
combinations associated with Eq. (53), which we denote by
$D^{A_{g}}=\left|\boldsymbol{d}\right|^{2}$, $\boldsymbol{D}^{G_{g}}$,
$\boldsymbol{D}^{H_{g},1}$ and $\boldsymbol{D}^{H_{g},2}$. Note that the two
$H_{g}$-bilinears are degenerate. As explained in Appendix C, we choose a
representation where $\boldsymbol{D}^{H_{g},1}=\boldsymbol{D}^{j=2}$ , as
defined in Eq. (29), and

$\displaystyle\boldsymbol{D}^{H_{g},2}=$
$\displaystyle\mathcal{R}_{5}\left(-\phi_{0}\right)\left(\begin{smallmatrix}\tilde{d}_{1}^{2}-\tilde{d}_{2}^{2}+\frac{1}{2}\left(\tilde{d}_{3}^{2}+\tilde{d}_{4}^{2}-2\tilde{d}_{5}^{2}\right)\\\
-2\tilde{d}_{1}\tilde{d}_{2}-\frac{\sqrt{3}}{2}\left(\tilde{d}_{3}^{2}-\tilde{d}_{4}^{2}\right)\\\
\left(\tilde{d}_{1}-\sqrt{3}\tilde{d}_{2}\right)\tilde{d}_{3}+\sqrt{3}\tilde{d}_{4}\tilde{d}_{5}\\\
\left(\tilde{d}_{1}+\sqrt{3}\tilde{d}_{2}\right)\tilde{d}_{4}+\sqrt{3}\tilde{d}_{3}\tilde{d}_{5}\\\
-2\tilde{d}_{1}\tilde{d}_{5}+\sqrt{3}\tilde{d}_{3}\tilde{d}_{4}\end{smallmatrix}\right).$
(56)

Here, we introduced
$\tilde{\boldsymbol{d}}=\mathcal{R}_{5}(\phi_{0})\boldsymbol{d}$ with angle
$\phi_{0}=\arccos\left(-1/4\right)$ and the particular rotation matrix

$\displaystyle\mathcal{R}_{5}(\phi)$
$\displaystyle=\left(\begin{smallmatrix}\cos\phi&\sin\phi&0&0&0\\\
-\sin\phi&\cos\phi&0&0&0\\\ 0&0&1&0&0\\\ 0&0&0&1&0\\\
0&0&0&0&1\end{smallmatrix}\right).$ (57)

Correspondingly, the two $H_{g}$-bilinears are related through the equation
$\boldsymbol{D}^{H_{g},2}=\mathcal{R}_{5}(-\phi_{0})\boldsymbol{D}^{H_{g},1}\big{|}_{\boldsymbol{d}\rightarrow\tilde{\boldsymbol{d}}}$,
and thus have the same magnitude
$\left|\boldsymbol{D}^{H_{g},1}\right|=\left|\boldsymbol{D}^{H_{g},2}\right|=\left|\boldsymbol{d}\right|^{2}$.
For our purposes, we do not need $\boldsymbol{D}^{G_{g}}$, although it can be
derived in a straightforward way. With the aid of the two $H_{g}$-bilinears,
the five $A_{g}$ invariants in Eqs. (53)-(55) can be constructed as
$|\boldsymbol{d}|^{2}$, $|\boldsymbol{d}|^{4}$,
$\boldsymbol{d}\cdot\boldsymbol{D}^{H_{g},1}$,
$\boldsymbol{d}\cdot\boldsymbol{D}^{H_{g},2}$, and
$\boldsymbol{D}^{H_{g},1}\cdot\boldsymbol{D}^{H_{g},2}$. The resulting nematic
Landau action becomes

$\displaystyle\mathcal{S}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{0}\left|\boldsymbol{d}\right|^{2}+\left|\boldsymbol{d}\right|^{3}\Big{(}g_{1}\hat{\boldsymbol{d}}\cdot\hat{\boldsymbol{D}}^{H_{g},1}+g_{2}\hat{\boldsymbol{d}}\cdot\hat{\boldsymbol{D}}^{H_{g},2}\Big{)}$
$\displaystyle\qquad\qquad\qquad+\left|\boldsymbol{d}\right|^{4}\Big{(}u_{1}+u_{2}\,\hat{\boldsymbol{D}}^{H_{g},1}\cdot\hat{\boldsymbol{D}}^{H_{g},2}\Big{)}\Big{\\}},$
(58)

with the unit vectors $\hat{\boldsymbol{d}}=\boldsymbol{d}/|\boldsymbol{d}|$
and
$\hat{\boldsymbol{D}}^{H_{g},1/2}=\boldsymbol{D}^{H_{g},1/2}/|\boldsymbol{d}|^{2}$,
the cubic coefficients $g_{1,2}$, and the quartic coefficients $u_{1,2}$,
which are restricted to $0\leq|u_{2}|<u_{1}$ in order for the action to be
bounded. Since the action (58) has not been explored in the literature, we
first establish the corresponding mean-field ground-state phase diagram. The
details of the derivation are given in Appendix C; here, we focus on the
numerically-obtained phase diagram shown in Fig. 3 in the
$\left(g_{2}/g_{1},\,u_{2}/u_{1}\right)$ parameter-space. Clearly, the Landau
expansion (58) has two distinct ground-state phases, which we dub
$C_{3}$-nematic and $C_{5}$-nematic, based on the residual rotational symmetry
of each state.

Figure 3: Nematic mean-field phase diagram of an icosahedral quasicrystal
obtained by a mean-field minimization of the action (58). The parameters
$u_{1},\,u_{2}$ and $g_{1},\>g_{2}$ are the quartic and cubic Landau
coefficients, respectively. Each state is characterized by the angle $\alpha$
denoting either a tensile nematic state ($\alpha=\pi/3)$ or a compressive
nematic state ($\alpha=0$), see Eq. (24). In this figure, we set $g_{1}>0$.
The case $g_{1}<0$ is recovered by switching
$(\alpha=0)\leftrightarrow(\alpha=\pi/3)$.

We first discuss the $C_{3}$-nematic phase, by focusing on the parameter
regime $\mathrm{sign}\,g_{1}=\mathrm{sign}\,g_{2}$ and $u_{2}<0$ (bottom right
quadrant in Fig. 3). In this region, the condition
$\hat{\boldsymbol{D}}^{H_{g},1}=\hat{\boldsymbol{D}}^{H_{g},2}=-\mathrm{sign}\,g_{1}\hat{\boldsymbol{d}}$
minimizes not only the two cubic terms of Eq. (58), but also the anisotropic
quartic term, i.e. all three direction-dependent terms are simultaneously
minimized. This condition is satisfied for $10$ directions of the nematic
order parameter $\hat{\boldsymbol{d}}$, which are associated with the $10$
threefold symmetry axes $\mathcal{V}_{3}^{\mathrm{ico}}$ in Eq. (44). Indeed,
in the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation (24), the
$C_{3}$-nematic ground state is given by

$\displaystyle g_{1}$ $\displaystyle>0:$
$\displaystyle\boldsymbol{d}_{3,i}^{\mathrm{ico}}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}\right|,\alpha=\pi/3,\boldsymbol{n}_{i},\boldsymbol{m}_{i}=\mathcal{V}_{3,i}^{\mathrm{ico}},\boldsymbol{l}_{i}\right],$
(59) $\displaystyle g_{1}$ $\displaystyle<0:$
$\displaystyle\boldsymbol{d}_{3,i}^{\mathrm{ico}}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}\right|,\alpha=0,\boldsymbol{n}_{i}=\mathcal{V}_{3,i}^{\mathrm{ico}},\boldsymbol{m}_{i},\boldsymbol{l}_{i}\right],$
(60)

with $i=1,\ldots,10$. Since $\alpha=\pi/3$ or $\alpha=0$, both ground states
are uniaxial, see Eqs. (36)-(37), with the state described by Eq. (59) being a
tensile nematic state ($\alpha=\pi/3$) and the state described by Eq. (60)
being a compressive nematic state ($\alpha=0$). In Fig. 4(a), we plot the
Fermi surface distortion associated with the tensile $C_{3}$-nematic ground
state (59) using Eq. (41). Note that while strictly speaking the crystal
momentum is not a good quantum number for a quasicrystal, plotting Eq. (42) is
useful to visualize, even if perturbatively only, the effect of the broken
rotational symmetry on the otherwise spherically symmetric electronic charge
distribution. To visualize the allowed directions of the nematic director
$\boldsymbol{m}$, we draw in the same figure an icosahedron concentric to the
Fermi surface, which consists of $20$ identical equilateral triangles.
Clearly, the long axis $\boldsymbol{m}$ points towards the center of one these
equilateral triangles. Since there are $20$ such triangles, and because
$\pm\boldsymbol{m}$ are identical nematic states, there are indeed $10$
degenerate states. Because $\boldsymbol{n}_{i}$ or $\boldsymbol{m}_{i}$ can
only point along 10 symmetry-related directions, the broken symmetry is
discrete. Note that the resulting point group in this case is
$\mathsf{D_{3d}}$. While it belongs to the crystallographic trigonal system,
this does not imply that the quasicrystal will become a crystal rather than an
incommensurate crystal. Nevertheless, it is interesting that a nematic
transition changes the point group from non-crystallographic to
crystallographic.

To discuss the $C_{5}$-nematic phase, we consider the parameter range
$\mathrm{sign}\,g_{1}=-\mathrm{sign}\,g_{2}$ and $u_{2}>0$, corresponding to
the top left quadrant in the phase diagram of Fig. 3. In this case, the
simultaneous minimization of the three anisotropic terms of the action (58) is
achieved by the condition
$\hat{\boldsymbol{D}}^{H_{g},1}=-\hat{\boldsymbol{D}}^{H_{g},2}=-\mathrm{sign}\,g_{1}\hat{\boldsymbol{d}}$.
This condition, in turn, is satisfied for six directions of
$\hat{\boldsymbol{d}}$ that are associated with the $6$ fivefold symmetry axes
$\mathcal{V}_{5}^{\mathrm{ico}}$ in Eq. (43). In analogy to the
$C_{3}$-nematic phase, the $C_{5}$-nematic ground state is uniaxial and
parametrized as:

$\displaystyle g_{1}$ $\displaystyle>0:$
$\displaystyle\boldsymbol{d}_{5,i}^{\mathrm{ico}}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}\right|,\alpha=\pi/3,\boldsymbol{n}_{i},\boldsymbol{m}_{i}=\mathcal{V}_{5,i}^{\mathrm{ico}},\boldsymbol{l}_{i}\right],$
(61) $\displaystyle g_{1}$ $\displaystyle<0:$
$\displaystyle\boldsymbol{d}_{5,i}^{\mathrm{ico}}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}\right|,\alpha=0,\boldsymbol{n}_{i}=\mathcal{V}_{5,i}^{\mathrm{ico}},\boldsymbol{m}_{i},\boldsymbol{l}_{i}\right],$
(62)

where $i=1,\ldots,6$. We visualize the tensile Fermi surface distortion
associated with the tensile state (61) in Fig. 4(b). In the $C_{5}$-nematic
phase, the nematic director $\boldsymbol{m}$ aligns with the corners of the
icosahedron. Since there are $12$ such corners, we find indeed $6$ degenerate
$C_{5}$-nematic states. The residual symmetries inside the nematic phase
result in the point group $\mathsf{D_{5d}}$, which is also a non-
crystallographic point group.

The nematic phase diagram in Fig. 3 shows that the $C_{3}$-nematic state (59)
and the $C_{5}$-nematic state (61) actually occupy the entire phase space,
including regions where not all three anisotropic terms in (58) can be
satisfied simultaneously. As shown in the top-right and bottom-left corner, a
tensile ground state can become a compressive one upon traversing the phase
diagram.

Figure 4: Distorted Fermi surface in the nematic phase, as given by Eq. (42),
plotted together with the nematic axes
$\boldsymbol{n}$,$\boldsymbol{m}$,$\boldsymbol{l}$. (a) Icosahedral
$C_{3}$-nematic state (59); (b) icosahedral $C_{5}$-nematic state (61); (c)
cubic triplet ($T_{2g}$) state (78); (d) cubic doublet ($E_{g}$) state (88).
In all cases, the $\alpha=\pi/3$ (tensile) state is plotted, and the nematic
director is parallel to $\boldsymbol{m}$. The unit vectors
$\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}$ are rescaled
by
$2.2\,k_{F}\left(\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}\right)$,
respectively.

###  III.2 Cubic crystals

The cubic crystal system consists of two octahedral point groups
($\mathsf{O_{h}}$, $\mathsf{O}$) and three tetrahedral point groups
($\mathsf{T_{h}}$, $\mathsf{T_{d}}$, $\mathsf{T}$). In all five cases, the
reduced symmetry with respect to the icosahedral point group leads to a
symmetry decomposition of the five-component nematic order parameter
$\boldsymbol{d}$ (17) according to:

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{cubic}}{\longrightarrow}\quad\left(\boldsymbol{d}^{e},\,\boldsymbol{d}^{t}\right)^{T},$
(63)

where the doublet $\boldsymbol{d}^{e}$ and the triplet $\boldsymbol{d}^{t}$
nematic vectors are given by:

$\displaystyle\boldsymbol{d}^{e}$
$\displaystyle=\left(\begin{array}[]{c}d_{1}\\\
d_{2}\end{array}\right)=\left(\begin{array}[]{c}d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}\\\
d_{x^{2}-y^{2}}\end{array}\right)=\left|\boldsymbol{d}^{e}\right|\left(\begin{array}[]{c}\cos\gamma_{e}\\\
\sin\gamma_{e}\end{array}\right),$ (70) $\displaystyle\boldsymbol{d}^{t}$
$\displaystyle=\left(d_{3},d_{4},d_{5}\right)^{T}=\left(d_{2yz},d_{2xz},d_{2xy}\right)^{T}.$
(71)

In terms of the IRs of the five groups ($\mathsf{O_{h}}$, $\mathsf{O}$,
$\mathsf{T_{h}}$, $\mathsf{T_{d}}$, $\mathsf{T}$), $\boldsymbol{d}^{t}$ always
transforms as a three-dimensional IR ($T_{2g}$, $T_{2}$, $T_{g}$, $T_{2}$, and
$T$, respectively) whereas $\boldsymbol{d}^{e}$ transforms as a two-
dimensional IR ($E_{g}$, $E$, $E_{g}\oplus\bar{E}_{g}$, $E$, and
$E\oplus\bar{E}$, respectively). Our goal is to derive the Landau expansions
for both types of nematics and then determine the mean-field ground states in
the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation (24). For
concreteness, we focus on the most symmetric cubic point group,
$\mathsf{O_{h}}$, and comment on the extension to the other four cubic point
groups in the end of this section.

We begin with the nematic triplet order parameter $\boldsymbol{d}^{t}$ (71),
which transforms as the $T_{2g}$ IR of $\mathsf{O_{h}}$. The decomposition of
the symmetrized product for each expansion order gives

$\displaystyle\big{[}\otimes_{j=1}^{2}T_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{g}\oplus T_{2g},$ (72)
$\displaystyle\big{[}\otimes_{j=1}^{3}T_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus T_{1g}\oplus 2T_{2g},$ (73)
$\displaystyle\big{[}\otimes_{j=1}^{4}T_{2g}\big{]}_{s}$
$\displaystyle=2A_{1g}\oplus 2E_{g}\oplus T_{1g}\oplus 2T_{2g},$ (74)

which implies the existence of $4$ invariants in the Landau expansion up to
fourth order. To systematically construct them, it is convenient to exploit
the bilinears associated with the quadratic decomposition, Eq. (72). We find
that $D^{A_{1g}}=|\boldsymbol{d}^{t}|^{2}$ and

$\displaystyle\boldsymbol{D}^{E_{g}}$
$\displaystyle=\left(\begin{smallmatrix}\frac{1}{\sqrt{3}}\left(2d_{5}^{2}-d_{3}^{2}-d_{4}^{2}\right)\\\
\left(d_{3}^{2}-d_{4}^{2}\right)\end{smallmatrix}\right),$
$\displaystyle\boldsymbol{D}^{T_{2g}}$
$\displaystyle=\left(\begin{smallmatrix}2d_{4}d_{5}\\\ 2d_{3}d_{5}\\\
2d_{3}d_{4}\end{smallmatrix}\right).$ (75)

In terms of these bilinears, the four invariants become
$\left|\boldsymbol{d}^{t}\right|^{2}$ to quadratic order,
$\boldsymbol{d}^{t}\cdot\boldsymbol{D}^{T_{2g}}$ to cubic order, and
$\left|\boldsymbol{d}^{t}\right|^{4}$ and
$\boldsymbol{D}^{T_{2g}}\cdot\boldsymbol{D}^{T_{2g}}$ to quartic order. Note
that the quartic term could equally be written in terms of
$\boldsymbol{D}^{E_{g}}\cdot\boldsymbol{D}^{E_{g}}$ by using the Fierz
identity
$\boldsymbol{D}^{E_{g}}\cdot\boldsymbol{D}^{E_{g}}=\frac{4}{3}\left|\boldsymbol{d}^{t}\right|^{4}-\boldsymbol{D}^{T_{2g}}\cdot\boldsymbol{D}^{T_{2g}}$.
Therefore, the resulting Landau expansion is

$\displaystyle\mathcal{S}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{0}\left|\boldsymbol{d}^{t}\right|^{2}+g\,d_{3}d_{4}d_{5}+u_{1}\left|\boldsymbol{d}^{t}\right|^{4}$
$\displaystyle+u_{2}\left(d_{3}^{2}d_{4}^{2}+d_{3}^{2}d_{5}^{2}+d_{4}^{2}d_{5}^{2}\right)\Big{\\}},$
(76)

with the cubic and quartic Landau coefficients $g$, $u_{1}$, $u_{2}$. The
Landau expansion (76) is the the same as that of the 4-state Potts model
($Z_{4}$-Potts) [103]. It is well established that the the upper critical
dimension of this model is below $3$, $d_{u}<3$ [104]. Consequently, a mean-
field solution of Eq. (76) is appropriate, resulting in a first-order
transition and in a fourfold degenerate ground-state manifold given by
$\boldsymbol{d}_{0,i}^{t}=\mathrm{sign}\left(-g\right)\left|\boldsymbol{d}^{t}\right|\mathcal{V}_{i}^{Z_{4}}$
with $i=1,\ldots,4$ and:

$\displaystyle\mathcal{V}^{Z_{4}}$
$\displaystyle=\frac{1}{\sqrt{3}}\left\\{\left(\begin{smallmatrix}1\\\ 1\\\
1\end{smallmatrix}\right),\left(\begin{smallmatrix}\text{-}1\\\ \text{-}1\\\
1\end{smallmatrix}\right),\left(\begin{smallmatrix}1\\\ \text{-}1\\\
\text{-}1\end{smallmatrix}\right),\left(\begin{smallmatrix}\text{-}1\\\ 1\\\
\text{-}1\end{smallmatrix}\right)\right\\}.$ (77)

This fourfold degeneracy is associated with the $4$ threefold symmetry axes
$\mathcal{V}^{111}$ defined in Eq. (45). Indeed, in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation of Eq. (24),
these four states are given by

$\displaystyle g$ $\displaystyle>0:$ $\displaystyle\boldsymbol{d}_{0,i}^{t}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}^{t}\right|,\alpha=\pi/3,\boldsymbol{n}_{i},\boldsymbol{m}_{i}=\mathcal{V}_{i}^{111},\boldsymbol{l}_{i}\right],$
(78) $\displaystyle g$ $\displaystyle<0:$
$\displaystyle\boldsymbol{d}_{0,i}^{t}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}^{t}\right|,\alpha=0,\boldsymbol{n}_{i}=\mathcal{V}_{i}^{111},\boldsymbol{m}_{i},\boldsymbol{l}_{i}\right].$
(79)

Since $\alpha=\pi/3$ or $\alpha=0$, the nematic states (78)-(79) are uniaxial
[recall Eqs. (36)-(37)], and the nematic director is determined by either
$\boldsymbol{m}$ or $\boldsymbol{n}$, respectively. In Fig. 4(c), we show the
distorted Fermi surface associated with the tensile (i.e. $\alpha=\pi/3$)
nematic state, Eq. (78). Note that the nematic director $\boldsymbol{m}$
points towards the corners of the cube that is concentric to the Fermi
surface. Since there are $8$ such corners, the number of degenerate states is
$4$, since $\pm\boldsymbol{m}$ correspond to the same state.

Next, we consider the nematic doublet order parameter $\boldsymbol{d}^{e}$
(70), which transforms as the $E_{g}$ IR of $\mathsf{O_{h}}$. Following the
same procedure as above, we first perform the decomposition of the symmetrized
product for each Landau-expansion order:

$\displaystyle\big{[}\\!\otimes_{j=1}^{2}E_{g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{g},$ (80)
$\displaystyle\big{[}\\!\otimes_{j=1}^{3}E_{g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus A_{2g}\oplus E_{g},$ (81)
$\displaystyle\big{[}\\!\otimes_{j=1}^{4}E_{g}\big{]}_{s}$
$\displaystyle=A_{1g}\\!\oplus\\!2E_{g},$ (82)

which reveals the existence of one $A_{1g}$ invariant per expansion order. The
bilinear combinations in Eq. (80) can be readily obtained:

$\displaystyle D^{A_{1g}}$
$\displaystyle=\left|\boldsymbol{d}^{e}\right|^{2},$
$\displaystyle\boldsymbol{D}^{E_{g}}$
$\displaystyle=\left|\boldsymbol{d}^{e}\right|^{2}\left(\begin{array}[]{c}\cos
2\gamma_{e}\\\ -\sin 2\gamma_{e}\end{array}\right).$ (85)

Using these bilinears, we can express the three invariants as
$\left|\boldsymbol{d}^{e}\right|^{2}$, $\left|\boldsymbol{d}^{e}\right|^{4}$
and $\boldsymbol{d}^{e}\cdot\boldsymbol{D}^{E_{g}}$, such that the nematic
Landau expansion becomes:

$\displaystyle\mathcal{S}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{0}\left|\boldsymbol{d}^{e}\right|^{2}+g\left|\boldsymbol{d}^{e}\right|^{3}\,\cos\left(3\gamma_{e}\right)+u\left|\boldsymbol{d}^{e}\right|^{4}\Big{\\}}.$
(86)

The nematic action (86) has the same Landau expansion as the $3$-state Potts
model ($Z_{3}$-Potts) [104]. We note that the $3$-state Potts model is in the
same universality class as the $3$-state clock model ($Z_{3}$-clock), so we
will use these terms interchangeably in the remainder of the paper. We will
discuss further the differences between $Z_{q}$-Potts and $Z_{q}$-clock models
in Sec. IV. The key point is that the upper critical dimension of the
$Z_{3}$-Potts/clock model is $d_{u}<3$. Therefore, as in the case of the
triplet nematic, a mean-field solution of Eq. (86) is appropriate. We find a
threefold degenerate ground-state manifold given by
$\boldsymbol{d}_{0,i}^{e}=\mathrm{sign}\left(-g\right)\left|\boldsymbol{d}^{e}\right|\mathcal{V}_{i}^{Z_{3}}$
with $i=1,2,3$ and

$\displaystyle\mathcal{V}^{Z_{3}}$
$\displaystyle=\left\\{\left(\begin{smallmatrix}1\\\
0\end{smallmatrix}\right),\frac{1}{2}\left(\begin{smallmatrix}-1\\\
\sqrt{3}\end{smallmatrix}\right),\frac{1}{2}\left(\begin{smallmatrix}-1\\\
-\sqrt{3}\end{smallmatrix}\right)\right\\}.$ (87)

This threefold degeneracy is related to the $3$ fourfold rotation axes
$\mathcal{V}^{100}$ in Eq. (51), which are parallel to the three coordinate
axes. This can be directly seen upon writing the nematic ground state in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation (24):

$\displaystyle g$ $\displaystyle>0:$ $\displaystyle\boldsymbol{d}_{0,i}^{e}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}^{e}\right|,\alpha=\pi/3,\boldsymbol{n}_{i},\boldsymbol{m}_{i}=\mathcal{V}_{i}^{100},\boldsymbol{l}_{i}\right],$
(88) $\displaystyle g$ $\displaystyle<0:$
$\displaystyle\boldsymbol{d}_{0,i}^{e}$
$\displaystyle=\boldsymbol{d}\left[\left|\boldsymbol{d}^{e}\right|,\alpha=0,\boldsymbol{n}_{i}=\mathcal{V}_{i}^{100},\boldsymbol{m}_{i},\boldsymbol{l}_{i}\right].$
(89)

As in the case of the $T_{2g}$ (triplet) nematic, the $E_{g}$ (doublet)
nematic state is uniaxial with the nematic director $\boldsymbol{m}$ or
$\boldsymbol{n}$ pointing along one of the three coordinate axes. This is
illustrated in Fig. 4(d), where we plot the Fermi-surface distortion
associated with the tensile (i.e. $\alpha=\pi/3$) state given by Eq. (88).

While the analysis performed above focused on the $\mathsf{O_{h}}$ group, the
same general results hold for the doublet and triplet nematics of the other
four cubic point groups, as shown in Table 1. It is important to note,
however, that the case of the nematic doublet order parameter
$\boldsymbol{d}^{e}$ in the tetrahedral point groups $\mathsf{T_{h}}$ and
$\mathsf{T}$ needs to be treated slightly differently, because
$\boldsymbol{d}^{e}$ transforms according to a _complex_ IR, which allows for
additional invariants in the Landau expansion, which in turn make the nematic
state biaxial rather than uniaxial. A detailed derivation of this case is
presented in Appendix D.

##  IV Electronic nematicity in crystallographic axial point groups

In contrast to the family of polyhedral point groups, which encompasses the
seven point groups discussed in Sec. III, there is an infinite number of axial
groups. Their defining property, as the name indicates, is their underlying
cylindrical symmetry, which implies that the system is invariant under some
$n$-fold proper or improper rotation with respect to the $z$-axis. This
symmetry can be accompanied by additional symmetries related to inversion or
twofold rotations with respect to in-plane axes, resulting in a total of seven
series of axial groups, as shown in Fig. 1. They can be constructed from two
axial group series corresponding to the cyclic group $\mathsf{C_{n}}$ and the
group $\mathsf{S_{2n}}$:

$\displaystyle\mathsf{C_{n}}$
$\displaystyle=\left\\{E,C_{nz},C_{nz}^{2},\dots,C_{nz}^{n-1}\right\\},$
$\displaystyle n$ $\displaystyle\in[2,3,\dots,\infty],$ (90)
$\displaystyle\mathsf{S_{2n}}$
$\displaystyle=\left\\{E,S_{2nz},S_{2nz}^{2},\dots,S_{2nz}^{2n-1}\right\\},$
$\displaystyle n$ $\displaystyle\in[2,3,\dots,\infty].$ (91)

Here, $E$ is the identity operation, $C_{nz}$ corresponds to an $n$-fold
rotation with respect to the $z$ axis, and the improper rotation (or
rotoinversion) operation $S_{nz}$ is defined as a proper rotation $C_{nz}$
followed by a reflection with respect to the horizontal mirror,
$S_{nz}=IC_{2z}C_{nz}$. Note that we can also extend these definitions to
include the groups $\mathsf{C_{1}}=\left\\{E\right\\}$ and
$\mathsf{S_{2}}=\mathsf{C_{i}}=\left\\{E,I\right\\}$. The remaining five
infinite series of axial groups are obtained from $\mathsf{C_{n}}$ or
$\mathsf{S_{2n}}$ in the following way:

$\displaystyle\mathsf{D_{n}}$
$\displaystyle=\mathsf{C_{n}}\times\left\\{E,C_{2x}\right\\},$ (92)
$\displaystyle\mathsf{C_{nv}}$
$\displaystyle=\mathsf{C_{n}}\times\left\\{E,IC_{2x}\right\\},$ (93)
$\displaystyle\mathsf{C_{nh}}$
$\displaystyle=\mathsf{C_{n}}\times\left\\{E,IC_{2z}\right\\},$ (94)
$\displaystyle\mathsf{D_{nd}}$
$\displaystyle=\mathsf{S_{2n}}\times\left\\{E,C_{2x}\right\\},$ (95)
$\displaystyle\mathsf{D_{nh}}$
$\displaystyle=\mathsf{C_{n}}\times\left\\{E,IC_{2z}\right\\}\times\left\\{E,C_{2x}\right\\},$
(96)

where $I$ denotes inversion and $C_{2x}$, twofold rotation with respect to the
in-plane $x$-axis. Note that the compositions $IC_{2z}$ and $IC_{2x}$
correspond to reflections with respect to the horizontal and vertical mirrors
(denoted $\sigma_{h}$ and $\sigma_{v}$), whereas $IC_{2y}C_{4z}=S_{4z}C_{2x}$
corresponds to a reflection with respect to the diagonal mirror (denoted
$\sigma_{d}$).

Upon imposing the crystallographic restriction theorem, one finds $27$
crystallographic axial point groups corresponding to the following six crystal
systems: hexagonal, trigonal, tetragonal, orthorhombic, monoclinic, and
triclinic. Together with the $5$ crystallographic polyhedral point groups
discussed in Sec. III that form the cubic crystal system, one finds $32$
crystallographic point groups. While this section focuses on the $27$
crystallographic axial groups, we will discuss non-crystallographic axial
point groups that are relevant to quasicrystalline materials [105, 88, 86, 89,
91] (octagonal, decagonal, and dodecagonal) and to twisted quasicrystals [93,
94, 95, 106] in Sec. V.

###  IV.1 2D isotropic systems

To set the stage for the analysis of electronic nematicity in axial groups, it
is instructive to analyze the case in which the system has full in-plane
rotational symmetry. To describe such a system, one possibility would be to
consider the continuous $\mathsf{SO(2)}$ group, which describes 2D rotations
and is equivalent to $\mathsf{C_{\infty}}$ defined in Eq. (90). However, for
our purposes, it is more convenient to consider the continuous dihedral group
$\mathsf{D_{\infty h}}$, since all other axial groups are subgroups of
$\mathsf{D_{\infty h}}$. Formally, the elements of this group include
continuous proper and improper 2D rotations as well as a horizontal mirror,
i.e. $\mathsf{D_{\infty h}=SO(2)}\times\\{E,I\\}\times\\{E,C_{2x}\\}$.
Hereafter, we will simply refer to such a system as a 2D isotropic system.

In the 2D isotropic case, the five-dimensional nematic order parameter
$\boldsymbol{d}$ (17) decomposes into three separate channels,

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{in\text{-}plane}\,\,\mathrm{isotropic}}{\longrightarrow}\quad\Big{(}\,\underline{d_{1}},\,\boldsymbol{d}^{\mathrm{ip}},\,\boldsymbol{d}^{\mathrm{op}}\,\Big{)}^{T}.$
(97)

The first component $d_{1}=d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}$
transforms as the trivial IR $A_{1g}$ of $\mathsf{D_{\infty h}}$, which we
indicate by the underline in the expression above. The other components
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ are two-
component “vectors” (i.e. doublets) that transform as different two-
dimensional IRs of $\mathsf{D_{\infty h}}$, $E_{2g}$ and $E_{1g}$,
respectively. They are given by:

$\displaystyle\boldsymbol{d}^{\mathrm{ip}}$
$\displaystyle=\left(\begin{array}[]{c}d_{2}\\\
d_{5}\end{array}\right)=\left(\begin{array}[]{c}d_{x^{2}-y^{2}}\\\
d_{2xy}\end{array}\right)=\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left(\begin{array}[]{c}\cos\gamma_{\mathrm{ip}}\\\
\sin\gamma_{\mathrm{ip}}\end{array}\right),$ (104)
$\displaystyle\boldsymbol{d}^{\mathrm{op}}$
$\displaystyle=\left(\begin{array}[]{c}d_{3}\\\
d_{4}\end{array}\right)=\left(\begin{array}[]{c}d_{2yz}\\\
d_{2xz}\end{array}\right)=\left|\boldsymbol{d}^{\mathrm{op}}\right|\left(\begin{array}[]{c}\cos\gamma_{\mathrm{op}}\\\
\sin\gamma_{\mathrm{op}}\end{array}\right),$ (111)

where we introduced a polar parameterization in terms of the angles
$\gamma_{\mathrm{ip}}$ and $\gamma_{\mathrm{op}}$. The superscripts
$\mathrm{ip}$ and $\mathrm{op}$ are used to indicate that the condensation of
the order parameter promotes an in-plane or out-of-plane distortion of the
Fermi surface, respectively. Importantly, as explained in Sec. II, we are
using the transformation rule in Eq. (20) to define the nematic order
parameter for all point groups, which results in the definitions (104)-(111).
In the literature, one often uses a point-group-specific basis to define the
nematic order parameter, resulting in definitions such as
$\tilde{\boldsymbol{d}}^{\mathrm{ip}}=\big{(}d_{x^{2}-y^{2}},\,-d_{2xy}\big{)}^{T}$
and $\tilde{\boldsymbol{d}}^{\mathrm{op}}=\left(d_{2yz},\,-d_{2xz}\right)^{T}$
[107]. Further implications arising from different definitions are discussed
in Sec. VI. Of course, the final results do not depend on the basis used to
express the nematic order parameter.

We emphasize that, because the nematic component $d_{1}$ transforms trivially
in $\mathsf{D_{\infty h}}$, it will do so for all axial groups. Therefore,
$d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}$ should not be interpreted as a
nematic order parameter, since it is generically non-zero for any temperature
or tuning parameter. Physically, from Eq. (42), $d_{1}\neq 0$ corresponds to a
distortion of the isotropic Fermi surface along the $k_{z}$-axis, lowering its
symmetry from spherical to cylindrical – as expected for an axial group. As we
will argue below, while $d_{1}$ does not impact the critical properties of the
nematic instability, it needs to be included if one chooses to express the
nematic order parameter $\boldsymbol{d}$ in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation of Eq. (24).

We now proceed to derive the Landau expansions for
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$. Since they
transform according to the IRs $E_{2g}$ and $E_{1g}$ of the infinite dihedral
group $\mathsf{D_{\infty h}}$, respectively, the number of invariants for each
expansion order can be read off from the symmetrized-product decomposition:

$\displaystyle\big{[}\otimes_{j=1}^{2N}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus_{n=1}^{N}E_{(4n)g},$ (112)
$\displaystyle\big{[}\otimes_{j=1}^{2N-1}E_{2g}\big{]}_{s}$
$\displaystyle=\oplus_{n=1}^{N}E_{(4n-2)g},$ (113)
$\displaystyle\big{[}\otimes_{j=1}^{2N}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus_{n=1}^{N}E_{(2n)g},$ (114)
$\displaystyle\big{[}\otimes_{j=1}^{2N-1}E_{1g}\big{]}_{s}$
$\displaystyle=\oplus_{n=1}^{N}E_{(2n-1)g},$ (115)

with $N=1,2,3$ and where we used a tensor summation notation, e.g.
$\oplus_{n=1}^{2}E_{(4n)g}=E_{4g}\oplus E_{8g}$. The decompositions
(112)-(115) are carried out up to sixth order in order to emphasize that no
$A_{1g}$ invariants occur other than
$\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2N}$ and
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2N}$, i.e. the invariants only
occur at even order. Consequently, the Landau expansions for the in-plane and
out-of-plane nematic order parameters become

$\displaystyle\mathcal{S}_{\mathrm{ip}}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{0}|\boldsymbol{d}^{\mathrm{ip}}|^{2}+u|\boldsymbol{d}^{\mathrm{ip}}|^{4}\Big{\\}},$
$\displaystyle\mathcal{S}_{\mathrm{op}}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{0}|\boldsymbol{d}^{\mathrm{op}}|^{2}+u|\boldsymbol{d}^{\mathrm{op}}|^{4}\Big{\\}}.$
(116)

Therefore, since $\boldsymbol{d}^{\mathrm{ip}}$ and
$\boldsymbol{d}^{\mathrm{op}}$ are two-component order parameters, their
actions have the same Landau expansion as the well-known XY model. Their
condensation leads to a continuous symmetry breaking, since the nematic angles
$\gamma_{\mathrm{ip}},\gamma_{\mathrm{op}}\in[0,2\pi]$ are not constrained by
symmetry and can point in any direction.

It is illustrative to rewrite the fully-isotropic nematic order parameters
(104)-(111) in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation of Eq. (24). For
the description of the in-plane nematic order parameter
$\boldsymbol{d}^{\mathrm{ip}}$, it is convenient to introduce the three
orthonormal vectors in cylindrical coordinates

$\displaystyle\boldsymbol{e}_{\parallel}^{A}$
$\displaystyle=\left(\begin{smallmatrix}\cos\left(\gamma_{\mathrm{ip}}/2\right)\\\
\sin\left(\gamma_{\mathrm{ip}}/2\right)\\\ 0\end{smallmatrix}\right),$
$\displaystyle\boldsymbol{e}_{\parallel}^{B}$
$\displaystyle=\left(\begin{smallmatrix}-\sin\left(\gamma_{\mathrm{ip}}/2\right)\\\
\cos\left(\gamma_{\mathrm{ip}}/2\right)\\\ 0\end{smallmatrix}\right),$
$\displaystyle\boldsymbol{e}_{z}$ $\displaystyle=\left(\begin{smallmatrix}0\\\
0\\\ 1\end{smallmatrix}\right).$ (117)

The magnitude of the nematic order parameter is given by
$\left|\boldsymbol{d}\right|=\sqrt{\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}+\left(d_{1}\right)^{2}}$,
where $d_{1}$, as explained above, is always non-zero, and can thus be
considered as an intrinsic parameter characterizing the 2D isotropic system
even in the absence of symmetry-breaking nematic order. Depending on the value
of the ratio $-1<d_{1}/\left|\boldsymbol{d}\right|<1$, we find three different
regimes of parameters that characterize $\boldsymbol{d}$ in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation:

$\displaystyle\frac{\left|d_{1}\right|}{\left|\boldsymbol{d}\right|}$
$\displaystyle\\!\leq\\!\frac{1}{2}\,:$ $\displaystyle\\!\\!\alpha$
$\displaystyle=\frac{\pi}{6}+\arcsin\left(\frac{d_{1}}{\left|\boldsymbol{d}\right|}\right),$
$\displaystyle\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{\parallel}^{A},\;\boldsymbol{m}=\boldsymbol{e}_{\parallel}^{B},\;\boldsymbol{l}=\boldsymbol{e}_{z},$
$\displaystyle\frac{d_{1}}{\left|\boldsymbol{d}\right|}$
$\displaystyle\\!<\\!\frac{\text{-}1}{2}\\!:$ $\displaystyle\\!\\!\alpha$
$\displaystyle=\arcsin\left(\frac{\left|d_{1}\right|}{\left|\boldsymbol{d}\right|}\right)-\frac{\pi}{6},$
$\displaystyle\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{\parallel}^{A},\;\boldsymbol{m}=\boldsymbol{e}_{z},\;\boldsymbol{l}=\boldsymbol{e}_{\parallel}^{B},$
$\displaystyle\frac{d_{1}}{\left|\boldsymbol{d}\right|}$
$\displaystyle\\!>\\!\frac{1}{2}\,:$ $\displaystyle\\!\\!\alpha$
$\displaystyle=\frac{\pi}{2}-\arcsin\left(\frac{\left|d_{1}\right|}{\left|\boldsymbol{d}\right|}\right),$
$\displaystyle\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{z},\;\boldsymbol{m}=\boldsymbol{e}_{\parallel}^{B},\;\boldsymbol{l}=\boldsymbol{e}_{\parallel}^{A}.$
(118)

As for the out-of-plane nematic order parameter, we have
$\left|\boldsymbol{d}\right|=\sqrt{\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}+\left(d_{1}\right)^{2}}$
and:

$\displaystyle\alpha$
$\displaystyle=\frac{\pi}{6}-\arcsin\left[d_{1}/\left(2\left|\boldsymbol{d}\right|\right)\right],$
$\displaystyle\boldsymbol{n}$
$\displaystyle=\left(\begin{smallmatrix}\sin\gamma_{\mathrm{op}}\cos\eta_{\mathrm{op}}\\\
\cos\gamma_{\mathrm{op}}\cos\eta_{\mathrm{op}}\\\
\sin\eta_{\mathrm{op}}\end{smallmatrix}\right),\;\;\boldsymbol{m}=\left(\begin{smallmatrix}\sin\gamma_{\mathrm{op}}\sin\eta_{\mathrm{op}}\\\
\cos\gamma_{\mathrm{op}}\sin\eta_{\mathrm{op}}\\\
-\cos\eta_{\mathrm{op}}\end{smallmatrix}\right),\;\;\boldsymbol{l}=\left(\begin{smallmatrix}\cos\gamma_{\mathrm{op}}\\\
-\sin\gamma_{\mathrm{op}}\\\ 0\end{smallmatrix}\right),$ (119)

where we defined the polar angle:

$\displaystyle\eta_{\mathrm{op}}$
$\displaystyle=-\frac{1}{2}\mathrm{sign}\left(d_{1}\right)\arcsin\Bigg{(}\left|\boldsymbol{d}^{\mathrm{op}}\right|\Big{/}\sqrt{\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}+\frac{3}{4}\left(d_{1}\right)^{2}}\;\Bigg{)}$
$\displaystyle+\frac{\pi}{2}\frac{1+\mathrm{sign}\left(d_{1}\right)}{2}.$
(120)

Both the in-plane (118) and the out-of-plane (119) nematic states are
generically biaxial, since $\alpha\neq 0,\,\pi/3$. Note that, if the
contribution from the trivial nematic component is negligible,
$\left|d_{1}\right|\ll\left|\boldsymbol{d}^{\mathrm{ip}}\right|,\left|\boldsymbol{d}^{\mathrm{op}}\right|$,
both states are “maximally” biaxial with $\alpha=\pi/6$ [cf. Fig. 2(a)]. In
this regard, a non-zero $d_{1}$ brings the states (118) and (119) closer to
the uniaxial regime. Indeed, in the limit
$\left|d_{1}\right|\gg\left|\boldsymbol{d}^{\mathrm{ip}}\right|,\left|\boldsymbol{d}^{\mathrm{op}}\right|$,
$\alpha\rightarrow 0,\,\pi/3$ in both cases, with the nematic director
$\boldsymbol{n}$ or $\boldsymbol{m}$ pointing along the $k_{z}$-axis.
Interestingly, the nematic in-plane order parameter (118) also establishes a
uniaxial state when
$\left|d_{1}\right|=\left|\boldsymbol{d}^{\mathrm{ip}}\right|/\sqrt{3}$,
corresponding to $\left|d_{1}\right|\big{/}\left|\boldsymbol{d}\right|=1/2$.
In this fine-tuned parameter regime, there is a swap between the nematic axes
that point along the $k_{z}$-axis, since for
$\left|d_{1}\right|>\left|\boldsymbol{d}^{\mathrm{ip}}\right|/\sqrt{3}$ either
the long or the short nematic axis ($\boldsymbol{m}$ or $\boldsymbol{n}$)
points out-of-plane, while for
$\left|d_{1}\right|<\left|\boldsymbol{d}^{\mathrm{ip}}\right|/\sqrt{3}$ the
long and short nematic axes both lie within the plane. Since
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ also emerge
in the description of hexagonal and tetragonal lattices, we will defer
plotting the Fermi surface distortion patterns triggered by the condensation
of these nematic order parameters to the subsequent subsections.

An alternative representation of the 2D isotropic in-plane nematic order
parameter $\boldsymbol{d}^{\mathrm{ip}}$ has also been widely employed in the
literature [37, 48, 61, 49, 50, 51, 62, 63], building on the analogy with the
classical _two-dimensional_ tensorial nematic order parameter
$Q_{\mu\mu^{\prime}}=\left\langle
a_{\mu}a_{\mu^{\prime}}-\frac{1}{2}\delta_{\mu\mu^{\prime}}\boldsymbol{a}^{2}\right\rangle$,
with director $\boldsymbol{a}=\left(a_{x},\,a_{y}\right)$. In terms of the
quadrupolar order parameters, the electronic nematic order parameter is given
by:

$Q=\left(\begin{array}[]{cc}d_{x^{2}-y^{2}}&d_{2xy}\\\
d_{2xy}&-d_{x^{2}-y^{2}}\end{array}\right)=\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{\tau}^{\mathrm{ip}},$
(121)

where $\boldsymbol{\tau}^{\mathrm{ip}}=(\tau^{3},\tau^{1})$ with Pauli
matrices $\tau^{i}$. This form establishes a straightforward connection
between the two-component “vector”
$\boldsymbol{d}^{\mathrm{ip}}=\left(d_{2},\,d_{5}\right)^{T}=\big{(}d_{x^{2}-y^{2}},\,d_{2xy}\big{)}^{T}$
and $Q$ similar to the 3D case in Eq. (18), with the symmetric Pauli matrices
replacing the symmetric Gell-Mann matrices. To preserve rotational invariance,
the Landau expansion can only depend on traces of powers of $Q$:

$\mathcal{S}_{\mathrm{ip}}=\int_{\mathsf{x}}\sum_{n}\frac{\kappa_{n}}{2}\,\mathrm{tr}\left(Q^{n}\right),$
(122)

where $\kappa_{n}$ are Landau coefficients. Using the facts that
$\left(\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{\tau}^{\mathrm{ip}}\right)^{2n}=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2n}\tau^{0}$,
and thus
$\mathrm{tr}\left[(\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{\tau}^{\mathrm{ip}})^{2n+1}\right]=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2n}\mathrm{tr}\left[\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{\tau}^{\mathrm{ip}}\right]=0$,
we find:

$\displaystyle\mathcal{S}_{\mathrm{ip}}$
$\displaystyle=\int_{\mathsf{x}}\sum_{n}\kappa_{2n}\,|\boldsymbol{d}^{\mathrm{ip}}|^{2n}=\int_{\mathsf{x}}\Big{\\{}r_{0}|\boldsymbol{d}^{\mathrm{ip}}|^{2}+u|\boldsymbol{d}^{\mathrm{ip}}|^{4}+\mathcal{O}\left(|\boldsymbol{d}^{\mathrm{ip}}|^{6}\right)\Big{\\}},$
(123)

which is identical to Eq. (116). When comparing to the Landau expansion of the
3D isotropic case, Eq. (31), the main difference is the absence of the cubic
term. Ultimately, this can be traced back to the fundamental differences in
the $\mathrm{SU(2)}$ and $\mathrm{SU(3)}$ Lie algebras, reflected in the fact
that different Pauli matrices necessarily anti-commute whereas different Gell-
Mann matrices do not necessarily anti-commute.

As discussed above in the context of Eq. (116), electronic nematicity in the
2D isotropic system (either in-plane or out-of-plane) belongs to the XY
universality class, implying the existence of a Goldstone mode in the
nematically ordered state. As discussed in Refs. [48, 108], the coupling
between this Goldstone mode and low-energy electronic degrees of freedom can
promote interesting phenomena, such as non-Fermi-liquid behavior. As we will
see below, the main effect of the explicit breaking of the full rotational
symmetry in a lattice described by an axial point group is the emergence of an
anisotropic term in the Landau expansion:

$\mathcal{S}_{\mathrm{axial}}=\int_{\mathsf{x}}\Big{\\{}r_{0}|\boldsymbol{d}^{\mathrm{\alpha}}|^{2}+u|\boldsymbol{d}^{\mathrm{\alpha}}|^{4}+h_{q}|\boldsymbol{d}^{\mathrm{\alpha}}|^{q}\cos\left(q\gamma_{\alpha}\right)\Big{\\}},$
(124)

where $q$ is a positive integer that is the same within the same crystal
system (hexagonal, trigonal, tetragonal, orthorhombic, monoclinic, and
triclinic), $\alpha$ can refer to $\mathrm{ip}$ or $\mathrm{op}$, and $h_{q}$
is a Landau coefficient. For the (crystallographic and non-crystallographic)
axial groups that we investigated here, we found the following general result:
if the group has $(2n)$-fold symmetry, then $q=n$ for $\alpha=\mathrm{ip}$ and
$q=2n$ for $\alpha=\mathrm{op}$. On the other hand, if the axial group has
$(2n+1)$-fold symmetry, $q=2n+1$ for both $\alpha=\mathrm{ip}$ and
$\alpha=\mathrm{op}$.

We note that Eq. (124) corresponds to the Landau expansion of the $q$-state
“soft” clock model ($Z_{q}$-clock), whose classical critical properties are
well understood [109]. In three dimensions (and non-zero temperatures), which
is the case of most relevance here, the character of the transition is (of
course, $q=1$ implies explicitly broken symmetry): 3D-Ising for $q=2$, first-
order for $q=3$, and 3D-XY for $q\geq 4$. The latter result is a consequence
of the fact that $h_{q}$ is a dangerously irrelevant perturbation (in the
renormalization-group sense) for $q\geq 4$, i.e. it is irrelevant at the
nematic critical point but relevant inside the nematically ordered state [110,
111, 112]. Thus, even though the transition is XY-like, the ordered state does
not display a Goldstone mode, but a gapped pseudo-Goldstone mode.

It is important to point out that crystals described by an axial point group
can be very anisotropic and display behavior intermediate between 2D and 3D.
Moreover, as we will discuss in Sec. V.4, twisted quasicrystals are 2D
systems. In this regard, it is interesting to note that the $Z_{q}$-clock
model has rather unique properties in 2D [109]. While the $q=2$ model
undergoes a standard 2D-Ising transition, the $q=3$ model undergoes a second-
order transition despite the existence of a cubic term in the Landau
expansion. More surprisingly, for $q\geq 5$, the 2D $Z_{q}$-clock model
undergoes two Berezinskii-Kosterlitz-Thouless (BKT) transitions: a higher-
temperature BKT transition towards a critical phase with quasi-long-range
order, analogous to that displayed by the XY model at non-zero temperatures,
and a lower-temperature BKT transition towards a state with long-range order
where the discrete clock symmetry is broken. Interestingly, for $q=4$, these
two BKT transitions merge into a single line. As a result, the critical
behavior of the 2D $Z_{4}$-clock model, which maps onto the 2D Ashkin-Teller
model, is only weakly-universal [113], since while the anomalous critical
exponent is fixed to $\eta=1/4$, the other exponents depend on the value of
$h_{4}$ [114]. The crossover from 2D to 3D $Z_{q}$-clock behavior has been
recently discussed numerically in Ref. [112] and also indirectly in the
context of quantum critical points in 2D $Z_{q}$-clock systems [115, 57, 116,
58, 59]. Finally, we emphasize that the critical properties of the
$Z_{q}$-clock model are generically different from those of the $Z_{q}$-Potts
model, except for the special case $q=3$, where they share the same
universality class [104]. For this reason, as explained in the previous
section, we use the terms $Z_{3}$-clock and $Z_{3}$-Potts interchangeably.

Finally, we will show in the Appendices E-G that, in the cases where the
nematic order parameter transforms as a complex IR, which reflects the lack of
twofold in-plane rotation axes in the crystal, the clock term in the action
(124) acquires an offset angle, which can be traced back to the existence of
another anisotropic term in the Landau expansion,
$\tilde{h}_{q}|\boldsymbol{d}^{\mathrm{\alpha}}|^{q}\sin\left(q\gamma_{\alpha}\right)$,
(see also Ref. [102]):

$\mathcal{S}_{\mathrm{axial}}=\int_{\mathsf{x}}\Big{\\{}r_{0}|\boldsymbol{d}^{\mathrm{\alpha}}|^{2}+u|\boldsymbol{d}^{\mathrm{\alpha}}|^{4}+h_{q}|\boldsymbol{d}^{\mathrm{\alpha}}|^{q}\cos\left(q\gamma_{\alpha}-\delta_{0}\right)\Big{\\}}.$
(125)

To distinguish it from Eq. (125), we denote the Landau expansion above a
$Z_{q}^{*}$-clock model. An offset angle is also found when the two nematic
doublets are “degenerate,” i.e., transform as the same IR, as is the case in
trigonal point groups.

###  IV.2 Hexagonal crystals

The hexagonal crystal system consists of seven crystallographic axial point
groups: $\mathsf{D_{6h}}$, $\mathsf{D_{6}}$, $\mathsf{D_{3h}}$,
$\mathsf{C_{6v}}$, $\mathsf{C_{6h}}$, $\mathsf{C_{3h}}$, and $\mathsf{C_{6}}$.
In all cases, the nematic order parameter $\boldsymbol{d}$ in Eq. (17)
decomposes into the three channels in the same way as in the 2D isotropic
case,

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{hexagonal}}{\longrightarrow}\quad\Big{(}\,\underline{d_{1}},\,\boldsymbol{d}^{\mathrm{ip}},\,\boldsymbol{d}^{\mathrm{op}}\,\Big{)}^{T}.$
(126)

We recall that an underline means that the corresponding nematic component
transforms as the trivial IR of the point group. The most notable difference
with respect to the isotropic case is, of course, the limited number of
$C_{nz}$ or $IC_{nz}$ symmetry elements, which reduces the continuous XY-
degeneracy of the nematic angles $\gamma_{\mathrm{ip}}$ and
$\gamma_{\mathrm{op}}$ down to either threefold or sixfold degeneracies, as we
demonstrate below. Within the hexagonal crystal system, one needs to
distinguish between the sets of point groups {$\mathsf{D_{6h}}$,
$\mathsf{D_{6}}$, $\mathsf{D_{3h}}$, $\mathsf{C_{6v}}$} and
{$\mathsf{C_{6h}}$, $\mathsf{C_{3h}}$, $\mathsf{C_{6}}$}, since only the
former possess the characteristic $6$ in-plane symmetry axes

$\displaystyle\mathcal{V}_{1}^{\mathrm{hex}}$
$\displaystyle=\left\\{\left(\begin{smallmatrix}1\\\ 0\\\
0\end{smallmatrix}\right),\left(\begin{smallmatrix}\text{-}\frac{1}{2}\\\
\frac{\sqrt{3}}{2}\\\
0\end{smallmatrix}\right),\left(\begin{smallmatrix}\text{-}\frac{1}{2}\\\
\text{-}\frac{\sqrt{3}}{2}\\\ 0\end{smallmatrix}\right)\right\\},$
$\displaystyle\mathcal{V}_{2}^{\mathrm{hex}}$
$\displaystyle=\left\\{\left(\begin{smallmatrix}0\\\ 1\\\
0\end{smallmatrix}\right),\left(\begin{smallmatrix}\text{-}\frac{\sqrt{3}}{2}\\\
\text{-}\frac{1}{2}\\\
0\end{smallmatrix}\right),\left(\begin{smallmatrix}\frac{\sqrt{3}}{2}\\\
\text{-}\frac{1}{2}\\\ 0\end{smallmatrix}\right)\right\\},$ (127)

see Eqs. (90),(94). Mathematically, the absence of in-plane symmetry axes is
manifested in the nematic doublets $\boldsymbol{d}^{\mathrm{ip}}$ and
$\boldsymbol{d}^{\mathrm{op}}$ transforming according to complex IRs, which
causes additional Landau invariants to emerge when compared to the generic
$Z_{q}$-clock Landau expansion of Eq. (124). As we show in detail in Appendix
E, this additional Landau invariant for the groups {$\mathsf{C_{6h}}$,
$\mathsf{C_{3h}}$, $\mathsf{C_{6}}$} can be recast as an offset angle in the
clock term, see Eq. (125). In this section, we focus instead on the first set
of groups, {$\mathsf{D_{6h}}$, $\mathsf{D_{6}}$, $\mathsf{D_{3h}}$,
$\mathsf{C_{6v}}$}, for which the nematic doublets
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ transform as
real two-dimensional IRs. For concreteness, we show the derivation for the
point group $\mathsf{D_{6h}}$, where $\boldsymbol{d}^{\mathrm{ip}}$ and
$\boldsymbol{d}^{\mathrm{op}}$ transform according to the IRs $E_{2g}$ and
$E_{1g}$, respectively. The results equally apply to the other groups of the
set, {$\mathsf{D_{6}}$, $\mathsf{D_{3h}}$, $\mathsf{C_{6v}}$}.

Consider first the in-plane nematic order parameter
$\boldsymbol{d}^{\mathrm{ip}}=\left(d_{2},\,d_{5}\right)^{T}=\big{(}d_{x^{2}-y^{2}},\,d_{2xy}\big{)}^{T}$.
The symmetrized product decomposition for each expansion order is given by

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{2g},$ (128)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus A_{2g}\oplus E_{2g},$ (129)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\\!\oplus 2E_{2g}.$ (130)

To construct the three $A_{1g}$ Landau invariants, we derive the bilinear
combinations associated with Eq. (128):

$\displaystyle D_{\mathrm{ip}}^{A_{1g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2},$
$\displaystyle\boldsymbol{D}_{\mathrm{ip}}^{E_{2g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left(\cos\left(2\gamma_{\mathrm{ip}}\right),\text{-}\sin\left(2\gamma_{\mathrm{ip}}\right)\right)^{T}.$
(131)

Then, the three invariants become
$\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}$,
$\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{4}$, and
$\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{D}_{\mathrm{ip}}^{E_{2g}}$,
since both $\boldsymbol{d}^{\mathrm{ip}}$ and
$\boldsymbol{D}_{\mathrm{ip}}^{E_{2g}}$ transform as $E_{2g}$. The
corresponding Landau expansion,

$\displaystyle\mathcal{S}_{\mathrm{ip}}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{0}\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}+g_{\mathrm{ip}}\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{3}\cos\left(3\gamma_{\mathrm{ip}}\right)+u\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{4}\Big{\\}},$
(132)

has the same form as that of the $Z_{3}$-clock model, Eq. (124), as derived
elsewhere [41, 74, 57, 79]. Correspondingly, the threefold-degenerate mean-
field ground state is given by

$\displaystyle\gamma_{\mathrm{ip}}^{0}$
$\displaystyle=\frac{\pi}{3}\left(\frac{1+\mathrm{sign}\,g_{\mathrm{ip}}}{2}\right)+\frac{2\pi}{3}n,$
$\displaystyle n$ $\displaystyle=\left\\{0,1,2\right\\}.$ (133)

In the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation (24),
this ground state is parametrized according to Eq. (118); recall that $d_{1}$
is the symmetry-conforming nematic component, which is non-zero at any
temperature or tuning parameter range. We plot the corresponding Fermi surface
distortion in Fig. 5(a) by using the general expression (42). In this figure,
we chose $g_{\mathrm{ip}}<0$ and three values of the trivial component
$d_{1}/\left|\boldsymbol{d}\right|=\left\\{0.1,\,0.5,\,0.9\right\\}$. In the
top panel where
$\left|d_{1}\right|<\left|\boldsymbol{d}^{\mathrm{ip}}\right|/\sqrt{3}$, the
short and long nematic axes $\boldsymbol{n}$ and $\boldsymbol{m}$ align,
respectively, with the in-plane symmetry axes $\mathcal{V}_{1}^{\mathrm{hex}}$
and $\mathcal{V}_{2}^{\mathrm{hex}}$ of Eq. (127), highlighted in purple and
light-blue in the figure. The three degenerate states correspond to
$\boldsymbol{m}$ being aligned with one of the three light-blue axes,
corresponding to $\mathcal{V}_{2}^{\mathrm{hex}}$ in Eq. (127). For the
opposite sign $g_{\mathrm{ip}}>0$, the long nematic axis $\boldsymbol{m}$
aligns with the purple axes, which correspond to
$\mathcal{V}_{1}^{\mathrm{hex}}$ in Eq. (127). The biaxial nature of the
nematic order parameter is clear from the shape of the distorted Fermi
surface. Upon increasing $\left|d_{1}\right|/\left|\boldsymbol{d}\right|$, we
demonstrate in the second panel of Fig. 5(a) how the nematic state passes
through a fine-tuned uniaxial point at
$\left|d_{1}\right|=\left|\boldsymbol{d}\right|/2$, according to Eq. (118). In
the parameter range $\left|d_{1}\right|>\left|\boldsymbol{d}\right|/2$ (third
panel), two nematic axes are interchanged, and the nematic state becomes
gradually more in-plane isotropic.

We now discuss the case of the out-of-plane nematic doublet
$\boldsymbol{d}^{\mathrm{op}}=\left(d_{3},\,d_{4}\right)^{T}=\left(d_{2yz},\,d_{2xz}\right)^{T}$,
which transforms as the $E_{1g}$ IR of $\mathsf{D_{6h}}$. The decomposition of
the symmetrized product gives:

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{2g},$ (134)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{1g}\big{]}_{s}$
$\displaystyle=B_{1g}\oplus B_{2g}\oplus E_{1g},$ (135)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\\!\oplus 2E_{2g},$ (136)
$\displaystyle\big{[}\otimes_{j=1}^{5}E_{1g}\big{]}_{s}$
$\displaystyle=B_{1g}\\!\oplus B_{2g}\oplus 2E_{1g},$ (137)
$\displaystyle\big{[}\otimes_{j=1}^{6}E_{1g}\big{]}_{s}$
$\displaystyle=2A_{1g}\\!\oplus A_{2g}\oplus 2E_{2g}.$ (138)

Here, we extended the expansion to sixth-order to include the leading
anisotropic invariant, i.e. the leading-order invariant that depends on the
angle $\gamma_{\mathrm{op}}$. Using the bilinears associated with Eq. (134),

$\displaystyle D_{\mathrm{op}}^{A_{1g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2},$
$\displaystyle\boldsymbol{D}_{\mathrm{op}}^{E_{2g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\left(\cos\left(2\gamma_{\mathrm{op}}\right),-\sin\left(2\gamma_{\mathrm{op}}\right)\right)^{T},$
(139)

the four invariants can be constructed as
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}$,
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}$,
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{6}$,
$(\boldsymbol{d}^{\mathrm{op}}\cdot\boldsymbol{D}_{\mathrm{op}}^{E_{2g}})^{2}$,
where we used the fact that
$\boldsymbol{d}^{\mathrm{op}}\cdot\boldsymbol{D}_{\mathrm{op}}^{E_{2g}}$
transforms as $B_{1g}$. The resulting nematic Landau expansion can be written
as

$\displaystyle\mathcal{S}_{\mathrm{op}}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{0}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}+u\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}+u_{6}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{6}+v_{6}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{6}\cos\left(6\gamma_{\mathrm{op}}\right)\Big{\\}},$
(140)

and corresponds to the $Z_{6}$-clock model of Eq. (124). Thus, we find an
important difference between the in-plane $\boldsymbol{d}^{\mathrm{ip}}$ and
out-of-plane $\boldsymbol{d}^{\mathrm{op}}$ nematic order parameters in the
hexagonal lattice: whereas the former undergoes a first-order transition, the
latter is expected to undergo a continuous 3D-XY transition, at which the
discrete nature of the nematic angle is irrelevant (in the renormalization
group sense).

The mean-field ground state of $\boldsymbol{d}^{\mathrm{op}}$ is therefore
sixfold degenerate and characterized by the nematic angles

$\displaystyle\gamma_{\mathrm{op}}^{0}$
$\displaystyle=\frac{\pi}{6}\left(\frac{1+\mathrm{sign}\,v_{6}}{2}\right)+\frac{2\pi}{6}n,$
$\displaystyle n$ $\displaystyle\in\\{0,1,\dots,5\\}.$ (141)

The distorted Fermi surface associated with this ground state is shown in Fig.
5(b). In terms of the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation derived in Eq.
(119), we see that, regardless of the value of the trivial nematic component
$d_{1}$, the nematic axis $\boldsymbol{l}$ is always in-plane and aligned with
either $\mathcal{V}_{1}^{\mathrm{hex}}$ (purple) or
$\mathcal{V}_{2}^{\mathrm{hex}}$ (light-blue) defined in Eq. (127), depending
on the sign of the Landau coefficient $v_{6}$. The long ($\boldsymbol{m}$) and
short ($\boldsymbol{n}$) nematic axes are rotated out of the plane by the tilt
angle $\eta_{\mathrm{op}}$ (120), which is exactly $\pm 45^{\circ}$ in the
limiting case $d_{1}=0$. To visualize the sixfold degeneracy of this state, it
is useful to invoke the symmetry elements $C_{2z}$ or $IC_{2z}$, one of which
is present in every hexagonal point group. They imply the existence of a non-
identical ground state configuration that emerges from Fig. 5(b) through a
$180^{\circ}$ rotation about the $k_{z}$-axis. The key point is that, under
such a rotation, $\boldsymbol{l}\rightarrow-\boldsymbol{l}$ but
$\boldsymbol{m}$, $\boldsymbol{n}$ are not mapped onto $-\boldsymbol{m}$,
$-\boldsymbol{n}$. This implies that, for each of the three choices for the
$\boldsymbol{l}$ axis in $\mathcal{V}_{1}^{\mathrm{hex}}$ or
$\mathcal{V}_{2}^{\mathrm{hex}}$, there are two different sets of possible
$\boldsymbol{m}$, $\boldsymbol{n}$ values, resulting in a sixfold degeneracy.
In contrast, such a sixfold degeneracy is absent for the in-plane state
$\boldsymbol{d}^{\mathrm{ip}}$ of Fig. 5(a), since $\boldsymbol{l}$ remains
invariant and $\boldsymbol{m}$, $\boldsymbol{n}$ are mapped onto
$-\boldsymbol{m}$, $-\boldsymbol{n}$ under a $C_{2z}$ or $IC_{2z}$ operation.

As we mentioned earlier, the same results derived here for $\mathsf{D_{6h}}$
hold for {$\mathsf{D_{6}}$, $\mathsf{D_{3h}}$, $\mathsf{C_{6v}}$} and
{$\mathsf{C_{6h}}$, $\mathsf{C_{3h}}$, $\mathsf{C_{6}}$}, as shown in Table 1.
The main difference in the latter set is that an offset angle appears in the
anisotropic terms of Eqs. (132) and (140), see Appendix E, which we indicate
as $Z_{3}^{*}$-clock and $Z_{6}^{*}$-clock models in Table 1.

Figure 5: Distorted Fermi surface in the nematic phase, as given by Eq. (42),
plotted together with the nematic axes
$\boldsymbol{n}$,$\boldsymbol{m}$,$\boldsymbol{l}$ for the cases of: (a)
hexagonal in-plane doublet state (133) for three values of
$d_{1}/\left|\boldsymbol{d}\right|$, and (b) out-of-plane doublet state (141);
(c)-(e) trigonal nematic state (179) for three values of $K_{r}$ (175). The
unit vectors $\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}$
are rescaled by
$2.5\,k_{F}\left(\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}\right)$,
respectively. Unless specified otherwise, we chose
$d_{1}/\left|\boldsymbol{d}\right|=0.1$.

###  IV.3 Trigonal crystals

The axial point groups $\mathsf{D_{3}}$, $\mathsf{D_{3d}}$, $\mathsf{C_{3v}}$,
$\mathsf{S_{6}}$ and $\mathsf{C_{3}}$ form the trigonal crystal system. In
what concerns nematicity, its key distinction with respect to the hexagonal
crystal system, Eq. (126), is that the two nematic doublets, in-plane
$\boldsymbol{d}^{\mathrm{ip}}=\left(d_{2},\,d_{5}\right)^{T}=\big{(}d_{x^{2}-y^{2}},\,d_{2xy}\big{)}^{T}$
and out-of-plane
$\boldsymbol{d}^{\mathrm{op}}=\left(d_{3},\,d_{4}\right)^{T}=\left(d_{2yz},\,d_{2xz}\right)^{T}$,
transform as the same two-dimensional IR. We thus represent the decomposition
of the five-component nematic order parameter $\boldsymbol{d}$ (17) as

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{trigonal}}{\longrightarrow}\quad\Big{(}\,\underline{d_{1}},\,\big{\\{}\boldsymbol{d}^{\mathrm{ip}},\,\boldsymbol{d}^{\mathrm{op}}\big{\\}}\,\Big{)}^{T},$
(142)

where the curly brackets are used to indicate degeneracy, i.e. nematic
components that transform as the same IR. Analogously to the hexagonal case,
we must also distinguish the trigonal group sets {$\mathsf{D_{3}}$,
$\mathsf{D_{3d}}$, $\mathsf{C_{3v}}$} and {$\mathsf{S_{6}}$,
$\mathsf{C_{3}}$}, as the latter lacks the in-plane symmetry directions
encompassed by $\mathcal{V}_{1}^{\mathrm{hex}}$,
$\mathcal{V}_{2}^{\mathrm{hex}}$ in Eq. (127). This case, for which
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ transform as
a complex IR, is discussed in Appendix F. Here, as a representative of the
group set {$\mathsf{D_{3}}$, $\mathsf{D_{3d}}$, $\mathsf{C_{3v}}$}, we focus
our analysis on the point group $\mathsf{C_{3v}}$, for which
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ transform as
the IR $E$. Repeating the same steps as in the other sections, we perform the
symmetrized product decomposition:

$\displaystyle\big{[}\otimes_{j=1}^{2}E\big{]}_{s}$ $\displaystyle=A_{1}\oplus
E,$ $\displaystyle\big{[}\otimes_{j=1}^{3}E\big{]}_{s}$
$\displaystyle=A_{1}\oplus A_{2}\oplus E,$ (143)
$\displaystyle\big{[}\otimes_{j=1}^{4}E\big{]}_{s}$ $\displaystyle=A_{1}\oplus
2E.$ (144)

Clearly, this decomposition is analogous to that of the $E_{2g}$ IR in the
hexagonal case, see Eqs. (128)-(132). This suggests that, individually, each
doublet $\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$
would behave as a $Z_{3}$-clock order parameter. Here, however, we must focus
on the combined four-component vector
$\boldsymbol{d}^{E}=\left(\boldsymbol{d}^{\mathrm{ip}},\boldsymbol{d}^{\mathrm{op}}\right)$,
which transforms as $\left(E\oplus E\right)$. Its symmetrized decomposition,
in turn, gives $13$ invariants up to fourth-order in the Landau expansion:

$\displaystyle\big{[}\otimes_{j=1}^{2}\left(E\oplus E\right)\big{]}_{s}$
$\displaystyle=3A_{1}\oplus A_{2}\oplus 3E,$ (145)
$\displaystyle\big{[}\otimes_{j=1}^{3}\left(E\oplus E\right)\big{]}_{s}$
$\displaystyle=4A_{1}\oplus 4A_{2}\oplus 6E,$ (146)
$\displaystyle\big{[}\otimes_{j=1}^{4}\left(E\oplus E\right)\big{]}_{s}$
$\displaystyle=6A_{1}\oplus 3A_{2}\oplus 13E.$ (147)

To systematically construct them, we use the $7$ bilinears that appear in the
decomposition (145). These include the $4$ bilinears that transform as one-
dimensional IRs

$\displaystyle D_{\mathrm{ip}}^{A_{1}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2},$ $\displaystyle
D_{\mathrm{io}}^{A_{1}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|\cos\left(\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}\right),$
$\displaystyle D_{\mathrm{op}}^{A_{1}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2},$ $\displaystyle
D_{\mathrm{io}}^{A_{2}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|\sin\left(\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}\right),$
(148)

and the $3$ bilinears that transform as two-dimensional IRs

$\displaystyle\boldsymbol{D}_{\mathrm{ip}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left(\begin{smallmatrix}\cos\left(2\gamma_{\mathrm{ip}}\right)\\\
-\sin\left(2\gamma_{\mathrm{ip}}\right)\end{smallmatrix}\right),$
$\displaystyle\boldsymbol{D}_{\mathrm{io}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|\left(\begin{smallmatrix}\cos\left(\gamma_{\mathrm{ip}}+\gamma_{\mathrm{op}}\right)\\\
-\sin\left(\gamma_{\mathrm{ip}}+\gamma_{\mathrm{op}}\right)\end{smallmatrix}\right),$
$\displaystyle\boldsymbol{D}_{\mathrm{op}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\left(\begin{smallmatrix}\cos\left(2\gamma_{\mathrm{op}}\right)\\\
-\sin\left(2\gamma_{\mathrm{op}}\right)\end{smallmatrix}\right).$ (149)

It is now straightforward to construct the $13$ invariants, each associated
with a Landau coefficient. The three quadratic bilinears are given by
$D_{\mathrm{ip}}^{A_{1}}$, $D_{\mathrm{op}}^{A_{1}}$, and
$D_{\mathrm{io}}^{A_{1}}$, whereas the four cubic invariants are
$\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{D}_{\mathrm{ip}}^{E}$,
$\boldsymbol{d}^{\mathrm{op}}\cdot\boldsymbol{D}_{\mathrm{op}}^{E}$,
$\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{D}_{\mathrm{io}}^{E}$, and
$\boldsymbol{d}^{\mathrm{op}}\cdot\boldsymbol{D}_{\mathrm{io}}^{E}$. As for
the quartic invariants, they are given by the squares of the quadratic
invariants $\big{(}D_{\mathrm{ip}}^{A_{1}}\big{)}^{2}$ and
$\big{(}D_{\mathrm{op}}^{A_{1}}\big{)}^{2}$, the cross terms
$D_{\mathrm{ip}}^{A_{1}}D_{\mathrm{op}}^{A_{1}}$,
$D_{\mathrm{ip}}^{A_{1}}D_{\mathrm{io}}^{A_{1}}$ and
$D_{\mathrm{op}}^{A_{1}}D_{\mathrm{io}}^{A_{1}}$, as well as
$\big{(}D_{\mathrm{io}}^{A_{1}}\big{)}^{2}-\big{(}D_{\mathrm{io}}^{A_{2}}\big{)}^{2}$.
Writing the action as
$\mathcal{S}=\mathcal{S}_{2}+\mathcal{S}_{3}+\mathcal{S}_{4}$, we thus have:

$\displaystyle\mathcal{S}_{2}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{\mathrm{ip}}\left(\boldsymbol{d}^{\mathrm{ip}}\right)^{2}+r_{\mathrm{op}}\left(\boldsymbol{d}^{\mathrm{op}}\right)^{2}+r_{\mathrm{io}}\left(\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{d}^{\mathrm{op}}\right)\Big{\\}},$
(150) $\displaystyle\mathcal{S}_{3}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}g_{\mathrm{ip}}\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{3}\mathsf{c}_{3\gamma_{\mathrm{ip}}}+g_{1}\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|\mathsf{c}_{2\gamma_{\mathrm{ip}}+\gamma_{\mathrm{op}}}$
$\displaystyle+g_{\mathrm{op}}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{3}\mathsf{c}_{3\gamma_{\mathrm{op}}}+g_{2}\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\mathsf{c}_{\gamma_{\mathrm{ip}}+2\gamma_{\mathrm{op}}}\Big{\\}},$
(151) $\displaystyle\mathcal{S}_{4}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}u_{\mathrm{ip}}\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{4}+\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|\left[u_{1}\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}+u_{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\right]\mathsf{c}_{\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}}$
$\displaystyle+u_{\mathrm{op}}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}+\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\left[u_{\mathrm{io}}^{0}+u_{\mathrm{io}}^{\mathrm{c}}\mathsf{c}_{2\gamma_{\mathrm{ip}}-2\gamma_{\mathrm{op}}}\right]\Big{\\}}.$
(152)

To avoid cumbersome notations, we defined
$\mathsf{c}_{\gamma}\equiv\cos\gamma$. As expected, the degeneracy between
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ is
manifested in the quadratic mixing term
$\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{d}^{\mathrm{op}}$.

To determine the mean-field ground state, we first diagonalize the quadratic
part $\mathcal{S}_{2}$ (150) to determine the combination of
$\left(\boldsymbol{d}^{\mathrm{ip}},\boldsymbol{d}^{\mathrm{op}}\right)$ that
orders first. We find:

$\displaystyle\mathcal{S}_{2}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}\lambda_{+}\left|\boldsymbol{d}^{+}\right|^{2}+\lambda_{-}\left|\boldsymbol{d}^{-}\right|^{2}\Big{\\}},$
(153)

with the eigenvalues

$\displaystyle\lambda_{\pm}$
$\displaystyle=\frac{1}{2}\left(r_{\mathrm{ip}}+r_{\mathrm{op}}\pm\sqrt{\left(r_{\mathrm{ip}}-r_{\mathrm{op}}\right)^{2}+r_{\mathrm{io}}^{2}}\,\right),$
(154)

and the eigenvectors

$\displaystyle\left(\begin{array}[]{c}\boldsymbol{d}^{+}\\\
\boldsymbol{d}^{-}\end{array}\right)$
$\displaystyle=U^{T}\left(\begin{array}[]{c}\boldsymbol{d}^{\mathrm{ip}}\\\
\boldsymbol{d}^{\mathrm{op}}\end{array}\right)=\left(\begin{array}[]{c}\mathrm{sign}\left(r_{\mathrm{io}}\right)\beta_{+}\boldsymbol{d}^{\mathrm{ip}}+\beta_{-}\boldsymbol{d}^{\mathrm{op}}\\\
-\mathrm{sign}\left(r_{\mathrm{io}}\right)\beta_{-}\boldsymbol{d}^{\mathrm{ip}}+\beta_{+}\boldsymbol{d}^{\mathrm{op}}\end{array}\right).$
(161)

In these expressions, we defined

$\displaystyle\beta_{\pm}$
$\displaystyle=\frac{1}{\sqrt{2}}\sqrt{1\pm\frac{r_{\mathrm{ip}}-r_{\mathrm{op}}}{\sqrt{\left(r_{\mathrm{ip}}-r_{\mathrm{op}}\right)^{2}+r_{\mathrm{io}}^{2}}}}\,,$
(162)

and the unitary matrix

$\displaystyle U$
$\displaystyle=\left(\begin{smallmatrix}\beta_{+}\mathrm{sign}r_{\mathrm{io}}&0&-\beta_{-}\mathrm{sign}r_{\mathrm{io}}&0\\\
0&\beta_{+}\mathrm{sign}r_{\mathrm{io}}&0&-\beta_{-}\mathrm{sign}r_{\mathrm{io}}\\\
\beta_{-}&0&\beta_{+}&0\\\ 0&\beta_{-}&0&\beta_{+}\end{smallmatrix}\right).$
(163)

By construction, $\lambda_{-}<\lambda_{+}$, which implies that the order
parameter combination $\boldsymbol{d}^{-}$ is the one that condenses. The
order parameter combination $\boldsymbol{d}^{+}$, on the other hand, is a
fluctuating field that primarily renormalizes the action for
$\boldsymbol{d}^{-}$. Therefore, to proceed, we consider only terms that are
quadratic or linear in $\boldsymbol{d}^{+}$ (Gaussian approximation).
Introducing the parametrization
$\boldsymbol{d}^{\pm}=\left|\boldsymbol{d}^{\pm}\right|\left(\cos\gamma_{\pm},\sin\gamma_{\pm}\right)$,
we can rewrite the action as $\mathcal{S}=\mathcal{S}_{-}+\mathcal{S}_{+-}$
with

$\displaystyle\mathcal{S}_{-}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}\lambda_{-}\left|\boldsymbol{d}^{-}\right|^{2}+g_{-}\left|\boldsymbol{d}^{-}\right|^{3}\cos\left(3\gamma_{-}\right)+u_{-}\left|\boldsymbol{d}^{-}\right|^{4}\Big{\\}},$
(164) $\displaystyle\mathcal{S}_{+-}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}\lambda_{+}\left|\boldsymbol{d}^{+}\right|^{2}+\tilde{g}_{2}\left|\boldsymbol{d}^{+}\right|\left|\boldsymbol{d}^{-}\right|^{2}\cos\left(\gamma_{+}+2\gamma_{-}\right)\Big{\\}}.$
(165)

Here, the Landau coefficients $g_{-}$, $\tilde{g}_{2}$, and $u_{-}$ can in
principle be expressed in terms of the original Landau coefficients of Eqs.
(150)-(152) by applying the unitary transformation $U$ of Eq. (163). Since the
$\boldsymbol{d}^{+}$ action is Gaussian, the minimization of $\mathcal{S}_{+}$
is straightforward and gives $\boldsymbol{d}^{+}$ in terms of
$\boldsymbol{d}^{-}$:

$\displaystyle\left|\boldsymbol{d}^{+}\right|$
$\displaystyle=\frac{\left|\tilde{g}_{2}\right|}{2\lambda_{+}}\left|\boldsymbol{d}^{-}\right|^{2},$
$\displaystyle\gamma_{+}$
$\displaystyle=-2\gamma_{-}+\left(\frac{1+\mathrm{sign}\,\tilde{g}_{2}}{2}\right)\pi.$
(166)

Substituting it back in $\mathcal{S}$, we find the Landau expansion in terms
of the ordering field $\boldsymbol{d}^{-}$ alone:

$\mathcal{S}=\int_{\mathsf{x}}\Big{\\{}\lambda_{-}\left|\boldsymbol{d}^{-}\right|^{2}+g_{-}\left|\boldsymbol{d}^{-}\right|^{3}\cos\left(3\gamma_{-}\right)+\tilde{u}_{-}\left|\boldsymbol{d}^{-}\right|^{4}\Big{\\}},$
(167)

with the reduced quartic coefficient
$\tilde{u}_{-}=u_{-}-\frac{\tilde{g}_{2}^{2}}{4\lambda_{+}}$.

The central result is that the nematic action (167) has the same form as the
Landau expansion of the $Z_{3}$-clock model, see Eq. (124) and Ref. [107].
Therefore, in all trigonal crystals, the nematic transition belongs to the 3D
$Z_{3}$-Potts/clock “universality class”, which actually corresponds to a
mean-field first-order transition. Indeed, minimization of Eq. (167) reveals
three degenerate ground states:

$\displaystyle\gamma_{-}^{0}$
$\displaystyle=\frac{\pi}{3}\left(\frac{1+\mathrm{sign}\,g_{-}}{2}\right)+\frac{2\pi}{3}n,$
$\displaystyle n$ $\displaystyle\in\\{0,1,2\\}.$ (168)

Substitution in Eq. (166) shows that, in the ground state,
$\boldsymbol{d}^{+}$ and $\boldsymbol{d}^{-}$ are collinear:

$\left(\begin{array}[]{c}\cos\gamma_{+}^{0}\\\
\sin\gamma_{+}^{0}\end{array}\right)=\mathrm{sign}\left(g_{-}\tilde{g}_{2}\right)\left(\begin{array}[]{c}\cos\gamma_{-}^{0}\\\
\sin\gamma_{-}^{0}\end{array}\right).$ (169)

Therefore, using Eq. (161), we find that the nematic in-plane and out-of-plane
doublets are also collinear, and thus can be parametrized as:

$\displaystyle\boldsymbol{d}^{\mathrm{ip}}$
$\displaystyle=\left|\boldsymbol{d}^{E}\right|\cos\delta_{E}\left(\begin{array}[]{c}\cos\gamma_{-}^{0}\\\
\sin\gamma_{-}^{0}\end{array}\right),$
$\displaystyle\boldsymbol{d}^{\mathrm{op}}$
$\displaystyle=\left|\boldsymbol{d}^{E}\right|\sin\delta_{E}\left(\begin{array}[]{c}\cos\gamma_{-}^{0}\\\
\sin\gamma_{-}^{0}\end{array}\right),$ (174)

where $\delta_{E}\in[0,\pi]$ and
$\left|\boldsymbol{d}^{E}\right|=\sqrt{\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}+\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}}$
are determined by the Landau coefficients.

To express the trigonal nematic ground state (174) in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation (24), it is
convenient to define two new quantities (recall the $d_{1}$ is the non-zero
symmetry-preserving nematic component)

$\displaystyle K_{r}$
$\displaystyle=-\frac{\sqrt{3}\,\mathrm{sign}(g_{-})}{2}\cos\left(\delta_{E}\right)\frac{\left|\boldsymbol{d}^{E}\right|}{\left|\boldsymbol{d}\right|}-\frac{d_{1}}{2\left|\boldsymbol{d}\right|}$
$\displaystyle\in$ $\displaystyle\left[-1,1\right],$ (175) $\displaystyle
K_{r}^{\prime}$
$\displaystyle=-\frac{\mathrm{sign}(g_{-})}{2}\cos\left(\delta_{E}\right)\frac{\left|\boldsymbol{d}^{E}\right|}{\left|\boldsymbol{d}\right|}+\frac{\sqrt{3}}{2}\frac{d_{1}}{\left|\boldsymbol{d}\right|}$
$\displaystyle\in$ $\displaystyle\left[-1,1\right],$ (176)

as well as the tilt angle $\eta_{0}\in\left[0,\frac{\pi}{2}\right]$

$\displaystyle\eta_{0}$
$\displaystyle=\frac{\pi}{2}\frac{1+\mathrm{sign}K_{r}^{\prime}}{2}-\frac{\mathrm{sign}K_{r}^{\prime}}{2}\arccos\Bigg{(}\frac{\left|K_{r}^{\prime}\right|}{\sqrt{1-\left|K_{r}\right|^{2}}}\Bigg{)},$
(177)

and the three orthonormal vectors

$\displaystyle\boldsymbol{e}_{\parallel}$
$\displaystyle=\left(\begin{smallmatrix}\cos\gamma_{-}^{0}\\\
-\sin\gamma_{-}^{0}\\\ 0\end{smallmatrix}\right),$
$\displaystyle\boldsymbol{e}_{\perp}^{A}$
$\displaystyle=\left(\begin{smallmatrix}\sin\gamma_{-}^{0}\cos\eta_{0}\\\
\cos\gamma_{-}^{0}\cos\eta_{0}\\\ \sin\eta_{0}\end{smallmatrix}\right),$
$\displaystyle\boldsymbol{e}_{\perp}^{B}$
$\displaystyle=\left(\begin{smallmatrix}\sin\gamma_{-}^{0}\sin\eta_{0}\\\
\cos\gamma_{-}^{0}\sin\eta_{0}\\\ -\cos\eta_{0}\end{smallmatrix}\right).$
(178)

Note that $\boldsymbol{e}_{\parallel}$ points along one of the six high-
symmetry directions of the hexagonal point groups, given by
$\mathcal{V}_{1}^{\mathrm{hex}}$ and $\mathcal{V}_{2}^{\mathrm{hex}}$ in Eq.
(127). Indeed, the trigonal nematic angle $\gamma_{-}^{0}$ in Eq. (168)
assumes the same values as the in-plane hexagonal nematic angles
$\gamma_{\mathrm{ip}}^{0}$ in Eq. (133).

In terms of these quantities, the trigonal nematic state (174) can be
conveniently expressed in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-notation via

$\displaystyle\left|K_{r}\right|\leq\frac{1}{2}$ $\displaystyle\\!:$
$\displaystyle\\!\\!\alpha$
$\displaystyle=\frac{\pi}{6}+\arcsin\left(K_{r}\right),$
$\displaystyle\\!\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{\perp}^{A},\;\;\boldsymbol{m}=\boldsymbol{e}_{\perp}^{B},\;\;\boldsymbol{l}=\boldsymbol{e}_{\parallel},$
$\displaystyle K_{r}<\frac{\text{-}1}{2}$ $\displaystyle\\!:$
$\displaystyle\\!\\!\alpha$
$\displaystyle=\arcsin\left(\left|K_{r}\right|\right)-\frac{\pi}{6},$
$\displaystyle\\!\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{\perp}^{A},\;\;\boldsymbol{m}=\boldsymbol{e}_{\parallel},\;\;\boldsymbol{l}=\boldsymbol{e}_{\perp}^{B},$
$\displaystyle K_{r}>\frac{1}{2}$ $\displaystyle\\!:$
$\displaystyle\\!\\!\alpha$
$\displaystyle=\frac{\pi}{2}-\arcsin\left(\left|K_{r}\right|\right),$
$\displaystyle\\!\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{\parallel},\;\;\boldsymbol{m}=\boldsymbol{e}_{\perp}^{B},\;\;\boldsymbol{l}=\boldsymbol{e}_{\perp}^{A}.$
(179)