Quantum transitions of nematic phases in a spin-1 bilinear-biquadratic model and their implications for FeSe

Wen-Jun Hu,1,2 Hsin-Hua Lai,2 Shou-Shu Gong ,3 Rong Yu,4 Elbio Dagotto,1,5 and Qimiao Si2,* 1Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA 2Department of Physics and Astronomy & Rice Center for Quantum Materials, Rice University, Houston, Texas 77005, USA 3Department of Physics, Beihang University, Beijing 100191, China 4Department of Physics, Renmin University of China, Beijing, 100872, China 5Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA


I. INTRODUCTION
Understanding the iron-based superconductors (FeSCs) has been a subject of extensive research in recent years [1][2][3][4][5]. The initial interest started with the discovery of superconductivity in the iron pnictides. More recently, iron chalcogenides have provided considerable material variety to this intriguing field and reached the new record of superconducting transition temperature (T c ) in FeSCs. These include the potassium iron selenides and other intercalated FeSe systems [6], as well as the single-layer FeSe built on substrates [7,8]. Because all these record-breaking materials involve FeSe as a building block, it is important to understand the physics of bulk FeSe [9,10]. Indeed, there is a vast current interest in this system, which possesses the simplest structure among the FeSCs. In contrast to the standard case of the iron pnictides, where a tetragonal-to-orthorhombic structural phase transition microscopic ingredients for the normal state of the FeSCs as well as to elucidate the degree to which the mechanism for superconductivity is universal across the many varieties of FeSCs.
In this paper we address this issue by exploring the quantum phases and phase transitions related to the nematic phase of FeSe. We focus on studying a spin-1 bilinear-biquadratic model on the square lattice, which has been considered before to understand the exotic magnetism and nematic order of the iron-chalcogenide superconductors [24,25,[31][32][33][34][35][36] but has not been systematically studied to understand the quantum phases and phase transitions of FeSe under pressure. We implement both the site-factorized wave-function analysis and large-scale DMRG method on this model. In general, models with different active microscopic degrees of freedom will have different types of phases in their phase diagrams and, thus, different kinds of quantum phase transitions. In our case, we find four stable spin dipolar and quadrupolar phases, including the Néel antiferromagnetic order, the (π, 0) collinear antiferromagnetic phase (CAFM), the (π/2, π ) antiferromagnetic phase (AFM*), and the (π, 0) antiferroquadrupolar phase (AFQ), and obtain the ground-state phase diagrams. Furthermore, we apply our theoretical results for the understanding of the low-temperature phase diagram that has been indicated by recent experiments of the NMR [37] and x-ray scattering [38] measurements in pressurized FeSe. These experiments have demonstrated that lowering temperature induces a tetragonal to orthorhombic (OR) transition, which accompanies a magnetic transition.
It is important to clarify that the actual values of the bilinear and biquadratic couplings, used in our study as free parameters to construct the phase diagrams, could be fixed by analyzing a higher level, and far more difficult, multiorbital Hubbard models (see Ref. [39] and references therein) and mapping the low-energy states into the spin-1 model used here. While this study will be carried out in the future, we note that in the bad-metal regime of such a multiorbital setting, the biquadratic interaction is expected to be sizable [40]. Thus, in the present effort generic phase diagrams of the spin-1 bilinear-biquadratic model varying couplings in a broad range will be presented. While future work can clarify with precision where each particular FeSC material is located in our phase diagram, here we focus on whether the overall phase diagram hosts nematic phases and their transitions that pertain to the properties of bulk FeSe.
Our paper is organized as follow. In Sec. II we introduce the spin-1 bilinear-biquadratic model and describe the computational details. Section III contains the results of the site-factorized wave-function analysis. In Secs. IV and V we present our DMRG results of the spin and quadrupolar structure factors, as well as the nematic order parameters. Finally, we provide our discussions and conclusions in Sec. VI.

II. MODEL AND METHODS
Our starting point is a spin-1 bilinear-biquadratic model on the square lattice [24,25,[31][32][33]35], which is defined as Here S i is a spin-1 operator at site i, which also forms the quadrupolar operator Q i , with five independent components: The biquadratic term in Eq. (1) can be re-expressed as In addition, J i j and K i j are, respectively, the bilinear and biquadratic couplings between the spins at sites i and j, with the pair i j denoting distinct bonds on the square lattice. The consideration of biquadratic interaction is quite necessary for spin-1 systems. Note that the ab initio method based on density functional theory (DFT) has been used to extract the biquadratic couplings [26]. However, this is a challenging task given that (i) FeSe is strongly correlated and (ii) quadrupoles, being rank-2 objects, do not efficiently couple to the singleparticle degrees of freedom [35] that come into the DFT approach. (For related reasons, the DFT simulation may only study different magnetic orders, which cannot explain the nonmagnetic phase in the FeSe.) The fact that J 2 is comparable with J 1 but K 2 does not appear in the DFT results illustrates this difficulty [26]. For a spin system with furtherneighbor dipolar interactions, it is reasonable to consider the quadrupolar interactions as well.
Following the idea of proposing the (π, 0) antiferroquadrupolar order phase as the candidate of the nonmagnetic phase in the FeSe [35], here we consider interactions beyond nearest neighbor up to the third neighbors. For a minimal model without loss of generality, we consider the nearestneighbor bilinear interaction J 1 = 1 as the energy unit, with varying second-neighbor interaction J 2 . In addition, we consider the first three neighbors of the biquadratic interactions to have the same strength for the purpose of simplifying the model and reducing the number of parameters, −K 1 = K 2 = −K 3 = K > 0. We will demonstrate the robustness of our results by studying the cases with variations of these parameters. For convenience, the present model will be referred to as the J 1 -J 2 -K model. In our site-factorized wave-function analysis and DMRG calculations, we can consider both magnetic and nonmagnetic phases in our model Eq. (1).
For a spin-1 model possibly harboring purely magnetic order, purely quadrupolar order, or coexisting magnetic and quadrupolar orders, it is convenient to choose the timereversal invariant basis of the SU(3) fundamental representation, namely where we abbreviate |S z = ±1 ≡ | ± 1 (|S z = 0 ≡ |0 ) and |1 ≡ | − 1 . Within this basis, the site-factorized wave functions at each site i which characterize any possible ordered state with short-ranged correlations can be expressed as where d x i , d y i , d z i are complex numbers and can be re-expressed in the vector form with the basis {|x , |y , |z } as It is convenient to separate the real and imaginary parts of d i as d i = u i + iv i . The normalization of the wave function leads to the constraint d i ·d i = 1, or equivalently, u 2 i + v 2 i = 1, and the overall phase can be fixed by requiring In a pure quadrupolar state, d will take either a real or imaginary value, but not both, and the associated director is parallel to the director vector d. This is to be contrasted with a magnetic order, for which d contains both real and imaginary components, thus yielding a dipolar magnetic moment. Within this framework, we can determine the spin operator from S i = 2u i × v i . In terms of the components of d, the spin and quadrupolar operators can be written as In addition to the site-factorized wave-function analysis, we also study the model Eq. (1) by the density matrix renormalization group (DMRG) with spin rotational SU(2) symmetry. We perform the DMRG simulations on L × 2L cylindrical systems with L = 4, 6, 8 in the y direction. The cylinder geometry used here has open boundary conditions along the x direction and periodic boundary conditions along the y direction. We keep up to 4000 SU(2) DMRG states. In the Néel AFM and the (π, 0) collinear antiferromagnetic phase (CAFM), the truncation error is around 10 −6 , while in the (π/2, π ) AFM * and the (π, 0) antiferroquadrupolar phase (AFQ) the truncation error is around 10 −5 . These small truncation errors ensure us to obtain accurate DMRG results.

III. SITE-FACTORIZED WAVE-FUNCTION ANALYSIS
To explore the possible quantum phases of the model Eq. (1), we start from an analysis based on a site-factorized wave-function analysis [24,33,35,41]. In this framework we can re-express the model Hamiltonian as In order to obtain the variational phase diagram, we should numerically minimize the Hamiltonian above. In the present analysis we consider four stable spin dipolar or quadrupolar ordered phases as illustrated in Fig. 1, including the CAFM, Néel AFM, and (π, 0) AFQ, as well as a newly discovered magnetic phase dubbed (π/2, π ) AFM * . Note that here we do not consider the ferroquadrupolar order discussed in Ref. [33] as a candidate for the nonmagnetic phase but only consider the (π, 0) AFQ order. From the perspective of purely theoretical explorations, the two types of quadrupolar states are both intriguing. However, the ferroquadrupolar order itself does not generate a nematic order.
The Néel AFM in Fig. 1(b) shows the conventional spin pattern with the nearest-neighbor spins antiparallel to each other. In the CAFM phase, shown in Fig. 1(c), the nearestneighbor spins are parallel to each other along one direction, while they are antiparallel to each other along the other direction. The (π, 0) AFQ with staggered quadrupolar order along one direction, as shown in Fig. 1(d), can be characterized by having mutually orthogonal nearest-neighbor directors, i.e., The novel (π/2, π ) AFM * , shown in Fig. 1(e), is a new phase where the spin direction along the x axis rotates with a commensurate period of four sites while its period along the y axis is still two sites reflecting the wave vector (π/2, π ). This AFM * state is nematic since it spontaneously breaks the lattice C 4 symmetry, by choosing between two degenerate wave vectors q = (π/2, π ) and (π, π/2). As an illustration, the spin configuration of the (π/2, π ) AFM * state along the x axis may take the 4-site periodic pattern as {|S z = 1 , |S x = 1 , |S z = −1 , |S x = −1 }, while the spin orientation still takes the conventional staggered pattern along the y axis. Importantly, in the (π/2, π ) AFM * phase, the q = (π, 0) quadrupolar order parameter is nonzero as well, which is responsible for the stability of this phase.
Within the site-factorized wave-function studies, the energy per site of each phase can be obtained (see Appendix B) as follows: For all energies we have neglected the constant term K i j in Eq. (7). Note that the SU(3) analysis was considered earlier, in particular in Ref. [35]. It is also worth emphasizing that 023359-3 these results capture the quantum fluctuations inherent to the S = 1 case, as discussed in some detail in Appendix B. Using these energies, we can analytically determine the boundaries as follows: (1) Phase boundary between Néel AFM and (π, 0) AFQ: (2) Phase boundary between Néel AFM and (π/2, π ) AFM * : K − 8J 2 + 4 = 0.
(5) Phase boundary between CAFM and (π, 0) AFQ: Using these equations, by employing the site-factorized wave-function analysis, we obtain the variational phase diagram shown in Fig. 1(a).
We would like to mention that within the site-factorized wave-function analysis, a phase with coexistent magnetic and quadrupolar orders at different real-space sites, dubbed AFMQ, has been found by minimizing the energy. To our best understanding, the AFMQ shows the same period as that of (π/2, π ) AFM * , but is an inhomogeneous phase with finite magnetic and quadrupolar orders at different realspace columns (rows). The spin pattern in AFMQ presents a staggered pattern between magnetically ordered sites, while the quadrupolar pattern is ferroquadrupolar (FQ), i.e., the quadrupolar directors are all parallel to each other. The sitefactorized wave functions between the magnetically ordered sites and the quadrupolar sites are orthogonal to each other. Since this regime with coexisting magnetic and quadrupolar orders at different columns (rows) appears between the (π/2, π ) AFM * and (π, 0) AFQ, whose period is consistent with both (π/2, π ) AFM * and (π, 0) AFQ, it is likely that this phase is just a transition regime between the purely magnetic phase and the purely quadrupolar phase. Its existence reflects the first-order nature of the transition between the two phases, and it is expected to be destabilized by quantum fluctuations; this is confirmed by our DMRG calculations later. We have therefore ignored this regime in the phase diagram displayed in Fig. 1(a).

IV. DMRG PHASE DIAGRAMS
Our analysis so far has been semiclassical. In order to explore the role of full quantum fluctuations and analyze the model of Eq. (1) in an unbiased way, we have also carried out large-scale density matrix renormalization group (DMRG) calculations [42]. First of all, we selected four points of the four stable phases shown in the phase diagram Fig. 1(a): J 2 = 0.4 and K = 0.6 for the Néel AFM, J 2 = 1.5 and K = 0.5 for the CAFM, J 2 = 1.5 and K = 0.7 for the (π, 0) AFQ, and J 2 = 0.8 and K = 0.35 for the (π/2, π) AFM * . Next, we computed the spin-spin ( S i · S j ) and quadrupolarquadrupolar ( Q i · Q j ) correlation functions for these four points by using DMRG on a L = 8 cylindrical geometry, which contain the real-space spin and quadrupolar configurations displayed in Fig. 2. The spin pattern of the Néel AFM state is in Fig. 2(a); the CAFM state in Fig. 2(b) automatically chooses the antiparallel configuration along the y direction and the parallel configuration along the x direction, due to the cylindrical geometry; the AFQ phase in Fig. 2(c) has the antiferroquadrupolar configuration along the x direction and the ferroquadrupolar configuration along the y direction. In the (π/2, π ) AFM * phase, along the y direction the spin configuration is antiparallel, while along the x direction the spin pattern in Fig. 2(d) indicates that the spins are orthogonal between two nearest-neighbor sites. This selection of otherwise degenerate nematic states is induced by the small symmetry-breaking geometry of the cylinders used in DMRG.
Next we consider the evolutions of the zero-temperature phases as a function of the biquadratic K coupling to obtain the ground-state phase diagram for fixed J 2 . In order to identify the phases we encounter, we have calculated the static spin and quadrupolar structure factors defined as where i, j are restricted to be only partially summed over the L × L sites in the middle of the cylinder so that the finite-size effects are reduced [43]. Figure 3 displays the DMRG results for J 2 = 1.5. Following the scaling behavior in magnetic order states which can be obtained from spin-wave theory, we show the order parameters as a function of 1/L. Note that, given the numerically intensive nature of the DMRG calculations for an extended parameter space, we only crudely extrapolate the data for the finite-size scaling, which is sufficient to show whether the order vanishes or not. The evolutions of the spin (m 2 S ) and quadrupolar (m 2 Q ) structure factors is shown in Fig. 3(a). Here, in the range K < 0.65, the peak at momentum q = (0, π ) of the spin structure factor suggests the presence of the CAFM phase, while the quadrupolar structure factor has its peak at ) and quadrupolar (m 2 Q ) structure factors obtained from the DMRG calculations on the 8 × 16 cylinders for J 2 = 1.5. Both structure factors display dramatic changes at K 0.65, indicating a phase transition from the (0, π ) CAFM to the (π, 0) AFQ. In the (π, 0) AFQ, m 2 Q exhibits a characteristic peak at (π, 0). Finite-size scaling for the spin (b) and quadrupolar (c) structure factors at different values of K for J 2 = 1.5. For the spin structure factor (b), the highest peak of m 2 S in its momentum distribution is shown. For the quadrupolar structure factors (c), the intensity at q = (π, 0) is plotted. According to the scaling, the value K = 0.65 is close to the phase boundary. Lines are guides to the eye. q = (0, 0). On the other hand, for K > 0.65, the quadrupolar structure factor develops a clear peak at q = (π, 0), while the peak in the spin structure factor has melted, indicating the (π, 0) AFQ order. Due to the cylindrical geometry, the CAFM state automatically selects the configuration with q = (0, π ), whereas the AFQ phase selects q = (π, 0). Furthermore, we examine the finite-size scalings of the spin (m 2 S ) and quadrupolar (m 2 Q ) order parameters at different values of K in Figs. 3(b) and 3(c). Upon increasing K, the spin order parameter m 2 S (0, π ) decreases and vanishes when K > 0.65. Instead, the quadrupolar order parameter m 2 Q (π, 0) develops for K > 0.65. The behavior of both these order parameters indicate a direct phase transition from CAFM to (π, 0) AFQ, and K = 0.65 is roughly the location of the phase boundary based on the finite-size scalings. Beyond the discovery of the (π, 0) AFQ phase as a genuine ground state in [35], here we show that the CAFM and (π, 0) AFQ appear as nearby phases through the variation of the biquadratic couplings.
We also performed DMRG calculations at fixed J 2 = 0.8, with results in Fig. 4. At small and large values of K, the system is in the CAFM and (π, 0) AFQ phases, respectively. In the intermediate region there is a new magnetic phase emerging, with a new peak in the spin structure factor m 2 S developing at q = (π/2, π ) as shown in Fig. 4(a), for example at K = 0.35. Also, we find the peak at q = (π, 0) for the quadrupolar structure factor m 2 Q . These results are consistent with our site-factorized wave-function analysis for the (π/2, π ) AFM * order in Sec. III, which suggests the coexistence of the (π/2, π ) magnetic and (π, 0) quadrupolar orders. Figures 4(b) and 4(c) contain the finite-size scalings -π π 0 -π π 0 -π π 0 -π π 0 -π π 0 -π S displays a transition from (0, π ) CAFM, through (π/2, π ) AFM * , and finally to (π, 0) AFQ. In both (π/2, π ) AFM * and (π, 0) AFQ, m 2 Q exhibits a characteristic peak at (π, 0). Finite-size scalings for the spin (b) and quadrupolar (c) structure factors at different values of K and at fixed J 2 = 0.8. For the spin structure factor (b), the highest peak of m 2 S in its momentum distribution is shown. For the quadrupolar structure factors (c), the intensity at q = (π, 0) is plotted. The momentum (π/2, π ) is not an allowed lattice vector on the 6 × 6 cluster, which is responsible for the apparent nonmonotonic dependence of m 2 S vs 1/L in the range 0.25 < K < 0.65. Lines are guides to the eye. of spin m 2 S and quadrupolar m 2 Q order parameters at different values of K. Although the momentum (π/2, π ) is not an allowed lattice vector on the 6 × 6 cluster, and creates the apparently nonmonotonic dependence of m 2 S vs 1/L in the range 0.25 < K < 0.65 shown in Fig. 4(b), still the results clearly demonstrates the nonzero values of both order parameters m 2 S and m 2 Q at K ∼ 0.6 after the finite-size scaling. This signature is less clear for smaller K, but we believe that both orders already coexist there. As inferred from the finite-size scalings, we find two quantum phase transitions, with the first one happening around K 0.25 and the second around K 0.65.
We would like to mention that the particular parameter cut in the model is for presentation purposes only. The nematic phases actually span over large parameter regions, and the general features of the nematic phase diagram remain the same if we choose other fixed parameters (see Appendix A).

V. NEMATICITY
To characterize the nematicity in the different phases, we introduce two nematic order parameters σ S B1g and σ Q B1g defined as wherex andŷ denote the unit length vectors along the x and y directions, respectively, and N m is the number of sites of the two columns in the middle of the cylinder. Analyses of these nematic order parameters have been shown efficient to determine the lattice rotational symmetry breaking in the DMRG calculations on the cylinder geometry [44]. The absolute value of the nematic order parameters as a function of K at J 2 = 1.5 and 0.8 are, respectively, presented in Figs. 5(a) and 5(b). Comparing σ S B1g and σ Q B1g clarifies whether the antiferromagnetic or antiquadrupolar fluctuations dominate the contributions to the nematic order. We find σ S B1g dominating over σ Q B1g inside the CAFM phase [K 0.65 (0.25) for J 2 = 1.5 (0.8)], and vice versa inside the (π, 0) AFQ phase (K 0.65 for both values of J 2 ). We also notice that the crossings of the two nematic order parameters occurs exactly at the location where the quadrupolar order at q = (π, 0) develops. This reflects the different types of fluctuations that are responsible for the nematic order on the two sides of the quantum phase transition. For J 2 = 1.5, the crossing is at the boundary between CAFM and (π, 0) AFQ, while for J 2 = 0.8 this crossing occurs at the boundary between CAFM and (π/2, π ) AFM * .
Finally, we show the finite-size scaling for the nematic order parameters σ S B1g and σ Q B1g of the four phases in Fig. 6. In the Néel AFM state, although there are small values of nematicity, due to the cylindric geometry used in the DMRG calculations, both σ S B1g and σ Q B1g decay fast and vanish with increasing size. For the other three phases, the finite-size scaling clearly indicates the presence of nonzero nematic orders in the thermodynamic limit.

VI. DISCUSSIONS AND CONCLUSIONS
We now discuss the implications of our results for the iron chalcogenides. Our work leads to a possible understanding of the properties of FeSe based on the presented phase diagram in Fig. 1(a). For clarity, we show a schematic phase diagram of the nematic phases in the inset of Fig. 7. Note that applying pressure increases the kinetic energy without affecting the local interactions as much and, thus, amounts to increasing w, the coherent electron spectral weight. Qualitatively comparing the ambient-pressure FeSe with the pressurized FeSe is similar to comparing the ambient-pressure FeSe with the typical iron arsenides. Because the ambient-pressure FeSe is more strongly correlated than the latter, it is expected to be more frustrated and, correspondingly, having a larger K/J ratio. In other words, under pressure, K/J should decrease, which is the parameter trajectory we have proposed in Fig. 7 for FeSe as a function of pressure. More microscopically, a nonperturbative procedure for calculating the effective exchange interactions in the bad metal regime has been developed using the slave-boson-type approach [45]. Here the bilinear spinexchange interaction J is given by a two-boson process and turns out to be (1 − w)J c , where w denotes the percentage of the physical electron spectral weight that resides in the coherent part near the Fermi energy and J c is the exchange interaction at the delocalization-localization transition (i.e., when w → 0 + ). A similar procedure for the biquadratic interaction implies that it is given by a four-boson process, and will be on the order of (1 − w) 2 K c . Thus, K/J is expected to be proportional to (1 − w). Because applying pressure enhances  Fig. 1(a)]. There are two possible sequences of quantum phase transitions from the (π, 0) AFQ phase presumably stable at ambient pressure towards the high pressure CAFM (π, 0) phase, as illustrated by the arrows in both the main panel and the inset. 023359-6 the coherent electron spectral w, it is expected to lead to a decrease in K/J.
Thus, pressurizing FeSe may amount to taking a horizontal cut in this phase diagram: we propose two such cuts as candidates for the parameter tuning, which are also illustrated in the inset of Fig. 7. The resulting phase diagram is illustrated in the main panel of Fig. 7, with the system undergoing either a direct transition between the (π, 0) AFQ and CAFM states, or a transition between them through an intermediate (π/2, π ) AFM * regime with coexistence of magnetic and quadrupolar orders. This phase diagram is qualitatively similar to that inferred from recent experiments. While the presence of AFM order at pressures of the order of 2 GPa had been indicated before [46], recent NMR measurements [37] have provided strong evidence that the order achieved by increasing pressure breaks the C 4 symmetry and has a (π, 0) wave vector. The x-ray scattering experiments [38] have also provided evidence that a C 4 symmetry breaking accompanies the magnetic ordering. There are indications in the existing experiments for two stages of phase transitions under pressure [46][47][48][49], with the onset of AFM order around p 1 ≈ 0.8 GPa and a change of the magnetic structure around p 2 ≈ 1.2 GPa [46]. Additional NMR and neutron scattering measurements in the intermediate pressure range 0.8 P 1.7 GPa are especially needed to clarify this issue and ascertain which of the two proposed sequences applies. We reiterate that the set of model parameters we choose is for the purpose of illustrating the phase transitions in Fig. 7. The nematic phases and the transitions actually span over a large parameter regime in the overall phase diagram, which underscores the fact that the proposed physical picture is robust instead of fine tuned.
Regardless of which of the two phase transition sequences is realized, our results have important implications for the single-electron excitations. The (π, 0) AFQ state contains two order parameters. The rank-2 AFQ order parameter does not efficiently couple with the (coherent) conduction electrons near the Fermi surface and, therefore, will not cause a reconstruction of the Fermi surface. However, the nematic order parameter, σ S B1g and σ Q B1g defined in Eqs. (10) and (11), will linearly mix with the occupancy difference in the 3d xz and 3d yz orbitals, thereby generating a splitting of the electronic bands and a distortion of the Fermi surface. All these features are consistent with the observations by the ARPES experiments [16][17][18]23]. Likewise, the (0, π ) CAFM state contains two order parameters. While the nematic order parameter acts similarly as in the AFQ case, distorting the Fermi surface, the AFM order is very different from the AFQ order: it provides a spatially modulated and spin-dependent potential for the conduction electrons, thereby reconstructing the Fermi surface. Thus, our proposed quantum phase transitions will be accompanied by drastic changes in the geometry of the Fermi surfaces. This is consistent with the dramatic evolution of the Fermi surface recently reported in the Shubnikov-de Haas (SdH) oscillation measurements on FeSe [50].
To summarize, the present work has advanced two key results. First, by considering the interplay between the frustrated magnetic interactions, we establish a quantum phase diagram in which the (nonmagnetic) (π ,0) antiferroquadrupolar order is robustly located near the (π, 0) collinear antiferromagnetic order and an intervening (π/2, π ) antiferromagnetic phase, all of which break the C 4 symmetry and thus promote a nematic order. Second, this theoretical result provides the basis to understand qualitatively the quantum phase transition of FeSe under pressure. The similarity of the quantum phase transitions we have identified in the frustrated bilinearbiquadratic model with the experimental observations provides evidence that a similar type of spin physics is important for the emergence of superconductivity in both iron chalcogenides and iron pnictides. This is not to say that the orbital degrees of freedom are decoupled. As discussed above, the nematic order of the spin quadrupolar or dipolar orders will be coupled to the orbital order. Nonetheless, the interactions among the spin degrees of freedom, as described in Eq. (1), will give rise to superconducting pairing in FeSe-and by extension in other iron chalcogenides-in a similar way as they do in the iron pnictides. Thus, our results not only contribute to the understanding of recent experiments in FeSe, but also provide evidence for a common origin of superconductivity across the extensive material classes of iron-based superconductors. More generally, our findings suggest the importance of correlation-induced short-range spin exchange interactions for both the normal state and superconductivity in the iron chalcogenides. This provides a new linkage between the superconductivity of the highest T c iron-based families with that arising in a broad array of strongly correlated electron systems, including the cuprates and heavy fermion metals. shown in Fig. 8 and for (K 1 = −K, K 2 = 0.7K, K 3 = −0.5K ) shown in Fig. 9. In both cases we have chosen J 1 = 1 and J 2 = 1.5 for the DMRG calculations. From the finitesize scalings of the spin (m 2 S ) and quadrupolar (m 2 Q ) order parameters, we find direct phase transitions between the CAFM and (π, 0) AFQ phases for both cases, and the transition points are around K = 0.8 for the parameter set (K 1 = −K, K 2 = 0.9K, K 3 = −0.7K ) and K = 1.1 for (K 1 = −K, K 2 = 0.7K, K 3 = −0.5K ). These results suggest that the nematic phases actually span over large parameter regions in the phase diagram, and the general features of the nematic phase diagram and relation to experiments remain qualitatively the same as described in our conclusions of the main text.

APPENDIX B: GROUND-STATE ENERGY IN THE SITE-FACTORIZED WAVE-FUNCTION APPROXIMATION
Within the site-factorized wave-function approximation of the SU(3) representation, the ground-state energy per site of a certain ordered phase in the S = 1 model can be readily determined. Here we show how this works for the (π, 0) CAFM state, with particular emphasis on the quantum contributions. The generalization to other states is straightforward.
Denote the two sublattices of the CAFM state to be A and B, respectively. Without losing generality, we assume the local wave functions on these two sublattices to be The corresponding directors are, respectively, This gives for the antiferromagnetically coupled bond, and for the ferromagnetic coupled bond. For the CAFM state, a site connects to one AFM and one FM nearest-neighbor bonds, two AFM next-nearest-neighbor bonds, and two FM third-nearest-neighbor bonds on average. Following Eq. (7), one gets the energy per site for the CAFM phase to be E CAFM = −2J 2 + 2J 3 + K 1 + 2K 2 = −2J 2 + K, (B5) when neglecting the constant term K i j in Eq. (7). The contribution from the constant term is 2(K 1 + K 2 + K 3 ) = −2K, and this shifts the energy to E CAFM = −2J 2 + 2J 3 + 3K 1 + 4K 2 + 2K 3 = −2J 2 − K.
Note that this energy is higher than the energy per site of the CAFM state in the classical limit, E c CAFM = −2J 2 + 2J 3 + 2K 1 + 2K 2 + 2K 3 = −2J 2 − 2K. The reason is as follows. In the classical limit, (S A · S B ) 2 = (S z A S z B ) 2 = 1 for an AFM bond. But in the SU(3) representation of the S = 1 model, one can show that for the AFM bond (S A · S B ) 2 = (S z A S z B ) 2 + S + A S − B S − A S + B /4 = 2. The larger value comes from the transverse correlation S + A S − B S − A S + B /4 and reflects the inherent quantum mechanical nature of the AFM state. Note that there are one nearest neighbor and two next-nearestneighbor AFM bonds in the CAFM state, therefore the energy difference between the S = 1 case and the classical limit is E CAFM − E c CAFM = K 1 + 2K 2 = K. By contrast, for any FM bond, the transverse corrections vanish.