Two-step melting of three-sublattice order in $S=1$ easy-axis triangular lattice antiferromagnets

We consider $S=1$ triangular lattice Heisenberg antiferromagnets with a strong single-ion anisotropy $D$ that dominates over the nearest-neighbour antiferromagnetic exchange $J$. In this limit of small $J/D$, we study low temperature ($T \sim J \ll D$) properties of such magnets by employing a low-energy description in terms of hard-core bosons with nearest neighbour repulsion $V \approx 4J + J^2/D$ and nearest neighbour unfrustrated hopping $t \approx J^2/2D$. Using a cluster Stochastic Series Expansion (SSE) algorithm to perform sign-problem-free quantum Monte Carlo (QMC) simulations of this effective model, we establish that the ground-state three-sublattice order of the easy-axis spin-density $S^z(\vec{r})$ melts in zero field ($B=0$) in a {\em two-step} manner via an intermediate temperature phase characterized by power-law three-sublattice order with a temperature dependent exponent $\eta(T) \in [\frac{1}{9}, \frac{1}{4}]$. For $\eta(T) < \frac{2}{9}$ in this phase, we find that the uniform easy-axis susceptibility of an $L \times L$ sample diverges as $\chi_L \sim L^{2-9 \eta}$ at $B=0$, consistent with a recent prediction that the thermodynamic susceptibility to a uniform field $B$ along the easy axis diverges at small $B$ as $\chi_{\rm easy-axis}(B) \sim B^{-\frac{4-18\eta}{4-9\eta}}$ in this regime.