Spin nematics in frustrated spin-dimer systems with bilayer structure

TL;DR

This paper investigates spin nematic phases in frustrated bilayer spin-1/2 dimer systems, revealing conditions under which these phases emerge and their relation to other magnetic orders through theoretical and numerical methods.

Contribution

It demonstrates that frustrated two-spin exchange interactions alone can induce spin nematic phases, and explores phase transitions and emergent phases from an SU(4) symmetric point.

Findings

Spin nematic phase exists over a wide parameter range.

Four-spin interactions are not necessary for spin nematic phases.

Continuous phase transitions occur between spin nematic and antiferromagnetic phases.

Abstract

We study frustrated spin-1/2 dimer systems in two dimensions with a bilayer structure, where spins are ferromagnetically coupled in dimers. Our model includes frustrated two-spin exchange interactions as well as four-spin interaction. We pay particular attention to the spin nematic phase, which does not exhibit any magnetic (spin-dipole) order but has a spin-quadrupolar long-range order. Employing a perturbation calculation, a mean-field approximation, and a numerical many-variable variational Monte Carlo method, we determine ground-state phase diagrams on various two-dimensional lattices. It is found that the model exhibits the spin nematic phase with ferro-quadrupolar order in a wide parameter region, in addition to conventional magnetically-ordered phases. In particular, it is shown that even when the four-spin interactions are absent, frustrated two-spin exchange interactions can realize the spin nematic phase as a result of strong interdimer correlations. It is also found that the phase transitions between the spin nematic phase and antiferromagnetic phases can be continuous. Furthermore, we present some exact arguments that various phases including the spin nematic phase and the vector chiral (p-type nematic) phase emerge from an SU(4) symmetric point in the model by the addition of appropriate perturbative interactions. The spin nematic phase generated from the SU(4) point is connected with the spin nematic phase found numerically in the system with only two-spin interactions.

Figures

Tables

Table 1. Table 1: Parameter ranges or regions of various phases in the ground state of the S = 1 𝑆 1 S=1 bilinear-biquadratic models ( ℋ bb subscript ℋ bb \mathcal{H}_{\rm bb} and ℋ ( 1 ) superscript ℋ 1 \mathcal{H}^{(1)} ) on (a) square and honeycomb lattices, and on (b) triangular and kagome lattices. FM and FQ respectively denote the ferromagnetic and ferro-quadrupolar phases, and AFM and AFQ respectively the antiferromagnetic and antiferro-quadrupolar phases. PVBS and TVBC respectively denote plaquette valence-bond solid and trimerized valence-bond crystal. Both AFM3 and AFQ3 represent phases with three-sublattice structure, and AFM120 ∘ does the AFM phase with 120 ∘ superscript 120 120^{\circ} structure. The results for ℋ bb subscript ℋ bb \mathcal{H}_{\rm bb} are taken from Refs.  Tóth  et al. , 2012 ; Niesen and Corboz, 2017 ; Zhao  et al. , 2012 ; Läuchli  et al. , 2006 ; Liu  et al. , 2015 .

(a)
	Square		Honeycomb
Phase	${\mathscr{H}}_{\mathrm{bb}}$ (Refs. Tóth et al., 2012; Niesen and Corboz, 2017)	${\mathscr{H}}^{\left(1\right)}$	${\mathscr{H}}_{\mathrm{bb}}$ (Ref. Zhao et al., 2012)	${\mathscr{H}}^{\left(1\right)}$
FM
FQ
Néel		$J_{\parallel }+J_{\times }>\max (0.49J_4,-\frac{J_4}{4})$		$J_{\parallel }+J_{\times }>\max (0.49J_4,-\frac{J_4}{4})$
Haldane		$0.49J_4>J_{\parallel }+J_{\times }>0.37J_4$	-	-
AFM3		$0.37J_4>J_{\parallel }+J_{\times }>\frac{J_4}{4}$	-	-
AFQ3			-	-
PVBS	-	-
(b)
	Triangular		Kagome
Phase	${\mathscr{H}}_{\mathrm{bb}}$ (Ref. Läuchli et al., 2006)	${\mathscr{H}}^{\left(1\right)}$	${\mathscr{H}}_{\mathrm{bb}}$ (Ref. Liu et al., 2015)	${\mathscr{H}}^{\left(1\right)}$
FM
FQ
AFM120∘		$J_{\parallel }+J_{\times }>\max (\frac{J_4}{4},-1.6J_4)$	-	-
AFQ3			-	-
TVBC	-	-		$J_{\parallel }+J_{\times }>\max (-0.03J_4,-4.2J_4)$
AFQ	-	-

Equations

$$\mathcal{H}$$

$${\mathcal{H}}_{\mathrm{d}}$$

$${\mathcal{H}}_{\parallel }$$

$${\mathcal{H}}_{\times }$$

$${\mathcal{H}}_4$$

$$T_j^{\alpha }$$

$$N_j^{\alpha }$$

$$Q_j^{(1)}$$

$$Q_j^{(2)}$$

$$Q_j^{(3)}$$

$$Q_j^{(4)}$$

$$Q_j^{(5)}$$

$$J_{\parallel }=J_4/4,J_{\times }=0$$

$$J_{\parallel }-\frac{J_4}{4}+J_{\times }=0$$

$$J_{\parallel }=J_4/4,J_{\times }=J_{\mathrm{d}}=0,$$

$${\mathcal{H}}^{(1)}$$

$${\mathcal{H}}_{\mathrm{bb}}=J_{\mathrm{bb}}\sum _{⟨j,j^{\prime }⟩}\left[\cos \theta {\tilde{\mathbf{S}}}_j\cdot {\tilde{\mathbf{S}}}_{j^{\prime }}+\sin \theta {\left({\tilde{\mathbf{S}}}_j\cdot {\tilde{\mathbf{S}}}_{j^{\prime }}\right)}^2\right]$$

$${\mathcal{H}}^{(2)}=J^{(2)}\sum _{⟨j,j^{\prime }⟩}\left[{\left({\tilde{S}}_j\cdot {\tilde{S}}_{j^{\prime }}\right)}^2-1\right]$$

$$J^{(2)}=-\frac{(J_{\parallel }-J_{\times }{)}^2}{8|J_{\mathrm{d}}|}=-\frac{J_{\parallel }^2}{2|J_{\mathrm{d}}|}=-\frac{J_{\times }^2}{2|J_{\mathrm{d}}|}.$$

$$|{\Phi }_{\mathrm{DP}}⟩=\prod _j|\phi {⟩}_j,$$

$$|\phi {⟩}_j=|\phi {⟩}_{\Lambda }=\sum _{{\sigma }_1,{\sigma }_2}{a}_{\Lambda ,{\sigma }_1{\sigma }_2}|{\sigma }_1{\sigma }_2⟩$$

$${\mathcal{H}}_{j{j}^{\prime }}$$

$${\mathbf{T}}_{\Lambda }^{\mathrm{MF}}$$

$${\mathbf{N}}_{\Lambda }^{\mathrm{MF}}$$

$${\mathbf{Q}}_{\Lambda }^{\mathrm{MF}}$$

$$J_{\parallel }=-1,$$

$${\mathbf{T}}_{\mathrm{A}}^{\mathrm{MF}}={\mathbf{T}}_{\mathrm{B}}^{\mathrm{MF}},|{\mathbf{T}}_{\Lambda }^{\mathrm{MF}}|=1,|{\mathbf{N}}_{\Lambda }^{\mathrm{MF}}|=0$$

$$E_{\mathrm{FM}}=\frac{1}{2}\left(\frac{J_{\mathrm{d}}}{z}+J_{\parallel }+J_{\times }\right)+\frac{1}{16}{J}_4.$$

$${\mathbf{N}}_{\mathrm{A}}^{\mathrm{MF}}={\mathbf{N}}_{\mathrm{B}}^{\mathrm{MF}},|{\mathbf{T}}_{\Lambda }^{\mathrm{MF}}|=0,|{\mathbf{N}}_{\Lambda }^{\mathrm{MF}}|>0$$

$${\mathbf{T}}_{\mathrm{A}}^{\mathrm{MF}}=-{\mathbf{T}}_{\mathrm{B}}^{\mathrm{MF}},|{\mathbf{T}}_{\Lambda }^{\mathrm{MF}}|=1,|{\mathbf{N}}_{\Lambda }^{\mathrm{MF}}|=0$$

$$E_{\text{C-AFM}}^{\mathrm{MF}}=\frac{1}{2}\left(\frac{J_{\mathrm{d}}}{z}-J_{\parallel }-J_{\times }\right)+\frac{1}{16}{J}_4.$$

$$|{\mathbf{T}}_{\Lambda }^{\mathrm{MF}}|=0,|{\mathbf{N}}_{\Lambda }^{\mathrm{MF}}|=0,$$

$${\mathbf{Q}}_{\mathrm{A}}^{\mathrm{MF}}={\mathbf{Q}}_{\mathrm{B}}^{\mathrm{MF}},|{\mathbf{Q}}_{\Lambda }^{\mathrm{MF}}|=\sqrt{\frac{4}{3}}$$

$$E_{\mathrm{SNf}}^{\mathrm{MF}}=\frac{1}{2}\frac{J_{\mathrm{d}}}{z}+\frac{3}{16}{J}_4.$$

$$|{\mathbf{T}}_{\Lambda }^{\mathrm{MF}}|>0,|{\mathbf{N}}_{\Lambda }^{\mathrm{MF}}|=0$$

$${\mathbf{T}}_{\Lambda }^{\mathrm{MF}}\cdot {\mathbf{T}}_{{\Lambda }^{\prime }}^{\mathrm{MF}}=-\frac{1}{2}|{\mathbf{T}}_{\Lambda }^{\mathrm{MF}}||{\mathbf{T}}_{{\Lambda }^{\prime }}^{\mathrm{MF}}|$$

Full text

I Introduction

The search for spin nematic states has been under active investigation in the last decade. The spin nematic state is characterized by the absence of any magnetic long-range order (except a trivial magnetization in an applied magnetic field) and the spontaneous symmetry breaking of spin rotation accompanied with a spin-quadrupolar long-range order.Andreev and Grishchuk Because of these peculiar properties, the spin nematic state is a novel, intriguing non-magnetic state with a hidden order.

In theoretical studies, there are already many proposals of the spin models that exhibit the spin nematic order at low temperatures. A typical example in such models is the spin-1 bilinear-biquadratic model,Blume and Hsieh (1969); Chen and Levy (1971); Papanicolaou (1984) which includes bilinear (two-spin) and biquadratic (four-spin) exchange interactions. The appearance of the spin nematic phase on the cubic lattice was shown for sufficiently large biquadratic interactions by a rigorous proof Tanaka et al. (2001) and quantum Monte Carlo simulations.Harada and Kawashima (2002) It has been also shown that the spin-1 bilinear-biquadratic model exhibits the spin-nematic ground-state phases with ferro- or antiferro-quadrupolar order on various two-dimensional lattices.Tóth et al. (2012); Niesen and Corboz (2017); Zhao et al. (2012); Läuchli et al. (2006); Liu et al. (2015) Another example of the models that show the spin nematic phase is a family of spin-1/2 frustrated ferromagnets which include ferromagnetic (FM) exchange interactions and competing antiferromagnetic (AFM) ones.Shannon et al. (2006); Shindou et al. (2011); Chubukov (1991); Hikihara et al. (2008); Sudan et al. (2009); Zhitomirsky and Tsunetsugu (2010); Ueda and Momoi (2013); Janson et al. (2016) In this case, the instability leading to the spin nematic phase appears in the saturated state in an external field, where two magnons form a bound state. When the external field decreases below the saturation field, those bound magnon pairs condense, which leads to spin nematicity.Momoi and Shannon (2005); Shannon et al. (2006) There are also a few other examples which show the spin-nematic ground state in anisotropic spin modelsSchulz (1986); Chen et al. (2003); Damle and Senthil (2006), the Heisenberg model on Shastry-Sutherland lattice,Momoi and Totsuka (2000); Wang and Batista (2018) and models with multiple-spin ring exchanges.Momoi and Shannon (2005); Hikihara and Yamamoto (2008); Momoi et al. (2012); Yokoyama and Hotta (2018)

In real materials, however, candidate materials for the spin nematic state are rather limited. The biquadratic interaction in spin-1 systems is not very strong in general, but the realization of the spin nematic state in the spin-1 systems requires relatively strong biquadratic interaction comparable to or larger than the bilinear exchange interaction. We note that there are some proposals to enhance the ratio of the biquadratic interaction to the bilinear one.Mila and Zhang (2000); Tanaka et al. (2018) In the case of spin-1/2 frustrated ferromagnets, the spin nematic phase can appear only in a narrow parameter range at zero field and it appears in a wider parameter space in a high magnetic field,Sindzingre et al. (2009) which is not easy to access in experiments. Nevertheless, active studies on spin-1/2 frustrated ferromagnets are on going, for example, in the quasi-one-dimensional spin-1/2 frustrated ferromagnets LiCuVO4Nawa et al. (2017); Orlova et al. (2017) and Rb2Cu2Mo3O12,Matsui et al. (2017) and a two-dimensional kagome compound, volborthite Cu3V2O7(OH)${}_{2}\cdot$2H2O.Yoshida et al. (2017); Kohama et al. (2019) For stimulating further studies, it is desirable to search still more theoretical models for describing the spin nematic phase at zero field in the systems with only two-spin interactions, which are easier to access experimentally.

Theoretical studies on spin nematics in spin-1/2 systems often face technical difficulties. The appearance of spin nematic phases is usually understood from the bi-magnon instability in the fully polarized state,Shannon et al. (2006); Kecke et al. (2007); Momoi et al. (2012) where some of the induced spin nematic phases are expected to remain down to low magnetization regime. Nevertheless, even in the spin-1/2 square-lattice $J_{1}$-$J_{2}$ model with ferromagnetic $J_{1}$, which is thought to be a typical model for spin nematic ordering,Shannon et al. (2006) there still exists a certain controversy if the spin nematic order remains at zero field.Richter et al. (2010); Shindou et al. (2011) For a better understanding of spin nematic ordering at zero field, it is hence helpful to find some specific spin-1/2 models which definitely show spin nematic phases at zero magnetic field. In this context, the mechanism of spin nematics in spin-1 systems is instructive. In the spin-1 bilinear-biquadratic model, it is known that an SU(3) symmetry plays an important role in the emergence of spin nematic phases at zero field.Papanicolaou (1988); Penc and Läuchli (2011) Referring to this, we also take into account the effect of four-spin interactions in our spin-1/2 systems, to elucidate mechanism of spin nematic ordering. We will see that an SU(4) symmetry appearing in our model at a certain parameter point serves as an origin of the spin-nematic phase.

The aim of this work is twofold: One is to propose a spin-1/2 model that requires neither the four-spin interaction nor magnetic field for realizing the spin nematic state. The other is to clarify the mechanism of this spin nematic ordering by considering four-spin interactions as well. For these purposes, we introduce in this paper a spin-1/2 frustrated ferromagnetic model in two dimensions. The model consists of ferromagnetically coupled dimers of $S=1/2$ spins forming a bilayer structure and contains frustrated bilinear Heisenberg exchange interactions and a type of four-spin interaction [see Eq. (1) and Fig. 1 in Sec. II]. The parameter space of the model includes the effective spin-1 bilinear-biquadratic model and an SU(4) symmetric model, which enables us to show the mechanism of spin nematics. Using analytical and numerical techniques, we obtain the following results.

First, we analyze the strong ferromagnetic dimer limit. Performing a perturbative calculation, we map the model to the spin-1 bilinear-biquadratic model. With the help of former studies on spin-1 systems, we can determine the ground-state phase diagrams on various lattices, which include the spin nematic phases in wide parameter regions. In particular, we show that, even when the four-spin interaction is absent, a second-order perturbation yields an effective biquadratic interaction, which leads to the appearance of the spin nematic phase with ferro-quadrupolar order in a parameter regime where all the first-order perturbations are canceled to each other.

Second, we study how the spin nematic phase emerging at the strong dimer limit is affected by decreasing the ferromagnetic intradimer coupling from infinite to finite values. We use a mean-field approximation with product-state ansatz and the many-variable variational Monte Carlo (mVMC) method.Misawa et al. (2019); mVM In the mean-field approximation, we determine the ground-state phase diagram as a function of two-spin and four-spin interactions. In addition to conventional magnetically-ordered phases, the spin nematic phase with the ferro-quadrupolar order appears in a wide parameter region. In this mean-field approximation, the spin nematic phase vanishes when the four-spin interaction is deleted. The result also reveals that this spin nematic phase connects with the SU(4) symmetric point in our model. To further take account of the effects of interdimer correlations, we study the model with only the two-spin interactions using mVMC method. The obtained phase diagrams for the square and triangular lattices show that, for large but finite intradimer interactions, the spin nematic phase emerges in a finite parameter range. This result is consistent with our aforementioned discussion from the second-order perturbation and thereby confirms that our model exhibits the spin nematic phase even when it does not include any four-spin interaction or external field.

Lastly, we focus on the model around the SU(4) symmetric point. From an exact symmetry argument, we find that our model at the SU(4) point has degenerate ground states including the spin nematic state and the vector chiral (also known as p-type nematicAndreev and Grishchuk ) state. It also shows that our mean-field approximation indeed presents exact ground states in the SU(4) model. Thus, the SU(4) model is a source of various exotic phases and one can realize each of them by adding perturbative interactions. As an example of such a perturbation, we show that an appropriate set of Ising-type interactions can yield the spin nematic phase.

The rest of the paper is organized as follows: In Sec. II, we introduce the model Hamiltonian, and also discuss symmetric properties of the model and the order parameters studied. In Sec. III, we present the results of the perturbation calculations from the strong ferromagnetic dimer limit. We further present the results of the mean-field approximation with a product state and the mVMC method in Secs. IV and V, respectively. In Sec. VI we discuss the exact arguments on emerging phases in the vicinity of the SU(4)-symmetric model and also on the effect of perturbations of Ising interactions. Section VII is devoted to summary and concluding remarks. Details of the mean-field calculation, SU(4) transformation on dimers, and nontrivial degeneracies in the mean-field solutions are discussed in Appendices A, B, and C, respectively.

II Model

Here we introduce the spin model studied in this paper. We also present all local observables which concern us and symmetries inherent in our model.

II.1 Hamiltonian

We study the spin-1/2 frustrated quantum magnet consisting of spin-dimer units, which form a two-dimensional bilayer structure of spins. The model Hamiltonian has the form

[TABLE]

where ${\bm{S}}_{l,j}$ is the spin-1/2 operator of the $l$th spin ($l=1,2$) in the $j$th dimer. The dimer sites, labeled with $j$ or $j^{\prime}$, form two-dimensional lattices and the sum $\sum_{\langle j,j^{\prime}\rangle}$ is taken for the nearest-neighboring sites in the lattices. We consider several lattices including the square, honeycomb, triangular, and kagome lattices. Schematic pictures of the model are shown in Fig. 1. Throughout this paper we consider the case that the intradimer exchange interaction is ferromagnetic or zero, $J_{\rm d}\leq 0$, where two $S=1/2$ spins in each dimer dominantly form a spin triplet.

II.2 Local observables

We explore local order parameters measured with the following local observables: For a parallel spin order on a dimer bond, we use the total spin operators on each dimer,

[TABLE]

with $\alpha=x,y,z$ and, for an antiparallel spin order, Néel-spin operators

[TABLE]

For a spin quadrupolar order on a dimer bond, we use the five-component spin quadrupolar operators

[TABLE]

which act on two $S=1/2$ spins on each dimer.

These quadrupolar operators are a natural extension of the on-site quadrupolar operators in spin-1 systems; one can obtain the above operators by inserting the total spin operators $T^{\alpha}_{j}$ into the spin-1 operators in the quadrupolar operators of spin-1 systems. (For the definition of the spin-1 quadrupolar operators, see for example Ref. Penc and Läuchli, 2011.) This derivation of Eq. (4) readily concludes that the commutation relations between the total-spin operators $T_{j}^{\alpha}$ and the quadrupolar operators $Q_{j}^{(n)}$ are the same as those in spin-1 systems. Hence the combined set of $T_{j}^{\alpha}$ ($\alpha=x,y,z$) and $Q_{j}^{(n)}$ ($n=1,\cdots,5$) forms SU(3) group; $T_{j}^{\alpha}$ and $Q_{j}^{(n)}$ are eight-dimensional generators of su(3) algebra.

Incidentally, one can obtain the quadrupolar operators (4) with a minus sign by substituting the Néel-spin operators $N^{\alpha}_{j}$ into the $S=1$ spin operators in the quadrupolar operators of spin-1 systems. From this fact, it also follows that the combined set of $N_{j}^{\alpha}$ ($\alpha=x,y,z$) and $-Q_{j}^{(n)}$ ($n=1,\cdots,5$) also forms another SU(3) group.

II.3 Symmetries

In addition to the apparent SU(2) symmetry, the model (1) has the following higher symmetries in specific parameter spaces.

*SU(3) symmetry: * In the parameter space defined by

[TABLE]

with any $J_{\rm d}$, the total Hamiltonian (1) has a global SU(3) symmetry. The eight generators of SU(3) rotation are given by $\sum_{j}T_{j}^{\alpha}$ ($\alpha=x,y,z$) and $\sum_{j}Q_{j}^{(n)}$ ($n=1,\cdots,5$). In particular, in the $J_{\rm d}\to-\infty$ limit, the model (1) in the broader parameter space

[TABLE]

also has the same SU(3) symmetry. We note that this space (6) includes the aforementioned space (5). In this limit, as the singlet state is gapped out in each dimer, the model is reduced to a spin-1 model. The model in the space (6) is mapped to the SU(3) symmetric spin-1 bilinear-biquadratic model (see Sec. III).

SU(4) symmetry: In the case of $J_{\rm d}=0$ among the SU(3) symmetric space (5), i.e.,

[TABLE]

the total Hamiltonian has a global SU(4) symmetry; $\mathcal{H}$ commutes with the operators $\sum_{j}S^{\alpha}_{1,j}$, $\sum_{j}S^{\alpha}_{2,j}$ ($\alpha=x,y,z$), and $\sum_{j}S^{\alpha}_{1,j}S^{\beta}_{2,j}$ ($\alpha,\beta=x,y,z$), which are known as the fifteen generators of SU(4) group. The SU(3) group mentioned above is a subgroup of this SU(4) group. This SU(4) symmetry also contains another SU(3) symmetry generated by $\sum_{j}N_{j}^{\alpha}$ and $-\sum_{j}Q_{j}^{(n)}$. We present our analysis on the SU(4) symmetric model with ferromagnetic coupling ($J_{\parallel}=J_{4}/4<0$) in Sec. VI and further describe SU(4) transformation in Appendix B.

We will see in the following sections that the above high symmetric models in a ferromagnetic coupling regime are on an exact phase boundary or a multiple point in the ground-state phase diagram. The SU(3) symmetric model in the space (5) with ferromagnetic couplings $J_{\rm d}<0$ and $J_{\parallel}<0$ is on an exact phase boundary between the ferromagnetic phase and the spin nematic phase with ferroquadrupolar order. Similarly, in the $J_{d}\rightarrow-\infty$ limit, the SU(3) model with the parameters (6) under the ferromagnetic condition $J_{\parallel}+J_{\times}+\frac{J_{4}}{4}<0$ is on the same phase boundary, which is shown in the next section. The SU(4) symmetric model given by Eq. (7) in the ferromagnetic case $J_{\parallel}<0$ is at a multiple point where many phases coexist, which is further discussed in Sec. VI.

III Strong ferromagnetic-dimer limit

In this section, we study the model (1) in the limit of strong ferromagnetic intradimer coupling, $J_{\rm d}\to-\infty$. Treating the intradimer exchange term $\mathcal{H}_{\rm d}$ as a unperturbed Hamiltonian and the rest of terms as a perturbation, we derive an effective Hamiltonian, which enables us to see the mechanism of spin nematic ordering.

In the ground state of the unperturbed Hamiltonian $\mathcal{H}_{\rm d}$, two $S=1/2$ spins in each dimer form a spin triplet and the ground states are $3^{N}$-fold degenerate, where $N$ is the number of dimers in the system. The first-order perturbation induces state transitions between degenerate ground states, whose matrix elements are written with the effective Hamiltonian

[TABLE]

where $\tilde{\bm{S}}_{j}$ denote the spin-1 operators acting on the $S=1$ triplet sector on the $j$th dimer. This first-order perturbation Hamiltonian is nothing but the spin-1 bilinear-biquadratic model.

The spin-1 bilinear-biquadratic model has been extensively studied on various lattices. In Table 1, we summarize the obtained ground-state phases. Here, defining the Hamiltonian

[TABLE]

with the parameter $\theta$ and $J_{\rm bb}>0$, we describe the phase diagrams as functions of $\theta$. For the square lattice the phase diagram contains at least five phases, i.e., the ferromagnetic (FM), ferro-quadrupolar (FQ), Néel, three-sublattice antiferromagnetic (AFM3), and three-sublattice antiferro-quadrupolar (AFQ3) phases.Tóth et al. (2012) The emergence of a quasi-one-dimensional Haldane phase in a narrow region between the Néel and AFM3 phases was also reported.Niesen and Corboz (2017) The phase diagram for the honeycomb lattice includes the FM, FQ, Néel, and plaquette valence-bond-solid phases.Zhao et al. (2012) For the triangular lattice, the phase diagram contains the FM, FQ, 120∘-structure antiferromagnetic (120∘-AFM), and AFQ3 phases.Läuchli et al. (2006) For the kagome lattice, the phase diagram was found to include the FM, FQ, antiferro-quadrupolar, and trimerized valence-bond-crystal phases.Liu et al. (2015) The parameter ranges of $\theta$ for these phases are shown in Table 1 for each lattice. It is noteworthy that the FQ phases in the geometrically-frustrated (triangular and kagome) lattices appear in wider regions than those in the bipartite (square and honeycomb) lattices because of the suppression of antiferromagnetic ordering in the former lattices.

From these results, we can derive the phase diagram for the first-order perturbation Hamiltonian ${\cal H}^{(1)}$ using the relation $J_{\rm bb}\cos\theta=(J_{\parallel}+J_{\times})/2+J_{4}/8$ and $J_{\rm bb}\sin\theta=J_{4}/4$. The resultant regions of each phase are also presented in Table 1 for the lattices considered. Again, we find that the regions of the FQ phases in the triangular and kagome lattice systems are wider than those in the square and honeycomb lattice systems. The phase diagrams in the first-order perturbation for the square and triangular lattices are shown in the $J_{\parallel}+J_{\times}$ versus $J_{4}$ plane in Fig. 2. The quadrupolar phases appear in between the ferromagnetic phase and antiferromagnetic phases. The parameter regions of the FQ and AFQ3 phases appearing for negative and positive $J_{4}$, respectively, shrink with decreasing $|J_{4}|$ and vanish at $J_{4}=0$ (within the first-order perturbation).

In the case of $J_{4}=0$, the first-order perturbation Hamiltonian contains only the two-spin exchange interactions and they vanish at $J_{\parallel}+J_{\times}=0$. In this situation, the second-order perturbation, which induces further effective interactions, becomes relevant. From a standard procedure of perturbation theory, the second-order perturbation Hamiltonian in the case of $J_{\parallel}+J_{\times}=J_{4}=0$ turns out to have a rather simple formHikihara and Starykh (2010)

[TABLE]

with

[TABLE]

Note that the coupling constant $J^{(2)}$ is always negative.

The second-order perturbation Hamiltonian (10) leads us to an important conclusion for the case of $J_{4}=0$. It is natural to expect that, for strong ferromagnetic $J_{\rm d}$, the biquadratic interaction in Eq. (10) is still dominant over other interactions in a finite parameter region around $J_{\parallel}+J_{\times}=0$. As shown in Table 1, the pure-biquadratic model (10), which is equivalent to the spin-1 bilinear-biquadratic model (9) with $\theta=-\pi/2$, has the ferro-quadrupolar ground state on the triangular and kagome lattices. We hence conclude that, in our model (1) in the vicinity of $J_{\parallel}+J_{\times}=0$ with strong ferromagnetic $J_{\rm d}$, the dominant effective biquadratic interaction leads to the ferro-quadrupolar ground state on the geometrically frustrated lattices even in the case of $J_{4}=0$, in which the original model includes only the bilinear exchange terms. We will confirm this conclusion numerically in Sec. V.

For the case of the square and honeycomb lattices, the spin-1 biquadratic Hamiltonian (10) is just on the phase boundary between the ferro-quadrupolar and Néel ordered phases as shown in Table 1, where the ferro-quadrupolar phase spreads to a finite region with ferromagnetic bilinear interactions. If we slightly shift the couplings $J_{\parallel}$ and $J_{\times}$ from the phase boundary $J_{\parallel}+J_{\times}=0$ into the region $J_{\parallel}+J_{\times}<0$, they yield a ferromagnetic bilinear interaction between the effective $S=1$ spins $\tilde{\bm{S}}_{j}$ and $\tilde{\bm{S}}_{j^{\prime}}$ due to the first-order perturbation. We hence expect that even when $J_{4}=0$, our original model (1) on these bipartite lattices realizes the ferro-quadrupolar phase in a finite parameter region in $J_{\parallel}+J_{\times}<0$. We will confirm in Sec. V that this is also the case.

IV Mean-Field Approximation with Product-State Ansatz

In this section, we employ a mean-field approximation with product-state ansatz to determine the ground-state phase diagram of the model (1). We consider the case of $J_{4}\leq 0$ and $J_{\rm d}\leq 0$. Some details of the method and results are presented also in Appendix A.

IV.1 Method

We employ the approximation in which the ground-state wave function is expressed by a direct product of dimer states,

[TABLE]

where the dimer states $|\varphi\rangle_{j}$ can take an arbitrary state spanned with the dimer bases. We further assume that the wave function has two- and three-sublattice structures, respectively, for the bipartite (square and honeycomb) lattice systems and the triangular lattice system; namely, the dimers in the same sublattice are in the same state,

[TABLE]

for $j\in\Lambda$, where $\Lambda={\rm A,B}$ $({\rm A,B,C})$ denotes the two sublattices (three sublattices) and $|\sigma_{1}\sigma_{2}\rangle$ denotes the dimer state with the eigenvalues $S^{z}_{1,j}=\sigma_{1}$ and $S^{z}_{2,j}=\sigma_{2}$. Optimizing the coefficients $a_{\Lambda,\sigma_{1}\sigma_{2}}$ in Eq. (13) variationally, we obtain the lowest-energy mean-field solution.

To obtain the ground state, we minimize the expectation value of the bond Hamiltonian of the model (1),

[TABLE]

where $z$ is the coordination number. We note that the coordination number $z$ is taken into account only through the coupling constant $J_{\rm d}/z$ in Eq. (14). Using the resultant ground state, we calculate the expectation values of local observables defined in Sec. II.2,

[TABLE]

for any $j\in\Lambda$ for each sublattice $\Lambda$.

The similar mean-field approximation with site-decoupled wave functions has been applied to the spin-1 bilinear-biquadratic model (9) on the square, honeycomb, and triangular lattices. This approximation provides pretty accurate results for $\pi/2<\theta<2\pi$: For the squareTóth et al. (2012) and honeycombZhao et al. (2012) lattices, the mean-field approximation yields the phase diagram for $\pi/2<\theta<2\pi$ which is completely the same as those obtained by other numerical approaches such as exact diagonalization and tensor renormalization group technique. For the triangular lattice,Läuchli et al. (2006) the phase diagram obtained by the mean-field approximation is essentially the same as that by the exact diagonalization; the only discrepancy appears in the phase boundary between the ferro-quadrupolar and 120∘-AFM phases, where the ferro-quadrupolar phase region becomes narrower in the mean-field approximation. On the other hand, for $0<\theta<\pi/2$, the mean-field approximation is rather unreliable at least for the bipartite lattices since the direct-product wave function is not able to describe the Haldane phase on the square lattice and the plaquette valence-bond-crystal phase on the honeycomb lattice, in which the entanglement between different dimers is essential. In our calculation, we hence restrict ourselves to explore the parameter region of $J_{4}\leq 0$, which corresponds, in the limit $J_{\rm d}\to-\infty$, to the region of $\pi\leq\theta\leq 2\pi$, where the mean-field approximation is expected to be reliable. In the following, setting

[TABLE]

we determine the ground-state phase diagram in $J_{\times}$ versus $J_{4}$ planes with $J_{4}\leq 0$ for several fixed values of $J_{\rm d}/z$.

IV.2 Two-sublattice case

First, we discuss the case of the two-sublattice structure. From the expectation values of local observables we studied, we found four distinct phases. These phases are characterized with the expectation values ${\bm{T}}^{\rm MF}_{\Lambda}$, ${\bm{N}}^{\rm MF}_{\Lambda}$, and ${\bm{Q}}^{\rm MF}_{\Lambda}$ on two sublattices $\Lambda={\rm A},{\rm B}$, which are summarized as follows:

• (i)

Ferromagnetic (FM) phase: All spins are fully polarized, pointing to the same direction,

[TABLE]

for $\Lambda={\rm A},{\rm B}$. The ground-state energy per bond of this phase is given by

[TABLE]

• (ii)

A-type antiferromagnetic (A-type AFM) phase: Two spins in each dimer are antiparallel to each other, and all of the staggered moments $\langle{\bm{N}}_{j}\rangle$ are in the same direction,

[TABLE]

for $\Lambda={\rm A},{\rm B}$. This state can be also regarded as two ferromagnetic layers whose moments are antiparallel to each other. The spin moments $\langle{\bm{S}}_{l,j}\rangle$ may shrink, i.e., $|{\bm{N}}^{\rm MF}_{\Lambda}|\leq 1$, due to quantum fluctuation.

• (iii)

C-type antiferromagnetic (C-type AFM) phase: Two spins in each dimer point to the same direction, and the total spin moments ${\bm{T}}^{\rm MF}_{\rm A}$ and ${\bm{T}}^{\rm MF}_{\rm B}$ are antiparallel, forming the Néel-type magnetic order in $\langle{\bm{T}}_{j}\rangle$,

[TABLE]

for $\Lambda={\rm A},{\rm B}$. Each local spin is fully polarized. The ground-state energy per bond of this phase (in the mean-field approximation) is

[TABLE]

• (iv)

Spin nematic phase with ferro-quadrupolar order: All the spin-dipole moments vanish, while the spin-quadrupolar moments take the same finite value on all dimers,

[TABLE]

for $\Lambda={\rm A},{\rm B}$. The quadrupolar moments are saturated. The ground-state energy per bond of this phase (in the mean-field approximation) is

[TABLE]

In the limit $J_{\rm d}\to-\infty$, this phase corresponds to the ferro-quadrupolar phase of the spin-1 bilinear-biquadratic model.

Schematic illustration of spin structures representing these phases are shown in Figs. 3(a)–(d). We note that the fully-saturated nature of the C-type AFM and spin nematic phases is an artifact of the approximation with product-state ansatz. Indeed, we will show in Sec. V that, in the mVMC calculations, the quantum reduction in the magnetic and spin-quadrupolar moments is observed also in these phases.

We have determined the phase diagrams for several values of $J_{\rm d}/z$. Figure 4 shows the results for some typical values of $J_{\rm d}/z$.

At $J_{\rm d}=0$, we find three regions; two are the FM and A-type AFM phases and the other corresponds to the boundary between the spin nematic phase, which appears for $J_{\rm d}<0$, and dimer-singlet phase for $J_{\rm d}>0$. The quadruple point where these four phases coexist is the SU(4) symmetric point, given by $J_{\parallel}=-1$, $J_{4}=-4$, and $J_{\rm d}=J_{\times}=0$. We further show in Sec. VI, using exact symmetry arguments on a generic model, that the spin nematic state is naturally generated by SU(4) symmetry and the SU(4) model is on the multiple point surrounded by, at least, five phases including the spin nematic phase. The phase boundary between the FM and A-type AFM phases is on the line $J_{\times}=0$ for $J_{4}>-4$, while the region of the phase boundary between the spin-nematic and dimer-singlet phases is surrounded by the two lines $J_{4}=\pm 4J_{\times}-4$. On the boundary between the spin-nematic and dimer-singlet phases, there exists additional non-trivial degeneracy in the mean-field solution, which is a remnant of ${\rm SU}(2)\times{\rm SU}(2)$ symmetry in the case of $J_{\rm d}=J_{\times}=0$. (See the Appendix C.2 for more details.)

When $J_{\rm d}$ is negative, the spin nematic state with ferro-quadrupolar order is selected from the non-trivially degenerate ground states in the degenerate region, resulting in the spin nematic phase defined in Eq. (22). Another signature of high SU(4) symmetry remains on the boundary between the FM and spin-nematic phases for $J_{\rm d}<0$ in the mean-field approximation, where the boundary line $J_{4}=4J_{\times}-4$ is obtained from the condition $E_{\rm FM}=E^{\rm MF}_{\rm SNf}$. On this boundary, the mean-field solution of the ground state has non-trivial SU(3) degeneracy, which we further explain in Appendix C.1.

As $|J_{\rm d}|/z$ increases, the region of the spin nematic phase enlarges toward smaller $|J_{4}|$ regime, and the boundary between the FM and A-type AFM phases also moves toward large $J_{\times}$ [see Figs. 4(b) and (c)]. Then, at $J_{\rm d}/z\sim-1.15$, the C-type AFM phase enters the parameter space calculated [see Fig. 4(d)]. The appearance of the C-type AFM phase is due to the competition between the exchange interactions $J_{\rm d}$ and $J_{\parallel}$. Indeed, at the limit of $J_{\times}\to\infty$ and $J_{4}=0$, the interactions $J_{\rm d}$ and $J_{\parallel}$ lead to the A-type AFM phase for $|J_{\rm d}|/z<|J_{\parallel}|=1$ and the C-type AFM phase for $|J_{\rm d}|/z>|J_{\parallel}|$.

With further increasing $|J_{\rm d}|/z$, we see that the A-type AFM phase shrinks while the other phases enlarge. [See Fig. 4(e).] The A-type AFM phase eventually vanishes at $J_{\rm d}/z\simeq-2$ and the spin nematic phase touches with the C-type AFM phase. [See Fig. 4(f).] The spin nematic phase extends to the line of $J_{4}=0$, but touches to the line $J_{4}=0$ only at the single point $J_{\times}=1$. Hence in the case of $J_{\rm d}/z\leq-2$ and $J_{4}=0$, the phase diagram contains only the FM and C-type AFM phases. For $J_{\rm d}/z\leq-2$, the phase diagram, which is unaltered at least down to $J_{\rm d}/z=-10$ in our calculation, is the same as the one obtained for the limit $J_{\rm d}\to-\infty$ in Sec. III. The boundary lines of the spin nematic phase are $J_{4}=\pm 4(J_{\times}-1)$, which are obtained from the conditions $E_{\rm FM}=E^{\rm MF}_{\rm SNf}$ and $E^{\rm MF}_{\text{C-AFM}}=E^{\rm MF}_{\rm SNf}$.

Among the transitions occurring in the present parameter space, the transition between the spin-nematic and A-type AFM phases is continuous. In contrast, the other transitions are always of first order, accompanied with a jump in ${\bm{T}}^{\rm MF}_{\Lambda}$. (See the Appendix A.2 for detailed analysis of the order of the transitions.) In Fig. 5, we plot $J_{\times}$-dependence of the order parameters on the parameter lines given by $J_{4}=-0.5$ and $-4.0$ in the plane with $J_{\rm d}/z=-1.2$ [dotted lines of Fig. 4(d)] and the parameter line given by $J_{4}=-3.0$ in the plane with $J_{\rm d}/z=-2.0$ [dotted line in Fig. 4(f)]. In the case of $J_{\rm d}/z=-1.2$ and $J_{4}=-0.5$ [Fig. 5(a)], the system undergoes two successive transitions from the FM phase to the A-type AFM phase, and then, to the C-type AFM phase as $J_{\times}$ increases. At both transitions, the order parameters exhibit finite jumps, showing first-order phase transitions. In the intermediate A-type AFM phase, the Néel-spin moment ${\bm{N}}^{\rm MF}_{\Lambda}$ shrinks from the saturated value due to the zero-point quantum reduction, while the total spin ${\bm{T}}^{\rm MF}_{\Lambda}$ on each dimer is fully polarized in the FM and C-type AFM phases. When negative $J_{4}$ becomes stronger [see Fig. 5(b)], the system exhibits four phases, the FM, spin-nematic, A-type AFM, and C-type AFM phases. It is found that the transition between the spin-nematic and A-type AFM phases is continuous while the other two transitions are of first order. The Néel-spin moment exhibits the quantum reduction in the A-type AFM phase and vanishes continuously at the boundary to the spin nematic phase. For $J_{\rm d}/z=-2.0$ and $J_{4}=-3.0$ [Fig. 5(c)], there occur two first-order transitions from the FM phase to the spin nematic phase, and then, to the C-type AFM phase with increasing $J_{\times}$. Finally, we note that the transition between the FM and C-type AFM phases occurring on the line of $J_{4}=0$ is of first order.

IV.3 Three-sublattice case

Next, we discuss the case of three-sublattice structure on the triangular lattice. Figure 6 presents the phase diagrams we obtained for various $J_{\rm d}/z$. We found four distinct phases; three of them are the FM, A-type AFM, and spin-nematic phases, which are translationally invariant and defined in the same manners as those for the two-sublattice case, i.e., Eqs. (17), (19), and (22), respectively. The other phase is as follows:

• (iii’)

C-type 120∘-structure antiferromagnetic (C-type 120∘-AFM) phase: Two spins in each dimer are parallel to each other,

[TABLE]

for different sublattices $\Lambda$ and $\Lambda^{\prime}$, which is the well-known spin structure in the triangular-lattice antiferromagnet.Miyashita and Shiba (1984); Kawamura and Miyashita (1984) The local spins shrink due to quantum fluctuation.

A schematic picture representing this phase is shown in Fig. 3(e).

For several values of $J_{\rm d}/z$, we have determined the ground-state phase diagrams, some of which are shown in Fig. 6. For small $|J_{\rm d}|/z$, the resultant diagrams are the same as those for the two-sublattice case. At $J_{\rm d}=0$, the diagram includes the FM phase, the A-type AFM phase, and the region with the degenerate ground states which is the phase boundary between the spin-nematic and dimer-singlet phases. The phase boundaries are given by $J_{4}=\pm 4J_{\times}-4$ ($J_{4}<-4$) and $J_{\times}=0$ ($J_{4}>-4$). The SU(4) point ($J_{\parallel}=-1$, $J_{4}=-4$, $J_{\rm d}=J_{\times}=0$) is on the quadruple point, where the spin nematic phase is generated by SU(4) symmetry as shown in Sec. VI.1. When a negative $J_{\rm d}$ is introduced, the spin nematic phase with ferro-quadrupolar order replaces the region of the degenerate ground states. As $|J_{\rm d}|/z$ increases, the FM and spin-nematic phases enlarge, while the A-type AFM phase becomes smaller. Then, at $J_{\rm d}/z\sim-1.7$, the C-type 120∘-AFM phase appears from the region around $1\lesssim J_{\times}\lesssim 2$ and $J_{4}=0$. The appearance of this phase can be understood from the strong coupling limit $J_{\rm d}\rightarrow-\infty$, where 120∘-AFM phase appears in large $J_{\times}$ regime. As $|J_{\rm d}|/z$ further increases, the spin nematic and C-type 120∘-AFM phases extend and eventually, around $J_{\rm d}/z\sim-3.5$, cover the A-type AFM phase region. For $J_{\rm d}/z\leq-3.5$, the phase diagram in our scope ($-2\leq J_{\times}\leq 4$ and $-6\leq J_{4}\leq 0$) does not depend on $J_{\rm d}/z$, at least down to $J_{\rm d}/z=-10$ in our calculation. This phase diagram coincides with the diagram obtained by the mean-field approximationLäuchli et al. (2006) to the spin-1 bilinear-biquadratic model (8) derived with the first-order perturbation in Sec. III on the triangular lattice. The spin nematic phase boundary touches to the $J_{4}=0$ line on the single point $J_{\times}=1$. The two boundary lines are given by $J_{4}=4(J_{\times}-1)$ and $J_{4}=-2(J_{\times}-1)$.

We also show the order of phase transitions in Fig. 6. The transitions from the FM phase to the other phases are of first order, accompanied with jumps in ${\bm{T}}^{\rm MF}_{\Lambda}$. The transitions from the spin nematic phase to the two types of AFM phases, i.e., A-type AFM and C-type 120∘-AFM phases, are continuous while the one between these two AFM phases is of first order. (See also Appendices A.2–A.4 for detailed results.) Figure 7 presents the $J_{\times}$-dependence of the order parameters on the parameter lines given by $J_{4}=-1.0$ and $-4.0$ in the plane with $J_{\rm d}/z=-2.3$ [dotted lines in Fig. 6(d)] and the line given by $J_{4}=-4.0$ in the plane with $J_{\rm d}/z=-3.5$ [dotted line in Fig. 6(f)]. For $J_{\rm d}/z=-2.3$ and $J_{4}=-1.0$ [Fig. 7(a)], there appear four phases, the FM, spin-nematic, C-type 120∘-AFM, and A-type AFM phases. The magnetic moments exhibit the quantum reduction in the C-type 120∘-AFM and A-type AFM phases and vanish in the spin nematic phase. The order parameters exhibit jumps at the phase transitions between the FM and spin-nematic phases and between the C-type 120∘-AFM and A-type AFM phases, while they change continuously at the transition between the spin-nematic and C-type 120∘-AFM phases. The continuous nature of the latter transition is also confirmed clearly in the results for $J_{\rm d}/z=-3.5$ and $J_{4}=-4.0$ in Fig. 7(c). On the parameter line of $J_{\rm d}/z=-2.3$ and $J_{4}=-4.0$ [Fig. 7(b)], the system undergoes two transitions from the FM phase to the spin nematic phase, and then to the A-type AFM phase. We find again that the transition between the FM and spin-nematic phases is of first order, while the transition between the spin-nematic and A-type AFM phases is continuous.

IV.4 Comments

Three comments are in order here. First, in our result by the mean-field approximation with the product-state ansatz, we find no finite region of the spin nematic phase in the case of $J_{4}=0$, contrary to the expectation from the perturbation theory in Sec. III. This result is attributed to the fact that the product-state ansatz, which completely ignores the entanglement between different dimers, is not able to include the effect of the second-order perturbation sufficiently. Indeed, as shown in the subsequent section, the mVMC method provides a result that the spin nematic phase emerges in a finite region in the $J_{4}=0$ case.

Second, it has been shown that, when applied to the spin-1 bilinear-biquadratic model (9) on the triangular lattice, the mean-field approximation yields an inaccurate result on the phase boundary between the spin-nematic and 120∘-AFM phases; the mean-field approximation gives the critical point $\theta_{\rm c}=\arctan(-2)\sim-0.35\pi$, while the exact diagonalization method, which can take the entanglement between different sites into account, gives $\theta_{\rm c}\sim-0.11\pi$.Läuchli et al. (2006) Hence, in our model (1) for large $|J_{\rm d}|/z$, the phase boundary between the spin-nematic and C-type 120∘-AFM phases is also expected to move from the mean-field line $J_{4}=-2(J_{\times}-1)$ (corresponding to $\theta_{\rm c}\sim-0.35\pi$ at $J_{\rm d}/z\to-\infty$) toward the line $J_{4}=-0.61(J_{\times}-1)$ (corresponding to $\theta_{\rm c}\sim-0.11\pi$), in the direction to enlarge the spin nematic phase.

Third, we must be careful when applying the results for the three-sublattice case to the kagome lattice system. For the spin-1 bilinear-biquadratic model in the kagome lattice, it was found that the 120∘-AFM phase is not present in the ground-state phase diagram and, instead, the trimerized valence-bond-crystal phase appears.Liu et al. (2015) Therefore, in our result, the C-type 120∘-AFM phase should be replaced by a phase corresponding to the trimerized valence-bond-crystal phase, and the boundary lines of the phase as well as the nature of the phase transitions to the phase may also be different from the mean-field results. We expect that our results for the FM, A-type AFM, and spin nematic phases remain valid semi-quantitatively for the kagome-lattice case, since the ground states of these ordered phases are discribed rather well by the mean-field approximation.

V Many-variable Variational Monte Carlo calculation

In this section, we focus on the case of $J_{4}=0$ and large negative $J_{\rm d}$. We numerically explore the ground-state phase diagram using mVMC method, and establish the emergence of the spin nematic phase with ferro-quadrupolar order suggested from the second-order perturbation calculation.

V.1 Method

To analyze the Hamiltonian defined in Eq. (1), we use the mVMC method,Misawa et al. (2019); mVM which can include the spatial correlations and quantum fluctuations beyond the direct product of dimer states. The variational wave function in mVMC is defined as

[TABLE]

Here, $\mathcal{P}_{\rm G}$ is the Gutzwiller factors defined as

[TABLE]

where $n_{l,j,\sigma}$ is the number operator of electron at the $l$th site in the $j$th dimer with spin $\sigma$. By taking $g\to\infty$, we completely exclude the doubly occupied states and express the localized spin-1/2 systems at half filling. The pair-product part $|\phi_{\rm pair}\rangle$ is the generalized pairing wave function defined as

[TABLE]

where $F_{l,j,\sigma,l^{\prime},j^{\prime},\sigma^{\prime}}$ denotes the variational parameters, $c_{l,j,\sigma}^{\dagger}$ represents the creation operator of electron at the $l$th site in the $j$th dimer with spin $\sigma$, $N_{\rm s}=2N$ is the number of spins (electrons) in the system, and $|0\rangle$ is the vacuum of electrons. In this form, we can express spin nematic states by using spin-triplet pairing wave functions.Shindou and Momoi (2009); Shindou et al. (2011) In our calculations, we have imposed $2\times 2$ ($3\times 3$) sublattice structure in the pair-product part for the square lattice (triangular lattice) to express the C-type AFM (C-type 120∘-AFM) state. All the variational parameters are simultaneously optimized by using the stochastic reconfiguration method.Sorella (2001); Sorella et al. (2007)

We note that, if one takes $F_{l,j,\sigma,l^{\prime},j^{\prime},\sigma^{\prime}}$ as

[TABLE]

the wave function becomes

[TABLE]

This result shows that the pair-product state $|\phi_{\rm pair}\rangle$ includes the dimer-product state defined in Eq. (12) as a special case. Although the entanglement between dimers is completely ignored in the dimer-product state, the mVMC can include such entanglement.

V.2 Square Lattice

Here, by using mVMC, we examine stability of the spin nematic phase around $J_{\parallel}+J_{\times}=0$ in the square lattice. As a typical case, we take $J_{\rm d}/z=-2~{}(J_{\rm d}=-8)$, $J_{\parallel}=-1$, and $J_{4}=0$.

To make the initial states of the mVMC calculations, we first impose the external fields that induce the candidate states of the ground state. In this calculation, we consider the FM, C-type AFM, and spin nematic states and take the external fields defined as

[TABLE]

where ${\bm{\pi}}=(\pi,\pi)$ is the wave vector of the C-type AFM ordering. We typically take the amplitude of the external field as unity, i.e., $h_{\rm FM}=1$ for example. We first optimize the variational parameters using the stochastic reconfiguration method under the external fields. Then, by turning off the external fields, we again optimize the variational parameters and obtain the FM, C-type AFM, and spin nematic states. We have also checked that the A-type AFM state, whose initial state can be prepared by imposing the external field $H_{\rm ex,AAFM}=-h_{\rm AAFM}\sum_{j}N_{j}^{z}$, is unstable and becomes one of the other states after the optimization without the external field. We hence omit the result of the A-type AFM state in the following.

Using the optimized wave functions of the FM, C-type AFM, and spin nematic states, we compute the energies, the local total-spin moment $T$ defined by

[TABLE]

and the local spin-quadrupolar moment $Q$ defined by

[TABLE]

for each state. The calculation was performed for finite systems with $N=L\times L$ sites under the periodic boundary condition. We found that the finite-size effects are small for the FM and spin-nematic states, so that we could achieve convergence to thermodynamic-limit values with the data for the systems with up to $L=10$. For the C-type AFM state, however, the system-size dependences of $T$ and $Q$ are large. We therefore performed the calculation for the systems with up to $L=14$ and extrapolated the data of $T(L)$ and $Q(L)$ using the least-square fitting with linear functions of $1/L$, such as $T(L)=T(L=\infty)+a/L$. We note that, in our calculation, only $Q^{(2)}$ becomes finite in the collinear magnetic ordered phases (FM phase and C-type AFM phase) and the spin nematic phase.

Figure 8 (a) shows $J_{\times}$ dependence of the energies for the FM, C-type AFM, and spin nematic states, while Fig. 8 (b) shows the local total-spin moment $T$ and the ferro-quadrupolar moment $Q$ in the ground state. The non-zero $T$ in Fig. 8 (b) indicates the appearance of the magnetic ordered phase such as the FM or C-type AFM phase. The spin nematic phase is characterized by finite spin-quadrupolar moment ($Q>0$) and absence of the magnetic order ($T=0$).

As shown in Fig. 8 (a), for $J_{\times}\leq 0.88$, we find that the FM state is the ground state. Its local moment is still fully polarized ($T=1$) even when we seriously take into account the interdimer correlations.

For $0.89\leq J_{\times}<1.00$, we find that the spin nematic state becomes the ground state even when $J_{4}=0$. In the spin nematic state, we confirm that no spontaneous polarization occurs in the spin degrees of freedom ($T=0$) and the spin-quadrupolar moment $Q$ becomes finite as shown in Fig. 8 (b). This result shows that effects of the interdimer correlations included in mVMC actually stabilize the spin nematic phase. The spin nematic phase widely extends for $J_{\times}<1$ while it does not for $J_{\times}>1$. This is consistent with the result of the second-order perturbation theory, which indicates the stability of the spin nematic phase for $J_{\parallel}+J_{\times}<0$.

At $J_{\times}=1.00$, the energy of the spin nematic state is slightly smaller than that of the C-type AFM state within the system size treated, but they are almost equal, suggesting that the transition point between the phases is in the range $1.00\leq J_{\times}<1.01$. We note that the C-type AFM state is not stable for $J_{\times}<1$, i.e., even if we choose the C-type AFM state as an initial state, the final state after optimization becomes a spin nematic state, which has no spin order. We therefore conclude that the spin nematic phase exists at least for $J_{\times}<1$.

For $J_{\times}\geq 1.01$, the C-type AFM state is the ground state. In contrast to the FM state, the C-type AFM state is largely affected by the interdimer correlations. Due to the quantum fluctuations, the energy of the C-type AFM state obtained by the mVMC method is significantly lower than that of the direct product of the dimer states. The local spin moment $T$ in the mVMC result is also reduced from the saturated value.

From the results above, we conclude that the system with $J_{4}=0$ and sufficiently large negative $J_{\rm d}$ exhibits the spin nematic phase in addition to the FM and C-type AFM phases. Both transitions between the FM and spin-nematic phases and between the spin-nematic and C-type AFM phases are of first order accompanied by a jump of the magnetic moment. These transition properties are the same as the results of the mean-field approximation for those transitions occuring at $J_{4}<0$.

V.3 Triangular Lattice

For the triangular lattice, we perform basically the same calculations as the case of the square lattice. As a typical case, we take $J_{\rm d}/z=-2~{}(J_{\rm d}=-12)$, $J_{\parallel}=-1$, and $J_{4}=0$.

In the triangular lattice, it is expected that the C-type 120*∘-AFM state becomes the ground state in addition to the FM and spin nematic states. To prepare the initial state of the C-type 120∘*-AFM state, we impose the external field defined as

[TABLE]

where $\phi({\bm{r}}_{j})=\phi(x_{j},y_{j})=2\pi x_{j}/3+2\pi y_{j}/3$. We have checked by the mVMC that the A-type AFM state is unstable for $J_{\rm d}/z=-2$.

We performed the calculation for finite systems with $N=6\times 6,12\times 6$, and $12\times 12$ sites under the periodic boundary condition. As in the case of the square lattice, the system-size dependences of the data are small for the FM and spin-nematic states, so that we could obtain a good convergence in the data for the systems with up to $N=12\times 12$ sites. For the C-type 120∘-AFM state, sizable system-size dependences still remain in the results of $T$ and $Q$. We hence extrapolated the data of $T(N)$ and $Q(N)$ by using least-square fitting to linear functions of $1/\sqrt{N}$, e.g., $T(N)=T(N=\infty)+a^{\prime}/\sqrt{N}$. We note that in our calculation, only $Q^{(2)}$ is finite in the FM and spin-nematic phases while $Q^{(1)}$ and $Q^{(2)}$ become finite in the C-type 120*∘*-AFM phase with the coplanar magnetic order.

In Fig. 9(a), we show $J_{\times}$ dependence of the energies for the FM, C-type 120∘-AFM, and spin nematic states, while we show in Fig. 9(b) the local total-spin moment $T$ and the spin-quadrupolar moment $Q$ in the ground state.

We find that the FM state is the ground state for small $J_{\times}$ ($J_{\times}\leq 0.9$). With increasing $J_{\times}$, the system undergoes a first-order transition into the spin nematic phase, accompanied with a level cross of the ground states, around $J_{\times}\sim 0.9$. At the other side of the spin nematic phase, it has turned out that the C-type 120∘-AFM state is unstable for $J_{\times}\leq 1.1$. We can thus conclude safely the appearance of the spin nematic phase for $0.9\lesssim J_{\times}\lesssim 1.1$. In contrast to the square lattice, the spin nematic phase emerges for both $J_{\times}>1$ and $J_{\times}<1$. This is consistent with the result of the perturbation theory, which indicates that the spin nematic phase extends to both sides of $J_{\parallel}+J_{\times}=0$.

At $J_{\times}=1.1$, the energies of the spin nematic and C-type 120∘-AFM states are almost degenerate, and their slopes also seem to be equal. Furthermore, the local spin moment $T$ in the C-type 120∘-AFM state for $J_{\times}>1.1$ decreases with decreasing $J_{\times}$ and seems to gradually vanish at $J_{\times}\simeq 1.1$. These results indicate that the continuous phase transition occurs between the C-type 120∘-AFM and spin nematic phases, which is consistent with the conclusion obtained by the mean-field approximation.

VI Vicinity of the SU(4) symmetric point

In this section, we present an exact analysis on the SU(4) symmetric point of the model (1) with ferromagnetic coupling. We show emergence of various phases including the spin nematic phase and the vector chiral (p-type nematic) phase, by adding perturbations to this symmetric point. All of these phases are transformed to each other through SU(4) rotation. We further argue that the spin nematic state can be stabilized, out of the ground-state manifold at the SU(4) point, by adding a set of Ising interactions as a perturbation.

VI.1 Degenerate ground states of the SU(4) symmetric model

The SU(4) symmetric model with ferromagnetic coupling is given by $J_{4}=4J_{\parallel}<0$ and $J_{\rm d}=J_{\times}=0$, as mentioned in Sec. II. The Hamiltonian reads

[TABLE]

where we set $J_{\rm su4}=-J_{\parallel}=-J_{4}/4>0$. In this paper we call Eq. (36) ferromagnetic SU(4) model. We note that the overall coupling constant is negative opposed to the SU(4) spin-orbital model for the two-orbital Hubbard model at quarter filling,Kugel’ and Khomskii ; Li et al. (1998) whose quantum criticality has been extensively investigated in one dimension.Yamashita et al. (1998); Azaria et al. (1999, 2000); Momoi et al. (2003); Lecheminant and Totsuka (2005, 2006) We rigorously show in the following that the model (36) has various degenerate ground states which are transformed to each other through the SU(4) rotation. One of them is a spin nematic state with ferro-quadrupolar order.

To analyze the SU(4) model, we use the following fifteen local generators of SU(4) group on each dimer. We first adopt the spin dipole operators $T_{j}^{\alpha}$ and $N_{j}^{\alpha}$ ($\alpha=x,y,z$), and the quadrupolar operators $Q_{j}^{(n)}$ ($n=1,\cdots,5$). We further use the vector chiral operators (or p-type nematic operatorsAndreev and Grishchuk )

[TABLE]

($\alpha=x,y,z$), where $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita tensor, and the Heisenberg exchange operator

[TABLE]

All of these fifteen operators, $T_{j}^{\alpha}$, $N_{j}^{\alpha}$, $\chi_{j}^{\alpha}$ ($\alpha=x,y,z$), $Q_{j}^{(n)}$ ($n=1,\cdots,5$), and ${\cal O}_{j}$, are the generators of SU(4) group on dimers.Lecheminant and Totsuka (2006) This is a useful representation of SU(4) generators for our analysis of emerging phases. We note that ${\cal O}_{j}$ is also the generator of the spin-chirality dual transformation on dimers.Momoi et al. (2003) See also Appendix B for a convenient definition of SU(4) generators.

Let us start our argument from the fact that the SU(4) Hamiltonian can be expressed in terms of permutation operators,

[TABLE]

From this, it immediately follows that eigenenergies of the ferromagnetic SU(4) model are lower-bounded by $-3J_{\rm su4}N_{\rm b}/4$, where $N_{\rm b}\equiv zN/2$ is the number of nearest-neighbor dimer pairs.

A trivial ground state with the energy $-3J_{\rm su4}N_{\rm b}/4$ is the state in which all spins are down,

[TABLE]

Due to SU(4) symmetry, any state obtained through global SU(4) rotation to Eq. (40) also belongs to the ground states. We summarize below five typical states of degenerate ground states and their order parameters. Among these states, the antiferromagnetic state and the vector chiral state are transformed to each other through the spin-chirality dual relation,Hikihara et al. (2003); Momoi et al. (2003) whereas the rest of states are invariant, i.e., self dual. We note that, since all of these states have product wave functions, the mean-field approximation presented in Sec. IV gives exact ground states in the ferromagnetic SU(4) model.

VI.1.1 Ferromagnetic state

In the ferromagnetic state $|{\rm FM}\rangle$, all spins are ferromagnetically ordered. The order parameter is given by $\sum_{j}{\bm{T}}_{j}$. In this state, spin SU(2) and time-reversal symmetries are broken.

VI.1.2 Antiferromagnetic state

Applying the SU(4) rotation $U=\prod_{j}(-i)\exp(i\pi T_{j}^{x}/2)\exp(i\pi N_{j}^{x}/2)$ to $|{\rm FM}\rangle$, we obtain the A-type antiferromagnetic state

[TABLE]

The spins have an antiferromagnetic order detected with $\sum_{j}{\bm{N}}_{j}$. In this state, spin SU(2) and time-reversal symmetries are broken.

VI.1.3 Spin nematic state

Applying the rotation $U=\prod_{j}\exp(i\pi Q_{j}^{(3)}/4)$ to $|{\rm FM}\rangle$, we obtain the spin nematic state with ferro-quadrupolar order,

[TABLE]

This state does not have any spin order

[TABLE]

while has a ferro-quadrupolar order

[TABLE]

Thus only spin SU(2) symmetry is broken.

VI.1.4 Vector chiral state

Applying the duality transformationMomoi et al. (2003) $U=\prod_{j}\exp[i\pi(-\sqrt{6}{\cal O}_{j}+1)/8]$ to $|\mbox{A-AF}\rangle$, we obtain the vector chiral (p-type nematic) state

[TABLE]

Each dimer state is a linear combination of the spin singlet state and a spin triplet state with the complex coefficients. The vector chiral state does not have any spin order

[TABLE]

while it has a vector chiral order

[TABLE]

Thus SU(2) and reflection symmetries are broken. This vector chiral order is also accompanied with a quadrupolar order, which corresponds to a quadrupolar moment of the vector chirality.

VI.1.5 Dimer singlet state

The dimer singlet state

[TABLE]

is also obtained, for example, through the SU(4) transformation $\prod_{j}e^{-i3\pi/4}\exp(i\pi N_{j}^{z}/4)|{\rm VC}\rangle$. No symmetry is broken in this state.

VI.2 Perturbations to lift the degeneracy

A weak perturbation to the ferromagnetic SU(4) model can stabilize one of these degenerate ground states, making its phase realize in a finite parameter region. For example, the ferromagnetic phase is realized in the vicinity of the ferromagnetic SU(4) model by adding a weak ferromagnetic Heisenberg interaction on interdimer bonds,

[TABLE]

with $\lambda>0$. The fully polarized state (40) still belongs to the lowest-energy eigenstates of this perturbed Hamiltonian. Another example of exact results is the stability of the dimer singlet phase. One can show that the dimer singlet state (48) is the lowest-energy eigenstate of the perturbed Hamiltonian

[TABLE]

with $\lambda^{\prime}>0$. The other states are, however, not eigenstates usually, if we add SU(2) symmetric perturbations. In the case of the spin nematic state (42), we infer that the spin nematic phase with the ferro-quadrupolar order is stabilized by adding a perturbation of biquadratic interactions,

[TABLE]

with $\lambda^{\prime\prime}>0$, though the state (42) is not an exact eigenstate of this perturbed Hamiltonian any more.

In the rest of this section, we consider an Ising anisotropic perturbation and prove that adding some Ising couplings to the ferromagnetic SU(4) model can stabilize the spin nematic phase. We use the Hamiltonian of the ferromagnetic SU(4) model with additional Ising couplings

[TABLE]

where

[TABLE]

with $J^{z}_{\rm d}>0$. This Hamiltonian conserves the number of dimers having the state $|\sigma_{1}\sigma_{2}\rangle_{j}$, which we denote by $N_{\sigma_{1}\sigma_{2}}$, and hence is block-diagonalized into subspaces characterized by the quantum numbers $\{N_{\sigma_{1}\sigma_{2}}\}=\{N_{\uparrow\uparrow},N_{\downarrow\downarrow},N_{\uparrow\downarrow},N_{\downarrow\uparrow}\}$. In the following, we consider the case that the number of dimers, $N=\sum_{\sigma_{1}\sigma_{2}}N_{\sigma_{1}\sigma_{2}}$, is even.

Our argument has some analogies to the derivation of the so-called $\eta$-pairing superconductivity in an extended Hubbard model.Yang (1989); Essler et al. (1992, 1993) We introduce $\eta$ operators of spins

[TABLE]

The operators $\eta^{\alpha}_{j}$ satisfy the commutation relation

[TABLE]

and form an SU(2) group. The two-dimensional fundamental representation of this SU(2) is spanned with two states $|\uparrow\uparrow\rangle_{j}$ and $|\downarrow\downarrow\rangle_{j}$. These two states carry pseudo-spin-1/2 degrees of freedomHikihara et al. (2017)

[TABLE]

and $\eta^{+}_{j}\equiv\eta^{x}_{j}+i\eta^{y}_{j}$ ($\eta^{-}_{j}\equiv\eta^{x}_{j}-i\eta^{y}_{j}$) is the raising (lowering) operator between them, while the rest of states, $|\uparrow\downarrow\rangle_{j}$ and $|\downarrow\uparrow\rangle_{j}$, are doubly degenerate singlet states,

[TABLE]

Since $\eta_{j}^{\alpha}$ ($\alpha=x,y,z$) are parts of SU(4) generators, $\sum_{j}\eta^{\alpha}_{j}$ commutes with $\mathcal{H}_{\rm su4}$. For later use, we rewrite $\mathcal{H}_{\rm su4}$ and $\mathcal{H}_{\rm int}^{z}$ with SU(4) generators and $\eta$ operators introduced in Eqs. (54),

[TABLE]

The degenerate ground states of the ferromagnetic SU(4) model come from all subspaces of distinct $\{N_{\sigma_{1}\sigma_{2}}\}$. We first consider the case that only the intradimer Ising coupling term $\mathcal{H}^{z}_{\rm d}$ is added to the SU(4) model. Since $\mathcal{H}_{\rm su4}$ is block diagonalized for each subspace of $\{N_{\sigma_{1}\sigma_{2}}\}$ and the eigenvalue of $\mathcal{H}^{z}_{\rm d}$ is $-(J^{z}_{\rm d}/4)(N_{\uparrow\uparrow}+N_{\downarrow\downarrow}-N_{\uparrow\downarrow}-N_{\downarrow\uparrow})$, this Ising term partially lifts the ground-state degeneracy of the SU(4) model. Hence, only the states in the subspace with $N_{\uparrow\downarrow}=N_{\downarrow\uparrow}=0$ have the lowest energy among the degenerate ground states of the SU(4) model. All other states in the subspaces with $N_{\uparrow\downarrow}+N_{\downarrow\uparrow}>0$ acquire finite energy costs of $J^{z}_{\rm d}(N_{\uparrow\downarrow}+N_{\downarrow\uparrow})/2$ compared to the lowest-energy states with $N_{\uparrow\downarrow}=N_{\downarrow\uparrow}=0$ and are gapped out of the ground-state manifold. We note that the role of $\mathcal{H}^{z}_{\rm d}$ in our argument is only to wipe out the states with $N_{\uparrow\downarrow}+N_{\downarrow\uparrow}>0$ from the ground states. Hence, $J^{z}_{\rm d}$ does not have to be small and can be comparable to or larger than $J_{\rm su4}$.

To see the nature of the ground states of the Hamiltonian $\mathcal{H}_{\rm su4}+\mathcal{H}^{z}_{\rm d}$, we consider the projection ${\cal P}$ to the subspace with $N_{\uparrow\downarrow}=N_{\downarrow\uparrow}=0$, to which the ground states belong. In this projected space, the terms containing ${\bm{\eta}}$ spins in the Hamiltonian give non-trivial operators and the rest of terms give constants. As a result, the Hamiltonian $\mathcal{H}_{\rm su4}+\mathcal{H}^{z}_{\rm d}$ in this subspace reduces to a “ferromagnetic” Heisenberg model of ${\bm{\eta}}$ spins,

[TABLE]

except for a constant. Thus the ground states of $\mathcal{H}_{\rm su4}+\mathcal{H}^{z}_{\rm d}$ are the perfectly ferromagnetic states of ${\bm{\eta}}$ spins.

We next add an interdimer Ising coupling term $\mathcal{H}^{z}_{\rm int}$. Note that $\mathcal{H}^{z}_{\rm int}$ also preserves the quantum numbers $\{N_{\sigma_{1}\sigma_{2}}\}$. In a subspace with $N_{\uparrow\downarrow}+N_{\downarrow\uparrow}>0$, the change of the lowest energy induced by $\mathcal{H}^{z}_{\rm int}$ (compared to that in the subspace with $N_{\uparrow\downarrow}=N_{\downarrow\uparrow}=0$) is, at most, of order $\mathcal{O}(N_{\uparrow\downarrow}+N_{\downarrow\uparrow})$. Therefore, if the intradimer Ising term $\mathcal{H}^{z}_{\rm d}$ is sufficiently strong compared to the interdimer Ising term $\mathcal{H}^{z}_{\rm int}$, i.e., $|J^{z}_{\parallel}|,|J^{z}_{\times}|\ll J^{z}_{\rm d}$, the ground state of total Hamiltonian $\mathcal{H}^{\prime}$ [in Eq. (52)] still belongs to the subspace with $N_{\uparrow\downarrow}=N_{\downarrow\uparrow}=0$.

To consider the ground state in the thermodynamic limit, a careful treatment of low-lying excited states is needed here. With increasing the system size, quasi-degenerate low-lying states (also known as Anderson tower states) in finite-size systems whose excitation energies decay in the form of ${\cal O}(1/N)$ also join to the lowest-energy state, forming a symmetry broken ground state in the thermodynamic limit.Anderson (1952); Bernu et al. (1992); Koma and Tasaki (1994) In our model, if the coupling $J_{\rm d}^{z}$ is sufficiently strong so that the energy gap from the ground states to the lowest excited states with $N_{\uparrow\downarrow}+N_{\downarrow\uparrow}>0$ is of order unity, these excitations do not contribute to the formation of ground states in the thermodynamic limit and hence the ground states are properly reproduced by the Hamiltonian in the projected space.

In the projected space with $N_{\uparrow\downarrow}=N_{\downarrow\uparrow}=0$, only $\eta$ operators remain non-trivial and the effective total Hamiltonian is written as

[TABLE]

with

[TABLE]

We thus find that the ground state of $\mathcal{H}^{\prime}$ is described with ${\bm{\eta}}$ spins and the effective Hamiltonian is equivalent to the ferromagnetic XXZ model. When $|\Delta_{\eta}|<1$, i.e., $0<J^{z}_{\parallel}+J^{z}_{\times}<2J_{\rm su4}$, the anisotropy is of easy-plane type, and in the ground state all ${\bm{\eta}}$ spins point to the same direction in the $xy$ plane of the $\eta$ spin space. This ground state is indeed a spin nematic state: ${\bm{\eta}}$ spins have an order $\langle{\bm{\eta}}_{j}\rangle=(\rho\cos\vartheta,\rho\sin\vartheta,0)$ with finite positive constant $\rho$, which corresponds to the ferro-quadrupolar order,

[TABLE]

and $\langle S^{z}_{1,j}+S^{z}_{2,j}\rangle=0$. This order is accompanied with the spontaneous breaking of the global spin-rotation symmetry around the spin $z$ axis. Furthermore, since the ground state belongs to the subspace with $N_{\uparrow\downarrow}=N_{\downarrow\uparrow}=0$, the expectation values of $S^{x}_{l,j}$, $S^{y}_{l,j}$ ($l=1,2$), and $N^{z}_{j}=S^{z}_{1,j}-S^{z}_{2,j}$ are always zeros and hence the ground state does not have any spin (dipole) order,

[TABLE]

In this argument, the ferromagnetic Ising couplings $J^{z}_{\rm d}$ on dimers are important and the choice of the other couplings $J^{z}_{\parallel}$ and $J^{z}_{\times}$ is relatively free inside the region $0<J^{z}_{\parallel}+J^{z}_{\times}<2J_{\rm su4}$ (and $|J^{z}_{\parallel}|,|J^{z}_{\times}|\ll J^{z}_{\rm d}$).

To summarize, we have shown that the ferromagnetic SU(4) model [Eq. (36)] has the degenerate various ground states including the spin nematic state with ferro-quadrupolar order and the vector chiral state. Each phase extends to a finite parameter region of a generalized model around the SU(4) symmetric point. Furthermore, adding some appropriate Ising couplings to the SU(4) model can also stabilize the spin nematic phase out of the degenerate ground states of the SU(4) model.

VII Concluding remarks

In this paper, we have studied the frustrated spin-1/2 dimer model (1) composed of ferromagnetic dimers on two-dimensional lattices with a bilayer structure, to explore the spin nematic phase. We have used various approaches — perturbative calculations, a mean-field approximation, mVMC method, and exact arguments. We have succeeded to show the appearance of the spin nematic phase with ferroquadrupolar order in a wide parameter range, which includes the model with only two-spin interactions and also the SU(4) symmetric model. All the spin nematic states found in this paper are adiabatically connected to each other on each lattice, forming the single spin nematic phase.

The appearance of the spin nematic phase in our model can be understood from two mechanisms, which are active in two distinct regimes in the phase. One is attributed to the effective biquadratic interactions between spin-triplet states in dimers, which couple the spin-quadrupolar degrees of freedom on neighboring dimers. This mechanism is elucidated by the perturbative analysis in the strong ferromagnetic dimer-coupling regime. In this limit, our spin-1/2 model is mapped to the effective spin-1 bilinear-biquadratic model. The effective biquadratic interaction derived perturbatively comes from the four-spin exchange interaction in the original model in the first-order process and also from the two-spin interactions in the second-order process. When this effective biquadratic interaction dominates the effective bilinear interaction, the spin nematic phase emerges.

This mechanism opens a possibility to realize the spin nematic phase in the model with only two-spin interactions. When there is no four-spin interaction in the original model, the effective biquadratic interaction is usually not strong as it comes from the second-order perturbation. However, when the interdimer two-spin interactions in the original model have a strong competition between ferromagnetic and antiferromagnetic couplings, the effective bilinear interaction becomes very weak because of the cancellation in perturbations and hence the effective biquadratic interaction becomes the strongest one. The spin nematic phase thus appears in between the ferromagnetic and antiferromagnetic phases, even if the original Hamiltonian does not contain any four-spin interaction. We note that a similar mechanism to surpress the bilinear interaction was argued for the multi-orbital Hubbard model.Mila and Zhang (2000); Tanaka et al. (2018) To confirm the emergence of the spin nematic phase in the model with only two-spin interactions, we have performed large-size numerical mVMC calculations for our model (1) with only two-spin interactions ($J_{4}=0$) and moderately strong ferromagnetic $J_{\rm d}$. The resultant phase diagrams for the square and triangular lattices indeed exhibit the ferro-quadrupolar spin-nematic phases in finite parameter ranges in between the ferromagnetic and antiferromagnetic phases.

The other mechanism to realize the spin nematic phase is found in the ferromagnetic SU(4) symmetric model. Applying exact arguments, we have shown that a spin-nematic product state becomes one of the exact ground states on the ferromagnetic SU(4) model, because the spin nematic order parameters are given by the generators of SU(4) symmetry. This result suggests that the mean-field approximation with product-state ansatz becomes exact on the ferromagnetic SU(4) model. Due to the SU(4) symmetry, various phases are generated in the vicinity of this high symmetric point. In addition to the spin nematic phase, the vector chiral (p-type nematic) phase as well as the conventional ferromagnetic and antiferromagnetic phases and the dimer singlet phase appear in the vicinity of the SU(4) symmetric point if appropriate perturbative interactions are added. Using an argument similar to the one on the $\eta$-pairing superconductivity, we have also proven that our model at the SU(4) point with appropriate Ising couplings can exhibit the spin nematic phase.

The spin nematic state found in this paper is stable in a wide parameter region of the SU(2) symmetric model without any magnetic field. This is because both the mechanisms for spin nematic ordering we showed are valid at zero field. This is a clear contrast to the spin nematic state found in spin-1/2 frustrated ferromagnets,Shannon et al. (2006); Sindzingre et al. (2009); Hikihara et al. (2008); Sudan et al. (2009); Ueda and Momoi (2013) which is caused by the two-magnon instabilityShannon et al. (2006); Chubukov (1991); Kecke et al. (2007) at the saturation field and hence usually more stable in a strong external magnetic field. Furthermore, the spin nematic phase in our model remains to exist even without any four-spin interaction in contrast to the spin nematic state in the spin-1 bilinear-biquadratic model, which requires a strong biquadratic coupling.

Our analysis also revealed that the phase transition between the ferroquadrupolar spin nematic and antiferromagnetic phases is continuous in many cases. For example, on the triangular lattice, the transition between the spin nematic phase and the C-type antiferromagnetic phase with $120^{\circ}$ structure is always continuous, as far as we studied. This is consistent with the former analysis in the $J_{\rm d}\rightarrow-\infty$ limit.Läuchli et al. (2006)

Though our model is a toy model, there are a few candidate materials which might capture some features of our model. We can find possible candidates in organic magnets, which realize spin systems with high flexibility in the control of exchange interactions.Yamaguchi et al. (2018a, b) In particular, several organic biradical molecules are known to contain dimer structure of two $S=1/2$ spins coupled ferromagnetically.Shiomi et al. (1993); Hosokoshi et al. (1999); Iwase et al. (2013) Arranging these dimers in a two-dimensional lattice with a bilayer structure may provide a playground for searching for the spin nematic state. Another candidate is an SU(4)-symmetric system. It was recently proposed that a Coulomb-impurity lattice on a graphene substrate can realize an SU(4)-symmetric spin-orbital model with a tunable coupling constant.Dou et al. (2016)

Even after a material exhibiting the spin nematic phase is prepared, the direct experimental detection of the phase is not easy, but we can manage to find it by combining various experiments. One of the most striking features of the spin nematic phase is the absence of the spin order down to zero temperature, which we can observe by neutron scattering and nuclear magnetic resonance measurements. In addition, spin-wave analyses revealed that the ferroquadrupolar state in the spin-1 bilinear-biquadratic models has gapless excitations with a linear dispersion.Papanicolaou (1988); Läuchli et al. (2006) These low-energy excitations result in algebraic temperature dependence in thermodynamic properties, e.g., $T^{2}$ dependence in specific heat data in two dimensions. Those properties are also expected to appear in the spin nematic phase in our spin-1/2 dimer models. As for the detection of the spin-quadrupolar order, although it is a challenging task, a few theories were recently proposed for inelastic light scattering,Michaud et al. (2011) resonant inelastic x-ray scattering,Savary and Senthil (2015) and electron spin resonance.Furuya and Momoi (2018) We hope that our study stimulates further the search for the spin nematic phase in real materials.

Acknowledgements.

We thank to Yuko Hosokoshi, Chisa Hotta, Hosho Katsura, Kouichi Okunishi, Nic Shannon, and Hironori Yamaguchi for fruitful discussions. This work was partially supported by JSPS KAKENHI Grant Numbers JP15K05198, JP16H06345, JP16K17746, JP17H02931, and JP16K05425. This work was also supported by Building of Consortia for the Development of Human Resources in Science and Technology from the MEXT of Japan. Our calculation was partly carried out at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo.

Appendix A Numerical results of the mean-field approximation

In this Appendix, we present details of the numerical procedure in the analysis of the mean-field approximation discussed in Sec. IV as well as the results of phase transitions in the approximation.

A.1 Details of the calculation

The mean-field solution was obtained by optimizing numerically the complex coefficients $\{a_{\Lambda,\sigma_{1}\sigma_{2}}\}$ in the product state $|\Phi_{\rm DP}\rangle$ defined in Eqs. (12) and (13) so that the state has the minimum expectation value of the bond Hamiltonian $\mathcal{H}_{jj^{\prime}}$ [Eq. (14)]; The minimized function is $\langle\Phi_{\rm DP}|\mathcal{H}_{jj^{\prime}}|\Phi_{\rm DP}\rangle$ ($j\in{\rm A},j^{\prime}\in{\rm B}$) for the two-sublattice case and $\langle\Phi_{\rm DP}|(\mathcal{H}_{jj^{\prime}}+\mathcal{H}_{j^{\prime}j^{\prime\prime}}+\mathcal{H}_{j^{\prime\prime}j})|\Phi_{\rm DP}\rangle$ ($j\in{\rm A},j^{\prime}\in{\rm B},j^{\prime\prime}\in{\rm C}$) for the three-sublattice case. In the calculation, the normalization condition $\sum_{\sigma_{1}\sigma_{2}}|a_{\Lambda,\sigma_{1}\sigma_{2}}|^{2}=1$ for each sublattice $\Lambda$ was imposed. From the arbitrariness of the global phase factor of $|\Phi_{\rm DP}\rangle$, we further imposed the constraint that $a_{\Lambda,\downarrow\downarrow}$ is real for each $\Lambda$, without loss of generality.

The minimization was achieved by using the steepest descent method. Since the optimization process often becomes slow and is trapped in a local minimum in the steepest descent method, we performed 1000 calculations starting from randomly-prepared initial states for each parameter point. The calculations were continued until the optimization converged or the method reached 10000 iterations. We note that we achieved the convergence for all the 1000 runs for most of the parameter points treated and for 165 runs even at the worst case. We then adopted the state giving the lowest energy as the ground state.

For determining the phase diagrams, we performed the calculation on $J_{\times}$ versus $J_{4}$ planes with several values of $J_{\rm d}/z$, varying $J_{\times}$ and $J_{4}$ with intervals of $0.1$. (Note that $J_{\parallel}$ is fixed to be $J_{\parallel}=-1$.) Figures 4 and 6 show typical examples of the obtained phase diagrams. In order to explore the nature of the phase transitions, we performed the calculation on several parameter lines with fixed $J_{\rm d}/z$ and $J_{4}$ ($J_{\times}$), varying $J_{\times}$ ($J_{4}$) with intervals of $0.01$ or $0.001$.

A.2 Transition between spin-nematic and A-type AFM phases

The phase transition between the spin nematic phase and the A-type AFM phase appears in both two- and three-sublattice cases. We find that this transition is continuous for all parameter lines studied. Figure 10 shows the $J_{\times}$-dependence of the Néel-spin moment around the transition in the parameter lines $J_{4}=-6.0$ and $J_{\rm d}/z=0,-0.01,-0.10$ in the two-sublattice case. For $J_{\rm d}/z\lesssim-0.10$, the Néel-spin moment rises from zero continuously with a moderate slope. The slope is steeper as $|J_{\rm d}|/z$ is smaller, however, the transition is still continuous even for $J_{\rm d}/z=-0.01$. At $J_{\rm d}=0$, the spin nematic phase reaches to the boundary region between the spin-nematic and dimer-singlet phases with degenerate ground states, and the order parameters exhibit finite jumps when the system moves from this region into the A-type AFM phase. We note that the steep rise of the Néel-spin moment at $J_{\rm d}/z\to-0$ was observed in both two- and three-sublattice cases. For not too small $|J_{\rm d}|/z$, $|{\bm{N}}^{\rm MF}_{\Lambda}|^{2}$ rises linearly with a moderate slope, suggesting $|{\bm{N}}^{\rm MF}_{\Lambda}|\propto\sqrt{J_{\times}-J_{\times,{\rm c1}}}$, where $J_{\times,{\rm c1}}$ is the critical value.

A.3 Transition between spin-nematic and C-type 120∘-AFM phases

The phase transition between the spin nematic phase and the C-type 120∘-AFM phase occurs in the three-sublattice case for large $|J_{\rm d}|/z$. This transition is found to be continuous. Figure 11 presents the data of the squared total-spin moment in a dimer, $|{\bm{T}}^{\rm MF}_{\rm A}|^{2}$, for $J_{\rm d}/z=-3.5$ and $J_{4}=-0.1,-0.5,-1.5$. When $|J_{4}|$ is not too small, the total-spin moment rises from zero continuously with a finite slope, indicating ${\bm{T}}^{\rm MF}_{\Lambda}\propto\sqrt{J_{\times}-J_{\times,{\rm c2}}}$, where $J_{\times,{\rm c2}}$ is the critical value. The slope becomes steeper as $J_{4}$ approaches zero. At $J_{4}=0$, the spin nematic phase vanishes (within the mean-field approximation) and there occurs a direct transition between the FM and C-type 120∘-AFM phases via a special point with the degenerate ground states at $J_{\times}=-J_{\parallel}=1$. (See Appendix C.3.)

A.4 Transition between A-type AFM and C-type 120∘-AFM phases

The phase transition between the A-type AFM phase and the C-type 120∘-AFM phase occurs in the three-sublattice case for large $|J_{\rm d}|/z$. Figure 12 shows the $J_{4}$-dependence of the total-spin and Néel-spin moments on the parameter line with $J_{\rm d}/z=-3.0$ and $J_{\times}=4.0$. The order parameters exhibit a clear jump at the transition. The transition thus turns out to be the first-order one, occuring between two magnetically-ordered phases with distinct symmetries.

Appendix B SU(4) transformation on dimers

In this Appendix, we briefly describe SU(4) transformation on a dimer. The fifteen generators of SU(4) groupLecheminant and Totsuka (2006) on a dimer are given by the spin operators $T_{j}^{\alpha}$ and $N_{j}^{\alpha}$, the vector chiral operators $\chi_{j}^{\alpha}$ ($\alpha=x,y,z$), the quadrupolar operators $Q_{j}^{(n)}$ ($n=1,\cdots,5$), and the spin exchange operator ${\cal O}_{j}$. Hereafter we omit the dimer index $j$. Using these operators, we conveniently define the fifteen generators $\lambda_{n}$ ($n=1,\cdots,15$) of SU(4) group as follows:

[TABLE]

As for the orthonormal bases on a dimer, we use the following states

[TABLE]

Using these definitions, one can explicitly show that the $4\times 4$ matrix representation of the generators $\lambda_{n}$ on the basis vector $(|x\rangle,|y\rangle,|z\rangle,|0\rangle)$ coincides with the $4\times 4$ generalized Gell-Man matricesGreiner and Müller (1994) which generate SU(4) algebra,

[TABLE]

[TABLE]

Arbitrary SU(4) transformation is given by these generators parameterized by generalized Euler angles.Tilma et al. (2002)

Appendix C Non-trivial degeneracies in the mean-field approximation

We describe non-trivial degeneracies in the mean-field solutions in the two types of phase boundaries and a triple point found in Sec. IV. One boundary exists between the ferromagnetic and spin-nematic phases, which indicates emergent SU(3) symmetry. The other exists in the boundary between the spin-nematic and dimer-singlet phases.

C.1 Emergent SU(3) symmetry in the boundary between the FM and spin-nematic

phases

In the phase boundary between the FM phase and the spin nematic phase with ferro-quadrupolar order, the ground state manifold has non-trivial degeneracy corresponding to SU(3) rotation. We first briefly summarize SU(3) rotation on a dimer. The eight generators $\lambda_{n}$ ($n=1,\cdots,8$) of SU(3) rotation are given by $T^{\alpha}$ ($\alpha=x,y,z$) and $Q^{(n)}$ ($n=1,\cdots,5$) as in Eqs. (B). Using these generators, we can write arbitrary SU(3) rotation in the formByrd (1997)

[TABLE]

with parameters ($\alpha,\beta,\gamma,\theta,a,b,c,\phi$). This gives arbitrary unitary transformation among the spin triplet. As an initial state, we use the state $|z\rangle_{j}$. Since this state is invariant under the last (right) four rotations, the transformed state is simply written as

[TABLE]

In the mean-field approximation performed in Sec. IV, the FM state and the spin nematic state are degenerate in energy in the boundary between these two phases. We further find that, in the mean-field solutions, any state obtained by arbitrary global SU(3) rotation to the FM state, which is a translationally invariant product state of Eq. (136), also takes the exactly same energy. This is emergent non-trivial degeneracy associated with the global SU(3) rotation. As mentioned in Sec. II.3, the Hamiltonian has the exact global SU(3) symmetry and the corresponding degenerate ground states only in the parameter space (5), which exists inside of the mean-field phase boundary. This degeneracy thus remains in the whole phase boundary in the mean-field approximation even though the model Hamiltonian does not possess the SU(3) symmetry.

C.2 Boundary between the spin-nematic and dimer-singlet phases

In the mean-field approximation, the spin nematic phase touches with the dimer singlet phase in a finite parameter plane at $J_{\rm d}=0$, as shown in Figs. 4(a) and 6(a). In this phase boundary, the ground state manifold has non-trivial degeneracy corresponding to ${\rm SU(2)}\times{\rm SU(2)}$ rotation, even though the Hamiltonian does not have this symmetry except for the special case of $J_{\rm d}=J_{\times}=0$. Among the spin-nematic ground states, we consider the product state of $|z\rangle$ without loss of generality. All other states are related with the global SU(2) rotation. We next apply SU(2) rotation only to one spin of each dimer,

[TABLE]

where $\hat{{\bm{\alpha}}}$ ($\alpha=x,y,z$) denotes the unit vector parallel to $\alpha$-axis and $\hat{{\bm{\omega}}}$ the unit vector parallel to ${\bm{\omega}}$. The quadrupolar state ($|{\bm{\omega}}|=0$) is continuously transformed to the vector chiral states ($0<|{\bm{\omega}}|<\pi$) and the dimer singlet state ($|{\bm{\omega}}|=\pi$) when $\hat{\bm{\omega}}\parallel\hat{\bm{z}}$. By a straightforward calculation, one can show that the translationally invariant product state of these dimer bases also has the same energy as the ferro-quadrupolar state. Thus the ground state manifold has the same degrees of freedom as global ${\rm SU(2)}\times{\rm SU(2)}$ rotation. This degeneracy exists in the whole phase boundary between the spin-nematic and dimer-singlet phases in the mean-field approximation.

C.3 Triple point for the FM, spin nematic, and C-type AFM phases

For large $|J_{\rm d}|/z$, the triple point for the FM, spin nematic, and C-type (120∘-)AFM phases exists at $J_{\times}=-J_{\parallel}$ and $J_{4}=0$ in both two-sublattice and three-sublattice cases in the mean-field approximation. On this triple point, any product state in which each dimer independently take an arbitrary superposition of spin-triplet states has the same energy to form the massively-degenerate ground state manifold in the mean-field level.