Finite-temperature behavior of a classical spin-orbit-coupled model for ${\mathrm{YbMgGaO}}_4$ with and without bond disorder

We present the results of finite-temperature classical Monte Carlo simulations of a strongly spin-orbit-coupled nearest-neighbor triangular-lattice model for the candidate $\mathrm{U}\left(1\right)$ quantum spin liquid ${\text{YbMgGaO}}_4$ at large system sizes. We find a single continuous finite-temperature stripe-ordering transition with slowly diverging heat capacity that completely breaks the sixfold ground-state degeneracy, despite the absence of a known conformal field theory describing such a transition. We also simulate the effect of random-bond disorder in the model, and find that even weak bond disorder destroys the transition by fragmenting the system into very large domains—possibly explaining the lack of observed ordering in the real material. The Imry-Ma argument only partially explains this fragility to disorder, and we extend the argument with a physical explanation for the preservation of our system's time-reversal symmetry even under a disorder model that preserves the same symmetry.

Figure 1

The magnetic ${\text{Yb}}^{3+}$ ions form a triangular lattice with a lattice spacing of 3.4 Å. The nearest-neighbor bonds are oriented parallel to the three directions ${\mathit{a}}_i$. The four sublattices $A,B,C,D$ described in the main text are colored differently.

Figure 2

Cross-section with $J_{\pm }=0.915J_{zz}$ of the approximate $T=0$ phase diagram for the Hamiltonian (1) from the Luttinger-Tisza ansatz. In the strong SO regime $|J_{\pm \pm }|\gg 0$, there are six degenerate stripe-ordered ground states. In the “in-plane SO phase” with $J_{\pm \pm }\ll 0$, the spins lie in the $x\text{−}y$ plane and point along the stripes. In the “out-of-plane SO phase” with $J_{\pm \pm }\gg 0$ (and extending to weak negative $J_{\pm \pm }$ for large $J_{z\pm })$, they point perpendicular to the stripes and tilt out of the $x\text{−}y$ plane (and so appear shorter when viewed from above). For weak SO couplings $J_{\pm \pm },J_{z\pm }\ll J_{zz}$, the ground state shows in-plane ${120}^{\circ }$ order with an emergent continuous $U\left(1\right)$ symmetry in spin space. The two black dots in the SO phases represent the Hamiltonians whose finite-temperature behavior is investigated in this paper.

Figure 3

The complex order parameter $\mathrm{\Psi }$ defined in (7) equals a different sixth root of unity for each of the six stripe-ordered ground states in the SO phases. Values of $\mathrm{\Psi }$ depicted in the same color correspond to ground state pairs related by time-reversal symmetry.

Figure 4

Heat capacity per site for clean systems of various sizes in the (a) in-plane and (b) out-of-plane phases.

Figure 5

Order parameter (7) for clean systems of various sizes in the (a) in-plane and (b) out-of-plane phases.

Figure 6

Thermal order parameter $\mathrm{\Psi }$ distributions for the clean $46\times 46$ system in the in-plane phase at temperatures (a) $T=1.0$, (b) $T=1.5$, (c) the critical temperature $T_c=1.56$, (d) $T=1.7$, and (e) $T=2.0$, and for the clean $64\times 64$ system in the out-of-plane phase at temperatures (f) $T=1.0$, (g) $T=1.4$, (h) the critical temperature $T_c=1.49$, (i) $T=1.6$, and (j) $T=2.0$. The circle in each plot represents the unit circle of the complex plane.

Figure 7

Thermal energy distributions for the clean $128\times 128$ systems in (a) the in-plane and (b) the out-of-plane phase. The distributions are single-peaked as expected for a continuous transition, rather than double-peaked as expected for a first-order transition.

Figure 8

Spin structure factors $S\left(\mathit{k}\right)$ defined by (10), and connected correlation lengths in lattice spacings, for the clean $46\times 46$ systems in the in-plane phase [(a) and (b)] and the out-of-plane phase [(c) and (d)], for temperatures $T=1.3$ [(a) and (c)] and $T=1.7$ [(b) and (d)]. The spin structure factors at the critical temperatures look qualitatively similar to those for $T=1.7$, except that a much higher fraction of the total spectral weight lies exactly at the $M$ points. The spectral weight at the $M$ points indicated by black dots is far higher than the plotted scale: they contain $62\%$ and $56\%$ of the total spectral weight for temperature $T=1.3,32\%$ and $34\%$ at the critical temperatures (not shown), and $7\%$ for $T=1.7$.

Figure 9

Order parameters for ${\mathbb{Z}}_6,{\mathbb{Z}}_3$, and ${\mathbb{Z}}_2$ symmetry breaking for $128\times 128$ systems in the (a) in-plane and (b) out-of-plane phases. All three symmetry subgroups appear to break together. The angle brackets denote a spatial as well as thermal average. The differences between these curves and the corresponding curves for the smaller systems are smaller than the plot markers.

Figure 10

Heat capacity and order parameter collapses for the Ising critical exponents $\beta =1/8,\nu =1$ via (14) and (13) for the clean system in the [(a) and (b)] in-plane and [(c) and (d)] out-of-plane phases.

Figure 11

Heat capacity and order parameter collapses for the ${\mathbb{Z}}_3$ universality class critical exponents $\alpha =1/3,\beta =1/9,\nu =5/6$ via (12) and (13) for the clean system in the [(a) and (b)] in-plane and [(c) and (d)] out-of-plane phases.

Figure 12

Disorder-averaged heat capacity per site for (a) the in-plane-phase system at disorder strength $\mathrm{\Delta }=0.3$ and (b) the out-of-plane-phase system at disorder strength $\mathrm{\Delta }=0.2$.

Figure 13

Snapshots of part of the in-plane-phase system at temperature $T=1.0$ in the case of (a) no disorder and (b) strong disorder $\mathrm{\Delta }=0.8$. (a) displays visible long-range order with stripes running parallel to the ${\mathit{a}}_3$ direction. (b) shows the intersection of three domains whose stripes run parallel to the ${\mathit{a}}_1$ (right), ${\mathit{a}}_2$ (bottom left), and ${\mathit{a}}_3$ directions (top left), with the domain walls indicated by red dashed lines.