Kosterlitz-Thouless Melting of Magnetic Order in the Triangular Quantum Ising Material TmMgGaO$_4$

TL;DR

This study uses advanced many-body computational methods to accurately model the frustrated magnet TmMgGaO$_4$, revealing it as an ideal realization of the triangular-lattice Ising model with a unique two-stage Kosterlitz-Thouless melting process.

Contribution

It introduces a novel computational protocol for modeling frustrated quantum magnets and demonstrates the realization of KT physics in TmMgGaO$_4$, including detailed dynamical spectra predictions.

Findings

TMGO is an ideal realization of the triangular-lattice Ising model with transverse field.

Magnetic order in TMGO melts via two successive KT phase transitions.

Dynamical spectra show vortex-antivortex pairs resembling rotons.

Abstract

Frustrated magnets host the promises of material realizations of new paradigm of quantum matter, while direct comparison of unbiased model calculations with experimental measurements is still very challenging. Here, we design and implement a protocol of employing many-body computation methodologies for accurate model calculation -- both equilibrium and dynamical properties -- of a frustrated rare-earth magnet TmMgGaO$_4$ (TMGO), which perfectly 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 detections of its intriguing properties.

Figures

Equations

$$H_{\mathrm{TLI}}=J_1\sum _{⟨i,j⟩}{S}_i^z{S}_j^z+J_2\sum _{⟨⟨i,j⟩⟩}{S}_i^z{S}_j^z-\sum _i(\Delta S_i^x+h{g}_{\parallel }{\mu }_B{S}_i^z),$$

$$\psi =m_{\mathrm{A}}^z+{\mathrm{e}}^{i2\pi /3}{m}_{\mathrm{B}}^z+{\mathrm{e}}^{i4\pi /3}{m}_{\mathrm{C}}^z,$$

$$\hat{\rho }(\Delta \tau )=e^{-\Delta \tau \hat{H}}=\sum _{n=0}^{N_c}\frac{(-\Delta \tau {)}^n}{n!}{\hat{H}}^n.$$

$$G(\mathbf{q},\tau )={\int }_0^{\infty }\frac{\mathbf{d}\omega }{\pi }\frac{{\mathbf{e}}^{-\tau \omega }+{\mathbf{e}}^{-(\beta -\tau )\omega }}{1+{\mathbf{e}}^{-\beta \omega }}\mathbf{B}(\mathbf{q},\omega )$$

$$H={ϵ}_0+\sum _i\sum _{\delta }t(c_i^{†}{c}_{i+\delta }+h.c.)$$

$$\begin{array}{rl}Z& =\text{Tr}\left[e^{-\beta \mathcal{H}}\right]\\ & =\text{Tr}\left[e^{-\Delta \tau \mathcal{H}}{e}^{-\Delta \tau \mathcal{H}}\cdots e^{-\Delta \tau \mathcal{H}}{e}^{-\Delta \tau \mathcal{H}}\right].\end{array}$$

$$\begin{array}{rl}1& =\prod _{i=1}^N\left[\sum _{S_i^z=\pm 1}|S_i^z⟩⟨S_i^z|\right]\\ & \equiv \sum _{\{S_i^z\}}|S^z⟩⟨S^z|,\end{array}$$

$$\begin{array}{rl}Z=\sum _{\{S_{i,l}=\pm 1\}}& ⟨S_1^z|e^{-\Delta \tau \mathcal{H}}|S_L^z⟩⟨S_L^z|e^{-\Delta \tau \mathcal{H}}|S_{L-1}^z⟩⟨S_{L-1}^z|e^{-\Delta \tau \mathcal{H}}|S_{L-2}^z⟩\cdots \\ & \cdots ⟨S_3^z|e^{-\Delta \tau \mathcal{H}}|S_2^z⟩⟨S_2^z|e^{-\Delta \tau \mathcal{H}}|S_1^z⟩,\end{array}$$

$$\begin{array}{rl}⟨S_{l+1}^z|e^{-\Delta \tau \mathcal{H}}|S_l^z⟩& =⟨S_{l+1}^z|e^{-\Delta \tau {\mathcal{H}}_1-\Delta \tau {\mathcal{H}}_0}|S_l^z⟩\\ & =⟨S_{l+1}^z|e^{-\Delta \tau {\mathcal{H}}_1}{e}^{-\Delta \tau {\mathcal{H}}_0}|S_l^z⟩+\mathbf{O}[(\Delta \tau {)}^2],\end{array}$$

$$⟨S_{l+1}^z|e^{-\Delta \tau {\mathcal{H}}_1}{e}^{-\Delta \tau {\mathcal{H}}_0}|S_l^z⟩={\Lambda }^N{e}^{-\Delta \tau J_1{\sum }_{⟨i,j⟩}{S}_{i,l}^z{S}_{j,l}^z-\Delta \tau J_2{\sum }_{⟨⟨i,j⟩⟩}{S}_{i,l}^z{S}_{j,l}^z+\gamma {\sum }_i{S}_{i,l}^z{S}_{i,l+1}^z},$$

$$\omega \{S_{i,l}^z\}=\left(\prod _l\prod _{⟨i,j⟩}{e}^{-\Delta \tau J_1{S}_{i,l}^z{S}_{j,l}^z}\right)\left(\prod _l\prod _{⟨⟨i,j⟩⟩}{e}^{-\Delta \tau J_2{S}_{i,l}^z{S}_{j,l}^z}\right)\left(\prod _{\tau }\prod _{⟨l,l^{\prime }⟩}\Lambda e^{\gamma S_{i,l}^z{S}_{i,l^{\prime }}^z}\right).$$

$$⟨\hat{O}⟩\approx \frac{1}{N}\sum _p^N\hat{O}(\{S_i{\}}_p)$$

$$G(\mathbf{q},\tau )=\frac{1}{N}\sum _{i,j}{e}^{i\mathbf{q}\cdot {\mathbf{r}}_{ij}}⟨S_i^z(\tau )S_j^z(0)⟩,$$

Full text

Kosterlitz-Thouless (KT) physics bestows interesting mechanism of phase transition upon two-dimension (2D) interacting system with a continuous symmetry. Although such symmetry is not allowed to break spontaneously at any finite temperature Mermin and Wagner (1966), phase transition can still take place from the high temperature disordered phase to a KT phase with quasi-long-range order, which has a topological root in the binding of the vortex and antivortex pair Kosterlitz and Thouless (1973); Kosterlitz (1974). Experimentally, the KT transition has been observed in thin helium films Bishop and Reppy (1978) and ultracold 2D bose gases Hadzibabic et al. (2006); Desbuquois et al. (2012). Two distinct types of elementary excitations, i.e., phonons and rotons, play essential roles in the related superfluid phenomena Landau (1949); Feynman (1955), and they are important for the understanding of liquid helium thermodynamics Bendt et al. (1959). Besides interacting bosons in liquid and gas, there are also theoretical proposals of KT transitions in solid-state magnetic systems such as the 2D classical XY  Kosterlitz and Thouless (1973); Kosterlitz (1974) and the frustrated quantum Ising models Isakov and Moessner (2003). However, to date, the material realization of the KT transition in 2D magnets has rarely been reported.

In the mean time, the search of exotic quantum magnetic states in the triangular lattice spin systems – the motif of frustrated magnets – has attracted great attention over the decades. Experimentally, the triangular lattice quantum magnets have been synthesized only very recently, including compounds Ba3CoSb2O9 Zhou et al. (2012); Susuki et al. (2013); Ma et al. (2016); Ito et al. (2017), Ba8CoNb6O24 Rawl et al. (2017); Cui et al. (2018), and a rare earth oxide YbMgGaO4 – which has been suggested as a quantum spin liquid candidate Li et al. (2015, 2016); Shen et al. (2016); Li et al. (2019a), while an alternative scenario of glassy and disordered induced state has also been proposed recently Zhu et al. (2017); Kimchi et al. (2018); Ma et al. (2018). On the other hand, an Ising-type triangular antiferromagnet TmMgGaO4 (TMGO, with Yb3+ replaced by another rare earth ion Tm3+) Cevallos et al. (2018); Li et al. (2020); Shen et al. (2019), as shown in Fig. 1(a) and explained in details in this paper, is the successful material realization of a quantum magnet with strong Ising anisotropy.

In this work, we construct the microscopic model of TMGO and employ two state-of-the-art quantum many-body simulation approaches: the exponential tensor renormalization group (XTRG) Chen et al. (2018) and quantum Monte Carlo equipped with stochastic analytic continuation (QMC-SAC) Sandvik (2016); Qin et al. (2017); Shao et al. (2017); Huang et al. (2018); Sun et al. (2018), to calculate both the thermodynamic and dynamic properties. By scanning various parameters and fit our simulation results to the existing experimental data Cevallos et al. (2018); Li et al. (2020); Shen et al. (2019), we find TMGO realizes a triangular-lattice transverse field Ising model and determine accurately its model parameters. Based on this, we conclude that TMGO should host the celebrated KT phase, and further predict several prominent features to be observed in TMGO, inspired by the experimental measurements for detecting KT physics in a superfluid thin film Bishop and Reppy (1978). It is worthwhile to point out that our calculation of quantum fluctuations goes beyond the linear spin-wave approximation in Ref. Shen et al. (2019) and puts the system in clock ordered (later melted through KT transitions), rather than disordered, regime. Therefore, our methods and results do not only explain the experimental findings, but more importantly establish a protocol of acquiring equilibrium and dynamic experiments of strongly correlated quantum material, such as TMGO, in an unbiased manner.

Results

Microscopic spin model. Due to strong spin-orbit coupling and crystal electric field splitting, TMGO can be described as an effective spin-1/2 model with strong easy-axis anisotropy, i.e., a triangular lattice Ising model (TLI), as shown in Fig. 1(a). First-principle calculations of the TMGO material are performed based on the density functional theory (DFT) Kresse and Furthmüller (1996); Kresse and Joubert (1999), where we see a large easy-axis anisotropy of room-temperature energy scale. By comparing the DFT energies of antiferromagnetic spin configuration to a ferromagnetic one, one finds the former has a lower energy, and the coupling strength can be estimated as a few tenths of meV. Density distribution of 4f electrons in TMGO can also be obtained, where it is observed that the 4f electrons of Tm3+ are coupled via superexchange mediated by 2p electrons of O2- within the triangular lattice plane (Supplementary Note 1).

In Ref. Li et al., 2020, the authors took the lowest two levels $|\Phi_{\pm}\rangle$ of Tm3+ as non-Kramers doublet, and construct a classical TLI with both nearest-neighbor (NN) and next-nearest-neighbor (NNN) interactions to account for the absence of zero point entropy observed in experiments. Substantial randomness was also introduced to explain the smooth magnetization curves even at very low temperature. Later on, inelastic neutron scattering (INS) results of TMGO reveal a clear magnon band Shen et al. (2019), suggesting the influence of coupling randomness should be modest in TMGO, and an adequate modelling of the material shall include in the Hamiltonian non-commuting terms with quantum fluctuations. Since the Kramers theorem is absent in Tm3+ system with total angular momentum $J=6$, a small level splitting $\Delta$ between the quasi-doublet $|\Phi_{\pm}\rangle$ is involved, as shown in Fig. 1(a). Therefore, a quantum TLI model was proposed Shen et al. (2019), with spin-1/2 Hamiltonian

[TABLE]

where $\langle,\rangle$ \big{(}\langle\langle,\rangle\rangle\big{)} stands for NN (NNN) couplings $J_{1}$ ($J_{2}$), $\Delta$ the energy splitting between $|\Phi_{\pm}\rangle$ (i.e., the intrinsic transverse field), and $h$ is the external magnetic field. $g_{\parallel}=2Jg_{J}$ constitutes the effective spin-1/2 $g$ factor, with $g_{J}$ the Landé factor.

The phase diagram of quantum TLI has been studied intensively with analytic and numeric methods in the past Isakov and Moessner (2003); Wang et al. (2017); Powalski et al. (2013); Fey et al. (2019) and is schematically shown in Fig. 1(b). We indicate the TMGO model parameter with the vertical arrow (the determination of parameters is given below). From high to low temperatures, the system first goes through a KT transition at $T_{h}$ from the paramagnetic phase to a KT phase with power-law (algebraic) spin correlations. At a lower temperature $T_{l}$, the system enters the clock phase with a true long-range order depicted in Fig. 1(a). This three-sublattice clock order breaks the discrete lattice point group as well as the $Z_{2}$ spin symmetries, giving rise to a low but finite transition temperature $T_{l}$.

Increasing the next-nearest neighbor couplings $J_{2}$, say, at $J_{2}/J_{1}=0.2$, we find the static magnetic structure factor develops a stripe order Metcalf (1974) with structure factor peak at M point of the Brillouin zone (BZ) (see Supplementary Notes 2,3). This magnetic stripe order, as shown in Fig. 1(c), has been observed previously in TLI material AgNiO2, where $J_{2}$ coupling is relatively strong ($\sim 0.15J_{1}$, along with other interactions) and the exotic KT physics is absent there Wheeler et al. (2009). In TMGO, however, the clock order wins over the stripe order as $J_{2}/J_{1}\simeq 0.05$ is relatively small in this material. Nevertheless, as will be shown below, a ghost of the latter order – the M rotonlike modes – remains in the spin spectrum Shen et al. (2019), which turns out to be related to a vortex-antivortex pair excitation in the topological language (see Discussion section).

Thermodynamics and parameter fittings. The model parameters in Eq. (Results) can be accurately determined through fitting the available experimental data of TMGO Cevallos et al. (2018); Li et al. (2020); Shen et al. (2019), from which we find $J_{1}=0.99$ meV, $J_{2}=0.05J_{1}$, $\Delta=0.54J_{1}$ and $g_{\parallel}=13.212$. We present in Fig. 2 the calculated thermodynamic quantities and their experimental counterparts Li et al. (2020); Shen et al. (2019), where excellent agreements are seen. In Fig. 2(a), at high temperatures $T>30$ K, the magnetic entropy $S_{\rm{m}}$ approaches $R\ln{2}$, corresponding to the paramagnetic phase with effective spin-1/2. As temperature decreases, $S_{\rm{m}}$ gradually releases throughout the intermediate-temperature regime and approaches zero below $T_{l}^{*}\simeq 1$ K, as the long-range clock order develops. In Fig. 2(b) the very good agreement in magnetic specific heat $C_{\rm{m}}$ extends from high $T$ ($\sim 30$ K) all the way down to low temperature $\sim$1 K.

In the fittings of the magnetic entropy and specific heat curves in Figs. 2(a,b), we rescale the $T$-axes (in the unit of $J_{1}$) to lay the model calculations on top of the experimental measurements, and in this way we find the optimized $J_{1}=$ 0.99 meV. In both plots, the $y$-axes scaling ratios are associated with the ideal gas constant $R=8.313$ J$\cdot$mol${}^{-1}\cdot$K*-1* and thus fixed. In Fig. 2(c), we fit the magnetic susceptibility $\chi(h)=\frac{M(h)}{h}$, with $M(h)$ the uniform magnetization (per Tm3+), under external fields of a small value $h=1$ kOe and a larger one $10$ kOe. As shown in Fig. 2(c), by setting the effective Landé factor $g_{\parallel}=13.212$, we fit both susceptibility curves very well. This completes the model parameters in the Hamiltonian Eq. (Results) (see more fitting details in Supplementary Note 4), and note this parameter set also leads to accurate entropy results at 5 kOe when put in a direct comparison to the experimental line in Fig. 2(a).

With the parameters $J_{1}$, $J_{2}$, $\Delta$, and $g_{\parallel}$ determined from above fittings, we can compute the magnetization curve $M(h)$ and compare it directly with the independently measured experimental curves in Fig. 2(d). Note there exists a turning point at about 1/3 magnetization, under a magnetic field around 20-25 kOe. It becomes clearer as $T$ decreases further down to 40 mK [see the inset of Fig. 2(d)], suggesting the existence of field-induced quantum phase transition in TMGO. This sharp change of behaviors can also be witnessed in the specific heat curves under various magnetic fields (see Supplementary Note 6), which can also be understood very well within the set of parameters obtained above.

Lastly, we briefly discuss the scalings in the uniform susceptibility $\chi$ in Fig. 2(c) as well as Fig. 3, where susceptibility curves under more external fields are computed theoretically in the later plot, as a complementary. These $\chi$ data reflects the two-step establishment of magnetic order as $T$ lowers. As pointed out in Refs. Damle, 2015; Biswas and Damle, 2018, a universal scaling $\chi(h)=h^{-\alpha}$ appears for small fields $h$ with $\alpha=\frac{4-18\eta(T)}{4-9\eta(T)}$ in the KT regime, where $\eta(T)\in[\frac{1}{9},\frac{2}{9}]$ is the anomalous dimension exponent of the emergent XY order parameter varying with temperature. For $T=T_{l}$, $\eta=\frac{1}{9}$ and $\chi(h)\sim h^{-2/3}$, which diverges as $h$ approaches zero; while for some higher temperature (slightly) below $T_{h}$, $\eta=\frac{2}{9}$ and $\chi(h)$ remains a constant vs. $h$. Therefore, at small external field, the increase of $\chi$ at intermediate $T$ reflects the decrease of $\eta(T)$ vs. $T$, and such enhancement becomes less prominent for a relatively larger field, say, $h=10$ kOe. This salient difference indeed can be noticed in the experimental as well as our numerical curves in Fig. 2(c). Moreover, in Fig. 3 we calculate more susceptibility curves under various magnetic fields $h$ between 6.5 kOe and 9 kOe, where the power-law scaling is shown explicitly and the anomalous exponent $\eta$ can be extracted therein.

Spin spectra and magnetic structure factors. Spin frustration can lead to strong renormalization effects, which in turn gives rise to interesting spectrum. Here we employ the QMC-SAC approach Sandvik (2016); Qin et al. (2017); Shao et al. (2017); Huang et al. (2018); Sun et al. (2018) to compute the spin spectra $S(\mathbf{q},\omega)$ from $S^{z}$ spin correlations, at various temperatures (see Methods). The obtained spectra, with model parameters determined from equilibrium data fittings, are plotted in Fig. 4 and compared directly to the INS results Shen et al. (2019). Figure 4(a) depicts the spin spectrum inside the clock phase, at a low temperature $T=0.5$ K. Since the clock phase is of discrete symmetry breaking, the S($$\mathbf{q}$$,\omega=0) signals the Bragg peak of the clock order and there exists a small gap $\sim 0.1$ meV between the $\omega=0$ and finite $\omega$ spectra, consistent with the INS result. The rotonlike modes are also clearly present in the QMC-SAC results, with an energy gap $\simeq 0.4$ meV, in quantitative agreement with that in Ref. Shen et al. (2019).

Figure 4(b) is the QMC-SAC spectrum calculated according to the parameters (${\Delta}/J_{1}\simeq 1.36$, $J_{2}/J_{1}\simeq 0.046$) given in Ref. Shen et al. (2019). As mentioned in the introduction, we find, via spin structure factor calculations, that such set of parameters actually put the model in the disordered paramagnetic phase with $\Delta>\Delta_{\rm{c}}\sim 0.8J_{1}$. It is possible that the fitting scheme adopted in Ref. Shen et al. (2019) is based on mean-field treatment and cannot capture the quantum fluctuations inherent to the quantum TLI model as well as the material TMGO. This is a clear sign that the unbiased quantum many-body calculation scheme in our work is the adequate approach to explain the experimental results.

We continue with the correct parameter set and rise the temperature to $T=1.45$ K in Fig. 4(c). It is interesting to see that the dispersion still resembles that in the clock phase of Fig. 4(a) but with a vanishing gap at the K point, as well as softened M roton modes. To show it more clearly, we plot the intensity at M in Fig. 4(d), where the roton gap gets reduced as $T$ increases, with substantially broadened linewidths. Since M rotonlike excitations can be related to vortex-pair excitation (see Discussions), this softening of M roton is consistent with the scenario of vortex proliferations near the upper KT transition. Such remarkable spectra constitute a nontrivial prediction to be confirmed in future INS experiments.

Besides, the static magnetic structure factor $S(\mathbf{q})=\sum_{i,j}$ e${}^{i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}\langle S_{i}^{z}\,S_{j}^{z}\rangle$ are also simulated, where $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$ run throughout the lattice. Fig. 5(a) shows the temperature dependence of $S_{\textup{K}}$ and $S_{\textup{M}}$, where one observes an enhancement of $S_{\textup{M}}$ at intermediate temperature, signifying its closeness towards the stripe order. At $T<T_{l}^{*}$, the enhancement of $S_{\textup{M}}$ vanishes and instead the $S_{\textup{K}}$ intensity becomes fully dominant. Figures 5 (b,c,d) show the $S(\mathbf{q})$ results at low ($T=0.57$ K), intermediate ($T=2.2$ K), and high temperaure ($T=4.5$ K). In the clock phase, $S(\mathbf{q})$ evidently peaks at the K point, the ordering wavevector of the three-sublattice clock phase, while in the intermediate temperature regime, notably there exists a “ghost” peak at the M point, manifesting the existence of short-range stripe order selected by thermal fluctuations. These interesting features are gone at higher $T=4.5$ K, where strong fluctuations considerably weaken the structure factor peaks.

Discussion

As postulated by Landau, the phonon-roton spectrum play an important role in understanding low-$T$ thermodynamics and critical velocities of superfluidity, etc, in liquid helium, where the rotons are believed to be related to local vorticity of the fluid Landau (1949). Roton constitutes a minima in the spectrum at finite momentum and energy, and has been regarded as a quantum analog of hydrodynamic vortex ring, as coined “the ghost of vanished vortex ring” by Onsager Donnelly (1974). On the other hand, as derived from a trial wavefunction introduced by Feynman Feynman (1955), the roton excitation has energy $\epsilon_{\mathbf{q}}\sim\hbar\mathbf{q}^{2}/[2m\,S(\mathbf{q})]$, with $m$ is the helium atom mass and $S(\mathbf{q})$ the structure factor peak, and it is therefore also associated with an incipient crystalization Bragg peak competing with superfluidity Nozières (2004). Rotonlike excitations are also found in thin helium films Carneiro et al. (1976) and frustrated triangular Heisenberg (TLH) magnets. In the latter case, M point rotonlike modes were predicted theoretically Starykh et al. (2006); Chernyshev and Zhitomirsky (2006); Zheng et al. (2006a, b) and confirmed in recent experiments Ma et al. (2016); Ito et al. (2017), whose nature is under ongoing investigations Verresen et al. (2019); Ferrari and Becca (2019). Notably, it has been proposed that the M rotons in TLH can be softened by further enhancing spin frustration (and thus quantum fluctuations) Ferrari and Becca (2019) or thermal fluctuations Chen et al. (2019), which melts the long-range or incipient semi-classical 120∘ order, driving the system into liquid-like spin states.

In the frustrated magnet TMGO, as mentioned earlier in Fig. 4, there exists rotonlike modes with large density of states, which becomes softened even at low temperatures, melting the clock order, and strongly influences thermodynamics of the system. Like in liquid helium, roton in TMGO has also a topological origin: we demonstrate below that the rotonlike modes represent bound states of topological vortex excitations, via a pseudo-spin mapping of the spin stripe order. As shown in Fig. 1(c), the stripe order constitutes a proximate competing order to the clock state. Although it eventually gets perished in TMGO at low $T$, the stripe order leaves a “ghost image” in the excitation spectrum, i.e., rotonlike dip along the $\Gamma$-M-$\Gamma$ path in Fig. 4(a). Correspondingly, there exists an incipient $S_{\textup{M}}$ peak in the static structure factor at intermediate $T$ [Fig. 5(c)], i.e., $\epsilon_{\mathbf{q}}\sim 1/S(\mathbf{q})$, similar to rotons in the superfluid helium discussed above Feynman (1955); Nozières (2004).

On top of the spin stripe order, the M rotonlike excitation can be related to a locally bounded vortex pair, some form of “rotational motion” happening in TMGO. We perform a pseudo-spin mapping

[TABLE]

where $\psi=|\psi|{\rm e}^{i\theta}$ is the complex order parameter Wang et al. (2017), i.e., the pseudo spin. As shown in Fig. 1(c), $\psi$ is located in the center of each triangle, with emergent XY degree of freedom $\theta$. In Eq. (1), $m_{\gamma}^{z}=\pm 1/2$ represents the spin-up(-down) of corresponding spin $S^{z}$ components at $\gamma$-sublattice ($\gamma={\rm A,B,C}$). This mapping helps establishing a Landau-Ginzburg theory of TLI Isakov and Moessner (2003); Wang et al. (2017), and the clock order shown in Fig. 1(a) corresponds to a ferromagnetic order of pseudo spins.

As shown in Fig. 1(c), we create a vortex pair by applying $S^{x}$ operator on two adjacent sites to simultaneously flip their spin orientations. We note that any closed loop enclosing only the vortex defect [red dot in Fig. 1(c)] leads to a winding number 1 (modulo $2\pi$), while those around the anti-vortex (black dot) to $-1$. Zero winding number can be counted when a pair of defects (or no defects at all) are enclosed by the loop. Moreover, one can further move the vortex on the triangular lattice, like in a “tight-binding” model, by flipping spins on further neighboring sites, which naturally leads to a quadratic-type low-energy dispersion near M point along $\Gamma$-M-$\Gamma$ (see Supplementary Note 7). Since the ghost peak in Fig. 5(c) only suggests a short-range stripe correlation, the vortex pair can thus move only within a small cluster with incipient stripe order, i.e., they are bounded. The vortex pair only unbinds at the upper KT transition $T_{h}$ where vortices are proliferated Kosterlitz and Thouless (1973); Kosterlitz (1974), as seen by the “softening” of M rotonlike mode in dynamical spectrum in Fig. 4.

To conclude, in this work, we have established a protocol of understanding and explaining experiments of frustrated magnets in an unbiased manner, with XTRG and QMC-SAC machinery. The thermodynamic and dynamic results of TMGO are captured to great accuracy, thus allowing comprehensive studies of KT physics therein. At intermediate temperature, the KT phase of TMGO realizes a magnetic analog of 2D superfluid phase with several intriguing properties: (i) There emerges a spin XY symmetry and correspondingly complex order parameter $\psi$, which bears quasi-long-range correlation and phase coherence; (ii) The finite-$T$ spin spectrum contains the long-wave-length magnon and competing gapped rotonlike modes near the BZ boundary with energy signifying the binding of vortex-antivortex pair, which plays a key role in determining finite-$T$ phase of the system; (iii) The quasi-long-range XY order melts and the TMGO becomes paramagnetic as $T$ is above $T_{h}$, driven by the proliferation of vortex excitations in analogy to the superfluid transition in a helium thin film Bishop and Reppy (1978); Minnhagen (1987). We note that, different from the liquid helium, TMGO has two-temperature scales that outline the intermediate-temperature KT phase. Similar separation of scales has been seen in other quantum magnetic materials, by temperature Banerjee et al. (2016); Chen et al. (2019) or spatial dimensions Lake et al. (2005), but the KT phase is the first time to be seen. Furthermore, the $T_{l}^{*},T_{h}^{*}$ in the present work constitute tentative estimated of two KT transition temperatures roughly from thermodynamics, which still needs to be precisely determined both numerically and experimentally in the future.

The extraction of anomalous exponents $\eta(T)$ constitutes another interesting future study. The exponent $\eta(T)$ of the KT phase appears in the algebraic correlations in $\langle S^{z}_{i}S^{z}_{j}\rangle\sim|r_{i}-r_{j}|^{-\eta(T)}$ and seems rather indirect to measure in solid-state experiements. Nevertheless, as discussed above in Fig. 3, $\eta(T)$ can be determined from the uniform susceptibility $\chi$, a routine magnetic quantity in experiments. Moreover, it would be interesting to check several distinct predictions of KT physics in this 2D magnetic material. One renowned phenomenon is the universal jump in superfluid density at the KT transition Nelson and Kosterlitz (1977), as observed experimentally in helium film Bishop and Reppy (1978). Through calculations of the $q$-clock model ($q=5,6$), people have revealed universal jumps in the spin stiffness at both upper and lower KT transitions Chatelain , which certainly is an interesting prediction to check in TMGO. Dynamically, nuclear magnetic resonance measurements of relaxation time 1/$T_{1}$ can be conducted, which probe signals of low-energy magnetic dissipations at the KT transition where vortices proliferate. Besides, non-equilibrium thermodynamics such as thermal transport, would also be very worthwhile to explore in the TMGO magnet.

Methods

Quantum manybody computations. In this work, we combined two many-body numerical approaches: quantum Monte Carlo (QMC) Wang et al. (2017) and the exponential tensor renormalization group (XTRG), the latter method is recently introduced based on matrix product operators (MPOs) and logarithmic temperature scales Chen et al. (2018). XTRG is employed to simulate the TLI down to temperatures $T<0.5$ K on YC $W\times L$ geometries up to width $W=9$ with various lengths up to $L=12$. Both dynamical and equilibrium properties are simulated, with the purpose of fitting the experimental data and obtaining the right parameters as well as to make some predictions. The QMC is performed in the space-time lattice of $L\times L\times L_{\tau}$, where $L=36$ and $L_{\tau}=\beta/\Delta\tau$ with $\Delta\tau=0.05$ and $\beta\equiv 1/T$. The space-time configuration is written in the $S^{z}_{i,\tau}$ basis with both local and Wolff-cluster updates to overcome the long autocorrelation time. Since the QMC method is standard, we will only introduce the stochastic analytic continuation (SAC) scheme below and leave the QMC itself to the Supplementary Note 8.

Exponential thermal tensor network method. For the calculations of equilibrium properties, we start from an high-$T$ density matrix $\hat{\rho}(\Delta\tau)=e^{-\Delta\tau H}$, whose MPO representation can be obtained conveniently and accurately (up to machine precision), at a small $\Delta\tau\sim 10^{-3\sim-4}$. One way to obtain such an accurate MPO representation is to exploit the series-expansion thermal tensor network technique Chen et al. (2017) via the expansion

[TABLE]

Given the $\hat{\rho}(\Delta\tau)$ representation, traditionally one evolves $\hat{\rho}(\beta)$ linear in $\beta$ to reach various lower temperatures, i.e., $\beta=L_{\tau}\Delta\tau$ increases by a small value $\Delta\tau$ after each step by multiplying $\hat{\rho}(\Delta\tau)$ to the density matrix Li et al. (2011). However, this linear scheme is not optimal in certain aspects, and encounters challenges in generalization to 2D. Instead, recent study shows that the block entanglement entropy of MPO is bound by $S_{E}\leq a\ln\beta+\mathrm{const.}$ at a conformal critical point, with $a$ an universal coefficient proportional to the central charge Chen et al. (2018). This suggests an exponential procedure of performing cooling procedure. Based on this idea, we have developed the XTRG method, which turns out to be highly efficient in simulating both 1D critical quantum chains and various 2D lattice systems Chen et al. (2018, 2019); Li et al. (2019b).

In XTRG, we cool down the system by multiplying the thermal state by itself, i.e., $\rho_{0}\equiv\rho(\Delta\tau)$, $\rho_{1}\equiv\rho_{0}\cdot\rho_{0}=\rho(2\Delta\tau)$, thus $\rho_{n}\equiv\rho_{n-1}\cdot\rho_{n-1}=\rho(2^{n}\Delta\tau)$, and reach the low-$T$ thermal states exponentially fast. Efficient compression of MPO bonds is then required to maintain the cooling procedure, where a truncation scheme optimizing the free energy, and in the mean time maintaining the thermal entanglement, is involved. One advantage of XTRG is the convenience and high efficiency to deal with long-range interactions after the quasi-1D mapping. For the TLI model with nearest- ($J_{1}$) and next-nearest-neighbor ($J_{2}$) interactions considered in this work, we map the 2D lattice into a quasi-1D geometry following a snake-like path. The Hamiltonian thus contains “long-range” interactions and has an efficient MPO representation with geometric bond dimension $D_{H}=2W+2$, with $W$ the width of the lattice. In XTRG calculations, the computational costs scale with power $O(D^{4})$, with $D$ the retained bond dimension in MPO, which is chosen as large as $500$-$600$ in the present study, assuring accurate thermodynamical results down to sub-Kelvin regime.

QMC-SAC approach. We exploit the path integral QMC Wang et al. (2017), equipped with SAC approach, to compute the dynamical properties. The time displaced correlated function, defined as $G(\tau)=\langle S^{z}(\tau)S^{z}(0)\rangle$, for a set of imaginary times $\tau_{i}\ (i=0,1,\cdots,L_{\tau})$ with statistical errors can be obtained from QMC simulations. By SAC method Sandvik (2016); Qin et al. (2017); Huang et al. (2018); Sun et al. (2018), the corresponding real-frequency spectral function $S(\omega)$ can be obtained via $S(\tau)=\int_{-\infty}^{\infty}d\omega S(\omega)K(\tau,\omega)$, where the kernel $K(\tau,\omega)$ depends on the type of the spectral function, i.e., fermionic or bosonic, finite or zero temperature. The spectra at positive and negative frequencies obey the relation of $S(-\omega)=e^{-\beta\omega}S(\omega)$ and we are restricted at the positive frequencies and the kernel can therefore be written as $K(\tau,\omega)=\frac{1}{\pi}(e^{-\tau\omega}+e^{-(\beta-\tau)\omega})$. In order to work with a spectral function that is itself normalized to unity, we further modify the kernel and the spectral function and arrive at the transformation between the imaginary time Green’s function $G(\bf{q},\tau)$ and real-frequency spectral function $B(\bf{q},\omega)$

[TABLE]

where $B(\bf{q},\omega)=S(\bf{q},\omega)(1+e^{-\beta\omega})$.

In the practical calculation, we parametrize the $B(\bf{q},\omega)$ with a large number of equal-amplitude $\delta$-functions sampled at locations in a frequency continuum as $B(\omega)=\sum_{i=0}^{N_{\omega}-1}a_{i}\delta(\omega-\omega_{i})$. Then the relationship between Green’s function obtained from Eq. (3) and from QMC can be described by the goodness of fit $\chi^{2}$, i.e. $\chi^{2}=\sum_{i=1}^{N_{\tau}}\sum_{j=1}^{N_{\tau}}(G_{i}-\bar{G}_{i})C_{ij}^{-1}(G_{j}-\bar{G}_{j})$, where $\bar{G}_{i}$ is the average of QMC measurement and $C_{ij}$ is covariance matrix $C_{ij}=\frac{1}{N_{B}(N_{B}-1)}\sum_{b=1}^{N_{B}}(G_{i}^{b}-\bar{G}_{i})(G_{j}^{b}-\bar{G}_{j})$, with $N_{B}$ the number of bins. Then we update the series of $\delta$-functions in a Metropolis process, from $(a_{i},\omega_{i})$ to $(a^{\prime}_{i},\omega^{\prime}_{i})$, to get a more probable configuration of $B(\bf{q},\omega)$. The weight for a given spectrum follows the Boltzmann distribution $P(B)\propto\exp(-\chi^{2}/2\Theta)$, with $\Theta$ a fictitious temperature chosen in an optimal way so as to give a statistically sound mean $\chi^{2}$ value, while still staying in the regime of significant fluctuations of the sampled spectra so that a smooth averaged spectral function is obtained. The resulting spectra will be collected as an ensemble average of the Metropolis process within the configurational space of $\{a_{i},\omega_{i}\}$, as detailed in Refs. Sandvik, 2016; Qin et al., 2017; Huang et al., 2018; Sun et al., 2018.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability

All numerical codes in this paper are available upon request to the authors.

Acknowledgements

The authors are indebted to Ziyu Chen, Wentao Jin, Andreas Weichselbaum, Yao Shen, and Yue-Sheng Li for stimulating discussions. This work was supported by the Ministry of Science and Technology of China through the National Key Research and Development Program (Grant No. 2016YFA0300502), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB28000000), the National Science Foundation of China (Grant No. 11421092, 11574359, 11674370, 11974036, 11874115, and 11834014) and Research Grants Council of Hong Kong Special Administrative Region of China through 17303019 and the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. B.-B.C. was supported by the German Research foundation, DFG WE4819/3-1. W.L. was supported by the Fundamental Research Funds for the Central Universities. We thank the Center for Quantum Simulation Sciences at Institute of Physics, Chinese Academy of Sciences, the Computational Initiative at the Faculty of Science at the University of Hong Kong and the Tianhe-1A platform at the National Supercomputer Center in Tianjin and Tianhe-2 platform at the National Supercomputer Center in Guangzhou for their technical support and generous allocation of CPU time.

Author contributions

H.L. and Y.D.L. contributed equally to this work. W.L., Y.Q. and Z.Y.M initiated the work. H.L. and B.B.C. performed the calculations of thermodynamics, X.T.Z. and X.L.S. performed first-principle calculations of the band structure, and Y.D.L. performed the dynamical property simulations. All authors contributed to the analysis of the results. Z.Y.M. and W.L. supervised the project.

Additional information

Supplementary Information is available in the online version of the paper.

Competing interests: The authors declare no competing interests.

I Supplementary Information for:

Kosterlitz-Thouless Melting of Magnetic Order in the Triangular Quantum Ising Material TmMgGaO4

Supplementary Note 1: The electron density distributions and magnetic interactions in TMGO.

Density-functional theory (DFT) calculations of TMGO can be performed via the Vienna $ab$ $initio$ simulation package, with the projector augmented wave method Kresse and Furthmüller (1996); Kresse and Joubert (1999). In the DFT calculations, we take the TMGO lattice parameters $a=b=3.4260$ Å  and $c=25.1690$ Å  as determined from experiments.

When Coulomb interaction is switched off ($U=0$), the band structure of TMGO indicates a metallic state [Supplementary Fig. S1(a)], while a finite $U$ ($=3$ eV) opens up a Mott gap as shown in Supplementary Fig. S1(b). The partial electron density relevant for magnetic exchange interactions include the contributions from 4f electrons of Tm3+ and 2p electrons of O2-, as seen in the density of states in Supplementary Figs. S1(a,b). The Tm-O-Tm superexchange paths are visualized as electron clouds overlap in Supplementary Fig. S1(c), and the Tm3+ ions form a triangular lattice [lower plot in Supplementary Fig. S1(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 factor $S(\mathbf{q})$ from $\langle S_{0}^{z}S_{\rm{r}}^{z}\rangle$ correlations, and show the results in Supplementary Fig. S2. From the real-space correlations on YC6 geometry, a three-sublattice order can be clearly identified, which translates into the $S(\mathbf{q})$ peak at the K point [see, e.g., Fig. 5(b) in the main text, or Supplementary Fig. S3], signifying the presence of clock order.

However, as $J_{2}/J_{1}$ increases to 0.2, $S(\mathbf{q})$ changes its peak to the M point [Supplementary Fig. S3(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. Metcalf (1974). Here we replot the possible classical spin configurations in Supplementary Fig. S4, along with their computed static structure factor $S(\mathbf{q})$. One can see clearly that only orders in Supplementary Figs. S4(c,d) have peaks at M point. Supplementary Fig. S4(d) corresponds to the stripy order, in agreement with our simulated data (with large $J_{2}$), while the one in (c) has a $\Gamma$ peak that is absent in our results. Therefore, the only classical Ising configuration that corresponds to our structure factor data is Supplementary Fig. S4(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 $\pi$ phase shift, i.e., antiferromagnetic correlation along primitive vectors $\mathbb{a}$ and $\mathbb{b}$ (see primitive vector in Fig. 1 of the main text), while the correlation is ferromagnetically along $\mathbb{a}+\mathbb{b}$ and $\mathbb{a}-\mathbb{b}$. This clearly corresponds to a magnetic stripy phase, in the context of TLI model.

Furthermore, we can also argue that large $J_{2}>0$ favors a stripe order: In the presence of $J_{2}$, the stripe Ising configuration [Supplementary Fig. S4(d)] leads to an energy estimate of $-J_{1}-J_{2}$, while the up-up-down (UUD) order [A,B sublattice spin up and C sublattice down, see Supplementary Fig. S4(b)] $-J_{1}+3J_{2}$, and the order in Supplementary Fig. S4(c) [math] 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 AgNiO2, the existence of stripe order was observed, which is ascribed to the relatively large $J_{2}$ in the compound Wheeler et al. (2009).

Supplementary Note 4: TLI parameter fittings.

The parameter fitting workflow is as follows: we scan the parameters $(J_{1},J_{2},\Delta)$ to fit the specific heat $C_{\rm{m}}(T)$, magnetic entropy $S_{\rm{m}}(T)$, as well as the susceptibility $\chi(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_{\rm{m}}$ at finite fields $h$, and compare directly to experimental data Li et al. (2020); Shen et al. (2019), so as to ensure that the parameter set is adequate and precise to model the material.

To be concrete, we show in Supplementary Fig. S5 part of simulation data in our scanning. In Supplementary Fig. S5(a), we start with $J_{2}=0$ and scan various $\Delta$ values. It is found that $C_{\rm{m}}$ curves are sensitive, in terms of the peak height as well as the overall shape, to different $\Delta$ values. By tuning $\Delta$ (while keeping $J_{2}=0$), we find $\Delta/J_{1}=0.4$ and $J_{1}=1.1$ meV can produce results in agreement with experimental $C_{\rm{m}}$ curves. However, with this set of parameters (as well as other $\Delta$ values) we clearly miss the experimental susceptibility line, as plotted in Supplementary Fig. S5(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 $\Delta/J_{1}=0.54$ with $J_{2}/J_{1}=0.05$ can well reproduce both the specific heat curve in Supplementary Fig. S5(c) and the magnetic susceptibility data in (d). To show how sensitive the fittings are with respect to $J_{2}$, we also provide in Supplementary Figs. S5(c,d) the simulated data with $J_{2}=0.01$ to 0.09, from which we see that $J_{2}$ considerably influences both $C_{m}(T)$ and $\chi(T)$ curves. Regarding error bar of the fitted parameter $J_{2}$, the most probable regime where $J_{2}/J_{1}$ resides is between 0.03 and 0.05 (while 0.05 is still more preferable).

With this parameter set $J_{1}=0.99$ meV, $\Delta/J_{1}=0.54$, and $J_{2}/J_{1}=0.05$ (as well as $g_{\parallel}=13.212$), we have computed the magnetization curves at two different temperatures and entropy $S_{\rm{m}}$ at finite magnetic fields. We compare them to experimental data in Fig. 2 of the main text and observe excellent agreement. 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 $\omega(k)$ with same parameter set on a $36\times 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. S6 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 $\Theta_{\rm{W}}\simeq 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\simeq$ 210.1 cm3K/mol leads to an estimate of the effective $g_{\parallel}=\frac{1}{\sqrt{S(S+1)}\mu_{\rm{B}}}\frac{\sqrt{3k_{\rm{B}}C}}{\sqrt{N_{\rm{A}}\mu_{0}}}\simeq 13.212$, where $S=1/2$ represents the effective spin-1/2, which is in excellent agreement with the various thermodynamic fittings till low temperatures in the main text and thus constitutes a self-consistency check. Remind this $g_{\parallel}$ value obtained is in reasonably good agreement with the ideal Landé factor $g_{J}=7/6$ (ideally $J_{\parallel}=2Jg_{\rm{J}}$).

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, $\Delta/J_{1}=0.54$, $J_{2}/J_{1}=0.05$, and $g_{\parallel}=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) earlier in the main text, the magnetization curves $M(h)$ show quasi-plateau structures at $M/M_{\rm{sat}}\simeq 1/3$, where $M_{\rm{sat}}$ is the saturation magnetization along the $z$ direction. In Supplementary Ref. Li et al. (2020), the specific heat curves have also been measured under various external magnetic fields. In Supplementary Fig. S7(a), we redrawn the experimental data, and compare them, side by side, to the simulated $C_{\rm{m}}(h,T)$ curves. It can be seen that in both Supplementary Figs. S7(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\geq 30$ kOe.

It is remarkable that these highly nontrivial and non-monotonic behaviors of $C_{\rm{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_{\rm{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_{\rm{m}}$, despite that the peaks in $C_{\rm{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_{\rm{m}}$ peak moves back to lower $T$ side with a decreasing height. Eventually, for $h\geq 30$ kOe, the system gradually polarizes into a ferromagnetic spin configuration, and the hump in $C_{\rm{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 $\Delta$ limit), we consider a “tight-binding” model

[TABLE]

on the triangular lattice, where $i$ labels the lattice site, and $|i\rangle$ labels a state with a spin flipped at site $i$. $\delta=\mathbb{a},\mathbb{b},\mathbb{a}-\mathbb{b}$ denotes nearest neighboring sites ($\mathbb{a},\mathbb{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\rangle$ and $|j\rangle$, so $t\sim\Delta^{2}$ (in the small $\Delta$ limit), given $i$ and $j$ constitute a pair of nearest neighboring sites.

For the sake of simplicity, we set $\epsilon_{0}=0,t=1$, and take Fourier transformation of Eq. (S1). The resulting dispersion $\epsilon(k)=2[\cos(k\cdot\mathbb{a})+\cos(k\cdot\mathbb{b})+\cos(k\cdot(\mathbb{a}-\mathbb{b}))]$ is plotted in Supplementary Fig. S8(a). By cutting along the $\Gamma$-M-$\Gamma$ path, we observe a quadratic low-energy dispersion near the minimum as shown in Supplementary Fig. S8(b). On the other hand, as shown in Supplementary Fig. S8(c) the M point constitutes a maximal along K-M-K path. Therefore, it is evident that the M point indeed constitutes a saddle point in the dispersion.

Supplementary Note 8: A brief introduction of path integral QMC for TLI.

The Hamiltonian of quantum TLI is $\mathcal{H}=J_{1}\sum_{\langle i,j\rangle}S_{i}^{z}S_{j}^{z}+J_{2}\sum_{\langle\langle i,j\rangle\rangle}S_{i}^{z}S_{j}^{z}-h\sum_{i}S_{i}^{x}$, where the external magnetic field part is omitted for clearance of narrative. QMC evaluates the partition function $Z=\text{Tr}\left[e^{-\beta\mathcal{H}}\right]$, where discretized imaginary time is used and $\beta=L_{\tau}\Delta\tau$ ($\Delta\tau=0.05$). It maps the problem onto a (2+1)D classical system as follows,

[TABLE]

By inserting a complete set of $S^{z}_{\rm{i}}$ eigenstates between each pair of exponentials, i.e.,

[TABLE]

we can then rewrite the partition function as,

[TABLE]

where $l$ indices the time slice $\tau=l\cdot\Delta\tau$.

Next, we employ the Trotter-Suzuki decomposition,

[TABLE]

where $\mathcal{H}_{0}=J_{1}\sum_{\langle i,j\rangle}S_{i}^{z}S_{j}^{z}+J_{2}\sum_{\langle\langle i,j\rangle\rangle}S_{i}^{z}S_{j}^{z}$ and $\mathcal{H}_{1}=-h\sum_{i}S_{i}^{x}$. Therefore the corresponding matrix element reads

[TABLE]

with $\gamma=-\frac{1}{2}\ln\tanh(\Delta\tau h)$, and $\Lambda^{2}=\sinh(\Delta\tau h)\cosh(\Delta\tau h)$. For a certain configuration $\{S^{z}_{i,l}\}$, the configurational weight is,

[TABLE]

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

[TABLE]

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,

[TABLE]

which should be calculated in prior to the spin spectrum $S(\mathbf{q},\omega)$. The latter can be obtained after a stochastic analytical continuation, as detailed in Methods section of the main text.