Quantum dynamics of topological strings in a frustrated Ising antiferromagnet

We investigate the quantum dynamics of the transverse field Ising model on the triangular lattice through large-scale quantum Monte Carlo simulations and stochastic analytic continuation. At weak transverse field, we capture for the first time the excitations related to topological quantum strings, which exhibits continuum features described by XY chain along the strings and those in accord with"Luttinger string liquid"in the perpendicular direction. The continuum features can be well understood from the perspective of topological strings. Furthermore, we identify the contribution of strings from the excitation spectrum. Our study provides characteristic features for the experimental search for string-related excitations and proposes a new theoretical method to pinpoint topological excitations in the experimental spectra.

Excitations violating such local constraints are frequently referred to as spinons, which are the gauge charges in the emergent gauge theory [10][11][12]. Spinons are point-like topological defects since a single spinon cannot be created/annihilated locally and such processes must be involved in pairs. Meanwhile, gauge theories also support string-like excitations, which connect a pair of spinons with opposite charges and mediate the effective interaction between them. While the spinon pair is annihilated along a non-contractible closed trajectory, such string becomes topological in nature, which we dub as 'topological string' [13,14].
In contrast to spinon excitations that are typically located at higher energies, quantum strings are basic degrees of freedom participating in low-energy physics. Their dynamics reflects the fluctuations of the background gauge field [15,16]. Quantum strings are 1d objects in space that can be described as twodimensional worldsheets in spacetime [17]. As emergent topological defects with finite dimensions, strings have much more internal degrees of freedom (corresponding to their vibration) and are thus expected to bring much richer physics. Strings appear in broad contexts in condensed matter physics, such as in hightemperature superconductors [18][19][20][21][22], quantum spin ice [23][24][25] and incommensurability [18,26,27]. They can even be directly observed via the ultra-cold atoms in the optical lattice [28,29]. However, in contrast to well-studied point-like excitations, strings are paid much less attention.
The quest for exploring string dynamics is also motivated by the rapid progress on frustrated rareearth materials. Recently, a number of rare-earth pyroclore, Kagomé and triangular lattice magnets have been evidenced to be described by frustrated quantum Ising models [30][31][32][33][34][35] which exhibit emergent gauge structures. For the triangular lattice, a large class of materials are described by the transverse field Ising model (TFIM) [31,34], a typical model that can be understood from the perspective of topological strings [15,16,29]: String picture can be utilised to understand its phase diagram [15,16] as well as to depict the quantum fluctuations therein [36][37][38][39]. Among these materials, a typical one is TmMgGaO 4 (TMGO). The spin excitation spectra was recently discovered in the neutron scattering experiment [31], and the so-called 'roton' excitation was claimed to relate to the BKT physics [33,40,41]. However, direct signatures of string excitations are absent from the neutron spectrum.
In this manuscript, we study a generic antiferromagnetic TFIM on the triangular lattice relevant to triangular lattice rare-earth compounds. Through large-scale quantum Monte Carlo (QMC) simulations [42], we focus on the dynamical signatures of quantum string excitations. At weak transverse field, we find continuum excitations at low energies that accords with the internal vibration and mutual interactions of quantum strings. In particular, the internal vibration matches well with the continuous spectrum of spin-1/2 XY chain, while the string interactions induce excitations exhibiting features of Luttinger liquid [43]. It should be noted that these continuum spectroscopic properties go beyond the effective field theory, where for the latter only coherent spin waves were predicted at low energies. Our work provides insight into understanding the dynamical properties of triangular quantum Ising antiferromagnet such as TMGO.

Starting from rare-earth compounds
We begin with the nearest-neighbour spin-1/2 TFIM relevant to triangular rare-earth compounds Here is the effective spin-1/2 moment acting on the lowest two crystal field levels of the non-Kramers ion.
denotes the nearest neighbour antiferromagnet coupling and the Δ is the energy splitting between the lowest crystal-field levels. The origin of the crystal field splitting term Δ depends on the nature of crystal field levels.
For non-Kramers doublet systems, the crystal field degeneracy is protected by the 3 on-site point group symmetry, therefore Δ is induced by breaking such symmetry, e.g., by applying in-plane uniaxial pressure.
For TMGO and some other Tm-based materials, the lowest two crystal field levels form a nondegenerate dipole-multipole doublet. The lowest crystal fields exhibit intrinsic splitting Δ since they carry the singlet representation of the point group. To quantitatively compare with the TMGO experiment, we also introduce a next-nearest-neighbour interaction = perturbation term to the Hamiltonian [34]. The effective model is illustrated in Fig. 2c inset.
The best-fit parameters for TMGO were found in the former work [33] to be = 0.99 meV, = 0.05 meV and Δ = 0.54 meV. Hereafter we set the NN coupling as the energy unit = 1. We carry out a QMC simulation [42,44,45] combined with SAC technique [46][47][48] to measure the spin excitation spectrum. The simulation is performed on × lattice ( = 24) with periodic boundary conditions. The temperature is set = 1/4 = 1/96 to mimic the ground state result. In inelastic neutron scattering experiments, what is measured is the dynamical -correlator given by Our measured excitation spectrum is presented in Fig. 1a.
At the parameters for TMGO, the simulated result agrees well with the previous inelastic neutron scattering experiment [31] and the previous QMC work [33] (See Supplementary Figure 1). Notice that for TMGO Δ is large where the system has been close to the clock-paramagnetic transition point (Fig. 2f). To gain a better understanding of the physics deep in the clock phase, we also measure the excitation spectra with decreased Δ. We find that with decreasing Δ, the spectrum splits into two branches. The higher-lying and lower-lying branches lie at the energy scale of and Δ, respectively, see Fig. 1b,c.

The string description
To better understand the two branches of excitations, we first consider the classical Ising limit Δ = 0.
In this case, the ground state is macroscopically degenerate which must satisfy the '2-up-1-down' or '1up-2-down' triangle-rule within each elementary triangle. Such local constraints give rise to emergent (1) gauge structures at low energies which further lead to emergent topological strings and topological sectors. Instead of understanding strings from dimers and gauge theories, here we present a more intuitive way from an alternative domain wall perspective [49][50][51]. We choose the stripe state (Fig. 2d) Figure 4). Under periodic boundary conditions, strings wrapping around non-contractible loops cannot be created or annihilated within finite steps of operation, which strongly reflects its topological feature in nature. Therefore the number of topological strings s can be used to label different topological sectors. On the other hand, breaking of triangle-rule corresponds to spinon excitations at the energy scale of [15,16] (Fig. 2a, string 4) . This is related to the high energy excitations.
Once the transverse field Δ is tuned on, quantum dynamics are present. When Δ , the low energy excitations do not mix up with the spinon excitations at a higher energy level. As shown above, the low energy quantum excitations are completely provided by the quantum strings. Here we first consider the excitations where only one string is involved, i.e., the internal vibration of a single string. The shape of a string is specified by its segments {s 1 , . . . , s }, which take on the values r and l (shown in Fig. 2a). Upon the action of the transverse field, a pair of antiparallel nearest-neighbour spins are flipped, i.e., | . . . rl . . . ↔ | . . . lr . . . (Fig. 2b). According to the first-order perturbation theory, these fluctuations can be mapped onto a spin-1/2 The effective model is exactly solvable under Jordan-Wigner transformation. Therefore, one can predict a continuous spectrum due to effective Jordan-Wigner fermions excitations corresponding to kink-antikink pairs in the language of strings. This vibration also brings an energy favour k = −Δ / to each string, stablizing it energetically against perturbations.
When multiple strings are considered, their interaction becomes significant to the quantum excitations.
The non-crossing condition first imposes a hard-core repulsion on strings. The vibration of strings then turns this hard-core repulsion into a dynamic one [49]. Specifically, when two strings come close to each other (e.g., the string 2 and 3 in Fig. 2a), their vibration will be constrained, reducing the energy gain from internal vibration (See Supplementary Note 1). This loss of vibration energy can be equivalently regarded as the repulsion interaction between them. As former work [27] has identified, the dynamic repulsion can be described by a power-law long-range potential. In other words, the internal vibration of strings can also generate mutual dynamic interaction among them.

Dilute string limit
The vibration inside the strings, together with the repulsion interaction between them, determines the low energy excitations of TFIM. To study these two kinds of excitations separately, we turn to the dilute string limit where the strings are far apart and the interaction between them is weak, so that the internal string vibration is not affected by string interaction. This can be parametrised by a quantity called the string density s = s / . The string density is determined by measuring the nearest-neighbor correlation s = 1 2 1 − 4 r r+x . In particular, the stripe phase and clock phase, whose configurations and string correspondence are depicted in Fig. 2d,e, have averaged string density s stripe = 0 and s clock = 2/3. The string density can be tuned through introducing various perturbations, like the nearest-neighbour exchange anisotropy = ( − 1) [27]. Here denote the NN bonds in -direction (Fig. 2c inset).
The string density s ( , ) as a function of exchange anisotropy and NNN interaction are measured through a Monte Carlo simulation, shown in Fig. 2c (See Supplementary Note 2). By increasing , the system enters the stripe phase from the clock phase through a first-order transition; by decreasing , the string density s decreases from 2/3 to 0 continuously.
We take = 0 and tune 0 < < 1 to fine-tune the string density to a comparatively low lever s = 1/3 and 1/2. The energy hierarchy in the spectrum persists (Fig. 3a). The spinon branch and the string branch stay at the energy scale of and Δ, respectively. To look more carefully into the string excitations, we scrutinise the low energy spectrum along ΓM and ΓKM lines (Fig. 3b-e).
The spectrum along the ΓM line corresponds to the internal string vibration. In the case of nearly independent strings at s = 1/3 (Fig. 3b), the spectrum is continuous and enveloped by two dome-like curves, most extensive at Γ point and shrinking when moving towards M point. We find these features in good agreement quantitatively with the continuous spectrum of spin-1/2 XY chain (Fig. 3b,

String and spinon spectra
Having figured out the relation between different string-related excitations and their features in the spectrum, we now return to the isotropic case (i.e., the exchange anisotropy is tuned to = 1). To locate the excitations of different topological defects in the spectrum, we try to isolate the contributions from string and spinon excitations. We define a kink operator, which singles out the flippable spins at the kinks of the strings (See Supplementary Note 4), and a spinon operator, which singles out the spins on triangles with the same spin polarization. As the kink operator is only sensitive to string excitations, and the spinon operator only to spinon excitations, their dynamic correlations isolate the contribution of two respective excitations -At low field Δ = 0.2 (Fig. 4a,d), the string spectrum has only low energy part and the spinon spectrum have only high energy part; while for higher Δ (Fig. 4c,f), the energy scales of two excitations nearly coincide, and they mix up into what we see in the spin spectrum measured by neutron scattering.
Our work focuses on the dynamical signature of quantum string excitations. We simulate the spectrum of frustrated TFIM on the triangular lattice, which is relevant to a large class of rare-earth magnets including the recently discovered TMGO. At low transverse field, we find continuum excitations at low energies.
We show that the features of this continuum are in good agreement with the internal vibration and mutual interactions of quantum strings. In particular, the internal vibration (spectrum along ΓM line) matches well with the continuous spectrum of spin-1/2 XY chain, and the string interaction (spectrum along ΓKM line) induces excitations which exhibit the features of Luttinger liquid.
Several trials can be considered in the experiment to observe the string-related excitations in the excitation spectra. If we can find materials with appropriate Δ/ , we expect well separation of string-related excitations from spinons in excitation spectrum [31,33,52,53]. In experiments, we can also study anisotropic models by applying in-plane uniaxial pressure. On the one hand, for non-Kramers doublet systems, the crystal field splitting Δ is generated and can be tuned by the applied pressure. This enables discovering string Our work not only provides insight into understanding the dynamical properties of triangular quantum Ising antiferromagnet such as TMGO, but also provides iconic features in the spectrum for the experimental search for emergent strings, as well as deepening the understanding of emergent quantum strings in general.
In SSE, the evaluation of partition function is done by a Taylor expansion, and the trace is taken by summing over a complete set of suitably chosen basis. By writing the Hamiltonian as the sum of a set of operators whose matrix elements are easy to calculate = − and truncating the Taylor expansion at a sufficiently large cutoff , we can further obtain

Stochastic analytical continuation (SAC)
We adopted the SAC method to obtain the spectral function ( ) from the imaginary time correlation ( ) measured from QMC [46][47][48]. The spectral function ( ) is connected to the imaginary time Green's is the kernal function . For spin systems, ( , ) = ( − + −( − ) )/ . To ensure the normalization of spectral function, we further modify the transformation, where ( ) = ( ) (1 + − ) is the renormalised spectral function. is the covariance matrix. A Metropolis process is utilised to update the series in sampling.
The weight for a given spectrum is taken to follow a Boltzmann distribution with Θ a virtue temperature. All the sampled spectral functions are then averaged to obtain the final result.
For more details, see the Supplementary Method 2. As an evidence of the validity of SAC, we also compare the spectral function obtained from SAC with the lowest excitation frequency extracted directly from the imaginary time correlations (see the Supplementary Figure 3).

DATA AVAILABILITY
The data that support the findings of this study are available from the authors upon reasonable request.

CODE AVAILABILITY
The computer codes that support the findings of this study are available from the authors upon reasonable request. [9] Yan, Z. Improved sweeping cluster algorithm for quantum dimer model. Preprint at https://arxiv.org/abs/2011. 08457 (2020).

SUPPLEMENTARY FIGURE 1. COMPARISON WITH NEUTRON SPECTRUM
In Supplementary Figure 1, we show that the QMC-SAC simulation on the antiferromagnetic Ising model with next-neareast-neighbour coupling and transverse field can reproduce the spectrum measured by inelastic neutron scattering [52] on TmMgGaO 4 . As indicated by Ref. [33], the best fit parameters are found to be = 0.99meV, = 0.05meV and Δ = 0.54meV. In the Fig. 1 of main text we have plotted the numerical results of the spectrum at such parameters in logarithm scale coloring so that the two modes can both be seen clearly. To compare this QMC-SAC spectrum with the neutron scattering result, we plotted the spectrum in linear scale coloring and take meV as the energy unit. The spectrum from numerical simulation in Panel a and that from neutron scattering in Panel b agrees considerably well.

SUPPLEMENTARY FIGURE 2. SPECTRA IN LINEAR SCALE
In Supplementary Figure 2, we replot the the spectra at different transverse field Δ in Fig. 1 Here , denotes the effective spins of the string and ± are the ladder operators. The denotes that the sum is taken over all the nearest neighbours. This model is well studied and is exactly solvable under Jordan-Wigner transformation. Its excitations are free fermions in the form of effective "spinons" on the effective spin chain and "kinks" in the spin representation on the original lattice.
The discussion above only considers single string case. When multiple strings appear in one configuration at the same time, they have to conform with the non-crossing condition. This condition is equivalent with the triangle rule in spin representation and the one-dimer-per-site constraint in dimer representation. Such condition restricts the Hilbert space of the strings. Same as the spin-1/2 chain, the size of a single string's Hilbert space is 2 . However, when two strings get close (e.g., the string 2 and 3 in Panel i), certain configurations with crossing strings are forbidden. So that the Hilbert space shrinks and an effective longrange interaction emerges (Panel i). Such interaction is found repulsive and should be algebraic to the distance between strings [27]. On the contrary, violation of this rule results in the triangles with three parallel spins. Such triangle hosts a spinon excitation, which is also known as "vortex" as the winding number of effective (1) spin is non-zero [33]. The vortices are connected by either a string loop or a non-directed string (Panel d).

SUPPLEMENTARY NOTE 2. INCOMMENSURATE ORDER
By decreasing from 1, the string density s , defined as which effectively counts the number of horizontal bonds connecting opposite spins, drops from 2/3 and finally reach zero (Panel a) as the system enters the stripe phase through a continuous phase transition.
Accordingly, the peak of magnetic structure factor

SUPPLEMENTARY NOTE 3. RESEMBLANCE TO THE LUTTINGER LIQUID
In this section, we demonstrate that the features in the spectrum at low string density are reminiscent of the Luttinger liquid when we add an anisotropy As described previously, in the perpendicular direction, the strings can be viewed as 0D fermions moving on 1D line. There exists a long-range power-law repulsive interaction between adjacent fermions. As indicated in Ref.
[? ], this may result in a Luttinger liquid phase at certain parameter region. It is non-zero only on the sites which is envelopped by a corner of the string, or to written explicitly, where e is the displacement from a site to its six neighbours (Panel c). N.b., this definition is equivalent to singling out the "flippable" spins.

SUPPLEMENTARY NOTE 5. COMPARISON WITH SPIN WAVE METHODS
In this section, we compare our string description with the spin wave descriptions that have been formerly used to study the spectrum of rare-earth compounds.
First, the spin wave theories are semi-quantitative theories. Particularly, in frustrated magnets, quantum fluctuation is largely neglected in spin waves. As a matter of fact, the predicted parameters, especially the transverse field Δ, for TMGO by spin wave theory deviates a lot from those extracted by more methods such as QMC. Therefore, tuning parameters away from their physical value is needed to produce spectra that agree with numerics or experiments.
By contrast, our string picture can produce quite simple and effective quantitative predictions in a wide range of system parameters. The most representative example is the envelop of continuum excitations.
Along ΓM direction, we predict that the excitations correspond to kinks inside each string, and the spectra are enveloped by the sine curves 1 = Δ cos √ 3 /2 and 2Δ cos √ 3 /4, which agrees considerably well with numerical results (Fig. 3b in main text). Along ΓKM direction, we predict the excitations correspond to Luttinger liquids above the energy scale of the gap. The nearly gapless point should lie at = 2 − and = 2 − , and the dispersions in vicinity should be linear, which also agree with our numerical results. As far as we know, our string picture is the simplest theory that explains these features altogether.
Meanwhile, some features of our spectrum are hard to be reproduced by spin wave theories. For example, at low transverse field Δ = 0.2, in vicinity of Γ point, we find three instinct peaks in the spectrum, centered at = 0.038, 0.337 and 1.302 respectively, composing 37%, 58% and 5% of the spectrum weight. These three peaks are explained in our string picture as the excitation due to interstring interaction, kink-antikink pair on single string and the spinons, respectively. On the other hand, spin wave theory only predicts the existence of two peaks -two of the three modes in the linear spin wave are degenerate at Γ point required by degeneracy, as shown in Fig. 9 of Ref. [35]. The two modes predicted by LSW corresponds to the highest frequency and lowest frequency peak in the spectrum we find. The peak that makes up the largest portion of weight is missed in the linear spin wave theory. This peak is unlikely to come from the non-linear corrections as well, since its weight is larger than the LSW contribution.
Our string picture becomes more rigorous and provide richer physics when we consider anisotropic NN interactions. In this case, the spatial anisotropy can stabilize the incommensurate phase out of the commensurate clock phase. The incommensurate order wave vector exhibits a linear relation with the string density. In the string picture, by minimizing the string energy, we are able to predict a string densityanisotropy relation that agrees with the numerical results (See Supplementary Note 1, 2). We are also able to clearly visualise the strings in the spin configuration snapshot (See Supplementary Figure 4). Furthermore, at low string density with weak inter-string interactions, the spectra exhibit clearer continuous features with agree perfectly with the effective model.

SUPPLEMENTARY METHOD 1. STOCHASTIC SERIES EXPANSION (SSE)
For the numerical works in this paper, we use a quantum Monte Carlo (QMC) method with stochastic series expansion (SSE) algorithm [42,44,45] to calculate the ground state properties and imaginary time Green function. This method will be briefly introduced below.
In quantum statistics, the measurement of observables is closely related to the calculation of partition function = tr ( exp(− )) / , = tr exp(− ) where = 1/ is the inverse temperature, is the Hamiltonian of the system and is an arbitrary observable. Typically, in order to evaluate the ground state property, one takes a sufficiently large such that ∼ , where is the system scale and is the dynamical exponent. In SSE, such evaluation of is done by a Taylor expansion of the exponential and the trace is taken by summing over a complete set of suitably-chosen basis.
We then write the Hamiltonian as the sum of a set of operators whose matrix elements are easy to calculate.
In practice we truncate the Taylor expansion at a sufficiently large cutoff and it is convenient to fix the sequence length by introducing in identity operator 0 = 1 to fill in all the empty positions despite it is not part of the Hamiltonian.
Diagonal update, where diagonal operators are inserted into and removed from the operator sequence, and cluster update, where diagonal and off-diagonal operates convert into each other, are adopted in update strategy.
In transverse field Ising model where a constant is added into the Hamiltonian for convenience. For the non-local update, a branching cluster update strategy is constructed [44], where a cluster is formed in ( + 1)-dimensional by grouping spins and operators together. Each cluster terminates on site operators and includes bond operators (Panel a). All the spins in each cluster is flipped together at a probability 1/2 after all clusters are identified.

SUPPLEMENTARY METHOD 2. STOCHASTIC ANALYTICAL CONTINUATION (SAC)
For the spectra in this paper we adopted a stochastic analytical continuation (SAC) [46][47][48] method to obtain the spectral function ( ) from the imaginary time correlation ( ) measured from QMC, which is generally believed a numerically unstable problem. This method will be briefly introduced below.
The spectral function ( ) is connected to the imaginary time Green's function ( ) through an integral equation where ( , ) is the kernel function depending on the temperature and the statistics of the particles. We where ( ) = ( ) (1 + − ) is the renormalized spectral function.
In practice, ( ) for a set of imaginary time ( = 1, · · · ) is measured in QMC simulation together with the statistical errors. The renormalized spectral function is parametrized into large number of equalamplitude -functions whose positions are sampled (b) Then the fitted Green's functions˜from Eq. 22 and the measured Greens functions¯are compared by the fitting goodness where the covariance matrix is defined as with the number of bins, the measured Green's functions of each .
A Metropolis process is utilized to update the series in sampling. The weight for a given spectrum is taken to follow a Boltzmann distribution with Θ a virtue temperature to balance the goodness of fitting 2 and the smoothness of the spectral function.
All the spectral functions of sampled series { , } is then averaged to obtain the spectrum as the final result.