Emergence of nematic paramagnet via quantum order-by-disorder and pseudo-Goldstone modes in Kitaev magnets

The appearance of nontrivial phases in Kitaev materials exposed to an external magnetic field has recently been a subject of intensive studies. Here, we elucidate the relation between the field-induced ground states of the classical and quantum spin models proposed for such materials, by using the infinite density matrix renormalization group (iDMRG) and the linear spin wave theory (LSWT). We consider the $K\mathrm{\Gamma }{\mathrm{\Gamma }}^{\prime }$ model, where $\mathrm{\Gamma }$ and ${\mathrm{\Gamma }}^{\prime }$ are off-diagonal spin exchanges on top of the dominant Kitaev interaction $K$. Focusing on the magnetic field along the [111] direction, we explain the origin of the nematic paramagnet, which breaks the lattice-rotational symmetry and exists in an extended window of magnetic field, in the quantum model. This phenomenon can be understood as the effect of quantum order-by-disorder in the frustrated ferromagnet with a continuous manifold of degenerate ground states discovered in the corresponding classical model. We compute the dynamical spin structure factors using a matrix operator based time evolution and compare them with the predictions from LSWT. We, thus, provide predictions for future inelastic neutron scattering experiments on Kitaev materials in an external magnetic field along the [111] direction. In particular, the nematic paramagnet exhibits a characteristic pseudo-Goldstone mode, which results from the lifting of a continuous degeneracy via quantum fluctuations.

Figure 1

Phase diagram of the $K\mathrm{\Gamma }{\mathrm{\Gamma }}^{\prime }$ model in a magnetic field along the [111] axis when (a) ${\mathrm{\Gamma }}^{\prime }=0$ and (b) $-0.06$. In both cases, the Kitaev spin liquid (KSL) is observed only in a small region near $\varphi \to 0$ and $h\to 0$. A large fraction of the phase diagram is occupied by a nematic paramagnet (NP) that spontaneously breaks the $C_3$ lattice-rotation symmetry. The inset illustrates the extension of NP down to zero field and right to the pure $\mathrm{\Gamma }$ model. A long-range zigzag (ZZ) ordered phase appears at $\mathrm{\Gamma }/\left|K\right|\gtrsim 1$ at intermediate fields. The orange dashed line indicates the transition between the fully polarized phase and frustrated ferromagnet based on ansatz (3) in the corresponding classical model. A finite ${\mathrm{\Gamma }}^{\prime }<0$, i.e., (b) ${\mathrm{\Gamma }}^{\prime }=-0.06$, induces ZZ at low fields such that NP occurs at intermediate fields between ZZ and the longitudinally polarized paramagnet (P) that is adiabatically connected to the fully polarized limit. The phase diagram (a,b) is obtained on the rhombic-2 cylinder geometry with $L_{\mathrm{circ}}=10$ (see Fig. 3). [(c)–(f)] Dynamical spin structure factors $\mathcal{S}(k,\omega )$ of P and NP near the upper critical field at two opposite limits of the $K\mathrm{\Gamma }$ model as illustrated by the orange hexagons in (a). In both limits, the magnon band at low energies, which is apparent in P, becomes diffuse in NP, whereas the continuum at higher energies persists across the transition P to NP. The characteristic feature of NP is the increased spectral weight at $\mathrm{\Gamma }$ and at low energies, which is indicative of the pseudo-Goldstone mode. Remarkably, the continuum exhibits a very different structure in both limits.

Figure 2

(a) Canting of the classical spin away from the [111] field, its U(1) manifold (orange line) of degenerate states, and the two sets of azimuthal angles selected by quantum fluctuations (empty and filled circles). (b) Lifting of the classical U(1) degeneracy via quantum order-by-disorder, and the azimuthal dependence of energy, Eq. (4), at $\varphi /\pi =0.17$, ${\mathrm{\Gamma }}^{\prime }=0$ and $h\in \{1.43,1.44,1.45\}$. The average energy $\overline{E}=1/2\pi {\int }_0^{2\pi }d\phi E\left(\phi \right)$ has been subtracted. As the field is tuned, a transition from one set of azimuthal angles to the other occurs. (c) Frustrated ferromagnet (FF, blue) in the phase diagram of the classical $K\mathrm{\Gamma }$ model (${\mathrm{\Gamma }}^{\prime }=0$) subjected to a magnetic field in [111] direction. FF is stabilized between the fully polarized state and other nontrivial ordering patterns, e.g., ZZ and 18-site order [30]. Upon incorporating quantum fluctuations into the classical FF state via LSWT, Eq. (4), FF turns into the nematic paramagnet (NP). The subscripts $\gamma$ and $\alpha \beta$ indicate that quantum fluctuations lift the classical U(1) degeneracy differently. The orange line signifies the transition between ${\mathrm{NP}}_{\gamma }$ and ${\mathrm{NP}}_{\alpha \beta }$.

Figure 3

In order to enable the use of iDMRG to obtain the ground-state wave function in the form of matrix product state (MPS), the two-dimensional honeycomb lattice is mapped to an infinitely long cylinder. The following geometries are used: (a) rhombic with a circumference of $L_{\mathrm{circ}}=6$ and 12 (not shown) sites and (b) its allowed momenta in reciprocal space, [(c) and (d)] rectangular with $L_{\mathrm{circ}}=8$ and 12 (not shown), and [(e) and (f)] rhombic-2 with $L_{\mathrm{circ}}=10$. The geometries are chosen such that they (1) capture the $K$ point in reciprocal space, at which gapless Majorana-Dirac cones of the Kitaev model with isotropic coupling strength occur [7], and (2) the flux decoration of the classical $\mathrm{\Gamma }$ model [47] is not frustrated.

Figure 4

Order parameters ${\mathcal{O}}_{C_3}^{\gamma }$ for each bond $\gamma \in \{x,y,z\}$ at $\varphi /\pi =0.1$ on (a) the rectangular geometry with $L_{\mathrm{circ}}=12$ and (b) the rhombic geometry with $L_{\mathrm{circ}}=12$. The nonzero ${\mathcal{O}}_{C_3}^{\gamma }$ indicates the broken $C_3$ symmetry of NP. Depending on the geometry both sets of nematic states occur. Namely, (a) the rectangular geometry selects ${\mathrm{NP}}_{\gamma }$, while (b) the rhombic geometry selects ${\mathrm{NP}}_{\alpha \beta }$ in a wide range of fields. Both geometries exhibit an intermediate NP phase with states which neither correspond to ${\mathrm{NP}}_{\gamma }$ nor ${\mathrm{NP}}_{\alpha \beta }$.

Figure 5

In the limit $h⟶0$, the quantum ground state obtained by iDMRG depends on the cylinder geometry and the bond dimension $\chi$ of the MPS. The solid line represents a quadratic fit to the energies of MPS that belong to NP, while for dashed lines the MPS exhibits a long-ranged magnetically ordered state. More specifically, we find the following states. At (a) $\varphi /\pi =0.1$ and (b) 0.25, the ground state is either ${\mathrm{NP}}_{\alpha \beta }$ for rhombic unit cells (rho) or ${\mathrm{NP}}_{\gamma }$ for rectangular unit cells (rec). The rhombic-2 (rho-2) geometry cannot be associated with either of them, as it does not exhibit the corresponding anisotropy in bond energies, namely $E_{\alpha }\simeq E_{\beta }<E_{\gamma }$ (${\mathrm{NP}}_{\alpha \beta }$) or $E_{\gamma }<E_{\alpha }\simeq E_{\beta }$ (${\mathrm{NP}}_{\gamma }$) with $a,b,c\in \{x,y,z\}$. Zigzag (ZZ) order is found close in energy for rho-2. (c) In the pure $\mathrm{\Gamma }>0$ model, i.e., $\varphi /\pi =0.5$, cylinders with $L_{\mathrm{circ}}=12$ ($3\mathrm{x}{L}_y$ rec, and $6\mathrm{x}{L}_y$ rho) exhibit transitions from states with finite magnetic moment and long-range order at small $\chi$ into either of the NP phase at large $\chi$. In particular, “5:2 rho-2” exhibits a ZZ order for $\chi \le 1200$, as does “6x2 rho” for $\chi \le 450$, “6x3 rho” exhibits a 12-site order for $\chi \le 256$, “3x1 rec” and “3x2 rec” show a six-site order for $\chi \le 640$ “3x3 rec” exhibits an 18-site order for $\chi \le 640$. This indicates that ${\mathrm{NP}}_{\gamma }$, ${\mathrm{NP}}_{\alpha \beta }$ and several magnetically ordered states observed in the classical model [30] are close in energy.

Figure 6

Dynamical spin structure factor $\mathcal{S}(\mathbf{k},\omega )$ at $\varphi /\pi =0.1$. (a) $h=2.0$ and (b) $h=1.2$ in the high-field paramagnet (P) phase, (c) $h=0.8$ in P right before the transition into the nematic paramagnet (NP), and (d) $h=0.4$ within NP. In (a)–(c), P exhibits two magnon bands, with the lower one remaining sharp but the upper one being obscured by the multimagnon continuum upon decreasing the field. Comparing to the dispersion obtained by LSWT, the lower magnon band (solid line) is shifted downwards in energy and hits zero at the $\mathrm{\Gamma }$ and ${\mathrm{\Gamma }}^{\prime }$ point at a higher critical field than that predicted classically. Within NP, as in (d), the lower magnon band becomes slightly diffusive. A broad continuum appears, ranging from (almost) zero frequency to about $\omega =1.6$. The spectral weight is mostly concentrated at low frequencies near the $\mathrm{\Gamma }$ point, indicating long wavelength excitations corresponding to a pseudo-Goldstone mode. LSWT spectrum in the FF phase, or NP in the quantum limit, is gapless by the Nambu-Goldstone theorem. Some of the spectral weight is redistributed to the $M$-point, signaling that a small perturbation, e.g., ${\mathrm{\Gamma }}^{\prime }$, may be sufficient to drive the system into the zigzag (ZZ) order.

Figure 7

In order to verify the nearly degenerate U(1) manifold of ground states in the quantum model (1), we tilt the magnetic field slightly from the [111] axis by $\vartheta =1^{\circ }$ and then rotate it by varying $\phi$. Plotted are the ground-state energy as obtained with iDMRG for two different geometries with 12 sites circumference each, 6x1 rho [see Fig. 3] and 3x1 rec [see Fig. 3]. The average energy ${\overline{E}}_{\mathrm{GS}}=1/2\pi {\int }_0^{2\pi }d\phi E_{\mathrm{GS}}\left(\phi \right)$ has been subtracted. The magnetic moment follows the rotating field smoothly with a slight deviation up to $7^{\circ }$. As discussed in Sec. 3, the classical U(1) degeneracy is lifted by quantum fluctuations. The induced gap is of the order of $\varepsilon \approx {10}^{-4}$ and increases upon lowering the field. The minima are located at $\phi =0$, $\frac{2}{3}\pi$, and $\frac{4}{3}\pi$ corresponding to ${\mathrm{NP}}_{\alpha \beta }$.

Figure 8

Dynamical spin structure factors $\mathcal{S}(\mathbf{k},\omega )$ complementing Figs. 1–1 and 6 in the main text. Superimposed are the single-magnon bands (solid line) and the onset of the two-magnon continuum (dashed line) as obtained from LSWT. Near the Kitaev limit, $\varphi /\pi =0.05$ (a,b), there exists a broad continuum akin to the KSL and an almost flat band at low energies within P. Its dispersion increases upon tuning up $\varphi$ and reaches a maximum around $\varphi /\pi =0.25$ (e,f). In the pure $\mathrm{\Gamma }$ model, $\varphi /\pi =0.5$ (g,h), the dispersion is still present although reduced. In the entire region $\varphi /\pi =\left(0,0.5\right]$ the lower band attains its minimum at the Brillouin zone center, i.e., the $\mathrm{\Gamma }$ and ${\mathrm{\Gamma }}^{\prime }$ points, and reaches zero frequency at the transition between P and NP. Within NP the main spectral weight is located at low frequencies and long wavelengths, which is consistent with the existence of a pseudo-Goldstone mode. Some spectral weight is distributed around the M point, where the band bends towards lower energies. Thus a small perturbation, e.g., ${\mathrm{\Gamma }}^{\prime }<0$, may trigger a transition into the ZZ order.

Figure 9

The magnon dispersions in the FF phase calculated with LSWT for the $K\mathrm{\Gamma }$ model at $\varphi /\pi =0.1$ and (a) $h=0.42$, (b) 0.4, and (c) 0.389. Upon lowering the field, the excitation gap at the M point decreases and approaches zero, as there is a nearby phase transition to the ZZ order in the classical model.