Quantum Electrodynamics in 2+1 Dimensions as the Organizing Principle of a Triangular Lattice Antiferromagnet

Quantum electrodynamics in 2+1 dimensions (QED 3 ) has been proposed as a critical field theory describing the low-energy effective theory of a putative algebraic Dirac spin liquid or of quantum phase transitions in two-dimensional frustrated magnets. We provide compelling evidence that the intricate spectrum of excitations of the elementary but strongly frustrated J 1 - J 2 Heisenberg model on the triangular lattice is in one-to-one correspondence to a zoo of excitations from QED 3 , in the quantum spin liquid regime. This includes a large manifold of explicitly constructed monopole and bilinear excitations of QED 3 , which is thus shown to serve as an organizing principle of phases of matter in triangular lattice antiferromagnets and their low-lying excitations. Moreover, we observe signatures of emergent valence bond solid (VBS) correlations. This can be interpreted either as evidence of critical VBS fluctuations of an emergent Dirac spin liquid or as a transition from the 120 ◦ N´eel order to a VBS whose quantum critical point is described by QED 3 . Our results are obtained by comparing ansatz wave functions from a parton construction to exact eigenstates obtained using large-scale exact diagonalization up to N = 48 sites.


I. INTRODUCTION
The emergence of collective equations of motion from seemingly unrelated microscopic interactions is one of the most fascinating aspects of many-body physics.Strongly correlated electrons can realize intriguing quantum field theories (QFT) as their low-energy effective description which otherwise are used to describe the fundamental laws of elementary particles.Quantum spin liquids (QSL) in frustrated magnetism are a particularly exciting instance of emergent QFTs [1].Topological QFTs describe certain gapped QSLs which have been shown both analytically and numerically to emerge in local spin models.This includes the emergent Z 2 lattice gauge theory in the toric code or the Kitaev's honeycomb model [2] as well as Chern-Simons theories realized in chiral spin liquids [3,4], which have been discovered in simple Heisenberg-like Hamiltonians on the triangular and kagome lattice [5][6][7][8][9][10][11].
The arguably most widely known QFT is quantum electrodynamics (QED).While on one hand, it is the fundamental theory of fermions coupled to a U(1) gauge field describing the physics of elementary electrons and photons, it has also been discussed as an emergent field theory in frustrated magnets.Remarkably, QED in three spatial dimensions can be realized in pyrochlore spin ice compounds [12][13][14].Condensed matter systems also allow for the realization of QED in less than 3 spatial dimensions.The physics of QED in 2 + 1 dimensions (QED 3 ) is considered to be more strongly coupled than its 3+1 dimensional counterpart while exhibiting a richer phenomenology than the confining 1 + 1 dimensional QED, also referred to as the Schwinger model.However, the physics of QED 3 is until today still a subject of intense research.
Early on, QED 3 has been suggested as an effective field theory for so-called algebraic or Dirac spin liquids (DSL) in quantum magnets [15][16][17].While model wave functions of these states have been studied [18][19][20] and signatures of gapless spin liquids have been detected in certain spin models [21,22], no realization of Dirac spin liquid has until today been unambiguously confirmed.However, several numerical studies have recently highlighted the relevance of QED 3 for elementary frustrated quantum magnets, such as the kagome Heisenberg antiferromagnet or the J 1 -J 2 Heisenberg model on the triangular lattice.
In this work, we demonstrate that QED 3 serves as on organizing principle of the physics of triangular lattice Heisenberg antiferromagnets by exact numerical calculations.Specifically, we consider the spin-1/2 Heisenberg model on the triangular lattice, with nearest-neighbor J 1 and next-nearestneighbor J 2 antiferromagnetic (AF) couplings, We construct a "vacuum" of QED 3 and a plethora of low-lying excitations directly on a lattice by means of Gutzwiller projection of a DSL parton ansatz.These model states of both the vacuum and the low-lying excitations are then directly compared to numerically exact ground state and low-energy eigenstates of an N = 36 sites simulation cluster.Moreover, we present the complete low-energy spectrum with space group and spin-parity resolution for the N = 48 sites cluster, obtained using large-scale Exact Diagonalization (ED).The N = 48 ED spectrum used approximately 25 million CPU hours which was possible thanks to an overall 40 million CPU hours grant from PRACE.The ground state phase diagram of Eq. ( 1) features a 120 • Néel ordered state for J 2 /J 1 ≲ 0.09 [23][24][25][26] and a stripy antiferromagnetic state for J 2 /J 1 ≳ 0.14.In between, a paramag- netic regime is stabilized whose nature has been subject of intense debate [20,[27][28][29][30][31].In particular, several references have pointed out the possibility of a Dirac spin liquid [20,22,27] or a gapped Z 2 spin liquid [28][29][30]32] being stabilized in this regime.Quite interestingly, the paramagnetic regime is highly sensitive to further perturbations: for instance, when adding a further neighbor AF coupling J 3 , a gapless chiral spin liquid (CSL), which spontaneously breaks time and lattice reflection symmetry, was reported using DMRG [33]; similarly, adding an explicit tiny chiral term in the Hamiltonian (J χ S i • (S j × S k )) leads to a gapped CSL [9,34,35].

II. QED3 AND ITS TORUS SPECTRUM
Based on a parton approach, it has been advocated that QED 3 , i.e. the quantum field theory of N f flavors of massless Dirac fermions coupled to a compact U (1) gauge field in 2+1D, could be an appropriate description of the low-energy physics of certain spin liquid phases, in particular for the kagome and the triangular lattice Heisenberg models.In these examples, the number of fermion flavors is N f = 4.
In our work, we want to study the low-energy spectrum of finite-size spin systems and investigate whether the spectrum is in agreement with the low-energy spectrum of QED 3 in finite volume.In order to execute this program we need to review what is known about QED 3 and how the finite size low-energy spectrum of a quantum field theory on a torus is structured.
It is useful to start understanding QED 3 in the limit of N f → ∞.In that limit, the gauge fluctuations are suppressed and the theory has a simpler structure.The theory is a conformal field theory (CFT) with a variety of gapless modes.On the one hand, the bare massless fermions remain gapless, as well as the gauge field in the form of photons.Another set of excitations, the monopoles -which will play an important role later -, are found at a large energy scale which is proportional to N f in the limit N f → ∞.In the correspondence to the spin model, only the gauge invariant and charge neutral excitations of QED 3 are allowed to appear, therefore only neutral fermion excitations are visible (e.g.vacuum, neutral bilinears c.f. Fig. 1(b)) in addition to the gauge field excitations and the monopoles, see Fig. 1(c).As N f decreases towards the value of interest N f = 4, the monopoles also become gapless lowenergy excitations.The role of the monopoles is also crucial for the stability of the CFT window of QED 3 itself.In the extreme limit N f = 0 it is known since Polyakov [36] that a pure compact U (1) gauge theory in 2+1D is confining in the presence of a UV cutoff.Analytical and numerical results on the N f extent of the conformal window of QED 3 did not yet reach a consensus, but it is conceivable that N f = 4 still belongs to the conformal window [37][38][39][40].
As QED 3 is a CFT for large enough values of N f , one can characterize the low energy excitations based on their scaling dimensions.The current best estimates for the scaling dimensions at N f = 4 are summarized in Ref. [40].For 1+1D systems described by a CFT, the energy spectrum on a circle is equivalent to the operator content on the 2D (space-time) torus and harbors thus the spectrum of scaling operators and their descendants, arranged into Verma modules.In 2+1D CFTs an analogous correspondence only holds for Hamiltonians quantized on a spatial sphere, while the lattice models studied in condensed matter physics more naturally live on a spatial torus.Recently the torus energy spectrum for a series CFTs, such as the Wilson-Fisher and Gross-Neveu-Yukawa theories have been studied using a combination of numerical and analytical results [41][42][43][44][45].The spectrum is understood to collapse as 1/L (L being the linear extent of the system) as expected for a relativistic theory.The low-energy spectrum ∆ i = E i − E GS multiplied by L then forms a fingerprint ξ i ≡ ∆ i × L of the conformal field theory governing the low- energy spectrum [46].While the torus spectrum of these theories is known to be different from the sphere spectrum [47], in the above works a phenomenological reminiscence among some of the low-lying states on the sphere and on the torus has been observed.So we expect the low-energy part of the torus spectrum of QED 3 to be formed by the torus doppelgänger of the vacuum, the neutral bilinears, the monopoles, and gaugefield excitations.
The torus spectrum of QED 3 on a (square) torus has been studied in the limit N f → ∞ in Ref. [43].This includes the bilinears and the photon excitations.We discuss the basic structure of the neutral fermionic sector adapted to the hexagonal Brillouin zone torus in App.F. A hallmark of the fermionic bilinear spectrum is the massive degeneracy of levels already at the first fermionic bilinear excitation above the ground state (larger than their number N 2 f on the sphere), combined with a rather soft photon excitation which additionally boosts the number of low-lying levels.So even in the absence of monopoles, the low-energy torus spectrum of QED 3 is much denser than for example the torus spectrum of a 2+1D Ising CFT [41].
The monopole excitations in lattice systems carry quantum numbers of the space and the spin symmetry group.In Refs.[48,49] the quantum numbers of the first q = ±1 monopoles for N f = 4 have been studied for various lattices.For the triangular lattice of interest here, spin-singlet monopoles at each of the six X points in the Brillouin zone have been identified together with spin-triplet monopoles at each of the two K points.Altogether these form twelve monopole states in total.We expect precisely these states to be present at low energy in the torus spectrum.Note that the SU (4) symmetry of QED 3 for N f = 4 predicts the energies of the 12 monopole states to be degenerate, also on the torus.Since this symmetry is not exact on the lattice but only emerges in the IR, one would expect some finite-size splitting to lift the degeneracy.
In the absence of an analytical result for the torus spectrum for N f = 4 we proceed now by a Gutzwiller projection procedure to generate model state wave functions for several classes of torus excited states of QED 3 on the triangular lattice.This way we will be able to describe the vacuum, the bilinears, and the monopoles, see Sec.IV.While we cannot predict the ξ i values of the torus spectrum this way, we can still check whether the significant overlaps of the Gutzwiller projected states cover relevant low-energy states of the ED spectrum on a given system size.In the next step, we can then interpret the relative excitation energies of the different classes of excitations, and possibly gain some insights about the relevance of QED 3 for the low energy physics of the triangular lattice J 1 − J 2 Heisenberg model, and perhaps even about the stability of QED 3 itself.
Note that we are not able to generate excitations of the gauge field sector in this way.In order to shed some light on this sector, we take a different approach and consider the spectrum of the Rokhsar-Kivelson quantum dimer model [50] on the same finite-size clusters, see Sec.IV E, for comparison.

III. SPECTRUM OF THE J1-J2 TRIANGULAR ANTIFERROMAGNET
We now turn to the main problem of interest, to understand the low-energy spectrum of the J 1 − J 2 Heisenberg model on the triangular lattice and the nature of the corresponding phases.The low-energy spectrum in many phases of matter is well understood, for example in magnetically ordered phases [51], or in topologically ordered phases with their torus ground state degeneracy.As discussed in the previous section, torus spectra for conformal field theories in 2+1D are however a topic of ongoing research.
We study the energy spectra of Eq. ( 1) obtained from largescale exact diagonalization [52] where we resolve the spec-trum w.r.t. the irreducible representations (irreps) of translational, point group (PG), and spin-flip symmetries in Fig. 2. States in the even (odd) spin-flip representation have even (odd) total spin S. We denote the irrep of these states by k.ρ, where k labels the momentum and ρ labels the point group irrep.For N = 36 Fig. 2(a) we cover the range of The ground state for J 2 /J 1 ≲ 0.09 is magnetically ordered with a 120 • Néel pattern and ordered with a stripy pattern for J 2 /J 1 ≳ 0.14 [53].In between, a quantum paramagnetic regime is stabilized, whose nature is the main subject of our discussion.Far outside the paramagnetic regime, e.g. for J 2 /J 1 = −0.1 in the 120 • Néel phase or at J 2 /J 1 = 0.2 in the collinear stripy phase, the main structure of the low-lying energy spectrum can be well understood by tower-of-states (TOS) analysis [51,53].In the 120 • Néel phase, TOS analysis predicts an S = 0 ground state in the Γ.A1 sector whereas the first two quasi-degenerate excited states are in the K.A1 and Γ.B1 sector with S = 1, precisely what is observed in the J 1 -J 2 model in Fig. 2(a).Similarly, TOS analysis in the stripy phase predicts three quasidegenerate S = 0 states, one at Γ.A1 and two forming a Γ.E2 irrep.For details on the TOS analysis and further predictions see App.D and Ref. [51].In the intermediate regime, however, the spectrum is rather dense and low-lying excitations belong to various irreps of the symmetry group.As shown in the subsequent sections, many of these levels can be identified with non-trivial excitations on top of the QED 3 vacuum.In particular, the prominent singlet levels with momentum k = X and k = M , which are neither part of the TOS for the 120 • Néel nor stripy order, will be related to monopole and bilinear excitations of QED 3 .Fig. 2(b) shows the low-energy spectrum of the N = 48 cluster at J 2 /J 1 = 0.125 in the paramagnetic regime.Again, we observe several low-lying singlet excitations (filled symbols) at non-trivial momenta and PG irreps.In particular, two singlet Γ.E2, one Γ.A2, and two levels at X 1 .A and X 2 .A which are sixfold degenerate can be found below the spin gap, and lower than their N = 36 occurrence.The lowest lying triplets can be found in the K.A1, Z.A and M .A1 sectors, i.e. along the Brillouin zone boundary at roughly the same energy.We refer to App.B for a detailed comparison of the two system sizes studied here.
We will provide an attempt at understanding the presence of these low-energy states from the QED 3 perspective in the following.

IV. FROM QED3 TO THE SPECTRUM OF THE J1-J2 MODEL
To elucidate the structure of the energy spectrum in the intermediate paramagnetic regime, we will now relate them to systematically constructed wave functions of excitations of the π-flux ansatz (cf, Fig. 1(a, e)) for the Dirac spin liquid.Our figure of merit to compare ansatz wave functions |ψ GW ⟩ with exact eigenstates |ψ n ⟩ of Eq. ( 1) is the overlap Note that o n is expected to approach zero exponentially fast as N → ∞ due to the orthogonality catastrophe.But on a given system size the behavior of o n as a function of J 2 /J 1 is insightful.The overlap calculations have been performed in the full Hilbert space without employing any symmetry.

A. Parton construction
Ansatz wave functions of QSLs can be systematically constructed using the parton construction [54,55].The original spin operator S i is decomposed into fermionic parton operators f iα , (α =↑, ↓) acting on an enlarged fermionic Hilbert space, where σ = (σ x , σ y , σ z ) denote the Pauli matrices.To impose the single occupancy constraint n iα = f † iα f iα = 1, the partons can be coupled to a dynamical gauge field [55].Rewriting a spin model such as Eq. ( 1) into fermionic operators yields a Hamiltonian with quartic interactions in the parton operators which are coupled to a dynamical gauge field.Applying a Hubbard-Stratonovich transformation and assuming a static gauge field, the resulting parton Hamiltonian is simplified into a mean-field model of the form, where χ i,j denote the hopping amplitudes of the particular mean-field ansatz.The π-flux ansatz of the DSL on the triangular lattice we employ is shown in Fig. 1(e) [20,56].A magnetic flux of π is implemented through the "up" triangles while zero flux is chosen through the "down" triangles.The bandstructure is shown in Fig. 1(a) and features two gapless Dirac cones with linear dispersion, cf.App. C. To explicitly construct numerical ansatz wave functions |ψ⟩, we perform a Gutzwiller projection of parton Slater determinants |ψ free ⟩, The Gutzwiller projection operator imposes exactly the single occupancy constraint.
Past studies focused mostly on variational studies based on the Gutzwiller vacuum, but a systematic and exhaustive study of Gutzwiller projected excited states has not been attempted before.In Ref. [57] a similar basis of Gutzwiller projected excited states was used to model dynamical response functions by projecting the Hamiltonian into the Gutzwiller excitation subspace, while we investigate the overlap of the Gutzwiller subspace with the exact eigenstates.

B. Vacuum as the filled Dirac sea
To construct a canonical ground state ansatz, we fill all single-particle energy levels of the lower band and perform a Gutzwiller projection.In our ansatz, the boundary conditions can be twisted, leading to a shift of the momentum grid in reciprocal space.The twist corresponds to changing the ansatz for the gauge field (i.e. the hopping phases) in mean-field Hamiltonian Eq. ( 4) [58].We find that the energy of the meanfield ansatz is minimized whenever the two Dirac nodes, indicated as red dots in Fig. 1(f), are shifted to the center of three resolved momenta of the finite size cluster.We refer to this choice as the centered boundary conditions, cf.App.E The resolved single particle momenta on the N = 36 site cluster with centered boundary conditions are shown as purple crosses in Fig. 1(f).We find that, also after Gutzwiller projection, the state with centered boundary conditions has lower variational energies than the standard periodic boundary conditions.Since our original model is real, we only consider the real part of the wave function after Gutzwiller projection (the centered Gutzwiller wave function is genuinely complex, while for the more common periodic/antiperiodic choices it would be real), where |ψ center D.S. ⟩ denotes the filled Dirac sea in the centered gauge [59].
We computed overlaps of |ψ vac ⟩ with low-lying energy states [60] using a Krylov technique which avoids calculating all low-lying eigenfunctions explicitly [61], see appendix A for details.Fig. 3(a,d) show substantial overlaps with the ED ground state of up to o 0 = 0.923 at J 2 /J 1 = 0.12.This maximum is attained in the paramagnetic regime and decays in the 120 • Néel ordered and stripy phase.In the stripy phase, the state with maximal overlap is at low energy but distinct from the ground state.

C. Monopoles
In the version of QED 3 we consider, the U (1) gauge field is compact and as a consequence, there exist charge neutral, topological excitations known as monopoles.If these defects proliferate, the DSL will become unstable and will enter a confined phase, in our context e.g. the familiar 120 • magnetic Néel order (for triplet monopoles) or a 12-site unit cell valence-bond solid (VBS) that breaks lattice symmetries (for singlet monopoles) [48].Conversely, in the limit of a large number of flavors N f , monopoles will be present but are irrelevant for the low energy physics.Hence, their precise role for a given N f is quite topical and still unclear.
Using approximately known scaling dimensions, on the triangular lattice the only symmetry-allowed three-monopole operator should be irrelevant [40,48].The quantum numbers of the single monopole operators on the torus are nontrivial and have been obtained in [48,49]: the singlet monopoles are in the X.A irrep while the triplet monopoles have momentum K and are even under reflection and rotation (A1).
Hence, a possible signature of the triplet monopole would be a low-energy mode at momentum K as advocated in some numerics [22,57,[62][63][64] Note however that this "exotic" monopole is evolving into a TOS level in the 120 • phase, and therefore its presence at low energy is not as much as a surprise as for the singlet.In this work, we want to be more precise and we specifically construct a microscopic monopole wave function by fixing a flux pattern such that a ±2π flux is distributed over the whole lattice (i.e. each triangular plaquette has an additional flux ±2π/(2N ) w.r.t. to the π-flux Ansatz), see App.G.By computing the tight-binding dispersion, we find an exact two-fold degenerate zero-energy state.Thus, we can generate several monopole wave functions by filling up all negative energy states and putting 2 electrons (up FIG. 4. Overlaps of low-energy levels at J2/J1 = 0.125 with the the QED3 ansatz wave functions from Gutzwiller projection.For the vacuum and monopole states, we show on, whereas for the bilinear excitations o B 0,1 n is shown.The diameter of the colored circles is proportional to the overlap with the state at the center of the circle.We observe that almost every eigenstate in the dense low-energy ED spectrum has significant overlap with only one of the various excitation ansatz types. or down) in these levels, and then Gutzwiller projecting them.In the end, we can build six singlet and two triplet monopole wave functions from which we can compute exact overlaps with the many-body eigenstates.
Results are presented in Figs.3(b,d) for the triplet and singlet monopole states.Quite remarkably, we do measure a significant overlap with the lowest singlet state found in the irrep X.A and the lowest triplet in the irrep K.A1, in perfect agreement with the expected quantum numbers.Both monopole excitations maximize o n in the paramagnetic regime, with overlaps of up to o n ≈ 0.65 for the singlet and o n ≈ 0.67 for the triplet monopole, Fig. 3(d).Moreover, we have also computed corresponding overlaps for two-monopole states, which are discussed in appendix H.

D. Bilinear fermionic excitations
The next set of excitations of QED 3 we investigate is particle-hole excitations of the Dirac sea represented as bilinears in the fermionic operators.Here, we construct them as particle-hole excitations of the parton ansatz, as shown in Fig. 1(b) with centered boundary conditions.First, we consider excitations of partons in the Dirac sea at single particle momenta closest to the Dirac node.For simplicity, we focus on S z = 0 wave functions at momenta Γ and M , i.e. we can construct 24 states which are all found to be linearly independent after Gutzwiller projection.These excitations are constructed by exciting a parton of either spin σ =↑, ↓ on one of the six momenta closest to the two Dirac nodes, either to the exact same momentum (resulting in a k = Γ or one of the k = M point excitation) or the corresponding momentum at the other Dirac node (resulting in the other two k = M excitations).App.F explains the counting of the degeneracies of the bilinear excitations in further details.
These 24 states span a large subspace which we compare to low-lying excitations of the spin Hamiltonian.To do so, we orthogonalize the states after Gutzwiller projection such that we can compute the overlap of states from ED with this 24 dimensional subspace we name B 0 .To compare ED eigenstates to a subspace B, we define the overlaps o B n , Our definition is equivalent to measuring the norm of the projected state onto the space B, since which yields o B n = 1 if ψ n is an element of the space B, and o B n = 0 if the state is orthogonal.By construction, non-zero overlaps are only found with states at momenta Γ or M due to the momentum conservation of the Gutzwiller projection.An overview of all states with significant overlap is shown as yellow or green symbols in Fig. 4. Prominently, the bilinear excitations have a sizeable overlap of up to o B0 n ≈ 0.74 with the low lying S = 1 Γ.B1 level belonging to the tower of states of the 120 • Néel state along with the S = 1 K.A1 level, which we previously found to have large overlaps with the triplet monopole excitation.Interestingly, also the lowest lying S = 0 level at M.B2 has sizeable overlaps of up to o B0 n ≈ 0.54.Regarding the point-group quantum numbers, some quantitative predictions have been made in Ref. [48]: the bilinear excitations are expected at Γ or M and they are all odd under reflection.We do confirm large overlaps with odd states (Γ.E2 or M .B2), but we also find states which are even under reflection.As Ref. [48] considered the N 2 f = 16 neutral bilinear field operators, the torus geometry with its more complex excitation spectrum explains this discrepancy in the number of bilinear states, as emphasized in Ref. [43].
Furthermore, we have considered bilinear excitations with the minimal momentum transfer distinct from Γ and M , which are degenerate in energy before projection with the already considered set.This yields another set of 48 S z = 0 ansatz wave functions.Again, we find these states to be linearly independent, spanning a space we call B 1 .Many lowlying exact eigenstates are found to have significant overlaps o B1 n with this space, cf.As a remarkable result of this construction, we find that almost the entire complex and dense low-energy spectrum in the paramagnetic regime has significant overlaps, with either the vacuum state, monopoles or bilinear excitations, as is beautifully visible in Fig. 4. Again, the overlaps of the bilinear excitations are maximized in the paramagnetic regime, while sharply dropping off in the stripy phase.Quite interestingly, even in the 120 • Néel phase, some overlaps are still significant.We interpret this as a sign that the 120 • Néel is a natural descendant of the DSL, while the collinear stripy phase has not been identified as an instability of the DSL.

E. Quantum dimer model
While we have a way to construct "model states" for the vacuum, the monopoles, and the fermion bilinear excitations, we are not aware of a correspondingly simple way to construct the gauge field states of QED 3 on the torus.At N f → ∞ the photon modes have been calculated for the square torus [43], and the softest gauge field excitation lies even below the first bilinear level on the square torus.At N f = 4 there are no corresponding results available.
We, therefore, pursue a rather different avenue here, by investigating the energy spectrum of the hardcore quantum dimer model (QDM) on the triangular lattice and comparing it to the J 1 − J 2 Heisenberg model on the same finite-size clusters.There are two different points of view on this procedure.The first one is that there is a long history in frustrated quantum magnetism to investigate QDMs as effective models for the singlet subspace of S = 1/2 Heisenberg Hamiltonians in magnetically disordered phases, such as VBS or spin liquid phases [50,[65][66][67][68].Most of these applications were for square and kagome lattices, but none for the triangular lattice so far [69].The other point of view is to consider the QDM as a quantum link model for a pure gauge theory possibly with static background charges [70].In this context however, the triangular lattice QDM would be expected to describe a Z 2 gauge theory and not a U (1) theory.The tension between these two points of view needs to be clarified in future work in view of our findings below.
The triangular lattice QDM was shown [71] to host an intriguing Z 2 spin liquid for values of V /t ∈ [0.8, 1], i.e in the vicinity of the Rokhsar-Kivelson point V /t = 1.For smaller V /t values a valence bond solid with a 12-site unit cell is realized with an estimated parameter extent in V /t ∈ [−0.75 ± 0.25, 0.8] [71,72].Our own exact diagonalization spectra shown in Fig. 11 suggest that V /t = −1.0 is still within the VBS phase for the considered clusters.For even more negative values of V /t a columnar VBS phase is found.
Interestingly we find that many singlet levels of the low energy spectrum of the Heisenberg model at J 2 /J 1 = 0.125 are surprisingly well reproduced by the low-energy spectrum of the QDM at V /t = −1.0,see Fig. 5(a).We observe that the first excited state in both cases is a singlet at Γ.E2.The next low-lying excitations are two excitations at X 1 .A and X 2 .B (b) for both the QDM and the spin model.These excitations are then followed by two singlet excitations at M.B1 and M.B2.Moreover, we observe a low-lying singlet state at K.E, which cannot be described as a monopole or bilinear excitation, see Fig. 4.These momenta are exactly what is expected for a lozenge VBS, cf.table I.However, we could not match all of the respective point group irreps, which can differ in a spin model, where non-trivial phases of the resonating dimers can occur.This unexpected (but not perfect) similarity of the lowenergy spectrum of the two models raises the possibility that the QED 3 region of the triangular lattice is actually unstable and flows to a confining phase in the IR.The known QDM phase diagram suggests that this confining phase could be a 12-site unit cell VBS.In order to probe for this possibility we have calculated the dimer-dimer correlations of both models in Fig. 5(b,c).While the qualitative agreement between the correlations of the two models is remarkable, in the spin model the correlations decay faster with distance [73].It will be interesting for future research to explore whether there is a small but finite VBS order parameter in the spin liquid region, or whether these correlations are ultimately just the critical VBS fluctuations expected in the QED 3 Dirac spin liquid.

V. DISCUSSION / CONCLUSION
In the above, we revealed a rather compelling one-to-one correspondence between the elementary excitations of QED 3 and the excitations of the triangular J 1 -J 2 Heisenberg model, see Fig. 4.Even though the structure of low-energy excitations in the paramagnetic regime is rather complex, we have demonstrated that monopole and bilinear excitations of the πflux ansatz have comprehensive overlaps with almost all lowlying eigenstates.From this, we can compellingly conclude QED 3 to be the organizing principle of the phases of matter in and close to the paramagnetic regime.Moreover, we have pointed out a close resemblance between the dimer correlations and low-lying energy spectrum of the J 1 -J 2 model in the paramagnetic regime and the quantum dimer model in the valence bond solid phase on both the N = 36 and N = 48 site clusters.This could constitute evidence for a valence bond solid being realized, also found in recent large-scale DMRG simulations pointing out quasi-long range dimer order [32].In light of our findings, we would like to discuss possible scenarios for the phase diagram of the J 1 -J 2 model.
A transition from the 120 • Néel phase to the 12-site valence bond solid appears to be a possible scenario.From a field theoretical perspective, such a transition can be described by a deconfined quantum critical point [74][75][76].As pointed out in Ref. [77] this precise critical point on the triangular lattice would be described by a N f = 4 QED 3 with an emergent P SU (4) = SU (4)/Z 4 symmetry.As a consequence of this enhanced symmetry, an exact degeneracy between the singlet and triplet monopoles in the spectrum should be observed.Indeed, we find the energy levels of the two types of monopoles close to being degenerate at the critical point J 2 /J 1 = 0.09, while a larger splitting between these levels is observed throughout the remaining paramagnetic regime.In the case of an extended DSL region, the same degeneracy between singlet and triplet monopoles would be expected, while our results show growing energy splitting between them as J 2 /J 1 ≳ 0.1.
A second scenario is that the paramagnetic regime indeed realizes a stable Dirac spin liquid phase.Even though we do not observe the expected degeneracy between the two monopole excitations to be realized, we cannot rule out that the remaining energy splitting is still a finite-size effect.Recent work is concerned with the stability of a DSL as a function of N f [37,38,40], and indications exist for the DSL to be a stable phase on the triangular lattice [78].Empirical signatures of bilinear and monopole excitations had also been reported in a previous DMRG study [22].
The data obtained from our exact diagonalization on small clusters does not allow to unambiguously distinguish between these two scenarios.From our data, both of the above scenarios are equally plausible.However, we think that evidence of either critical or long-range lozenge VBS correlations constitutes a new insight into the physics of the intermediate paramagnetic regime.We would like to point out, that a recent study has found evidence for a similar VBS state to be stabilized upon coupling the triangular lattice to phonons [79].
Finally, from previous work, we also notice that a chiral spin liquid (CSL) is in close vicinity to the paramagnetic regime when adding small scalar chirality interactions on the triangles of the form S i •(S j ×S k ) [9].A DSL can be unstable towards a magnetically ordered phase or a CSL phase through a QED3-Gross-Neveu quantum critical point [80].While we do not observe direct evidence of a gapped CSL being realized, such a transition could be in close vicinity to the paramagnetic regime.Another point to clarify is whether the fact that the lowest parton ground state is reached for the centered boundary conditions with its complex (i.e.time-reversal symmetry breaking) nature, indicates an instability towards a spontaneous chiral spin liquid formation.This effect would not appear in other geometries, such as the square torus studied in Ref. [43], where the time-reversal invariant antiperiodic boundary conditions minimize the energy.
While in our study we directly compared ansatz wave functions with eigenstates from ED on a finite cluster, we would like to point out that this technique could also be applied in complementary numerical techniques, for example, to construct excited states in tensor network simulations [81][82][83] or other variational Monte Carlo techniques.This would allow also these techniques to study the non-trivial monopole and bilinear excitations revealed in the above on cylinder geometries or larger periodic tori.An application of the technology developed in this work to the kagome lattice S = 1/2 Heisenberg model is the natural next step in the quest for quantum spin liquids.the main manuscript, the centered boundary conditions with a shift vector, have been used, centering the Dirac nodes on the left-pointing triangles.When using a square torus, the grid in momentum space is also a square lattice, and then the analogous shift vector procedure would result in antiperiodic boundary conditions, centering the Dirac nodes in a square plaquette in momentum space.

Appendix F: Construction of bilinear excitations
Bilinear excitations are constructed by performing particlehole excitations out of the filled Dirac sea with centered boundary conditions.Gutzwiller projection conserves momentum and hence the momentum of these bilinear excitations can be derived from the momentum of the parton wave function.We only focus on excitations created from removing one parton from the lower band and inserting it into the upper band in one of the six degenerate energy levels with minimal energy.When exciting a parton from one momentum in the lower band to the same momentum in the upper band a net k = Γ = (0, 0) state is obtained.If a parton is excited from a momentum in one valley to the corresponding momentum in the other valley, a net k = M state is obtained.Further excitations with momenta close to the k = Γ and k = M can be obtained by exciting a parton to the other four momenta.We summarize these possibilities in Fig. 8.

Appendix G: Construction of monopole excitations
As sketched in Fig. 1(c), putting a monopole inside the torus amounts to having a single flux quantum threading the whole system.By inspection, we have found one possible solution on the N = 36 cluster that we consider, see Fig. 9(a) where we specify the phase factors appearing in the tightbinding model.Such a phase pattern provides a uniform flux of Φ = 2π/(2N ) per plaquette (on top of the π-flux Ansatz required for the DSL which amounts to having π flux on all up triangles and zero flux on the down ones).The corresponding eigenspectrum of the tight-binding model is also given in Fig. 9(b), showing the existence of two zero-energy modes.Indeed, if the Dirac points are available for Φ = 0, then it is known that the monopole spectrum will exhibit N f degenerate states at zero energy (as found for instance in the quantum Hall effect of graphene).Now in order to construct spin wave functions, we need to fill up all negative energy states and also put two fermions in the zero-energy modes, thus leading to 3 singlet and one triplet states [49].Similarly, we can consider "antimonopoles", i.e. a negative flux Φ = −2π/(2N ) in order to get six additional states.In total, after Gutzwiller projection, we can generate six singlet and two triplet monopole wave functions.We have checked that they are linearly independent after projection.

Appendix H: Construction and overlaps of two-monopoles
In addition to the overlaps obtained from the single monopole excitations presented in the main text, we have also performed overlap calculations with ansatz wave functions for a two-monopole excitation at J 2 /J 1 = 0.125.The results are shown in Fig. 10.These states are on average higher in energy than the single monopole and bilinear excitations.However, several low-lying states with energies E < J 1 still have significant overlap with these two-monopole excitations.The ansatz for these excitations is created analogously to the single-monopole excitations in Fig. 9, where instead of having 2π flux through the system, we choose 4π flux.
Note that in such a situation, the tight-binding model exhibits 4 modes at approximately zero energy so that we can construct 2 × 6 2 = 72 wave functions with a total S z = 0 for a positive or negative flux.

FIG. 1 .
FIG. 1.Quantum electrodynamics in 2 + 1 dimensions implemented on the triangular lattice.(a) The QED3 vacuum is constructed by filling the lower band of the π-flux ansatz shown in (e) and performing a Gutzwiller projection.The red dots indicate the Dirac nodes with linear dispersion.(b) Bilinear excitations are obtained by creating particle-hole excitations on the Dirac-sea.(c) Monopole excitations are created by distributing a unit of magnetic flux across the full torus.(d) Gauge field (photon) excitations are studied using an effective quantum dimer model.(e) The π-flux ansatz is implemented on a two-site unit cell and features π-flux through up-triangles and 0-flux through down-triangles.An example gauge is shown with the dashed and solid lines which indicate a phase of π and 0 on the hopping, respectively.(f) Brillouin zone (BZ) of the triangular lattice and folded BZ of the π-flux ansatz.The high symmetry points Γ, K, M, X are marked with orange symbols.The grey stars show the momentum resolution of the N = 36 simulation cluster.The purple crosses indicate the shifted momenta minimizing the kinetic energy of the filled Dirac sea.

FIG. 2 .
FIG. 2. Energy spectrum of the J1-J2 model on the triangular lattice for N = 36 as a function J2/J1 in (a) and for N = 48 at J2/J1 = 0.125 in (b).Symbol colors and shapes indicate the momentum and point group representation of the eigenstate.Filled (open) symbols indicate even (odd) total spin S. The 120 • Néel state persists for J2/J1 ≲ 0.09, the paramagnetic regime is realized for 0.09 ≲ J2/J1 ≲ 0.14, and the stripy magnet is stabilized for J2/J1 ≳ 0.14.The vertical dashed line at J2/J1 = 0.125 in panel (a) highlights the spectrum in order to compare with the N = 48 data shown in panel (b).

FIG. 3 .
FIG. 3. Overlaps of low-energy levels of the J1-J2 model with various ansatz wave functions constructed out of various excitations above the QED3 vacuum.For the vacuum and monopole states, we show on, whereas for the bilinear excitations o B 0,1 n is shown.The diameter of the colored circles is proportional to the overlap with the state at the center of the circle.(a) The vacuum state has significant overlap (up to ≈ 0.92) with the ground state.(b) The singlet monopole has significant overlap with the low-lying X.A level (≈ 0.65) while the triplet monopole has sizeable overlap with the low-lying triplet K.A1 (≈ 0.69) level.(c) Among the bilinear-excitations B0 there is significant overlap with the low-lying Γ.B1 and M.B2 levels.The Y0.B energy level has significant overlaps (≈ 0.67) with the bilinear excitations B1 (d) numerical values of the maximal overlap of the aforementioned states with eigenstates from ED.For every excitation, the maximum is attained in the paramagnetic regime.

Fig. 4 .
Most notably, the low-lying Y 0 .B state has overlaps of up to o B1 n ≈ 0.67.We would like to emphasize the large number of low-lying bilinear excitations, 24 (k = Γ, M ) + 48 (k ̸ = Γ, M ) = 72 states with S z = 0 and another 72 states with total S z = ±1 expected for the DSL on the torus.Thus, there are a total of 144 lowlying bilinear excitations of the parton ansatz, cf.appendix F. The described counting is valid, whenever each Dirac point has three symmetric neighboring momenta, whose quasiparticle energy bands have exactly the same energy.This holds for both the N = 36 and N = 48 site clusters in this manuscript as well as more generically 6N × 6N clusters for N ≥ 1.

FIG. 5 .
FIG. 5. (a) Comparison of the low-lying energy spectrum of the J1-J2 model at J2/J1 = 0.125 and the QDM model at V /t = −1.0 on the N = 48 cluster.The filled (open) circles denote even (odd) spin levels.The QDM is only expected to reproduce filled levels because it cannot describe states with total spin S > 0. (b) Connected dimer correlations ⟨(S z 0 S z 1 )(S z i S z j )⟩c of the ground state from ED of the J1-J2 model at J2/J1 = 0.125 on the N = 48 cluster.(c) Connected dimer correlations of the QDM in the VBS phase at V /t = −1.0.We observe a close resemblance between these patterns.

FIG. 7 .
FIG. 7. Comparison between momentum resolution with periodic and centered boundary conditions for the N = 36 site simulation cluster.In the centered boundary conditions, the position of the Dirac nodes shown as red dots is at maximal distance from all resolved momenta.

deg: 6 (FIG. 8 .
FIG. 8. Bilinear excitations as particle-hole excitations of the parton ansatz with degeneracies.(a) excitation with k = 0 by exciting one parton to the upper band with exactly the same momentum.(b) excitation with k = M by exciting one parton to the corresponding momentum at the Dirac valley.(c,d) excitation by moving one parton to another momentum either at the same (c) valley or the other valley(d).

FIG. 9 .
FIG. 9. Top: flux pattern chosen to distribute one flux quantum over the N = 36 triangular lattice so that all triangles contain the same additional flux Φ wrt the π flux Ansatz.By convention, all fluxes point upwards or to the right.Bottom: corresponding tight-binding dispersion energies.

FIG. 11 .
FIG. 11.Excitation energies of the quantum dimer model at different momemta k = Γ, M, K, X.The top row shows the spectra for the topologically trivial sector while the lower row shows the spectra for the other three topological sectors.Filled and empty symbols indicate the data from the N = 36 and N = 48 site clusters, respectively.

FIG. 12 .
FIG. 12. (a) Comparison of the low-lying energy spectrum of the J1-J2 model at J2/J1 = 0.125 and the QDM model at V /t = −1.0 on the N = 36 cluster.The filled (open) circles denote even (odd) spin levels.(b) Connected dimer correlations ⟨(S0 • S1)(Si • Sj)⟩c of the ground state from ED of the J1-J2 model at J2/J1 = 0.125 on the N = 36 cluster.(c) Connected dimer correlations of the QDM in the VBS phase at V /t = −1.0.This figure is the analog of Fig. 5 with N = 36 in the main text.