Topological phase transition and single/multi anyon dynamics of Z2 spin liquid

Among the quantum many-body models that host anyon excitation and topological orders, quantum dimer models (QDM) provide a suitable playground for studying the relation between single-anyon and multi-anyon continuum spectra. However, as the prototypical correlated system with local constraints, the generic solution of QDM at different lattice geometry and parameter regimes is still missing due to the lack of controlled methodologies. Here we obtain, via sweeping cluster quantum Monte Carlo algorithm, the excitation spectra in different phases of the triangular lattice QDM. Our results reveal the single vison excitations inside the Z2 quantum spin liquid (QSL) and its condensation towards the 12×12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{12}\times \sqrt{12}$$\end{document} valence bond solid (VBS), and demonstrate the translational symmetry fractionalization and emergent O(4) symmetry at the QSL-VBS transition. We find the single vison excitations, whose convolution qualitatively reproduces the dimer spectra, are not free but subject to interaction effects throughout the transition. The nature of the VBS with its O(4) order parameters are unearthed in full scope. Our approach opens the avenue for generic solution of the static and dynamic properties of QDMs and has relevance towards the realization and detection of fractional excitations in programmable quantum simulators.


INTRODUCTION
Fractionalized anyon excitations are among the most important features of topologically ordered phases, a class of phases beyond the Landau paradiam of classifying phases with symmetry breaking 1 . The fractionalized nature of these anyon excitations renders that they cannot be created or annihilated individually by physical probes. This phenomenon is both a blessing and a curse: it reflects the topological nature of the excitations and the phase, but also obscures any direct detection of single anyon excitations. Instead, they can only be observed indirectly from a multi-particle continuum of spectral functions. For example, a continuum in inelastic neutron scattering spectrum is often used as a signature to detect quantum spin liquids with fractionalized spin excitations, which is considered as a two-spinon continuum [2][3][4][5][6] . Consequently, understanding the relation between physical spectra and underlying single-anyon excitation is an essential question in the study of topologically ordered phases.
In the theoretical study of topologically ordered phases including the quantum spin liquid (QSL) 7,8 , one usually relies on approximate tools to model the fractionalized excitations because they cannot be directly accessed in experiments and numerical simulations. In simple mean-field theories of QSL, as a physical probe excites a pair of fractionalized excitations, the corresponding spectrum is given by the convolution of spectra of the underlying anyons. However, in realistic systems, this simple relation is modified by interactions between anyons, it is therefore important to know how much change has happened due to the interaction effect.
Quantum dimer models (QDM) 9,10 provide a suitable playground for studying the relation between single-anyon and twoanyon spectra in QSLs. Originally proposed to model the resonant valence bond state in high-T c superconductors 11 and frustrated magnets, it realizes a gapped Z 2 QSL at the exactly-solvable Rokhsar-Kivelson (RK) point 10 if put on a nonbipartite lattice such as the triangle and the kagome [12][13][14] . Comparing to other models of QSLs, the QDMs are suitable as the spinful excitations are absent in the Hilbert space due to the one-dimer-per-site constraint. This means that the spinon excitations in the Z 2 spin liquid are absent, leaving the visons as the only low-energy anyon excitations. As a result, the spectrum of vison excitations can be directly measured in numerical simulations. This feature of QDM allows one to compare the spectra of both the fractionalized single-vison excitations and the physical dimer-dimer correlations, which involves a pair of visons. Although the ground state of QDM is exactly known at the RK point, the excited states are not exactly solvable due to interactions among the visons.
Furthermore, away from the RK point, the QDM on the triangular lattice can be tuned into a ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p valence bond solid (VBS) phase 12,13 . The phase transition is conjectured to be continuous and of the O(4) universality, driven by the condensation of visons [15][16][17][18] . Therefore, the QDM near this transition is an ideal system to study the spectral properties of anyon condensation, if there exist controlled theoretical and numerical methods.
Recently, a quantum Monte Carlo (QMC) scheme, the sweeping cluster method, is invented by the author [19][20][21] . The method is able to keep track of the strict local constraint of dimer covering and at the same time perform Markov chain Monte Carlo (MC) in the space-time path integral such that both static and dynamic properties of the QDM can be obtained, only subject to finite system sizes. Therefore, it is different from the projection QMC employed in the previous literatures [15][16][17][18] , where the interplay of quantum and thermal fluctuations of the QDM models is not present, and the computation complexity has been reduced such that larger system sizes can now be accessed (as will show later, the largest system size is three times larger than that in previous literature). The method has been applied to the square lattice QDM and a mixed phase separating columnar phase at strong dimer attraction and staggered phase at strong dimer repulsion are found 20 . In this work, we further develop the method to study the static and dynamic properties of triangular lattice QDM.
The problem has a long and interesting history. From the work of Moessner and Sondhi 12 , one knows that from the mapping to frustrated Ising model on honeycomb lattice, the problem is in principle solvable via MC simulations on the frustrated Ising model, and a ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p VBS and a Z 2 QSL are suggested. Then in a series of works with zero temperature Green's function MC [15][16][17][18] , the transition from the QSL to VBS, with the notion that the gap of topological vison excitations is closed at the transition is revealed, although the numerical method therein only work close to the RK point and zero temperature. Later, the dynamical dimer correlations at the RK point is presented in ref. 22 , taking the advantage that at the RK point, the quantum mechanics in imaginary time among the equally weighted dimer coverings is equivalent to a classical stochastic process 23 . Despite of these important progresses, the complete spectra of both dimer and vison excitations, not only the gap but also the spectral weight, and the complete understanding of the transition from Z 2 QSL to the ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p VBS in terms of symmetry fractionalization of topological order, and the nature of the complex O(4) order parameter of the VBS, are not revealed. Here we try to answer these questions with unbiased QMC and symmetry analysis.

Model and measurements
We study the QDM on triangular lattice with the Hamiltonian, where the sum runs over all plaquettes including the three possible orientations. The kinetic term, controlled by t, flips the two dimers on every flippable plaquette, i.e., on every plaquette with two parallel dimers, while the potential term V describes interactions between nearest-neighbor dimers. Throughout the paper, we set t = 1 as the energy unit and the inverse temperature β = 1/T with temperature scale also measured according to t.
Before the sweeping cluster QMC [19][20][21] , one commonly employs the projector approaches to study QDMs, which includes the Green's function [15][16][17][18] and diffusion MC schemes 24,25 . These projector methods obey the geometric constraints, but are not effective away from RK point 26 and only work at T = 0. Also, there exists no cluster update for the projector methods to reduce the computational complexity. On the contrary, the sweeping cluster algorithm is based on path-integral in the world-line MC configurational space of all finite temperatures and features efficient cluster update for constrained systems. It is an general extension of the directed-loop algorithm 27,28 for the D dimension classic dimer model 29 to the quantum dimension of (D + 1). Since our QMC works at finite temperature, we can also access the imaginary time correlation functions. And from here, we employ the stochastic analytic continuation (SAC) method [30][31][32][33][34][35][36][37][38][39] to obtain the real frequency excitation spectra from their QMC imaginary time correspondance. The reliability of such QMC-SAC scheme has been extensively tested in quantum many-body systems, ranging from 1D Heisenberg chain 33 compared with Bethe ansatz, 2D Heisenberg model 36,38 compared with exact diagonalization, field theoretical analysis and neutron scattering of square lattice quantum magnet, Z 2 quantum spin liquid model with fractionalized spectra 35,39 compared with anyon condensation theory to quantum Ising model with direct comparison with neutron scattering and NMR experiments 40,41 .
We compute three dynamical correlation functions. The first one is dimer correlation. The dimer operator D i = 1 or 0 when there is a/no dimer on the link i. The dimer correlation function is defined as C d ðr i;j ; τÞ ¼ P i;j hD i ðτÞD j ð0Þi À hD i i 2 , and C d (q, τ) through the Fourier transformation, then the excitation spectrum The second one is vison correlation. Visons (V i ) live in the centre of triangle plaquettes and they must arise in pairs, as shown in the right inset of Fig. 1c. The correlation function is defined as Pij is the number of dimers along the path P ij we chose between plaquettes i and j as shown in Fig. 1c. It is clear that the value of V i V j is path dependent. In order to eliminate this dependence, one can choose a reference configuration, and follow the same path P ij again to obtain another N 00 Pij and then the observable Pij is path independent. Then we redefine C v ðr i;j ; τÞ ¼ hV i ðτÞV j ð0Þi ¼ hV i ðτÞV i ð0ÞV i ð0ÞV j ð0Þi ¼ hðÀ1Þ NH t þNP ij i where N Ht means the number of the t-term in Eq. (1) between V i (τ) and V i (0). We choose the reference configuration as the columnar VBS shown in Fig. 1a, which doubles the unit cell and the corresponding BZ under this reference (gauge choice) is the dashed rectangle with high symmetry points A, B and C in Fig. 1b.
The last one is another "dimer", i.e., the vison-convolution correlation function. We denote this "dimer" -the visonconvolution (VC) operator -as D vc The idea is that if two visons are closest to each other, sharing the same link, then the D vc i on link i can be represented as the product of these two vison operators, with i 1 and i 2 the triangle plaquettes closest to the link i. d i = ±1 when there is no/one dimer on link i of the reference 2 can be computed using Wick's theorem as the convolution of two vison operators, (2) Here, d i is constant for link i under the gauge choice, and can be taken outside the brackets. The spectrum C vc d ðq; ωÞ, which we refer to as the vison-convolution spectrum, gives rise to the convolution of two vison excitations. It is therefore of great importance to compare it with the dimer spectrum C d (q, ω), where the difference will reveal the interaction effects among the visons in different regions of the phase diagram. And we emphasize that although the bottom of the dimer dispersion has been discussed in the refs. 17,18 , the full numerical calculation of the C d (q, ω), C v (q, ω) and C vc d ðq; ωÞ dynamical correlation functions, both in the frequency and momentum axes, are being presented here and they provide the well-characterised example of the dynamics of a Z 2 spin liquid and a phase transition driven by condensation of fractional excitations.

Spectra of dimer, vison and vison-convolution
In the Z 2 QSL phase, the visons are the emergent and fractionalized elementary excitation with no spin and charge quantum numbers 42 . As discussed in the introduction, this is a suitable advantage of the QDM that single vison spectrum can be measured unambiguously, as usually the vison excitations have to be constructed in mean-field as built-in without knowing the unbiased physics 43 , or measured indirectly in lattice models of frustrated magnets 35,39,[44][45][46][47] .
We therefore measure the correlation functions of C d (q, τ), C v (q, τ) and C vc d ðq; τÞ in QMC and then using SAC [33][34][35] to generate the real frequency spectra C d (q, ω), C v (q, ω) and C vc d ðq; ωÞ. These results are presented in Fig. 2. Inside the Z 2 QSL phase with V = 1, all the spectra are gapped. The vison spectra acquire the smallest gap at the order of ω~0.1 at B point of BZ. And the dimer and VC correlations are also gapped with their minimal at M point. It is interesting to notice that the VC spectral gap at M point is higher than the dimer gap at the same point, suggesting that actually visons have a binding energy in forming the dimer correlation and consequently their interaction effect is attractive and gives rise to a bound state with lower energy than the naive convolution. In addition to SAC, we also fit the excitation gaps directly from the imaginary time correlation functions, as shown in Supplemental Notes.
As V is reduced from 1 to 0.9 and 0.8, a QSL-VBS transition is expected at V c~0 .85 [15][16][17][18] , and previous works from the gap measurements and field analytical analysis 12 have proposed emergent O(4) symmetry at the transition. But how the entire spectra change across the transition has not been shown due to the lack of access to finite temperature fluctuation effects. With our QMC + SAC scheme, we observe that the vison gap closes at the B point and the dimer and VC spectrum gap close at X and M points of the BZ (subject to finite size effect of the QMC simulation), as shown in Fig. 2 for V = 0.8 and 0.9. The minimal at Γ and K of the VC spectra come from the allowed momentum convolution of single vison spectra which has minimal at B. Such gap closing process is a manifestation of the symmetry fractionalization mechanism of anyon condensation in Z 2 topological order 35,[48][49][50] . That is, since here the Z 2 gauge field is odd in nature (see the discussion in Supplemental Notes), the visons carry π-flux throughout the lattice. As the QSL-VBS critical point is approached, the vison gap will close and the entire vison spectral weight will condensed at a finite momentum point. In a similar manner, the dimer spectra, generated from the vison bound states, will also close at a finite momentum point. This is different from the usual Bose condensation from disorder symmetric state to ordered symmetry-breaking state, where the condensation of the low-lying bosons usually close gap at the Γ point. Since in our case the disordered state has intrinsic topological order with elementary excitations (visons) carrying finite momentum (π-flux), the condensation gap manifests finite momentum closing. Similar translation symmetry fractionalization process, has also been observed in π-flux Z 2 spin liquid realized in the Kagome lattice model 35,[48][49][50] , which is proposed to be used as an experimental

Dimer
Vison-convolution Vison signature of quantum spin liquid in neutron scattering for Kagome antiferromagnet [51][52][53] . Also, one sees that at V = 0.9 and 0.8, there are more higher energy spectral weights in dimer, VC and vison spectra, coming from the enhanced quantum critical fluctuations of the QSL-VBS transition.
Emergent O(4) symmetry and order parameter of VBS Next we discuss the nature of the QSL-VBS transition and the symmetry breaking pattern of the ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p phase. As explained in the Supplemental Notes, it is expected theoretically 13,54 that this transition is driven by the condensation of visons, which is described by a four-component order parameter {ϕ i }, i = 0, 1, 2, 3 constructed from the Fourier transformation vison configuration at momenta B, i.e. ±ð π 6 ; π 6 Þ and ±ðÀ π 6 ; 5π 6 Þ in Fig. 1b. The order parameter transforms as a 4D representation under the lattice wallpaper-group symmetries, and the matrix form of group actions are summarized in the Supplemental Notes.
In order to numerically confirm that the order parameter ϕ i indeed captures the QSL-VBS transition, we perform a principal component analysis (PCA) on the vison correlation function C v to extract the condensing mode near the transition. PCA diagonalizes the 4 × 4 matrix of the momentum-space vison correlation function at the B point, and identifies the eigenvectors with the largest eigenvalues corresponds to the modes represented by the order parameter ϕ i . We list the ratio of the first largest eigenvalue over the second at V = 0.5-1 in Table 1. Since the largest eigenvalue always dominate, it shows that the principal component of the VBS structure is indeed the expected ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p order. The theoretical analysis further predicts that, at the QSL-VBS critical point, the transition point acquires an emergent O(4) symmetry, as O(4)-symmetry-breaking terms become irrelevant.
In other words, the order parameter lives homogeneously on a four-dimensional sphere 12 .
To reveal such emergent O(4) symmetry at the QSL-VBS critical point and its breaking inside the ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p VBS phase. We prepare the order parameter histogram in Fig. 3. By reorganizing the order parameters into Since the order parameter is four dimensional and hard to visualize, we draw two-dimensional projected histogram (w, x) and (y, z) of the 4D order parameter near the phase transition point at V = 0.85 and deep inside the VBS phase at V = 0. Figure 3a and b are the two independent projections of the 4D (w, x, y, z) space and clearly an emergent O(4) symmetry is present. The Inset shows the modulus distribution of the 4D sphere (with arbitrary unit) which means the order indeed lives homogeneously on a four-dimensional sphere 55 . Figure 3c and d are the same analysis inside the VBS phase, and here clearly distinctive points in the two projected phase are present, which are in full consistency with the symmetry analysis in Supplemental Notes, i.e. the ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p VBS breaks the O(4) symmetry.

DISCUSSION
Via the sweeping cluster QMC algorithm, supplemented with SAC scheme to obtain the real-frequency data and symmetry analysis of the VBS order parameter, we reveal the excitation spectra in different phases of the triangular lattice QDM, in particular, the single vison excitations inside the Z 2 QSL and its condensation towards the ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p VBS with the translational symmetry fractionalization. We found the vison-convolution spectrum is different from the dimer spectrum due to the vison interaction effect, and we also unearth the emergent O(4) symmetry at the QSL-VBS transition and the nature of the ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p VBS with its O(4) order parameter and symmetry breaking. We note that our results not only confirm expectations on triangular lattice QDM by previous works 12,13,[15][16][17][18] , but more importantly, move forward by directly and reliably characterising the single particle dynamics of fractional excitations using controlled numerics, and demonstrating their condensation towards symmetry-breaking phase. We believe our work provide the well-characterised example of the dynamics of a Z 2 spin liquid and opens an avenue for generic solution of the static and dynamic properties of QDMs and other strict constrained systems, such as those in programmable quantum simulators based on Rydberg atom arrays [56][57][58] and superconducting qubits 59,60 where geometry frustration and dynamics of quantum Ising models have been proposed and partially realized.

Sweeping cluster algorithm
This is a quantum Monte Carlo method developed by author which can work well in constrained spin models [19][20][21] . The key idea of sweeping cluster algorithm is to sweep and update layer by layer along the imaginary time direction, so that the local constraints (gauge field) are recorded by update-lines. Via this way, all the samplings are done in the restricted Hilbert space, i.e. the low-energy space. In this article, we can measure the information of single vison because in a strictly constrained space, the energy gap of other quasi-particles such as spinon, becomes infinite large and thus these quasi-particles does not exist in the restricted Hilbert space. We also note that due to the reduced computational complexity with global updates, the system sizes simulated here is three times larger than those simulated with the projection methods in previous works [16][17][18] .

Stochastic analytic continuation
The main idea of this method [30][31][32] is to obtain the optimal solution of the inverse Laplace transform via sampling depend on importance of goodness. From sweeping cluster method, we can obtain a set of imaginary time correlation functions G(τ). The real-frequency spectral function and the imaginary time correlation function have the following transformation relationship as GðτÞ ¼ 1 π R 1 0 dωðe Àτω þ e ÀðβÀτÞω ÞSðωÞ. In order to inversely solve this equation, we must fit a better spectral function. Let the spectral function has a general form as S(ω) = ∑ i a i δ(ω − ω i ). By sampling according to the importance of goodness of fit, we can finally get the spectral function numerically. The reliability of such QMC-SAC scheme has been extensively tested in quantum many-body systems, ranging from 1D Heisenberg chain 33 compared with Bethe ansatz, 2D Heisenberg model 36,38 compared with exact diagonalization, field theoretical analysis and neutron scattering spectra in real square lattice quantum magnets, deconfined quantum critical point 36,37 and deconfined U(1) spin liquid phase with emergent photon excitations 61 , Z 2 quantum spin liquid model with fractionalized spectra 35,39 compared with anyon condensation theory, to quantum Ising model with direct comparison with neutron scattering and NMR experiments 40,41 .