The similarity to the in-plane nematic state of the 2D isotropic system, Eq.
(118), is apparent – indeed, one can verify that upon setting $\delta_{E}=0$
(pure in-plane nematicity) one recovers the exact same equations. Conversely,
upon setting $\delta_{E}=\pi/2$ (pure out-of-plane nematicity), the out-of-
plane isotropic nematic solution (119) is recovered. Therefore, the trigonal
nematic state is generally biaxial with $\alpha\neq 0,\,\pi/3$, except for the
special point $\left|K_{r}\right|=\frac{1}{2}$, where it becomes uniaxial. In
Figs. 5(c)-(e), we plot the Fermi surface distortion in the trigonal nematic
state (179) for three representative values of $K_{r}$, encompassing each of
the three regimes above. For $\left|K_{r}\right|<\frac{1}{2}$, shown in Fig.
5(e), the nematic distortion is qualitatively similar to that caused by an
out-of-plane hexagonal order parameter $\boldsymbol{d}^{\mathrm{op}}$ [Fig.
5(b)], since the nematic axis $\boldsymbol{l}=\boldsymbol{e}_{\parallel}$
aligns with one of the three in-plane symmetry axes, represented by the purple
lines in the figure [corresponding to $\mathcal{V}_{1}^{\mathrm{hex}}$ in Eq.
(127)]. The key difference is that the trigonal groups do not possess the
symmetry elements $C_{2z}$ or $IC_{2z}$, such that a $180^{\circ}$ degree
rotation about the $k_{z}$-axis does not lead to a degenerate state. As a
result, what used to be a $Z_{6}$-degenerate state in the hexagonal case [Fig.
5(b)] splits into two sets of $Z_{3}$-degenerate states in the trigonal case,
one for each sign of $g_{-}$ [Eq. (168)]. In the regime
$\left|K_{r}\right|>\frac{1}{2}$, plotted in Fig. 5(c), the trigonal nematic
state has either the long ($\boldsymbol{m}$) or the short ($\boldsymbol{n}$)
nematic axis aligned with an in-plane symmetry axis
$\boldsymbol{e}_{\parallel}$ , which for $g_{-}<0$ corresponds to
$\mathcal{V}_{1}^{\mathrm{hex}}$. This state resembles the in-plane nematic
state of the hexagonal lattice depicted in Fig. 5(a), albeit slightly rotated
out-of-plane. Exactly at $\left|K_{r}\right|=\frac{1}{2}$, $\alpha=0,\,\pi/3$
and the nematic state is uniaxial, as shown in Fig. 5(d). The properties of
the nematic state of all trigonal groups are summarized in Table 1.

###  IV.4 Tetragonal crystals

In any of the seven tetragonal point groups, the five-component nematic order
parameter $\boldsymbol{d}$ (17) is decomposed as:

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{tetragonal}}{\longrightarrow}\quad\Big{(}\,\underline{d_{1}},\,d_{2},\,d_{5},\,\boldsymbol{d}^{\mathrm{op}}\,\Big{)}^{T}.$
(180)

Hence, in contrast to the 2D isotropic system (97), as well as to the
hexagonal and trigonal lattices, in tetragonal crystals the nematic in-plane
doublet $\boldsymbol{d}^{\mathrm{ip}}=\left(d_{2},d_{5}\right)^{T}$ is further
decomposed into the two channels $d_{2}=d_{x^{2}-y^{2}}$ and $d_{5}=d_{2xy}$,
which transform as non-trivial one-dimensional IRs. Thus, the tetragonal
crystal system is the highest-symmetry crystal for which a non-trivial nematic
component exists that does not transform as a multi-dimensional IR.

There are seven different axial point groups in the tetragonal crystal system:
{$\mathsf{D_{4h}}$, $\mathsf{D_{4}}$, $\mathsf{C_{4v}}$, $\mathsf{D_{2d}}$},
which possess in-plane twofold rotational symmetry axes, and
{$\mathsf{C_{4h}}$, $\mathsf{C_{4}}$, $\mathsf{S_{4}}$}, which do not possess
such axes. As in the previous subsections, the two sets need to be treated
slightly differently, since in the latter the out-of-plane nematic order
parameter $\boldsymbol{d}^{\mathrm{op}}$ transforms as a complex IR and
$\left\\{d_{2},d_{5}\right\\}$ transform as the same one-dimensional IR. In
this section we will focus on the first set, illustrating the results for the
point group $\mathsf{D_{4h}}$, while leaving the discussion of the second set
of point groups to Appendix G.

In $\mathsf{D_{4h}}$, the single-component order parameters
$d_{2}=d_{x^{2}-y^{2}}$ and $d_{5}=d_{2xy}$ transform according to the IRs
$B_{1g}$ and $B_{2g}$, respectively, while the out-of-plane doublet
$\boldsymbol{d}^{\mathrm{op}}=\left(d_{3},\,d_{4}\right)^{T}=\left(d_{2yz},\,d_{2xz}\right)^{T}$
transforms as the IR $E_{g}$. As single-component real-valued order
parameters, both $d_{2}$ and $d_{5}$ are Ising variables that undergo a
$Z_{2}$-Ising transition described by the action

$\displaystyle\mathcal{S}_{\nu}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{\nu}\left|d_{\nu}\right|^{2}+u_{\nu}\left|d_{\nu}\right|^{4}\Big{\\}},$
(181)

with $\nu=\\{2,5\\}$. We note that this Landau expansion could also be recast
in terms of the in-plane nematic doublet
$\boldsymbol{d}^{\mathrm{ip}}=\left(d_{2},\,d_{5}\right)^{T}=\big{(}d_{x^{2}-y^{2}},\,d_{2xy}\big{)}^{T}$
as an effective $Z_{2}$-clock model, see Eq. (124):

$\mathcal{S}_{\mathrm{ip}}=\int_{\mathsf{x}}\Big{\\{}r_{0}|\boldsymbol{d}^{\mathrm{ip}}|^{2}+h_{2}|\boldsymbol{d}^{\mathrm{ip}}|^{2}\cos\left(2\gamma_{\mathrm{ip}}\right)+\mathcal{O}\left(|\boldsymbol{d}^{\mathrm{ip}}|^{4}\right)\Big{\\}}.$
(182)

The second term is nothing but $h_{2}\big{(}d_{2}^{2}-d_{5}^{2}\big{)}$, which
implies that the transitions toward $d_{2}$ and $d_{5}$ electronic nematic
order take place at different values of the control parameter, $r_{0}=-h_{2}$
for $d_{2}$ and $r_{0}=h_{2}$ for $d_{5}$. The symmetry here also allows for
an additional biquadratic term $d_{2}^{2}d_{5}^{2}$. One could in principle
integrate out the fluctuations of the sub-leading channel, whose effect on the
leading channel will only be important if $h_{2}$ is small. Thus, the Landau
coefficients $r_{\nu}$ and $u_{\nu}$ in Eq. (181) should be understood as
renormalized Landau coefficients of the “original” in–plane nematic action.

The $d_{2}=d_{x^{2}-y^{2}}$ and $d_{5}=d_{2xy}$ Ising-nematic order parameters
are the types of nematic order most widely studied in the literature [64, 65,
66, 67, 68, 69, 70, 71, 39, 72, 73]. Whereas one often describes them in terms
of an effective 2D nematic order parameter, see Eq. (121), it is
straightforward to describe them in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation by employing Eq.
(118) with the additional constraint of $\gamma_{\mathrm{ip}}=\\{0,\pi\\}$ for
$d_{2}$ [corresponding to $h_{2}<0$ in Eq. (182)] and
$\gamma_{\mathrm{ip}}=\\{\frac{\pi}{2},\frac{3\pi}{2}\\}$ for $d_{5}$
[corresponding to $h_{2}>0$ in Eq. (182)]. The corresponding Fermi surface
distortions are shown in Figs. 6(a),(b). For the $d_{2}=d_{x^{2}-y^{2}}$
state, the in-plane nematic axes align with the in-plane principal axes
$\mathcal{V}_{1}^{\mathrm{tet}}$ (light-blue) while for $d_{5}=d_{2xy}$ they
align with the in-plane diagonal axes $\mathcal{V}_{2}^{\mathrm{tet}}$
(purple), with

$\displaystyle\mathcal{V}_{1}^{\mathrm{tet}}$
$\displaystyle=\left\\{\left(\begin{smallmatrix}1\\\ 0\\\
0\end{smallmatrix}\right),\left(\begin{smallmatrix}0\\\ 1\\\
0\end{smallmatrix}\right)\right\\},$
$\displaystyle\mathcal{V}_{2}^{\mathrm{tet}}$
$\displaystyle=\left\\{\frac{1}{\sqrt{2}}\left(\begin{smallmatrix}1\\\ 1\\\
0\end{smallmatrix}\right),\frac{1}{\sqrt{2}}\left(\begin{smallmatrix}1\\\
-1\\\ 0\end{smallmatrix}\right)\right\\}.$ (183)

Therefore, the $Z_{2}$ character of the Ising-nematic order parameter is
associated with the $C_{4z}$ operation ($90^{\circ}$ rotation about the
$k_{z}$-axis). Recall that the alignment of the $\boldsymbol{n}$,
$\boldsymbol{m}$, $\boldsymbol{l}$ axes in Eq. (118), which determines which
nematic axes are in-plane, depends on the ratio
$d_{1}/\left|\boldsymbol{d}\right|$ between the trivial (i.e. symmetry-
preserving) nematic component
$d_{1}=d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}$ and
$\left|\boldsymbol{d}\right|=\sqrt{d_{\nu}^{2}+d_{1}^{2}}$, with
$\nu=\\{2,5\\}$.

To derive the Landau invariants for the out-of-plane nematic order parameter
$\boldsymbol{d}^{\mathrm{op}}$ we compute, once again, the decomposition of
the symmetrized product for each expansion order,

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus B_{1g}\oplus B_{2g},$ (184)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{g}\big{]}_{s}$
$\displaystyle=2E_{g},$ (185)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{g}\big{]}_{s}$
$\displaystyle=2A_{1g}\\!\oplus A_{2g}\oplus B_{1g}\oplus B_{2g}.$ (186)

From the bilinears associated with Eq. (184),
$D_{\mathrm{op}}^{A_{1g}}=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}$, and

$\displaystyle D_{\mathrm{op}}^{B_{1g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\cos\left(2\gamma_{\mathrm{op}}\right),$
$\displaystyle D_{\mathrm{op}}^{B_{2g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\sin\left(2\gamma_{\mathrm{op}}\right),$
(187)

we can readily identify the three Landau invariants as
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}$,
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}$, and
$\big{(}D_{\mathrm{op}}^{B_{1g}}\big{)}^{2}-\big{(}D_{\mathrm{op}}^{B_{2g}}\big{)}^{2}$.
We then obtain the Landau expansion

$\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}+v_{4}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}\cos\left(4\gamma_{\mathrm{op}}\right)\Big{\\}},$
(188)

which corresponds to the $Z_{4}$-clock model of Eq. (124). The fourfold
degenerate ground-state is parametrized by

$\displaystyle\gamma_{\mathrm{op}}^{0}$
$\displaystyle=\frac{\pi}{4}\left(\frac{1+\mathrm{sign}v_{4}}{2}\right)+\frac{2\pi}{4}n,$
$\displaystyle n$ $\displaystyle\in\\{0,1,2,3\\}.$ (189)

Therefore, the out-of-plane nematic transition in a tetragonal crystal belongs
to the 3D-XY universality class; we emphasize that such a critical behavior
may only emerge very close to the transition, depending on the degree of out-
of-plane anisotropy of the system.

The corresponding representation of Eq. (189) in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-formalism can be directly read
off from Eq. (119). Fig. 6(c) shows the corresponding Fermi surface
distortion, with the nematic axis $\boldsymbol{l}$ aligned with the in-plane
high-symmetry axes of either $\mathcal{V}_{1}^{\mathrm{tet}}$ (light-blue) or
$\mathcal{V}_{2}^{\mathrm{tet}}$(purple) [see Eq. (183)], depending on the
sign of $v_{4}$. The fourfold degeneracy of the ground state arises from the
two axes in either $\mathcal{V}_{1}^{\mathrm{tet}}$ or
$\mathcal{V}_{2}^{\mathrm{tet}}$ with each $\boldsymbol{l}$ can align,
combined with the two sets of directions of $\boldsymbol{n},\,\boldsymbol{m}$
that are related by a $C_{2z}$ symmetry operation (which is present in all
tetragonal groups). A summary of the results for all tetragonal crystal groups
is presented in Table 1, combining the analysis of this section with that of
Appendix G.

Figure 6: 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)-(b)
tetragonal in-plane $d_{2}=d_{x^{2}-y^{2}}$ and $d_{5}=d_{2xy}$ states (181),
(c) tetragonal out-of-plane doublet state $\boldsymbol{d}^{\mathrm{op}}$
(189), (d)-(f) orthorhombic nematic states $d_{5}=d_{2xy}$, $d_{3}=d_{2yz}$
and $d_{4}=d_{2xz}$. The unit vectors
$\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}$ are rescaled
by
$2.25\,k_{F}\left(\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}\right)$,
respectively. For the tetragonal and orthorhombic cases we chose the trivial
components as $d_{1}/\left|\boldsymbol{d}\right|=0.1$ and
$\left(\left|\boldsymbol{d}^{e}\right|/\left|\boldsymbol{d}\right|,\gamma_{e}\right)=\left(0.3,1.1\right)$,
respectively.

###  IV.5 Orthorhombic, Monoclinic, and Triclinic crystals

In the remaining three crystal systems – orthorhombic, monoclinic, and
triclinic – all the point groups are Abelian, and as such do not admit multi-
dimensional IRs (real or complex). As a result, all nematic components
transform as one-dimensional IRs, implying that the non-trivial nematic
components behave as $Z_{2}$-Ising order parameters. Another key difference
between these axial groups and the other ones previously analyzed in this
section is that at least one of the components of
$\boldsymbol{d}^{\mathrm{ip}}=\left(d_{2},\,d_{5}\right)^{T}=\big{(}d_{x^{2}-y^{2}},\,d_{2xy}\big{)}^{T}$
and
$\boldsymbol{d}^{\mathrm{op}}=\left(d_{3},\,d_{4}\right)^{T}=\left(d_{2yz},\,d_{2xz}\right)^{T}$
transform as the trivial IR, like
$d_{1}=d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}$.

We start with the orthorhombic crystal system, described by the $3$ point
groups {$\mathsf{D_{2h}}$, $\mathsf{D_{2}}$, $\mathsf{C_{2v}}$}. In terms of
the IR of these groups, the five nematic components $\boldsymbol{d}$ (17) are
decomposed according to

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{orthorombic}}{\longrightarrow}\quad\Big{(}\,\underline{d_{1},\,d_{2}},\,d_{3},\,d_{4},\,d_{5}\,\Big{)}^{T}.$
(190)

The first two components $\left\\{d_{1},d_{2}\right\\}$, highlighted by the
underline, transform trivially under the groups symmetry operations and thus
are generically non-zero. The remaining three components $d_{3}=d_{2yz}$,
$d_{4}=d_{2xz}$ and $d_{5}=d_{2xy}$ correspond to shear distortions of the
Fermi surface, and transform according to three different one-dimensional IRs,
as shown explicitly in Table 1. As such, they all undergo a $Z_{2}$-Ising
transition described by an action analogous to Eq. (181).

These three Ising-nematic states can be expressed in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation provided that
one also includes the non-zero trivial components
$\boldsymbol{d}^{e}=\left(d_{1},d_{2}\right)=\left|\boldsymbol{d}^{e}\right|\left(\cos\gamma_{e},\sin\gamma_{e}\right)$.
In the case of the nematic order parameter $d_{5}=d_{2xy}$, one can directly
apply the in-plane isotropic parametrization of Eq. (118) with
$\left|\boldsymbol{d}\right|=\sqrt{d_{5}^{2}+\left|\boldsymbol{d}^{e}\right|^{2}}$
and $\tan\gamma_{\mathrm{ip}}=d_{5}/d_{2}$. The resulting Fermi surface
distortion, calculated via Eq. (42), is shown in Fig. 6(d). Since $d_{2}$ is
always non-zero, the in-plane nematic axes $\boldsymbol{e}_{\parallel}^{A}$
and $\boldsymbol{e}_{\parallel}^{B}$ defined in Eq. (117) are never aligned
with $\mathcal{V}_{2}^{\mathrm{tet}}$ in Eq. (183), which is consistent with
the fact that $\mathcal{V}_{2}^{\mathrm{tet}}$ are not high-symmetry
directions in orthorhombic crystals. Note that one of the nematic axes is
aligned with the $k_{z}$-axis. The $Z_{2}$ degeneracy is associated with the
$C_{2x}$ ($C_{2y}$) or $IC_{2x}$ ($IC_{2y}$) symmetry operations,
corresponding to in-plane twofold rotations or vertical mirror reflections of
the nematic axes in Fig. 6(d).

As for the nematic order parameters $d_{3}$ and $d_{4}$, one cannot just apply
the out-of-plane parametrization of Eq. (119), since $d_{2}$, which is a
component of $\boldsymbol{d}^{\mathrm{ip}}$, is always non-zero. To express
$d_{3}$ and $d_{4}$ in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation, it is
convenient to define the auxiliary variables.

$\displaystyle c_{k}^{e}$
$\displaystyle=\frac{\left|\boldsymbol{d}^{e}\right|}{\left|\boldsymbol{d}\right|}\cos\left(\gamma_{e}-\frac{k\pi}{3}\right),$
$\displaystyle s_{k}^{e}$
$\displaystyle=\frac{\left|\boldsymbol{d}^{e}\right|}{\left|\boldsymbol{d}\right|}\sin\left(\gamma_{e}-\frac{k\pi}{3}\right),$
(191)

with $k=\\{2,4\\}$ and
$\left|\boldsymbol{d}\right|=\sqrt{\left|\boldsymbol{d}^{e}\right|^{2}+(d_{3})^{2}+(d_{4})^{2}}$.
We then define the angles

$\displaystyle\eta_{3}^{\pm}$
$\displaystyle=\frac{1}{2}\frac{d_{3}}{\left|d_{3}\right|}\arccos\Bigg{(}\frac{\left|s_{2}^{e}\right|}{\sqrt{1-\left(c_{2}^{e}\right)^{2}}}\Bigg{)}\pm\frac{1+\mathrm{sign}\left(s_{2}^{e}\right)}{2}\frac{\pi}{2},$
(192) $\displaystyle\eta_{4}^{\pm}$
$\displaystyle=\frac{1}{2}\frac{d_{4}}{\left|d_{4}\right|}\arccos\Bigg{(}\frac{\left|s_{4}^{e}\right|}{\sqrt{1-\left(c_{4}^{e}\right)^{2}}}\Bigg{)}\pm\frac{1+\mathrm{sign}\left(s_{4}^{e}\right)}{2}\frac{\pi}{2},$
(193)

and, from them, the four unit vectors

$\displaystyle\boldsymbol{e}_{A}^{yz}$
$\displaystyle=\left(0,\sin\eta_{3}^{+},\cos\eta_{3}^{-}\right)^{T}\\!\\!,$
$\displaystyle\boldsymbol{e}_{B}^{yz}$
$\displaystyle=\left(0,\cos\eta_{3}^{-},\text{-}\sin\eta_{3}^{+}\right)^{T}\\!\\!,$
$\displaystyle\boldsymbol{e}_{A}^{zx}$
$\displaystyle=\left(\cos\eta_{4}^{-},0,\sin\eta_{4}^{+}\right)^{T}\\!\\!,$
$\displaystyle\boldsymbol{e}_{B}^{zx}$
$\displaystyle=\left(\text{-}\sin\eta_{4}^{+},0,\cos\eta_{4}^{-}\right)^{T}\\!\\!.$
(194)

Now, the nematic order parameter $d_{3}=d_{2yz}$ can be conveniently expressed
in the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation via:

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

Analogously, the nematic order parameter $d_{4}=d_{2xz}$ becomes:

$\displaystyle\left|c_{4}^{e}\right|\leq\frac{1}{2}$ $\displaystyle:$
$\displaystyle\\!\alpha$
$\displaystyle=\frac{\pi}{6}+\arcsin\left(c_{4}^{e}\right),$
$\displaystyle\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{A}^{zx},\;\boldsymbol{m}=\boldsymbol{e}_{B}^{zx},\;\boldsymbol{l}=\boldsymbol{e}_{y},$
$\displaystyle c_{4}^{e}<\frac{\text{-}1}{2}$ $\displaystyle:$
$\displaystyle\\!\alpha$
$\displaystyle=\arcsin\left(\left|c_{4}^{e}\right|\right)-\frac{\pi}{6},$
$\displaystyle\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{A}^{zx},\;\boldsymbol{m}=\boldsymbol{e}_{y},\;\boldsymbol{l}=\boldsymbol{e}_{B}^{zx},$
$\displaystyle c_{4}^{e}>\frac{1}{2}$ $\displaystyle:$
$\displaystyle\\!\alpha$
$\displaystyle=\frac{\pi}{2}-\arcsin\left(\left|c_{4}^{e}\right|\right),$
$\displaystyle\\!\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{y},\;\boldsymbol{m}=\boldsymbol{e}_{B}^{zx},\;\boldsymbol{l}=\boldsymbol{e}_{A}^{zx}.$
(196)

Note that the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representations
(195) and (196) are very similar to the in-plane isotropic
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation, Eq. (118). In
fact, the latter can be recovered from Eq. (191) by setting $k=0$ as
$c_{0}^{e}=d_{1}/\left|\boldsymbol{d}\right|$, and the three cases support a
fine-tuned uniaxial nematic order parameter for
$\left|c_{4}^{e}\right|=\left|c_{2}^{e}\right|=\left|c_{0}^{e}\right|=\frac{1}{2}$.
In Figs. 6(e)-(f), we show the Fermi surface distortions corresponding,
respectively, to the condensation of the nematic order parameters
$d_{3}=d_{2yz}$ and $d_{4}=d_{2xz}$. In the former, one of the nematic in-
plane axis always aligns with $\boldsymbol{e}_{x}$ while for $d_{4}=d_{2xz}$,
it aligns with $\boldsymbol{e}_{y}$, in accordance with Eqs. (195) and (196).
In both cases, the twofold degeneracy is associated with the symmetry element
$C_{2z}$, i.e. a $180^{\circ}$ rotation about the $k_{z}$-axis.

We now analyze the case of the monoclinic crystal system, encompassed by the
$3$ axial point groups {$\mathsf{C_{2h}}$, $\mathsf{C_{1h}}=\mathsf{C_{s}}$,
$\mathsf{C_{2}}$}. In all cases, the decomposition of the five nematic
components $\boldsymbol{d}$ (17) is

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{monoclinic}}{\longrightarrow}\quad\Big{(}\,\underline{d_{1},\,d_{2},\,d_{5}},\,\big{\\{}d_{3},\,d_{4}\big{\\}}\,\Big{)}^{T}.$
(197)

As indicated by the underline, the three components
$d_{1}=d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}$, $d_{2}=d_{x^{2}-y^{2}}$,
and $d_{5}=d_{2xy}$ transform as the trivial IR and are thus always non-zero.
On the other hand, the nematic order parameters $d_{3}=d_{2yz}$ and
$d_{4}=d_{2xz}$ transform according to the same one-dimensional IR, as
indicated by the curly brackets. Consequently, the analysis of the Landau
expansion is analogous to that performed in Appendix G for the tetragonal
groups {$\mathsf{C_{4h}}$, $\mathsf{C_{4}}$, $\mathsf{S_{4}}$}. Indeed, we
obtain the same form of the nematic action
$\mathcal{S}=\mathcal{S}_{2}+\mathcal{S}_{4}$:

$\displaystyle\mathcal{S}_{2}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{1}\left(d_{3}\right)^{2}+r_{2}\left(d_{4}\right)^{2}+r_{3}d_{3}d_{4}\Big{\\}},$
(198) $\displaystyle\mathcal{S}_{4}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}u_{1}\left(d_{3}\right)^{4}+u_{2}\left(d_{4}\right)^{4}+u_{3}\left(d_{3}\right)^{2}\left(d_{4}\right)^{2}$
$\displaystyle+u_{4}d_{3}\left(d_{4}\right)^{3}+u_{5}\left(d_{3}\right)^{3}d_{4}\Big{\\}}.$
(199)

The outcome, as explained in Appendix G, is that a linear combination of
$d_{3}$ and $d_{4}$, which is enforced by the Landau coefficients, undergoes a
$Z_{2}$-Ising transition. The visualization of this nematic order parameter in
the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation offers
little insight, as the nematic axes $\boldsymbol{n}$, $\boldsymbol{m}$,
$\boldsymbol{l}$ can point anywhere in space.

The last crystal system is the triclinic one, described by the $2$ axial point
groups {$\mathsf{C_{1}}$, $\mathsf{C_{i}}=\mathsf{S_{2}}$}. The symmetry of
these crystals is so low that all five nematic components $\boldsymbol{d}$
transform as the trivial IR, which is a direct consequence of the absence of
rotational symmetry axes:

$\displaystyle\boldsymbol{d}$
$\displaystyle\quad\overset{\mathrm{triclinic}}{\longrightarrow}\quad\Big{(}\,\underline{d_{1},\,d_{2},\,d_{3},\,d_{4},\,d_{5}}\Big{)}^{T}.$
(200)

Therefore, nematic phase transitions cannot occur in triclinic crystals.

##  V Electronic nematicity in quasicrystalline axial point groups

While a large number of known quasicrystals is described by the icosahedral
(i.e., polyhedral) point groups {$\mathsf{I_{h}}$, $\mathsf{I}$}, whose
electronic nematic properties were analyzed in Sec. III.1, there are also
quasicrystalline materials with eightfold, tenfold, and twelvefold symmetry
[86]. In contrast to the icosahedral quasicrystals, the latter are
quasiperiodic in two directions and periodic along an axial direction;
consequently, they are described by non-crystallographic axial point groups
[88, 89, 105, 91]. In this subsection, we investigate the properties of the
electronic nematic order parameter in the class of octagonal, decagonal, and
dodecagonal point groups. Importantly, these point groups describe not only
quasicrystalline materials, but also artificial quasicrystals obtained from
twisting two crystalline 2D materials. The latter will be discussed in more
depth in Subsection V.4.

In all cases studied here, the symmetry-decomposition of the five nematic
components of $\boldsymbol{d}$ (17) have the same form as that for 2D
isotropic systems:

$\boldsymbol{d}\quad\overset{\mathrm{axial\>\>quasicrystals}}{\longrightarrow}\quad\Big{(}\,\underline{d_{1}},\,\boldsymbol{d}^{\mathrm{ip}},\,\boldsymbol{d}^{\mathrm{op}}\,\Big{)}^{T}.$
(201)

Following the discussion in Sec. IV.1 for 2D isotropic systems, our goal is to
obtain the first invariant in the Landau expansion of
$\boldsymbol{d}^{\alpha}$ that is not isotropic, i.e. that is not of the form
$\left|\boldsymbol{d}^{\alpha}\right|^{2n}$ with integer $n$ (here, $\alpha$
refers to either $\mathrm{ip}$ or $\mathrm{op}$). To accomplish this in a
systematic way, we evaluate the symmetrized decomposition of the
$N^{\mathrm{th}}$-order product
$\big{[}\otimes_{j=1}^{N}E_{\Gamma}\big{]}_{s}$, where $E_{\Gamma}$ is the
two-dimensional IR according to which $\boldsymbol{d}^{\alpha}$ transforms.
The leading-order anisotropic term is obtained from the product with the
smallest $N$ for which the decomposition gives either two invariants (if $N$
is even) or one invariant (if $N$ is odd). We note that Ref. [106] performed a
related analysis for the $d_{2xy}$ and $d_{x^{2}-y^{2}}$ superconducting order
parameters in $\mathsf{D_{n}}$ point groups.

A summary of all the results presented in this section, as well as in Sec.
III.1 for the icosahedral quasicrystal, is contained in Table 2.

###  V.1 Dodecagonal quasicrystals

The non-crystallographic point groups that possess dodecagonal symmetry are
{$\mathsf{D_{12h}}$, $\mathsf{D_{12}}$, $\mathsf{C_{12v}}$, $\mathsf{D_{6d}}$,
$\mathsf{C_{12h}}$, $\mathsf{C_{12}}$, $\mathsf{S_{12}}$}. For concreteness,
hereafter we consider $\mathsf{D_{12h}}$, for which
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ transform as
the $E_{2g}$ and $E_{1g}$ IR, respectively. We will discuss later how the
results generalize to the other dodecagonal groups.

To derive the Landau expansion of the doublets, we follow the procedure
outlined in the beginning of this section. Focusing first on the in-plane
doublet $\boldsymbol{d}^{\mathrm{ip}}$, we compute the symmetrized-product
decomposition of the corresponding IR:

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{4g},$ (202)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{2g}\big{]}_{s}$
$\displaystyle=B_{1g}\oplus B_{2g}\oplus E_{2g},$ (203)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\\!\oplus 2E_{4g},$ (204)
$\displaystyle\big{[}\otimes_{j=1}^{5}E_{2g}\big{]}_{s}$
$\displaystyle=B_{1g}\oplus B_{2g}\oplus 2E_{2g},$ (205)
$\displaystyle\big{[}\otimes_{j=1}^{6}E_{2g}\big{]}_{s}$
$\displaystyle=2A_{1g}\\!\oplus A_{2g}\oplus 2E_{4g}.$ (206)

Therefore, since there are two $A_{1g}$ invariants in (206), whereas all other
even orders have only one invariant, we conclude that the leading order
anisotropic term in the Landau expansion of $\boldsymbol{d}^{\mathrm{ip}}$
appears at sixth-order. It can be constructed from the product of the $B_{1g}$
and $B_{2g}$ trilinears,

$\displaystyle B_{1g}:$
$\displaystyle\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{3}\cos\left(3\gamma_{\mathrm{ip}}\right),$
$\displaystyle B_{2g}:$
$\displaystyle\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{3}\sin\left(3\gamma_{\mathrm{ip}}\right),$
(207)

and thus, it becomes
$\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{6}\cos\left(6\gamma_{\mathrm{ip}}\right)$,
resulting in the action:

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

The Landau expansion (208) is equivalent to that of the $Z_{6}$-clock model,
displaying a sixfold-degenerate ground state

$\displaystyle\gamma_{\mathrm{ip}}^{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\\}.$ (209)

As for the out-of-plane doublet $\boldsymbol{d}^{\mathrm{op}}$, which
transforms as the $E_{1g}$ IR, the symmetrized-product decomposition up to
sixth-order yields only isotropic invariants of the form
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2n}$:

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{2g},$ (210)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{1g}\big{]}_{s}$
$\displaystyle=E_{1g}\oplus E_{3g},$ (211)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{2g}\oplus E_{4g},$ (212)
$\displaystyle\big{[}\otimes_{j=1}^{5}E_{1g}\big{]}_{s}$
$\displaystyle=E_{1g}\oplus E_{3g}\oplus E_{5g},$ (213)
$\displaystyle\big{[}\otimes_{j=1}^{6}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus B_{1g}\oplus B_{2g}\oplus E_{2g}\oplus E_{4g}.$
(214)

It turns out that an anisotropic invariant can only be constructed at twelfth-
order, from the product of the $B_{1g}$ and $B_{2g}$ IRs that appear in the
sixth-order product decomposition:

$\displaystyle B_{1g}:$
$\displaystyle\left|\boldsymbol{d}^{\mathrm{op}}\right|^{6}\cos\left(6\gamma_{\mathrm{op}}\right),$
$\displaystyle B_{2g}:$
$\displaystyle\left|\boldsymbol{d}^{\mathrm{op}}\right|^{6}\sin\left(6\gamma_{\mathrm{op}}\right).$
(215)

We thus obtain the nematic action:

$\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}+v_{12}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{12}\cos\left(12\gamma_{\mathrm{op}}\right)\Big{\\}},$
(216)

which has the same form as the Landau expansion of the $Z_{12}$-clock model.
The ground-state is parametrized by

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

The in-plane and out-of-plane nematic states given by (209),(217) can be
expressed in the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation
in a straightforward way via the relationships (118) and (119). We show the
corresponding Fermi surface distortions in Figs. 7(a),(b). For
$\boldsymbol{d}^{\mathrm{ip}}$, both in-plane nematic axes point along the
high-symmetry directions associated with twofold in-plane rotations, whereas
the third axis is parallel to the $k_{z}$-axis. In contrast, for
$\boldsymbol{d}^{\mathrm{op}}$, only the nematic axis $\boldsymbol{l}$ aligns
with a high-symmetry direction.

The same results obtained for $\mathsf{D_{12h}}$ also hold for the dodecagonal
groups {$\mathsf{D_{12}}$, $\mathsf{C_{12v}}$, $\mathsf{D_{6d}}$}. On the
other hand, for {$\mathsf{C_{12h}}$, $\mathsf{C_{12}}$, $\mathsf{S_{12}}$},
due to the lack of in-plane symmetry axes, the nematic doublets transform as
complex IRs, resulting in an offset angle in the clock term of the type (125),
similarly to what is shown in Appendices E-G for crystalline axial point
groups.

###  V.2 Decagonal quasicrystals

The are seven non-crystallographic point groups that possess tenfold symmetry,
{$\mathsf{D_{10h}}$, $\mathsf{D_{10}}$, $\mathsf{D_{5h}}$, $\mathsf{C_{10v}}$,
$\mathsf{C_{10h}}$, $\mathsf{C_{5h}}$, $\mathsf{C_{10}}$}. Similarly to the
previous subsection, we need to distinguish the groups {$\mathsf{D_{10h}}$,
$\mathsf{D_{10}}$, $\mathsf{D_{5h}}$, $\mathsf{C_{10v}}$} from the groups
{$\mathsf{C_{10h}}$, $\mathsf{C_{5h}}$, $\mathsf{C_{10}}$}, as the latter do
not have in-plane symmetry axes. The only effect of this lack of in-plane
rotational symmetry is that the clock term of the nematic action acquires an
offset angle, see Eq. (125).

We consider here the case of $\mathsf{D_{10h}}$. Starting with the in-plane
doublet $\boldsymbol{d}^{\mathrm{ip}}$, which transforms as the $E_{2g}$ IR,
we obtain the symmetrized-product decomposition:

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{4g},$ (218)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{2g}\big{]}_{s}$
$\displaystyle=E_{2g}\oplus E_{4g},$ (219)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\\!\oplus E_{2g}\oplus E_{4g},$ (220)
$\displaystyle\big{[}\otimes_{j=1}^{5}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\\!\oplus A_{2g}\oplus E_{2g}\oplus E_{4g}.$ (221)

According to what was explained in the introduction of this section, the
fifth-order invariant must be the leading-order term in the Landau expansion
that is not proportional to the square of the doublet. Directly employing a
fifth-order transformation condition, as an extension to Eq. (266), one finds
this invariant to be
$\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{5}\cos\left(5\gamma_{\mathrm{ip}}\right)$,
and thus the nematic action becomes:

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

which maps onto the $Z_{5}$-clock model. Minimization of the clock term gives
with the fivefold-degenerate ground state

$\displaystyle\gamma_{\mathrm{ip}}^{0}$
$\displaystyle=\frac{\pi}{5}\left(\frac{1+\mathrm{sign}v_{5}}{2}\right)+\frac{2\pi}{5}n,$
$\displaystyle n$ $\displaystyle\in\\{0,1,2,3,4\\}.$ (223)

Moving on to the out-of-plane doublet $\boldsymbol{d}^{\mathrm{op}}$, we note
that it transforms as the $E_{1g}$ IR, whose symmetrized product decomposition
is:

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{2g},$ (224)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{1g}\big{]}_{s}$
$\displaystyle=E_{1g}\oplus E_{3g},$ (225)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{2g}\oplus E_{4g},$ (226)
$\displaystyle\big{[}\otimes_{j=1}^{5}E_{1g}\big{]}_{s}$
$\displaystyle=B_{1g}\oplus B_{2g}\oplus E_{1g}\oplus E_{3g},$ (227)
$\displaystyle\big{[}\otimes_{j=1}^{6}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{2g}\oplus 2E_{4g}.$ (228)

Similarly to the case of the dodecagonal point group $\mathsf{D_{12h}}$, only
powers of $\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}$ occur up to sixth-
order. Following the same steps as in that case, we obtain the clock term from
the product of the fifth-order product-decomposition elements $B_{1g}$ and
$B_{2g}$:

$\displaystyle B_{1g}:$
$\displaystyle\left|\boldsymbol{d}^{\mathrm{op}}\right|^{5}\cos\left(5\gamma_{\mathrm{op}}\right),$
$\displaystyle B_{2g}:$
$\displaystyle\left|\boldsymbol{d}^{\mathrm{op}}\right|^{5}\sin\left(5\gamma_{\mathrm{op}}\right),$

yielding:

$\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}+v_{10}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{10}\cos\left(10\gamma_{\mathrm{op}}\right)\Big{\\}}.$
(229)

We therefore find $\boldsymbol{d}^{\mathrm{op}}$ to be a $Z_{10}$-clock
nematic order parameter, whose tenfold degenerate ground states are given by:

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

In Figs. 7(c)-(d), we show the distorted Fermi surfaces obtained from
expressing Eqs. (223),(230) in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation. The Fermi
surface properties are analogous to the dodecagonal case, except that the in-
plane high-symmetry directions are those associated with tenfold rotational
symmetry.

Figure 7: 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
dodecagonal in-plane and out-of-plane nematicity (209),(217) [panels (a)-(b)],
decagonal in-plane and out-of-plane nematicity (223),(230) [panels (c)-(d)],
and octagonal in-plane and out-of-plane nematicity (236),(242) [panels
(e)-(f)]. The unit vectors
$\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}$ are rescaled
by
$2.25\,k_{F}\left(\left\\{\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}\right\\}\right)$,
respectively. We chose the trivial component such that
$d_{1}/\left|\boldsymbol{d}\right|=0.1$.

###  V.3 Octagonal quasicrystals

The analysis of the octagonal point groups {$\mathsf{D_{8h}}$,
$\mathsf{D_{8}}$, $\mathsf{C_{8v}}$, $\mathsf{D_{4d}}$, $\mathsf{C_{8h}}$,
$\mathsf{C_{8}}$, $\mathsf{S_{8}}$} mirrors the analyses of the previous two
subsections for the dodecagonal and decagonal quasicrystals. Considering the
point group $\mathsf{D_{8h}}$ for concreteness, we note that
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ transform as
$E_{2g}$ and $E_{1g}$, respectively. Performing the symmetrized product
decomposition for $E_{2g}$ gives:

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{2g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus B_{1g}\oplus B_{2g},$ (231)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{2g}\big{]}_{s}$
$\displaystyle=2E_{2g},$ (232)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{2g}\big{]}_{s}$
$\displaystyle=2A_{1g}\\!\oplus A_{2g}\oplus B_{1g}\oplus B_{2g}.$ (233)

Using the bilinears

$\displaystyle D^{B_{1g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\cos\left(2\gamma_{\mathrm{ip}}\right),$
$\displaystyle D^{B_{2g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\sin\left(2\gamma_{\mathrm{ip}}\right),$
(234)

it is straightforward to construct the three invariants up to quartic order,
$\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}$,
$\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{4}$ and
$(D^{B_{1g}})^{2}-(D^{B_{2g}})^{2}$. We therefore obtain the Landau expansion
of the $Z_{4}$-clock model:

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

which has the fourfold-degenerate ground state:

$\displaystyle\gamma_{\mathrm{ip}}^{0}$
$\displaystyle=\frac{\pi}{4}\left(\frac{1+\mathrm{sign}v_{4}}{2}\right)+\frac{2\pi}{4}n,$
$\displaystyle n$ $\displaystyle\in\\{0,1,2,3\\},$ (236)

see also Eq. (188).

Considering now the out-of-plane doublet $\boldsymbol{d}^{\mathrm{op}}$, the
symmetrized product decomposition of $E_{1g}$ to sixth-order only gives
invariants that are powers of $\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}$:

$\displaystyle\big{[}\otimes_{j=1}^{2}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus E_{2g},$ (237)
$\displaystyle\big{[}\otimes_{j=1}^{3}E_{1g}\big{]}_{s}$
$\displaystyle=E_{1g}\oplus E_{3g},$ (238)
$\displaystyle\big{[}\otimes_{j=1}^{4}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus B_{1g}\oplus B_{2g}\oplus E_{2g},$ (239)
$\displaystyle\big{[}\otimes_{j=1}^{5}E_{1g}\big{]}_{s}$
$\displaystyle=E_{1g}\oplus 2E_{3g},$ (240)
$\displaystyle\big{[}\otimes_{j=1}^{6}E_{1g}\big{]}_{s}$
$\displaystyle=A_{1g}\oplus B_{1g}\oplus B_{2g}\oplus 2E_{2g}.$ (241)

Analogously to the case of the dodecagonal and decagonal quasicrystals, it is
straightforward to find the leading-order anisotropic term, this time from the
product of the quartic combinations:

$\displaystyle B_{1g}:$
$\displaystyle\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}\cos\left(4\gamma_{\mathrm{op}}\right),$
$\displaystyle B_{2g}:$
$\displaystyle\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}\sin\left(4\gamma_{\mathrm{op}}\right).$

Hence, we find the nematic action

$\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}+v_{8}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{8}\cos\left(8\gamma_{\mathrm{op}}\right)\Big{\\}},$
(242)

corresponding to the $Z_{8}$-clock model. The directions of
$\boldsymbol{d}^{\mathrm{op}}$ that minimize the action are eightfold
degenerate:

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

The same results hold for {$\mathsf{D_{8}}$, $\mathsf{C_{8v}}$,
$\mathsf{D_{4d}}$}, whereas in the groups {$\mathsf{C_{8h}}$,
$\mathsf{C_{8}}$, $\mathsf{S_{8}}$}, the clock terms in the nematic actions of
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{d}^{\mathrm{op}}$ acquire an
offset angle of the form (125).

The $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation of the in-
plane and out-of-plane nematic doublets (236) and (243) is represented in
Figs. 7(e)-(f) via the corresponding distortions of the Fermi surface. The
high-symmetry in-plane directions associated with the eightfold symmetry of
the quasicrystal are highlighted in the figure.

###  V.4 Twisted quasicrystals

The rapid advances in the field of twistronics opened a new path to
investigate the properties of emergent quasicrystals created by twisting 2D
materials. In some cases, the twisted structure is an incommensurate lattice
that, however, retains a crystallographic point group. This is the case of
twisted bilayer graphene with a non-commensurate twist angle, although the
effects of quasi-periodicity are dramatically enhanced when three graphene
layers are twisted by two different angles, as shown recently in Ref. [85].
Interestingly, this work reported the emergence of robust superconductivity in
this quasi-periodic structure. In what concerns electronic nematicity, as long
as the point group is the same as a crystallographic one, the results derived
in Sec. IV can be directly applied.

An alternative route to realize twisted quasicrystals, as shown in Refs. [94,
95], is to start with two identical 2D materials described by the (planar)
crystallographic point group $\mathsf{D_{n}}$ (with $n=2,3,4,\,6$). Upon
twisting the two crystals by a relative angle of $\pi/n$ with respect to their
common $z$-axis, the resulting structure becomes invariant under an $n$-fold
rotation about $z$ followed by a horizontal mirror reflection (which maps one
layer onto the other). This is nothing but the improper rotation element
$S_{2nz}$. Comparing Eqs. (92) and (95), we conclude that the twisted bilayer
is actually a quasi-periodic “lattice” with point group $\mathsf{D_{nd}}$,
i.e. $\mathsf{D_{n}}\overset{\mathrm{twist}}{\longrightarrow}\mathsf{D_{nd}}$.
Experimentally, both twisted bilayer graphene with a twist angle of
$30^{\circ}$ and twisted bilayer cuprate with a twist angle of $45^{\circ}$
have been realized [92, 96].

While much of the recent interest in these constructions have focused on
realizing time-reversal symmetry-breaking superconductivity [93, 94, 95], it
has been recently pointed out that electronic nematicity can be strongly
impacted by the enhanced symmetry of the twisted bilayer [59]. Here, we
further explore this idea by discussing the properties of the in-plane nematic
order parameter $\boldsymbol{d}^{\mathrm{ip}}$ in the four cases of
$\pi/n$-twisted bilayers ($n=2,3,4,\,6$).

The point group in the cases $n=2$ and $n=3$ remain crystallographic,
$\mathsf{D_{2}\overset{\mathrm{twist}}{\longrightarrow}\mathsf{D_{2d}}}$ and
$\mathsf{D_{3}\overset{\mathrm{twist}}{\longrightarrow}\mathsf{D_{3d}}}$.
There is, however, one important difference. While $\mathsf{D_{3}}$ and
$\mathsf{D_{3d}}$ both belong to the trigonal crystal system, $\mathsf{D_{2}}$
and $\mathsf{D_{2d}}$ belong to different systems – monoclinic and
orthorhombic, respectively. In terms of the components of
$\boldsymbol{d}^{\mathrm{ip}}=\left(d_{2},\,d_{5}\right)^{T}$, this means that
while in $\mathsf{D_{2}}$ only the shear component $d_{5}$ transforms non-
trivially, in $\mathsf{D_{2d}}$ both components transform non-trivially as two
different one-dimensional IRs – see Table 1. Therefore, twisting two
monoclinic layers by $90^{\circ}$ opens up a new nematic instability channel.

For two tetragonal layers twisted by $45^{\circ}$ ($n=4$), the twisted bilayer
has a non-crystallographic octagonal point group,
$\mathsf{D_{4}\overset{\mathrm{twist}}{\longrightarrow}\mathsf{D_{4d}}}$ [93].
Using the results above, we see that the two components of the in-plane
nematic doublet $\boldsymbol{d}^{\mathrm{ip}}$ change from transforming as the
one-dimensional IRs $B_{1}$ and $B_{2}$ of $\mathsf{D_{4}}$ in each layer to
transforming as a single two-dimensional IR $E_{2}$ of $\mathsf{D_{4d}}$ in
the twisted bilayer. Therefore, the character of the nematic transition
changes from $Z_{2}$-Ising for the individual tetragonal layers to
$Z_{4}$-clock for the coupled $45^{\circ}$-twisted bilayer. As discussed in
the end of Sec. IV.1, the 2D $4$-state clock model (which has the same
properties as the 2D Ashkin-Teller model) undergoes a transition from the
disordered to the ordered phase that can be understood as the merging of the
two BKT transitions of the $q$-state clock model as $q\rightarrow 4^{+}$.
Because of its unique character, the exponents of this transition are non-
universal and depend on $g_{\mathrm{ip}}$, except for the anomalous exponent
$\eta=1/4$. Since electronic nematicity is observed in tetragonal cuprates and
iron pnictides/chalcogenides [43, 44], this setting offers an interesting path
to realize Ashkin-Teller nematicity.

Finally, the case of two hexagonal layers twisted by $30^{\circ}$ ($n=6$)
corresponds to
$\mathsf{D_{6}\overset{\mathrm{twist}}{\longrightarrow}\mathsf{D_{6d}}}$ and
was investigated in Ref. [59]. The outcome, which follows from the nematic
properties of dodecagonal quasicrystals, is that while the uncoupled layers
undergo a second-order nematic transition in the 2D $3$-state clock/Potts
universality class (see Table 1), the coupled $30^{\circ}$-twisted bilayer
undergoes two BKT transitions: the higher one from the disordered to the
critical phase (where there is only quasi-long-range nematic order) and the
lower one from the critical phase to the long-range ordered nematic phase.

cubic | IR | nem. basis | universality | res. PG | transf. matrix
---|---|---|---|---|---
$\mathsf{O_{h}}$ | $E_{g}$ | $\boldsymbol{d}^{e}$ |  $Z_{3}$-Potts | $\mathsf{D_{4h}}$ |
| $T_{2g}$ | $\boldsymbol{d}^{t}$ |  $Z_{4}$-Potts | $\mathsf{D_{3d}}$ | $\mathcal{R}^{A_{2u}}_{g}\mathcal{R}^{(3)}_{g}$
$\mathsf{O}$ | $E$ | $\boldsymbol{d}^{e}$ |  $Z_{3}$-Potts | $\mathsf{D_{4}}$ |
| $T_{2}$ | $\boldsymbol{d}^{t}$ |  $Z_{4}$-Potts | $\mathsf{D_{3}}$ | $\mathcal{R}^{A_{2}}_{g}\mathcal{R}^{(3)}_{g}$
$\mathsf{T_{d}}$ | $E$ | $\boldsymbol{d}^{e}$ |  $Z_{3}$-Potts | $\mathsf{D_{2d}}$ |
| $T_{2}$ | $\boldsymbol{d}^{t}$ |  $Z_{4}$-Potts | $\mathsf{C_{3v}}$ | $\mathcal{R}^{(3)}_{g}$
$\mathsf{T_{h}}$ | $\bar{E}_{g}$ | $\boldsymbol{d}^{e}$ |  $Z_{3}^{*}$-Potts | $\mathsf{D_{2h}}$ |
| $T_{g}$ | $\boldsymbol{d}^{t}$ |  $Z_{4}$-Potts | $\mathsf{S_{6}}$ | $\mathcal{R}^{A_{u}}_{g}\mathcal{R}^{(3)}_{g}$
$\mathsf{T}$ | $\bar{E}$ | $\boldsymbol{d}^{e}$ |  $Z_{3}^{*}$-Potts | $\mathsf{D_{2}}$ |
| $T$ | $\boldsymbol{d}^{t}$ |  $Z_{4}$-Potts | $\mathsf{C_{3}}$ | $\mathcal{R}^{(3)}_{g}$
hexagonal | IR | nem. basis | universality | res. PG | transf. matrix
$\mathsf{D_{6h}}$ | $E_{2g}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{3}$-clock | $\mathsf{D_{2h}}$ | $\mathcal{R}_{g}^{B_{1u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
| $E_{1g}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{6}$-clock | $\mathsf{C_{2h}}$ | $\mathcal{R}_{g}^{A_{1u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{D_{6}}$ | $E_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{3}$-clock | $\mathsf{D_{2}}$ | $\mathcal{R}_{g}^{B_{1}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
| $E_{1}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{6}$-clock | $\mathsf{C_{2}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{D_{3h}}$ | $E^{\prime}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{3}$-clock | $\mathsf{C_{2v}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
| $E^{\prime\prime}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{6}$-clock | $\mathsf{C_{2}}$ | $\mathcal{R}_{g}^{A_{1}^{\prime\prime}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{6v}}$ | $E_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{3}$-clock | $\mathsf{C_{2v}}$ | $\mathcal{R}_{g}^{B_{2}}\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$
| $E_{1}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{6}$-clock | $\mathsf{C_{s}}$ | $\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$
$\mathsf{C_{6h}}$ | $\bar{E}_{2g}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{3}^{*}$-clock | $\mathsf{C_{2h}}$ | $\mathcal{R}_{g}^{B_{u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
| $\bar{E}_{1g}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{6}^{*}$-clock | $\mathsf{C_{i}}$ | $\mathcal{R}_{g}^{A_{u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{3h}}$ | $\bar{E}^{\prime}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{3}^{*}$-clock | $\mathsf{C_{s}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
| $\bar{E}^{\prime\prime}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{6}^{*}$-clock | $\mathsf{C_{1}}$ | $\mathcal{R}_{g}^{A^{\prime\prime}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{6}}$ | $\bar{E}_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{3}^{*}$-clock | $\mathsf{C_{2}}$ | $\mathcal{R}_{g}^{B}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
| $\bar{E}_{1}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{6}^{*}$-clock | $\mathsf{C_{1}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
trigonal | IR | nem. basis | universality | res. PG | transf. matrix
$\mathsf{D_{3d}}$ | $E_{g}$ | $\left\\{\boldsymbol{d}^{\mathrm{ip}},\boldsymbol{d}^{\mathrm{op}}\right\\}$ |  $Z_{3}$-clock | $\mathsf{C_{2h}}$ | $\mathcal{R}_{g}^{A_{1u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{D_{3}}$ | $E$ | $\left\\{\boldsymbol{d}^{\mathrm{ip}},\boldsymbol{d}^{\mathrm{op}}\right\\}$ |  $Z_{3}$-clock | $\mathsf{C_{2}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{3v}}$ | $E$ | $\left\\{\boldsymbol{d}^{\mathrm{ip}},\boldsymbol{d}^{\mathrm{op}}\right\\}$ |  $Z_{3}$-clock | $\mathsf{C_{s}}$ | $\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$
$\mathsf{S_{6}}$ | $\bar{E}_{g}$ | $\left\\{\boldsymbol{d}^{\mathrm{ip}},\boldsymbol{d}^{\mathrm{op}}\right\\}$ |  $Z_{3}^{*}$-clock | $\mathsf{C_{i}}$ | $\mathcal{R}_{g}^{A_{u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{3}}$ | $\bar{E}$ | $\left\\{\boldsymbol{d}^{\mathrm{ip}},\boldsymbol{d}^{\mathrm{op}}\right\\}$ |  $Z_{3}^{*}$-clock | $\mathsf{C_{1}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$

tetragonal IR nem. basis universality res. PG transf. matrix $\mathsf{D_{4h}}$
$B_{1g}$ $d_{x^{2}\text{-}y^{2}}$ $Z_{2}$-Ising $\mathsf{D_{2h}}$ $B_{2g}$
$d_{2xy}$ $Z_{2}$-Ising $\mathsf{D_{2h}}$ $E_{g}$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{4}$-clock $\mathsf{C_{2h}}$
$\mathcal{R}_{g}^{A_{1u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{D_{4}}$ $B_{1}$ $d_{x^{2}\text{-}y^{2}}$ $Z_{2}$-Ising
$\mathsf{D_{2}}$ $B_{2}$ $d_{2xy}$ $Z_{2}$-Ising $\mathsf{D_{2}}$ $E$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{4}$-clock $\mathsf{C_{2}}$
$\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$ $\mathsf{C_{4v}}$ $B_{1}$
$d_{x^{2}\text{-}y^{2}}$ $Z_{2}$-Ising $\mathsf{C_{2v}}$ $B_{2}$ $d_{2xy}$
$Z_{2}$-Ising $\mathsf{C_{2v}}$ $E$ $\boldsymbol{d}^{\mathrm{op}}$
$Z_{4}$-clock $\mathsf{C_{s}}$ $\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$
$\mathsf{D_{2d}}$ $B_{1}$ $d_{x^{2}\text{-}y^{2}}$ $Z_{2}$-Ising
$\mathsf{D_{2}}$ $B_{2}$ $d_{2xy}$ $Z_{2}$-Ising $\mathsf{C_{2v}}$ $E$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{4}$-clock $C_{2}/\mathsf{C_{s}}$
$\mathcal{R}_{g}^{(2)}$ $\mathsf{C_{4h}}$ $B_{g}$
$\big{\\{}d_{x^{2}\text{-}y^{2}},d_{2xy}\big{\\}}$ $Z^{*}_{2}$-Ising
$\mathsf{C_{2h}}$ $\bar{E}_{g}$ $\boldsymbol{d}^{\mathrm{op}}$
$Z_{4}^{*}$-clock $\mathsf{C_{i}}$
$\mathcal{R}_{g}^{A_{u}}\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$
$\mathsf{C_{4}}$ $B$ $\big{\\{}d_{x^{2}\text{-}y^{2}},d_{2xy}\big{\\}}$
$Z^{*}_{2}$-Ising $\mathsf{C_{2}}$ $\bar{E}$ $\boldsymbol{d}^{\mathrm{op}}$
$Z_{4}^{*}$-clock $\mathsf{C_{1}}$ $\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$
$\mathsf{S_{4}}$ $B$ $\big{\\{}d_{x^{2}\text{-}y^{2}},d_{2xy}\big{\\}}$
$Z^{*}_{2}$-Ising $\mathsf{C_{2}}$ $\bar{E}$ $\boldsymbol{d}^{\mathrm{op}}$
$Z_{4}^{*}$-clock $\mathsf{C_{1}}$ $\mathcal{R}_{g}^{(2)}$ orthorombic IR nem.
basis universality res. PG transf. matrix $\mathsf{D_{2h}}$ $B_{1g}$ $d_{2xy}$
$Z_{2}$-Ising $\mathsf{C_{2h}}$ $B_{2g}$ $d_{2xz}$ $Z_{2}$-Ising
$\mathsf{C_{2h}}$ $B_{3g}$ $d_{2yz}$ $Z_{2}$-Ising $\mathsf{C_{2h}}$
$\mathsf{D_{2}}$ $B_{1}$ $d_{2xy}$ $Z_{2}$-Ising $\mathsf{C_{2}}$ $B_{2}$
$d_{2xz}$ $Z_{2}$-Ising $\mathsf{C_{2}}$ $B_{3}$ $d_{2yz}$ $Z_{2}$-Ising
$\mathsf{C_{2}}$ $\mathsf{C_{2v}}$ $A_{2}$ $d_{2xy}$ $Z_{2}$-Ising
$\mathsf{C_{2}}$ $B_{1}$ $d_{2xz}$ $Z_{2}$-Ising $\mathsf{C_{s}}$ $B_{2}$
$d_{2yz}$ $Z_{2}$-Ising $\mathsf{C_{s}}$ monoclinic IR nem. basis universality
res. PG transf. matrix $\mathsf{C_{2h}}$ $B_{g}$
$\left\\{d_{2xz},d_{2yz}\right\\}$ $Z_{2}$-Ising $\mathsf{C_{i}}$
$\mathsf{C_{1h}}\equiv\mathsf{C_{s}}$ $A^{\prime\prime}$
$\left\\{d_{2xz},d_{2yz}\right\\}$ $Z_{2}$-Ising $\mathsf{C_{1}}$
$\mathsf{C_{2}}$ $B$ $\left\\{d_{2xz},d_{2yz}\right\\}$ $Z_{2}$-Ising
$\mathsf{C_{1}}$

Table 1: Symmetry properties of the allowed nematic order parameters in all
crystallographic point groups (PG). The columns list the corresponding
irreducible representation (IR), the relevant nematic basis vector (see main
text for the definition), the universality class of the corresponding Landau
expansion, the residual PG, and the transformation matrices (see main text).
Note that the $Z_{3}$-clock and $Z_{3}$-Potts models have the same properties.
The $Z_{q}^{*}$-clock model has the same form as the $Z_{q}$-clock model,
except that the clock term has a non-universal offset angle. The overbar over
the IR denotes a complex IR.

octagonal | IR | nem. basis | universality | res. PG | transf. matrix
---|---|---|---|---|---
$\mathsf{D_{8h}}$ | $E_{2g}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{4}$-clock | $\mathsf{D_{2h}}$ |
| $E_{1g}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{8}$-clock | $\mathsf{C_{2h}}$ | $\mathcal{R}_{g}^{A_{1u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{D_{8}}$ | $E_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{4}$-clock | $\mathsf{D_{2}}$ |
| $E_{1}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{8}$-clock | $\mathsf{C_{2}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{8v}}$ | $E_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{4}$-clock | $\mathsf{C_{2v}}$ |
| $E_{1}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{8}$-clock | $\mathsf{C_{s}}$ | $\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$
$\mathsf{D_{4d}}$ | $E_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{4}$-clock | $\mathsf{D_{2}}/\mathsf{C_{2v}}$ |
| $E_{3}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{8}$-clock | $\mathsf{C_{2}}/\mathsf{C_{s}}$ | $\mathcal{R}_{g}^{B_{1}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{8h}}$ | $\bar{E}_{2g}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{4}^{*}$-clock | $\mathsf{C_{2h}}$ |
| $\bar{E}_{1g}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{8}^{*}$-clock | $\mathsf{C_{i}}$ | $\mathcal{R}_{g}^{A_{u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{8}}$ | $\bar{E}_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{4}^{*}$-clock | $\mathsf{C_{2}}$ |
| $\bar{E}_{1}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{8}^{*}$-clock | $\mathsf{C_{1}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{S_{8}}$ | $\bar{E}_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{4}^{*}$-clock | $\mathsf{C_{2}}$ |
| $\bar{E}_{3}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{8}^{*}$-clock | $\mathsf{C_{1}}$ | $\mathcal{R}_{g}^{B}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
decagonal | IR | nem. basis | universality | res. PG | transf. matrix
$\mathsf{D_{10h}}$ | $E_{2g}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{5}$-clock | $\mathsf{D_{2h}}$ |
| $E_{1g}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{10}$-clock | $\mathsf{C_{2h}}$ | $\mathcal{R}_{g}^{A_{1u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{D_{10}}$ | $E_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{5}$-clock | $\mathsf{D_{2}}$ |
| $E_{1}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{10}$-clock | $\mathsf{C_{2}}$ | $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{D_{5h}}$ | $E_{2}^{\prime}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{5}$-clock | $\mathsf{C_{2v}}$ |
| $E_{1}^{\prime\prime}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{10}$-clock | $\mathsf{C_{2}}/\mathsf{C_{s}}$ | $\mathcal{R}_{g}^{A_{1}^{\prime\prime}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{10v}}$ | $E_{2}$ | $\boldsymbol{d}^{\mathrm{ip}}$ |  $Z_{5}$-clock | $\mathsf{C_{2v}}$ |
| $E_{1}$ | $\boldsymbol{d}^{\mathrm{op}}$ |  $Z_{10}$-clock | $\mathsf{C_{s}}$ | $\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$
|  |  |  |  |

$\mathsf{C_{10h}}$ $\bar{E}_{2g}$ $\boldsymbol{d}^{\mathrm{ip}}$
$Z_{5}^{*}$-clock $\mathsf{C_{2h}}$ $\bar{E}_{1g}$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{10}^{*}$-clock $\mathsf{C_{i}}$
$\mathcal{R}_{g}^{A_{u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{5h}}$ $\bar{E}_{2}^{\prime}$ $\boldsymbol{d}^{\mathrm{ip}}$
$Z_{5}^{*}$-clock $\mathsf{C_{s}}$ $\bar{E}_{1}^{\prime\prime}$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{10}^{*}$-clock $\mathsf{C_{1}}$
$\mathcal{R}_{g}^{A^{\prime\prime}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{10}}$ $\bar{E}_{2}$ $\boldsymbol{d}^{\mathrm{ip}}$
$Z_{5}^{*}$-clock $\mathsf{C_{2}}$ $\bar{E}_{1}$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{10}^{*}$-clock $\mathsf{C_{1}}$
$\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$ dodecagonal IR nem. basis
universality res. PG transf. matrix $\mathsf{D_{12h}}$ $E_{2g}$
$\boldsymbol{d}^{\mathrm{ip}}$ $Z_{6}$-clock $\mathsf{D_{2h}}$ $E_{1g}$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{12}$-clock $\mathsf{C_{2h}}$
$\mathcal{R}_{g}^{A_{1u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{D_{12}}$ $E_{2}$ $\boldsymbol{d}^{\mathrm{ip}}$ $Z_{6}$-clock
$\mathsf{D_{2}}$ $E_{1}$ $\boldsymbol{d}^{\mathrm{op}}$ $Z_{12}$-clock
$\mathsf{C_{2}}$ $\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{12v}}$ $E_{2}$ $\boldsymbol{d}^{\mathrm{ip}}$ $Z_{6}$-clock
$\mathsf{C_{2v}}$ $E_{1}$ $\boldsymbol{d}^{\mathrm{op}}$ $Z_{12}$-clock
$\mathsf{C_{s}}$ $\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$ $\mathsf{D_{6d}}$
$E_{2}$ $\boldsymbol{d}^{\mathrm{ip}}$ $Z_{6}$-clock
$\mathsf{D_{2}}/\mathsf{C_{2v}}$ $E_{5}$ $\boldsymbol{d}^{\mathrm{op}}$
$Z_{12}$-clock $\mathsf{C_{2}}/\mathsf{C_{s}}$
$\mathcal{R}_{g}^{B_{1}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{12h}}$ $\bar{E}_{2g}$ $\boldsymbol{d}^{\mathrm{ip}}$
$Z_{6}^{*}$-clock $\mathsf{C_{2h}}$ $\bar{E}_{1g}$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{12}^{*}$-clock $\mathsf{C_{i}}$
$\mathcal{R}_{g}^{A_{u}}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$
$\mathsf{C_{12}}$ $\bar{E}_{2}$ $\boldsymbol{d}^{\mathrm{ip}}$
$Z_{6}^{*}$-clock $\mathsf{C_{2}}$ $\bar{E}_{1}$
$\boldsymbol{d}^{\mathrm{op}}$ $Z_{12}^{*}$-clock $\mathsf{C_{1}}$
$\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$ $\mathsf{S_{12}}$ $\bar{E}_{2}$
$\boldsymbol{d}^{\mathrm{ip}}$ $Z_{6}^{*}$-clock $\mathsf{C_{2}}$
$\bar{E}_{5}$ $\boldsymbol{d}^{\mathrm{op}}$ $Z_{12}^{*}$-clock
$\mathsf{C_{1}}$
$\mathcal{R}_{g}^{B}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$ icosahedral IR
nem. basis universality res. PG transf. matrix $\mathsf{I_{h}}$ $H_{g}$
$\boldsymbol{d}$ $\mathsf{D_{3d}}/\mathsf{D_{5d}}$ $\mathsf{I}$ $H$
$\boldsymbol{d}$ $\mathsf{D_{3}}/\mathsf{D_{5}}$

Table 2: Symmetry properties of the allowed nematic order parameters in non-
crystallographic point groups (PG) describing quasicrystals. The columns
represent the same quantities as those in Table 1.

##  VI Discussion and summary

In summary, in this paper we derived the properties of electronic nematic
order in all 32 crystallographic point groups and in the non-crystallographic
point groups associated with quasicrystals. We expressed the Fermi surface
distortion patterns caused by nematic order in terms of a general three-
dimensional nematic order parameter expressed in tensorial form in Eq. (21)
and in five-component vectorial form in Eq. (24). We also established the
critical properties of the allowed nematic transitions in these point groups.
These latter results are summarized in two tables: Table 1, for
crystallographic point groups, and Table 2, for the point groups of
quasicrystals.

The crystallographic point groups shown in Table 1 are organized in six blocks
corresponding to the cubic, hexagonal, trigonal, tetragonal, orthorhombic and
monoclinic crystal systems; we do not show the triclinic crystal system since
it does not support any nematic transition . For each point group in the first
column, we list all nematic order parameters that transform as a non-trivial
irreducible representation (IR, second column), together with the
corresponding basis for that order parameter (third column). The possible
basis are: the full five-component nematic order parameter
$\boldsymbol{d}=\left(d_{1},\,d_{2},\,d_{3},\,d_{4},\,d_{5}\right)^{T}$, whose
elements correspond, respectively, to the usual charge quadrupolar order
parameters $d_{\frac{1}{\sqrt{3}}(2z^{2}-x^{2}-y^{2})}$, $d_{x^{2}-y^{2}}$,
$d_{2yz}$, $d_{2xz}$, and $d_{2xy}$; the $E_{g}$ and $T_{2g}$ cubic bases
$\boldsymbol{d}^{e}=\left(d_{1},\,d_{2}\right)^{T}$ and
$\boldsymbol{d}^{t}=\left(d_{3},\,d_{4},\,d_{5}\right)^{T}$; and the in-plane
and out-of-plane axial bases
$\boldsymbol{d}^{\mathrm{ip}}=\left(d_{2},\,d_{5}\right)^{T}$, and
$\boldsymbol{d}^{\mathrm{op}}=\left(d_{3},\,d_{4}\right)^{T}$. Curly brackets
indicate that the enclosed nematic order parameters transform as the same IR.
The fourth column lists the universality class of the respective nematic
Landau expansion: $Z_{2}$-Ising [e.g., Eq. (181)]; $Z_{3}$-clock, which is
equivalent to $Z_{3}$-Potts [e.g., Eq. (86) and (132)]; $Z_{4}$-Potts [e.g.,
Eq. (76)]; $Z_{4}$-clock [e.g., Eq. (188)]; and $Z_{6}$-clock [e.g., Eq.
(140)]. The asterisk notation in $Z_{q}^{*}$-clock indicates that the clock
term has an offset angle [e.g., Eq. (328)], which is always the case when the
group lacks twofold in-plane rotation axes and the nematic order parameter
transforms as a complex IR (which is indicated by an overbar in the second
column). The fifth column shows the residual point group (PG) after the onset
of nematic order, which indicates the set of symmetry elements that remain
intact after nematicity is established.

Finally, in the last column of Table 1, we list the transformation matrices of
the nematic order parameters expressed in terms of the transformation matrices
of the coordinate vector. To understand what this entails, consider first the
five-component nematic order parameter $\boldsymbol{d}$. It transforms
according to the IR $\Gamma_{j=2}^{+}$ of the fully-isotropic orthogonal group
$\mathsf{O(3)}$, and therefore transforms with the $5\times 5$ matrices
$\mathcal{R}_{+,j=2}(g)=\mathcal{R}_{j=2}(\vartheta,\hat{\boldsymbol{\ell}})=\exp\left(-\mathsf{i}\vartheta\,\boldsymbol{J}^{(2)}\cdot\hat{\boldsymbol{\ell}}\right)$,
where $g=\big{(}\vartheta,\hat{\boldsymbol{\ell}},\mathcal{I}\big{)}$ are the
symmetry elements of $\mathsf{O(3)}$, see Appendix A for more details. Moving
from isotropic space to point groups has two effects: first, it restricts the
symmetry elements $g$ to the discrete set of operations that characterize the
group. Second, it causes the nematic order parameter to be decomposed into
different irreducible channels. Accordingly, in the point groups, the $5\times
5$ matrices $\mathcal{R}_{+,j=2}(g)$ assume a block-diagonal form, where each
block is labeled by the IR characterizing the nematic order parameter. For
one-dimensional IRs, one simply recovers the characters of the respective IR.
For multi-dimensional IRs, the transformation matrices appear as blocks
constructed from $\mathcal{R}_{+,j=2}(g)$. In particular, for the four multi-
component nematic bases $\boldsymbol{d}^{e}$, $\boldsymbol{d}^{t}$,
$\boldsymbol{d}^{\mathrm{ip}}$, and $\boldsymbol{d}^{\mathrm{op}}$ used here,
we have:

$\displaystyle\mathcal{R}_{\boldsymbol{d}^{e}}(g)$
$\displaystyle=\left(\begin{array}[]{cc}[\mathcal{R}_{+,j=2}(g)]_{11}&[\mathcal{R}_{+,j=2}(g)]_{12}\\\
{}[\mathcal{R}_{+,j=2}(g)]_{21}&[\mathcal{R}_{+,j=2}(g)]_{22}\end{array}\right),$
(246) $\displaystyle\mathcal{R}_{\boldsymbol{d}^{t}}(g)$
$\displaystyle=\left(\begin{smallmatrix}[\mathcal{R}_{+,j=2}(g)]_{33}&[\mathcal{R}_{+,j=2}(g)]_{34}&[\mathcal{R}_{+,j=2}(g)]_{35}\\\
{}[\mathcal{R}_{+,j=2}(g)]_{43}&[\mathcal{R}_{+,j=2}(g)]_{44}&[\mathcal{R}_{+,j=2}(g)]_{45}\\\
{}[\mathcal{R}_{+,j=2}(g)]_{53}&[\mathcal{R}_{+,j=2}(g)]_{54}&[\mathcal{R}_{+,j=2}(g)]_{55}\end{smallmatrix}\right),$
(247) $\displaystyle\mathcal{R}_{\boldsymbol{d}^{\mathrm{ip}}}(g)$
$\displaystyle=\left(\begin{array}[]{cc}[\mathcal{R}_{+,j=2}(g)]_{22}&[\mathcal{R}_{+,j=2}(g)]_{25}\\\
{}[\mathcal{R}_{+,j=2}(g)]_{52}&[\mathcal{R}_{+,j=2}(g)]_{55}\end{array}\right),$
(250) $\displaystyle\mathcal{R}_{\boldsymbol{d}^{\mathrm{op}}}(g)$
$\displaystyle=\left(\begin{array}[]{cc}[\mathcal{R}_{+,j=2}(g)]_{33}&[\mathcal{R}_{+,j=2}(g)]_{34}\\\
{}[\mathcal{R}_{+,j=2}(g)]_{43}&[\mathcal{R}_{+,j=2}(g)]_{44}\end{array}\right).$
(253)

In principle, we could simply use these matrices to define the transformation
properties of the nematic order parameters listed in the table. However, to
construct invariants involving the coupling between the nematic order
parameter and other physical quantities, it is often more convenient to
express these four nematic transformation matrices (246)-(253) in terms of
transformation matrices of the coordinate vector. The latter transforms
according to the IR $\Gamma_{j=1}^{-}$ of $\mathsf{O(3)}$. As a result, its
transformation matrices for any symmetry element
$g=\big{(}\vartheta,\hat{\boldsymbol{\ell}},\mathcal{I}\big{)}$ of
$\mathsf{O(3)}$ are given by:

$\displaystyle\mathcal{R}_{g}^{(3)}$ $\displaystyle=\mathcal{R}_{-,j=1}(g),$
(254) $\displaystyle\mathcal{R}_{g}^{(2)}$
$\displaystyle=\left(\begin{array}[]{cc}[\mathcal{R}_{-,j=1}(g)]_{11}&[\mathcal{R}_{-,j=1}(g)]_{12}\\\
{}[\mathcal{R}_{-,j=1}(g)]_{21}&[\mathcal{R}_{-,j=1}(g)]_{22}\end{array}\right),$
(257)

where
$\mathcal{R}_{-,j=1}(g)=\mathcal{I}\mathcal{R}_{j=1}(\vartheta,\hat{\boldsymbol{\ell}})$
and
$\mathcal{R}_{j=1}(\vartheta,\hat{\boldsymbol{\ell}})=\exp\left(-\mathsf{i}\vartheta\,\boldsymbol{J}^{(1)}\cdot\hat{\boldsymbol{\ell}}\right)$,
see Eq. (261) in Appendix A. Here, we introduced two different matrices to
distinguish between the cases of a 3D-vector, which transforms with
$\mathcal{R}_{g}^{(3)}$ (254) and is relevant for the polyhedral point groups,
and a 2D in-plane vector, which transforms with $\mathcal{R}_{g}^{(2)}$ (257)
and is relevant for the axial point groups. The key point illustrated by the
last column of Table 1 is that, for most of the nematic order parameters, the
transformation matrices (246)-(253) can be directly related to either
$\mathcal{R}_{g}^{(3)}$ or $\mathcal{R}_{g}^{(2)}$. For instance, for the
nematic order parameter that transforms as the $T_{2g}$ IR of
$\mathsf{O_{h}}$, the last entry on the row means that
$\mathcal{R}_{\boldsymbol{d}^{t}}(g)=\mathcal{R}_{g}^{A_{2u}}\mathcal{R}_{g}^{(3)}$,
where $\mathcal{R}_{g}^{A_{2u}}$ are the characters of the $A_{2u}$ IR.
Similarly, for the nematic order parameter transforming as the $E$ IR of
$\mathsf{C_{4v}}$, we have
$\mathcal{R}_{\boldsymbol{d}^{\mathrm{op}}}(g)=\sigma^{x}\mathcal{R}_{g}^{(2)}\sigma^{x}$,
where $\sigma^{j}$ are Pauli matrices.

To illustrate the usefulness of this representation of the transformation
matrices, let us construct invariants involving the coupling between the
nematic order parameters and the electric polarization $\boldsymbol{P}$, which
transforms like a vector. To make the example more transparent, we consider
the hexagonal point group $\mathsf{C_{3h}}$ and the nematic in-plane order
parameter $\boldsymbol{d}^{\mathrm{ip}}$. This doublet transforms according to
the IR $E^{\prime}$ and via the transformation matrix
$\mathcal{R}_{\boldsymbol{d}^{\mathrm{ip}}}(g)=\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}$,
as shown in Table 1. The in-plane polarization
$\boldsymbol{P}_{\mathrm{ip}}=(P_{x},P_{y})$ , on the other hand, transforms
according to the same IR $E^{\prime}$, but via the matrices
$\mathcal{R}^{\boldsymbol{P}_{\mathrm{ip}}}(g)=\mathcal{R}_{g}^{(2)}$ that
characterize a vector. Indeed, $\mathcal{R}^{\boldsymbol{P}_{\mathrm{ip}}}(g)$
is related to $\mathcal{R}_{\boldsymbol{d}^{\mathrm{ip}}}(g)$ through a
similarity transformation, the Pauli matrix $\sigma^{z}$. Since both
$\boldsymbol{d}^{\mathrm{ip}}$ and $\boldsymbol{P}_{\mathrm{ip}}$ transform as
$E^{\prime}$, there must exist a bilinear invariant. However, due to the
difference in transformation matrices, this invariant is not simply
$(\boldsymbol{d}^{\mathrm{ip}})^{T}\cdot\boldsymbol{P}_{\mathrm{ip}}$, but
instead it is
$(\boldsymbol{d}^{\mathrm{ip}})^{T}\sigma^{z}\boldsymbol{P}_{\mathrm{ip}}$,
where the additional Pauli matrix $\sigma^{z}$ is needed to compensate for the
similarity transformation relating
$\mathcal{R}^{\boldsymbol{P}_{\mathrm{ip}}}(g)$ and
$\mathcal{R}_{\boldsymbol{d}^{\mathrm{ip}}}(g)$. As a second example, consider
the coupling between the nematic out-of-plane doublet
$\boldsymbol{d}^{\mathrm{op}}$ and the out-of-plane polarization $P_{z}$
within the same point group $\mathsf{C_{3h}}$. The nematic order parameter
transforms according to the IR $E^{\prime\prime}$ and via the transformation
matrix
$\mathcal{R}_{\boldsymbol{d}^{\mathrm{op}}}(g)=\mathcal{R}_{g}^{A^{\prime\prime}}\big{(}\sigma^{z}\mathcal{R}_{g}^{(2)}\sigma^{z}\big{)}$,
see Table 1. The out-of-plane polarization $P_{z}$, on the other hand,
transforms as $A^{\prime\prime}$ and via the $1\times 1$ “matrix”
$\mathcal{R}^{P_{z}}(g)=\mathcal{R}_{g}^{A^{\prime\prime}}$. Using the fact
that $\boldsymbol{P}_{\mathrm{ip}}$ transforms according to
$\mathcal{R}_{g}^{(2)}$, we can readily construct the invariant
$P_{z}\,\left(\boldsymbol{d}^{\mathrm{op}}\right)^{T}\sigma^{z}\boldsymbol{P}_{\mathrm{ip}}$.

Table 2 has the same layout as Table 1, but refers to the quasicrystalline
point groups considered in this work. They correspond to octagonal, decagonal,
dodecagonal, and icosahedral quasicrystals.

The group-theoretical classification presented here can be valuable in
searching for new nematic systems in material databases. More broadly, our
results offer interesting insights into which types of lattice may realize
exotic electronic nematic phenomena. Indeed, for any of the nematic doublets
that behave as a $Z_{q}$-clock order parameter with $q\geq 4$, the nematic
transition in 3D belongs to the XY universality class, which in turn also
describes the nematic transition of a 2D isotropic model. On the other hand,
in a 2D system with a $Z_{q}$-clock nematic order parameter, a critical
nematic phase with quasi-long-range order precedes the onset of long-range
order for $q\geq 5$, whereas for $q=4$ one obtains an Ashkin-Teller nematic
model, whose critical properties are described by non-universal critical
exponents. From Tables 1 and 2, we see that this condition is satisfied by
$\left(d_{2yz},\,d_{2xz}\right)$ nematicity in hexagonal and tetragonal
crystals, as well as $\left(d_{2yz},\,d_{2xz}\right)$ and
$\big{(}d_{x^{2}-y^{2}},\,d_{2xy}\big{)}$ nematicity in octagonal, decagonal,
and dodecagonal quasicrystals. In this regard, it will be interesting to
investigate the properties of the Landau expansion (58) of the icosahedral
nematic order parameter, since the mean-field ground states are either sixfold
or tenfold degenerate.

Experimentally, these properties should be directly accessible via probes that
measure the nematic susceptibility in different channels, such as elasto-
resistance [9] and Raman spectroscopy [117]. Theoretically, one expects that
nematic fluctuations will be enhanced near the 3D-XY nematic transition, since
there will be not only soft longitudinal (i.e., amplitude) fluctuations but
also soft transverse (i.e., phase) fluctuations associated with the emergent
continuous symmetry. Given that nematic fluctuations have been proposed as
potential drivers of non-Fermi liquid behavior [49, 50, 51, 69] and pairing
[118, 119], the investigation of emergent XY-nematicity could provide further
insights into these problems as well. As for the exotic 2D behavior associated
with $q\geq 4$ clock in-plane nematicity, the most promising realizations
would be in $30^{\circ}$-twisted hexagonal bilayers ($Z_{6}$-clock) and
$45^{\circ}$-twisted tetragonal bilayers ($Z_{4}$-clock). Both settings have
been recently realized experimentally by twisting graphene [92] and monolayers
of cuprate BSCCO [96]. Interestingly, nematic order has been observed in
underdoped BSCCO [120].

One important ingredient that was not included in this paper and that deserves
further investigation is the role of strain, both extrinsic and intrinsic. For
instance, it was recently shown that the impact of uniaxial strain on in-plane
$Z_{3}$-Potts/clock nematics is fundamentally different from the well-studied
case of Ising-nematics [81]. This raises the broader question of how external
strain can be used to probe and modify the behavior of $Z_{q}$-clock nematic
order parameters. Moreover, long-wavelength quantized strain fluctuations,
which are manifested as acoustic phonons, are known to significantly change
the critical properties of $Z_{2}$ and $Z_{3}$ electronic nematic transitions
in tetragonal, hexagonal, and trigonal crystals by promoting long-range
nematic interactions [121, 122, 123, 124, 74, 107]. It will be important to
perform similar analyses in cubic crystals as well as quasicrystals, which
also host phason modes on top of phonon modes [125, 126]. Finally, crystalline
defects such as dislocations, vacancies, and interstitials create
inhomogeneous strain distributions, which act as random nematic conjugate
fields [64, 127, 128]. While most studies have focused on a random-field
Ising-model description of this rich phenomenon, it will be important to
develop models that can capture the long-range and correlated nature of the
various components of the strain fields generated by these defects.

###### Acknowledgements.

We thank T. Birol and J. Viñals for fruitful discussions. This work was
primarily supported by the U. S. Department of Energy, Office of Science,
Basic Energy Sciences, Materials Sciences and Engineering Division, under
Award No. DE-SC0020045 (M.H., A.R., and R.M.F.). D.F.A. was supported by
Department of Energy, Office of Basic Energy Science, Division of Materials
Sciences and Engineering under Award No. DE-SC0021971.

##  Appendix A Generators of SO$(3)$ and parametrization of $\boldsymbol{d}$

In this Appendix, we give further details about the generators of
$\mathsf{SO(3)}$. Within the full rotation group $\mathsf{SO(3)}$, the
transformation matrices associated with an irreducible representation (IR)
$\Gamma_{j}$ are given by
$\mathcal{R}_{j}(\vartheta,\hat{\boldsymbol{\ell}})=\exp\left(-\mathsf{i}\vartheta\,\boldsymbol{J}^{(j)}\cdot\hat{\boldsymbol{\ell}}\right)$,
parametrized by a rotation angle $\vartheta$ around a unit-length rotation
axis $\hat{\boldsymbol{\ell}}$. The three generators $J_{\mu}^{(j)}$ satisfy
the Lie algebra

$\displaystyle\left[J_{\mu_{1}}^{(j)},J_{\mu_{2}}^{(j)}\right]$
$\displaystyle=\mathsf{i}\,\epsilon_{\mu_{1}\mu_{2}\mu_{3}}J_{\mu_{3}}^{(j)},$
(258)

with the anti-symmetric Levi-Civita tensor $\epsilon_{\mu_{1}\mu_{2}\mu_{3}}$.
In principle, the matrices $J_{\mu}^{(j)}$ can be defined in the “Lie basis”
using the standard relationships

$\displaystyle\hat{J}_{z}^{(j)}\,|j,m\rangle$ $\displaystyle=m\,|j,m\rangle,$
(259) $\displaystyle\hat{J}_{\pm}^{(j)}\,|j,m\rangle$
$\displaystyle=\sqrt{j(j+1)-m(m\pm 1)}\,|j,m\pm 1\rangle,$ (260)

where $\hat{J}_{\pm}^{(j)}=\hat{J}_{x}^{(j)}\pm\mathsf{i}\hat{J}_{y}^{(j)}$.
However, for the transformation matrices to be orthogonal, i.e.
$\mathcal{R}_{j}^{T}(\vartheta,\hat{\boldsymbol{\ell}})=\mathcal{R}_{j}(-\vartheta,\hat{\boldsymbol{\ell}})=\mathcal{R}_{j}^{-1}(\vartheta,\hat{\boldsymbol{\ell}})$,
one needs to impose the requirement
$\big{(}J_{\mu}^{(j)}\big{)}^{T}=\big{(}J_{\mu}^{(j)}\big{)}^{*}=-J_{\mu}^{(j)}$,
i.e. the generators need to be anti-symmetric, and thus, a similarity
transformation has to be applied on Eqs. (259)-(260). We now list the
resulting generators. For $j=0$ the three generators are
$J_{x}^{(0)}=J_{y}^{(0)}=J_{z}^{(0)}=0$, which ensures that a $\Gamma_{j=0}$
scalar is invariant under all rotations. For $j=1$ the three generators are
the well-known anti-symmetric Gell-Mann matrices

$\displaystyle J_{x}^{(1)}$
$\displaystyle=\mathsf{i}\left(\begin{smallmatrix}0&0&0\\\ 0&0&-1\\\
0&1&0\end{smallmatrix}\right),$ $\displaystyle J_{y}^{(1)}$
$\displaystyle=\mathsf{i}\left(\begin{smallmatrix}0&0&1\\\ 0&0&0\\\
-1&0&0\end{smallmatrix}\right),$ $\displaystyle J_{z}^{(1)}$
$\displaystyle=\mathsf{i}\left(\begin{smallmatrix}0&-1&0\\\ 1&0&0\\\
0&0&0\end{smallmatrix}\right),$ (261)

which can also be expressed in terms of the anti-symmetric Levi-Cevita tensor
$\big{(}J_{\mu_{1}}^{(1)}\big{)}_{\mu_{2}\mu_{3}}=-\mathsf{i}\epsilon_{\mu_{1}\mu_{2}\mu_{3}}$.
The three generators for $j=2$ are given by

$\displaystyle J_{x}^{(2)}$
$\displaystyle=\mathsf{i}\left(\begin{smallmatrix}0&0&\sqrt{3}&0&0\\\
0&0&1&0&0\\\ -\sqrt{3}&-1&0&0&0\\\ 0&0&0&0&1\\\
0&0&0&-1&0\end{smallmatrix}\right),$ $\displaystyle J_{y}^{(2)}$
$\displaystyle=\mathsf{i}\left(\begin{smallmatrix}0&0&0&-\sqrt{3}&0\\\
0&0&0&1&0\\\ 0&0&0&0&-1\\\ \sqrt{3}&-1&0&0&0\\\
0&0&1&0&0\end{smallmatrix}\right),$ $\displaystyle J_{z}^{(2)}$
$\displaystyle=\mathsf{i}\left(\begin{smallmatrix}0&0&0&0&0\\\ 0&0&0&0&-2\\\
0&0&0&1&0\\\ 0&0&-1&0&0\\\ 0&2&0&0&0\end{smallmatrix}\right),$ (262)

and the ones for $j=4$ are

$\displaystyle J_{x}^{(4)}$
$\displaystyle=\frac{\mathsf{i}}{\sqrt{2}}\left(\begin{smallmatrix}0&0&0&0&2\sqrt{5}&0&0&0&0\\\
0&0&2&0&0&0&0&0&0\\\ 0&-2&0&-\sqrt{7}&0&0&0&0&0\\\ 0&0&\sqrt{7}&0&3&0&0&0&0\\\
-2\sqrt{5}&0&0&-3&0&0&0&0&0\\\ 0&0&0&0&0&0&3&0&0\\\
0&0&0&0&0&-3&0&-\sqrt{7}&0\\\ 0&0&0&0&0&0&\sqrt{7}&0&2\\\
0&0&0&0&0&0&0&-2&0\end{smallmatrix}\right),$ (263) $\displaystyle J_{y}^{(4)}$
$\displaystyle=\frac{\mathsf{i}}{\sqrt{2}}\left(\begin{smallmatrix}0&0&0&0&0&-2\sqrt{5}&0&0&0\\\
0&0&0&0&0&0&0&2&0\\\ 0&0&0&0&0&0&\sqrt{7}&0&-2\\\ 0&0&0&0&0&3&0&-\sqrt{7}&0\\\
0&0&0&0&0&0&-3&0&0\\\ 2\sqrt{5}&0&0&-3&0&0&0&0&0\\\
0&0&-\sqrt{7}&0&3&0&0&0&0\\\ 0&-2&0&\sqrt{7}&0&0&0&0&0\\\
0&0&2&0&0&0&0&0&0\end{smallmatrix}\right),$ (264) $\displaystyle J_{z}^{(4)}$
$\displaystyle=\mathsf{i}\,\left(\begin{smallmatrix}0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&-4\\\ 0&0&0&0&0&0&0&3&0\\\ 0&0&0&0&0&0&-2&0&0\\\
0&0&0&0&0&1&0&0&0\\\ 0&0&0&0&-1&0&0&0&0\\\ 0&0&0&2&0&0&0&0&0\\\
0&0&-3&0&0&0&0&0&0\\\ 0&4&0&0&0&0&0&0&0\end{smallmatrix}\right).$ (265)

Clearly, the generators are never unique in the sense that they can always be
rotated via an orthogonal matrix $O^{(j)}$ according to
$[O^{(j)}]^{T}\boldsymbol{J}^{(j)}O^{(j)}$, which corresponds to a basis
rotation. Throughout this work, the generators are fixed according to Eqs.
(261)-(265).

Next, we determine the bilinears $\boldsymbol{D}^{j=2}$ and
$\boldsymbol{D}^{j=4}$ that occur in the decomposition (26) and which were
used in the derivation of the Landau expansion of the 3D isotropic nematic
system. The transformation conditions are structurally similar to Eq. (5). For
example, for the $\Gamma_{j=4}^{+}$ bilinear
$D^{j=4,\nu}=d_{i}\Lambda_{i,i^{\prime}}^{j=4,\nu}d_{i^{\prime}}$, one needs
to solve the equation

$\displaystyle\mathcal{R}_{+,j=2}^{T}(g)\Lambda^{j=4,\nu}\mathcal{R}_{+,j=2}(g)$
$\displaystyle=\mathcal{R}_{+,j=4}(g)_{\nu\nu^{\prime}}\Lambda^{j=4,\nu^{\prime}},$
(266)

where $\Lambda^{j=4,\nu}$ are $5\times 5$-dimensional matrices with
$i,i^{\prime}=1,\dots,5$ and $\nu,\nu^{\prime}=1,\dots,9$. Recall that the
$\Gamma_{j=2}^{+}$ nematic order parameter $\boldsymbol{d}$ transforms via the
matrices
$\mathcal{R}_{+,j=2}(g)=\mathcal{R}_{j=2}(\vartheta,\hat{\boldsymbol{\ell}})$
for each element
$g=\big{(}\vartheta,\hat{\boldsymbol{\ell}},\mathcal{I}\big{)}$.

A straightforward but tedious calculation then yields
$D^{j=0}=\left|\boldsymbol{d}\right|^{2}$, $\boldsymbol{D}^{j=2}$ given by Eq.
(29), and

$\displaystyle\boldsymbol{D}^{j=4}$
$\displaystyle=\frac{\sqrt{5}}{6}\left(\begin{smallmatrix}\frac{1}{\sqrt{5}}\left[6d_{1}^{2}-4\left(d_{3}^{2}+d_{4}^{2}\right)+\left(d_{2}^{2}+d_{5}^{2}\right)\right]\\\
\sqrt{7}\left(d_{2}^{2}-d_{5}^{2}\right)\\\
\sqrt{14}\left(d_{2}d_{3}+d_{4}d_{5}\right)\\\
2\sqrt{3}d_{1}d_{2}-2\left(d_{3}^{2}-d_{4}^{2}\right)\\\
\sqrt{2}\left[\left(2\sqrt{3}d_{1}+d_{2}\right)d_{3}-d_{4}d_{5}\right]\\\
\sqrt{2}\left[\left(2\sqrt{3}d_{1}-d_{2}\right)d_{4}-d_{3}d_{5}\right]\\\
2\sqrt{3}d_{1}d_{5}+4d_{3}d_{4}\\\
\sqrt{14}\left(d_{2}d_{4}-d_{3}d_{5}\right)\\\
2\sqrt{7}d_{2}d_{5}\end{smallmatrix}\right).$ (267)

Importantly, all three bilinears have the same magnitude,
$\left(D^{j=0}\right)^{2}=\left(\boldsymbol{D}^{j=2}\right)^{2}=\left(\boldsymbol{D}^{j=4}\right)^{2}=\left|\boldsymbol{d}\right|^{4}$.
From the condition (266) we see that a basis rotation $O^{(j=4)}$ on the
generators $\boldsymbol{J}^{(j=4)}$ would rotate the basis state
$\boldsymbol{D}^{j=4}$ (267) into $O^{(j=4)}\boldsymbol{D}^{j=4}$.

To finish this Appendix, we further discuss the parameterization of the five-
component nematic order parameter within the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation (24). Recall
that:

$\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{\\}}.$
(268)

The three angles describing the three orthonormal eigenvectors
$\boldsymbol{n}$, $\boldsymbol{m}$, $\boldsymbol{l}$ can be chosen in several
ways. The parameterization employed in this work makes use of the spherical
unit vectors

$\displaystyle\hat{\boldsymbol{e}}_{r}$
$\displaystyle=\left(\begin{smallmatrix}\cos\varphi\sin\theta\\\
\sin\varphi\sin\theta\\\ \cos\theta\end{smallmatrix}\right),$
$\displaystyle\hat{\boldsymbol{e}}_{\varphi}$
$\displaystyle=\left(\begin{smallmatrix}\sin\varphi\\\ -\cos\varphi\\\
0\end{smallmatrix}\right),$ $\displaystyle\hat{\boldsymbol{e}}_{\theta}$
$\displaystyle=\left(\begin{smallmatrix}\cos\varphi\cos\theta\\\
\sin\varphi\cos\theta\\\ -\sin\theta\end{smallmatrix}\right).$ (269)

One of the eigenvectors, e.g. $\boldsymbol{l}$, can always be chosen to be
aligned with the radial vector $\hat{\boldsymbol{e}}_{r}$, i.e.
$\boldsymbol{l}=\hat{\boldsymbol{e}}_{r}$. Conversely, the remaining two
vectors can be arbitrarily rotated about this axis by an angle $\eta$:

$\displaystyle\boldsymbol{l}$ $\displaystyle=\hat{\boldsymbol{e}}_{r},$
$\displaystyle\boldsymbol{n}$
$\displaystyle=\cos\eta\,\hat{\boldsymbol{e}}_{\varphi}+\sin\eta\,\hat{\boldsymbol{e}}_{\theta},$
$\displaystyle\boldsymbol{m}$
$\displaystyle=-\sin\eta\,\hat{\boldsymbol{e}}_{\varphi}+\cos\eta\,\hat{\boldsymbol{e}}_{\theta}.$
(270)

Thus, within the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation
the five degrees of freedom are encoded in the magnitude $|\boldsymbol{d}|$
and the four angles $\alpha,\varphi,\theta,\eta$. To ensure a one-to-one
mapping between the two representations, (19) and (21), the parameter ranges
have to be restricted according to 111For the limiting values of $\theta$ the
rotation angles are actually more restricted: For $\theta=0$ we constrain
$\varphi\in[0,0]$ and $\eta\in[0,\pi]$, and for $\theta=\frac{\pi}{2}$ we
constrain $\varphi\in[0,2\pi]$ and $\eta\in[0,\frac{\pi}{2}]$.

$\displaystyle|\boldsymbol{d}|$ $\displaystyle\geq 0,$ $\displaystyle\alpha$
$\displaystyle\in\left[0,\frac{\pi}{3}\right],$ $\displaystyle\theta$
$\displaystyle\in\left[0,\frac{\pi}{2}\right],$ $\displaystyle\varphi$
$\displaystyle\in\left[0,2\pi\right],$ $\displaystyle\eta$
$\displaystyle\in\left[0,\pi\right].$ (271)

##  Appendix B Symmetrized products

In this Appendix, we discuss in more detail the product decomposition
associated with a generic order parameter
$\boldsymbol{\eta}=(\eta_{1},\dots,\eta_{\dim n})$ that transforms according
to an IR $n$ of some group. Because the order parameter components commute,
i.e. $\eta_{i}\eta_{j}=\eta_{j}\eta_{i}$, the number of product components is
reduced. For example, for a two-component order parameter
$\boldsymbol{\eta}=(\eta_{1},\eta_{2})$ one would naively expect $2\times 2=4$
product components, while in fact, there are only $3$: $\eta_{1}^{2}$,
$\eta_{2}^{2}$ and $\eta_{1}\eta_{2}$, since the antisymmetric combination
$\eta_{1}\eta_{2}-\eta_{2}\eta_{1}$ vanishes. Here, we use group theory to
derive expressions for the product decompositions in the symmetric channel,
see e.g. Ref. [130]. Doing so automatically excludes vanishing contributions,
such as the antisymmetric combination $\eta_{1}\eta_{2}-\eta_{2}\eta_{1}$ in
the above example, and removes redundancies associated with double-counting.
To illustrate the last point, consider the symmetrized bilinears in the above
example: $N_{1}=\eta_{1}^{2}+\eta_{2}^{2}$ and
$\boldsymbol{N}_{2}=\left(\eta_{1}^{2}-\eta_{2}^{2},2\eta_{1}\eta_{2}\right)$.
Then, a straightforward decomposition would suggest that two invariants exist
in fourth order, $N_{1}^{2}$ and $\boldsymbol{N}_{2}\cdot\boldsymbol{N}_{2}$,
while in fact they are identical since
$N_{1}^{2}=\boldsymbol{N}_{2}\cdot\boldsymbol{N}_{2}=\left(\eta_{1}^{2}+\eta_{2}^{2}\right)^{2}$.
Such double-countings are removed by using the symmetrized decomposition. As
we are interested in Landau expansions, we consider products up to sixth
order, but extension to higher orders is straightforward.

Let us begin with the second-order product $P_{2,ij}=\eta_{i}\eta_{j}$. Using
the order parameter transformation
$\eta_{i}^{\prime}=\mathcal{R}_{n}(g)_{ii^{\prime}}\eta_{i^{\prime}}$ under a
symmetry element $g$, one finds the product to transform as

$\displaystyle P_{2,ij}^{\prime}$
$\displaystyle=\mathcal{R}^{\eta^{2}}(g)_{ij,i^{\prime}j^{\prime}}\,P_{2,i^{\prime}j^{\prime}},$
(272)

where
$\mathcal{R}^{\eta^{2}}(g)_{ij,i^{\prime}j^{\prime}}=\mathcal{R}_{n}(g)_{ii^{\prime}}\mathcal{R}_{n}(g)_{jj^{\prime}}$.
The transformation matrices $\mathcal{R}^{\eta^{2}}(g)$ themselves form a
representation of the group, i.e. it holds that
$\mathcal{R}^{\eta^{2}}(g_{1}g_{2})=\mathcal{R}^{\eta^{2}}(g_{1})\mathcal{R}^{\eta^{2}}(g_{1})$,
see e.g. Ref. [130] for details. Correspondingly, the characters of this
representation result from the trace,

$\displaystyle\chi^{\eta^{2}}(g)$
$\displaystyle=\sum_{i,j}\mathcal{R}^{\eta^{2}}(g)_{ij,ij}=\sum_{i,j}\mathcal{R}_{n}(g)_{ii}\mathcal{R}_{n}(g)_{jj}=\chi_{n}^{2}(g).$
(273)

The formula (273) contains all symmetry channels, including the ones that
eventually vanish due to the commutativity of the components of
$\boldsymbol{\eta}$. To have these removed, we symmetrize the transformation
matrix $\mathcal{R}^{\eta^{2}}(g)$ before computing its character.
Symmetrization involves the addition of all permutations in Eq. (272) with
respect to $(i^{\prime},j^{\prime})$, divided by the number of permutations,
i.e.

$\displaystyle\mathcal{R}_{s}^{\eta^{2}}(g)_{ij,i^{\prime}j^{\prime}}$
$\displaystyle=\frac{1}{2!}\left[\mathcal{R}^{\eta^{2}}(g)_{ij,i^{\prime}j^{\prime}}+\mathcal{R}^{\eta^{2}}(g)_{ij,j^{\prime}i^{\prime}}\right].$
(274)

Then, the characters for the symmetrized product (274) become

$\displaystyle\chi_{s}^{\eta^{2}}(g)$
$\displaystyle=\sum_{i,j}\mathcal{R}_{s}^{\eta^{2}}(g)_{ij,ij}=\frac{1}{2}\left(\chi_{n}^{2}(g)+\chi_{n}(g^{2})\right).$
(275)

Thus, the decomposition of a symmetrized product (275) only contains a subset
of the original symmetry channels associated with the bare product (273). The
symmetrized decomposition is accomplished by using
$\mathcal{R}_{s}^{\eta^{2}}(g)$ rather than $\mathcal{R}^{\eta^{2}}(g)$.

For the higher-order products such as $P_{3,ijk}=\eta_{i}\eta_{j}\eta_{k}$ and
$P_{4,ijkl}=\eta_{i}\eta_{j}\eta_{k}\eta_{l}$, one proceeds analogously. As in
the case above, the corresponding transformation matrices

$\displaystyle\mathcal{R}^{\eta^{3}}(g)_{ijk,i^{\prime}j^{\prime}k^{\prime}}$
$\displaystyle=\mathcal{R}_{n}(g)_{ii^{\prime}}\mathcal{R}_{n}(g)_{jj^{\prime}}\mathcal{R}_{n}(g)_{kk^{\prime}},$
(276)
$\displaystyle\mathcal{R}^{\eta^{4}}(g)_{ijkl,i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$
$\displaystyle=\mathcal{R}_{n}(g)_{ii^{\prime}}\mathcal{R}_{n}(g)_{jj^{\prime}}\mathcal{R}_{n}(g)_{kk^{\prime}}\mathcal{R}_{n}(g)_{ll^{\prime}},$
(277)

form representations of the group with characters

$\displaystyle\chi^{\eta^{3}}(g)$
$\displaystyle=\sum_{ijk}\mathcal{R}^{\eta^{3}}(g)_{ijk,ijk}=\chi_{n}^{3}(g),$
(278) $\displaystyle\chi^{\eta^{4}}(g)$
$\displaystyle=\sum_{ijkl}\mathcal{R}^{\eta^{4}}(g)_{ijkl,ijkl}=\chi_{n}^{4}(g).$
(279)

The symmetrized transformation matrices are

$\displaystyle\mathcal{R}_{s}^{\eta^{3}}\\!\\!(g)_{ijk,i^{\prime}j^{\prime}k^{\prime}}$
$\displaystyle=\frac{1}{3!}\left[\mathcal{R}^{\eta^{3}}\\!\\!(g)_{ijk,i^{\prime}j^{\prime}k^{\prime}}+5\left(i^{\prime}\leftrightarrow
j^{\prime}\leftrightarrow k^{\prime}\right)\right],$ (280)
$\displaystyle\mathcal{R}_{s}^{\eta^{4}}\\!\\!(g)_{ijkl,i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$
$\displaystyle=\frac{1}{4!}\Big{[}\mathcal{R}^{\eta^{4}}\\!\\!(g)_{ijkl,i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$
$\displaystyle\qquad+23\left(i^{\prime}\leftrightarrow
j^{\prime}\leftrightarrow k^{\prime}\leftrightarrow l^{\prime}\right)\Big{]},$
(281)

with the number of added permutation terms indicated explicitly [cf. Eq.
(274)]. One finds the characters associated with these symmetrized products to
be

$\displaystyle\chi_{s}^{\eta^{3}}\\!\\!(g)$
$\displaystyle=\\!\sum_{ijk}\\!\mathcal{R}_{s}^{\eta^{3}}\\!\\!(g)_{ijk,ijk}\\!$
$\displaystyle=\big{[}\chi_{n}^{3}(g)+2\chi_{n}(g^{3})+3\chi_{n}(g)\chi_{n}(g^{2})\big{]}\big{/}6,$
(282) $\displaystyle\chi_{s}^{\eta^{4}}\\!\\!(g)$
$\displaystyle=\\!\sum_{ijkl}\mathcal{R}_{s}^{\eta^{4}}\\!\\!(g)_{ijkl,ijkl}=\frac{1}{24}\Big{[}\chi_{n}^{4}(g)+6\chi_{n}^{2}(g)\chi_{n}(g^{2})$
$\displaystyle\quad+8\chi_{n}(g)\chi_{n}(g^{3})+3\chi_{n}^{2}(g^{2})+6\chi_{n}(g^{4})\Big{]}.$
(283)

Similarly, one finds the fifth- and sixth-order expressions

$\displaystyle\chi_{s}^{\eta^{5}}\\!\\!(g)$
$\displaystyle=\frac{1}{5!}\Big{[}\chi_{n}^{5}(g)+10\chi_{n}^{3}(g)\chi_{n}(g^{2})+20\chi_{n}^{2}(g)\chi_{n}(g^{3})$
$\displaystyle+24\chi_{n}(g^{5})+30\chi_{n}(g)\chi_{n}(g^{4})+20\chi_{n}(g^{2})\chi_{n}(g^{3})$
$\displaystyle+15\chi_{n}(g)\chi_{n}^{2}(g^{2})\Big{]},$ (284)
$\displaystyle\chi_{s}^{\eta^{6}}\\!\\!(g)$
$\displaystyle=\frac{1}{6!}\Big{[}15\chi_{n}^{4}(g)\chi_{n}(g^{2})+120\chi_{n}(g)\chi_{n}(g^{2})\chi_{n}(g^{3})$
$\displaystyle+\chi_{n}^{6}(g)+15\chi_{n}^{3}(g^{2})+40\chi_{n}^{2}(g^{3})+45\chi_{n}^{2}(g)\chi_{n}^{2}(g^{2})$
$\displaystyle+144\chi_{n}(g)\chi_{n}(g^{5})+40\chi_{n}^{3}(g)\chi_{n}(g^{3})+120\chi_{n}(g^{6})$
$\displaystyle+90\chi_{n}(g^{2})\chi_{n}(g^{4})+90\chi_{n}^{2}(g)\chi_{n}(g^{4})\Big{]}.$
(285)

The presented derivation is not restricted to $\eta$ transforming according to
a single IR. In fact, we can equally assume that $\eta$ transforms according
to the representation $n_{1}+n_{2}+\ldots+n_{\nu}$. For example, in the case
of a complex IR, we choose $n_{1}=E$ and $n_{2}=\bar{E}$; in the case of a
degeneracy such as in a trigonal system, we have $n_{1}=E_{g}$ and
$n_{2}=E_{g}$, or even combinations thereof such as in $\mathsf{C_{3}}$ with
$n_{1}=n_{3}=E$ and $n_{2}=n_{4}=\bar{E}$. For clarity, we let
$\tilde{\mathcal{R}}_{n}(g)$ be the transformation matrices associates with
the IR $n$, and $\tilde{\chi}_{n}(g)$ the corresponding character. The nematic
order parameter $\boldsymbol{\eta}=(\eta_{1},\eta_{2}\dots,\eta_{\dim
n_{1}+\ldots+\dim n_{\nu}})$, on the other hand, transforms with the matrices
$\mathcal{R}_{n}(g)=\tilde{\mathcal{R}}_{n_{1}}(g)\oplus\ldots\oplus\tilde{\mathcal{R}}_{n_{\nu}}(g)$.
Clearly, its characters are given by
$\chi_{n}(g)=\tilde{\chi}_{n_{1}}(g)+\ldots+\tilde{\chi}_{n_{\nu}}(g)$. More
importantly, since the transformation matrices $\mathcal{R}_{n}(g)$ are block-
diagonal, we find

$\displaystyle\left[\mathcal{R}_{n}(g)\right]^{\ell}$
$\displaystyle=\left[\tilde{\mathcal{R}}_{n_{1}}(g)\right]^{\ell}\oplus\ldots\oplus\left[\tilde{\mathcal{R}}_{n_{\nu}}(g)\right]^{\ell},$
(286)

with $\ell=\\{0,1,2,\dots\\}$. As a result, the characters for any power
$\ell$ become

$\displaystyle\chi_{n}(g^{\ell})$
$\displaystyle=\tilde{\chi}_{n_{1}}(g^{\ell})+\ldots+\tilde{\chi}_{n_{\nu}}(g^{\ell}).$
(287)

We now note that the above formalism equally applies to the “mixed”
$\boldsymbol{\eta}=(\eta_{1},\eta_{2}\dots,\eta_{\dim n_{1}+\ldots+\dim
n_{\nu}})$ with transformation matrices $\mathcal{R}_{n}(g)$. The only
difference is that the characters (287) have to be inserted into the formulae
(275)-(285).

##  Appendix C Minimization of the icosahedral nematic Landau expansion

In this Appendix, we derive the mean-field phase diagram of the nematic
instability in an icosahedral group, which is presented in Fig. 3. Starting
from the symmetrized product decomposition (53)-(55) we first determine the
three bilinear combinations $D^{A_{g}}$, $\boldsymbol{D}^{H_{g},1}$ and
$\boldsymbol{D}^{H_{g},2}$. The trivial one is given by
$D^{A_{g}}=|\boldsymbol{d}|^{2}$. For the two degenerate $H_{g}$-bilinears we
choose a representation where both bilinears have equal amplitude
$\left|\boldsymbol{D}^{H_{g},1}\right|=\left|\boldsymbol{D}^{H_{g},2}\right|=\left|\boldsymbol{d}\right|^{2}$.
This condition indeed fixes the two bilinears and one finds
$\boldsymbol{D}^{H_{g},1}=\boldsymbol{D}^{j=2}$ to be identical to the
$\mathsf{O(3)}$ bilinear (29). With the amplitudes being equal, it is clear
that the obtained $\boldsymbol{D}^{H_{g},2}$ must be related to
$\boldsymbol{D}^{H_{g},1}$ through a rotation of the kind
$\boldsymbol{D}^{H_{g},2}=\mathcal{R}^{A}\tilde{\boldsymbol{D}}^{H_{g},1}$
where
$\tilde{\boldsymbol{D}}^{H_{g},1}=\boldsymbol{D}^{H_{g},1}\big{|}_{\boldsymbol{d}\rightarrow\tilde{\boldsymbol{d}}=\mathcal{R}^{B}\boldsymbol{d}}$
with rotation matrices $\mathcal{R}^{A,B}$. Clearly, this choice preserves the
magnitude as
$(\boldsymbol{D}^{H_{g},2})^{T}\boldsymbol{D}^{H_{g},2}=(\tilde{\boldsymbol{D}}^{H_{g},1})^{T}\tilde{\boldsymbol{D}}^{H_{g},1}=\left((\tilde{\boldsymbol{d}})^{T}\tilde{\boldsymbol{d}}\right)^{2}=\left((\boldsymbol{d})^{T}\boldsymbol{d}\right)^{2}$.
We find this relation to be satisfied for
$\mathcal{R}^{A}=\mathcal{R}_{5}(-\phi_{0})$ and
$\mathcal{R}^{B}=\mathcal{R}_{5}(\phi_{0})$, as presented in the main text in
Eqs. (56)-(57). Introducing 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}$,
we rewrite the Landau expansion (58) as

$\displaystyle\mathcal{S}$
$\displaystyle=\int_{\mathsf{x}}|\boldsymbol{d}|^{2}\Big{\\{}r_{0}+|\boldsymbol{d}|f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(3)}+|\boldsymbol{d}|^{2}f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(4)}\Big{\\}},$
(288)

with the direction-dependent functions

$\displaystyle f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(3)}$
$\displaystyle=g_{1}\hat{\boldsymbol{d}}\cdot\hat{\boldsymbol{D}}^{H_{g},1}+g_{2}\hat{\boldsymbol{d}}\cdot\hat{\boldsymbol{D}}^{H_{g},2},$
(289) $\displaystyle
f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(4)}$
$\displaystyle=u_{1}+u_{2}\,\hat{\boldsymbol{D}}^{H_{g},1}\cdot\hat{\boldsymbol{D}}^{H_{g},2}\,.$
(290)

Within a mean-field analysis, the cubic term in the action (288) triggers a
first-order transition at a reduced temperature $r_{0}>0$. To derive this
solution, one first solves the Landau equations for $|\boldsymbol{d}|$:

$\displaystyle|\boldsymbol{d}|_{0}$
$\displaystyle=\frac{3\big{|}f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(3)}\big{|}}{8f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(4)}}\left[1+\sqrt{1-\frac{32r_{0}f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(4)}}{9[f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(3)}]^{2}}}\right].$
(291)

where we implicitly assumed that
$f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(3)}<0$; we later
verified that this is indeed the case for the mean-field solution. The first-
order phase transition occurs at the reduced temperature $r_{0}^{c}$ where the
expansion (288) evaluated at the solution (291) vanishes. We find:

$\displaystyle r_{0}^{c}$
$\displaystyle=\max_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}\frac{[f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(3)}]^{2}}{4f_{\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}}^{(4)}},$
(292)

where the maximum function determines the optimal direction of the nematic
direction parameters $\alpha,\boldsymbol{n},\boldsymbol{m},\boldsymbol{l}$.
This optimization with respect to the four degrees of freedom has been
conducted numerically to derive the phase diagram in Fig. 3.

An analytical solution is readily available in two regions of the
$(g_{2}/g_{1},\,u_{2}/u_{1})$ parameter-space. First, it is convenient to
parameterize the rotated nematic vector
$\tilde{\boldsymbol{d}}=\mathcal{R}_{5}(\phi_{0})\boldsymbol{d}$ as

$\displaystyle\tilde{\boldsymbol{d}}$
$\displaystyle=\frac{|\boldsymbol{d}|}{\sqrt{3}}\Big{\\{}\cos\left(\tilde{\alpha}\right)\left(\tilde{n}_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}\tilde{n}_{\mu^{\prime}}\right)+\cos\left(\tilde{\alpha}+\frac{2\pi}{3}\right)\left(\tilde{m}_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}\tilde{m}_{\mu^{\prime}}\right)$
$\displaystyle\quad+\cos\left(\tilde{\alpha}+\frac{4\pi}{3}\right)\left(\tilde{l}_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}\tilde{l}_{\mu^{\prime}}\right)\Big{\\}},$
(293)

similar to the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation
(24). Here,
($\tilde{\alpha}$,$\tilde{\boldsymbol{n}}$,$\tilde{\boldsymbol{m}}$,$\tilde{\boldsymbol{l}}$)
are functions of ($\alpha$,$\boldsymbol{n}$,$\boldsymbol{m}$,$\boldsymbol{l}$)
defined implicitly through
$\tilde{\boldsymbol{d}}=\mathcal{R}_{5}(\phi_{0})\boldsymbol{d}$. This allows
us to rewrite the cubic terms in Eq. (288) as

$\displaystyle\boldsymbol{d}\cdot\boldsymbol{D}^{H_{g},1}$
$\displaystyle=|\boldsymbol{d}|^{3}\cos\left(3\alpha\right),$ (294)
$\displaystyle\boldsymbol{d}\cdot\boldsymbol{D}^{H_{g},2}$
$\displaystyle=\tilde{\boldsymbol{d}}\cdot\tilde{\boldsymbol{D}}^{H_{g},1}=|\boldsymbol{d}|^{3}\cos\left(3\tilde{\alpha}\right).$
(295)

The two parameter regimes where we can analytically derive the nematic ground
state are spanned by (i) $\mathrm{sign}g_{1}=\mathrm{sign}g_{2}$ and
$u_{2}<0$; and by (ii) $\mathrm{sign}g_{1}=-\mathrm{sign}g_{2}$ and $u_{2}>0$.
In both regimes, all direction-dependent terms in the expansion (288) can be
simultaneously minimized.

In the first parameter regime (i), the cubic terms are minimized by
$\hat{\boldsymbol{d}}\cdot\hat{\boldsymbol{D}}^{H_{g},1}=\hat{\boldsymbol{d}}\cdot\hat{\boldsymbol{D}}^{H_{g},2}=-\mathrm{sign}g_{1}$,
which corresponds to
$\cos\left(3\alpha\right)=\cos\left(3\tilde{\alpha}\right)=-\mathrm{sign}g_{1}$,
see Eqs. (294)-(295). For concreteness, we choose $g_{1}<0$, such that the
minimization gives $\alpha=\tilde{\alpha}=0$; a similar procedure can be
carried out for $g_{1}>0$. Technically, one still needs to demonstrate that
eigenvectors $\boldsymbol{n}$,$\boldsymbol{m}$,$\boldsymbol{l}$ exist for
these values. Setting $\alpha=\tilde{\alpha}=0$ in the definition
$\tilde{\boldsymbol{d}}=\mathcal{R}_{5}(\phi_{0})\boldsymbol{d}$ and in
expression (293), we find two equations

$\displaystyle\tilde{\boldsymbol{d}}$
$\displaystyle=|\boldsymbol{d}|\frac{\sqrt{3}}{2}\;n_{\mu}\left(\mathcal{R}_{5}(\phi_{0})\boldsymbol{\lambda}^{j=2}\right)_{\mu\mu^{\prime}}n_{\mu^{\prime}},$
(296) $\displaystyle\tilde{\boldsymbol{d}}$
$\displaystyle=|\boldsymbol{d}|\frac{\sqrt{3}}{2}\;\tilde{n}_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}\tilde{n}_{\mu^{\prime}},$
(297)

which have to be identical. Equating the two lines (296) and (297) leads to
the five equations

$\displaystyle 2\tilde{n}_{z}^{2}-\tilde{n}_{x}^{2}-\tilde{n}_{y}^{2}$
$\displaystyle=\frac{1+3\sqrt{5}}{4}n_{x}^{2}+\frac{1-3\sqrt{5}}{4}n_{y}^{2}-\frac{1}{2}n_{z}^{2},$
(298) $\displaystyle\tilde{n}_{x}^{2}-\tilde{n}_{y}^{2}$
$\displaystyle=\frac{-1+\sqrt{5}}{4}n_{x}^{2}+\frac{1+\sqrt{5}}{4}n_{y}^{2}-\frac{2\sqrt{5}}{4}n_{z}^{2},$
(299) $\displaystyle\tilde{n}_{x}\tilde{n}_{z}$ $\displaystyle=n_{x}n_{z},$
(300) $\displaystyle\tilde{n}_{y}\tilde{n}_{z}$ $\displaystyle=n_{y}n_{z},$
(301) $\displaystyle\tilde{n}_{x}\tilde{n}_{y}$ $\displaystyle=n_{x}n_{y},$
(302)

that need to be simultaneously satisfied by the four degrees of freedom
comprised in $\boldsymbol{n}$ and $\tilde{\boldsymbol{n}}$ (recall that each
unit vector has two independent degrees of freedom). The set of equations
(298)-(302) has a total of $10$ solutions. The first four solutions correspond
to the equal amplitude condition
$\left|n_{\mu}\right|=\left|\tilde{n}_{\mu}\right|=1/\sqrt{3}$, which
automatically satisfies Eqs. (298)-(299). Eqs. (300)-(302) are solved by
appropriately choosing the relative signs between the components of each
vector, yielding
$\boldsymbol{n}_{i}=\tilde{\boldsymbol{n}}_{i}=\mathcal{V}_{i}^{111}$ with
$\mathcal{V}^{111}$ given by Eq. (45) and $i=1,2,3,4$. Recall that changing
the sign of $\boldsymbol{n}_{i}$ or $\tilde{\boldsymbol{n}}_{i}$ does not
change the corresponding nematic order parameter $\boldsymbol{d}$, which
allows us to discard the solutions
$\boldsymbol{n}_{i}=-\tilde{\boldsymbol{n}}_{i}$.

The other six solutions are obtained upon setting one component $n_{\mu}=0$ to
zero. This necessarily requires $\tilde{n}_{\mu}=0$ in order to solve two of
the Eqs. (300)-(302). The remaining three equations can then be solved in a
straightforward way. For example, for $n_{z}=0$ we find the two solutions

$\displaystyle\boldsymbol{n}_{5}$
$\displaystyle=\left(\alpha_{-}^{(3)},\alpha_{+}^{(3)},0\right),$
$\displaystyle\tilde{\boldsymbol{n}}_{5}$
$\displaystyle=\left(\alpha_{+}^{(3)},\alpha_{-}^{(3)},0\right),$ (303)
$\displaystyle\boldsymbol{n}_{6}$
$\displaystyle=\left(\alpha_{-}^{(3)},-\alpha_{+}^{(3)},0\right),$
$\displaystyle\tilde{\boldsymbol{n}}_{6}$
$\displaystyle=\left(\alpha_{+}^{(3)},-\alpha_{-}^{(3)},0\right),$ (304)

with $\alpha_{\pm}^{(3)}=\frac{1}{\sqrt{6}}\sqrt{3\pm\sqrt{5}}$. The other
four solutions are found analogously and are presented in Eqs. (60), (44).
Note that the $\tilde{\boldsymbol{n}}$ directions in Eqs. (303)-(304) do not
carry any information, they are merely a mathematical construct used to find
the solution.

In the second parameter regime (ii), the cubic terms are minimized by
$\hat{\boldsymbol{d}}\cdot\hat{\boldsymbol{D}}^{H_{g},1}=-\hat{\boldsymbol{d}}\cdot\hat{\boldsymbol{D}}^{H_{g},1}=-\mathrm{sign}g_{1}$
or, equivalently,
$\cos\left(3\alpha\right)=-\cos\left(3\tilde{\alpha}\right)=-\mathrm{sign}g_{1}$.
Setting again $g_{1}<0$ for concreteness, we search for solutions with
$\alpha=0$ and $\tilde{\alpha}=\pi/3$. As a result, one obtains Eq. (296)
while the second equation (297) is replaced by

$\displaystyle\tilde{\boldsymbol{d}}$
$\displaystyle=-\frac{\sqrt{3}}{2}\,|\boldsymbol{d}|\;\tilde{m}_{\mu}\boldsymbol{\lambda}_{\mu\mu^{\prime}}^{j=2}\tilde{m}_{\mu^{\prime}}.$
(305)

Equating Eqs. (305) and (296) gives once again five equations

$\displaystyle\tilde{m}_{x}^{2}+\tilde{m}_{y}^{2}-2\tilde{m}_{z}^{2}$
$\displaystyle=\frac{1+3\sqrt{5}}{4}n_{x}^{2}+\frac{1-3\sqrt{5}}{4}n_{y}^{2}-\frac{1}{2}n_{z}^{2},$
(306) $\displaystyle-\tilde{m}_{x}^{2}+\tilde{m}_{y}^{2}$
$\displaystyle=\frac{-1+\sqrt{5}}{4}n_{x}^{2}+\frac{1+\sqrt{5}}{4}n_{y}^{2}-\frac{2\sqrt{5}}{4}n_{z}^{2},$
(307) $\displaystyle-\tilde{m}_{x}\tilde{m}_{z}$ $\displaystyle=n_{x}n_{z},$
(308) $\displaystyle-\tilde{m}_{y}\tilde{m}_{z}$ $\displaystyle=n_{y}n_{z},$
(309) $\displaystyle-\tilde{m}_{x}\tilde{m}_{y}$ $\displaystyle=n_{x}n_{y}.$
(310)

In contrast to the previous case, Eqs. (306)-(310) do not allow for solutions
with equal amplitudes due to the relative minus signs in Eqs. (308)-(310). The
six existing solutions are obtained by imposing a vanishing component
$n_{\mu}=0$, which directly implies $\tilde{m}_{\mu}=0$, leading to three
remaining equations. Those can be solved in a similar fashion as the previous
case. For instance, for $n_{z}=0$, we find the two solutions

$\displaystyle\boldsymbol{n}_{1}$
$\displaystyle=\left(\alpha_{+}^{(5)},-\alpha_{-}^{(5)},0\right),$
$\displaystyle\tilde{\boldsymbol{m}}_{1}$
$\displaystyle=\left(\alpha_{-}^{(5)},\alpha_{+}^{(5)},0\right),$ (311)
$\displaystyle\boldsymbol{n}_{2}$
$\displaystyle=\left(\alpha_{+}^{(5)},\alpha_{-}^{(5)},0\right),$
$\displaystyle\tilde{\boldsymbol{m}}_{2}$
$\displaystyle=\left(\alpha_{-}^{(5)},-\alpha_{+}^{(5)},0\right).$ (312)

A similar procedure can be applied for $n_{x}=0$ and $n_{y}=0$, resulting in
the solutions presented in Eqs. (62),(43).

##  Appendix D Landau expansion for the tetrahedral groups $\mathsf{T_{h}}$
and $\mathsf{T}$

In this Appendix, we derive the Landau expansion of the nematic doublet order
parameter $\boldsymbol{d}^{e}$ in the case of the cubic point groups
{$\mathsf{T_{h}}$, $\mathsf{T}$}. The main difference with respect to the
derivation presented in Sec. III.2, which applied to the cubic point groups
{$\mathsf{O_{h}}$, $\mathsf{O}$, $\mathsf{T_{d}}$}, is that here the nematic
doublet $\boldsymbol{d}^{e}$ transforms according to a complex IR. We
emphasize that, for the {$\mathsf{T_{h}}$, $\mathsf{T}$} point groups, the
nematic triplet order parameter $\boldsymbol{d}^{t}$ still transforms as a
real IR, so the results in Sec. III.2 apply directly to those groups as well.

To keep the notation transparent, we focus on the point group
$\mathsf{T_{h}}$; the results apply equally to the group $\mathsf{T}$. Within
$\mathsf{T_{h}}$, the complex combination
$\Delta^{e}=\left|\boldsymbol{d}^{e}\right|e^{\mathsf{i}\gamma_{e}}$
transforms as the $E_{g}$ IR while its complex conjugate $\bar{\Delta}^{e}$
transforms as the $\bar{E}_{g}$ IR, such that the two-component doublet
$\boldsymbol{d}^{e}$ transforms according to the representation
$\left(E_{g}\oplus\bar{E}_{g}\right)$. To derive the corresponding Landau
expansion, we first compute the decomposition of the symmetrized product for
each expansion order, following the procedure outlined in Appendix B:

$\displaystyle\big{[}\\!\otimes_{j=1}^{2}\left(E_{g}\oplus\bar{E}_{g}\right)\big{]}_{s}$
$\displaystyle=A_{g}\oplus\left(E_{g}\oplus\bar{E}_{g}\right),$ (313)
$\displaystyle\big{[}\\!\otimes_{j=1}^{3}\left(E_{g}\oplus\bar{E}_{g}\right)\big{]}_{s}$
$\displaystyle=2A_{g}\oplus\left(E_{g}\oplus\bar{E}_{g}\right),$ (314)
$\displaystyle\big{[}\\!\otimes_{j=1}^{4}\left(E_{g}\oplus\bar{E}_{g}\right)\big{]}_{s}$
$\displaystyle=A_{g}\oplus 2\left(E_{g}\oplus\bar{E}_{g}\right).$ (315)

We note the existence of one additional cubic invariant in this case, when
compared to the cases in which $\boldsymbol{d}^{e}$ transforms as a real IR,
see Eq. (81). To construct the four $A_{g}$ invariants in (313)-(315), we
first determine the bilinears associated with Eq. (313):

$\displaystyle D^{A_{g}}$ $\displaystyle=\left|\boldsymbol{d}^{e}\right|^{2},$
$\displaystyle D^{E_{g}}$
$\displaystyle=\left|\boldsymbol{d}^{e}\right|^{2}e^{-\mathsf{i}2\gamma_{e}},$
$\displaystyle D^{\bar{E}_{g}}$
$\displaystyle=\left|\boldsymbol{d}^{e}\right|^{2}e^{\mathsf{i}2\gamma_{e}}.$
(316)

Then, the four invariants can be written as
$\left|\boldsymbol{d}^{e}\right|^{2}$ , $\left|\boldsymbol{d}^{e}\right|^{4}$,
$\bar{\Delta}^{e}D^{E_{g}}$ and $\Delta^{e}D^{\bar{E}_{g}}$. From the latter
two, we construct two real combinations
$\left(\bar{\Delta}^{e}D^{E_{g}}+\Delta^{e}D^{\bar{E}_{g}}\right)$ and
$\mathsf{i}\left(\bar{\Delta}^{e}D^{E_{g}}-\Delta^{e}D^{\bar{E}_{g}}\right)$,
such that the Landau expansion becomes

$\displaystyle\mathcal{S}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}r_{0}\left|\boldsymbol{d}^{e}\right|^{2}+\left|\boldsymbol{d}^{e}\right|^{3}\Big{[}g^{\mathrm{c}}\cos\left(3\gamma_{e}\right)+g^{\mathrm{s}}\sin\left(3\gamma_{e}\right)\\!\Big{]}+u\left|\boldsymbol{d}^{e}\right|^{4}\\!\Big{\\}}.$
(317)

The presence of the second cubic term proportional to
$\sin\left(3\gamma_{e}\right)$ can be traced back to the fact that the cubic
groups {$\mathsf{T_{h}}$, $\mathsf{T}$} lack an axis of (proper or improper)
fourfold rotational symmetry compared to the cubic groups {$\mathsf{O_{h}}$,
$\mathsf{O}$, $\mathsf{T_{d}}$}. To proceed, it is instructive to rewrite the
cubic terms as

$\displaystyle
g^{\mathrm{c}}\cos\left(3\gamma_{e}\right)+g^{\mathrm{s}}\sin\left(3\gamma_{e}\right)$
$\displaystyle=g\,\cos\left(3\gamma_{e}-\delta_{0}\right),$ (318)

with

$\displaystyle g$
$\displaystyle=\mathrm{sign}\left(g^{\mathrm{c}}\right)\sqrt{(g^{\mathrm{c}})^{2}+(g^{\mathrm{s}})^{2}},$
$\displaystyle\delta_{0}$
$\displaystyle=\arctan\left(g^{\mathrm{s}}\big{/}g^{\mathrm{c}}\right).$ (319)

Then, the action (317) 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}-\delta_{0}\right)+u\left|\boldsymbol{d}^{e}\right|^{4}\Big{\\}},$
(320)

illustrating the emergence of an offset angle $\delta_{0}$ for the nematic
director that is non-universal, i.e. it is determined by the Landau
coefficients, which in turn depend on the microscopic model. As discussed
above, this is due to the lack of proper/improper fourfold rotational
symmetry. This Landau expansion thus behaves effectively as a $Z_{3}$-Potts
model whose threefold-degenerate ground-state directions,

$\displaystyle\gamma_{e}^{0}$
$\displaystyle=\frac{1}{3}\delta_{0}+\frac{m\pi}{3},$ (321)

are offset by an angle $\delta_{0}/3$. Here, $m=\left\\{1,3,5\right\\}$ for
$g>0$ and $m=\left\\{0,2,4\right\\}$ for $g<0$. In Table 1, we signal this
modified form of the $Z_{3}$-Potts model via $Z_{3}^{*}$.

To further visualize the effect of the offset angle, we rewrite the order
parameter
$\boldsymbol{d}^{e}=\left|\boldsymbol{d}^{e}\right|\left(\cos\gamma_{e},\sin\gamma_{e}\right)$
in the $(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation of Eq.
(24). To accomplish that, we define
$n_{e}=\mathrm{div}\left(\gamma_{e},\frac{\pi}{3}\right)$ and
$\tilde{\gamma}_{e}=\mathrm{mod}\left(\gamma_{e},\frac{\pi}{3}\right)=\gamma_{e}-n_{e}\frac{\pi}{3}\in[0,\frac{\pi}{3}]$
with $\mathrm{div}$, $\mathrm{mod}$ denoting integer division and modulo,
respectively. Since the definition in Eq. (319) implies
$\delta_{0}\in\left[-\pi/2,\,\pi/2\right]$, it follows that
$n_{e}=m-\left(\frac{1-\mathrm{sign}(\delta_{0})}{2}\right)$ and
$\tilde{\gamma}_{e}=\frac{\delta_{0}}{3}+\frac{\pi}{3}\left(\frac{1-\mathrm{sign}(\delta_{0})}{2}\right)\in\left[0,\,\frac{\pi}{3}\right]$.
Then, for even $n_{e}$, the $\boldsymbol{d}^{e}$ order parameter in the
$(\boldsymbol{n}\boldsymbol{m}\boldsymbol{l})$-representation of Eq. (24)
becomes

$\displaystyle\alpha$ $\displaystyle=\tilde{\gamma}_{e},$
$\displaystyle\boldsymbol{n}$
$\displaystyle=\boldsymbol{e}_{\tilde{n}_{e}+2},$
$\displaystyle\boldsymbol{m}$
$\displaystyle=\boldsymbol{e}_{\tilde{n}_{e}+1},$
$\displaystyle\boldsymbol{l}$ $\displaystyle=\boldsymbol{e}_{\tilde{n}_{e}},$
(322)

where $\tilde{n}_{e}=\frac{1}{2}n_{e}+1\in\\{1,2,3\\}$ and the summation in
(322) is understood as modulo $3$. Conversely, one obtains for odd $n_{e}$,

$\displaystyle\alpha$ $\displaystyle=\frac{\pi}{3}-\tilde{\gamma}_{e},$
$\displaystyle\boldsymbol{n}$ $\displaystyle=\boldsymbol{e}_{\tilde{n}_{e}},$
$\displaystyle\boldsymbol{m}$
$\displaystyle=\boldsymbol{e}_{\tilde{n}_{e}+1},$
$\displaystyle\boldsymbol{l}$
$\displaystyle=\boldsymbol{e}_{\tilde{n}_{e}+2},$ (323)

with $\tilde{n}_{e}=\frac{1}{2}n_{e}+\frac{1}{2}\in\\{1,2,3\\}$.
Interestingly, in this case of a $Z_{3}^{*}$-Potts transition, the offset
angle only affects $\alpha$ but leaves the nematic axes aligned with the
coordinate axes $\boldsymbol{e}_{1}$, $\boldsymbol{e}_{2}$, and
$\boldsymbol{e}_{3}$. Thus, within the point groups $\mathsf{T_{h}}$ and
$\mathsf{T}$, the nematic doublet state is generically biaxial with
$\alpha\neq\\{0,\pi/3\\}$. Only for $\gamma_{e}=n\frac{\pi}{3}$
[$n\in\mathbb{N}$], which is a symmetry-enforced condition in the point groups
{$\mathsf{O_{h}}$, $\mathsf{O}$, $\mathsf{T_{d}}$}, the state is uniaxial with
$\alpha=\\{0,\pi/3\\}$.

##  Appendix E Landau expansion for the hexagonal groups $\mathsf{C_{6h}}$,
$\mathsf{C_{3h}}$ and $\mathsf{C_{6}}$

In this Appendix, we derive the nematic Landau expansion for the hexagonal
point groups {$\mathsf{C_{6h}}$, $\mathsf{C_{3h}}$, $\mathsf{C_{6}}$}. In
contrast to the cases presented in Sec. IV.2, the nematic in-plane and out-of-
plane doublets $\boldsymbol{d}^{\mathrm{ip}}$ and
$\boldsymbol{d}^{\mathrm{op}}$ transform according to complex IRs. For
concreteness, we focus on the point group $\mathsf{C_{6h}}$, for which the
complex in-plane and out-of-plane nematic order parameters,
$\Delta^{\mathrm{ip}}=\left|\boldsymbol{d}^{\mathrm{ip}}\right|e^{\mathsf{i}\gamma_{\mathrm{ip}}}$
and
$\Delta^{\mathrm{op}}=\left|\boldsymbol{d}^{\mathrm{op}}\right|e^{\mathsf{i}\gamma_{\mathrm{op}}}$,
transform according to the IRs $E_{2g}$ and $E_{1g}$, respectively, whereas
their complex conjugates transform as $\bar{E}_{2g}$ and $\bar{E}_{1g}$. The
same results hold for the other two groups.

To derive the Landau expansion for the in-plane nematic doublet, we first
compute the decomposition of the symmetrized product for each expansion order,

$\displaystyle\big{[}\otimes_{j=1}^{2}\left(E_{2g}\oplus\bar{E}_{2g}\right)\big{]}_{s}$
$\displaystyle=A_{g}\oplus\left(E_{2g}\oplus\bar{E}_{2g}\right),$ (324)
$\displaystyle\big{[}\otimes_{j=1}^{3}\left(E_{2g}\oplus\bar{E}_{2g}\right)\big{]}_{s}$
$\displaystyle=2A_{g}\oplus\left(E_{2g}\oplus\bar{E}_{2g}\right),$ (325)
$\displaystyle\big{[}\otimes_{j=1}^{4}\left(E_{2g}\oplus\bar{E}_{2g}\right)\big{]}_{s}$
$\displaystyle=A_{g}\oplus 2\left(E_{2g}\oplus\bar{E}_{2g}\right).$ (326)

The four invariants are written in terms of the bilinears in Eq. (324),

$\displaystyle D^{A_{g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2},$ $\displaystyle
D^{E_{2g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}e^{-\mathsf{i}2\gamma_{\mathrm{ip}}},$
$\displaystyle D^{\bar{E}_{2g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}e^{\mathsf{i}2\gamma_{\mathrm{ip}}}.$
(327)

resulting in $\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}$,
$\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{4}$ ,
$\Delta^{\mathrm{ip}}D^{\bar{E}_{2g}}$ and
$\bar{\Delta}^{\mathrm{ip}}D^{E_{2g}}$. Thus, enforcing the invariants to be
real-valued, we obtain the nematic 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}}-\delta_{0}\right)+u\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{4}\Big{\\}},$
(328)

with an angular-shifted cosine term obtained from the relationship (318). The
Landau expansion (328) is that of a $Z_{3}$-clock model with an offset angle
$\delta_{0}$, whose minimization gives

$\displaystyle\gamma_{\mathrm{ip}}^{0}$
$\displaystyle=\frac{1}{3}\delta_{0}+\frac{\pi}{3}\left(\frac{1+\mathfrak{s}_{g}}{2}\right)+\frac{2\pi}{3}n,$
$\displaystyle n$ $\displaystyle=\left\\{0,1,2\right\\},$ (329)

where $\mathfrak{s}_{g}=\mathrm{sign}g_{\mathrm{ip}}$. In the nematic state,
the Fermi surface is similar to that shown in Fig. 5(a), but arbitrarily
rotated about the $k_{z}$-axis according to the offset angle $\delta_{0}$,
which is a Landau coefficient. To see this, we start from the parametrization
of $\boldsymbol{d}^{\mathrm{ip}}$ in Eq. (117) and employ the identities

$\displaystyle\left(\begin{smallmatrix}\cos\frac{\gamma_{\mathrm{ip}}^{0}}{2}\\\
\sin\frac{\gamma_{\mathrm{ip}}^{0}}{2}\end{smallmatrix}\right)$
$\displaystyle=\left(\text{-}1\right)^{n}R_{z}\left(\frac{\delta_{0}}{6}\right)\Bigg{[}\left(\begin{smallmatrix}\cos\tilde{\gamma}_{\mathrm{ip}}^{0}\\\
\text{-}\sin\tilde{\gamma}_{\mathrm{ip}}^{0}\end{smallmatrix}\right)\frac{1-\mathfrak{s}_{g}}{2}+\left(\begin{smallmatrix}\sin\tilde{\gamma}_{\mathrm{ip}}^{0}\\\
\cos\tilde{\gamma}_{\mathrm{ip}}^{0}\end{smallmatrix}\right)\frac{1+\mathfrak{s}_{g}}{2}\Bigg{]},$
(330)
$\displaystyle\left(\begin{smallmatrix}\text{-}\sin\frac{\gamma_{\mathrm{ip}}^{0}}{2}\\\
\cos\frac{\gamma_{\mathrm{ip}}^{0}}{2}\end{smallmatrix}\right)$
$\displaystyle=\left(\text{-}1\right)^{n}R_{z}\left(\frac{\delta_{0}}{6}\right)\Bigg{[}\left(\begin{smallmatrix}\sin\tilde{\gamma}_{\mathrm{ip}}^{0}\\\
\cos\tilde{\gamma}_{\mathrm{ip}}^{0}\end{smallmatrix}\right)\frac{1-\mathfrak{s}_{g}}{2}-\left(\begin{smallmatrix}\cos\tilde{\gamma}_{\mathrm{ip}}^{0}\\\
\text{-}\sin\tilde{\gamma}_{\mathrm{ip}}^{0}\end{smallmatrix}\right)\frac{1+\mathfrak{s}_{g}}{2}\Bigg{]},$
(331)

with
$\tilde{\gamma}_{\mathrm{ip}}^{0}=\frac{\pi}{3}\frac{1+\mathfrak{s}_{g}}{2}+\frac{2\pi}{3}n$
and the rotation matrix

$\displaystyle R_{z}\left(\delta\right)$
$\displaystyle=\left(\begin{smallmatrix}\cos\delta&-\sin\delta\\\
\sin\delta&\cos\delta\end{smallmatrix}\right).$ (332)

The relationships in Eqs. (330)-(331) explicitly show that the in-plane
nematic axes are rotated by an angle $\delta_{0}/6$ away from the hexagonal
high-symmetry directions listed in Eq. (127).

Repeating the same steps for the out-of-plane doublet order parameter
$\boldsymbol{d}^{\mathrm{op}}$, we obtain the symmetrized-product
decompositions

$\displaystyle\big{[}\otimes_{j=1}^{2}\left(E_{1g}\oplus\bar{E}_{1g}\right)\big{]}_{s}$
$\displaystyle=A_{g}\oplus\left(E_{2g}\oplus\bar{E}_{2g}\right),$ (333)
$\displaystyle\big{[}\otimes_{j=1}^{3}\left(E_{1g}\oplus\bar{E}_{1g}\right)\big{]}_{s}$
$\displaystyle=2B_{g}\oplus\left(E_{1g}\oplus\bar{E}_{1g}\right),$ (334)
$\displaystyle\big{[}\otimes_{j=1}^{4}\left(E_{1g}\oplus\bar{E}_{1g}\right)\big{]}_{s}$
$\displaystyle=A_{g}\oplus 2\left(E_{2g}\oplus\bar{E}_{2g}\right),$ (335)
$\displaystyle\big{[}\otimes_{j=1}^{5}\left(E_{1g}\oplus\bar{E}_{1g}\right)\big{]}_{s}$
$\displaystyle=2B_{g}\oplus 2\left(E_{1g}\oplus\bar{E}_{1g}\right),$ (336)
$\displaystyle\big{[}\otimes_{j=1}^{6}\left(E_{1g}\oplus\bar{E}_{1g}\right)\big{]}_{s}$
$\displaystyle=3A_{g}\oplus 2\left(E_{2g}\oplus\bar{E}_{2g}\right),$ (337)

with the bilinears

$\displaystyle D^{A_{g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2},$ $\displaystyle
D^{E_{2g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}e^{\mathsf{i}2\gamma_{\mathrm{op}}},$
$\displaystyle D^{\bar{E}_{2g}}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}e^{-\mathsf{i}2\gamma_{\mathrm{op}}}.$
(338)

Compared to the case of $\boldsymbol{d}^{\mathrm{ip}}$, here the cubic terms
$\Delta^{\mathrm{op}}D^{E_{2g}}$ and
$\bar{\Delta}^{\mathrm{op}}D^{\bar{E}_{2g}}$ do not transform as $A_{g}$ but
as $B_{g}$. Thus, the five invariants in Eqs. (333)-(337) are
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}$,
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}$,
$\left|\boldsymbol{d}^{\mathrm{op}}\right|^{6}$,
$\left(\Delta^{\mathrm{op}}D^{E_{2g}}\right)^{2}$ and
$\big{(}\bar{\Delta}^{\mathrm{op}}D^{\bar{E}_{2g}}\big{)}^{2}$, yielding the
Landau expansion

$\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}$
$\displaystyle+v_{6}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{6}\cos\left(6\gamma_{\mathrm{op}}-\delta_{0}\right)\Big{\\}},$
(339)

with an offset angle $\delta_{0}$ due to the relation (318). Minimization
leads to

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

which, from Eq. (119), corresponds to a rotation of the nematic axes by
$\delta_{0}/6$ about the $k_{z}$-axis.

##  Appendix F Landau expansion for the trigonal groups $\mathsf{C_{3}}$ and
$\mathsf{S_{6}}$

This Appendix presents the derivation of the Landau expansion for the trigonal
groups {$\mathsf{S_{6}}$, $\mathsf{C_{3}}$}, which lack the in-plane symmetry
directions $\mathcal{V}_{1}^{\mathrm{hex}}$, $\mathcal{V}_{2}^{\mathrm{hex}}$
in Eq. (127), thus complementing the analysis shown in Sec. IV.3 for the
trigonal groups {$\mathsf{D_{3}}$, $\mathsf{D_{3d}}$, $\mathsf{C_{3v}}$}. We
consider the point group $\mathsf{C_{3}}$ for concreteness; in this case, the
combinations
$\Delta^{\mathrm{ip}}=\left|\boldsymbol{d}^{\mathrm{ip}}\right|e^{\mathsf{i}\gamma_{\mathrm{ip}}}$
and
$\Delta^{\mathrm{op}}=\left|\boldsymbol{d}^{\mathrm{op}}\right|e^{\mathsf{i}\gamma_{\mathrm{op}}}$
transform according to the complex IR $E$, whereas their complex conjugates
$\bar{\Delta}^{\mathrm{ip}}$, $\bar{\Delta}^{\mathrm{op}}$ transform according
to $\bar{E}$. Since we are interested in the Landau expansion of the total
nematic order parameter
$\boldsymbol{d}^{E}=\big{(}\boldsymbol{d}^{\mathrm{ip}},\boldsymbol{d}^{\mathrm{op}}\big{)}$,
we need the symmetrized product decomposition of
$\left(E\oplus\bar{E}\right)\oplus\left(E\oplus\bar{E}\right)$:

$\displaystyle\big{[}\otimes_{j=1}^{2}\left(\left(E\oplus\bar{E}\right)\oplus\left(E\oplus\bar{E}\right)\right)\big{]}_{s}$
$\displaystyle=4A\oplus 3\left(E\oplus\bar{E}\right),$ (341)
$\displaystyle\big{[}\otimes_{j=1}^{3}\left(\left(E\oplus\bar{E}\right)\oplus\left(E\oplus\bar{E}\right)\right)\big{]}_{s}$
$\displaystyle=8A\oplus 6\left(E\oplus\bar{E}\right),$ (342)
$\displaystyle\big{[}\otimes_{j=1}^{4}\left(\left(E\oplus\bar{E}\right)\oplus\left(E\oplus\bar{E}\right)\right)\big{]}_{s}$
$\displaystyle=9A\oplus 13\left(E\oplus\bar{E}\right).$ (343)

The bilinear combinations obtained from Eq. (341) are

$\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|e^{\mathsf{i}\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{ip}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}e^{-\mathsf{i}2\gamma_{\mathrm{ip}}},$
$\displaystyle D_{\mathrm{op}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}e^{-\mathsf{i}2\gamma_{\mathrm{op}}},$
$\displaystyle D_{\mathrm{io}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|e^{-\mathsf{i}\left(\gamma_{\mathrm{ip}}+\gamma_{\mathrm{op}}\right)},$
(344)

as well as the complex conjugates of the bilinears that are not real-valued.
The $8$ cubic invariants transforming as $A$ are given by

$\displaystyle\bar{\Delta}^{\mathrm{ip}}D_{\mathrm{ip}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{3}e^{-\mathsf{i}3\gamma_{\mathrm{ip}}},$
$\displaystyle\bar{\Delta}^{\mathrm{ip}}D_{\mathrm{io}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|e^{-\mathsf{i}\left(2\gamma_{\mathrm{ip}}+\gamma_{\mathrm{op}}\right)},$
$\displaystyle\bar{\Delta}^{\mathrm{op}}D_{\mathrm{op}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{3}e^{-\mathsf{i}3\gamma_{\mathrm{op}}},$
$\displaystyle\bar{\Delta}^{\mathrm{op}}D_{\mathrm{io}}^{E}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}e^{-\mathsf{i}\left(\gamma_{\mathrm{ip}}+2\gamma_{\mathrm{op}}\right)},$
(345)

plus their complex conjugates. Similarly, the $9$ quartic invariants are

$\displaystyle D_{\mathrm{ip}}^{A}D_{\mathrm{ip}}^{A}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{4},$ $\displaystyle
D_{\mathrm{ip}}^{A}D_{\mathrm{io}}^{A}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{3}\left|\boldsymbol{d}^{\mathrm{op}}\right|e^{\mathsf{i}\left(\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}\right)},$
$\displaystyle D_{\mathrm{op}}^{A}D_{\mathrm{op}}^{A}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4},$ $\displaystyle
D_{\mathrm{op}}^{A}D_{\mathrm{io}}^{A}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|^{3}e^{\mathsf{i}\left(\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}\right)},$
$\displaystyle D_{\mathrm{ip}}^{A}D_{\mathrm{op}}^{A}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2},$
$\displaystyle D_{\mathrm{io}}^{A}D_{\mathrm{io}}^{A}$
$\displaystyle=\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}e^{\mathsf{i}\left(2\gamma_{\mathrm{ip}}-2\gamma_{\mathrm{op}}\right)},$
(346)

plus the complex conjugates of the terms on the right column. We can now write
down the Landau expansion in terms of the real-valued combinations of these
invariants. The resulting action becomes
$\mathcal{S}=\mathcal{S}_{2}+\mathcal{S}_{3}+\mathcal{S}_{4}$, with

$\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}1}\left(\boldsymbol{d}^{\mathrm{ip}}\cdot\boldsymbol{d}^{\mathrm{op}}\right)$
$\displaystyle+r_{\mathrm{io}2}(\boldsymbol{d}^{\mathrm{ip}})^{T}\left(-\mathsf{i}\sigma^{y}\right)\boldsymbol{d}^{\mathrm{op}}\Big{\\}},$
(347) $\displaystyle\mathcal{S}_{3}$
$\displaystyle=\int_{\mathsf{x}}\Big{\\{}\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{3}\left[g_{\mathrm{ip}}^{\mathrm{c}}\mathsf{c}_{3\gamma_{\mathrm{ip}}}+g_{\mathrm{ip}}^{\mathrm{s}}\mathsf{s}_{3\gamma_{\mathrm{ip}}}\right]+\left|\boldsymbol{d}^{\mathrm{op}}\right|^{3}\left[g_{\mathrm{op}}^{\mathrm{c}}\mathsf{c}_{3\gamma_{\mathrm{op}}}+g_{\mathrm{op}}^{\mathrm{s}}\mathsf{s}_{3\gamma_{\mathrm{op}}}\right]$
$\displaystyle+\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\left[g_{2}^{\mathrm{c}}\mathsf{c}_{\gamma_{\mathrm{ip}}+2\gamma_{\mathrm{op}}}+g_{2}^{\mathrm{s}}\mathsf{s}_{\gamma_{\mathrm{ip}}+2\gamma_{\mathrm{op}}}\right]$
$\displaystyle+\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|\left[g_{1}^{\mathrm{c}}\mathsf{c}_{2\gamma_{\mathrm{ip}}+\gamma_{\mathrm{op}}}+g_{1}^{\mathrm{s}}\mathsf{s}_{2\gamma_{\mathrm{ip}}+\gamma_{\mathrm{op}}}\right]\Big{\\}},$
(348) $\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|^{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}\left[u_{\mathrm{io}}^{\mathrm{c}}\mathsf{c}_{2\gamma_{\mathrm{ip}}-2\gamma_{\mathrm{op}}}+u_{\mathrm{io}}^{\mathrm{s}}\mathsf{s}_{2\gamma_{\mathrm{ip}}-2\gamma_{\mathrm{op}}}\right]$
$\displaystyle+u_{\mathrm{io}}^{0}\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{2}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{2}+\left[u_{2}^{\mathrm{c}}\mathsf{c}_{\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}}+u_{2}^{\mathrm{s}}\mathsf{s}_{\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}}\right]\left|\boldsymbol{d}^{\mathrm{ip}}\right|\left|\boldsymbol{d}^{\mathrm{op}}\right|^{3}$
$\displaystyle+u_{\mathrm{op}}\left|\boldsymbol{d}^{\mathrm{op}}\right|^{4}+\left[u_{1}^{\mathrm{c}}\mathsf{c}_{\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}}+u_{1}^{\mathrm{s}}\mathsf{s}_{\gamma_{\mathrm{ip}}-\gamma_{\mathrm{op}}}\right]\left|\boldsymbol{d}^{\mathrm{ip}}\right|^{3}\left|\boldsymbol{d}^{\mathrm{op}}\right|\Big{\\}},$
(349)

with $23$ invariants. For brevity, we defined
$\mathsf{c}_{\gamma}\equiv\cos\gamma$ and
$\mathsf{s}_{\gamma}\equiv\sin\gamma$; the Pauli matrix $i\sigma^{y}$ in the