Kosterlitz-Thouless melting of magnetic order in the triangular quantum Ising material TmMgGaO4

Frustrated magnets hold the promise of material realizations of exotic phases of quantum matter, but direct comparisons of unbiased model calculations with experimental measurements remain very challenging. Here we design and implement a protocol of employing many-body computation methodologies for accurate model calculations—of both equilibrium and dynamical properties—for a frustrated rare-earth magnet TmMgGaO4 (TMGO), which explains the corresponding experimental findings. Our results confirm TMGO is an ideal realization of triangular-lattice Ising model with an intrinsic transverse field. The magnetic order of TMGO is predicted to melt through two successive Kosterlitz–Thouless (KT) phase transitions, with a floating KT phase in between. The dynamical spectra calculated suggest remnant images of a vanishing magnetic stripe order that represent vortex–antivortex pairs, resembling rotons in a superfluid helium film. TMGO therefore constitutes a rare quantum magnet for realizing KT physics, and we further propose experimental detection of its intriguing properties.

. The partial electron density relevant for magnetic exchange interactions include the contributions from 4f electrons of Tm 3+ and 2p electrons of O 2− , as seen in the density of states in Supplementary   Figs. 1(a,b). The Tm-O-Tm superexchange paths are visualized as electron clouds overlap in Supplementary Fig. 1(c), and the Tm 3+ ions form a triangular lattice [lower plot in Supplementary Fig. 1(c)]. The two dimensionality of the material TMGO is manifested, as the electron density at the interlayer regime is negligibly small compared to that in the Tm-O-Tm path. Supplementary Note 2: Spin ordering and phase diagram of the J 1 -J 2 TLI.
As mentioned in the main text, by increasing the ratio of J 2 /J 1 we can drive the system from the clock order phase to a stripe order phase. To verify this, firstly we calculate the static magnetic structure Supplementary Figure 2. Real-space correlation S z 0 · S z r and a three-sublattice order. Each solid circle represents the real-space S z 0 · S z r correlation related to the central site 0, the size of circle denotes the magnitude of correlation, and the red(blue) color for the positive(negative) sign. The clock order pattern with an enlarged unit cell can be clearly seen. The XTRG calculation is performed with parameters J 1 = 0.99 meV, J 2 = 0.05J 1 , and ∆ = 0.54J 1 at T 0.57 K, on the YC geometry also specified in the plot.
Supplementary Figure 3. Schematic phase diagram of the J 1 -J 2 TLI and static structure factors. As J 2 increases, there exists a quantum phase transition between the clock and stripe phases, taking place at J 2c /J 1 ∼ 0.1 (for ∆/J 1 = 0.54). In (a) we collect the M-point intensity and plot it vs. J 2 /J 1 at four different temperatures, from which we see clearly that S M in small J 2 regime is continuously connected to that in the relatively large J 2 regime, i.e., the stripy phase. The contour plots of S(q) are shown in panels (b-e), with various J 2 /J 1 values (0 to 0.2). In the calculations, we fix the parameter as ∆/J 1 = 0.54, J 1 = 0.99 meV, and compute S(q) at T = 2.24 K for (b-e). factor S(q) from S z 0 S z r correlations, and show the results in Supplementary Fig. 2. From the real-space correlations on YC6 geometry, a three-sublattice order can be clearly identified, which translates into the S(q) peak at the K point [see, e.g., Fig. 5(b) in the main text, or Supplementary Fig. 3], signifying the presence of clock order.

Supplementary
However, as J 2 /J 1 increases to 0.2, S(q) changes its peak to the M point [ Supplementary Fig. 3(e)], consistent with a two-sublattice stripe-order pattern. Therefore, there must be a quantum phase transition between the clock and stripe phases, probably of first order. As the stripe order (with structure factor peak at M point) is in close proximity to the clock phase (a small J 2 drives the phase transition), we relate the finite energy M-roton excitations in the dynamic spin spectra with instability towards the stripe order.

Supplementary Note 3: Classical spin orders and their static structure factors.
The classical spin orders of J 1 −J 2 TLI has been throughly investigated in Supplementary Ref. [3]. Here we replot the possible classical spin configurations in Supplementary Fig. 4, along with their computed static structure factor S(q). One can see clearly that only orders in Supplementary Figs. 4(c,d) have peaks at M point. Supplementary Fig. 4(d) corresponds to the stripy order, in agreement with our simulated data (with large J 2 ), while the one in (c) has a Γ peak that is absent in our results. Therefore, the only classical Ising configuration that corresponds to our structure factor data is Supplementary Fig. 4(d), i.e., stripy order.
The other way around, in the context of our TLI model, Bragg peak at M = (1/2, 1/2) corresponds to a π phase shift, i.e., antiferromagnetic correlation along primitive vectors a and b (see primitive vector in  Supplementary Fig. 4(c) 0 energy. Therefore, it is clear that the stripy configuration is energetically more favorable in large J 2 limit. Note that in a recent study of triangular lattice magnet AgNiO 2 , the existence of stripe order was observed, which is ascribed to the relatively large J 2 in the compound [4].

Supplementary Note 4: TLI parameter fittings.
The parameter fitting workflow is as follows: we scan the parameters (J 1 , J 2 , ∆) to fit the specific heat C m (T ), magnetic entropy S m (T ), as well as the susceptibility χ(T ) (at a small magnetic field h=1 kOe), and find the optimal parameters. Given that, we compute the magnetization curves (at different T ) and magnetic entropy S m at finite fields h, and compare directly to experimental data [5,6], so as to ensure that the parameter set is adequate and precise to model the material.
To be concrete, we show in Supplementary Fig. 5 part of simulation data in our scanning. In Supplementary Fig. 5(a), we start with J 2 = 0 and scan various ∆ values. It is found that C m curves are sensitive, in terms of the peak height as well as the overall shape, to different ∆ values. By tuning ∆ (while keeping J 2 = 0), we find ∆/J 1 = 0.4 and J 1 = 1.1 meV can produce results in agreement with experimental C m curves. However, with this set of parameters (as well as other ∆ values) we clearly miss the experimental susceptibility line, as plotted in Supplementary Fig. 5(b). It therefore suggests that a finite J 2 should be involved in the fittings.
After a thorough scanning in the parameter space, we find ∆/J 1 = 0.54 with J 2 /J 1 = 0.05 can well Remind that there we push the calculations to YC9 lattice with width W = 9, and the fittings are equally good, suggesting the robustness of fittings vs. system sizes. Beyond equilibrium properties, we have also computed dynamical properties ω(k) with same parameter set on a 36 × 36 lattice, which also show excellent agreement with experimental data (e.g., the overall dispersion line and gap values). These direct comparisons lead us unambiguously to the conclusion that the above parameters of TLI can describe the material TMGO precisely.
Supplementary Note 5: The Curie-Weiss Fitting of high-T susceptibility.
As a complementary of the thermodynamic fittings in the main text, in Supplementary Fig. 6 we compare the experimental and simulated high-temperature magnetic susceptibility. It is found that the XTRG data lie on top of two experimental curves, with fitted Θ W 19.3 K, in very good agreement with the estimates in experimental works (e.g., 18.9 -19.1 K as indicated in the plot). Besides, the fitted constant C To conclude, through large-scale simulations of both equilibrium and dynamical properties, we pinpoint the model parameters of TMGO as J 1 = 0.99 meV, ∆/J 1 = 0.54, J 2 /J 1 = 0.05, and g = 13.212, which can be used to fit virtually all available experimental data, including the magnetic specific heat, entropy, susceptibility (both high-and low-temperature parts), and dynamical spin spectrum, etc.

Supplementary Note 6: Specific heat curves under external fields.
As noted in Fig. 2(d)  Ref. [5], the specific heat curves have also been measured under various external magnetic fields. In Supplementary Fig. 7(a), we redrawn the experimental data, and compare them, side by side, to the simulated C m (h, T ) curves. It can be seen that in both Supplementary Figs. 7(a) and (b), a low-T shoulder gradually appears under small fields, e.g., at h = 5 kOe, and then prominent peaks show up at around h = 10-20 kOe.
By further increasing the magnetic fields, the peak height gets declined and at the same time its position moves towards higher temperatures for h ≥ 30 kOe. It is remarkable that these highly nontrivial and non-monotonic behaviors of C m (h, T ) curves can also be understood within the TLI model. As shown in Fig. 2(d) of the main text, there exists a quasi-plateau structure at about 1/3 M sat in the curve, which suggests that the system undergoes a transition from the clock order to an UUD spin state upon increasing fields to h = 10-20 kOe. The forming of the UUD structure releases entropy and gives rise to the prominent low-T peak in C m , despite that the peaks in C m between 1 and 2 K are less pronounced in the simulated results than experiments. As the field strength further increases, the UUD order becomes weaker and accordingly the C m peak moves back to lower T side with a decreasing height. Eventually, for h ≥ 30 kOe, the system gradually polarizes into a ferromagnetic spin configuration, and the hump in C m moves towards higher T as h enhances.

Supplementary Note 7: Saddle point in the triangular tight-binding model.
When restricted in a subspace of configurations with only one pair of spins flipped (while others remain in the classical stripy order, in the small ∆ limit), we consider a "tight-binding" model on the triangular lattice, where i labels the lattice site, and |i labels a state with a spin flipped at site i. δ = a, b, a − b denotes nearest neighboring sites (a, b are primitive vectors shown in Fig. 1 of the main text). Remember that a second-order process related to S x terms in the Hamiltonian can actually tunnel between |i and |j , so t ∼ ∆ 2 (in the small ∆ limit), given i and j constitute a pair of nearest neighboring sites.
For the sake of simplicity, we set 0 = 0, t = 1, and take Fourier transformation of Eq. (1). The resulting dispersion (k) = 2[cos(k · a) + cos(k · b) + cos(k · (a − b))] is plotted in Supplementary Fig. 8(a). By cutting along the Γ-M-Γ path, we observe a quadratic low-energy dispersion near the minimum as shown in Supplementary Fig. 8(b). On the other hand, as shown in Supplementary Fig. 8(c)  Supplementary Note 8: A brief introduction of path integral QMC for TLI.

The Hamiltonian of quantum TLI is
By inserting a complete set of S z i eigenstates between each pair of exponentials, i.e., we can then rewrite the partition function as, where l indices the time slice τ = l · ∆τ .
Next, we employ the Trotter-Suzuki decomposition, where H 0 = J 1 i,j S z i S z j + J 2 i,j S z i S z j and H 1 = −h i S x i . Therefore the corresponding matrix element reads S z l+1 |e −∆τ H 1 e −∆τ H 0 |S z l = Λ N e −∆τ J 1 i,j S z i,l S z j,l −∆τ J 2 i,j S z i,l S z j,l +γ i S z i,l S z i,l+1 , with γ = − 1 2 ln tanh(∆τ h), and Λ 2 = sinh(∆τ h) cosh(∆τ h). For a certain configuration {S z i,l }, the configurational weight is, Now, the 2D quantum problem becomes a (2+1)D classical Ising model, which can be solved by local or global update schemes, both adopted in our practical Monte Carlo samplings.
Within this framework, the physical observables can be evaluated as where {S i } p denotes the spin configurations in which the measurement is performed at time p of the Markov chain. Besides, we are also interested in the imaginary time spin-spin correlation function, which should be calculated in prior to the spin spectrum S(q, ω). The latter can be obtained after a stochastic analytical continuation, as detailed in Methods section of the main text.