# Dynamical Signature of Fractionalization at the Deconfined Quantum Critical Point

###### Abstract

Deconfined quantum critical points govern continuous quantum phase transitions at which fractionalized (deconfined) degrees of freedom emerge. Here we study dynamical signatures of the fractionalized excitations in a quantum magnet (the easy-plane J-Q model) that realizes a deconfined quantum critical point with emergent O(4) symmetry. By means of large-scale quantum Monte Carlo simulations and stochastic analytic continuation of imaginary-time correlation functions, we obtain the dynamic spin structure factors in the and channels. In both channels, we observe broad continua as expected from the collective fluctuation of the deconfined excitations. We also provide field-theoretical calculations at the mean field level that explain the overall shapes of the computed spectra, while also pointing to the importance of interactions and gauge fluctuations to explain some aspects of the spectral-weight distribution. We make further comparisons with the conventional Landau paradigmatic O(2) transition in a different quantum magnet, at which no signature of fractionalization are observed. The distinctive spectral signatures of the deconfined quantum critical point suggest the feasibility of its experimental detection in neutron scattering and nuclear magnetic resonance experiments.

## I Introduction

The deconfined quantum critical point (DQCP), which separates the Néel antiferromagnetic (AFM) and spontaneously dimerized valence bond solid (VBS) phases in (2+1)D quantum magnets, was proposed as an example of continuous quantum phase transition outside the conventional Landau-Ginzburg-Wilson (LGW) paradigm Senthil et al. (2004a, b). The AFM and VBS order parameters both vanish continuously and simultaneously at the DQCP. This scenario is generically not expected within the standard LGW description, where such a case should be realizable only by fine tuning two separate transitions to coincide at special multi-critical points. Multiple field theory descriptions Motrunich and Vishwanath (2004); Senthil et al. (2004a, b); Hermele et al. (2005); Ran and Wen (2006); Senthil and Fisher (2006); Grover and Vishwanath (2013); Lu and Lee (2014); Xu and You (2015); Karch and Tong (2016); Hsin and Seiberg (2016); Cheng and Xu (2016); Wang et al. (2017); Metlitski and Thorngren (2017); You et al. (2017) have been proposed for the DQCP which are believed to be equivalent (or dual) to each other at low energy, including the non-compact CP (NCCP) theory Senthil et al. (2004a, b) and some versions of the quantum electrodynamics (QED) and quantum chromodynamics (QCD) theories Wang et al. (2017). In contrast to the LGW description which formulates the critical theory in terms of the order parameters directly, these gauge theory descriptions of the DQCP are formulated in terms of deconfined degrees of freedom (fractionalized particles and emergent gauge fields). The order parameters on either side of the DQCP can be expressed as different compositions of the fractionalized particles or gauge fluctuations within the same theoretical framework. This mechanism captures the intertwinement of the AFM and VBS orders and provides a natural route beyond the LGW paradigm to a non-fine-tuned quantum critical point between the two distinct symmetry-breaking phases. Based on the physical picture of deconfinement of the experimentally accessible spin excitation into two spinons at the DQCP, a broad continuum is expected in the spectral function. This is in sharp contrast to an LGW transition of the AFM state into a nondegenerate (trivial) quantum paramagnet, where spinwaves (magnons) survive as elementary excitations at the critical point (albeit being highly damped) Chubukov et al. (1994). The aim of this paper is to present a comprehensive numerical study of the signature of magnon fractionalization in the dynamic spin structure factor of a (2+1)D square-lattice spin model hosting a DQCP.

Following the DQCP proposal, intensive theoretical and numerical efforts have been invested in the possibility of unambiguously observing such critical points in lattice models. In the traditional frustrated quantum spin models that exhibit VBS phases, sign problems in quantum Monte Carlo (QMC) and other technical difficulties in methods such as the density matrix renormalization group and tensor product states prohibit studies of the large system sizes needed in order to reliably characterize critical points. However, for generic and universal properties, other “designer hamiltonians” Kaul et al. (2013) can be constructed that do not suffer from QMC sign problems but host the desired phases. Many such studies have pointed to the existence of the DQCP in both two-dimensional (2D) quantum magnets Sandvik et al. (2002); Sandvik (2007); Lou et al. (2009); Sandvik (2010a); Pujari et al. (2013); Block et al. (2013); Shao et al. (2016); D’Emidio and Kaul (2017); Qin et al. (2017a); Zhang et al. (2017) and related (through the path integral) three-dimensional (3D) classical models Nahum et al. (2011, 2015a, 2015b). In these studies it has been observed, e.g., that the order parameters have unusually large anomalous dimensions Lou et al. (2009); Block et al. (2013); Nahum et al. (2015b); Shao et al. (2016); Qin et al. (2017a); Zhang et al. (2017), which is an important deviation from the common 3D Wilson-Fisher fixed point. More concrete evidence of deconfinement has been found by directly probing the length scale associated with the fractionalization process Tang and Sandvik (2013); Shao et al. (2016) and from thermodynamics Sandvik et al. (2011). However, the experimentally most direct signatures of a DQCP, the dynamic spin structure factor , have so far not been calculated in the case of electronic spins (while there are already some intriguing results for an SU(3) symmetric model Assaad and Grover (2016)). Historically, in quasi-1D systems, the experimentally observed spinon continuum, in agreement with calculations for the spin- Heisenberg chain, was crucial in establishing spinon deconfinement. Indications of fractionalized magnetic excitations in 2D quantum spin liquids have also been similarly observed Han et al. (2012); Fu et al. (2015); Feng et al. (2017); Wen (2017); Wei et al. (2017); Feng et al. (2017). Given that is detectable by multiple experimental techniques, including inelastic neutron scattering (INS), resonant inelastic X-ray scattering (RIXS) and nuclear magnetic resonance (NMR), identifying the distinct signatures of fractionalization in at the DQCP will provide a useful guide to experimental searches for DQCPs in magnetic materials. Moreover, due to a recently investigated duality relation between the DQCP and a certain bosonic topological transitions in fermion systems Wang et al. (2017); Qin et al. (2017a); Zhang et al. (2017), similar dynamical signature of fractionalization is also expected in interaction-driven topological phase transitions. Therefore our work also can impact the ongoing efforts in finding experimental accessible signatures on topological phase transitions in strongly correlated electron systems.

In this work, we will investigate a U(1) version of the DQCP on the square lattice, with the easy-plane - (EPJQ) model defined by the Hamiltonian

(1) |

were denotes the spin-1/2 operator on each site and is the singlet-projection operator on the link (between nearest-neighbor sites). The two- and six-spin terms are both illustrated in Fig. 1(a). For this is the previously studied SU(2) - model Lou et al. (2009); Sen and Sandvik (2010); Sandvik (2010a), which is an extension of the original - model (or - model) model Sandvik (2007), with two instead of three singlet projectors in the terms. With three singlet projectors we can go further into the VBS state while still keeping in sign-free QMC simulations. The term with introduces the easy-plane anisotropy that breaks the symmetry down to explicitly. It has been shown Qin et al. (2017a) that when , the EPJQ model exhibits a direct and continuous quantum phase transition between the AFXY and VBS phases, as illustrated in Fig. 1(a), realizing the easy-plane DQCP (while for larger anisotropy, such as , the transition becomes first-order). The XY order parameter has a U(1) rotational symmetry and the VBS order parameter exhibits an emergent U(1) symmetry as the DQCP is approached, and, as argued based on dualities Wang et al. (2017), the two U(1) symmetries combine to form an emergent higher O(4) symmetry exactly at the DQCP.

To make a comparison with the EPJQ model, we will also study an easy-plane - (EPJJ) model,

(2) |

where . The bonds and the bonds correspond to the thin black and the thick blue bonds in Fig. 1(b) respectively. Since the term explicitly breaks the lattice symmetry and pins a corresponding pattern of columnar singlets, the large phase will simply be a trivial, non-degenerate paramagnet. The transition out of the AFXY phase is then the conventional O(2) Wilson-Fisher transition in the 3D XY universality class, as illustrated in Fig. 1(b), which we will contrast with the DQCP. We here take for the anisotropy parameter.

By means of stochastic analytic continuation (SAC) of imaginary-time correlation functions calculated using large-scale QMC simulations, we extract the dynamic spin structure factor over a wide range of momentum () and energy () transfers in both the EPJQ and EPJJ models. The calculations are performed in all phases of the models as well as at the critical points for both the and spin channels.

The rest of the paper is organized as follows: In Sec. II, the theoretical background of the DQCP phenomenon and the consequently expected low energy spectral features are laid out. In Sec. III, we discuss in detail the dynamic spin structure factors in the EPJQ and EPJJ models, as they are driven through their phase transitions. This comparison reveals the distinct spectral features of the DQCP. In Sec. IV, we provide a theoretical calculation of the dynamic spin structure factor at the DQCP which nicely matches the numerical observations. Sec. V summarizes the significance of our findings and their relevance to bridging the DQCP to experimentally accessible information, and points out future directions. A detailed finite-size scaling analysis of the not previously studied quantum phase transition of the EPJJ model and additional discussion of the spin spectra of this model are given in Appendix A and B, respectively.

## Ii Theoretical Expectations of Low-Energy Spectral Features

Before presenting our numerical result, we would like to first provide a theoretical overview of the expected spectral features at low energy, as summarized in Tab. 1 and Fig. 2. Let us define the Néel AFM and VBS order parameters as

(3) |

With the easy-plane anisotropy, the XY order parameter in the AFXY phase is just the planar component of the Néel order parameter .

Deep in the AFXY phase of both the EPJQ and EPJJ models, the low-energy fluctuations are described by the XY model , where is the spin stiffness and is the spin-wave Goldstone mode, such that the XY order parameter can be written as . The XY spin correlation function near momentum is expected to follow

(4) |

where . The imaginary part of this correlation function (with ) is shown in Fig. 2 (a), demonstrating the well-defined magnon mode with linear dispersion. On the other hand, the spin fluctuation is gapped at due to the easy-plane anisotropy, but becomes gapless at . The excitation of corresponds to the spin density fluctuation , which can decay into two gapless magnon modes, each around , so that the total momentum is close to . Therefore, we expect near to be of the form

(5) |

The imaginary part of this correlation function (with ) is shown in Fig. 2 (b). Here the spectral weight of the linearly dispersing mode is suppressed as . As we will see, the low-energy spectral features of the dynamic spin structure factors and match our QMC-SAC results nicely [see (a) and (b) for both Fig. 3 and Fig. 4].

low energy | channel | mode | ||
---|---|---|---|---|

(a) | AFXY | |||

(b) | AFXY | |||

(c) | DQCP | |||

(d) | DQCP | |||

(e) | DQCP | |||

(f) | DQCP |

At the DQCP, the low-energy dynamic spin susceptibility around is expected to be of the from

(6) |

with large values of the anomalous dimension and , characterizing the complete breakdown of a well-defined magnon at the critical point. The imaginary part of these correlation functions are shown in Fig. 2 (c,d) respectively, with and taken from Ref. Qin et al. (2017a) for illustration purpose. Strictly speaking, any non-zero anomalous dimension would imply the breakdown of well-defined magnons, but compare to the small anomalous dimension Guida and Zinn-Justin (1998); Campostrini et al. (2001) at the 3D O(2) Wilson-Fisher transition, we expect to observe a much more prominent continuum at the DQCP in the EQJQ model [as clearly seen in Fig. 3 (b,e)]. This is in sharp contrast to the spectrum at the 3DXY critical point in the EPJJ model [as shown in Fig. 4(b)], where there is essentially no continuum in the gapless channel and that in the gapped channel is much less prominent (though there are also interesting features there that cannot be explained at the level of analysis discussed above).

Another important feature in the spectrum of the DQCP is the gapless excitations at momenta and [as well as by symmetry] in both and channels, with much weaker spectral weight, as shown in Fig. 3 (b,e). Theoretically, they correspond to the (generally non-conserved) SO(5) current fluctuations, where the SO(5) group rotates the Néel and VBS order parameters as a combined vector . The SO(5) current can be written in terms of the combined order parameter as

(7) |

with and . By matching the momentum and the SO(5) symmetry quantum numbers, it is straight forward to identify the and fluctuations around to and respectively, and identify those around to and respectively. The emergent O(4) symmetry at the easy-plane DQCP corresponds to the subgroup of SO(5) that rotates only (keeping invariant), so the currents and are emergent conserved currents at low-energy. Their correlation functions can be calculated based on the QCD theory or the QED theory , where the order parameters are fractionalized as and the current-current correlations are given by

(8) |

with being the SO(5) generator that rotates components. These spectral functions of the currents in the field theory correspond in the lattice model to the spin spectrum around and around . The imaginary part of these correlation functions are show in Fig. 2 (e,f) respectively.

The correlation function of the non-conserved currents and are expected to take a similar form with another anomalous dimension ,

(9) |

They correspond to around and around . As we will discuss in more detail in Sec. III and Sec. IV, all these expected spectral features are qualitatively observed in the QMC-SAC spectrum of the EPJQ model [see Fig. 3 (b,e)], consistent with the QCD or QED description of the DQCP.

In the VBS phase, all excitations (in both and channels) are gapped. There is no low-energy feature in the spectrum that can be reliably predicted at the field theory level. With our QMC-SAC numerics we can easily go into the VBS, however, and we will present results along with the results in the XY phase and DQCP in the next section.

## Iii Numerical Calculations of the Spin Spectra

We here present results for both the EPJQ and the EPJJ models. The key quantity computed in our QMC simulations with the stochastic seies expansion (SSE) method Sandvik (2010b) is the spin correlation function in the imaginary time domain (for ),

(10) |

where and the summation is over all sites of lattice. From the imaginary time data for a set of points, we reconstruct the corresponding real-frequency spectral function by performing a numerical analytic continuation using the SAC method Sandvik (1998); Beach (2004); Syljuåsen (2008); Sandvik (2016); Qin et al. (2017b); Shao and Sandvik (); Shao et al. (2017); Huang et al. (2017). With this method, we average over Monte Carlo importance-sampled spectral functions , from which the dynamic spin structure factor is later obtained as . The intermediate spectrum has the advantage of being normalized to when integrating over positive frequencies only. In the sampling procedure we thus fix the normalization and use the relationship

(11) |

to define the goodness-of-fit between this function and the SSE-computed result (including covariance among the SSE data for different ). The weight for a given spectrum is , with a fictitious temperature chosen in an optimal way so as to give a statistically sound mean value, while still staying in the regime of significant fluctuations of the sampled spectra so that a smooth evaraged spectral function is obtained. The most recent incarnation of the SAC method uses a parametrization with a large number of equal-amplitude -function sampled at locations in a frequency continuum and collected in a histogram, as explained in Refs. Qin et al., 2017b; Shao and Sandvik, ; Shao et al., 2017; Huang et al., 2017; Shu et al., 2017. We refer to these works for technical details.

We have extracted for the EPJQ and EPJJ models in both the and channels. The imaginary time correlations in these channels were independently calculated in different simulations with the stochastic series expansion QMC method Sandvik (1999), implemented in the basis of the and spin components, respectively (i.e., the operators used in the correlators are always diagonal). For the EPJQ model in Eq. (1), we set and define the ratio to be the driving parameter, and for the EPJJ model we define . All results presented here are for square lattices with and periodic boundary conditions, and the inverse temperature . While we estimate some remaining finite size effects on these lattices, by comparing with smaller lattices we have confirmed that the main features of the spectra are stable and should be close to the thermodynamic limit and .

For small , the EPJQ model essentially reduces to an XXZ model, which has an AFXY ground state that breaks the symmetry spontaneously. When is large, the dimer interaction favors a VBS (columnar-dimerized) ground state, which spontaneously breaks the lattice rotation symmetry. In this work, we set the anisotropy parameter to , where we have found the signature of an easy-plane version of the DQCP separating the AFXY and the VBS phases at , based on the finite-size analysis of the critical exponent in our previous work Qin et al. (2017a). The phase diagram of the EPJQ model Fig. 1(a) is similar to those of the - and - models Sandvik et al. (2002); Sandvik (2007); Lou et al. (2009); Sen and Sandvik (2010); Sandvik (2010a); Shao et al. (2016), but the DQCP is in a different universality class due to the lowered symmetry. In the EPJJ model, as increases there is a O(2) Wilson-Fisher transition from the AFXY phase to the trivial columnar dimer phase, as illustrated in Fig. 1(b). The critical point is at as determined in Appendix A. The O(3) vesion of this quantum phase transition has been investigated extensively with various statically dimerized Heisenberg Hamiltonians (see the review in Ref. Sandvik, 2010b), including also a recent calculation of dynamic spectral functions Lohöfer et al. (2015). The O(2) transition, however, has not been investigated in detail with 2D quanum spin models, as far as we are aware.

Turning now to the salient features of the spin spectra, in Fig. 3 and Fig. 4 we show our results for the two models along the high symmetry path of wavevectors ---. In both cases we present results both inside the two phases and at the critical point.

In the AFXY phase, which is common to both the EPJQ and the EPJJ models, we observe the gapless Goldstone mode at in the chancel, as shown in Fig. 3 (a) and Fig. 4 (a). In the channel, the fluctuations are gapped due to the easy-plane anisotropy. However, as seen in Fig. 3 (b) and Fig. 4 (b), the modes around are still gapless, but with vanishing spectral weight as , as expected due to the conserved total . These behaviors are consistent with theoretical expectations in Fig. 2 (a,b) based on the field theory of the XY model.

Now let us focus on the critical points of both models. Fig. 3 (b,e) show the EPJQ spectra at , close to the DQCP at . Fig. 4 (b,e) show the EPJJ spectra at , close to the 3DXY transition at . By comparison, several exotic features of the DQCP spectra can be unambiguously identified. First, we observe broad and prominent continua in both and at the DQCP, which reflects the expected magnon fractionalization and emergence of deconfined, essentially independently propagating spinons. In contrast, at the 3DXY transition, the gapless magnon mode remains sharp in around with very weak continuum due to the critical fluctuations.

In the case of the DQCP, we find that the lower edges of both spectral function can be well accounted for by a remarkably simple single-spinon dispersion relation, , which matches the dispersion relation of a deconfined fermionic parton in the square lattice -flux state Hermele et al. (2005); Senthil and Fisher (2006); Ran and Wen (2006); Assaad and Grover (2016). This points us to the QCD or the QED theory Wang et al. (2017); You et al. (2017) that were previously proposed to describe the DQCP. If indeed the broad spectral functions seen in Fig. 3 (b,e) are due to two independently propagated spinons, and contributions from four- or more spinon excitations can be neglected, the upper spectral bound is obtained by maximizing with . This indeed appears to be in reasonable agreement with the observed distribution of the main spectral weight, though some weight, presumably arising from states with more than two spinons, is also present at higher energies. As we will elaborate further in Sec. IV, the overall shape and weight distribution of the continuum can be nicely captured at the mean field level.

Note that the gapless continuum in around is present at the DQCP but is absent at the 3DXY transition. The excitations are simply the low-energy fluctuations of the field. In the NCCP description of the DQCPSenthil et al. (2004a, b), , the fluctuation corresponds to the two-spinon excitation . In the dual NCCP description, , the fluctuation corresponds to the dual gauge flux . Thus, either the fractionalized spinon or the deconfined gauge field will enter the response function in the channel and reflect their criticality in the gapless continuum around . However for the 3DXY transition, described by the XY order parameter in the LGW paradigm, the fluctuation is not associated with any critical fluctuation of the order parameter and thus remains gapped at the critical point. Moreover, at the DQCP, as a consequence of the criticality, the fluctuation also becomes gapless around , because this corresponds to spin density which can decay into the gapless continuum of both and .

Furthermore gapless continua are also observed at momentum [and at as well by symmetry] in both and at the DQCP only. Such exotic features are not observed at the 3DXY transition. The excitation around corresponds to the fluctuation of the conserved current associated to the emergent O(4) symmetry (in the XY-VBS rotation channel), which is an unique feature of the easy-plane DQCP. The gapless point also follows naturally, because the XY-VBS current can deday into the continuum at and the continuum at , such that the momenta add up to . A similar interpretation applies to the channel as well. The only difference is that the spin-VBS current there is not conserved, but is nevertheless still critical. These “shadow” continua allow us to probe the critical VBS fluctuation in the spin excitation spectrum, which is another remarkable hallmark of the DQCP.

As discussed in the Sec. I, the spectral features uncovered here are relatively easy to probe in INS or RIXS experiments, hence paving way for observation of the seeming ephemeral DQCP in real materials. Whereas the previous studies of DQCP mainly focused on the critical scaling and exponents from the theoretical perspective, these quantities are rather difficult to measure in experiments. Even if the DQCP turns out to be first-order (as expected if the anisotropy is strong) or becomes unstable against other intermediate phases at low temperature, its distinct spectral features over a large range of frequencies can still be robustly observed above the low energy scale at which the potentially other transitions of phases become manifest.

Finally, the spectra of the EPJQ model in the VBS phase is shown in Fig. 3 (c,f). Their EPJJ counterpart in the columnar singlet phase is shown in Fig. 4 (c,f). All spin excitations are gapped in both and for both models. For the EPJQ model, the spectra in the VBS phase still maintain broad continua above the gap, in contrast to the much sharper spectra of gapped magnons in the EPJJ columnar phase. This might be related to the two-length-scale phenomena, which is inherent to the DQCP, persisting in the VBS phase of the standard JQ model Shao et al. (2016), namely, the domain wall size of the VBS order may still remain large while the spin correlation length is small. The domain wall size of the VBS order is directly related to the confinement length scale of the spinons. This implies that although the spin correlation length is finite, the confinement length scale of the spinon can still be large, which leads to the large continuum above the spin gap in the spin excitation spectrum.

## Iv Parton Mean Field Theory for the DQCP Spectra

In this section, we provide theoretical account for the overall shape of the dynamic spin structure factors and observed at the DQCP. The easy-plane DQCP admits several candidate field theory descriptions, including the easy-plane NCCP theory Motrunich and Vishwanath (2004); Senthil et al. (2004a, b), the non-compact QED theory Senthil and Fisher (2006); Grover and Vishwanath (2013); Xu and You (2015); Karch and Tong (2016); Hsin and Seiberg (2016); Cheng and Xu (2016); Wang et al. (2017) and the QCD theoryRan and Wen (2006); Wang et al. (2017) (or its Higgs descendent compact QEDHermele et al. (2005); Senthil and Fisher (2006); Wang et al. (2017); You et al. (2017)) with additional anisotropy in the SO(5) symmetric tensor representation. Although all theories are believed to provide equivalent descriptions of the low-energy physics under proposed duality relations Wang et al. (2017), some of them are more convenient to handle by mean-field treatment than others. Among these theories, we found that the QCD (or QED) theory gives the best account for the overall spectral features at the mean field level. Because in these theories, both the AFM and VBS order parameters are treated on equal footing as fermionic parton bilinears, it is already possible to approximately capture both spin and dimer fluctuations at the free fermion level (ignoring gauge fluctuations and local interactions). Fig. 5 shows the comparison of the dynamics spin structure factors between numerics and theory, based on the parton mean field theory. The overall features match quite nicely.

Let us start with the parton construction on the square lattice Wen (2004), where the spin operator is fractionalized into fermionic partons at each site as

(12) |

An SU(2) gauge structure emerges in association with the above fractionalization scheme, but at the mean-field treatment we will ignore the SU(2) gauge fluctuation completely and place the fermionic parton in the square-lattice -flux state Wen (2004); Hermele et al. (2005); Ran and Wen (2006). Thus, we use the following mean-field Hamiltonian

(13) |

such that each plaquette hosts a -flux for the fermionic parton. Four Dirac fermions are obtained at low energy. The fermionic parton dispersion is simply given by

(14) |

It is interesting to find that the lower edge of the DQCP spectra follows this simple dispersion relation quite nicely without any adjustable parameters beyond an overall velocity, as shown in Fig. 5 (a,c), which justifies the -flux state as our starting point. The upper edge of the two-parton continuum can also be obtained from by adding up single-parton energies. This gives a rough estimate for the energy range of the parton continuum, which is also consistent with the numerical observation in Fig. 5 (a,c).

Given Eq. (12) and Eq. (13), it is straight forward to calculate the spin-spin correlation function

(15) |

on the free fermion ground state of the mean field Hamiltonian . Then we can obtain the dynamic spin susceptibility

(16) |

from which we obtain the dynamic spin structure factor

(17) |

graphed in Fig. 6. This spectral function was also calculated in Ref. Assaad and Grover, 2016 previously. One can see that already captures the gapless continua at momenta , , , and in all spin channels. Because the mean field Hamiltonian is symmetric under SU(2), there is no difference between and . The easy-plane anisotropy only enters the parton theory starting from four-fermion interactions, since it is expressed in the SO(5) symmetric tensor representation that can not be written down at the quadratic level. Therefore, the anisotropy is not manifest in the mean field approximation, where the interaction effects are ignored. This observation provides a natural explanation for the strikingly similar spectra of and seen in the numerical results in Sec. III at the DQCP, despite of the presence of a rather large anisotropy in the EPJQ model.

The gauge fluctuations are expected to further renormalize the spectrum and enhance the critical fluctuations around , which are not taken into account in the simple mean field theory presented in Fig. 6. While including the gauge interactions in the calculation is highly non-trivial and beyond the scope of this work, we next discuss a phenomenological model that captures the spectral weight enhancement, and leave the more extensive calculation to future work. Let us consider modeling the interaction effect phenomenologically by a random phase approximation (RPA) correction,

(18) |

where . The coupling parameterize the strength of the spin-spin interaction in the channel. We can introduce the easy-plane anisotropy simply by considering . We found that the fluctuation is indeed enhanced by the interaction . The resulting RPA corrected spectral functions are already shown in Fig. 5 (b,d), with tuned to the magnetic ordering critical point and ^{1}^{1}1Although such a Grose-Neveu critical point is different from the DQCP, we only use it to provide a rough estimate of the spectral features close to a magnetic ordered phase. We do not claim that the criticality of DQCP can be correctly understood by our mean field + RPA approach.. Compared to Fig. 6, the spin spectra in Fig. 5 (b,d) are much improved by the interaction effect. Our phenomenological study combined with the QMC-SAC result demonstrates that the -flux state fermionic parton with interaction accounts well for the overall features of the DQCP spectra in both and channels, which is consistent with the expectations from the QCD or QED theories. An interesting open problem is a systematic route to incorporating the effects of gauge fluctuations in calculating the spin excitation spectrum.

## V Discussion

In this work, we have demonstrated dynamical signatures of fractionalization at the DQCP in a planar, U(1), quantum magnet by computing both the in-plane and out-of-plane dynamic spin structure factors at low temperature. By contrasting with analogous results for a conventional LGW critical point, we explicitly observe how fractionalization of the critical magnon into two spinons is manifested by a large continuum, in sharp contrast to a much less prominent continuum due to conventional critical quantum fluctuations at the ordinary 3DXY transition. We also discovered several low-energy spectral features that are unique to the DQCP, notably the continuum in the channel and the continua in both channels. These features are missing in the 3DXY transition. They will provide us the smoking-gun evidences to guide the search for the DQCP on the experimental venue.

In a series of recent works Karch et al. (2017); Potter et al. (2016); Wang et al. (2017), it was shown that the easy-plane DQCP is dual to the quantum electrodynamics with fermionic matter fields, and it describes the bosonic topological transition (BTT) between a bosonic symmetry protected topological state and a trivial insulator state Grover and Vishwanath (2013); Lu and Lee (2014). In addition, it also exhibits a self-duality Xu and You (2015); Karch et al. (2017); Mross et al. (2016); Hsin and Seiberg (2016). Recently this bosonic topological transition has also been realized in a lattice model and the static properties of the phase transition has been simulated via determinantal QMC Slagle et al. (2015); He et al. (2016); Wu et al. (2016). It would be very interesting to measure the dynamic structure factor at the BTT, and compare the results with our current study.

The DQCP also has a natural large- generalization Block et al. (2013); D’Emidio and Kaul (2017), i.e. instead of two flavors of bosonic matter fields, there are -flavors of matter fields. The field theory is called the noncompact CP model. This model is understood very well in the large- limit, and for finite and large , a systematic expansion can be performed to understand the details of the NCCP. In the future a similar dynamic structure factor measurement of large- version of the DQCP and comparison with the theoretical calculation would be another very interesting research direction. We would also like to compare these calculations to neutron scattering on real magnetic materials that are on the verge of quantum disorder.

###### Acknowledgements.

The authors thank Fakher Assaad, Subir Sachdev, Yin-Chen He, Max Metlitski, T. Senthil, Chong Wang and Stefan Wessel for helpful discussions. NSM, GYS and ZYM are supported by the Ministry of Science and Technology of China under Grant No. 2016YFA0300502, the key research program of the Chinese Academy of Sciences under Grant No. XDPB0803, and the National Science Foundation of China under Grant Nos. 11421092, 11574359 and 11674370, as well as the National Thousand-Young Talents Program of China. NSM would like to thank Boston University for support under its Condensed Matter Theory Visitors Program. AV and YZY are supported by a Simons Investigator Grant. AWS is supported by the NSF under Grant No. DMR-1710170 and by a Simons Investigator Grant. CX is supported by the David and Lucile Packard Foundation and NSF Grant No. DMR-1151208. We thank the following institutions for allocation of CPU time: the Center for Quantum Simulation Sciences at the Institute of Physics, Chinese Academy of Sciences, and the Tianhe-1A platform at the National Supercomputer Center in Tianjin.## References

- Senthil et al. (2004a) T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Science 303, 1490 (2004a).
- Senthil et al. (2004b) T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. Fisher, Phys. Rev. B 70, 144407 (2004b).
- Motrunich and Vishwanath (2004) O. I. Motrunich and A. Vishwanath, Phys. Rev. B 70, 075104 (2004).
- Hermele et al. (2005) M. Hermele, T. Senthil, and M. P. A. Fisher, Phys. Rev. B 72, 104404 (2005).
- Ran and Wen (2006) Y. Ran and X.-G. Wen, eprint arXiv:cond-mat/0609620 (2006), cond-mat/0609620 .
- Senthil and Fisher (2006) T. Senthil and M. P. A. Fisher, Phys. Rev. B 74, 064405 (2006).
- Grover and Vishwanath (2013) T. Grover and A. Vishwanath, Phys. Rev. B 87, 045129 (2013).
- Lu and Lee (2014) Y.-M. Lu and D.-H. Lee, Phys. Rev. B 89, 195143 (2014).
- Xu and You (2015) C. Xu and Y.-Z. You, Phys. Rev. B 92, 220416 (2015).
- Karch and Tong (2016) A. Karch and D. Tong, Phys. Rev. X 6, 031043 (2016).
- Hsin and Seiberg (2016) P.-S. Hsin and N. Seiberg, Journal of High Energy Physics 2016, 95 (2016).
- Cheng and Xu (2016) M. Cheng and C. Xu, Phys. Rev. B 94, 214415 (2016).
- Wang et al. (2017) C. Wang, A. Nahum, M. A. Metlitski, C. Xu, and T. Senthil, Phys. Rev. X 7, 031051 (2017).
- Metlitski and Thorngren (2017) M. A. Metlitski and R. Thorngren, ArXiv e-prints (2017), arXiv:1707.07686 [cond-mat.str-el] .
- You et al. (2017) Y.-Z. You, Y.-C. He, A. Vishwanath, and C. Xu, ArXiv e-prints (2017), arXiv:1711.00863 [cond-mat.str-el] .
- Chubukov et al. (1994) A. V. Chubukov, S. Sachdev, and J. Ye, Phys. Rev. B 49, 11919 (1994).
- Kaul et al. (2013) R. K. Kaul, R. G. Melko, and A. W. Sandvik, Annual Review of Condensed Matter Physics 4, 179 (2013).
- Sandvik et al. (2002) A. W. Sandvik, S. Daul, R. R. P. Singh, and D. J. Scalapino, Phys. Rev. Lett. 89, 247201 (2002).
- Sandvik (2007) A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
- Lou et al. (2009) J. Lou, A. W. Sandvik, and N. Kawashima, Phys. Rev. B 80, 180414 (2009).
- Sandvik (2010a) A. W. Sandvik, Phys. Rev. Lett. 104, 177201 (2010a).
- Pujari et al. (2013) S. Pujari, K. Damle, and F. Alet, Phys. Rev. Lett. 111, 087203 (2013).
- Block et al. (2013) M. S. Block, R. G. Melko, and R. K. Kaul, Phys. Rev. Lett. 111, 137202 (2013).
- Shao et al. (2016) H. Shao, W. Guo, and A. W. Sandvik, Science 352, 213 (2016).
- D’Emidio and Kaul (2017) J. D’Emidio and R. K. Kaul, Phys. Rev. Lett. 118, 187202 (2017).
- Qin et al. (2017a) Y. Q. Qin, Y.-Y. He, Y.-Z. You, Z.-Y. Lu, A. Sen, A. W. Sandvik, C. Xu, and Z. Y. Meng, Phys. Rev. X 7, 031052 (2017a).
- Zhang et al. (2017) X.-F. Zhang, Y.-C. He, S. Eggert, R. Moessner, and F. Pollmann, ArXiv e-prints (2017), arXiv:1706.05414 [cond-mat.str-el] .
- Nahum et al. (2011) A. Nahum, J. T. Chalker, P. Serna, M. Ortuño, and A. M. Somoza, Phys. Rev. Lett. 107, 110601 (2011).
- Nahum et al. (2015a) A. Nahum, J. T. Chalker, P. Serna, M. Ortuño, and A. M. Somoza, Phys. Rev. X 5, 041048 (2015a).
- Nahum et al. (2015b) A. Nahum, P. Serna, J. T. Chalker, M. Ortuño, and A. M. Somoza, Phys. Rev. Lett. 115, 267203 (2015b).
- Tang and Sandvik (2013) Y. Tang and A. W. Sandvik, Phys. Rev. Lett. 110, 217213 (2013).
- Sandvik et al. (2011) A. W. Sandvik, V. N. Kotov, and O. P. Sushkov, Phys. Rev. Lett. 106, 207203 (2011).
- Assaad and Grover (2016) F. F. Assaad and T. Grover, Phys. Rev. X 6, 041049 (2016).
- Han et al. (2012) T. H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Nature 492, 406 (2012).
- Fu et al. (2015) M. Fu, T. Imai, T.-H. Han, and Y. S. Lee, Science 350, 655 (2015).
- Feng et al. (2017) Z. Feng, Z. Li, X. Meng, W. Yi, Y. Wei, J. Zhang, Y.-C. Wang, W. Jiang, Z. Liu, S. Li, F. Liu, J. Luo, S. Li, G. qing Zheng, Z. Y. Meng, J.-W. Mei, and Y. Shi, Chinese Physics Letters 34, 077502 (2017).
- Wen (2017) X.-G. Wen, Chinese Physics Letters 34, 90101 (2017).
- Wei et al. (2017) Y. Wei, Z. Feng, W. Lohstroh, C. dela Cruz, W. Yi, Z. F. Ding, J. Zhang, C. Tan, L. Shu, Y.-C. Wang, J. Luo, J.-W. Mei, Z. Y. Meng, Y. Shi, and S. Li, ArXiv e-prints (2017), arXiv:1710.02991 [cond-mat.str-el] .
- Feng et al. (2017) Z. Feng, Y. Wei, R. Liu, D. Yan, Y.-C. Wang, J. Luo, A. Senyshyn, C. dela Cruz, W. Yi, J.-W. Mei, Z. Y. Meng, Y. Shi, and S. Li, ArXiv e-prints (2017), arXiv:1712.06732 [cond-mat.str-el] .
- Sen and Sandvik (2010) A. Sen and A. W. Sandvik, Phys. Rev. B 82, 174428 (2010).
- Guida and Zinn-Justin (1998) R. Guida and J. Zinn-Justin, Journal of Physics A: Mathematical and General 31, 8103 (1998).
- Campostrini et al. (2001) M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev. B 63, 214503 (2001).
- Sandvik (2010b) A. W. Sandvik, AIP Conf. Proc. 1297, 135 (2010b).
- Sandvik (1998) A. W. Sandvik, Phys. Rev. B 57, 10287 (1998).
- Beach (2004) K. S. D. Beach, eprint arXiv:cond-mat/0403055 (2004), cond-mat/0403055 .
- Syljuåsen (2008) O. F. Syljuåsen, Phys. Rev. B 78, 174429 (2008).
- Sandvik (2016) A. W. Sandvik, Phys. Rev. E 94, 063308 (2016).
- Qin et al. (2017b) Y. Q. Qin, B. Normand, A. W. Sandvik, and Z. Y. Meng, Phys. Rev. Lett. 118, 147207 (2017b).
- (49) H. Shao and A. W. Sandvik, “Parametrization of the spectrum in stochastic analytic continuation,” In preparation.
- Shao et al. (2017) H. Shao, Y. Q. Qin, S. Capponi, S. Chesi, Z. Y. Meng, and A. W. Sandvik, Phys. Rev. X 7, 041072 (2017).
- Huang et al. (2017) C.-J. Huang, Y. Deng, Y. Wan, and Z. Y. Meng, ArXiv e-prints (2017), arXiv:1707.00099 [cond-mat.str-el] .
- Shu et al. (2017) Y.-R. Shu, M. Dupont, D.-X. Yao, S. Capponi, and A. W. Sandvik, ArXiv e-prints (2017), arXiv:1712.01701 [cond-mat.str-el] .
- Sandvik (1999) A. W. Sandvik, Phys. Rev. B 59, R14157 (1999).
- Lohöfer et al. (2015) M. Lohöfer, T. Coletta, D. G. Joshi, F. F. Assaad, M. Vojta, S. Wessel, and F. Mila, Phys. Rev. B 92, 245137 (2015).
- Wen (2004) X.-G. Wen, Quantum Field Theory of Many-body Systems : From the Origin of Sound to an Origin of Light and Electrons (Oxford Graduate Texts, 2004) Chap. 9.
- (56) Although such a Grose-Neveu critical point is different from the DQCP, we only use it to provide a rough estimate of the spectral features close to a magnetic ordered phase. We do not claim that the criticality of DQCP can be correctly understood by our mean field + RPA approach.
- Karch et al. (2017) A. Karch, B. Robinson, and D. Tong, Journal of High Energy Physics 2017, 17 (2017).
- Potter et al. (2016) A. C. Potter, C. Wang, M. A. Metlitski, and A. Vishwanath, ArXiv e-prints (2016), arXiv:1609.08618 [cond-mat.str-el] .
- Mross et al. (2016) D. F. Mross, J. Alicea, and O. I. Motrunich, Phys. Rev. Lett. 117, 016802 (2016).
- Slagle et al. (2015) K. Slagle, Y.-Z. You, and C. Xu, Phys. Rev. B 91, 115121 (2015).
- He et al. (2016) Y.-Y. He, H.-Q. Wu, Y.-Z. You, C. Xu, Z. Y. Meng, and Z.-Y. Lu, Phys. Rev. B 93, 115150 (2016).
- Wu et al. (2016) H.-Q. Wu, Y.-Y. He, Y.-Z. You, T. Yoshida, N. Kawakami, C. Xu, Z. Y. Meng, and Z.-Y. Lu, Phys. Rev. B 94, 165121 (2016).
- Hasenbusch and Török (1999) M. Hasenbusch and T. Török, Journal of Physics A: Mathematical and General 32, 6361 (1999).
- Meng and Wessel (2008) Z. Y. Meng and S. Wessel, Phys. Rev. B 78, 224416 (2008).

## Appendix A Quantum critical point of the EPJJ model

In this section, we discuss how the critical point of the EPJJ model is determined from finite-size scaling of QMC results. The physical observable used here is the spin stiffness

(19) |

where is the free energy and is the twisting angle between spins in two different columns.

is easy to measure in SSE QMC simulations Sandvik (2010b) and it has well defined finite size scaling form

(20) |

with the dynamic exponent and the correlation length exponent , and is the scaling function. The quantum phase transition from AFXY phase to the columnar phase in the EPJJ model belongs to D O(2) universality class Hasenbusch and Török (1999); Meng and Wessel (2008). In the simulation we fix such that the is a constant and we effectively have a single-parameter scaling function. Then the quantity is dimensionless and does not change with at critical point.

In the QMC simulation, is evaluated using winding number fluctuations as Sandvik (2010b)

(21) |

where and are the number of all operators that transport spin in positive and negative directions along the lattice direction, respectively, during propagation of the spin state in imaginary time. Since the EPJJ model is spatially anisotropic, one can calculate along the or directions of the lattice, and both of them satisfy the same finite size scaling form Eq. (20) with different scaling functions. We label the two quantities as and , and for the sake of simplicity, only show here.

Fig. 7 (a) depicts the raw data for four different system sizes, . In Fig. 7(b), we plot the scaled quantity against the control parameter. The the crossing point of the curves should drift toward the critical point , and from our data we obtain . Finally, in Fig. 7(c), we further rescale the -axis as , with and , thus obtaining the scaling function in Eq. (20) as the common curve onto which the data for different system sizes collapse. The good data collapse without any adjustable parameters supports the expected D O(2) universality class.

## Appendix B Spectra of EPJJ model

In this section, we provide more detailed information of the EPJJ spectra inside the AFXY phase and close to the 3DXY transition point. Fig. 8 shows a scan through the square-lattice Brilliune zone (BZ) with more points, along the path . In Fig. 8(a), is shown at . Here the left part is identical to Fig. 4(a), where the Goldstone mode at is seen, and the spectra at is gapped. The right part is slightly different, with the spectra at also gapped, but, due to the folding of the BZ coming from the doubling of the real space unit cell, i.e., twice along the direction of the square lattice, a “shadow” band originating from the gapless dispersion from of the AFXY phase, presents itself from . Fig. 8 (c) shows along the the same path. Since the excitations are gapped in the AFXY phase, except for the point, there is no obvious sign of the band folding close to .

Fig. 8(b),(d) show and close to the 3DXY critical point at . The Goldstone mode at still presents itself, as well as the band folding close to in the channel. In the channel, on the other hand, all the spectra are gapped except , and above the gap, some signature of the band-folding can also be seen close to .