Broken mirror symmetry, incommensurate spin correlations, and $B_{2g}$ nematic order in iron pnictides

TL;DR

This paper investigates the origin of $B_{2g}$ nematic order in hole-doped iron pnictides, linking it to incommensurate magnetic fluctuations and proposing a microscopic mechanism supported by theoretical modeling.

Contribution

It introduces a classification of nematic orders by mirror symmetries, constructs a Ginzburg-Landau theory connecting magnetic fluctuations to nematicity, and identifies incommensurate $(q,q)$ fluctuations as the driver of $B_{2g}$ nematic order.

Findings

Incommensurate $(q,q)$ magnetic fluctuations underlie $B_{2g}$ nematic order.

Microscopic calculations support the link between spin correlations and nematicity.

The proposed mechanism aligns with existing experimental evidence and suggests new experimental tests.

Abstract

Motivated by growing indications for a distinct form of nematic correlations in the extremely hole doped iron pnictide compounds $A$Fe$_2$As$_2$ ($A$=K,Rb,Cs), we consider spin-driven nematic order in the general case of incommensurate magnetic fluctuations. We classify the nematic order parameters by broken mirror symmetries of the tetragonal $D_{4h}$ point group, and use this scheme to construct a general Ginzburg-Landau theory that links the nematic order to spatial pattern of magnetic fluctuations. Our analysis points to incommensurate magnetic fluctuations of wavevector $(q,q)$ as underlying a $B_{2g}$ nematic order in $A$Fe$_2$As$_2$. We substantiate this idea by microscopic calculations of the nematic order based on 3-sublattice $(2\pi/3,2\pi/3)$ spin correlations in an extended bilinear-biquadratic Heisenberg model. We summarize the existing evidence in support of the proposed mechanism, suggest ways to further test it in future experiments, and discuss the implications of our results for iron-based high temperature superconductivity.

Tables

Table 1. Table 1: Symmetry classification of nematicity in FeSCs based on the broken mirror symmetries of the D 4 ​ h subscript 𝐷 4 ℎ D_{4h} group.

Nematicity	${\sigma }_{x/y}$	${\sigma }_{d/d^{\prime }}$	$C_4={\sigma }_{x/y}\times {\sigma }_{d/d^{\prime }}$
$B_{1g}$	1	-1	-1
$B_{2g}$	-1	1	-1
$A_{2g}$	-1	-1	1

Equations

$${\sigma }_{B1}$$

$${\sigma }_{B2}$$

$${\sigma }_{A2}$$

$$f_{B2}$$

$$f$$

$$f_{B1}$$

$$f_{12}$$

$$H$$

$$S$$

$$f\sim r_{\sigma }({\sigma }_{B1}^2+{\sigma }_{B2}^2+{\sigma }_{A2}^2)+b_{\sigma }{\sigma }_{B1}{\sigma }_{B2}{\sigma }_{A2}+O({\sigma }^4).$$

$$f$$

$$f_{B2}$$

$$f_{B1}$$

$$f_{12}$$

$$\frac{\partial f}{\partial |m_1|}$$

$$\frac{\partial f}{\partial |m_2|}$$

$$\frac{\partial f}{\partial |m_3|}$$

$$\frac{\partial f}{\partial |m_4|}$$

$$S=$$

$$S_2=$$

$$S_4=$$

$$-v_2\left[(m_1^2-m_2^2-m_3^2+m_4^2{)}^2+(m_1^2-m_2^2+m_3^2-m_4^2{)}^2+(m_1^2+m_2^2-m_3^2-m_4^2{)}^2\right]\},$$

$$-r\sum _i{m}_i^2-\frac{v_1}{N}(\sum _i{m}_i^2{)}^2$$

$$v_2(m_1^2-m_2^2-m_3^2+m_4^2{)}^2$$

$$v_2(m_1^2-m_2^2+m_3^2-m_4^2{)}^2$$

$$v_2(m_1^2+m_2^2-m_3^2-m_4^2{)}^2$$

$$S=$$

$$+{\sigma }_{B1}(m_1^2-m_2^2-m_3^2+m_4^2)+{\sigma }_{B2}(m_1^2-m_2^2+m_3^2-m_4^2)+{\sigma }_{A2}(m_1^2+m_2^2-m_3^2-m_4^2)\}$$

$$+\sum _{\mathbf{k}}c{\mathbf{k}}^2\sum _i{m}_i(\mathbf{k}{)}^2.$$

$$f=$$

$$+(\lambda -{\sigma }_{B1}+{\sigma }_{B2}-{\sigma }_{A2})m_3^2+(\lambda +{\sigma }_{B1}-{\sigma }_{B2}-{\sigma }_{A2})m_4^2+g(\lambda ,{\sigma }_{B1},{\sigma }_{B2},{\sigma }_{A2})$$

$$g(\lambda ,{\sigma }_{B1},{\sigma }_{B2},{\sigma }_{A2})=$$

$$+\ln (\lambda -{\sigma }_{B1}+{\sigma }_{B2}-{\sigma }_{A2}+c{\mathbf{k}}^2)+\ln (\lambda +{\sigma }_{B1}-{\sigma }_{B2}-{\sigma }_{A2}+c{\mathbf{k}}^2)$$

$$\frac{\partial f}{\partial \lambda }$$

$$\frac{\partial f}{\partial {\sigma }_{B1}}$$

$$\frac{\partial f}{\partial {\sigma }_{B2}}$$

$$\frac{\partial f}{\partial {\sigma }_{A2}}$$

$$\frac{\partial f}{\partial m_1}$$

$$\frac{\partial f}{\partial m_2}$$

$$\frac{\partial f}{\partial m_3}$$

$$\frac{\partial f}{\partial m_4}$$

$$\frac{r-\lambda }{v_1}+\frac{{\sigma }_{B1}+{\sigma }_{B2}+{\sigma }_{A2}}{v_2}+4{\int }_0^{\Lambda }\frac{d^{2}k}{(2\pi {)}^2}\frac{1}{\lambda +{\sigma }_{B1}+{\sigma }_{B2}+{\sigma }_{A2}+c{k}^2}=0,$$

$$\frac{r-\lambda }{v_1}+\frac{-{\sigma }_{B1}-{\sigma }_{B2}+{\sigma }_{A2}}{v_2}+4{\int }_0^{\Lambda }\frac{d^{2}k}{(2\pi {)}^2}\frac{1}{\lambda -{\sigma }_{B1}-{\sigma }_{B2}+{\sigma }_{A2}+c{k}^2}=0,$$

Full text

*Introduction. * Strongly correlated systems often involve multiple building blocks for their macroscopic properties. Iron-based superconductors (FeSCs)  Kamihara2008 ; Johnston ; Dai2015 ; NatRevMat:2016 ; Hirschfeld2016 ; FWang-science2011 provide a prototype example. Typically, the phase diagram contains an antiferromagnetic (AFM) order, pointing to the role of spins. It prominently features a nematic order, which may be driven by spin or other degrees of freedom. Understanding its origin and the associated fluctuations will likely shed light on the mechanism of high temperature superconductivity.

In the most common iron pnictides, an electronic nematic order  MYi:2011 ; IFisher:2012 accompanies AF order of wave vector $(\pi,0)$ [Fig. 1(a)]. It lowers the $C_{4}$ rotational symmetry of the tetragonal lattice to $C_{2}$ by making the tetragonal $a$ and $b$ axes inequivalent. According to the tetragonal lattice notation, the nematic order has a $B_{1g}$ symmetry. However, nematic order in the FeSCs has considerable variations. The bulk FeSe, for example, has a $B_{1g}$ nematic order which is not accompanied by any AF order Zhao:2015 . A great deal of efforts have recently devoted to study this nematic order of FeSe.

A new surprise has emerged from heavily hole doped (Rb,Cs)Fe2As2 Feng:2018 ; Shibauchi:2018 ; Wu:2016 . Recent scanning tunneling microscopy (STM) measurements observe a two-fold symmetric quasiparticle interference (QPI) pattern about the two diagonal directions of Fe lattice Feng:2018 . Elastoresistance data also reveal an anisotropy along this direction Shibauchi:2018 . Both experiments evidence that the nematic order here has a $B_{2g}$ symmetry, which corresponds to a pattern that is rotated from its $B_{1g}$ counterpart by $45^{\circ}$. Equally important, for a range of doping and temperature in RbxBa1-xFe2As2 ($x$ near $0.8$), the $B_{2g}$ and $B_{1g}$ nematic channels are nearly degenerate Shibauchi:2018 .

An important question is whether a universal origin exists for the variety of nematic orders. One candidate mechanism attributes the $B_{1g}$ nematicity to an Ising order that is constructed from AFM or antiferroquadrupolar (AFQ) fluctuations Dai_PNAS:2009 ; FangKivelson:2008 ; XuMullerSachdev:2008 ; Yu:2015 at wave vector $(\pi,0)$ or $(0,\pi)$. To consider the possibility of the $B_{2g}$ nematicity in this light, we are motivated to explore more general types of magnetic fluctuations. Indeed, the spin excitations of KFe2As2 [Fig. 1(c)] are peaked near wave vector $(q,0)$ with $q\approx 2\pi/3$ at low energies, and with increasing energy the wave vector saturates near $(q,q)$. Compared to BaFe2As2 Harriger:2011 and K0.5Ba0.5Fe2As2 (see Fig. S2 of SM SM ), the $(q,q)$ spin excitations occupy a large spectral weight in KFe2As2. In addition, AFe2As2 has been evidenced to move towards an AFM quantum critical point as one goes from A=K to A=(Rb,Cs) Eilers:2016 , making it likely that the $(q,q)$ spin excitations further soften and grow in spectral weight for the (Rb,Cs) cases.

In this manuscript, we are thus motivated to study the role of incommensurate magnetic fluctuations on the nematicity. To this end, we consider the electrons residing on the tetragonal lattice and classify the nematic orders in terms of a broken mirror symmetry to $B_{1g}$, $B_{2g}$, and $A_{2g}$. Building on this symmetry analysis, we propose a general Ginzburg-Landau theory and connect the various nematic orders with the underlying incommensurate magnetic fluctuations. This allows for a unified understanding for the nematicity in FeSCs. In particular, we demonstrate that incommensurate $(q,q)$ and $(-q,q)$ magnetic fluctuations lead to a $B_{2g}$ Ising nematic order. This result is further supported by calculations on a microscopic bilinear-biquadratic Heisenberg model, which find a $B_{2g}$ nematic order from 3-sublattice $(2\pi/3,2\pi/3)$ AFM correlations [Fig. 1(b)]. Finally, through the formulation of broken mirror symmetry, we advance a robust mechanism for a quasi-degeneracy between several nematic channels.

*Classification of nematicity. * The nematic order of interest breaks a $Z_{2}$ symmetry, and is characterized by an Ising variable or a scalar order parameter. It can be classified according to the one-dimensional (1D) irreducible representations of the tetragonal point group ($D_{4h}$). Since inversion symmetry is preserved, a nematic order should transform as $B_{1g}$, $B_{2g}$, or $A_{2g}$. Each of them is uniquely determined by examining its transformation under the mirror symmetries $\sigma_{x/y}$ and $\sigma_{d/d^{\prime}}$ [see Table 1 and Fig. 2(a)-(c)]. The usual $B_{1g}$ nematic order breaks the mirror plane passing through the diagonal directions ($\sigma_{d/d^{\prime}}$), but preserves the one through axes ($\sigma_{x/y}$). For the $B_{2g}$ nematic order, the roles of the two mirror planes are reversed. Finally, the $A_{2g}$ nematic order breaks both mirror symmetries $\sigma_{x/y}$ and $\sigma_{d/d^{\prime}}$ but preserves their product, which is the $C_{4}$ symmetry; it qualifies as a nematic state because the $C_{2}$ symmetry about either the $x$ or $y$ axis is broken.

*Construction of the Ginzburg-Landau theory. * We are led to a Ginzburg-Landau theory for the nematicity. Consider an incommensurate magnetic moment $\mathbf{m}_{1}=\mathbf{m}(q_{1},q_{2})$ with a generic wave vector $(q_{1},q_{2})$ and other three moments related by mirror symmetries, $\mathbf{m}_{2}=\mathbf{m}(-q_{2},q_{1})$, $\mathbf{m}_{3}=\mathbf{m}(q_{2},q_{1})$, and $\mathbf{m}_{4}=\mathbf{m}(-q_{1},q_{2})$. The Ising-nematic parameters are conventionally defined within each plaquette in real space to be (see Fig. 2(d)),

$\sigma_{B1}=(\mathbf{S}_{A}-\mathbf{S}_{D})\cdot(\mathbf{S}_{B}-\mathbf{S}_{C})$, and $\sigma_{B2}=\mathbf{S}_{A}\cdot\mathbf{S}_{D}-\mathbf{S}_{B}\cdot\mathbf{S}_{C}$. In momentum space, we can write

[TABLE]

Since it preserves $\sigma_{d/d^{\prime}}$, the $B_{2g}$ nematic order is naturally connected to magnetic moments $\mathbf{m}_{1}=\mathbf{m}(q,q)$ and $\mathbf{m}_{2}=\mathbf{m}(-q,q)$ that have the same symmetry [Fig. 2(b)]. We can then construct an effective Landau free energy as follows:

[TABLE]

This construction parallels that for $B_{1g}$ nematicity Yu:2017 . The ground state phase diagram [Fig. 2(e)] has an incommensurate AFM order at either $(q,q)$ or $(-q,q)$ when $u_{2}>0$ and $u_{2}>u_{3}$. Since $\mathbf{m}_{1}^{2}-\mathbf{m}_{2}^{2}\propto\sigma_{B2}$, this phase supports a $B_{2g}$ nematic order at finite temperature. There are two additional incommensurate double-Q phases, with $\mathbf{m}_{1}\parallel\mathbf{m}_{2}$ and $\mathbf{m}_{1}\perp\mathbf{m}_{2}$, respectively. They are analogies of the two double-Q AFM phases in the $B_{1g}$ case  Giovannetti_NC:2011 ; Yu:2017 , and are expected to have enhanced $B_{2g}$ nematic susceptibility.

We can construct a general free energy for both the $B_{1g}$ and $B_{2g}$ nematicity in terms of the relevant $(q,0)/(0,q)$ and $(q,q)/(-q,q)$ magnetic moments.

[TABLE]

where $\mathbf{m}_{3}=\mathbf{m}(q,0)$ and $\mathbf{m}_{4}=\mathbf{m}(0,q)$. The phase diagram is even richer (see SM SM ), containing single-Q AFM states with either $\mathbf{m}_{i}$ being ordered, which supports either $B_{1g}$ or $B_{2g}$ nematic order, and several double-Q AFM states with $C_{4}$ symmetry. The states with ordered moments $\mathbf{m}_{1/2}$ and $\mathbf{m}_{3/4}$ either are separated by a bicritical point or coexist, depending on the model parameters (see SM SM ).

*Bilinear-biquadratic Heisenberg model. * We now turn to a microscopic model. A bilinear-biquadratic Heisenberg model has successfully explained the $B_{1g}$ nematicity in iron pnictides and iron selenide Yu:2015 . Here we reexamine this model and explore the phase diagram. The Hamiltonian reads as

[TABLE]

where $\delta=1,2,3$, and the summation is up to the 3rd-nearest neighbors. We set $J_{1}=1$ as the energy unit. The frustrating interactions cause a rich phase diagram even in the classical spin limit (see SM). Fig. 3(a), shows the ground-state phase diagram for $J_{3}=0.1$, $K_{2}=K_{3}=0$ and varying $J_{2}$ and $K_{1}$. A $(\pi,0)/(0,\pi)$ AFM occurs when $J_{2}>J_{1}/2$ and $K_{1}<0$. A double-Q $C_{4}$ AFM state with $\mathbf{m}(\pi,0)\perp\mathbf{m}(0,\pi)$ is stabilized when $K_{1}>0$. We find that increasing $K_{1}$ while decreasing $J_{2}$ stabilizes a $(\pi,q)$ and further a $(q,q)$ AFM state. Here, the $(q,q)$ AFM state is stabilized due to the competition of $J_{1}$ and $K_{1}$, which is different from what happens in the classical $J_{1}-J_{2}-J_{3}$ model Moreo:1990 . The incommensurate $(q,q)$ state does support a $B_{2g}$ nematic order below a transition temperature as shown in the Monte Carlo result in Fig. 3(b).

DMRG study on the quantum $S=1$ model. We have in addition investigated the $S=1$ bilinear-biquadratic Heisenberg model, Eq.(8), by the density matrix renormalization group (DMRG) method. Including quantum fluctuations makes the phase diagram even richer, with several magnetic and quadrupolar phases, as well as a nematic spin liquid 2017arXiv171106523H . We find evidence for a robust 3-sublattice $(2\pi/3,2\pi/3)$ AFM phase, signaled by a clear peak at momentum $(2\pi/3,2\pi/3)$ of the spin structure factor [Fig. 4, inset]. This phase can be stabilized for $J_{2}=J_{3}=K_{2}=0$, $K_{1}=0.8$, and $-0.4\lesssim K_{3}\lesssim-0.1$, a parameter regime close to that in the classical model. This phase supports a $B_{2g}$ nematic order. As shown in Fig. 4 main panel, the $B_{2g}$ nematic order, $\sigma_{B2}$, scales to a nonzero value in the thermodynamic limit.

Degeneracy of nematic channels. The advantage of the mirror symmetry formulation is even clearer for the generic $q_{1}\neq q_{2}$ case. There are four moments, $\mathbf{m}_{1}=\mathbf{m}(q_{1},q_{2})$ and its mirror-symmetry related $\mathbf{m}_{2}$, $\mathbf{m}_{3}$, $\mathbf{m}_{4}$, defined earlier. The Ginzburg-Landau action with $D_{4h}$ symmetry reads

[TABLE]

where $\mathbf{m}_{i}(\mathbf{k})=\mathbf{m}(\mathbf{q}_{i}+\mathbf{k})$ with $\mathbf{q}_{i}=(\pm q_{1},\pm q_{2})$. A Hubbard-Stratonovich transformation (see SM SM ) yields

[TABLE]

Here, $r_{\sigma}$ comes from contributions of magnetic fluctuations and is dictated by symmetry to be identical in each nematic channel. When $r_{\sigma}>0$, the nematic fluctuations are exactly degenerate among the three channels, and when $r_{\sigma}<0$ a nematic order arises as shown in a large-$N$ calculation in SM SM . In the nematic phase, either the degeneracy is lifted by spontaneous ordering to one nematic channel, or the ordering takes place in all three channels with $\sigma_{B1}=\sigma_{B2}=\sigma_{A2}$. The latter case is due to the cubic term of $\sigma_{i}$ in Eq. (10), which reflects the discrete $D_{4h}$ symmetry. Note that the nature of the nematic transition is very different from the magnetic ordering in a Heisenberg or a XY model, where the spontaneous symmetry breaking can take place along any direction.

Recent elastoresistance measurement indeed reveals a quasi-degeneracy between the $B_{1g}$ and $B_{2g}$ nematic fluctuations in the intermediate hole doping regime of iron pnictides Shibauchi:2018 . Neutron scattering measurements Horigane:2016 show that, upon hole doping, enhanced incommensurate $(q_{1},q_{2})$ fluctuations appear in the low-energy spin excitation spectrum. Thus, this quasi-degeneracy is well understood in our theory.

When the $(q_{1},q_{2})$ type magnetic fluctuations couple to $(q,0)$ or $(q,q)$ fluctuations, the degeneracy among the three nematic channels can be lifted to a degree. In real materials, magnetic fluctuations couple to other degrees of freedom, such as orbital and lattice, which may also help break the exact degeneracy of the three nematic channels to stabilize a particular type of nematic order Dagotto:2016 ; Kontani:2018 . Nonetheless, our formulations reveals that the spin-driven nematicity naturally accounts for the observed quasi-degeneracy between the $B_{1g}$ and $B_{2g}$ fluctuations. It would also be interesting to explore the possibility of an $A_{2g}$ nematicity in FeSCs.

*Discussions and Conclusions. * We now note on several points. First, the proposed mechanism for a $B_{2g}$ nematicity well accounts for the observations by recent STM, elastoresistance, and NMR measurements in heavily hole doped iron pnictides.Feng:2018 ; Shibauchi:2018 ; Wu:2016 In our analysis the $B_{2g}$ nematic order is associated with the $(q,q)$-type incommensurate magnetic fluctuations, which are a large part of the spin spectral weight in KFe2As2 Horigane:2016 . The 3-sublattice $(2\pi/3,2\pi/3)$ AFM spin order is consistent with that part of the fluctuation spectrum [Fig. 1(c)]. Thermodynamic measurements have suggested that (Rb,Cs) replacement for K drives the system toward an AFM quantum critical point Eilers:2016 . It is thus likely that the $(q,q)$ AFM fluctuations will be enhanced in the (Rb,Cs) cases, thereby strengthening the $B_{2g}$ correlations. Inelastic neutron scattering measurements in (Rb,Cs)Fe2As2 are called for. We note in passing that the phase diagram of the bilinear-biquadratic model also contains a 3-sublattice $(2\pi/3,2\pi/3)$ AFQ order and a double-stripe $(\pi/2,\pi/2)$ AFM order, either of which may support the $B_{2g}$ nematicity Mila:2012 ; Yu:2015 ; Dagotto:2017 ; Lai:2016 ; Zhang_Fernandes:2017 .

Second, the softening of $(q,q)$ magnetic fluctuations with hole doping suggests reduced $J_{2}$ value from BaFe2As2, which drives the system from the $(\pi,0)$ to $(q,q)$ AFM order, as shown in Fig. 3(a). The strong electron correlations in (K,Rb,Cs)Fe2As2 make ab initio estimates of $J$’s and $K$’s difficult. Still, the $(q,q)$ AFM order is relevant for (K,Rb,Cs)Fe2As2, which has $N=5.5$ and corresponds to the strongly hole-doped counterpart of the $N=6$ BaFe2As2. The schematic zero-temperature phase diagram (Fig. 5)  Eilers:2016 ; Yu_COSSMS:2013 illustrates that two types of antiferromagnetic orders are respectively associated with the $N=6$ and $N=5-5.5$ regimes. Given that the half-filled $N=5$ case is expected to have a commensurate $(\pi,\pi)$ AFM order, it is natural for the $N=5.5$ to develop the $(q,q)$ AFM order. In this sense, the $(q,q)$ AFM order implicated by our work elucidates the microscopic physics of the FeSCs over an extended doping range.

Third, upon doping alkaline ions, experiments suggest that the low-temperature electronic states may evolve from a $B_{1g}$ nematic state to a double-Q $C_{4}$ state Avici:2014 ; Boehmer:2015 ; Allred_Osborn:2016 , and to a $B_{2g}$ nematic state. All these states appear in the phase diagrams of our Landau theory as well as in the proposed microscopic model with frustrated bilinear-biquadratic interactions. Thus, the proposed mechanism represents a unified description of this rich variety of nematic orders. Because this unified description involves the magnetic degrees of freedom, this overall understanding suggests the important role of spin interactions in promoting the emergent properties of the iron-based materials including their high temperature superconductivity.

Finally, classifying nematic order through broken rotational symmetry goes back to its liquid crystal root. Ours is the first to frame it via broken mirror symmetry, which is natural in crystalline settings. Our approach will likely be important in the context of electronic topology as well, where non-local symmetries such as mirror symmetry play an important role.

In conclusion, we have introduced a framework for nematic orders by broken mirror symmetries. Using this approach, we have advanced a mechanism for a $B_{2g}$ nematic order and for a robust understanding of quasi-degenerate nematic channels. The mechanism provides a unified description of nematicity in iron-based superconductors, and elucidates the physics of the FeSCs in the heavily hole-doped regime .

I

Supplemental Material for Broken mirror symmetry, incommensurate spin correlations, and $B_{2g}$ nematic order in iron pnictides

I.1 Ginzburg-Laudau theory for $B_{1g}$ and $B_{2g}$ nematicities with $(q,0)/(0,q)$ and $(q,q)/(q,-q)$ incommensurate magnetic fluctuations

Since $\mathbb{m}(q,0)$ and $\mathbb{m}(0,q)$ are related by the reflection operator $\sigma_{d/d^{\prime}}$: $\sigma_{d/d^{\prime}}\mathbb{m}(q,0)=\mathbb{m}(0,q)$, and $\mathbb{m}(q,q)$ and $\mathbb{m}(q,-q)$ are related by the reflection operator $\sigma_{x/y}$: $\sigma_{x/y}\mathbb{m}(q,q)=\mathbb{m}(q,-q)$, a general free energy should be invariant under these two reflection operations. Correspondingly, the Ginzburg-Landau free energy takes the following form:

[TABLE]

where $r_{B2}=\frac{r-g}{2}$,$r_{B1}=\frac{r+g}{2}$, and $\mathbb{m}_{1}=\mathbb{m}(q,q)$, $\mathbb{m}_{2}=\mathbb{m}(-q,q)$, $\mathbb{m}_{3}=\mathbb{m}(q,0)$, $\mathbb{m}_{4}=\mathbb{m}(0,q)$. Here we assume $u_{1}>u_{2},u_{1}>u_{3}$, $v_{1}>v_{2},v_{1}>v_{3}$ so that we can neglect higher order terms of the moments in the free energy expansion. The complete phase diagram of the above free energy is very complicated, and here we only show the results for all magnetic moments being in parallel, $\mathbb{m}_{1}\parallel\mathbb{m}_{2}\parallel\mathbb{m}_{3}\parallel\mathbb{m}_{4}$.

Taking the derivatives of the free energy with respect to $|\mathbb{m_{1}}|$, $|\mathbb{m_{2}}|$, $|\mathbb{m_{3}}|$ and $|\mathbb{m_{4}}|$, we obtain the following saddle-point equations:

[TABLE]

These equations lead to the following saddle-point solutions:

(1) a paramagnetic phase where all magnetic moments vanish;

(2) a $(q,q)/(-q,q)$ phase supporting $B_{2g}$ nematic order with $|\mathbb{m_{1/2}}|$=$\sqrt{-\frac{r-g}{u_{1}-u_{2}}}$ while other magnetic moments vanish;

(3) a $(q,0)/(0,q)$ phase supporting $B_{1g}$ nematic order with $|\mathbb{m_{3/4}}|$=$\sqrt{-\frac{r+g}{v_{1}-v_{2}}}$ while other magnetic moments vanish;

(4) a $C_{4}$ double-Q $(q,q)+(-q,q)$ phase with $|\mathbb{m_{1}}|$=$|\mathbb{m_{2}}|$=$\sqrt{-\frac{r-g}{2(u_{1}-u_{3})}}$ and $|\mathbb{m_{3}}|$=$|\mathbb{m_{4}}|$=0;

(5) a $C_{4}$ double-Q $(q,0)+(0,q)$ phase with $|\mathbb{m_{1}}|$=$|\mathbb{m_{2}}|$=0 and $|\mathbb{m_{3}}|$=$|\mathbb{m_{4}}|$=$\sqrt{-\frac{r+g}{2(v_{1}-v_{3})}}$;

(6) a $B_{1g}$ and $B_{2g}$ coexisting phase with $|\mathbb{m_{1}}|$ or $|\mathbb{m_{2}}|$=$\sqrt{-\frac{(v_{1}-v_{2})(r-g)-(w_{1}+w_{2})(r+g)}{(u_{1}-u_{2})(v_{1}-v_{2})-(w_{1}+w_{2})^{2}}}$ and $|\mathbb{m_{3}}|$ or $|\mathbb{m_{4}}|$= $\sqrt{-\frac{(u_{1}-u_{2})(r+g)-(w_{1}+w_{2})(r-g)}{(u_{1}-u_{2})(v_{1}-v_{2})-(w_{1}+w_{2})^{2}}}$;

(7) a mixture phase supporting $B_{2g}$ nematicity with $|\mathbb{m_{1}}|$ or $|\mathbb{m_{2}}|$ coexisting with $|\mathbb{m_{3}}|$ and $|\mathbb{m_{4}}|$ where $|\mathbb{m_{1}}|$ or $|\mathbb{m_{2}}|$= $\sqrt{-\frac{2(v_{1}-v_{3})(r-g)-(w_{1}+w_{2})(r+g)}{2((u_{1}-u_{2})(v_{1}-v_{3})-(w_{1}+w_{2})^{2})}}$ and $|\mathbb{m_{3}}|$ = $|\mathbb{m_{4}}|$=$\sqrt{-\frac{(u_{1}-u_{2})(r+g)-2(w_{1}+w_{2})(r-g)}{2((u_{1}-u_{2})(v_{1}-v_{3})-(w_{1}+w_{2})^{2})}}$;

(8) a mixture phase supporting $B_{1g}$ nematicity with $|\mathbb{m_{3}}|$ or $|\mathbb{m_{4}}|$ coexisting with $|\mathbb{m_{1}}|$ and $|\mathbb{m_{2}}|$ where $|\mathbb{m_{1}}|$=$|\mathbb{m_{2}}|$=$\sqrt{-\frac{(v_{1}-v_{2})(r-g)-2(w_{1}+w_{2})(r+g)}{2((u_{1}-u_{2})(v_{1}-v_{3})-(w_{1}+w_{2})^{2})}}$ and $|\mathbb{m_{3}}|$ or $|\mathbb{m_{4}}|$=$\sqrt{-\frac{2(u_{1}-u_{3})(r+g)-(w_{1}+w_{2})(r-g)}{2((u_{1}-u_{2})(v_{1}-v_{3})-(w_{1}+w_{2})^{2})}}$;

(9) a $C_{4}$ phase with coexisting $\mathbb{m_{1}}$, $\mathbb{m_{2}}$, $\mathbb{m_{3}}$, and $\mathbb{m_{4}}$ moments where $|\mathbb{m_{1}}|$=$|\mathbb{m_{2}}|$=$\sqrt{-\frac{(v_{1}-v_{3})(r-g)-(w_{1}+w_{2})(r+g)}{2((u_{1}-u_{3})(v_{1}-v_{3})-(w_{1}+w_{2})^{2})}}$ and $|\mathbb{m_{2}}|$=$|\mathbb{m_{3}}|$= $\sqrt{-\frac{(u_{1}-u_{3})(r+g)-(w_{1}+w_{2})(r-g)}{2((u_{1}-u_{3})(v_{1}-v_{3})-(w_{1}+w_{2})^{2})}}$.

Among these, solutions (1),(4),(5), and (9) preserve full $D_{4h}$ point group symmetry. Other solutions, however, have symmetries lower than $D_{4h}$. Solutions (2) and (7) support $B_{2g}$ nematicity by preserving mirror symmetry $\sigma_{d/d^{\prime}}$ and breaking $\sigma_{x/y}$ symmetry. Solution (3) and (8) support $B_{1g}$ nematicity by preserving mirror symmetry $\sigma_{x/y}$ and breaking $\sigma_{d/d^{\prime}}$ symmetry. Solution (6) has the lowest symmetry, breaking $\sigma_{d/d^{\prime}}$, $\sigma_{x/y}$, and their product, $C_{4}$, symmetries. It allows coexistence of $B_{1g}$ and $B_{2g}$ nematic orders.

These phases can be classified into two phase diagrams, which are shown in Fig.S1. The white region corresponds to the paramagnetic solution (1), the blue region refers to either solution (2) or (4), the orange region specifies either solution (3) or (5), and the green region corresponds to solution (6) if the blue region refers to (2) and the orange region refers to (3), to solution (7) if the blue region refers to (2) and the orange one refers to (5), to solution (8) if the blue region refers to (4) and the orange refers to (3), and to (9) if the blue region refers to (4) and the orange one refers to (5). If $(u_{1}-u_{2})(v_{1}-v_{2}),(u_{1}-u_{2})(v_{1}-v_{3}),(u_{1}-u_{3})(v_{1}-v_{2})$ and $(u_{1}-u_{3})(v_{1}-v_{3})$ are all smaller than $(w_{1}+w_{2})^{2}$, the Ginzburg-Landau free energy has a bicritical point separating the blue and orange regions, as shown in Fig.1(a). In this case, the $B_{1g}$ and $B_{2g}$ nematic orders can not coexist. They are separated by a first-order transition. By contrast, in Fig.1(b), there is a coexistence region (green) for $B_{1g}$ and $B_{2g}$ nematic orders. All transitions in this case are second-order and there is a tetracritical point in the $r-g$ plane. Note that the two phase diagrams discussed here is similar to those of a two-component $\phi^{4}$ model in Ref. [19]. We summarize these results in Table S1.

Supplemental Table S1. Different conditions for the phase diagrams of FIG.S1

[TABLE]

I.2 Ginzburg-Laudau theory for nematic orders with generic $(q_{1},q_{2})$ incommensurate magnetic fluctuations

In this section, we perform a large-$N$ calculation [20,21] for the Ginzburg-Landau theory involving incommensurate magnetic moments with generic wave vectors $(q_{1},q_{2})$. Because the symmetry group $D_{4h}$ naturally connects states with wavevectors $\mathbb{Q}_{1}=(q_{1},q_{2}),\mathbb{Q}_{2}=(-q_{2},q_{1}),\mathbb{Q}_{3}=(q_{2},q_{1})$ and $\mathbb{Q}_{4}=(-q_{1},q_{2})$, we construct a Ginzburg-Laudau action in terms of these four magnetic states (here we assume all magnetic moments are in parallel for simplicity), $\mathbb{m}_{1}=\mathbb{m}(q_{1},q_{2}),\mathbb{m}_{2}=\mathbb{m}(-q_{2},q_{1}),\mathbb{m}_{3}=\mathbb{m}(q_{2},q_{1})$ and $\mathbb{m}_{4}=\mathbb{m}(-q_{1},q_{2})$:

[TABLE]

where $\mathbb{m}_{i}(\bm{k})\equiv\mathbb{m}_{i}(\mathbb{Q}_{i}+\bm{k})$. We define three nematic order parameters: $\sigma_{B1}=\mathbb{m}_{1}^{2}+\mathbb{m}_{4}^{2}-\mathbb{m}_{2}^{2}-\mathbb{m}_{3}^{2},\sigma_{B2}=\mathbb{m}_{1}^{2}+\mathbb{m}_{3}^{2}-\mathbb{m}_{2}^{2}-\mathbb{m}_{4}^{2}$ and $\sigma_{A2}=\mathbb{m}_{1}^{2}+\mathbb{m}_{2}^{2}-\mathbb{m}_{3}^{2}-\mathbb{m}_{4}^{2}$, which are respectively conserved under $\sigma_{d/d^{\prime}}$, $\sigma_{x/y}$, and $\sigma_{d/d^{\prime}}\times\sigma_{x/y}$. Then we rewrite the quartic term of the action as

[TABLE]

where $v_{1}=(u_{1}+3u_{2})/4$, $v_{2}=(u_{2}-u_{1})/4$. If $u_{2}>u_{1}$, then $v_{2}>0$, allowing for nematic orders.

In the large-$N$ limit, we rescale $v_{1}$ and $v_{2}$ to $v_{1}/N$ and $v_{2}/N$ and perform the Hubbard-Stratonovich transformation,

[TABLE]

We then arrive at

[TABLE]

Next we express $\mathbb{m}_{i}=(\sqrt{N}m_{i},\bm{\pi}_{i})$, where $m_{i}$ refers to the longitudinal ordered component and $\bm{\pi}_{i}$ are the transverse modes with $N-1$ components. We integrate out the transverse modes, treat $\sigma_{i},\lambda$, $\sigma_{B1}$, $\sigma_{B2}$ and $\sigma_{A2}$ at the saddle point level, and obtain the following free energy density (here we have redefined $i\lambda\longrightarrow\lambda$, such that $\lambda$ is real):

[TABLE]

with

[TABLE]

and the following saddle point equations:

[TABLE]

When $\sigma_{B1}=\sigma_{B2}=\sigma_{A2}=0$, we can immediately get $|m_{1}|=|m_{2}|=|m_{3}|=|m_{4}|\equiv m$ from Eqs.($\ref{N1},\ref{N2}$ and S23). If $m=0$, then the solution is a paramagnetic state. On the other hand, if $m\neq 0$, then the solution refers to a $C_{4}$ magnetic state which remains $C_{4}$ rational symmetry. This solution, however, is not physical in our classical model due to Mermin-Wagner theorem; $\partial g/\partial\lambda$ blows up in the state.

On the contrary, when $m_{1}=m_{2}=m_{3}=m_{4}=0$ is a solution of Eqs.(S24-S27), we rearrange Eqs.(S20-S23), and have the following equations:

[TABLE]

Since the above four equations have the same form, we only need to treat the following equation:

[TABLE]

where $x=\sigma_{B1}+\sigma_{B2}+\sigma_{A2},-\sigma_{B1}-\sigma_{B2}+\sigma_{A2},-\sigma_{B1}+\sigma_{B2}-\sigma_{A2}$ and $\sigma_{B1}-\sigma_{B2}-\sigma_{A2}$. This equation has maximally two solutions and one can verify that these solutions can be classified into three types:

1. $x=0$, which is disordered with $\sigma_{B1}=\sigma_{B2}=\sigma_{A2}=0$, and the value of $\lambda=\lambda_{0}(r,v_{1})$ can be determined:

[TABLE]

2. $x=\pm\sigma$, ordering of one out of the three nematic orders with $\sigma=\sigma_{B1}\neq 0/\sigma_{B2}\neq 0/\sigma_{A2}\neq 0$, we have:

[TABLE]

3. $x=-3\sigma$, or $x=\sigma$, corresponds to simultaneous ordering of all three nematic components with $\sigma=|\sigma_{B1}|=|\sigma_{B2}|=|\sigma_{A2}|$ (although $x=3\sigma,x=-\sigma$ is another set of solutions, the corresponding free energy is larger than the former case, so we neglect it.), we have:

[TABLE]

To understand the nature of these nematic phases, we expand the free energy in Eqs.(I.2-S19) and the related saddle point equations in Eqs.(S20-S23) to third order in $\{\sigma_{B1},\sigma_{B2},\sigma_{A2},\delta\lambda\}$, where $\delta\lambda=\lambda-\lambda_{0}$. We next substitute the solution of $\delta\lambda$ with respect to $\sigma_{B1},\sigma_{B2}$ and $\sigma_{A2}$ back into the free energy, and we have an effective free energy for the three nematic orders:

[TABLE]

where $r_{\sigma}=\frac{1}{4v_{2}}-\frac{1}{4\pi c}(\frac{1}{\lambda_{0}}-\frac{1}{\lambda_{0}+c\Lambda^{2}}),a=\frac{1}{2\pi c}(\frac{1}{\lambda_{0}^{2}}-\frac{1}{(\lambda_{0}+c\Lambda^{2})^{2}})$. The trilinear term in Eq.(S38) is the manifestation of the discrete symmetry of $D_{4h}$, since identical representation $A_{1g}=B_{1g}\times B_{2g}\times A_{2g}$. This term accounts for the last solution in the above list, and makes the physics of this model very different from that of a Heisenberg or XY model with a continuous symmetry.

I.3 Evolution of spin excitations with hole doping

As shown in Fig.S2, neutron scattering [22] on K doped BaFe2As2 compound shows that the spin excitations contain rich incommensurate magnetic fluctuations. The incommensurate $(q,q)$ fluctuations at high energies are considerably softened with increasing the K (hole doping) concentration.

I.4 Details on the numerical calculations

We determine the phase diagram of the classical bilinear-biquadratic model in Eq.(10) of the main text by using the Luttinger-Tisza method [23] and verified by Monte Carlo simulations at $T/J_{1}=0.01$ on lattices with size up to $64\times 64$. The model parameters used in the phase diagram in Fig.3(a) of the main text are $K_{1}/J_{1}=0.8$, $J_{3}=0.1$, $J_{3}=K_{2}=K_{3}=0$. To show that the $(2\pi/3,2\pi/3)$ state (labeled by the red line in Fig.3(a)) can indeed be stabilized for non-zero $K_{3}$ values, here we show the phase diagram of the model for $K_{3}/J_{1}=-0.2$ while keeping all the other parameters same as those in Fig.3(a) of the main text. From Fig. S3 one clearly sees that the phase diagram is similar to that in Fig.3(a) of the main text, and the $(2\pi/3,2\pi/3)$ state (labeled by the red dashed line) is stabilized for $K_{3}\neq 0$.

In the DMRG calculation, we choose two types of lattice geometries: both the rectangular (RC) and tilted (TC) cylinders, which are denoted as RC/TC$L_{y}\times L_{x}$, where $L_{x}(y)$ is the number of sites along the $x$ ($y$) direction, respectively. We performed DMRG simulations with 2000 $SU(2)$ DMRG states, and the truncation error is around $10^{-5}$ to ensure the accuracy of the results.