# Two-step melting of three-sublattice order in easy-axis triangular lattice antiferromagnets

## Abstract

We consider triangular lattice Heisenberg antiferromagnets with a strong single-ion anisotropy that dominates over the nearest-neighbour antiferromagnetic exchange . In this limit of small , we study low temperature () properties of such magnets by employing a low-energy description in terms of hard-core bosons with nearest neighbour repulsion and nearest neighbour unfrustrated hopping . 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 melts in zero field () in a two-step manner via an intermediate temperature phase characterized by power-law three-sublattice order with a temperature dependent exponent . For in this phase, we find that the uniform easy-axis susceptibility of an sample diverges as at , consistent with a recent prediction that the thermodynamic susceptibility to a uniform field along the easy axis diverges at small as in this regime.

###### pacs:

75.10.Jm## I Introduction

Controlled and quantitatively accurate theoretical calculations of the low temperature properties of geometrically frustrated magnets represent a major challenge, in part because it is difficult to devise analytical methods that can reliably capture the competition between different low-energy states(1), and in part because most such systems are not amenable to large-scale computations using numerically exact and unbiased quantum Monte Carlo (QMC) methods.(2). The latter difficulty has its origins in the presence of minus signs or phase factors in the statistical weights of different configurations that contribute to the partition function when it is sampled by such methods.

Triangular lattice antiferromagnets with moments at each site, with Hamiltonian

(1) |

in which the single-ion easy-axis anisotropy is much larger than the nearest-neighbour antiferromagnetic exchange between the moments , provide an example where interesting physics is accessible to controlled theoretical treatments in spite of the strongly frustrated nature of antiferromagnetic interactions. This is because the low-energy, low-temperature physics of in the limit can be mapped (3) on to that of hard-core bosons on the triangular lattice with large nearest-neighbour repulsion and parametrically small, unfrustrated boson hopping amplitude . The low-energy effective Hamiltonian is thus given by:

Here, is the chemical potential that controls deviations of the density from the half-filling, is the boson number at site , and is the corresponding boson annihilation operator.

The crucial minus sign in front of the hopping term ( is positive) renders the boson hopping unfrustrated and ensures that matrix elements of the boson Hamiltonian can, by addition of a suitable constant to Eq. LABEL:eq:H_b, be made strictly negative in the boson occupation-number basis. This implies that one can rewrite the partition function as a sum over configurations with positive weight in this basis. This representation can then be used to perform large-scale numerically exact quantum Monte Carlo studies(2) of the boson Hamiltonian, Eq. LABEL:eq:H_b, on finite lattices, from which one may draw conclusions about the low temperature physics of .

In order to describe the low temperature physics of using this approach, one sets(3) , and , , , and . This places the corresponding boson problem in a regime in which the frustrated nearest-neighbour repulsion between hard-core bosons dominates over the unfrustrated kinetic energy of these bosons. In this limit of , the interplay of frustration and quantum fluctuations is known to result in a “supersolid” ground state for , in which superfluidity coexists with three-sublattice (“”) density-wave order at wavevector (7); (8); (9); (10); (11); (12). For the easy-axis magnet, this implies(3) an unusual co-existence of three-sublattice spin-density order of and spin nematic order (6); (5); (4) for the transverse spin components ().

Here, our interest is in the zero field () melting of these coexisting low-temperature orders, which translates, in boson language, to the melting of supersolid order at half-filling (). In this boson language, which we use here and henceforth to make contact with earlier computational work on ground state properties of and other closely related bosonic systems, it is clear that the superfluid order will be lost via a Kosterlitz-Thouless transition at temperature . In this article, we confirm this expectation using detailed QMC computations that map out this line of Kosterlitz-Thouless transitions in the phase diagram of . How does three-sublattice density-wave order melt at ? Two very different answers have been given earlier: In Ref. (13) Boninsegni and Prokofiev have suggested that the melting proceeds even at half-filling via a single continuous transition in the universality class of the two-dimensional three-state Potts model. Guided by the analogy(14); (15) with the two dimensional six-state clock model(16); (17); (18); (19), the present authors(9) have instead argued that the melting will proceed at half-filling via a two-step melting process characterized by an intermediate phase exhibiting power-law three-sublattice () order.

As was emphasized recently(20), both scenarios in fact represent entirely consistent and generic possibilities for the melting of three-sublattice order in such systems. Either possibility can arise without fine-tuning of system parameters when the microscopic Hamiltonian has particle-hole symmetry in boson language (which translates to the Ising symmetry that flips the sign of the easy-axis component of the spin-density). Indeed, in the general phase diagram of such systems (for instance, with additional further-neighbour interactions), these two possibilities represent generic behaviours that are either separated from each other either by an intervening region of first order transitions, or connected to each other via the multicritical point whose universal properties were studied in Ref. (20). Thus, the ordered state of could, in principle, display any of these possibilities upon heating.

In this article, we study this question using QMC simulations of the effective Hamiltonian , and establish that three-sublattice order melts on heating in a two-step manner via an intermediate phase characterized by power-law three-sublattice order with a temperature dependent exponent . For in this phase, we find that the compressibility of an sample diverges as at . This implies that the uniform susceptibility of an sample of the easy-axis antiferromagnet will diverge as at , consistent with the prediction made in Ref. (20). By standard finite-size scaling arguments(20), this in turn implies that the thermodynamic susceptibility to a uniform field along the easy axis diverges at small as in the power-law three-sublattice ordered phase of . As emphasized in Ref. (20), this provides an interesting thermodynamic signature of the intermediate-temperature power-law ordered phase in such magnets.

The remainder of this article is organized as follows. In Sec. II, we describe the actual computational strategy used, as well as define the quantities that are measured in our QMC simulations. In Sec. III, we describe the results of these simulations and analyze these using standard finite-size scaling ideas. Finally, in Sec. IV, we conclude with a brief discussion of the potential experimental relevance of our results, and discuss some follow-up calculations that appear to be natural extensions of the work presented here.

## Ii Model and methods

In our numerical work, we set and study at on triangular lattices with periodic boundary conditions, with a multiple of up to . Using a customized “cluster” implementation(9); (21); (22) of the Stochastic Series Expansion (SSE) (23); (24); (25) QMC method to compute equilibrium averages at nonzero temperature , we focus on the finite-temperature phase diagram (Fig. 1) above the supersolid state that obtains for (7); (8); (9); (10); (13). In our cluster(9); (21); (22) implementation, is written as a sum comprising terms living on elementary triangles of the triangular lattice. Here the index distinguishes between the diagonal (in the number basis) term on each triangle and the three different hopping terms on the three links of each triangle. This cluster decomposition of the Hamiltonian allows the algorithm to distinguish between minimally frustrated triangular plaquettes (with total particle number or but not or ) and triangular plaquettes that pay a high interaction energy cost by having or particles on them. This turns out to be the key to obtaining reliable results at large values of , as was demonstrated in earlier work(9) on ground-state properties of the same model.

We characterize superfluidity of the low-temperature state by computing the superfluid stiffness . In order to characterize the three-sublattice density-wave order, we focus on , the component of the density operator at wavevector

and measure the equal-time correlation function as well as the corresponding fourth moment . In addition, we measure the static structure factor , where . We also measure the second moment of the SSE estimator ,whose average over the QMC run gives . Although is a basis-dependent quantity, is expected to scale in the same way as , whose computation involves technical difficulties that we side-step by focusing on . Finally, we also measure the compressibility , where .

## Iii Analysis of results

We expect superfluidity to be lost on heating via a Kosterlitz-Thouless transition at (16); (26); (27). In the vicinity of for fixed , we expect a finite-size scaling form (28)

(3) |

where and are fitting parameters which depends on temperature. At , we expect (16); (27); (28). Therefore, to pinpoint the critical point at each , we fit finite-size data for to this form and identify with the temperature at which the best-fit gives (29) (see Fig. 2). Using this procedure, we map out the superfluid transition shown in Fig 1.

If three-sublattice density-wave order is lost via a single continuous thermal phase transition at , one expects the Binder cumulants and to tend to deep in the ordered phase, and rise monotonically to a value of deep in the disordered phase. Additionally, in the vicinity of the continuous transition, one expects them to collapse on to scaling functions where , and is the density-wave correlation length exponent. Such a continuous transition can therefore be located from either Binder cumulant by looking for a crossing point of -dependent curves corresponding to various sizes . On the other hand, if long-range density-wave order is lost upon heating via a two-step melting process, with a power-law ordered intermediate phase for , one expects and to saturate at large in this power-law ordered phase to -independent limiting values that define the universal functions for . In other words, if there is an intermediate power-law ordered phase, we expect the Binder cumulant curves to stick over a finite range of rather than cross at one temperature.

With this in mind, we compare our QMC data for and with these competing predictions for their and dependence. We find clear evidence at all that the Binder cumulants stick rather than cross; a representative example is shown in Fig 3. When the Binder cumulants stick, we also find that , with increasing with temperature. An illustrative example is shown in Fig. 4. On the basis of this evidence, we conclude that three-sublattice density-wave order is lost via a two-step melting process with an intermediate temperature phase exhibiting power-law density-wave order.

It is not entirely straightforward to obtain the precise location of the upper and lower phase boundaries of this power-law density-wave ordered phase directly from the Binder cumulant data (for instance, the example displayed in Fig 3). Therefore, we exploit the fact that the long-distance correlations of in such a power-law ordered phase are expected(9); (14); (15); (20) to correspond to the long-distance properties of the order-parameter field in the analogous power-law ordered phase of the six-state clock model(16); (17); (18); (19), In such six-state clock models, the long-distance correlations of are controlled, in renormalization group (RG) language, by a line of fixed points (16), characterized by effective free-energy with , where , and is the phase of complex order parameter field that is subjected to a six-fold anisotropy. This range of corresponds to values of for which this fixed-line is stable to vortex excitations as well as six-fold anisotropy(16). For in this range, is characterized by power-law correlators with . This implies that the power-law exponent in such a power-law ordered phase lies in the range (16).

Reasoning by this analogy, we therefore expect the static structure factor of an system to obey , with in the phase with power-law three-sublattice order. With this in mind, we determine the upper (lower) transition temperature () at each by plotting () against for various and looking for a crossing point; a representative example is shown in Fig. 5. Using this procedure, we are able to map out the dependence of the phase boundaries and (Fig 1), and conclusively establish a persistent intermediate phase with power-law three-sublattice density-wave order for all studied here. For close to , power-law density-wave order co-exists with superfluidity, while at larger values of , this power-law ordered phase lies entirely above the superfluid transition . Our data gives no indication that the nature of the superfluid transition at is affected by the presence of power-law density-wave order, nor does it give any indication that the superfluidity affects the nature of the power-law density-wave ordered phase or transitions out of it. This may be understood by noting that the leading biquadratic term that couples the square of the modulus of the superfluid order parameter to the square of the modulus of the density-wave order parameter (in a Landau theory description of the transitions) represents an irrelevant perturbation of the decoupled critical points.

As is well-known, standard RG analysis (16) also predicts that for in the six-state clock model(18), with where and is a dimensionless constant. For the static structure factor in a system just above , we therefore expect(18) the analogous finite-size scaling form where is a universal function and is a constant. If we use the value of deduced from our earlier analysis, this scaling form has a single adjustible parameter at each , and the corresponding scaling collapse of data represents a stringent test of scaling; a representative example of such a scaling collapse is shown in Fig. 6. All of the foregoing provides clear evidence for the existence of an intermediate phase with power-law three-sublattice order at , with properties consistent with RG predictions and finite-size scaling (16); (28); (18).

The analysis of Ref. (20) predicts that , the slowly-varying (uniform) component of the density, will also display power-law correlations in the phase with power-law three-sublattice order: . For , these correlations decay slowly enough that they are predicted(20) to lead to a divergent contribution to the finite-size compressibility of an system. To check this prediction(20), we fit our finite-size compressibility data in the power-law ordered phase to the following form:

(4) |

for . In these fits, the value of is first obtained from a separate fit of the dependence of to a power-law form at each temperature , and the same value of is used in fitting to Eq. 4. An example of such fits is shown in Fig. 7, with the corresponding values of and displayed in Fig. 8. Based on this evidence, we conclude that our numerical results are completely consistent with the prediction made in Ref. (20).

As noted in Ref. (20), this result for a finite-size system at is equivalent, via finite-size scaling, to the statement that the thermodynamic compressibility at small has the singular form: . Thus, the numerical results displayed here for the finite-size scaling of , and the mapping between and discussed in the introduction, together establish that the power-law three-sublattice ordered intermediate temperature phase of is characterized by a divergent susceptibility to a uniform easy-axis field : for small . As noted in Ref. (20), this constitutes a striking thermodynamic signature of two-step melting of three-sublattice order in , which may be of potential experimental relevance.

## Iv Conclusion and Outlook

Our results thus conclusively establish that long-range three-sublattice order, present at low temperature in the phase diagram of when , melts in a two-step manner, via an intermediate phase [for ] with power-law three-sublattice order. On the other hand, superfluidity, present at low temperature in the phase diagram of for all , is lost via a transition at in the Kosterlitz-Thouless universality class. For very large values of , i.e. for , emerges as the smallest of these three transition temperatures. Conversely, for , is larger than both and , since these two transitions go rapidly to zero as approaches the threshold value of at which the ground state develops three-sublattice order. Finally, for , the Kosterlitz-Thouless transition occurs within the phase with power-law three-sublattice density-wave order.

Upon heating, frustrated magnets described by with large (which map to with ) are thus expected to first undergo a Kosterlitz-Thouless transition from a low-temperature state with co-existing transverse spin-nematic order () and longitudinal three-sublattice order () to a state with just longitudinal three-sublattice order of the easy-axis spin density. A second phase transition is expected to give rise to an intermediate state with power-law three-sublattice order of the easy-axis spin density. Finally, the system is expected to reaches a paramagnetic high-temperature state after a third transition involving loss of power-law three-sublattice order. In the intermediate phase with power-law three-sublattice order, the susceptibility to a uniform field along the easy-axis is predicted to diverge as in the limit, providing a thermodynamic signature of two-step melting of three-sublattice order in such magnets. Thus, this class of frustrated magnets is expected to exhibit several interesting physical phenomena, and it would be interesting to identify candidate materials whose magnetic properties are well-described by .

This divergent uniform susceptibility to a field along the easy-axis is perhaps easier to rationalize when the low-temperature three-sublattice ordered state is of the ferrimagnetic kind, i.e. characterized by a net easy-axis moment that accompanies spatial symmetry breaking. This is indeed the case in the low-temperature state of , as is clear upon noting that the spontaneous deviation from half-filing identified in the low temperature state of in Ref. (9) maps to a net easy-axis moment in the low-temperature state of . However, the arguments of Ref. (20) are independent of the nature of three-sublattice ordering, and predict that would also diverge in the power-law three-sublattice ordered phase associated with the two-step melting of antiferromagnetic three-sublattice order. A natural and interesting follow-up to our work would therefore be to test this stronger claim in other examples of easy-axis magnets which exhibit two-step melting of antiferromagnetic three-sublattice order. One such example is the transverse field Ising antiferromagnet on the triangular lattice,(30) in which this physics could perhaps be studied using a recently developed quantum cluster algorithm.(31)

## V Acknowledgements

We would like to thank L. Balents, and M. Boninsegni and N. Prokof’ev, for useful correspondence regarding the results of Ref. (8) and Ref. (13) respectively. We acknowledge use of computational resources of the Department of Theoretical Physics of the TIFR (KD), the Universtiy of Toronto (DH), and SISSA (DH), as well as computational resources funded by grant DST-SR/S2/RJN-25/2006 (KD). DH was the recipient of a fellowship for foreign scholars at the TIFR and a postdoctoral fellowship at SISSA (Italy) during the preliminary stages of this work. Another part of this work was completed while KD was a participant in the KITP Program on Frustrated Magnetism and Quantum Spin Liquids, funded in part by National Science Foundation grant NSF PHY11-25915. One of us would like to acknowledge the hospitality of ICTP (Trieste) and ICTS-TIFR (Bengaluru) during the writing of this manuscript.

### References

- R. Moessner and A. P. Ramirez, Phys. Today 59, 24 (2006).
- R. G. Melko in Strongly Correlated Systems, 185-206, Springer (Berlin 2013).
- K. Damle and T. Senthil, Phys. Rev. Lett. 97, 067202 (2006).
- P. Chandra and P. Coleman, Phys. Rev. Lett. 66, 100 (1991).
- A. F. Andreev and I. A. Grishchuk, Sov. Phys. JETP 60 (2), 267 (1984).
- H. H. Chen, and P. M. Levy, Phys. Rev. Lett. 27, 1383 (1971).
- G. Murthy, D. Arovas, and A. Auerbach, Phys. Rev. B 55, 3104 (1997).
- R. G. Melko, A. Paramekanti, A. A. Burkov, A. Vishwanath, D. N. Sheng, and L. Balents, Phys. Rev. Lett. 95, 127207 (2005).
- D. Heidarian and K. Damle, Phys. Rev. Lett. 95, 127206 (2005).
- S. Wessel and M. Troyer, Phys. Rev. Lett. 95, 127205 (2005).
- A. Sen, P. Dutt, K. Damle, and R. Moessner, Phys. Rev. Lett. 100, 147204 (2008).
- D. Heidarian and A. Paramekanti, Phys. Rev. Lett. 104, 015301 (2010).
- M. Boninsegni and N. Prokof’ev, Phys. Rev. Lett. 95, 237204 (2005).
- E. Domany, M. Schick, J. S. Walker, and R. B. Griffiths, Phys. Rev. B 18, 2209 (1978).
- E. Domany and M. Schick, Phys. Rev. B 20, 3828 (1979).
- J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson , Phys. Rev. B 16, 1217 (1977).
- J. Tobochnik, Phys. Rev. B 26, 6201 (1982).
- M. S. S. Challa and D. P. Landau, Phys. Rev. B 33, 437 (1986).
- E. Rastelli, S. Regina, and A. Tassi, Phys. Rev. B 69, 174407 (2004).
- K. Damle, Phys. Rev. Lett. 115, 127204 (2015).
- K. Damle, unpublished.
- K. Louis and C. Gros, Phys. Rev. B 70, 100410(R) (2004).
- O. F. Syljuasen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).
- A. W. Sandvik, Phys. Rev. B 59, R14157 (1999).
- A. W. Sandvik, J. Phys. A: Math. Gen. 25, 3667 (1992).
- J M Kosterlitz and D J Thouless, J. Phys. C: Solid State Phys. 5, 124 (1972); J. Phys. C: Solid State Phys. 6, 1181 (1973).
- D. R. Nelson and J. M. Kosterlitz, Phys. Rev. Lett. 39, 1201 (1977).
- H. Weber and P. Minnhagen, Phys. Rev. B 37, 5986 (1988).
- K. Harada and N. Kawashima, Phys. Rev. B 55, R11949 (1997).
- S. V. Isakov and R. Moessner, Phys. Rev. B 68, 104409 (2003).
- S. Biswas, G. Rakala, and K. Damle, arXiv:1512.00931 (unpublished).