Excitonic instability towards a Potts-nematic quantum paramagnet

Magnetic frustration can lead to peculiar magnetic orderings that break a discrete symmetry of the lattice in addition to the fundamental magnetic symmetries (i.e., spin rotation invariance and time-reversal symmetry). In this work, we focus on frustrated quantum magnets and study the nature of the quantum phase transition between a paramagnet and a magnetically ordered state with broken threefold ($Z_3$) crystal rotation symmetry. We predict the transition to happen in two stages, giving rise to an intermediate nematic phase in which rotation symmetry is broken but the system remains magnetically disordered. The nematic transition is described by the three-state Potts model. This prediction is based on an analysis of bound states formed from two-magnon excitations in the paramagnet, which become gapless while single-magnon excitations remain gapped. By considering three different lattice models, we demonstrate a generic instability towards two-magnon bound state formation in the Potts-nematic nematic channel. We present both numerical results and a general analytical perturbative formula for the bound state binding energy similar to BCS theory. We further discuss a number of different materials that realize key features of the model considered, and thus provide promising venues for possible experimental observation.


I. INTRODUCTION
The study of ordering phenomena is one of the central pillars of condensed matter physics.Landau's theory of phase transitions provides a general framework for understanding the structure of ordered phases and the nature of the transitions between different phases.The order parameter is a fundamental concept in the study of phase transitions and characterizes the ordered phase in terms of its symmetry.In the simplest and most common cases, the ordering transition is associated with the breaking of a fundamental symmetry, such as translation symmetry in the case of crystalline order; charge conservation in the case of superconductivity; and time-reversal combined with spin-rotation symmetry in the case of magnetic order.Ordering transitions with a much richer structure arise when, in addition to the fundamental symmetry, one or more secondary symmetries are broken in the ordered state, such as spatial rotation or inversion symmetry.This occurs, for instance, in unconventional non-s-wave superconductors and superfluids, where the internal structure of the Cooper pairs (i.e., their spin and orbital angular momentum) can lead to spontaneous spatial anisotropy associated with broken rotation symmetry [1,2].Similar phenomena can occur in magnets when the exchange interactions are (strongly) frustrated, which prohibits the formation of simple ferromagnetic or antiferromagnetic Néel order.Whereas the ferromagnet and antiferromagnet do not break the symmetries of the crystal lattice, the more complicated magnetic configurations that arise as a result of frustration generally do break one or more spatial symmetries.
In this paper, we study a particular class of such frustrated magnets: magnets which spontaneously break a discrete Z 3 (i.e., threefold) crystal rotation symmetry in the ordered state.Such phases naturally arise in simple Heisenberg spin models with frustrated exchange interactions, for instance in two-dimensional (or layered quasi-two-dimensional) magnets with triangular [3][4][5] or honeycomb lattice [6] geometry.The phase diagrams of triangular and honeycomb lattice spin models with further neighbor exchange couplings are known to include collinear single-Q magnetic states, characterized by staggered ferromagnetic stripes (see Fig. 1), which are stabilized by the "order by disorder" mechanism [3,[7][8][9][10].Since a single ordering wave vector is spontaneously selected from a set of three wave vectors related by crystal rotation symmetry, the latter symmetry is spontaneously broken in these single-Q states.This may be compared to Heisenberg antiferromagnets with tetragonal symmetry, such as frustrated square lattice magnets, which can support single-Q stripe order described by either one of two wave vectors (π, 0) and (0, π) [9,11].Such single-Q phases break a Z 2 Ising symmetry corresponding to the two possible orientations of the stripes, which in two dimensions implies a finite temperature Ising transition to an Ising nematic phase with no magnetic order [11].In three dimensions the Ising and magnetic transitions may still occur at different temperatures, in which case a so-called vestigial Ising nematic phase arises [12].The emergence of Ising nematicity is well-known example of the rich ordering phenomenology that can occur in complex magnets with multiple broken symmetries [13].
In contrast to the Ising case, the breaking of a discrete Z 3 symmetry is governed by the three-state Potts model [14] and thus provides a manifestly differentand much less studied-window into the nature of magnetic ordering in unconventional complex magnets [15][16][17].The purpose of this paper is to explore the nature of the phase transition in magnets that break an additional Z 3 symmetry, such as a threefold crystal rotation, and thus belong to the three-state Potts universality class.As far as the thermal phase transition is concerned, general arguments suggest that, similar to the Ising case, a finite temperate transition to a Potts-nematic phase is possible [16,17].
Rather than the thermal phase transition, our goal here is to study the nature of the zero temperature quantum phase transition between a quantum paramagnet and a magnetically ordered phase which breaks a discrete Z 3 crystal rotation symmetry.The question we specifically seek to answer is whether an intermediate Potts-nematic phase, in which the system remains paramagnetic but breaks the Z 3 rotation symmetry, can exist.
To address this question, we consider a minimal XXZ model of S = 1 quantum spins with single-ion anisotropy [18].In this minimal model, defined in Eq. ( 8), the single-ion anisotropy D controls the transition between the quantum paramagnet realized for sufficiently large D and the magnetically ordered state.The latter is controlled by the frustrated exchange interactions, which are chosen such that, at the classical level, the energy is minimized by a set of degenerate ordering wave vectors related by threefold rotation symmetry.Following the approach of Ref. 18, we start from the quantum paramagnet and study its instability towards the formation of two-magnon bound states as D is lowered, which may occur due to an attractive interaction between singlemagnon excitations.If the energy gap of two-magnon bound states closes at a value of D larger than the value at which the single-magnon gap closes, this indicates an instability towards bound state formation.Furthermore, if the internal structure of the bound state solution, defined by the relative angular momentum of the magnons, breaks the Z 3 rotation symmetry, this suggests that proliferation of two-magnon excitations gives rise to an intermediate Potts-nematic quantum paramagnetic phase.Such analysis is similar in spirit to the classic Cooper problem in the context of superconductivity [19], which addresses the instability of the Fermi sea to the formation of two-particle bound states (Cooper pairs).
We apply this analysis to three different models of frustrated quantum magnets: the triangular lattice, the honeycomb lattice, and the three-dimensional face-centered cubic (FCC) lattice.In all three cases we focus on a part of the phase diagram where rotation symmetry broken single-Q magnetic order is favored.We obtain a matrix Schrödinger equation for the two-magnon bound states, show that it decouples in channels of distinct symmetry quantum numbers, and solve it numerically.We demonstrate that when the exchange interactions favor magnetic order broken rotation symmetry, generically twomagnon bound states have largest (positive) binding energy in the nematic channel, suggesting a generic instability towards a Potts-nematic phase in this class of quantum magnets.The numerical analysis is supplemented + F q 0 5 J 5 s 5 h j 9 w P n 8 A i 5 m N U w = = < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " K W h 0 R L J 0 b w 8 e m / x 3 P U 2 c j X 7 5 q z e I a S q Z s l Q Q Y 7 q + l 9 g g I 9 p y K t i 0 1 E s N S w g d k y H r O q q I Z C b I 5 g d P 8 Z l T B j i K t S t l 8 V z 9 P Z E R a c x E h q 5 T E j s y y 9 5 M / M / r p j a 6 by a perturbative analytical approach, is independent of the underlying lattice structure and in two dimensions yields an expression for the width of the intermediate Potts-nematic phase similar to BCS theory.In particular, as expressed in Eq. ( 27), the binding energy, which is proportional to the width of the nematic phase, depends on the single-magnon density of states and the effective coupling constant in the nematic channel.
To make a connection with possible candidate materials, we identified a number of compounds that can host the proposed Potts-nematic vestigial phase under certain conditions.We point to two specific material classes, which we propose as a prospective platform to test our predictions.One of the two groups is the family of transition metal trichalcogenides (MP X 3 ), which have attracted much attention recently [20][21][22][23][24][25][26].The other is iron-intercalated transition metal dichalcogenide Fe x NbS 2 , which was recently reported to host a Z 3breaking magnetism, revealed in optical measurements [27,28].
The paper is organized as follows.In the next Section II we define our model and give general introduction to the physics of Z 3 -symmetry breaking magnets.We proceed with the Section III, in which we develop the main formalism used in our work and derive its perturbative limit.We then apply this machinery to the specific cases of triangular, honeycomb, and face-centered cubic lattices and present the results in Section IV.Finally, we propose concrete experimental implementations in Section V and conclude with a discussion (Section VI).The details of the derivation for a non-Bravais lattice are relegated to Appendix A. Throughout the paper we consider zero temperature and work in units of the lattice constant (i.e., a = 1).

II. ROTATION SYMMETRY BREAKING IN FRUSTRATED ANTIFERROMAGNETS
This work is focused on frustrated Heisenberg antiferromagnets described by the general Hamiltonian where J ij = J ji are the exchange coupling constants.
When the spins {S i } are considered as classical vectors, the exchange couplings describe interactions between pairs of spins which either favor alignment (J ij < 0) or anti-alignment (J ij > 0).Frustration arises when the interactions between different pairs of spins compete and the energetic requirement of perfect alignment or antialignment cannot be simultaneously satisfied for all pairs of spins.Frustration may originate from the geometry of the lattice, as is the case for the nearest neighbor triangular lattice antiferromagnet, or from the presence of multiple competing exchange constants, such as in J 1 -J 2 -J 3 models.
A standard approach to determining the classical spin configuration with the lowest energy is the Luttinger-Tisza method [29,30].It relies on expanding the spin variables in Fourier modes as S rα = q e iq•(r+ α) S qα , where α denotes the position of sublattice site α with respect to the lattice vector r, and writing the spin Hamiltonian in terms of the Fourier modes as H = (N/2) q (J q ) αβ S * qα • S qβ .Here J q is the Fourier transform of the exchange couplings, which in the case of non-Bravais lattices is matrix in sublattice space.Minimization of the classical energy then requires computing the eigenvalues of J q for all q and finding the minimum of lowest eigenvalue branch with respect to q.The minimum occurs at some specific wave vector Q or at a set of symmetry related (i.e., "degenerate") wave vectors {Q i }.The latter case is a consequence and signature of frustration.The eigenmodes corresponding to the minima of J q are generally referred to as the optimal Luttinger-Tisza modes and can be used to construct spin states that minimize the magnetic energy.This construction must be subject to the fixed-length constraint, i.e., |S i | 2 = S 2 , however, and this implies that, in the case of non-Bravais lattices, it may not be possible to construct ground states from optimal Luttinger-Tisza modes only.In this case, the Luttinger-Tisza method only provides a lower bound the classical energy.Note that in the case of primitive Bravais lattices, however, it is always possible to construct ground states from the optimal Luttinger-Tisza modes.
In this work we focus on a particular class of frustrated magnetic systems for which the minima of the eigenval-ues of J q occur at a set of three symmetry-related wave vectors {Q i } (with i = 1, 2, 3) in some part of the phase diagram.The degenerate wave vectors are related by crystal rotation symmetry and satisfy the two additive relations Note that the first of these relations implies that the Fourier modes are real.Prominent examples of Heisenberg magnets which exhibit such phases in the phase diagram are the triangular J 1 -J 2 model [3,10], the honeycomb lattice J 1 -J 2 -J 3 model [6], and the threedimensional face-centered cubic (FCC) lattice J 1 -J 2 model [8,31,32].These three models will be the subject of the remainder of this paper and will be discussed in more detail below.In the case of the triangular and honeycomb lattices, the wave vectors Q i correspond to the high symmetry M points of the Brillouin zone shown in Fig. 1(a), which are clearly related by threefold rotation.
The triangular lattice provides the simplest example of the magnetism of interest.The minima of J q are located at the three M points when J 1 and J 2 are both positive (and hence frustrated) and 1/8 < J 2 /J 1 < 1.
To understand the nature of the classical ground state manifold when the optimal Luttinger-Tisza modes derive from the wave vectors {Q i }, consider first the simple case of a Bravais lattice.The triangular and FCC lattices are examples.Up to a global spin rotation, the most general parametrization of the Fourier modes S Qi consistent with the fixed length requirement takes the form where S is the spin length and m = (m x , m y , m z ) ≡ (sin θ cos φ, sin θ sin φ, cos θ) is a unit vector defined in terms of the angles (θ, φ).The angles (θ, φ) parametrize a family of degenerate classical ground states, which can be represented as points on the Bloch sphere.Generically, such points correspond to triple-Q magnetic states, but there are special lines and points on the Bloch sphere correspond to double-Q and single-Q states.In particular, states which correspond to m = (1, 0, 0) and its equivalents define single-Q configurations.
A natural and convenient way to resolve the degeneracy of the classical ground states is to consider composite order parameters built from the primary magnetic Fourier modes.The distinct symmetry properties of the composite orders unambiguously expose the symmetries of magnetic configurations and can thus be used to sharply distinguish classical ground states.In the present context, we define the nematic composite order parameter n and the chiral order parameter χ as Here ω = e 2πi/3 is a cubic root of unity.It is clear from these definitions that both n and χ transform trivially under translations, but n is time-reversal even whereas χ is time-reversal odd.Most importantly, since nonzero n must originate from an unequal contribution of the three wave vectors, it signals the breaking of crystal rotation symmetry in the magnetic state.Substituting the Fourier mode representation (3) into ( 4) and ( 5) one finds n and chi in terms of m as This shows that n is zero when m 2 x = m 2 y = m 2 z , and that finite χ requires all components of m to be nonzero (i.e., noncoplanar order).Note in particular that single-Q ordered states are associated with n = 0.
Since any energy functional of n contains a cubic Pottsanisotropy ∝ n 3 + (n * ) 3 , thus giving rise to three distinct nematic ground states, n describes Pott-nematic order [].As a result, n is uniquely associated with the spontaneous breaking of threefold crystal rotation symmetry in the magnetically ordered phase.
The continuous degeneracy of the classical ground states parametrized by m does not survive the effect of fluctuations.This is due to the well known "orderby-disorder" mechanism, which tends to favor collinear spin alignment [3,[7][8][9][10][11].For instance, whereas the classical energy is independent of the angles (θ, φ), the energy associated with the quantum zero-point motion of the magnon excitations does depend on the orientation of m, and is lowest for the single-Q states such as m = (1, 0, 0) [3,5,10,31].

III. INSTABILITIES OF A QUANTUM PARAMAGNET
In this section we turn to the central part of this work.As outlined in the introduction, our goal is to determine the instability of a quantum paramagnet towards the formation of nematic magnon bound states.To achieve this, we start from an XXZ model for a system of S = 1 spins, introduced in Ref. 18, which is defined by the Hamiltonian The Hamiltonian is a sum of two terms.The first term is of the general form introduced in (1) and describes exchange couplings between the spins.Specifically, pairs of spins separated by a distance r − r = δ interact with coupling constant J δ .Note that compared to (1) here we have also introduced an XXZ anisotropy given by η δ ≡ J z δ /J δ , which denotes the ratio of the exchange couplings in the z direction and the xy plane.This is useful, since the parameter η δ can be interpreted as controlling the interactions between magnons.As far as the magnetic ground state is concerned, there is no difference with (1).The second term in ( 8) describes a single-ion anisotropy of strength D > 0.
For simplicity, in what follows we focus the discussion on the case of primitive Bravais lattices and leave the straightforward generalization to non-Bravais lattices to Appendix A, the results of which will be used when applying the analysis to the honeycomb lattice.Since we are interested in the quantum behavior of this model we consider the system at T = 0.
As argued in Ref. 18, in the limit where D is much larger than the exchange couplings J δ , the ground state is a quantum paramagnet with all spins in the |0 ≡ |S = 1, 0 state.We denote the paramagnetic ground state as |Ω = ⊗ j |0 j .Instead, when D → 0 one expects a magnetically ordered state, the nature of which is determined by the exchange couplings J δ .Crucially, in what follows we assume that the exchange couplings take values such that J q has minima at three degenerate wave vectors {Q i }, giving rise to a degenerate magnetic ground manifold at the classical level.This is precisely the part of the phase diagram considered in the previous section.In this case, D describes the transition from a quantum paramagnet to a magnetically ordered state which spontaneously breaks a discrete Z 3 rotation symmetry.Our goal is to examine this transition and establish whether or not the paramagnetic ground state |Ω is unstable towards the formation of nonmagnetic two-magnon bound states in the vicinity of the magnetic transition.

A. Schrödinger equation for magnon pairs
To address this question, we follow the approach of Ref. 18 and consider a general two-magnon state |ψ defined as where ψ rr is the wave function, which is a function of the positions of the two magnons.By projecting Schrödinger equation H |ψ = E |ψ into the two subspace of twomagnon states, we find the equation where J z δ ≡ η δ J δ .Note that, in accordance with (9), the two-magnon wavefunction is understood to vanish when r = r (i.e., ψ rr = 0) [33].The first key observation is that the second term on the right hand side, which is proportional to J z δ , generically describes an attractive interaction between the magnons, whereas the first term represents a kinetic term for the two magnons.To proceed, we expand the wave function ψ rr in Fourier modes as ψ rr = K,q e iK•(r+r )/2 e iq•(r−r ) ψ q (K), (11) where K is the momentum conjugate to the center of mass coordinate and q is conjugate to the relative coordinate r − r .Substituting the Fourier mode expansion into the Schrödinger equation ( 10) yields the equation where can be viewed as a kinetic energy of the two magnons, and B δ , defined as corresponds to the real-space wavefunction as a function of the relative coordinate r − r = δ.Note that in the case of zero center of mass momentum one has E q (K = 0) ≡ E q = 2ε q , where ε q = δ J δ cos q • δ is the singlemagnon dispersion.To solve Eq. ( 12), we divide by the kernel E − 2D − E q (K) and use the definition of B δ to obtain a matrix equation given by where the matrix M δδ is defined as Solutions to the Schrödinger equation are then determined by the condition In our case, we seek solutions corresponding to zero energy (E = 0), which means that the energy for creating two-magnon excitations vanishes.The general strategy for solving (17) after setting E = 0 is to diagonalize M and determine when the eigenvalues of M − 1 vanish.This is done as a function of D and the first eigenvalue which vanishes as D decreases corresponds to the solution of interest.Inspection of the structure of M reveals that it is possible to bring it in block diagonal form by making use of the symmetry properties of the underlying lattice and the center of mass momentum K. To illustrate how group theory machinery can be employed to block diagonalize M, it is useful to consider the example of the triangular lattice J 1 -J 2 model discussed in the previous section.
The triangular lattice has six nearest neighbor vectors a j and six next-nearest neighbor vectors a j [34], which we label j = 0, . . ., 5. As a result, δ, δ in ( 15) and ( 16) take values in this set of twelve vectors.Let us furthermore consider the case K = 0.The matrix M is then block diagonalized by a unitary matrix U of the form where the matrix elements of U are Here U transforms from the basis of (next)-nearest neighbor vectors to a basis of angular momentum labeled by l = −2, . . ., 3. Since the triangular lattice (as well as the center of mass momentum K = Γ) has full hexagonal symmetry, different angular momentum channels cannot mix, implying that U † MU is block diagonal.The blocks therefore correspond to distinct angular momentum channels, which we denote M l and are given by Here f l q and g l q are the nearest neighbor and next-nearest neighbor symmetry-adapted lattice harmonics of the triangular lattice and are defined as It is worth noting that the matrix in (20) can be decomposed as which leads to further simplifications within a weak coupling approach, see Sec.III B. We further point out that the l = 0 angular momentum channel is special and should be excluded, since it leads to a two-magnon wavefunction ψ rr that is nonzero when r = r, thus contradicting the assumptions that underlie (10).Equation ( 17) now reduces to a set of separate equations in each angular momentum channel given by Det(M l − 1) = 0.These equations can be solvedgenerally numerically-to obtain the critical value D * The angular momentum channel with the largest value D * c then determines the structure of excitons, much like the structure of Cooper pairs is determined by the relative angular momentum of the constituent electron pairs.
It is worth noting that symmetries generally imply degeneracies of angular momentum channels.For instance, in the case of the triangular lattice, l = ±2 (as well as l = ±1) are degenerate and form a single irreducible channel.
This detailed example of the triangular lattice serves to illustrate the general symmetry-based method for solving (17).In particular, the matrix M can always be block diagonalized such that each block M i corresponds to distinct symmetry quantum numbers (i.e., representation of the lattice symmetry group) and matrix elements within each block are products of symmetry-adapted lattice harmonics.So, Eq. ( 17) reduces to a set of decoupled equations in each symmetry channel of the form Det(M i − 1) = 0. Their solutions determine whether two-magnon bound states can form and if so, in what symmetry channel they form.This general method applies to finite center-of-mass momentum K as well, with the modification that the relevant symmetry group is the little group of K.

B. Weak coupling
Further analysis of Eq. ( 20) is possible when the interaction between the magnons is small.It follows from Eq. ( 10) that the interaction between magnons originates from the S z S z coupling and is parametrized by η δ = J z δ /J δ , which may be treated as a perturbative parameter.Similar to Ref. 18, this assumption will allow us to derive a compact expression for the exciton binding energy, which is shown to depend on: (i) an effective coupling constant and (ii) magnon density of states.
In Section III A, we used (11) to express the twomagnon problem to Schrödinger equation in terms of center-of-mass and relative coordinates.Under this transformation, the magnon interaction translates into quantum well potential for the relative coordinate r − r = δ.As a result, the weak coupling limit is equivalent to the shallow well approximation in quantum mechanics.This is why we expect the bound state binding energy ε b to be small compared to the exchange couplings J δ , while the size of the two-magnon bound state (i.e., the "exciton") should be large on the lattice scale.Consequently, the extent of the exciton in reciprocal space will be small, which allows us to expand the dispersion (13) in the vicinity of its minima.
To illustrate the weak coupling analysis, we will again consider the case of triangular lattice, and relegate the generalization to non-Bravais lattices to Appendix A. Following the scheme outlined above, we assume that the summand in (20) is concentrated in the vicinities of the wave vectors Q i , the locations of the minima of (13).Assuming further that K = 0, we expand E q in the de-nominator up to the second order in p = q − Q i and evaluate the lattice harmonics in the numerator at Q i .We also use the fact that the summand is invariant under rotation of q by 2π/3, which allows to consider just one out of three Q i points; the contributions from each wave vector must be equal.Then the matrix structure of ( 20) is determined by a single q point and has rank one.Evidently, for fixed angular momentum l it has just one nonzero eigenvalue corresponding to the eigenvector (f l Qi , g l Qi ) T .That leads to the following equation, where the summation was replaced by integration in the vicinity of Q i : Here V BZ is the area of the Brillouin zone, p is a small momentum reckoned from Q i , and The effective coupling constants λ l and the exciton binding energy ε b are defined as Integral ( 24) is most easily performed by switching to the integration over single-magnon energy ξ = p T A i p/2: where ω c is a cutoff on the order of the exchange energies J 1,2 and ν Qi is magnon density of states at one of the three minima.The latter is finite in two dimensions at a quadratic dispersion minimum and can be evaluated by A i diagonalization followed by a linear coordinate transform for the integration variables p x , p y .If ρ 1 and ρ 2 are the eigenvalues of A i (obviously, they are the same for each Q i ), then ν Qi ∝ 1/ √ ρ 1 ρ 2 .Note that flattening of the minima enhances the density of states and increases the binding energy.The value of the latter follows from ( 26) which is a key result of this section.The expression for the binding energy closely resembles the one for the Cooper problem of two attracting electrons above the two-dimensional Fermi surface.That is expected as both problems are in a weak coupling regime and effectively two spacial dimensions.
The energy of the two-magnon state is given by the familiar expression for two-particle bound states, i.e., E = 2(D + ε Qi ) − ε b , which implies that zero energy solutions exist for a critical value (23).Since the binding energy is positive for channels with nonzero (attractive) coupling constant, one has D * c > D c , demonstrating that the quantum paramagnet is unstable towards the formation of non-magnetic two-magnon excitons.The internal structure of the excitons, as defined by their symmetry + F q 0 5 J 5 s 5 h j 9 w P n 8 A i 5 m N U w = = < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " K W h 0 R L J 0 b w 8 e m / x 3 P U 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " o M Y c x 8 n A F v X 5 n n Y P e e d y K O z s j l 4 = " > + F q 0 5 J 5 s 5 h j 9 w P n 8 A j q O N V Q = = < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " H 0 y W J m P a y + a l T f G a V A I e x s i U N n q u / J z I S a T 2 J f N s Z E T P S y 9 5 M / M / r p S a 8 9 j I u k 9 Q w S R e L w l R g E + P Z 3 z j g i l E j J p Y Q q r i 9 FIG. 2. Two-magnon "exciton" bound states in the triangular lattice model.(a) Binding energy ε b defined in Eq. ( 23) as a function of η = J z 1,2 /J1,2 for different values of J2/J1.Solid (and colored) curves correspond to numerical solutions of Eq. ( 20) for l = 2. Dashed (and black) curves correspond to the binding energy computed using Eq. ( 27) (with prefactor fitted), showing excellent agreement for small η.(b) The binding energies as a function of J2/J1 computed numerically from Eq. ( 20) for the l = 2 (solid curves) and l = 3 (dashed curves) angular momentum channels.Black squares indicate where the nematic l = 2 channel becomes dominant.(c) The binding energies as a function of J2/J1 in the interval 1/8 < J2/J1 < 1 (bounded by vertical dashed lines).(d) Single-magnon density of states νM (in units of J −1 1 ) for magnon excitations on top of the paramagnetic phase evaluated at the M points.νM is given by Eq. ( 28) and explains the general behavior of ε b shown in (c).
quantum numbers, directly follows from comparing the coupling constants λ l : the angular momentum channel with the largest coupling gives rise to the largest value of D * c and determines the leading instability.As explained in the previous Section, zero angular momentum channel should be eliminated.
We expect the result (27) to be general for twodimensional materials with three equivalent minima of the magnon dispersion.Details of the model in each case are encoded into magnon density of states and effective coupling constants λ l .For non-Bravais lattices, Eq. ( 25) should be modified by introducing projectors to the lower magnon band as explained in Appendix A. We also comment on generalization to three spacial dimensions in the next Section.

IV. POTTS-NEMATIC MAGNON BOUND STATES
In this section, which forms the heart of the paper, we present a detailed application of the theory developed in the previous section to three particular lattice models that exhibit rotation symmetry broken magnetism.Specifically, we consider the triangular and honeycomb lattice in two dimensions, and the FCC lattice in three dimensions.We present both the outcome of the exact solution of the two-particle problem obtained by numerical integration in Eq. ( 20) and the prediction of the analytical perturbation theory given by (26).For the case of the FCC lattice we show how weak coupling approach, which in two dimensions gives rise to (27), still yields useful approximations to the numerically exact result.

A. Triangular lattice magnets
We begin by considering the triangular lattice model briefly introduced in Sec.II (see also Fig. 1).This model includes antiferromagnetic exchange couplings between nearest neighbors, given by J 1 > 0, and next-nearest neighbors, given by J 2 > 0, which define the XXZ part of (8).Since both couplings are antiferromagnetic the interactions are frustrated.We specifically focus on the part of the phase diagram where 1/8 < J 2 /J 1 < 1, since in this regime J q has minima at three commensurate M -point wave vectors (i.e., the centers of the edges of the Brillouin zone), thus implying a magnetically ordered state of collinear stripes at small D [35] [shown in Fig. 1(b)].Furthermore, since the dispersion of the single-magnon excitations in the quantum paramagnetic state is proportional to J q (i.e., ε q ∼ J q , see Sec.III A), its minima are also located at the M points.
Consider first the numerical solutions of Eq. ( 20).These are obtained straightforwardly by numerical evaluation of the sum and the key results are presented in Fig. 2. Most importantly, we find that in the regime of interest, i.e., 1/8 J 2 /J 1 < 1, the binding energy for two-magnon bound states is positive and largest in the l = 2 angular momentum channel, which is the doubly degenerate (Potts-)nematic channel.In Fig. 2(a) we show the binding energy ε b obtained from solving (20) for l = 2 as a function of η, which controls the strength of the attractive magnon interaction.For convenience here we have taken η 1,2 = η, where η i = J z i /J i , i.e., the ratio between J z and J is equal for nearest and next-nearest neighbor couplings.Different (solid and colored) curves correspond to different values of J 2 /J 1 , showing that the binding energy increases as J 2 /J 1 decreases.For comparison, the dashed black curves show the binding energy  calculated using Eq. ( 27), which was derived under the assumption that η is small.We see that the BCS-like formula ( 27) is in excellent agreement with the numerically exact result in this regime.
In panels (b) and (c) of Fig. 2 we show the binding energy, or equivalently D * c − D c , as a function of the strength of J 2 relative to J 1 for various values of η.As expected, panel (c) shows that the binding energy increases for increasing η.Furthermore, we find that the binding energy is a convex function of J 2 /J 1 and increases towards to boundaries of the interval 1/8 < J 2 /J 1 < 1.As discussed below, this behavior can be explained by (27), in particular the density of states ν M at the M points.
The exchange coupling ratio J 2 /J 1 = 1/8 is of special importance since the minima of J q change at this value, thus implying that the classical ground state changes.In particular, for J 2 /J 1 < 1/8 the minima of J q are located at the K points, giving rise to the well-known coplanar 120 • spiral state on the classical level.Note that in the vicinity of J 2 /J 1 = 1/8 both the M points and the K are local minima of J q ; it is the global minimum which changes at J 2 /J 1 = 1/8.
In the limit J 2 /J 1 = 0, Ref. 18 showed that twomagnon bound states can form in the (nondegenerate) l = 3 channel, which is the channel naturally associated with the K points.Excitons with l = 3 break an Ising Z 2 symmetry, which corresponds to the chirality of the 120 • classical order.To examine the behavior in the vicinity of J 2 /J 1 = 1/8 we compute the binding energy in the interval 0 < J 2 /J 1 < 1/4 for both channels and different values of η, and show the result in panel (b) of Fig. 2. As is clear, in the close vicinity of J 2 /J 1 = 1/8 the binding energy is positive in both channels, but is larger for the "chiral" l = 3 channel.At some value J 2 /J 1 1/8, which is marked by a black squares in Fig. 2(b) and increases with η, the "nematic" l = 2 channel becomes dominant, implying a transition from a chiral instability towards a Potts-nematic instability.Note that the curves exhibit a cusp (i.e., discontinuous derivative) at J 2 /J 1 = 1/8, which is due to the aforementioned abrupt change of the global minimum.
To further understand the formation of nematic twomagnon bound states, we follow the weak coupling approach described in Sec.III B and compute the magnon density of states ν M and effective coupling constant λ by expanding around the M points.The density of states is given by and the effective coupling constant is found as These two parameters enter the expression for the binding energy in Eq. ( 27) and have been used to compute the black dashed curves in Fig. 2(a).(The prefactor coming from the cutoff was fitted.)The density of states, which is shown in panel (d) of Fig. 2, is singular and diverges at J 2 = J 1 .This is caused by the fact that for J 2 > J 1 the minima of J q move away from M towards Γ, making the M points saddle points rather than minima.In contrast, ν M is regular at J 2 /J 1 = 1/8, because M remains a local minimum.Since the effective coupling ( 29) remains finite and does not change significantly in the interval 1/8 < J 2 /J 1 < 1, the main dependence of the binding energy ε b on J 2 /J 1 ratio comes from J 2 /J 1 dependence of ν M , which is illustrated by Figs.2(c) and (d).
As explained in Section III A, excitonic solution in the Potts-nematic channel correspond to l = ±2 values of the angular momentum and thus has a twofold degeneracy, which we observe both in numerical and weak coupling approaches.In the latter, it is straightforward to obtain an explicit expression for the wavefunctions in the real space, which up to a constant factor reads: + F q 0 5 J 5 s 5 h j 9 w P n 8 A i 5 m N U w = = < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " K W h 0 R L J 0 b w 8 e m / x 3 P U 2 + F q 0 5 J 5 s 5 h j 9 w P n 8 A j q O N V Q = = < / l a t e x i t > (c) FIG. 5. (a-1) Zigzag (blue) and striped (yellow) phases on the honeycomb lattice for antiferromagnetic nearest neighbor coupling (J1 < 0 case).Effective Potts-nematic coupling constant (31) is positive only at J2 > J1/2, so vestigial state can emerge only in the part of the zigzag domain (shown in darker blue).Lines demonstrate contours at which perturbative theory predicts the width of the intermediate phase to be (D * c − Dc)/J1 ∼ 10 −1 (solid), ∼ 10 −2 (dashed), and ∼ 10 −3 (dot-dashed).We assumed J z i = Ji when plotting this picture.(a-2) Similar plot for J1 > 0 case.In that case vestigial phase can appear for all points inside the zigzag and striped magnetic domains.(b) (c) The core real-space (ρ) structure of the Potts-nematic exciton on the honeycomb lattice in the perturbative regime for a) zigzag phase; b) striped phase.The spin-up magnon resides on the first sublattice, while the spin-down magnon is located at a site, shifted by ρ.The solution is doubly degenerate in each case.Slow exponential decay is not shown.
Here ρ is the vector connecting the two lattice sites occupied by magnons constituting the exciton and ω = e 2πi/3 is the cubic root of unity.The factor e −r/lε , where 1 appears in the next orders of the perturbation theory and makes the wavefunction normalizable.Eq. ( 30) represents the fact that the solution is a linear combination of exponentials e iMi•ρ with the wavevectors M i pointing to the minima of the single-magnon dispersion, so that the exciton is formed form the most lowlying excitations.The three exponentials are combined in a function with l = ±2 in the same way as the nematic order parameter (4) is constituted from the spin densities at M i points.Taking then real (imaginary) part of l = 2 wavefunction [Eq.(30)] gives the wavefunctions with the symmetry of d xy (d x 2 −y 2 ) orbitals.The third linearly independent combination of exponentials e iM1•ρ + e iM2•ρ + e iM3•ρ provides the rotation-symmetric l = 0 wavefunction, which contradicts assumptions made in (10) and thus should be discarded.
We show the bound state wavefunctions (30) graphically in Fig. 3 ignoring the slow decaying factor.Note that the radial structure is not completely trivial: one can see that the angular harmonics are not aligned on the neighboring shells but are rotated with respect to each other.That makes ψ(ρ) obey the approximate translational symmetry ρ → ρ + 2δ 1(2) with doubled period reflecting the fact the 2M i ∼ 0 in the reciprocal space.(The translation symmetry is not exact due to the decaying factor e −ρ/lε .)Note also that the excitonic wavefunction ψ(ρ) nullifies when the two magnons reside on the third-nearest sites.In accordance with this, we found that introducing third-nearest neighbor interaction J 3 does not change the value of the effective coupling (29).
Finally, to make sure that excitons with zero centerof-mass momentum K represent the leading instability, we run a numerical calculation for the case of a non-zero K, which still boils down to solving (17) with M given by Eq. (16).In this place we opted for fixing the binding energy at the value of ε = 0.05J 1 and evaluating the minimal interaction strength η required to create such bound state.The result is presented in Fig. 4, where we allowed the center-of-mass momentum K to take values across a path in the Brillouin zone, connecting the points Γ − K − M − Γ.Indeed, we observed that K = 0 requires the minimal amount of interaction for the bound state formation, which justifies considering K = 0 in the rest of the paper.Interestingly, we found that K = M i choice is also favorable for exciton formation, however, this instability remains subleading.

B. Honeycomb lattice magnets
Next, we examine Potts-nematic two-magnon bound state formation on the honeycomb lattice.The motivation for considering the honeycomb lattice is twofold.First, it has the same symmetry group as the triangular lattice, but contrary to the latter it has sublattice structure, thus providing a generalization of the analysis presented in Sec.III to non-Bravais lattices.Second, a variety of experimentally available candidate materials are realizations of the honeycomb lattice model which makes the honeycomb lattice highly relevant from an experimental perspective.The connection to experimental compounds will be discussed in more detail below in Sec.V.
As mentioned in Sec.II, as far as the exchange interactions included in Eq. ( 8) are concerned, we consider a frustrated J 1 -J 2 -J 3 honeycomb model with up to third-nearest neighbor couplings.The classical phase diagram of the isotropic J 1 -J 2 -J 3 Heisenberg model is known [6,36] and notably includes two magnetic phases with broken threefold rotation symmetry.The corresponding regions in the phase diagram are defined by property that the minima of J q are located at the M points.Note that since J q is a matrix in the space of the two sublattices, two distinct solutions exist, which are related by the relative orientation of the spins on the sublattices, i.e., either aligned or anti-aligned.The single-Q collinear orderings with M -point wave vectors that are selected by fluctuations are schematically shown in Fig. 5, and are generally referred to as stripe (sublatticeeven) and zigzag (sublattice-odd) orders.In our analysis we restrict to exchange couplings (J 1 , J 2 , J 3 ) which favor zigzag or stripe order on the classical level.
To account for the two sublattices in the system, we extend the general theory presented in Sec.III A to the case of a non-Bravais lattice (see Appendix A).The main occurring difference is the appearance of projectors to the Bloch bands in the expression for the effective matrix elements (A12).We present solutions of the modified equations obtained by numerical Brillouine zone integration complemented by the band summation in Fig. 6.For that plot, we picked values of the spin couplings inside the zigzag domain on the Fig. 5 (a 3) and varied the effective interaction strength η = η i = J z i /J i .We observed that bound state with the angular momentum l = ±2 has the largest binding energy and its dependence on η is captured well by the perturbation theory with a fitted prefactor (dashed line).As in the case of triangular lattice, that means that proliferation of excitations of this kind constitutes the leading instability, when D is lowered and that marks the transition to the Potts-nematic phase.
The weak coupling theory for a non-Bravais lattice in general follows the previous derivation given in Section III B. Again, the main contribution stems from the most low-lying excitations, which now also means that it is sufficient to consider the lower band.We expand around the minima of the latter (located at M i points) and, after some algebra, get expressions for the effective coupling constant in the Potts-nematic channel which appears to be the same for both zigzag and stripe M -orderings.The density of low-lying states at M points equals: ) where J 13 = J 2 1 + J 1 J 3 − 4J 2 3 and β = sign (J 1 − 3J 3 ).Eqs. (31) and (32) should be plugged into the general perturbative formula for the binding energy (27), which gives a result consistent with numerical computation (Fig. 6).
Perturbative wavefunction of l = ±2 excitions are again doubly degenerate and appear to be a linear combination of the exponentials e iMir as in the case of triangular lattice, but bear a sublattice structure, which slightly differs between the zigzag and stripe phase.We present it graphically in Fig. 5 (b) for underlying zigzag ordering an Fig. 5 (c) for the stripe one, choosing the d xy and d x 2 −y 2 basis and using the ρ plane, where ρ is the distance between the magnons in the pair.We assumed that the spin-up magnon resides on the first sublattice; the opposite case leads to similar pictures.The figure depicts the core of the exciton wavefunction with omitted slow decay factor e −ρ/lε .
With obtained analytical approximation, we explore exciton formation for all striped phases in the phase diagram of Fig. 5 (a).The fist important observation is that for the bound state to emerge, the effective coupling constant (31) needs to be positive.It turns out that at J 1 < 0 and inside the zigzag domain (J 3 > 0), λ > 0 only when J 2 /J 1 > 1/2, which defines a subregion potent to develop a vestigial phase (shown by darker tone and solid frame).Contrary to that, for the stripe phase (J 3 < 0) and for both phases at J 1 > 0, the effective coupling is positive everywhere inside the magnetic domain, so that in the latter cases magnetic ordering should always be preempted by a vestigial phase.
We now turn to the discussion of how prominent should be the effect of exciton formation across the honeycomb lattice phase diagram.As suggested by Eq. ( 27) and common knowledge, enhancement of the density of states at minima of the dispersion is beneficial for the formation of the bound state.Inspection of (32) implies that such enhancement happens in the vicinity of the right boundaries of zigzag and stripe domains.That can be explained by the fact that at these lines the system undergoes phase transition into the spiral ordering, characterized by a non-commensurate value of ordering wavevector Q.After the transition, the minima of the coupling matrix J q eigenvalues and, consequently, the minima of the single-magnon dispersion on top of paramagnet drift away from the M-points, turning the M-point itself into a saddle point rather than the minimum.Before the transition but near the boundary the effect is manifested by appearing of the soft direction at each M-point, leading to the increase of the single-magnon density of states.This phenomenon is completely analogous to the one, which happens on the triangular lattice in the vicinity of J 2 /J 1 = 1 point as described in Section IV A. To further illustrate this statement and make it quantitative, we plot contours, at which our perturbative theory predicts the binding energy (and hence the width of the Pottsnematic phase) to be of the order 10 −1 (solid line), 10 −2 (dashed line), and 10 −3 (dot-dashed line) in the units of |J 1 | [Fig.5 (a)].As can be seen in the picture, the most perceptible effect appears near the aforementioned boundary of the domains.

C. FCC lattice magnets
In the final part of this section we turn to an application of our theory to a model in three dimensions: the face-centric cubic (FCC) lattice.Our goal is to demonstrate that the analysis developed in Sec.III is not lim-ited to two dimensions, but extends to three dimensions as well.The three-dimensional case is different in one important way, however.Since it is well known that bound states do not exist for arbitrarily shallow potential in three dimensions, making it a threshold phenomenon, we expect the two-magnon bound states to occur only for finite interaction between the magnons.In particular, this means that Eq. ( 27) does not apply to the threedimensional case.Below we describe how the analysis of Sec.III B can be adapted to the case of the FCC lattice.
We consider a J 1 -J 2 model on the FCC lattice with antiferromagnetic J 1 and small ferromagnetic J 2 .For a simple nearest neighbor FCC antiferromagnet (J 2 = 0), the Fourier transformed exchange coupling J q has minima on lines in momentum space.The lines of minima reflects the frustration of the FFC lattice and are defined by the constraint (q x , q y ) = (2π, 0) with arbitrary q z (all symmetry equivalent constraints also define minima).Introducing ferromagnetic coupling between the next-nearest neighbors (J 2 < 0) lifts that degeneracy and results in minima of J q at the three wave vectors Q 1 = (2π, 0, 0), Q 2 = (0, 2π, 0), and Q 3 = (0, 0, 2π).These wave vectors satisfy Eq. ( 2).The fact that there are three equivalent minima points implies that Z 3 symmetry is broken in the ordered phase with the ordering wavevector provided by one of these minima.Such ordering is formed by stripes, which now can take one of the three directions parallel to the principal axes in the three-dimensional space.We thus expect that in this case also one might expect intermediate Potts-nematic state at a range of single-ion anisotropy parameter D between the ordered and paramagnetic phases.
We start with the demonstration of the solution for the binding energy ε b , obtained by numerical integration in (16).It is shown in Fig. 7, where ε b is plotted against the interaction strength η = η i = J z i /J i for a number of J 2 /J 1 values.As predicted above, we observe that the bound state is formed only at η > η c , which is the main distinction with respect to the two-dimensional case.At η η c , the binding energy has approximately quadratic dependence on η − η c .A more accurate formula (34) contains a logarithmic correction and is discussed below.We plot it in dashed in Fig. 7 and observe that it captures well the numerical dependence at small η −η c (with fitted parameters).
Having identified the threshold nature of the bound state formation, we studied the dependence of the critical interaction η c on the spin-spin couplings (Fig. 8).The resulting dependence is monotonous and drops to zero in a sharp nonanalytical way at J 2 → 0. This behavior is explained by the fact that at J 2 = 0 the singlemagnon dispersion reaches its minimum on a set of lines, which effectively renders the problem two-dimensional.From a quantitative perspective, it is interesting that for isotropic Heisenberg model with η = J z i /J i = 1, the bound states exist only at |J 2 |/J 1 0.1.
We proceed with the weak coupling theory for the three-dimensional case.Naive attempt at taking the in- tegral in (16) using expansion (24) in the vicinity of Q i points leads to a linearly divergent integral due to increased phase volume in three dimensions.The common resolution would be to subtract from the integrand of ( 24) its value at ε b = 0 (to be taken numerically over the whole Brillouin zone) and then proceed with the regularized integral.In the regime of small binding energies ε b that leads to the following implicit expression for ε b (η): where B is a constant determined by the density of states in the vicinity of M -points.The result indeed signifies that the bound state exists only at η > η c .However, we found that Eq. ( 33) is suitable to describe the asymptotics of the exact solution only at very small values of ε b , which correspond to a much smaller scale than the one in Fig. 7.This observation is explained by the following considerations.As discussed in the beginning of this subsection, at |J 2 | → 0 the dispersion acquires a very soft direction near the minimum, which eventually transforms into a line of degenerate minima.This process leads to the increase of the single-magnon density of states, which, as discussed in previous sections, enhances the effect.Given the smallness of relevant J 2 (in Fig. 7 |J 2 |/J 1 0.1) that means that at moderate ε b the value of ( 20) is not captured well by the quadratic expansion (24).Instead, one can expect logarithmical corrections to (33) due to the nascent crossover to 2D.Indeed, we find that a corrected formula captures the numerical dependence much better (see dashed lines in Fig. 7).
To conclude, our consideration predicts that the effect of exciton formation and appearance of the vestigial phase is weaker in three dimensions for nonspecific values of the spin-spin interactions.However, in the vicinity of special points in the phase diagram (J 2 = 0 in our case), the effect becomes more pronounced due to the enhanced density of states of relevant magnons.

V. CONNECTION TO EXPERIMENT
The analysis presented in the previous section indicates that the instability towards a Potts-nematic phase is rather general and does not depend on the specific realization of the proposed model on any specific lattice.This generality motivates the important question whether material realizations of Potts-nematic quantum magnets exist among experimentally known material systems.In pursuit of this question, in this section we highlight and discuss a number of promising candidate materials to which our theory applies or relates.
For our theory to directly apply, candidate materials should satisfy three requirements.First, the interactions between the spin must be such that magnetic order with broken Z 3 crystal symmetry is favored.This implies a tacit constraint on the symmetry of the crystal structure and limits the search for materials to hexagonal, trigonal, tetrahedral, or octahedral systems.In practice, the search should be focused on materials or families of materials in which rotation symmetry broken magnetic order is observed, which is clearest signature of the nature of the exchange couplings.Examples will be discussed below.The second requirement is that the candidate materials exhibit significant easy-plane single-ion anisotropy, which implies that spin-orbit coupling is important.In our model, the single-ion anisotropy controls the (quantum phase) transition between the paramagnet and the ordered state, and therefore would be responsible for driving a magnetically ordered state into a putative Potts-nematic phase.Closely related to this second requirement is the third requirement: the quantum spins must be S = 1 degrees of freedom, such that in the limit of large single-ion anisotropy the ground state is a perfect quantum paramagnet product state.In principle other integer-spin systems may be considered, such as S = 2, since these also admit a quantum paramagnetic ground state.
Given these requirements, one of the more promising material classes is the family of transition-metal phosphorous trichalcogenides (M PX 3 ) [20][21][22][23][24][25][26]37].These are layered van der Waals materials in which the transition metal M =(Fe, Ni, Mn, Co) sites are magnetic and form a honeycomb lattice within each layer.The X are occupied by chalcogen atoms, typically sulfur (S) or selenium (Se).The M PX 3 materials have attracted much attention recently as a versatile platform for exploring intrinsic two-dimensional magnetism in few-layer or monolayer systems.While some members of this class show antiferromagnetic Neél order, others were found to realize antiferromagnetic zigzag order-precisely the type of order discussed in Sec.IV B. Furthermore, strong easyaxis or easy-plane behavior was observed, such as Isinglike easy-axis behavior in FePS 3 [21,23,38] and XY-like easy-plane behavior in NiPS 3 [22,39].To explain the observed magnetic orderings across the family of M PX 3 compounds a J 1 -J 2 -J 3 honeycomb lattice spin model with single-ion anisotropy was proposed [21,22,[40][41][42], which, quite remarkably, is identical to the honeycomb lattice model studied in this work.It is for this reason that the trichalcogenide magnets offer a particularly compelling venue for Potts-nematicity associated with rotation symmetry broken magnetic order.
Within the family of M PX 3 materials NiPS 3 deserves special attention.Not only do the ordered moments of the observed zigzag phase lie in the plane [22,39], suggesting easy-plane anisotropy with D > 0, but the Ni 2+ ions also give rise to S = 1 spins on the honeycomb sites.This implies that NiPS 3 provides a realization of the model studied in this work.
It is important to point out one caveat concerning the bulk M PX 3 materials, in particular concerning the stacking of the constituent layers in the bulk structures.In the sulfur-based compounds M PS 3 with M =(Fe,Ni,Mn) the layers are stacked in monoclinic fashion, giving rise to space group C2/m.This means that the threefold rotation symmetry of the constituent honeycomb layers is already broken in the bulk structure, thus precluding the spontaneous breaking of rotation symmetry by the magnetic interactions.The coupling between the layers is believed to be very weak, however, and may be weak enough to prevent a selection of the nematic director by the monoclinic stacking direction.Indeed, encouraging evidence for this has recently been reported by linear and non-linear optical measurements [26], which have found three different magnetic domains below the Neél temperature within one homogeneous structural domain.
A second material of interest is the iron-intercalated transition metal dichalcogenide (TMD) Fe x NbS 2 with x = 1/3, which is a member of the class of intercalated TMDs M x T A 2 [43][44][45].Here M is a transition metal, T =(Ta, Nb), and A=(S, Se).The magnetic Fe atoms are intercalated between the adjacent van der Waals layers of 2H-NbS 2 and form a triangular lattice with √ 3 × √ 3 periodicity.The two neighboring triangular lattice layers of Fe atoms that make up the unit cell (and are located at c = 1/4 and c = 3/4) are shifted with respect to each other, effectively forming a honeycomb structure.The resulting space group of the intercalated structure is P 6 3 22.Importantly, both optical birefringence [27,28] and detailed neutron crystal diffraction [28] measurements have reported observation of rotational symmetry broken stripe order.Here, the two triangular lattice layers each form the ordering pattern shown in Fig. 1 and together form the honeycomb lattice stripy order shown in Fig. 5.As such, the observed magnetic order falls in the class of orderings considered in our work.(Note that Ref. 28 reported great sensitivity of the ordering pattern to small changes in the intercalation parameter x close to the commensurate value 1/3.)Furthermore, the Fe 2+ ions give rise to S = 2 spins [46,47].To understand the observed magnetic orderings, recent theory works proposed a Heisenberg spin with extended exchange interactions and single-ion anisotropy [48,49], essentially equivalent to our honeycomb lattice model.An open question is the nature of the single-ion anisotropy, however, with initial estimates suggesting it might be negative, thus implying an easy-axis anisotropy [48,49].
We conclude this section with a general remark.To obtain the nematic phase at low temperatures one needs a suitable tuning knob to weaken and eventually destroy the magnetic order.This can be provided by varying strain and pressure in the system.An alternative possibility is the formation of Potts-nematic phase at the thermal phase transition as function of temperature, rather than the quantum phase transition studied here.A theory of thermal transition to the Potts-nematic phase will be reported elsewhere.

VI. DISCUSSION AND CONCLUSION
The central result of this paper is the demonstration that a magnetically disordered Potts-nematic phase can exist as an intermediate phase at the quantum phase transition from a quantum paramagnet to a Z 3 threefold rotation symmetry broken antiferromagnet.Since such a phase is magnetically disordered, spin-rotation symmetry and time-reversal symmetry are preserved, yet a Z 3 crystal rotation symmetry is broken, thus giving rising to nematic order.The relevant order parameter for the nematic phase is expressed in terms of bilinears of the primary magnetic degrees of freedom, i.e., the spin variables.
Our conclusion is inferred from the analysis of twomagnon excitations in the paramagnetic state of an S = 1 quantum spin model, in which the transition between the paramagnet and the magnetic state is regulated by the value of easy-plane single-ion anisotropy.With both numerical and analytical methods we established the existence of two-magnon bound states-"excitons"-with the angular momentum l = 2.That makes the twoparticle gap close slightly before the single-particle gap as the transition is approached, leading to the proliferation of excitons.The latter can be described by a Jastrow ansatz [18,50].We derived a general equation for the excitonic bound state in terms of lattice harmonics [Eqs.(17) and (20)] and derived its analytical solution for the weak-coupling regime (27).The latter expresses the binding energy and, hence, the width of the vestigial phase, as a function of magnon density of states at the dispersion minima (in the paramagnetic state) and the effective interaction strength, which is a linear combination of the spin-spin coupling constants.Analytical perturbation theory is corroborated by a detailed numer-ical analysis for the cases of triangular, honeycomb, and face-centered cubic lattices (see Figs. 2-8).
It is important to describe the relation between the character of the vestigial phase determined by the twomagnon wavefunction and the properties of the singlemagnon dispersion above the paramagnetic phase.In the present paper we assumed that the dispersion has three degenerate minima in the Brillouin zone (Q 1 , Q 2 , and Q 3 momenta, satisfying 2Q i ∼ 0), which are related to each other by an in-plane rotation for the cases of triangular and honeycomb lattice.Under these assumptions, perturbative expression (25) for the effective coupling constants in different angular momenta channels (l) nullifies for l = 1 and l = 3 leaving only Potts-nematic channel with l = 2 and unphysical symmetric wavefunction with l = 0. Similarly, we checked that in the setup of Ref. 18, in which the dispersion had two inequivalent minima, the only physical channel arising in the perturbation theory is the chiral one (l = 3).This relation provides a universal scheme allowing to determine the nature of the vestigial phase given the single-particle spectrum.In the numerical calculation we observed that the effective couplings in the 'wrong' angular momentum channels are nonzero, but they are small and thus irrelevant deep inside the domain with a given magnetic ordering [see Fig. 2(b)].
While in the present work we focused on the specific case of an S = 1 spin model, our analysis straightforwardly generalizes to any integer spin system.Indeed, the ground state at large values of the single-ion anisotropy in that case is still a perfect paramagnet, given by a product of S z = 0 states on each site.A spin flip to a state with S z = s on any site costs an energy Ds 2 and therefore, the energetically lowest excitations in this case still correspond to s = ±1.It means that for the purpose of determining the instability of the paramagnet marked by a gap closure, one can restrict the consideration to these two branches only, which makes all further analysis the same as in the S = 1 scenario.In contrast, we expect that half-integer models behave qualitatively different, since they do not admit a paramagnetic product state.
The prediction of the formation of Potts-nematic phase at low temperatures in the vicinity of the quantum phase transition poses a natural question: does there exist a similar vestigial phase at the thermal phase transition between a striped antiferromagnet and the symmetry hightemperature phase?This question will be addressed in a separate paper.
We conclude by emphasizing that there exist a number of different materials which provide promising venues for experimental observation of Potts-nematicity in quantum magnets.As discussed in Sec.V, given its currently established magnetic properties, the M PX 3 family of transition-metal trichalcogenides is of particular interest.Within this class of materials NiPS 3 is most promising, as it is likely to provide a realization of the model considered in this work.Another possible candidate material is the Fe-intercalated dichalcogenide NbS 2 .

<
l a t e x i t s h a 1 _ b a s e 6 4 = " 2 R R wx L X l Y 8 T R O I o M 9 8 j 2 W c O j p r o = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o M Q L 2 F X g n o M e v E Y 0 T w g W c L s p D c Z M j u 7 z M w K I e Q T v H h Q x K t f 5 M 2 / c Z L s Q R M L G o q q b r q 7 g k R w b V z 3 2 8 m t r W 9 s b u W 3 C z u 7 e / s H x c O j p o 5 T x b D B Y h G r d k A 1 C i 6 x Y b g R 2 E 4 U 0 i g Q 2 A p G t z O / 9 Y R K 8 1 g + m n G C f k Q H k o e c U W O l h z I 9 7 x V L b s W d g 6 w S L y M l y F D v F b + 6 / Z i l E U r D B N W 6 4 7 m J 8 S d U G c 4 E T g v d V G N C 2 Y g O s G O p p B F q f z I / d Ur O r N I n Y a x s S U P m 6 u + J C Y 2 0 H k e B 7 Y y o G e p l b y b

a 0 1 FIG. 1 .
FIG. 1. Triangular lattice model.(a) Triangular lattice Heisenberg model with exchange couplings J1 between nearest neighbors (exemplified by red arrows) and J2 between next-nearest neighbor (blue arrow).When both J1,2 > 0 and 1/8 < J2/J1 < 1, the minima of the Fourier-transformed exchange coupling matrix Jq are located at the three M points in the Brillouine zone with the wave vectors Q1,2,3, indicated by blue dots on the right.The generating lattice vectors a1,2 are indicated in red.(b) Three possible realizations of triangular lattice single-Q ordering corresponding to the three M -point wave vectors.Black (white) solid dots represent classical moments pointing out of (into) the page.

c
at which magnon pairs can form with zero energy.Each angular momentum channel will yield a different value of D * c .The critical values D * c are then compared to D c , the value of D at which the gap to single magnon excitations closes, which is obtained by setting 2D c + E Qi = 2(D c + ε Qi ) = 0.It thus follows that D c is determined by the minimum of the magnon dispersion ε Qi .(Note that ε Qi is negative.)When D * c > D c for one or more l, two-magnon bound states can form with positive binding energy

12 <
l a t e x i t s h a 1 _ b a s e 6 4 = " j d f 3

1 <
3 h z x s 6 L 8 + 5 8 L F o L T j 5 z D H / g f P 4 A N O 6 S Z Q = = < / l a t e x i t > l a t e x i t s h a 1 _ b a s e 6 4 = " 2 R R w x L X l Y 8 T R O I o M 9 8 j 2 W c O j p r o = " > A A A B 6 n i c b

0 1 3 d x 2 y 2 < l a t e x i t s h a 1 _ b a s e 6 4 =
v w 3 b t s c t P X B w O O 9 G W b m e T F n S t v 2 t 1 V Y W 9 / Y 3 C p u l 3 Z 2 9 / Y P y o d H b R U l k t A W i X g k u x 5 W l D N B W 5 p p T r u x p D j 0 O O 1 4 k 9 u Z 3 3 m k U r F I P O g 0 p m 6 I R 4I F j G B t p J 4 / z J 4 G t Y t 0 U J s O y x W 7 a s + B V o m T k w r k a A 7 L X 3 0 / I k l I h S Y c K 9 V z 7 F i 7 G Z a a E U 6 n p X 6 i a I z J B I 9 o z 1 C B Q 6 r c b H 7 y F J 0 Z x U d B J E 0 J j e b q 7 4 k M h 0 q l o W c 6 Q 6 z H a t m b i f 9 5 v U Q H 1 2 7 G R J x o K s h i U Z B w p C M 0 + x / 5 T F K i e W o I J p K Z W x E Z Y 4 m J N i m V T A j O 8 s u r p F 2 r O p f V + n 2 9 0 r j J 4 y j C C Z z C O T h w B Q 2 4 g y a 0 g E A E z / A K b 5 a 2 X q x 3 6 2 P R W r D y m W P 4 A + v z B 7 7 w k O o = < / l a t e x i t > " n U Z 7 S w 1 U m a e V r U D O z 8 v P e Y r o k p I = " > A A A B 7 X i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R i 8 c K 9 g P a U D a b T b t 2 s x t 2 N 2 I I / Q 9 e P C j i 1 f / j z X / j t s 1 B W x 8 M P N 6 b Y W Z e k H C m j e t + O 6 W V 1 b X 1 j f J m Z W t 7 Z 3 e v u n / Q 1 j J V h L a I 5 F J 1 A 6 w p Z 4 K 2 D D O c d h N F c R x w 2 g n G N 1 O / 8 0 i V Z l L c m y y h f o y H g k W M Y G Ol d j j I n 7 L J o F p z 6 + 4 M a J l 4 B a l B g e a g + t U P J U l j K g z h W O u e 5 y b G z 7 E y j H A 6 q f R T T R N M x n h I e 5 Y K H F P t 5 7 N r J

FIG. 3 .
FIG.3.Real space structure of the two-magnon bound state wave function on the triangular lattice in the l = 2 angular momentum channel.Shown are the wave functions ψ1,2(ρ) given in Eq.(30), where ρ = r−r is the distance between the magnons.Doubly degenerate solution for the exciton wavefunction on the triangular lattice in the weak coupling regime, shown in real space [ψ(ρ) = (1/2π) dq ψqe iq•ρ , where ρ is the vector connecting the two magnons].The values slowly decay in space due to the finite exciton size (not shown).

FIG. 4 .
FIG.4.The value of the interaction parameter η = J z i /Ji required to create an exciton with the binding εq = 0.05J1 for different values of the center-of-mass momentum K, taking values along the Γ-K-M -Γ contour in the Brillouin zone (horizontal axis).The minimal η corresponds to the K = 0 case (Γ point), which makes it the leading instability channel and justifies considering zero center-of-mass momentum excitons in the rest of the paper.

FIG. 6 . 10 FIG. 7 .
FIG.6.The width of the Potts-nematic phase (D * c − Dc)/J1 obtained from exciton binding energy as a function of interaction strength η = J z i /Ji for the honeycomb lattice model.Solid line shows numerical solution, while the dashed one is for the analytical perturbative result, valid at η 1.The prefactor in the latter was fitted.

2 FIG. 8 .
FIG.8.Critical values of the interaction parameter ηc on the face-centered cubic (FCC) lattice as a function of spin couplings ratio |J2|/J1.Formation of Potts-nematic excitons and, hence, the vestigial phase occurs only at η = J z i /Ji > ηc due to the corresponding property of the three-dimensional potential well problem.