Dynamics of Compact Quantum Electrodynamics at Large Fermion Flavor

TL;DR

This study uses advanced quantum Monte Carlo simulations to explore the dynamic properties of large-flavor compact QED in (2+1)D, revealing phase transitions and spectral features that support the existence of a U(1) spin liquid and shed light on the phase transition's universality class.

Contribution

It provides the first detailed numerical investigation of the dynamical behaviors across the deconfinement-to-confinement transition at large fermion flavor in compact QED.

Findings

Identification of continua in spectral functions in the deconfined phase.

Observation of gapped spectra and translational symmetry breaking in the confined phase.

Evidence supporting the U(1) deconfined phase and insights into the phase transition nature.

Abstract

Thanks to the development in quantum Monte Carlo technique, the compact U(1) lattice gauge theory coupled to fermionic matter at (2+1)D is now accessible with large-scale numerical simulations, and the ground state phase diagram as a function of fermion flavor ($N_f$) and the strength of gauge fluctuations is mapped out~\cite{Xiao2018Monte}. Here we focus on the large fermion flavor case ($N_f=8$) to investigate the dynamic properties across the deconfinement-to-confinement phase transition. In the deconfined phase, fermions coupled to the fluctuating gauge field to form U(1) spin liquid with continua in both spin and dimer spectral functions, and in the confined phase fermions are gapped out into valence bond solid phase with translational symmetry breaking and gapped spectra. The dynamical behaviors provide supporting evidence for the existence of the U(1) deconfined phase and could shine light on the nature of the U(1)-to-VBS phase transition which is of the QED$_3$-Gross-Neveu chiral O(2) universality whose properties still largely unknown.

Figures

Tables

Table 1. Table 1: Emergent symmetry and conserved currents for the N f = 8 subscript 𝑁 𝑓 8 N_{f}=8 lattice QED model.

Monopole operator	irrelevant	relevant
Emergent symmetry	$\mathrm{PSU}\left(16\right)$	$\mathrm{PSU}\left(8\right)\times \mathrm{PSU}\left(8\right)\times {\mathbb{Z}}_4$
Conserved currents	$j_{i0}^{\mu },j_{0a}^{\mu },j_{ia}^{\mu }$	$j_{0a}^{\mu },j_{3a}^{\mu }$
Number of currents	255	126

Table 2. Table 2: Identification of operators at high symmetry points.

$𝒒$	$(0,0)$	$(\pi ,0)$	$(0,\pi )$	$(\pi ,\pi )$
$S_{𝒒}^a$	$-j_{0a}^0$	$j_{1a}^2$	$j_{2a}^1$	Néel
$D_{𝒒}^x$	$ℒ-{𝒯}^{11}$	$x$-VBS	$-i\overline{c}{\gamma }^0{\tau }^2{\partial }_{x}c$	$j_{30}^2$
$D_{𝒒}^y$	$ℒ-{𝒯}^{22}$	$-i\overline{c}{\gamma }^0{\tau }^1{\partial }_{y}c$	$y$-VBS	$-j_{30}^1$

Equations

$$\begin{array}{rl}H=& \frac{1}{2}J{N}_f\sum _{⟨i,j⟩}\frac{1}{4}{\hat{L}}_{ij}^2-t\sum _{⟨i,j⟩\alpha }\left(c_{i\alpha }^{†}{e}^{\mathrm{i}{\hat{\varphi }}_{ij}}{c}_{j\alpha }+\text{h.c.}\right)\\ & +\frac{1}{2}K{N}_f\sum _{\square }\cos \left(\text{curl}\hat{\varphi }\right),\end{array}$$

$$Z=\mathrm{Tr}{e}^{-\beta H}=\sum _{[\varphi ]}{\omega }_B[\varphi ]{\omega }_F[\varphi ],$$

$$\begin{split}\omega_{B}[\phi]=\prod_{\tau=M}^{1}&e^{-\frac{4}{JN_{f}\Delta\tau}\sum_{\langle ij\rangle}\Big{(}1-\cos\big{(}\phi_{ij}(\tau+1)-\phi_{ij}(\tau)\big{)}\Big{)}}\\ &e^{-\frac{1}{2}KN_{f}\Delta\tau\cos\big{(}\sum_{\langle i,j\rangle\in\Box}\phi_{i,j}(\tau)\big{)}},\end{split}$$

$${\omega }_F[\varphi ]={\left(\det \left[I+B^M{B}^{M-1}\cdots B^{\tau }\cdots B^2{B}^1\right]\right)}^{N_f}.$$

$$\begin{array}{rl}{B}^{\tau }& =\prod _{⟨ij⟩}\exp \left(\Delta \tau \left[\begin{array}{cccccc}0& & & & & \\ & 0& & e^{\mathrm{i}{\varphi }_{ij}(\tau )}& & \\ & & 0& & & \\ & e^{-\mathrm{i}{\varphi }_{ij}(\tau )}& & 0& & \\ & & & & 0& \end{array}\right]\right)\\ & =\prod _{⟨ij⟩}\left[\begin{array}{cccccc}1& & & & & \\ & \cosh \Delta \tau & 0& e^{\mathrm{i}{\varphi }_{ij}}\sinh \Delta \tau & & \\ & 0& 1& 0& & \\ & e^{-\mathrm{i}{\varphi }_{ij}}\sinh \Delta \tau & 0& \cosh \Delta \tau & & \\ & & & & 1& \end{array}\right].\end{array}$$

$$r=\frac{{\omega }_B[{\varphi }^{\prime }]}{{\omega }_B[\varphi ]}\cdot \frac{{\omega }_F[{\varphi }^{\prime }]}{{\omega }_F[\varphi ]}.$$

$$\displaystyle\begin{split}\frac{\omega_{B}[\phi^{\prime}]}{\omega_{B}[\phi]}=&\frac{\exp\left({\frac{4}{JN_{f}\Delta\tau}\left(\cos{\left(\phi_{ij}\left(\tau+1\right)-\phi^{\prime}_{ij}\left(\tau\right)\right)}+\cos\left(\phi_{ij}\left(\tau-1\right)-\phi^{\prime}_{ij}\left(\tau\right)\right)\right)}\right)}{\exp\left({\frac{4}{JN_{f}\Delta\tau}\left(\cos{\left(\phi_{ij}\left(\tau+1\right)-\phi_{ij}\left(\tau\right)\right)}+\cos\left(\phi_{ij}\left(\tau-1\right)-\phi_{ij}\left(\tau\right)\right)\right)}\right)}\cdot\frac{\exp\left(-\frac{1}{2}\Delta\tau KN_{f}\cos\left(\sum_{\langle ij\rangle\in\Box}\phi^{\prime}_{i,j}(\tau)\right)\right)}{\exp\left(-\frac{1}{2}\Delta\tau KN_{f}\cos\left(\sum_{\langle ij\rangle\in\Box}\phi_{i,j}(\tau)\right)\right)}\\ =&\exp\Bigg{[}-\frac{16}{JN_{f}\Delta\tau}\cdot\cos\left(\frac{\phi_{ij}(\tau+1)-\phi_{ij}(\tau-1)}{2}\right)\cdot\sin\left(\frac{\phi_{ij}(\tau+1)+\phi_{ij}(\tau-1)-\phi^{\prime}_{ij}(\tau)-\phi_{ij}(\tau)}{2}\right)\cdot\\ &\sin\left(\frac{\phi_{ij}(\tau)-\phi^{\prime}_{ij}(\tau)}{2}\right)-\frac{\Delta\tau KN_{f}}{2}\left(\cos\left(\sum_{\langle ij\rangle\in\Box}\phi^{\prime}_{ij}(\tau)\right)-\cos\left(\sum_{\langle ij\rangle\in\Box}\phi_{ij}\left(\tau\right)\right)\right)\Bigg{]},\end{split}$$

$$\begin{array}{rl}\frac{{\omega }_F[{\varphi }^{\prime }]}{{\omega }_F[\varphi ]}& =\frac{\det \left[I+B(\beta ,\tau )(1+\Delta )B(\tau ,0)\right]}{\det \left[I+B(\beta ,\tau )B(\tau ,0)\right]}\\ & =\det \{1+\Delta \left[1-G(\tau ,\tau )\right]\}.\end{array}$$

$$G^{\prime }(\tau ,\tau )=G(\tau ,\tau ){\left[1+\Delta (1-G(\tau ,\tau ))\right]}^{-1},$$

$$\Delta =\left[\begin{array}{cc}{\Delta }_{ii}& {\Delta }_{ij}\\ {\Delta }_{ji}& {\Delta }_{jj}\end{array}\right],$$

$$\begin{array}{rl}& {\Delta }_{ii}=1-e^{-\mathrm{i}\left({\varphi }_{ij}-{\varphi }_{ij}^{\prime }\right)}{\sinh }^2\Delta \tau ,\\ & {\Delta }_{jj}=1-e^{\mathrm{i}\left({\varphi }_{ij}-{\varphi }_{ij}^{\prime }\right)}{\sinh }^2\Delta \tau ,\\ & {\Delta }_{ij}=\left(-e^{\mathrm{i}{\varphi }_{ij}}+e^{\mathrm{i}{\varphi }_{ij}^{\prime }}\right)\sinh \Delta \tau \cosh \Delta \tau ,\\ & {\Delta }_{ji}=\left(-e^{-\mathrm{i}{\varphi }_{ij}}+e^{-\mathrm{i}{\varphi }_{ij}^{\prime }}\right)\sinh \Delta \tau \cosh \Delta \tau .\end{array}$$

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

$$H_{\text{MF}}=t\sum _i\mathrm{i}(c_{i+\hat{x}}^{†}{c}_i+(-{)}^x{c}_{i+\hat{y}}^{†}{c}_i)+h.c.$$

$$H_{\text{MF}}=-2t\sum _{\mathbf{k}}{c}_{\mathbf{k}}^{†}(\sin k_x{\sigma }^{10}+\sin k_y{\sigma }^{31})\otimes \text{mathds}{1}_{8\times 8}{c}_{\mathbf{k}},$$

$$\mathcal{L}=\bar{c}(\mathrm{i}{\gamma }^{\mu }{\partial }_{\mu })c,$$

$$S_i^a=\frac{1}{2}{c}_i^{†}{T}^a{c}_i,$$

$$\begin{array}{rl}{D}_i^x& =-(c_{i+\hat{x}}^{†}{e}^{\mathrm{i}{\varphi }_{i+\hat{x},i}}{c}_i+h.c.)/2\\ & =(\mathrm{i}{c}_{i+\hat{x}}^{†}{c}_i+h.c.)/2,\\ D_i^y& =-(c_{i+\hat{y}}^{†}{e}^{\mathrm{i}{\varphi }_{i+\hat{y},i}}{c}_i+h.c.)/2\\ & =(-{)}^x(\mathrm{i}{c}_{i+\hat{y}}^{†}{c}_i+h.c.)/2,\end{array}$$

$${\mathcal{O}}_{\mathbf{q}}=\sum _{\mathbf{k}}{c}_{\mathbf{k}}^{†}{v}_{\mathbf{q}}{c}_{\mathbf{k}+\mathbf{q}},$$

$$\begin{array}{rl}{S}_{\mathbf{q}}^a& =\frac{1}{2}\sum _{\mathbf{k}}{c}_{\mathbf{k}}^{†}{\sigma }^{00}{T}^a{c}_{\mathbf{k}+\mathbf{q}},\\ D_{\mathbf{q}}^x& =\sum _{\mathbf{k}}{c}_{\mathbf{k}}^{†}(\sin (k_x+\frac{q_x}{2})e^{-\mathrm{i}{q}_x/2}{\sigma }^{10}\text{mathds}1)c_{\mathbf{k}+\mathbf{q}},\\ D_{\mathbf{q}}^y& =\sum _{\mathbf{k}}{c}_{\mathbf{k}}^{†}(\sin (k_y+\frac{q_y}{2})e^{-\mathrm{i}{q}_y/2}{\sigma }^{31}\text{mathds}1)c_{\mathbf{k}+\mathbf{q}}.\end{array}$$

$${\chi }_{\mathcal{O}}(\mathbf{q},\omega )=\int dt{e}^{\mathrm{i}\omega t}⟨{\mathcal{O}}_{-\mathbf{q}}(t){\mathcal{O}}_{\mathbf{q}}(0)⟩,$$

$$\begin{array}{rl}{A}_{\mathcal{O}}(\mathbf{q},\omega )& =\sum _{m,n,\mathbf{k}}⟨m,\mathbf{k}|v_{-\mathbf{q}}|n,\mathbf{k}+\mathbf{q}⟩⟨n,\mathbf{k}+\mathbf{q}|v_{\mathbf{q}}|m,\mathbf{k}⟩\\ & \delta (\omega +{\xi }_{m,\mathbf{k}}-{\xi }_{n,\mathbf{k}+\mathbf{q}})(\Theta ({\xi }_{m,\mathbf{k}})-\Theta ({\xi }_{n,\mathbf{k}+\mathbf{q}}))\end{array}$$

$$\mathcal{L}=\bar{c}{\gamma }^{\mu }(\mathrm{i}{\partial }_{\mu }-a_{\mu })c+g(\mathcal{M}+{\mathcal{M}}^{†})+\cdots ,$$

$$\mathrm{SU}(8):c\to e^{\mathrm{i}{\theta }_a{T}^a}c,{\mathbb{Z}}_4:c\to e^{\mathrm{i}\frac{\pi }{4}{\tau }^3}c,$$

$$m=\frac{1}{2}{c}^{†}{\tau }^{3}c=0,\pm 1,\pm 2,\pm 3,\pm 4,$$

$$j_{i0}^{\mu }=\bar{c}{\gamma }^{\mu }{\tau }^{i}c,j_{0a}^{\mu }=\bar{c}{\gamma }^{\mu }{T}^{a}c,j_{ia}^{\mu }=\bar{c}{\gamma }^{\mu }{\tau }^i{T}^{a}c.$$

$$\frac{\mathrm{SU}(8)}{{\mathbb{Z}}_8}\times \frac{\mathrm{SU}(8)}{{\mathbb{Z}}_8}\times \frac{{\mathbb{Z}}_8}{{\mathbb{Z}}_2}=\mathrm{PSU}(8)\times \mathrm{PSU}(8)\times {\mathbb{Z}}_4.$$

$$j_{0a}^{\mu }=\bar{c}{\gamma }^{\mu }{T}^{a}c,j_{3a}^{\mu }=\bar{c}{\gamma }^{\mu }{\tau }^3{T}^{a}c.$$

$$\begin{array}{rl}{S}_{(0,0)}^a& \sim c^{†}{\sigma }^{00}{T}^{a}c=-\bar{c}{\gamma }^0{T}^{a}c=-j_{0a}^0,\\ S_{(\pi ,0)}^a& \sim c^{†}{\sigma }^{30}{T}^{a}c=\bar{c}{\gamma }^2{\tau }^1{T}^{a}c=j_{1a}^2,\\ S_{(0,\pi )}^a& \sim c^{†}{\sigma }^{03}{T}^{a}c=\bar{c}{\gamma }^1{\tau }^2{T}^{a}c=j_{2a}^1,\\ D_{(\pi ,\pi )}^x& \sim -c^{†}{\sigma }^{23}\text{mathds}1c=\bar{c}{\gamma }^2{\tau }^{3}c=j_{30}^2,\\ D_{(\pi ,\pi )}^y& \sim -c^{†}{\sigma }^{02}\text{mathds}1c=-\bar{c}{\gamma }^1{\tau }^{3}c=-j_{30}^1.\end{array}$$

$${\mathcal{T}}^{\mu \nu }=\frac{\partial \mathcal{L}}{\partial ({\partial }_{\mu }c)}{\partial }^{\nu }c-{\delta }^{\mu \nu }\mathcal{L}=\mathrm{i}\bar{c}({\gamma }^{\mu }{\partial }^{\nu }-{\delta }^{\mu \nu }{\gamma }^{\lambda }{\partial }_{\lambda })c.$$

$$\begin{array}{rl}{D}_{(0,0)}^x& \sim -\mathrm{i}{c}^{†}{\sigma }^{10}\text{mathds}1{\partial }_{x}c=-\mathrm{i}\bar{c}{\gamma }^1{\partial }_{x}c=\mathcal{L}-{\mathcal{T}}^{11},\\ D_{(0,0)}^y& \sim -\mathrm{i}{c}^{†}{\sigma }^{31}\text{mathds}1{\partial }_{y}c=-\mathrm{i}\bar{c}{\gamma }^2{\partial }_{y}c=\mathcal{L}-{\mathcal{T}}^{22}.\end{array}$$

Full text

I Introduction

For several decades, the topic of dynamical coupling between lattice gauge fields and fermionic matter fields have attracted considerable attention among physicists from high-energy Fiebig and Woloshyn (1990); Armour et al. (2011); Karthik and Narayanan (2016, 2018a, 2018b); Braun et al. (2014); Kotikov et al. (2016) and condensed matter Affleck and Marston (1988); Marston and Affleck (1989); Lee et al. (1998); Rantner and Wen (2001, ); Senthil and Fisher (2000); Senthil and Motrunich (2002); Herbut and Seradjeh (2003); Hermele et al. (2004a); Assaad (2005); Hermele et al. (2005, 2007); Nogueira and Kleinert (2008); He et al. (2017); Assaad and Grover (2016); Gazit et al. (2017, 2018); Prosko et al. (2017) communities. Previous works have quickly established the understanding at the large fermion flavor $(N_{f})$ limit Affleck and Marston (1988); Marston and Affleck (1989); Lee et al. (1998); Rantner and Wen (2001, ); Senthil and Fisher (2000); Senthil and Motrunich (2002); Herbut and Seradjeh (2003); Hermele et al. (2004a), as $1/N_{f}$ expansion is controlled thence, but left the physically most interesting cases of small $N_{f}$ – for example $N_{f}=2$ corresponds to the spin-1/2 case of electrons – unsolved. The very recent breakthrough of quantum Monte Carlo (QMC) simulations of Z2 Assaad and Grover (2016); Gazit et al. (2017, 2018); Chen et al. (2019) and U(1) Xu et al. (2019) gauge fields coupled to fermions provide the possibility of concrete investigations at the small $N_{f}$, and the expected deconfinement-to-confinement phase transitions and special properties of these phases are discovered. In such settings, the interactions between fermions are mediated via the fluctuating gauge bosons, which resemble the situation of fractionalized particles and emergent gauge fields in several prototypical strongly correlated systems including, but not limited to, the low-energy description of the high-temperature superconductors Lee et al. (1998); Senthil and Motrunich (2002); Hermele et al. (2005), frustrated magnets Hermele et al. (2005); Balents et al. (2002); Hermele et al. (2004b) and deconfined quantum criticalities Senthil et al. (2004); Sandvik (2007); Qin et al. (2017a); Ma et al. (2018); Li et al. (2019); Zhou et al. (2019). The quantum phases and phase transitions discovered are clearly beyond the Landau-Ginzburg-Wilson paradigm built upon the concepts of symmetry-breaking and local order parameters, and served as the building bulks of the new paradigm of quantum matter.

As for the discrete Z2 gauge field coupled to fermionic matter in (2+1)D Assaad and Grover (2016); Gazit et al. (2017, 2018), the deconfined phase with fractionalized fermionic excitations at weak gauge fluctuation and confined phase with symmetry-breaking at strong gauge fluctuation have been revealed. The (Z2) deconfinement-to-confinement transitions are continuous and associated with fermion gap opening in the excitation spectrum. Further developments that involve not only Z2 gauge but also Z2 matter fields to dynamically couple to Fermi surface (FS) give rise to the long-thought orthogonal metal phase which has metallic transport but no quasiparticles at the FS Chen et al. (2019), probably the simplest non-Fermi-liquid that can be generated without ambiguity Nandkishore et al. (2012) in (2+1)D lattice models.

As for the continuous U(1) gauge field coupled to fermionic matter at (2+1)D, such as the compact quantum electrodynamics (cQED3), there are fundamental physical questions awaiting for affirmative answer. The pure gauge theory at (2+1)D is known to be always confined Herbut and Seradjeh (2003); Polyakov (1977); Mandelstam (1976); Case et al. (2004), but whether the coupling to gapless fermionic matter could drive the system towards deconfinement have been debated Herbut and Seradjeh (2003); Polyakov (1977); Mandelstam (1976); Case et al. (2004); Song et al. (2018a). As mentioned above, large $N_{f}$ limits Hermele et al. (2004a, 2007); Nogueira and Kleinert (2008); Mithat (2008) demonstrated the existence of the U(1) deconfined phase, but previous QMC works at medium and small $N_{f}$ values are shown inconclusive Armour et al. (2011); Karthik and Narayanan (2016, 2018a, 2018b) due to the difficulties in effectively simulating the continuous gauge fields in the $(2+1)$ space-time with zero modes at the fermion spectra. It is only till very recent, that in Ref. Xu et al. (2019), with the help of fast updates and high level parallelization, that the phase diagram of U(1) gauge field coupled to fermion field at small $N_{f}$ has been mapped out, and the existence of U(1) deconfined phase, or algebraic quantum spin liquid in the condensed matter parlance Hermele et al. (2005), for $N_{f}=2,4,6,8$ are discovered with certainty. However, even the latest QMC simulations are still by no means easy, suffering long autocorrelation time in the critical phase of U(1) deconfined spin liquid, the transitions from U(1) deconfined phase to various confined phases (antiferromagnetic insulator (AFM) phase for $N_{f}=2$, and valence bond solid (VBS) phase and AFM for $N_{f}=4,6,8$) have not been able to investigated in detail, both statically and dynamically.

These transitions, dubbed QED3-Gross-Neveu chiral Heisenberg (QED3-GN-O(3)) or XY (QED3-GN-O(2)) transitions, are of high interests to both condensed matter and high-energy physicists, as the phase transition of algebraic quantum spin liquid to other magnetically order phases have experimental relevance. Inspired by the numerical work of Ref. Xu et al. (2019), there are recently several analytical works addressing the critical properties of them Gracey (2018); Ihrig et al. (2018); Zerf et al. (2019); Dupuis et al. (2019); Zerf et al. (2018); Boyack et al. (2019). And the conclusions drawn there are that the U(1)-to-AFM and U(1)-to-VBS phase transitions are indeed possible and the higher-order perturbative RG calculations performed also suggest the possible range of critical exponents of these QED3-GN transitions within $1/N_{f}$ and $4-\epsilon$ expansions Gracey (2018); Ihrig et al. (2018); Zerf et al. (2019); Dupuis et al. (2019); Zerf et al. (2018); Boyack et al. (2019).

While the QMC evaluation of the critical exponents are still difficult (currently the largest system accessed are $L=20$ due to the aforementioned computational complexity), the dynamical signatures of the transition would then provide guiding evidence for comprehensive understanding of them. Similar as the case of deconfined quantum critical point with emergent O(4) symmetry, where the coupling effects of fractionalized spinon and emergent U(1) gauge fields manifest in the spin spectral functions Ma et al. (2018), the unearthness of the QED3-GN dynamical signatures will provide similar physical understanding. In terms of quantum Monte Carlo simulations, the dynamical signatures can be obtained in two steps. One first measures the imaginary time correlation functions with good statistics, then performs analytic continuation to convert the correlation from imaginary to real frequencies. The recent developments of stochastic analytic continuation (SAC) scheme Sandvik (2016) is proven to be more reliable and could reveal non-trivial results in both unfrustrated and frustrated magnetic systems in 2D and 3D Ma et al. (2018); Shao et al. (2017); Huang et al. (2018); Qin et al. (2017b); Sun et al. (2018). Therefore the techniques for investigating the dynamical properties of the QED3-GN transitions are available.

Aware of the high interests and great difficulty in studying the U(1) deconfined to confined transition, and with the help of the state-of-art QMC methodology and SAC machinery, in this work, we nevertheless take the first step to investigate the dynamical signature of the QED3-GN transition at a large but finite fermion flavor of $N_{f}=8$. As a function of the strength of the gauge fluctuations, the deconfined phase (the algebraic quantum spin liquid) and the confined phase (VBS) are investigated in detail, and the dynamical sigature of their transition in the form of the spin and dimer spectra with continua are discovered. The physical meaning of such continua inside the U(1) deconfined phase and at the transition are addressed as well. These results set the stage for the further investigations of the smaller and physically more relevant $N_{f}$ and can be used to guide the experimental detection in inelastic neutron scattering and nuclear magnetic resonance for condensed matter materials which could host fractionalized excitations coupled with emergent gauge structures Feng et al. (2017); Wei et al. (2017). For example, the observation of the conserved current correlations in the spin and dimer spectra Ma et al. (2019); Lee et al. (2019); Huang et al. (2019) would be the decisive evidence for the deconfinement and emergent gauge fields.

With these thoughts in mind, we organize the rest of the paper as follows. In Sec. II, a quantum rotor model that describes the setting of compact QED3 coupled to fermions is introduced. The QMC and SAC methods employed to solve the model are also explained in a concise manner. The analysis of the theoretical interpretation of the continua and associated symmetry properties of the U(1) deconfined phase are given in Sec. III. In Sec. IV, QMC numerical results including the net gauge flux and most importantly, the spin and dimer spectra are presented, with the physical meaning of the continua therein thoroughly discussed, which serve as the dynamical signature of the U(1) phase and U(1)-to-VBS QED3-GN transition. Conclusion and outlook are presented in Sec. V.

II Model and Method

In this work we study a 2D quantum rotor model coupled to fermions considered in Ref. Xu et al. (2019), whose Hamiltonian can be written as

[TABLE]

where $c_{i\alpha}$($c_{i\alpha}^{\dagger}$) is the annihilation (creation) operator for a fermion with fermion flavor $\alpha$. The $\alpha$ runs from 1 to $N_{f}$ and here we focus on the case of $N_{f}=8$. As shown in Fig. 1 (b), the nearest hopping of fermions is associated with a phase $\phi_{ij}$, this phase inserts magnetic flux through each plaquette. The flux term with $K>0$ favors $\pi$-flux in each elementary plaquette $\Box$. Following the convention in Ref. Xu et al. (2019), we fixed $t=1$ and $K=1$ and scan the $J$-axis. $\hat{L}_{ij}$ is the canonical angular momentum operator and it satisfies the commutation relation of $[\hat{L}_{ij},e^{\pm i\hat{\phi}_{ij}}]=\pm e^{\pm i\hat{\phi}_{ij}}$, and $J$ is the strength of the gauge field fluctuations. The overall phase diagram of Eq. (1) is obtained in the previous QMC work Xu et al. (2019) and is adapted here in Fig. 1 (a).

In the quantum Monte Carlo simulation in Ref.Xu et al. (2019), the quantum critical points can be extracted by means of correlation ratio, which is defined as $r=1-\frac{\chi(\mathbf{X}+\delta{\mathbf{q}})}{\chi(\mathbf{X})}$, where $\mathbf{X}$ is the order wavevector for AFM ($\mathbf{X}=(\pi,\pi)$) or VBS ($\mathbf{X}=(\pi,0)$) on the square lattice and $\delta(\mathbf{q})=(\frac{2\pi}{L},0)$ is the smallest momentum away from $\mathbf{X}$. The $\chi$ is the correlation function of the corresponding order that one probes, for example, the AFM order is determined by the $\chi$ of spin-spin correlation function $\chi_{S}(\mathbf{k})=\frac{1}{L^{4}}\sum_{ij}\sum_{\alpha\beta}\langle S_{\beta}^{\alpha}(i)S_{\alpha}^{\beta}(j)\rangle e^{i\mathbf{k}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}$ where the spin operator $S_{\beta}^{\alpha}(i)=c_{i\alpha}^{\dagger}c_{i\beta}-\frac{1}{N_{f}}\delta_{\alpha\beta}\sum_{\gamma}c_{i\gamma}^{\dagger}c_{i\gamma}$, and the VBS order is determined by the $\chi$ of dimer-dimer correlation function $\chi_{D}(\mathbf{k})=\frac{1}{L^{4}}\sum_{ij}\left(\langle D_{i}D_{j}\rangle-\langle D_{i}\rangle\langle D_{j}\rangle\right)e^{i\mathbf{k}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}$ with the dimer operator $D_{i}=\sum_{\alpha\beta}S_{\beta}^{\alpha}(i)S_{\alpha}^{\beta}(i+\hat{x})$ is defined as dimer along the nearest-neighbor bond in $\hat{x}$ direction.

By monitoring the corresponding correlation ratios, Ref. Xu et al. (2019) gives that for $N_{f}=2$, the transition of U1D-to-AFM is at $J=1.6(2)$; for $N_{f}=4$, the transition of U1D-to-VBS is at $J=1.2(3)$, and the transition of VBS-to-AFM is at $J=17(4)$; for $N_{f}=6$, the transition of U1D-to-VBS is at $J=1.9(3)$; for $N_{f}=8$, the transition of U1D-to-VBS is at $J=2.5(1)$. The illustration of the model in Eq. (1) is shown in Fig. 1 (b).

This compact U(1) lattice gauge theory coupled to fermionic matter at (2+1)D is now accessible with large-scale QMC simulations. The Hamiltonian in Eq. (1) can be formulated in a coherent-state path integral. To simulate the above model with determinantal QMC method, we start with the partition function as detailed in Ref. Xu et al. (2019, 2019),

[TABLE]

where the configuration space of $[\phi]$ is comprised of the (2+1)D gauge field. The bosonic part of the partition function is

[TABLE]

and the fermionic part of the partition function is

[TABLE]

Here all the flavors of fermion are subject to the same gauge field configuration, so for every fermion flavor, the $B^{\tau}$ matrix in the fermionic weight $\omega_{F}[\phi]$ is given by

[TABLE]

Since the gauge field $\phi_{i,j}(\tau)$ are continuous variables at the $(2+1)$D space time, and matrix elements in $B^{\tau}$ are complex numbers, it is very important to use an efficient strategy to update the gauge field $[\phi]$ Xu et al. (2019). We update the U(1) gauge field on $l$-th imaginary-time slice at $ij$-th lattice bond from $\phi_{ij}$ to $\phi^{\prime}_{ij}$. The ratio which determines whether we accept the updating can be expressed as

[TABLE]

For the boson part, the ratio of the weight is

[TABLE]

and for the fermionic part, the ratio of the weight is

[TABLE]

If the update is accepted, we also need update equal-time Green’s function as

[TABLE]

with the $2\times 2$ matrix of $\Delta$

[TABLE]

where

[TABLE]

The more detail of DQMC method used in this work can be found in Ref. Xu et al. (2019). It is with such QMC methodology that the ground state phase diagram as a function of fermion flavor $N_{f}$ and the strength of gauge fluctuations $J$ is mapped out, as shown in Fig. 1 (a).

In this paper, we focus on the large fermion flavor case ($N_{f}=8$) by means of stochastic analytic continuation (SAC) of imaginary-time correlation functions obtained from DQMC, where the deconfine-confine phase transition is investigated in detail Xu et al. (2019), the $1/L$ extrapolation of the correlation ratio crossings estimates U1D-to-VBS transition point at $J_{c}=2.5(1)$ for $N_{f}=8$.

The time displaced correlated function (defined as $G_{i,j}(\tau)=\langle\hat{O}_{i}(\tau)\hat{O}_{j}(0)\rangle$ and $G(\mathbf{q},\tau)=\frac{1}{N^{2}}\sum_{i,j}e^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}_{i,j}}G_{i,j}(\tau)$) of an operator $\hat{O}$ for a set of imaginary times $\tau_{i}\ (i=0,1,\cdots,N_{\tau})$ with statistical errors can be obtained from DQMC simulations. By SAC method Sandvik (2016); Shao et al. (2017); Huang et al. (2018); Qin et al. (2017b), the corresponding real-frequency spectral function $A(\omega)$ can be obtained from them according to the relationship of $G(\tau)=\int_{-\infty}^{\infty}d\omega A(\omega)K(\tau,\omega)$, where the kernel $K(\tau,\omega)$ depends on the type of the spectral function, i.e., fermionic or bosonic, finite or zero temperature. The spectra at positive and negative frequencies obey the relation of $A(-\omega)=e^{-\beta\omega}A(\omega)$ and we are restricted at the positive frequencies and the kernel can then be written as $K(\tau,\omega)=\frac{1}{\pi}(e^{-\tau\omega}+e^{-(\beta-\tau)\omega})$. In order to work with a spectral function that is itself normalized to unity on the positive frequency axis, we modify the kernel and the spectral function and arrive at the transformation between the imaginary time Green’s function $G(\mathbf{q},\tau)$ and real-frequency spectral function $B(\mathbf{q},\omega)$

[TABLE]

where $B(\omega)=A(\omega)(1+e^{-\beta\omega})$.

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

III Field Theory Analysis of U(1) Deconfined Phase

III.1 $\pi$-Flux State Mean-Field Theory

Before presenting our numerical result, we first provide a mean-field study of the spin and dimer excitation spectra in the U1D phase ignoring the $\mathrm{U}(1)$ gauge fluctuation (or considering the $J\to 0$ limit of Eq. (1) model). The mean-field treatment will be asymptotically exact in the large $N_{f}$ limit. At the mean-field saddle point, the fermions experiences a $\pi$-flux (per plaquette) background, described by the following Hamiltonian on the square lattice

[TABLE]

where $c_{i},c_{i}^{\dagger}$ are the creation and annihilation operators for the fermions on site $i$ with $N_{f}=8$ internal flavors. The $N_{f}$ lattice fermions will give rise to $2N_{f}$ Dirac fermions at low energy following the fermion doubling theorem. To see this, we transform the Hamiltonian to the momentum space,

[TABLE]

where we have chosen the four-site unit cell (sublattices are arranged surrounding a plaquette) such that $c_{\bm{k}}=(c_{{\bm{k}}A},c_{{\bm{k}}B},c_{{\bm{k}}C},c_{{\bm{k}}D})^{\intercal}$ and each component $c_{{\bm{k}}a}$ further contains $N_{f}=8$ flavors. Here $\sigma^{ij\cdots}\equiv\sigma^{i}\otimes\sigma^{j}\otimes\cdots$ denotes the tensor product of Pauli matrices and $\sigma^{0}$ stands for $2\times 2$ identity matrix. $\mathds{1}_{8\times 8}$ is identity matrix with dimension $8\times 8$, we will use $\mathds{1}$ for short. The fermion dispersion is given by $\epsilon_{\bm{k}}=\pm 2t(\sin^{2}k_{x}+\sin^{2}k_{y})^{1/2}$, which is gapless at the momentum ${\bm{k}}=(0,0)$. Expand around the Dirac point and rescale the theory to eliminate $t$, the low-energy continuum model can be written in terms of the Lagrangian density as,

[TABLE]

where $\gamma^{\mu}=(\mathrm{i}\sigma^{21}\mathds{1},-\sigma^{31}\mathds{1},\sigma^{10}\mathds{1})$ and $\bar{c}=c^{\dagger}\gamma^{0}$. Note that the gamma matrices are of dimensions $4N_{f}\times 4N_{f}$. Given that the minimal dimension of gamma matrices for (2+1)D Dirac fermion is $2\times 2$, the above flavor counting confirms that the $\pi$-flux model contains $2N_{f}$ Dirac fermions at low energy.

III.2 Spin and Dimer Excitation Spectrum

Given $N_{f}=8$ fermions on each site, our model has the $\mathrm{SU}(8)$ flavor symmetry on the lattice level. The $\mathrm{SU}(8)$ spin operators are defined as

[TABLE]

where $T^{a}$ (with $a=1,\cdots,63$) are the $N_{f}^{2}-1=63$ generators of $\mathrm{SU}(8)$. The dimmer operators along $x$ and $y$ directions are defined as $D_{i}^{x}=\sum_{a}S_{i}^{a}S_{i+\hat{x}}^{a}$ and $D_{i}^{y}=\sum_{a}S_{i}^{a}S_{i+\hat{y}}^{a}$ respectively. One may expect to expand them to four-fermion operators by inserting Eq. (16). However, the following fermion bilinear operators

[TABLE]

are gauge invariant (where we have replaced the dynamic gauge connection $\phi_{ij}$ by its mean field value specified in Eq. (13)) and symmetry-wise equivalent to the dimmer operators. In the large $N_{f}$ limit, the fermion bilinear operators are generally more relevant at low energy. So under the renormalization group flow, the dimmer operators should be represented by the fermion bilinear operator in Eq. (17).

The Fourier transform for generic operator is defined via $\mathcal{O}_{\bm{q}}=\sum_{i}\mathcal{O}_{i}e^{-\mathrm{i}{\bm{q}}\cdot{\bm{r}}_{i}}$. For fermion bilinear operators, they take the general form of

[TABLE]

where $v_{\bm{q}}$ is the vertex function (matrix) that depends on the momentum transfer ${\bm{q}}$. When ${\bm{q}}$ goes beyond the fermion Brillouin zone, we apply the following rules to map $c_{{\bm{k}}+{\bm{q}}}$ back: $c_{{\bm{k}}+(\pi,0)}=\sigma^{30}\mathds{1}c_{\bm{k}}$, $c_{{\bm{k}}+(0,\pi)}=\sigma^{03}\mathds{1}c_{\bm{k}}$, $c_{{\bm{k}}+(\pi,\pi)}=\sigma^{33}\mathds{1}c_{\bm{k}}$. Applying to the spin and dimer operators, we explicitly have

[TABLE]

With these, we can evaluate the correlation function for spin or dimer operators

[TABLE]

from which the spectral function $A_{\mathcal{O}}({\bm{q}},\omega)=-2\operatorname{Im}\chi_{\mathcal{O}}({\bm{q}},\omega+\mathrm{i}0_{+})$ can be obtained. In the zero temperature limit, the spectral function is given by

[TABLE]

where $\lvert n,{\bm{k}}\rangle$ is the $n$-th eigenstate of the single-particle Hamiltonian Eq. (14) with momentum ${\bm{k}}$ and $\xi_{n,{\bm{k}}}$ is the corresponding eigen energy, and $\Theta$ denotes the step function. Given the $v_{\bm{q}}$ in Eq. (19), the above calculation will provide us the understanding of the overall shape of the spectral function for both spin and dimer correlations in the $J\to 0$ limit of the lattice QED model in Eq. (1), where $\mathrm{U}(1)$ gauge fluctuation is supressed. We will demonstrate the spectra in Fig. 3 and Fig. 5 in Sec.IV, and compare them with our QMC result involving gauge fluctuations. We find that the low-energy spectral features match nicely on the qualitative level between the free fermion and the QMC results (although the scaling dimensions will be altered by gauge fluctuations).

III.3 Emergent Symmetry and Conserved Currents

The mean-field understanding of the excitation spectrum helps us to identify signatures of emergent symmetry in the U1D phase. Let us restore the $\mathrm{U}(1)$ gauge fluctuation in the following discussion, and consider the compact QED theory with $16$ Dirac fermions in $(2+1)$D spacetime. The Lagrangian in Eq. (15) becomes

[TABLE]

where $\mathcal{M}$ and $\mathcal{M}^{\dagger}$ are the annihilation and creation operators of the $2\pi$ flux of the $\mathrm{U}(1)$ gauge field, also known as the monopole operator (event) in the spacetime.Song et al. (2018a, b) Such monople terms are generally allowed if not forbidden by the physical symmetry. Here the physical symmetry includes the spin $\mathrm{SU}(8)$ symmetry and the four-fold rotation symmetry $\mathbb{Z}_{4}$ of the square lattice, which act on the fermion field $c$ as

[TABLE]

where $\tau^{3}=\sigma^{12}\mathds{1}$ generates the rotation between $x$-VBS and $y$-VBS. Based on operator-state correspondence, the $\mathrm{U}(1)$ gauge monopole operator with charge 1 can be effectively mapped to the state on $S^{2}\times\mathbb{R}$ with 1 unit of background magnetic flux through $S^{2}$, and the states contain fermion zero mode guaranteed by Atiyah-Singer index theorem Borokhov et al. (2002). Therefore, when a $\mathrm{U}(1)$ monopole operator is inserted, each Dirac cone will contribute a zero mode, so there are totally 16 zero modes. Different ways of filling these zero modes leads to different monopole states that are degenerated in energy. A gauge neural monopole must have these fermion zero modes half-filled (i.e. filling 8 fermions on 16 zero modes), which results in $C_{16}^{8}=12870$ different monopole states. Among them, only $9$ states preserves the spin $\mathrm{SU}(8)$ symmetry. They can be labeled by the following quantum number

[TABLE]

which corresponds to the monopole angular momentum because $\tau^{3}/2$ is the generator of the $\mathbb{Z}_{4}$ lattice rotation symmetry. For example, the $m=4$ state is created by $\mathcal{M}^{\dagger}\sim\prod_{\alpha=1}^{8}\frac{1+\tau^{3}}{2}c_{\alpha}^{\dagger}$, which fills all the 8 fermions on the monopole modes of the same valley (of $\tau^{3}=+1$). Each fermion occupies a distinct $\mathrm{SU}(8)$ spin flavor, such that the monopole state is $\mathrm{SU}(8)$ symmetric. Further imposing the $\mathbb{Z}_{4}$ rotation symmetry to the monopole, its angular momentum must satisfy $m=0\mod 4$, which further singles out 3 monopole states labeled by $m=0,\pm 4$. These states span the Hilbert space of a single monopole that preserve all the $\mathrm{SU}(8)\times\mathbb{Z}_{4}$ physical symmetry, so their corresponding monopole operators are generally allowed to appear in the Lagrangian of the QED theory Eq. (22).

Depending on the relevance of the single monopole term $g$ at the large-$N_{f}$ fixed point, the lattice QED model in Eq. (1) can have different emergent symmetry in the U1D phase. The scaling dimension $\Delta_{\mathcal{M}}$ of the single monopole operator has been calculated in the large-$N_{f}$ limit in Ref. Borokhov et al., 2002; Dyer et al., 2013; Pufu and Sachdev, 2013, which reads $\Delta_{\mathcal{M}}=0.265\times(2N_{f})-0.0383+\mathcal{O}(1/N_{f})$ 111Note that the notion of $N_{f}$ in the lattice model is differed from that in the field theory by a factor of 2. With $N_{f}=8$, $\Delta_{\mathcal{M}}=4.2>3$, so the single monopole operator is irrelevant to the leading orders in $1/N_{f}$, nevertheless the conclusion may still change at higher orders of $1/N_{f}$. But if we accept that the single monopole term is irrelevant, the theory will flow to the QED3 fixed point, where the emergent symmetry is the full $\mathrm{SU}(16)/\mathbb{Z}_{16}=\mathrm{PSU}(16)$ flavor symmetry of the $16$ Dirac fermions, where the $\mathbb{Z}_{16}$ center of the $\mathrm{SU}(16)$ group should be quotient out because this subgroup is shared with the $\mathrm{U}(1)$ gauge group. The $\mathrm{SU}(16)$ generators can be enumerated as $\{{\bm{\tau}},T^{a},{\bm{\tau}}T^{a}\}$. Here ${\bm{\tau}}=(\sigma^{01}\mathds{1},\sigma^{13}\mathds{1},\sigma^{12}\mathds{1})$ are the generators of valley rotations. These generators are found by requiring them to commute with $\gamma^{\mu}$, such that the Lagrangian in Eq. (22) remains invariant under the fermion flavor rotation. Using the $\mathrm{SU}(16)$ generators, one can define the $\mathrm{SU}(16)$ currents (labeled by $i=1,2,3$ and $a=1,\cdots,63$)

[TABLE]

There are $255$ current operators in total (each current further contains 3 spacetime components labeled by $\mu=0,1,2$). All these currents are emergent conserved currents at low energy.

However, although unlikely, if the single monopole operator turns out to be relevant and if we assume the theory flows to another non-trivial conformal fixed point (when $N_{f}$ is within the conformal window), the emergent symmetry can be lowered by the non-vanishing monopole term. The single monopole term $g$ will break the emergent symmetry from $\mathrm{PSU}(16)$ to

[TABLE]

The above symmetry group is most easily seen for the $m=4$ monopole: the two $\mathrm{SU}(8)$ acts on the spin flavors in the two valleys ($\tau^{3}=\pm 1$) respectively, $\mathbb{Z}_{8}$ is the opposite 8-fold rotation of fermion phases in opposite valleys, and all the center subgroups must be quotient out as they are shared between the $\mathrm{U}(1)$ gauge group. More careful symmetry analysis for the other monopoles of different $m$ shows that the above symmetry is indeed the largest possible residual symmetry of a single monopole operator. In this case, the emergent conserved currents are reduced to

[TABLE]

In this case, there are $126$ emergent conserved currents in total. We summarize the above analysis in Tab. 1.

Our analysis shows that the relevance of the $\mathrm{U}(1)$ gauge monopole crucially affects the emergent symmetry and the emergent conserved currents that can be probed at low-energy. Identifying these current fluctuations in the spin and dimer excitation spectra will be the first step towards pinning down the emergent symmetry and studying the monopole effects in the lattice QED model. The analysis can be carried out on the mean-field level. According to Eq. (19) and $\bar{c}=c^{\dagger}\gamma^{0}$, $\gamma^{\mu}=(\mathrm{i}\sigma^{21}\mathds{1},-\sigma^{31}\mathds{1},\sigma^{10}\mathds{1})$, $\tau^{i}=(\sigma^{01}\mathds{1},\sigma^{13}\mathds{1},\sigma^{12}\mathds{1})$, we can identify the following spin and dimer operators to current operators

[TABLE]

They are summarized in Tab. 2. Among them, $j_{0a}^{0}$ is the conserved current of the physical spin $\mathrm{SU}(8)$ symmetry, and the remaining currents $j_{1a}^{2},j_{2a}^{1},j_{30}^{2},j_{30}^{1}$ are all emergent conserved current of $\mathrm{PSU}(16)$ but not of $\mathrm{PSU}(8)\times\mathrm{PSU}(8)$, see Tab. 1. By measuring the scaling dimension of these currents from the spin and dimer correlation functions, one can decide if they are conserved or not to further determine the emergent symmetry and the relevance of the single monopole operator Ma et al. (2019). At current stage, our numerical resolution is not sufficient to fully resolve the scaling dimension of these current fluctuations, nevertheless we will first map out the overall shape of the excitation spectra and identify these low-energy spectral features in this work and provide a road map for future study of the emergent conserved currents.

It is worth mentioning that the dimer fluctuation at $(0,0)$ momentum is gapless, but its spectral weight fades away much faster towards low energy as shown in Fig. 5 (a) and (b). This continuum corresponds to the energy-momentum tensor $T^{\mu\nu}$ which is the conserved current associated with the translation symmetry,

[TABLE]

Based on this definition, we can identify that

[TABLE]

The scaling dimensions of $\mathcal{L},\mathcal{T}^{\mu\nu}\sim\delta\mathcal{L}/\delta g_{\mu\nu}$ are 3, which follows from the fact that the action $S=\int\mathrm{d}^{3}x\mathcal{L}$ and the metric $g_{\mu\nu}$ must be dimensionless. Because of the relatively high scaling dimension of the $D_{(0,0)}^{x}$ fluctuation, it is much more irrelevant under renormalization compared to the current and order parameter fluctuations. Therefore the low-energy spectral weight of $D_{(0,0)}^{x}$ is expected to be much weaker compare to other continua (e.g. $D_{(\pi,0)}^{x}$, $D_{(\pi,\pi)}^{x}$) in the dimer excitation spectra in Fig. 5 in Sec. IV.

IV Quantum Monte Carlo Results

Here we present the QMC results, first begin with the definition of the physical observables.

IV.1 Physical observables

To understand the deconfine-confine phase transition, we focus on gauge-invariant dynamical structure factors obtained in QMC simulations, including the spin and dimer dynamical structure factor. They can be defined as the following formsXu et al. (2019); Zhou et al. (2018)

[TABLE]

[TABLE]

where the spin operator is $S_{\beta}^{\alpha}(i)=c_{i\alpha}^{\dagger}c_{i\beta}-\frac{1}{N_{f}}\delta_{\alpha\beta}\sum_{\gamma}c_{i\gamma}^{\dagger}c_{i\gamma}$ and the dimer operator is $D_{i}=\sum_{\alpha\beta}S_{\beta}^{\alpha}(i)S_{\alpha}^{\beta}(i+\hat{x})$. The dimer operator is defined along the nearest-neighbor bond in $\hat{x}$ direction.

As mentioned in Sec. II, the time displaced correlated functions $S({\bm{q}},\tau)$ and $D({\bm{q}},\tau)$ can be obtained in QMC for a set of imaginary times $\tau_{i}\ (i=0,1,\cdots,N_{\tau})$ with statistical errors. From which, the SAC will be further applied to extract the real-frequency spectral functions $S({\bm{q}},\omega)$ and $D({\bm{q}},\omega)$.

Another quantity, that has distinctively different behaviors in U1D and VBS confined phases, is the net fluctuation of flux in each time slices with Monte Carlo steps Xu et al. (2019); Kleinert et al. (2002). Flux in each plaquette can be written as $\sum_{b\in\Box}\phi_{b}=\Phi_{\Box}+2\pi m_{\Box}$ with $\Phi\in[0,2\pi)$ and $m_{\Box}$ an integer. The fluctuation of net flux in one time slice $M(\tau)$ is defined as a sum of $m_{\Box}$ of each plaquette at time slice $\tau$, $M(\tau)=\sum_{\Box}m_{\Box}(\tau)$. The evolution of $M(\tau)$ with Monte Carlo time series, both inside U1D and VBS phases at time slices $\tau_{1}$ and $\tau_{2}$, are shown in Fig. 2 (a) and (b), respectively. The parameters of calculation were given as $L=12$ and $\beta=2L$. Inside the U1D phase ($J=1.80<J_{c}$), as shown in Fig. 2 (a) the net flux favors $\pi$ flux in each plaquette, and the net fluctuation $M(\tau)$ at each time slice seldom changes, $M(\tau_{1})$ and $M(\tau_{2})$ follow closely to each other and their value only take the integers [math], $\pm 1$ and $\pm 2$; while in the VBS phase ($J=3.50>J_{c}$), as shown in Fig. 2 (b) the net fluctuations change almost randomly with more extended values, [math], $\pm 1,\pm 2,\pm 3,\pm 4,\pm 5,\pm 6$ (in unit of $2\pi$), and large deviation between different time slices $\tau_{1}$ and $\tau_{2}$ can all be seen. These large fluctuations in the net flux indicate the proliferate of monopoles in the confined VBS phase.

IV.2 Spectra and excitation gaps

In this part we present $S({\bm{q}},\omega)$ and $D({\bm{q}},\omega)$ inside both the U1D and VBS phases, these results are obtained from QMC-SAC simulations. We also show the corresponding spectra from the non-interacting $\pi$-flux model of Dirac fermions without gauge fluctuations, which are the spectra at the limit of $J=0$ derived in Sec. III.2.

IV.2.1 Spin Spectra in UID and VBS phases

Fig. 3 shows the features of the spin spectra through the U1D-VBS transition, the results are shown along the high-symmetry-path of $(0,0)-(\pi,0)-(\pi,\pi)-(0,0)$ . We present results for the non-interacting Dirac fermions corresponding to $J=0$ (Fig. 3 (a)), inside the U1D phase with $J=2.00<J_{c}$ (Fig. 3 (b)), close to the QED3-GN critical point at $J=3.10$ (Fig. 3 (c)) and inside the VBS confined phase at $J=4.50>J_{c}$ (Fig. 3 (d)).

The $\pi$-flux spinons, as discussed in Sec. III.2 with the U(1) gauge fluctuations suppressed, give rise to gapless spin spectra at momenta $(0,0)$, $(\pi,0)$ and $(\pi,\pi)$ in Fig. 3 (a). The situation persists as one introduces the U(1) gauge fluctuations, for example at $J=2$ in Fig. 3 (b). Of course on a finite size system $L=14$ for Fig. 3 (b), the spectra look gapped due to finite size effect, we have performed the $1/L$ extrapolation of the spin gaps at $(\pi,0)$ and $(\pi,\pi)$ with the gaps directly obtained from fitting the imaginary time decay of $S({\bm{q}},\tau)$ without SAC, the results are shown in Fig. 4 (a) and (d), and it is clear that in the U1D phase, the spin excitations at $(\pi,0)$ and $(\pi,\pi)$ are gapless in the thermodynamic limit. As discussed in Sec. III.3, the $(\pi,\pi)$ excitation corresponds to the spin $\mathrm{SU}(8)$ order parameter fluctuation, and the $(\pi,0)$ excitation corresponds to the current fluctuation whose charge operator generates the AFM-VBS rotation. If the emergent symmetry is $\mathrm{PSU}(16)$, the scaling dimension of the spin excitation at $(\pi,0)$ will be pinned at 2. However if the emergent symmetry is $\mathrm{PSU}(8)\times\mathrm{PSU}(8)\times\mathbb{Z}_{4}$, the scaling dimension will deviate from integer. More importantly, we also observed broad and prominent continuous spectra in Fig. 3 (b), which reflects the expected deconfinement and fractionalization of spinons and their interactions mediated by the fluctuating U(1) gauge field. Similar $S({\bm{q}},\omega)$, with gapless excitations at $(0,0)$, $(\pi,0)$ and $(\pi,\pi)$ and pronounced continua upto high energy, have also been seen at the deconfined quantum critical point with emergent O(4) symmetry Ma et al. (2018, 2019).

As one moves towards the critical point (actually slightly above it in Fig. 3 (c) with $J=3.1>J_{c}=2.5$), broad and prominent continuous spectra can still be observed, signifying the effects of gauge fluctuations. And the $1/L$ extrapolation of the spin gap are shown in Fig. 4 (b) and (e). It is clear now that once inside the VBS confined phase, the spin spectra are gapped due to the translational symmetry breaking of the VBS phase.

Deep inside the VBS phase, as shown in Fig. 3 (d) with $J=4.5$, the spin spectra are fully gapped and the continua above it also become less extended in frequency domain. This is expected as well since here both the gauge fields and the fermions are interacting at the length scale shorter than that associated with the excitation gaps. Below the gap, the system is an insulator with fermions forming singlets along the $(\pi,0)$ or $(0,\pi)$ directions, i.e., translational symmetry breaking. The corresponding $1/L$ extrapolation of the gaps at $(\pi,0)$ and $(\pi,\pi)$ are shown in Fig. 4 (c) and (f), respectively.

IV.2.2 Dimer Spectra in UID and VBS phases

Fig. 5 shows the dynamic dimer spectra through the U1D-to-VBS transition.

The dimer spectra of the non-interacting $\pi$-flux Dirac fermions are given in Fig. 5 (a). The spectra are gapless at momenta $(0,0)$, $(\pi,0)$ and $(\pi,\pi)$, similar to that of the spin in Fig. 3 (a), but the spectral weights have different distribution. The calculation details of the non-interacting spin and dimer spectra are given in Sec. III.2.

As one moves into the U1D phase at $J=2.0<J_{c}$ (Fig. 5 (b)), the gapless dimer spectra persist. Again the spectral gaps at the high symmetry points are due to the finite size effect. The interesting observation here is that the continua are very broad, extending all the way from $\omega=0$ to $\omega=8J$, i.e., beyond the upper boundary of the non-interacting spectra. This points to the importance of higher order continua of multi-spinon excitations due to the strong interaction effect mediated by the U(1) gauge field fluctuations. It is also interesting to notice that at low energy, $\omega/J\leq 2$, the spectral weight bear similar distribution with that of the non-interacting one, in particular, the weight is greatly reducing as one approaches momentum $(0,0)$, this is related with the fact that $D^{x}_{(0,0)}$ is the energy-momentum tensor with larger scaling dimensions, as pointed out in Sec. III.3. Also, the scaling dimension of $D^{x}_{(\pi,\pi)}$ could help with distinguishing the emergent symmetry of $\mathrm{PSU}(16)$ or $\mathrm{PSU}(8)\times\mathrm{PSU}(8)\times\mathbb{Z}_{4}$ in the U1D phase, given larger system sizes could be simulated.

Near and slightly above the critical point at $J=3.1$, as shown in Fig. 5 (c), broad and prominent continuous spectra can still be observed and there are gapless spectra at $(\pi,0)$. The gapless excitation close to $(\pi,0)$ are the critical fluctuations associated with the QED3-GN transition. With larger system sizes and lower temperature in the future QMC studies, one will be able to measure the anomalous dimension exponent $\eta$ from the momentum and frequency dependence of such critical fluctuation, and could compare with the predictions of QED3-GN transitions from the recent perturbative RG calculations Gracey (2018); Ihrig et al. (2018); Zerf et al. (2019); Dupuis et al. (2019).

Inside the VBS phase, the dimer spectra are gapped due to the $(\pi,0)$ or $(0,\pi)$ translational symmetry breaking. However, since the dimer order parameter will contribute a Bragg peak at $(\pi,0)$ and $\omega=0$, the analytic continuation is notoriously difficult for finding such not-smoothed spectra, i.e., one delta peak at $\omega=0$ followed with a gap and then continua above it. To solve this problem, we add a small pinning field to strengthen the VBS order Assaad and Herbut (2013) in the simulation of Fig. 5 (d). The pinning fields are added according to the pattern of Fig. 6 (a), with $t^{\prime}=1.05t$, and the simulation results are consistent with the expectation, in that, in Fig. 5 (d), the spectra looked gapped at low energy, however, a fixed momentum cut at $(\pi,0)$, as shown in Fig. 6 (b), indeed reveals that there is a Bragg peak at $\omega=0$ and a continuous spectra beyond a gap due to the break of discrete symmetry in VBS phase.

V Conclusions and discussion

In this work, we have performed both numerical and analytical analyses of the dynamics of a model realizing the compact QED3 at large fermion flavor ($N_{f}=8$). As mentioned in the introduction, the question of U(1) gauge field coupled to fermionic matter field at (2+1)D is of high interests to both condensed matter and high-energy physics communities. The U1D phase is a realization of the algebraic quantum spin liquid in which Dirac spinons dynamically coupled to the emergent U(1) gauge field. The transition of U1D-VBS is the deconfinement-to-confinement transition in the QED3 setting and is also the transition from symmetric quantum spin liquid to the symmetry-breaking phase that several potential quantum spin liquid compounds could have already realized by tuning the doping concentration, pressure and magnetic field Zheng et al. (2017); Kitagawa et al. (2018). The dynamical information of the U1D phase and U1D-VBS transition revealed here – the continua in spin and dimer spectra and their field theoretical interpretations in the emergent symmetries and conserved currents – provide the first piece of concrete evidence of the aforementioned exotic physical phenomena.

Looking forward, better algorithm in QMC simulations would certainly be desirable to access larger system sizes and lower temperatures. In particular, the critical properties of the U1D-VBS transition, that of the QED3-GN types, have already been discussed in the high-order perturbative RG calculations Gracey (2018); Ihrig et al. (2018); Zerf et al. (2019); Dupuis et al. (2019), but the system sizes in this work is too small to extract accurate values of the critical exponents. Further developments, in terms of algorithm improvement and more focus close to the QED3-GN critical points, are on-going.

From analytical perspective, calculation of the spectra with the fluctuations of the U(1) gauge fields included would be very useful, similar as the analysis in the Ref. Ma et al. (2018), the dynamical signature of the strongly correlated systems of fractionalized spinons and their coupling effects with the emergent gauge field could be revealed and provide clearer guidance for future numerical simulations and eventually to experiments.

ACKNOWLEDGMENTS

We thank Yang Qi, Lukas Janssen, Michael Scherer, Joseph Maciejko, John Gracey, Cenke Xu, Chong Wang, Yin-Chen He and Liujun Zou for helpful discussions. WW and ZYM acknowledge the supports from the Ministry of Science and Technology of China through the National Key Research and Development Program (Grant No. 2016YFA0300502), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB28000000), and the National Science Foundation of China (Grant No. 11421092, 11574359 and 11674370). X. Y. X. is thankful for the support of Hong Kong Research Grants Council (HKRGC) through C6026-16W, 16324216 and 16307117. We thank the Center for Quantum Simulation Sciences at Institute of Physics, Chinese Academy of Sciences, and the Tianhe-1A platform at the National Supercomputer Center in Tianjin for their technical support and generous allocation of CPU time.