Spin-orbit interaction driven dimerization in one dimensional frustrated magnets

TL;DR

This paper investigates how spin-orbit interaction influences one-dimensional frustrated magnets, revealing that it induces a bond density wave that breaks translational symmetry and opens an excitation gap.

Contribution

It demonstrates that spin-orbit coupling causes spontaneous bond density wave formation and symmetry breaking in 1D frustrated magnets with nematic fluctuations.

Findings

Spin-orbit interaction induces a bond density wave.

The bond density wave breaks translational symmetry.

A gap opens in the excitation spectrum due to this effect.

Abstract

We study the effect of spin-orbit interaction on one-dimensional U(1)-invariant frustrated magnets with dominant critical nematic fluctuations. The spin-orbit coupling explicitly breaks the U(1) symmetry of arbitrary global spin rotations about the high-symmetry axis down to $Z_2$ (invariance under a $\pi$-rotation). Given that the nematic order parameter is invariant under a $\pi$-rotation, it is relevant to ask if other discrete symmetries can be spontaneously broken. Here we demonstrate that the spin-orbit coupling induces a bond density wave that spontaneously breaks the translational symmetry and opens a gap in the excitation spectrum.

Figures

Equations

$$\mathcal{H}$$

$$+\alpha \sum _j{S}_j^x{S}_{j+1}^x,$$

$$\Re ⟨S_r^{+}{S}_{r+1}^{+}⟩=⟨S_r^x{S}_{r+1}^x-S_r^y{S}_{r+1}^y⟩\equiv ⟨{\mathcal{O}}^a(r)⟩,$$

$$\Im ⟨S_r^{+}{S}_{r+1}^{+}⟩=⟨S_r^x{S}_{r+1}^y+S_r^y{S}_{r+1}^x⟩,$$

$$\mathcal{F}=r\left(|{\Psi }_Q{|}^2+|{\Psi }_{\bar{Q}}{|}^2\right)+u\left({\Psi }_Q^{\ast }{\Psi }_{\bar{Q}}^{\ast }+{\Psi }_{\bar{Q}}{\Psi }_Q\right),$$

$${\mathcal{O}}^a$$

$${\mathcal{O}}^b$$

$${\mathcal{O}}^a(r)$$

$${\mathcal{O}}^b(r)$$

$$n_{pair}(r)$$

$$={\Psi }_0^2+4{\Psi }_{\pi }^2+4(-1{)}^r{\Psi }_0{\Psi }_{\pi }\cos (\theta ).$$

$$⟨S_r^z{S}_{r+1}^z⟩\simeq ⟨{\Psi }_r{⟩}^{\ast }⟨{\Psi }_r⟩+const..$$

$$S_j^{+}$$

$$S_j^{-}$$

$$S_j^z$$

$$\mathcal{H}={\mathcal{H}}_0+{\mathcal{H}}_{int},$$

$${\mathcal{H}}_0$$

$$+i\frac{\alpha }{4}\sum _k\sin k\left(c_k^{†}{c}_{-k}^{†}-c_{-k}{c}_k\right).$$

$${\mathcal{H}}_{int}$$

$$U(Q,q,p)$$

$$+2J_2\sin (2q)\sin (2p).$$

$$c_k$$

$$c_k^{†}$$

$$u_k$$

$$v_k$$

$${ϵ}_k=J_1\cos k+J_2\cos 2k-\left(h+J_1+J_2\right).$$

$${\mathcal{H}}_0=\sum _k{E}_k{\alpha }_k^{†}{\alpha }_k+E_0.$$

$$i{G}_0(k,\omega )$$

$${\mathcal{H}}_{int}^N$$

$${\Gamma }_Q^{(0)N}(q,p)$$

$$-4U(p+q,\frac{Q-p+q}{2},\frac{Q+p-q}{2})\left(u_{\frac{Q}{2}+p}{v}_{\frac{Q}{2}-p}{v}_{\frac{Q}{2}-q}{u}_{\frac{Q}{2}+q}+u_{\frac{Q}{2}-p}{v}_{\frac{Q}{2}+p}{v}_{\frac{Q}{2}+q}{u}_{\frac{Q}{2}-q}\right).$$

$${\mathcal{H}}_{int}^A$$

$${\Gamma }_Q^{(0)A}(q,p)$$

$$-2U(p+q,\frac{Q-p+q}{2},\frac{Q+p-q}{2})\left(u_{\frac{Q}{2}+p}{v}_{\frac{Q}{2}-p}{u}_{\frac{Q}{2}-q}{v}_{\frac{Q}{2}+q}+v_{\frac{Q}{2}+p}{u}_{\frac{Q}{2}-p}{v}_{\frac{Q}{2}-q}{u}_{\frac{Q}{2}+q}\right)$$

$$+2U(p-q,\frac{Q+p+q}{2},\frac{Q-p-q}{2})\left(v_{\frac{Q}{2}+p}{u}_{\frac{Q}{2}-p}{u}_{\frac{Q}{2}-q}{v}_{\frac{Q}{2}+q}+u_{\frac{Q}{2}+p}{v}_{\frac{Q}{2}-p}{v}_{\frac{Q}{2}-q}{u}_{\frac{Q}{2}+q}\right).$$

$${\Gamma }_{\omega Q}^N(q,p)$$

$${\Gamma }_{\omega Q}^A(q,p)$$

Full text

I Introduction

Frustrated magnetism is a continuous source of exotic states of matter that challenge the existing characterization probes ref:review1 ; ref:review2 . Once quantum fluctuations melt the traditional magnetic long-range order, it often happens that the remaining liquid or multipolar orderings do not couple directly to the usual experimental probes. A simple example is the spin nematic phase proposed to be the ground of the one-dimensional (1D) $J_{1}-J_{2}$ Heisenberg model near its saturation field  ref:nematic1 ; ref:nematic2 ; ref:nematic3 ; ref:nematic4 ; ref:nematic5 ; ref:nematic5b ; ref:nematic5c ; ref:nematic5d ; ref:nematic6 ; ref:nematic7 ; ref:nematic8 ; ref:nematic9 ; ref:nematic10 ; ref:nematic11 . This phase arises from a Bose-Einstein condensation of magnon pairs right below the saturation field $h_{\rm sat}$ ref:nematic1 ; ref:nematic2 ; ref:nematic3 ; ref:nematic4 ; ref:nematic5 ; ref:nematic5b ; ref:nematic5c ; ref:nematic5d . The attractive magnon-magnon interaction is generated by a ferromagnetic (FM) nearest neighbor (NN) exchange ($J_{1}<0$), which competes against an antiferromagnetic (AFM) next-nearest-neighbor (NNN) exchange $J_{2}>0$. $J_{2}$ must be bigger than $|J_{1}|/4$ for the zero field ground state not to be ferromagnetic.

Several quasi-1D materials are approximately described by the $J_{1}-J_{2}$ model with FM and AFM exchange interactions $J_{1}$ and $J_{2}$, respectively. Known examples include Rb2Cu2Mo3O12, ref:1DmaterialA1 ; ref:1DmaterialA2 , LiCuVO4, ref:1Dmaterial2 ; ref:1Dmaterial3 ; ref:1Dmaterial4 ; ref:1Dmaterial6 ; LiCuVO4 ; LiCuVO4b ; ref:1Dmaterial7 ; ref:1Dmaterial8 ; ref:1Dmaterial9 ; ref:1Dmaterial10 ; ref:1Dmaterial11 LiCuSbO4,  LiCusbO4 ; grafe2016 PbCuSO4(OH)2, PbCuSO4OH2 which span a wide range of $J_{2}/\rvert J_{1}\rvert$ values. However, a direct experimental observation of the predicted nematic ordering is rahter challenging. sato09 ; sato11 Given the symmetry of the order parameter, the nematic spin ordering is expected to induce a local quadrupolar electric moment via the always present spin-orbit coupling combined with the lattice anisotropy. However, * the spin anisotropy induced by this combination has the additional effect of breaking the global U(1) symmetry of spin rotations along the magnetic field axis down to a finite group*. For most of the known compounds, this group is not bigger than $Z_{2}$ for any direction of the applied magnetic field (only C2 rotation axes). This fact raises another concern because the nematic order parameter $\langle S_{r}^{+}S_{r+1}^{+}\rangle$ does not break the remaining $Z_{2}$ group. In other words, the nematic order parameter is invariant under $\pi$-rotations. This simple observation implies that if some form of ordering still exists right below the high field paramagnetic phase of these compounds, it should not be called “nematic ordering”. Nevertheless, the dominant nematic susceptibility of the U(1) invariant model may still induce additional symmetry breaking in the presence of spin anisotropy. If this is the case, it is important to identify those discrete symmetries.

In this work we investigate the relevant effect of spin-orbit interaction on the 1D frustrated $J_{1}-J_{2}$ model. Based on our previous considerations, there are two possible scenarios: i) The field-induced transition from a quantum paramagnet to the nematic phase is replaced by a crossover (no discrete symmetry breaking); ii) The field-induced transition from a quantum paramagnet to the nematic phase is replaced by a discrete symmetry breaking. We can anticipate that the problem under consideration belongs to the second case because it is known that the magnon-pairs condense at a finite wave vector $\pm Q$. In other words, the nematic ordering breaks the translational symmetry, which is not affected by the inclusion of the spin-orbit interaction. In addition, the system has an additional $Z_{2}$ symmetry besides the $\pi$-rotation about the $z$-axis. This symmetry is the product of two operations, ${\cal T}{\cal R}(\pi)$, where ${\cal T}$ is the time reversal operator and ${\cal R}(\pi)$ is a $\pi$-rotation operator about an axis perpendicular to the field direction. As we explain below, the real part of the nematic order parameter, ${\cal O}^{a}$, preserves this $Z_{2}$ symmetry, while the imaginary part, ${\cal O}^{b}$, does not. Fig. 1 (a) shows the the different broken symmetry states associate with the stabilization of the real or imaginary part of the nematic order parameter.

In the following sections, we derive a phenomenological Ginzburg-Landau theory that is complemented by microscopic analytical and numerical calculations. Moreover, we demonstrate that the combination of a divergent nematic susceptibility and spin anisotropy stabilizes the real component of the nematic order parameter, ${\cal O}^{a}$, which in turn produces bond dimerization in most of the phase diagram. This bond ordering is accompanied by an orthorhombic distortion of the surrounding oxygen octahedron, as it is schematically shown in the upper panel of Fig. 1 (b). In contrast, the imaginary component of the “nematic” order parameter, ${\cal O}^{b}$, does not produce bond dimerization. If stabilized by other mechanisms, this phase would produce a local orthorhombic distortion of the surrounding oxygen octahedron along diagonal directions, as illustrated in the lower panel of Fig. 1 (b).

Given that bond dimerization couples to the lattice via magneto-elastic interaction and lowers the space group of the material under consideration, the combined effect of high spin nematic susceptibility and spin anisotropy can in principle be detected with X-rays. In addition, the incommensurate bond-density wave expected for smaller values of $J_{2}/|J_{1}|$ should lead to a double-horn shape of the nuclear magnetic resonance (NMR) line. These conclusions can shed light on the search for the spin “nematic ordering” predicted on the basis of a U(1) invariant $J_{1}-J_{2}$ Heisenberg model.

The structure of this paper is as follows. In Sec. II, we introduce a simple model Hamiltonian in which the U(1) symmetry is reduced to Z2 via the inclusion of an Ising term (symmetric anisotropy). In Sec. III we consider a simple Ginzburg-Landau (GL) theory, which describes the possible scenarios. The phenomenological input parameters of the GL theory are calculated in Sec. IV by means of an analytical approach to the microscopic Hamiltonian. The results of this analytical approach are confirmed by numerical Density Matrix Renormalization Group (DMRG) calculations presented in Sec. V. The general implications of our results for experimental studies of unconventional magnetic ordering in quasi-1D frustrated compounds are discussed in Sec. VI.

II Model Hamiltonian

We consider a spin-$1/2$ anisotropic Heisenberg model on a 1D chain with ferromagnetic nearest-neighbor exchange, $J_{1}<0$, and antiferromagnetic next nearest neighbor exchange $J_{2}>0$:

[TABLE]

The last term is an Ising interaction that in real materials arises from the combined effect of spin-orbit coupling and lattice anisotropy. This term reduces the U(1) symmetry of continuous spin rotations about the field axis to $Z_{2}$. We note that the U(1) symmetry is restored if the magnetic field is applied along the $x$-direction. However, in real quasi-1D materials, the U(1) symmetry is not present for any direction of the magnetic field because of the interaction with other chains ref:nematic1 . In spite of these considerations, the U(1) invariant model has been invoked to describe various quasi-1D transition metal compounds. sato13

Based on the U(1) invariant model ($\alpha=0$), several authors proposed that nematic quasi-long range ordering should be observed right below the saturation field $h_{\rm sat}$ ref:nematic1 ; ref:nematic2 ; ref:nematic3 ; ref:nematic4 ; ref:nematic5 ; ref:nematic5b ; ref:nematic5c ; ref:nematic5d . $h_{\rm sat}$ is finite only for $J_{2}>|J_{1}|/4$ because the zero field ground state is ferromagnetic for $J_{2}\leq|J_{1}|/4$. The nematic ordering corresponds to a Bose-Einstein condensation of two-magnon bound states with a local order parameter $\langle S_{j}^{+}S_{j+1}^{+}\rangle$. The attractive magnon-magnon interaction is provided by the ferromagnetic nearest-neighbor exchange $J_{1}$. The ratio $J_{2}/\lvert J_{1}\rvert$ controls the total momentum $\pm Q$ of the two-magnon bound state. $Q$ is incommensurate for small values of $J_{2}/\lvert J_{1}\rvert$ and it is equal to $\pi$ for $J_{2}/\lvert J_{1}\rvert\geq 0.375$.

In general, the continuous SU(2) symmetry of the Heisenberg interaction is broken down to a discrete symmetry group in real materials ref:1DmaterialA2 ; ref:1Dmaterial2 . Even for an idealized single-chain system, the exchange interaction turns out to be anisotropic, instead of SU(2) invariant, once the spin-orbit interaction is included. This is so because an isolated chain has only one symmetry axis parallel to the chain direction ($x$-direction in our notation). In other words, the exchange interaction between spin components parallel to the chain direction is different from the exchange interaction between the spin components perpendicular to the chain direction, as it is clear from the $\alpha$-term in Eq. (1). Consequently, the pure 1D Hamiltonian has only discrete point group symmetries if the external magnetic field is not parallel to the chain direction. For the case under consideration (magnetic field perpendicular to the chain direction), the U(1) symmetry of the Heisenberg model is reduced to a discrete $Z_{2}$ symmetry corresponding to a $\pi$-rotation about the $z$-axis ${\cal R}_{z}(\pi):\;S_{r}^{z}\rightarrow S_{r}^{z},S_{r}^{x,y}\rightarrow-S_{r}^{x,y}$. Correspondingly, the spin components $S_{r}^{\pm}$ transform into $e^{\pm i\phi}S_{r}^{\pm}$ under a rotation by $\phi$ about the $z$-axis. This means that the nematic order parameter, $\langle S_{r}^{+}S_{r}^{+}\rangle$, transforms into $e^{i2\phi}\langle S_{r}^{+}S_{r}^{+}\rangle$, implying that it is invariant under $\pi$-rotations, as expected for a director. The inclusion of spin anisotropy then forces us to reconsider the problem because other symmetries (different from rotations) have to be invoked to characterize the phase that replaces the nematic quasi-long range ordering.

Besides the above mentioned $\pi$ rotation about the $z$-axis, the Hamiltonian of Eq. (1) is invariant under the product of the time reversal operation and a $\pi$-rotation about the $y$-axis, ${\cal T}{\cal R}_{y}(\pi)$, which changes the sign of the $y$ spin component: $S_{r}^{x,z}\rightarrow S_{r}^{x,z},S_{r}^{y}\rightarrow-S_{r}^{y}$. The real-part of the nematic order parameter,

[TABLE]

remains invariant under this transformation. In contrast, the imaginary part,

[TABLE]

changes sign. Finally, the nematic order parameter breaks the translational symmetry because the magnon-pairs condense at a finite momentum $\pm Q$. This symmetry is then expected to break spontaneously for $\alpha\neq 0$, as long as $Q$ is commensurate.

Based on this simple symmetry analysis, the spin anisotropy should stabilize a state that breaks the translational symmetry (in a strong or a weak sense) and select either the real or the imaginary part of the original nematic order parameter. Only one of these two components should be selected, as supposed to some linear combination, because they belong to different irreducible representations of the point group of $\mathcal{H}$.

III Ginzburg-Landau theory

The attractive interaction between magnons arising from the ferromagnetic nearest-neighbor interaction, leads to two-magnon bound states for $J_{2}/|J_{1}|>1/4$. The minimum energy of the two-magnon bound state is achieved for a finite value, $\pm Q$, of the center of mass momentum. The two-magnon bound states condense for $h<h_{c}$ (note that $h_{c}$ is higher than the field required to close the single-magnon gap). The two-magnon condensate is characterized by a two component complex order parameter $\Psi_{\pm Q}$ (macroscopic wave function of condensate) whenever $Q\neq-Q$. The spin-orbit interaction generates an effective coupling between these two components, as it can be inferred from the lowest order expansion of the Ginzburg-Landau free energy:

[TABLE]

where $\bar{Q}=-Q$. Due to the $Z_{2}$ symmetry restriction, the complex field $\Psi_{\pm Q}$ is fixed up to a phase factor $\pm 1$. We have also assumed that $u$ is real based on the underlying microscopic theory. We will first assume $Q=\pi$ which is the condensation wave vector for $J_{2}>J_{2c}$ ($J_{2c}\simeq 0.375\rvert J_{1}\rvert$ for $\alpha=0$). Given that $Q=\pi$ is invariant under spatial inversion, we have $\Psi_{\pi}=\Psi_{-\pi}$. Then, upon minimization of the free energy, we obtain a real order parameter $\Psi_{\pi}$ for $u<0$, and a purely imaginary order parameter for $u>0$. In the original spin language, we have $\Psi_{\pi}=\frac{1}{N}\sum_{r}e^{i\pi r}\langle S_{r}^{-}S_{r+1}^{-}\rangle$, whose real and imaginary parts are

[TABLE]

A real order parameter only breaks the translational symmetry, while an imaginary order parameter breaks additional symmetries, such as, ${\cal T}{\cal R}_{y}(\pi):\;S_{r}^{x,z}\rightarrow S_{r}^{x,z},S_{r}^{y}\rightarrow-S_{r}^{y}$. In both cases, the system should develop long-range ordering at $T=0$ because only discrete symmetries are broken. We note that the spin anisotropy corresponds to a uniform nematic field that couples linearly to the uniform component of the nematic order parameter $\Psi_{0}=\frac{1}{N}\sum_{r}\langle S_{r}^{-}S_{r+1}^{-}\rangle$, implying that $\Psi_{0}$ becomes finite for a finite $\alpha$. As we will see next, the interference between $\Psi_{0}$ and the $\pi$ component, $\Psi_{\pi}$, of the real part of the order parameter leads to a real space modulation (dimerization) of the expectation value of nearest-neighbor bond operators.

In general, the nematic order parameter is a complex number, $|\Psi_{\pi}|e^{i\theta}$, where $\theta=0$ and $\theta=\frac{\pi}{2}$ correspond to ${\cal O}^{a}$ and ${\cal O}^{b}$, respectively. The real space version of these order parameters is obtained via a Fourier transformation,

[TABLE]

which gives the real and imaginary parts of $\Psi_{r}=\langle S_{r}^{-}S_{r+1}^{-}\rangle$. It is clear that the amplitude of the real component, ${\cal O}^{a}(r)$, is modulated in real space for finite values of the spin-orbit interaction ($\Psi_{0}\neq 0$). In contrast, only the phase is modulated for $\alpha=0$. In other words, the spin-orbit coupling induces a dimerized bond ordering if the real component of the original nematic order parameter is selected. This interference between the [math] and $\pi$ components of ${\cal O}^{a}$ also leads to a magnon pair density wave:

[TABLE]

It follows that the long range ${\cal O}^{a}$ ordering is accompanied by another bond ordering associated with the longitudinal spin component

[TABLE]

This is just the usual bond dimerization that appears in spin-Peierls systems. Miller82 Indeed, Eq. (10) implies that the usual bond order parameter, $\langle{\bf S}_{r}\cdot{\bf S}_{r+1}\rangle$, must also exhibit dimerization.

The condensation wave vector becomes incommensurate ($Q<\pi$) for smaller values of $J_{2}/|J_{1}|$ (about $0.375$ for $\alpha=0$). In this case we need to consider a two component order parameter with phases $e^{i\theta_{\pm}}$ for $\pm Q$. Minimization of the free energy leads to $\theta_{+}+\theta_{-}=0\text{ mod }(2\pi)$ for $u<0$ and to $\phi_{+}+\phi_{-}=\pi\text{ mod }(2\pi)$ for $u>0$. The complex order parameter does not have a fixed phase because of the additional $U(1)$ phase factor $e^{\pm iQr}$, arising from translational symmetry. This $U(1)$ symmetry precludes long-range order for the single-chain problem. The free energy minimization also leads to the same amplitude for both components of the order parameter: $|\Psi_{Q}|=|\Psi_{\bar{Q}}|$, implying that the ground state must exhibit quasi-long range bond density wave order ${\cal O}^{a}(r)=\Psi_{0}+2|\Psi_{Q}|\cos(Qr)$ or ${\cal O}^{b}(r)=2|\Psi_{Q}|\sin(Qr)$.

IV Microscopic theory

Our discussion in the previous section indicates that two different kinds of bond order can be induced by spin anisotropy. To determine which order parameter is selected we need to consider the underlying microscopic theory. To this end, we use the Jordan-Wigner transformation to reformulate the spin Hamiltonian (1) as a model for interacting spinless fermions:

[TABLE]

where $c_{j}$ is the fermionic spinless operator which represents a spin flip on site $j$. The fermionic Hamiltonian is defined as

[TABLE]

where

[TABLE]

For $\alpha=0$, the single-particle spectrum, corresponding to single-magnon excitations, has a minimum at $k_{0}=\arccos(-\frac{J_{1}}{4J_{2}})$ if $4J_{2}>|J_{1}|$ and at $k_{0}=0$ if $4J_{2}<|J_{1}|$. Fig. 2 shows the evolution of the single-particle spectrum with increasing $J_{2}/\rvert J_{1}\rvert$. The spin-orbit interaction ($\alpha\neq 0$) breaks the $U(1)$ symmetry and the fermion number is no longer conserved. The interacting part of the Hamiltonian is

[TABLE]

where

[TABLE]

The effective attractive interaction between nearest neighbor sites, $J_{1}+2J_{2}\cos Q$, is maximized at $Q=\pi$. Therefore, the lowest energy “two-magnon” bound state is expected to have this momentum for large enough $J_{2}/|J_{1}|$.

In the isotropic limit, the nematic phase corresponds to a magnon pair condensate. This state can be approximated by a coherent state built with the two-particle wave function of the bound state. In presence of magnetic anisotropy, the particle number is not conserved because particles can be created or annihilated in pairs. The Hamiltonian eigenstates can then be grouped into two categories based on the particle number parity. The following analysis assumes that the magnetic anisotropy is weak ($\alpha\ll 1$), implying that the finite particle density induced by the spin anisotropy term above and at the critical field, $h_{c}$ can be made arbitrarily small for small enough values of $\alpha$. This condition also guarantees that two magnon condensation is still the dominant instability for $\alpha\neq 0$. For this reason, we will still refer to the bound state as a “two-magnon” bound states (adiabatic continuation) and we will use the two-particle Green’s function to compute the energy of the two-magnon modes. The two-particle Green’s function is obtained from the two-particle scattering amplitude by solving the corresponding Bethe-Salpeter equation.

IV.1 Bogoliubov representation

The above non-interacting fermionic Hamiltonian can be diagonalized with a Bogoliubov transformation,

[TABLE]

where

[TABLE]

with $E_{k}=\sqrt{\epsilon_{k}^{2}+4\Delta_{k}^{2}}$, $\Delta_{k}=\frac{\alpha}{4}\sin k$ and

[TABLE]

The diagonal Hamiltonian

[TABLE]

leads to the non-interacting Green’s function

[TABLE]

The next step is to write the interaction vertex in terms of the Bogoliubov quasi-particle operators. The normal interaction term is

[TABLE]

with a normal interaction vertex

[TABLE]

We can verify that $\Gamma_{Q}^{(0)N}(q,p)=\Gamma_{Q}^{(0)N}(-q,-p)$ due to the fermionic statistics. Furthermore, $\Gamma_{Q}^{(0)N}(q,p)=\Gamma_{\bar{Q}}^{(0)N}(-q,-p)=\Gamma_{\bar{Q}}^{(0)N}(q,p)$ because of inversion symmetry. The interaction vertex has been symmetrized with respect to the exchange of external lines. The anomalous interaction terms of the form $\alpha^{\dagger}\alpha^{\dagger}\alpha^{\dagger}\alpha^{\dagger}$ and $\alpha\alpha\alpha\alpha$ are

[TABLE]

with an anomalous interaction vertex

[TABLE]

We can verify that $\Gamma_{Q}^{(0)A}(q,p)=\Gamma_{\bar{Q}}^{(0)A}(p,q)=\Gamma_{Q}^{(0)A}(-q,-p)$ due to fermionic statistics and $\Gamma_{Q}^{(0)A}(q,p)=\Gamma_{\bar{Q}}^{(0)A}(-q,-p)=\Gamma_{\bar{Q}}^{(0)A}(q,p)$ due to inversion symmetry. This interaction vertex has also been symmetrized with respect to the exchange of external lines.

The remaining anomalous terms ($\alpha^{\dagger}\alpha^{\dagger}\alpha^{\dagger}\alpha$ and $\alpha\alpha\alpha\alpha^{\dagger}$) will not be considered because they give subdominant contributions (in an expansion in powers of $\alpha$) to the effective anomalous interaction vertex shown in Fig. 3. We note that the nature of the order parameter is determined by the effective anomalous interaction vertex because the effective normal vertex preserves the U(1) symmetry of ${\cal H}(\alpha=0)$.

Finally, we also note the existence of additional diagrams which are obtained by exchange of the external legs of the diagrams shown in Fig. 3. These diagrams are of the same order in an expansion in powers of $\alpha$. However, for the purpose of calculating the two magnon bound state, they can be neglected because they give contributions that remain regular in the proximity of the critical field $h=h_{c}$.

IV.2 Bethe-Salpeter equation

In the dilute limit, the scattering amplitude can be calculated by summing up the ladder diagrams shown in Fig. 3. This sum leads to the Bethe-Salpeter equation BSE :

[TABLE]

where $\Gamma_{Q}^{(0)N/A}$ is the bare scattering amplitude of the normal and anomalous type. The energy of the magnon pair bound state can be extracted from the poles of the scattering amplitude. For a wide range of $J_{2}/|J_{1}|$ values, the bound state dispersion has its minimum at $Q=\pi$. The energy of the two-magnon bound state increases with $\alpha$, implying that the critical field for the “two-magnon” condensation decreases relative to the saturation field of the isotropic ($\alpha=0$) Hamiltonian.

For the isotropic Heisenberg model, the condensate wave function has a $U(1)$ phase freedom. The spin anisotropy reduces this freedom to $Z_{2}$. The Ginzburg-Landau theory tells us that the phase of the macroscopic wave function can either be real or imaginary depending on the sign of the effective anomalous coupling parameter $u$. In this section we determine the phase of the wave function and also include a microscopic calculation of the parameter $u$. These properties are enclosed in the scattering amplitude obtained from the solution of the Bethe-Salpeter equation. Our analysis shows that wave function of the magnon pair condensate is real, implying that the dominant order parameter is ${\cal O}^{a}$.

IV.2.1 Wave function of the bound state

We start by introducing the two-magnon Green’s function, which can be easily obtained through the scattering amplitude $\Gamma_{\omega Q}^{N/A}(q,p)$:

[TABLE]

where $G_{0}^{(2)}$ is the non-interacting two particle Green’s function:

[TABLE]

The Lehmann representation shows explicitly that the two particle Green’s function has the following singular behavior near the pole of bound state, $\omega\sim\omega_{B}$:

[TABLE]

where the regular terms come from higher excited states and $\omega_{B}>0$ is the bound state energy relative to the ground state. The projection of the bound state wave function on the two-magnon sector is obtained from the residue of the pole,

[TABLE]

where $|B\rangle$ is the ket of the bound state and $|G\rangle$ the ground state. Due to the anomalous interaction arising from the spin anisotropy term, the bound state wave function is a linear combination of states with different particle number. The poles of the two particle Green’s function given by Eqs. (34) and (35) are obtained after inserting the scattering amplitude, $\Gamma_{\omega Q}^{N/A}(q,p)$, which results from the Bethe-Salpeter equation. According to Eqs. (39) and (40), the bound state wave function is then obtained by extracting the residue near the pole $\omega_{B}$. For $Q=\pi$, the bound state wave functions, $\psi_{\pi}(p)$ and $\phi_{\pi}(p)$, are even under spatial inversion. This result is found to be correct for all $J_{2}/|J_{1}|$ ratios and for any value of the bound state energy $\omega_{B}$, indicating that the broken symmetry state below the critical field must preserve the spatial inversion symmetry.

As the system approaches the critical field corresponding to the onset of the “two-magnon” condensate, $\omega_{B}\rightarrow 0$, the normal and the anomalous Green’s function become exactly the same (see Fig. 4). Consequently, the particle pair wave function, $\psi_{Q}(p)$, is exactly the same as the hole pair wave function $\phi_{\bar{Q}}^{*}(p)$. This is a manifestation of an emergent particle-hole symmetry at zero energy, which sets a constraint on the phase of the bound state wave function. If we describe the condensate state as a coherent state built with the two-body bound state wave function, the phases $\theta_{\pm}$ for the macroscopic components $\Psi_{\pm Q}$ are given by $\psi_{Q}(p)$ and $\phi_{\bar{Q}}(p)$, respectively. The relationship $\psi_{Q}(p)=\phi_{\bar{Q}}^{*}(p)$ indicates that $\theta_{+}+\theta_{-}=0$, which leads to a real order parameter $\langle S_{r}^{+}S_{r+1}^{+}\rangle$ in real space.

Beyond the condensation point, the new ground state, $|G\rangle$, is characterized by the order parameter $\Phi_{Q}(q)=\langle G|\alpha_{\frac{Q}{2}+q}\alpha_{\frac{Q}{2}-q}|G\rangle$. The coherent representation enables us to identify the bound state wave functions $\psi_{Q}(p)$ and $\phi_{\bar{Q}}(p)$ with the two component order parameter, $\Phi_{Q}(p)$ and $\Phi_{\bar{Q}}(p)$, that was discussed in Section III. This correspondence leads to the self-consistent equation for the order parameter based on the Bethe-Salpeter equation, from which one can straightforwardly confirm which order is favored. The analysis becomes more transparent by adopting an equivalent but more straightforward approach. We just introduce a small pairing field term into the Hamiltonian, which couples to the order parameter:

[TABLE]

The renormalized pairing fields $h_{Q}^{R},h_{\bar{Q}}^{R}$ are indicated by the ladder series of vertex corrections in Fig. 5. The pairing susceptibility diverges at the condensation point, implying that the order parameter develops spontaneously beyond this point, i.e., in absence of the pairing fields $h_{Q}$ and $h_{\bar{Q}}$. The order parameter in momentum space, $\Phi_{Q}(q)=\langle\alpha_{\frac{Q}{2}+q}\alpha_{\frac{Q}{2}-q}\rangle$, can be calculated as $\Phi_{Q}(q)=-h_{Q}^{R}(q)/\Omega_{q}$. Therefore, the ladder series of vertex corrections in Fig. 5 leads to the following self-consistent equation:

[TABLE]

where $\Omega_{Q}(q)=E_{{Q\over 2}+q}+E_{{Q\over 2}-q}$ is the energy of a two-magnon excitation and

[TABLE]

The order parameters $\Phi_{Q}^{\pm}(k)$ become finite when the corresponding matrix in Eq. (42) is singular. For $Q=\pi$, $\Phi_{\pi}^{+}(k)$ coincides with imaginary part of $\Phi_{\pi}(k)$, while $\Phi_{\pi}^{-}(k)$ is the real part. The numerical calculation shows that the order parameter $\Phi_{\pi}(k)$ is purely imaginary for $Q=\pi$ and that it satisfies $\Phi_{Q}(k)=\Phi_{-Q}^{*}(k)$ for $Q<\pi$.

To understand the meaning of this result in real space, we just need to consider the order parameter $\langle S_{r+1}^{-}S_{r}^{-}\rangle$, which is given by $\langle c_{r+1}c_{r}\rangle$ in terms of the Jordan-Wigner fermionic annihilation operators. By applying Fourier and Bogoliubov transformations, we find

[TABLE]

The relationship $\Phi_{Q}(q)=\Phi_{\bar{Q}}^{*}(q)$ implies that the real space order parameter is

[TABLE]

where $\phi_{+}$ is the phase of $\Phi_{Q}(q)$. Therefore, the order parameter $\langle S_{r+1}^{-}S_{r}^{-}\rangle$ is real.

IV.2.2 Microscopic calculation of phenomenological parameters

To provide a microscopic derivation of the phenomenological parameters $r$ and $u$ we express the Ginzburg-Landau free energy in its diagonal form:

[TABLE]

where $\Psi_{\pm}=\Psi_{Q}\pm\Psi_{\bar{Q}}^{*}$ and $\Psi_{Q}$ is the Fourier transform of the bond order parameter $\langle S_{r}^{-}S_{r+1}^{-}\rangle$. The order parameter can be expressed in terms of the fermionic Bogoliubov quasi-particles,

[TABLE]

where $\Phi_{\pm}(q)=\Phi_{Q}(q)\pm\Phi_{\bar{Q}}^{*}(q)$ and $\Phi_{Q}(q)=\langle\alpha_{\frac{Q}{2}+q}\alpha_{\frac{Q}{2}-q}\rangle$. From Eq. (46), we can identify the phenomenological parameters with the inverse of the corresponding static susceptibilities:

[TABLE]

In other words, $\chi_{\pm}$ are the response functions to pairing fields that couple linearly to the $\Psi_{\pm}$ order parameters given in Eq. (47):

[TABLE]

where $\Gamma_{0Q}^{N/A}(q,p)$ is the scattering amplitude at zero frequency, $\Omega_{Q}(q)=E_{\frac{Q}{2}+q}+E_{\frac{Q}{2}-q}$ and $B_{q}^{\pm}=u_{\frac{Q}{2}+q}u_{\frac{Q}{2}-q}\pm v_{\frac{Q}{2}+q}v_{\frac{Q}{2}-q}$. The second term is the non-interacting susceptibility of the Bogoliubov fermions, which is negligible near the critical point where both $\Gamma_{0Q}^{N}(q,p)$ and $\Gamma_{0Q}^{A}(q,p)$ diverge. The finite $u$ value arises from the non-zero anomalous scattering amplitude, $\Gamma_{0Q}^{A}(q,p)$, which forces the two susceptibilities $\chi_{Q}^{\pm}$ to be different:

[TABLE]

Numerically, we always find $u<0$ for different ratios of $J_{2}/J_{1}$, in agreement with our previous discussions.

V Numerical simulations

The above theoretical analysis for small $\alpha$ indicates that the spin anisotropy stabilizes a dimerized ground state. In this section we present DMRG calculations ref:dmrg1 ; ref:dmrg2 for the anisotropic 1D spin Hamiltonian ${\cal H}$, which confirm this analysis. The calculations have been done right below the critical field for chains of $L=160$ spins with open boundary conditions. We used up to 400 states and kept the truncation tolerance below $10^{-8}$ throughout the DMRG iterations. We did 6 full sweeps of finite algorithm of DMRG to get well converged observables. Fig. 6 shows the correlation functions for the real and imaginary parts of the nematic order parameter, ${\cal O}^{a}$, ${\cal O}^{b}$, for the “pair-density" operator $S_{i}^{z}S_{i+1}^{z}$, and for the bond operators ${\bf S}_{i}\cdot{\bf S}_{i+1}$, for different values of $\alpha$. The frustration ratio is taken as $J_{2}/J_{1}=0.9659$, which gives rise to a lowest energy “two-magnon” bound state with center of mass momentum equal to $\pi$. In agreement with the “two-magnon” calculation, the correlation functions of the order parameters ${\cal O}^{a}$ and ${\cal O}^{b}$ oscillate with wave vector $\pi$. The long wave length oscillations are just a consequence of the open boundary conditions and the large spin density wave susceptibility. Note, however, that the incommensurate nature of oscillations precludes the possibility of having long range spin density wave ordering (the incommensurate spin density wave ordering breaks the continuous U(1) symmetry group of translations).

The first column in Fig. 6 with $\alpha=0$ corresponds to the U(1) symmetric case for which all the correlators exhibit the expected power law behavior. The order parameters ${\cal O}^{a}$ and ${\cal O}^{b}$ are connected by a $\pi/4$ spin rotation about the $z$-axis, which is a symmetry of the Hamiltonian for $\alpha=0$. Correspondingly, both correlation functions exhibit an identical power-law decay. The bond and the pair density correlators also exhibit a power law decay with long wave length oscillations, which are magnified by the open boundary conditions ref:chiral2 .

The upper panels of Fig. 6 show that the real component of the nematic order parameter dominates over the imaginary part and develops long range ordering for nonzero $\alpha$. The amplitude of the order parameter increases with $\alpha$. As expected from the previous analysis, the pair density, $S_{i}^{z}S_{i+1}^{z}$, and the bond, ${\bf S}_{i}\cdot{\bf S}_{i+1}$, operators become dimerized as a consequence of the coexistence of uniform and staggered components of the ${\cal O}^{a}$ order parameter. We recall that the uniform component is directly induced by the $\alpha$ term of the Hamiltonian, while the staggered ($\pi$) component is spontaneously generated. The dimerization is identified by Fourier transforming the real space correlation functions shown in Figs. 6(e-h). This information is provided in the insets of Fig. 7, which show a clear peak at $Q=\pi$. It is also clear from Figs. 6(e-h) that the dimerization order parameter becomes more robust upon increasing $\alpha$. Moreover, the insets of Fig. 7 show that the Friedel oscillations of $\langle S_{i}^{z}\rangle$ are strongly suppressed for $\alpha=0.7$, i.e., for a large amplitude of the dimerization order parameter. Consistently with these results, Fig. 7 shows that the dimerized phase (shaded region) becomes stable over a larger window of magnetic field values as we increase $\alpha$. The same figure includes the field dependence of the magnetization curve for different values of $\alpha$. As expected, the slope of the magnetization curve above the critical field increases with $\alpha$ and the magnetization saturates asymptotically for $H\rightarrow\infty$.

VI Discussion and summary

Finally, we discuss the consequences of our theoretical study on the experimental search for nematic phases in various quasi-1D materials. We first note that two previous works studied the role of an XXZ magnetic anisotropy, which preserves the $U(1)$ symmetry of the spin Hamiltonian. Both works concluded that the effect of the U(1) invariant spin anisotropy is to stabilize the nematic phase. Similarly, we are finding that the effect of a U(1) symmetry breaking anisotropy is to stabilize the dimerized phase, which is a “direct descendant" of the nematic ordering (real part of the nematic order parameter).

The spin anisotropy of the quasi-1D compound LiCuVO4 is approximately $\alpha=0.1$. ref:Z2material1 ; ref:1Dmaterial2 ; ref:1Dmaterial3 ; ref:1Dmaterial6 ; LiCuVO4 ; LiCuVO4b ; observation The results of our simulations for this value of $\alpha$ [see Figs. 6 (b) and (f)] clearly show that the bond-bond correlation function develops a visible $\pi$ ordering (dimerization). This is a salient experimental signature that can be detected with X-rays. Rb2Cu2Mo3O12 is another quasi-1D $J_{1}$-$J_{2}$ frustrated magnet ref:1DmaterialA1 ; ref:1DmaterialA2 with $J_{2}/\lvert J_{1}\rvert\simeq 0.33$. This value falls in the “mutipolar" phase (quasi-condensate of $n$-magnon bound states with $n>2$) of the U(1) invariant spin model. Indeed, the nematic phase becomes unstable for $J_{2}/\rvert J_{1}\rvert<0.37$. Although we have not considered the situation of multipolar orderings of higher order than nematic, our simple GL analysis shows that the $\alpha$-term should once again select the real or the imaginary part of the multipolar order parameter.

A more recent experimental work reports evidence of a gapped “nematic ordering" for the quasi-1D $J_{1}-J_{2}$ compound LiCuSbO4. grafe2016 Nuclear magnetic resonance (NMR) measurements indicate the opening of a rather large spin gap above a critical magnetic field value $H_{c1}\simeq 13$T. Conventional explanations for the origin of this spin gap, such as saturation of the magnetization or the presence of a staggered Dzyaloshinskii-Moriya (DM) interaction have been ruled out by the experiments. Moreover, the absence of a kink in the low temperature magnetization vs. field curve is consistent with an XYZ magnetic anisotropy. Based on our analysis, this U(1) symmetry breaking anisotropy is responsible for the the spin gap that is inferred from the NMR measurements. According to Ref. grafe2016 the value of the spin anisotropy is $\alpha\simeq 0.08$. In addition, a $J_{2}/\rvert J_{1}\rvert\simeq 0.28$ ratio is obtained from a fit of the magnetization curve with a pure one-dimensional Hamiltonian. This ratio puts the material outside the region of stability of the nematic phase ref:chiral2 (higher multipolar orderings are expected for $J_{2}/\rvert J_{1}\lesssim 0.367$). However, Ref. LiCusbO4 reports a significantly larger ratio $J_{2}/\rvert J_{1}\rvert\simeq 0.45$, which falls in the nematic phase of the U(1) invariant Hamiltonian. ref:chiral2 . If this $J_{2}/\rvert J_{1}\rvert$ ratio is correct, the gapped high-field spin phase of this compound must correspond to a spin dimerized phase. However, further inspection of this material grafe2016 shows that the lattice is already dimerized (consecutive bonds are not equivalent) implying a continuous crossover between the dimerized spin phase and the high field paramagnetic state.

In summary, we have demonstrated that the spin-orbit interaction has important consequences for the field-induced spin nematic ordering of U(1) invariant frustrated models. The symmetry reduction of ${\cal H}$ due to the presence of the Ising term renders the nematic order parameter inapplicable. However, the real and imaginary parts of the nematic bond order parameter still break discrete symmetries, which can be directly related with observable quantities. Our analytical and numerical results demonstrate that the spin-orbit interaction stabilizes a bond density wave (bond dimerization for $Q=\pi$) which couples to the lattice via the magneto-elastic effect. These results are confirmed by our DMRG simulations. Given that the spin-orbit interaction is ubiquitous in nature and that continuous symmetries are never strictly present in real magnets, our study indicates that the proposed nematic ordering is likely to be replaced by bond dimerization in systems described by a $J_{1}$-$J_{2}$ model with $|J_{2}|/|J_{1}|\gtrsim 0.38$. Even in quasi-1D systems, which are approximately described by a U(1) invariant $XXZ$ model, the application of a magnetic field perpendicular to the chains should induce the dimerized state that we are proposing here.

Acknowledgements.

We thank Zhentao Wang, Shi-Zeng Lin and Yukitoshi Motome for helpful discussions. S. Z. and C.D.B. are supported by funding from the Lincoln Chair of Excellence in Physics and from the Los Alamos National Laboratory Directed Research and Development program. N. K. was supported by the National Science Foundation, under Grant No. DMR-1404375. E. D. was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Materials Sciences and Engineering Division.