Quantum Spin Supersolid as a precursory Dirac Spin Liquid in a Triangular Lattice Antiferromagnet

Based on the recent experiments on the triangular lattice antiferromagnet Na$_2$BaCo(PO$_4$)$_2$, we propose the easy-axis XXZ spin-1/2 model on the triangular lattice, that exhibits a quantum spin supersolid, to be a precursory Dirac spin liquid. Despite the presence of a three-sublattice magnetic order as a spin supersolid, we suggest that this system is close to a Dirac spin liquid by exploring its spectroscopic response. The physical consequence is examined from the spectroscopic response, and we establish the continuous spectra near the M point in addition to the K point excitation from the spinon continuum on top of the three-sublattice order. Moreover, the satellite peaks were predicted at the mid-points connecting the $\Gamma$ and K points. This proposal offers a plausible understanding of the recent inelastic neutron scattering measurement in Na$_2$BaCo(PO$_4$)$_2$ and could inspire further research in relevant models and materials, such as K$_2$Co(SeO$_3$)$_2$ and Rb$_2$Co(SeO$_3$)$_2$, and even more anisotropic magnets like PrMgAl$_{11}$O$_{19}$.


I. INTRODUCTION
The U(1) Dirac spin liquid (DSL) has been viewed as the mother state of the competing orders such as the Néel and valence bond solid states in 2D [1].By studying the symmetry properties of the monopole operators of the emergent quantum electrodynamics and ensuring the concomitant of the ordering-driven spontaneous mass generation and the confinement, recent theoretical work suggested the increased stability of the U(1) DSL on frustrated lattices [2].It inspires that the signatures of DSL may occur even in the nearby ordered phases on frustrated lattices.In this work, we attempt to identify the specific material and model that may reveal the signatures of DSL and competing orders on the triangular lattice.
Recent experiments on the triangular lattice antiferromagnet Na 2 BaCo(PO 4 ) 2 show several pieces of puzzling phenomena [3].While the system is clearly ordered antiferromagnetically with a three-sublattice order of the spin supersolid type, the magnetic excitations do not behave like typical spin-wave modes.Instead, broad continuous magnetic excitations appear in both the momentum and energy domains.Moreover, there exist strong continuous excitations at low energies around the M points that seem to be quite unexpected for an antiferromagnetic order at the K point.We interpret these continuous excitations as the signatures for the precursory U(1) DSL.
For the more material's aspect of the physics, it is the conventional belief that the 3d transition metal oxides should have more symmetric spin models.In reality, all the 3d transition metal ions have rather different physical properties, exhibiting rich behaviors with a diverse of models that depend sensitively on the structures and the crystal field environments [4][5][6][7].In many cases, the orbitals get involved and wildly change the physics.For the Co ion, the often-discussed scenario is actually about the selection of the high or low spin states due to the competition between the Hund's coupling and the crystal electric field splitting.More recently, however, the spin-orbit coupling of the Co ion has received quite some attention.This occurs when the electrons partially fill the t 2g shell [8].Due to the active spinorbit coupling and the orbitals, the local moment of the relevant Co ion is no longer pure spins, and is instead an entangled object of spin and orbital.The resulting effective model is expected to be anisotropic instead of Heisenberg-like.The well-known limiting examples were the quasi-1D Ising magnets CoNb 2 O 6 , BaCo 2 V 2 O 8 and SrCo 2 V 2 O 8 [9][10][11][12][13][14].For the same reason, the Kitaev interaction has been extended to the honeycomb cobaltates including Na 2 Co 2 TeO 6 and Na 2 Co 2 SbO 6 [15][16][17], and a similar form of Kitaev interaction is re-invoked for CoNb 2 O 6 [18,19].For our motivation of the triangular lattice antiferromagnet Na 2 BaCo(PO 4 ) 2 and other Co-based triangular lattice compounds, the model is expected to be the spin-1/2 easy-axis XXZ model [3,20].
Na 2 BaCo(PO 4 ) 2 develops a three-sublattice antiferromagnetic order below about 0.15K.The weak exchange couplings allow the external magnetic fields about 2T to polarize the spins for both in-plane and out-of-plane directions [21].The magnetic excitations of these polarized states are quantitatively captured by the linear spin-wave theory based on the spin-1/2 XXZ model with strong fields.The ratio of the exchange coupling for the xy components over the z-component coupling is ∼ 0.58, and both are antiferromagnetic and thus fullyfrustrated.The spin-1/2 XXZ model on the triangular lattice is a well-studied problem.The phase diagrams with and without the field are quite clear.The supersolid order, that breaks both the U(1) symmetry via the in-plane spin order and the lattice translation via the z-component order, was known to be present in the ferromagnetic xy and the antiferromagnetic z regime, and was later shown to persist to the fully-frustrated regime with both antiferromagnetic xy and z couplings until the Heisenberg point [22,23] and is entirely due to quantum origin.Therefore, the antiferromagnetic order in Na 2 BaCo(PO 4 ) 2 at zero field is a quantum spin supersolid at low temperatures.Due to the U(1) symmetry and the easy-axis coupling, the intermediate field along the z direction stabilizes a 1/3 magnetization plateau with a 'UUD' ("up-up-down") spin configuration in Na 2 BaCo(PO 4 ) 2 .This is consistent with the expectation of the XXZ model.What is not understood from the model is the presence of large continuous magnetic excitations at zero field.This feature is absent for Ba 3 CoSb 2 O 9 that is described by the same model but in the easy-plane regime [24,25].This difference is expected to arise from the enhanced frustration in the easy-axis regime of the XXZ model.
To understand the puzzling experiments in the triangular lattice antiferromagnet Na 2 BaCo(PO 4 ) 2 , we propose an interpretation from the precursory U(1) DSL within the framework of a mean-field theory.This mean-field theory synthesizes both the supersolid antiferromagnetic orders and the fractionalized excitations in the DSL.This mixed type of mean-field theory has been invoked for the proposed AF * state in the high-temperature cuprates, where the AF * state is a Z 2 topologically ordered state with the fermion bilinear antiferromagnetic order on the square lattice.The " * " in AF * refers to the Z 2 topological order with deconfinement and fractionalization in the context of Ref. 26.At the level of the mean-field analysis, the treatment and the spirit are quite similar.Our analysis in this work will be restricted to the meanfield theory.We consider the π-flux U(1) DSL on the triangular lattice with the supersolid antiferromagnetic order.At the mean-field level, we are able to reproduce the continuous excitations around the M and K points that were observed by the experiments.In addition, as the predicting power of the mean-field analysis, we further predict that, the satellite peaks of the dynamic spin structure factor appear at the midpoints between the Γ points and the K points.
The remaining parts of the paper are organized as follows.In Sec.II, we introduce the model Hamiltonian and set up the mean-field theory for both supersolid order and the spinons.In Sec.III, we analyze the mean-field theory for the spinons and the supersolid order.In Sec.IV, the spectroscopic properties for the spin dynamics are computed from the mean-field theory.Finally, in Sec.V, we discuss various aspects associated with the triangular lattice easy-axis quantum magnets.In the appendices, various details of the calculation from the main text are explained.

II. MODEL HAMILTONIAN AND MEAN-FIELD SETUP
The triangular lattice antiferromagnet Na 2 BaCo(PO 4 ) 2 is modeled as the spin-1/2 XXZ model with [21,27], where J ⊥ , J z > 0 are both antiferromagnetic and the last term is the Zeeman effect from an out-of-plane magnetic field B. With an easy-axis anisotropy 0 < J ⊥ /J z < 1, one finds that the magnetic phase undergoes through the Y-shape, UUD, and V-shape spin configurations with a √ 3 × √ 3 three-sublattice structure with increasing fields [28], and eventually becomes a fully polarized state above the saturation field The spin orientation on each sublattice is depicted for different ordered phases in Fig. 1.Throughout this work, we focus on the zero-field Y-shape state.Because the state breaks the U(1) symmetry and the lattice transition at the same time, this state is also known as the spin supersolid [29,30].The spin-wave study based on the spin supersolid cannot produce the low-energy continuous magnetic excitations near the M points, indicating the insufficiency of the spin wave theory for the spin dynamics.
We aim to find a phenomenological understanding of the zero-field magnetic state that contains the supersolid order and at the same time provides the continuous excitations near both the K and M points.Based on the recent progress about the U(1) DSL and the competing orders [2], we thus consider the possibility of the precursory DSL for Na 2 BaCo(PO 4 ) 2 and the related easy-axis XXZ model.To describe this exotic scenario, we first use the Abrikosov fermionic spinon representation for the spin operator with and the Hilbert space constraint where f † iα (f iα ) creates (annihilates) a fermionic spinon at the site i with the spin α.The supersolid antiferromagnetic order is captured by the fermionic spinon bilinear in ⟨S i ⟩.We use the Weiss mean-field theory to incorporate the supersolid order and rely on the slave-fermion mean-field theory to capture the U(1) DSL.For the first one, one expresses S µ i S µ j as ⟨S µ i ⟩S µ j where ⟨S µ i ⟩ captures the supersolid antiferromagnetic order.The continuous spectra in the dynamic spin structure factor would be captured by the particle-hole continuum of the fermionic spinons of the DSL at least at the mean-field level.It is known that, the π-flux DSL naturally gives the continuous excitations near the M point [31].
Combining these two points, we propose the following mean-field Hamiltonian for B = 0, where we fix the gauge with t ij = ±t (see Fig. 1(b)) such that a π-flux through each unit cell is presented, and we have used J⊥ and Jz to distinguish them from the original bare couplings.The chemical potential µ is introduced to enforce the fermion number occupation such that 1/2filling is obtained on average.The spinon hopping parameter t can be regarded as the energy unit.While the chemical potential µ can be fixed by the fermion filling, the parameters J⊥ and Jz can be determined selfconsistently by making sure that the mean-field theory yields the given antiferromagnetic supersolid orders.The procedure is briefly described here.From the mean-field Hamiltonian, one can obtain the spinon meanfield ground state by filling the fermionic spinons up to the spinon Fermi energy (E F ), and the ground state is given as where ϵ j refers to the spinon energy for the j-th spinon.
In reality, one should label the spinon with the crystal momentum, the band/sublattice indices, and the spin flavor.One can further perform the Gutzwiller projection to further improve the wavefunction.Here, we are content with the mean-field treatment and aim for a phenomenological description.The magnetic order is then computed as where the left are the input order parameters, and the right are the expectation values of the corresponding spin operators with respect to the mean-field ground state.
The operator hat symbol is introduced explicitly to be distinguished from the left order parameters.With the global U(1) symmetry, the above equations can be reduced to two equations with one for the in-plane component and the other for the z component.The two parameters J⊥ and Jz can thus be determined.Thus, at the mean-field level H MF describes a π-flux U(1) DSL with the antiferromagnetic supersolid order, i.e. the precursory DSL in the terminology of Sec.I.
The original XXZ spin model can be obtained from H MF by imposing the local Hilbert space constraint of having one spinon per site that leads to precisely the Hilbert space of the microscopic spin model.This can be effectively achieved by introducing a strong onsite Hubbard-U -like interaction.The spin model can be effectively recovered by doing a second-order perturbation of the spinon hopping [32].The resulting model is of the XXZ form with J ⊥ = J⊥ + 4t 2 /U and J z = Jz + 4t 2 /U .

III. MEAN-FIELD THEORY FOR THE SPINONS AND THE SUPERSOLID ORDER
In the above mean-field setup, we have schematically separated the spin into "two parts".One part is respon-sible for the antiferromagnetic supersolid order, and the other part is treated as the spinon bilinears.Each separately would not able to sufficiently capture the continuous excitations that were observed experimentally.We need these two ingredients together to modify the spin dynamics.It is illuminating to write down the spin correlation where the expectation is taken with respect to the precursory DSL that is described by the mean-field theory in the previous section.Here the first term gives the spin supersolid order, and the second term has the spin correlation from the spinons of the DSL.
In fact, separately computing the z-component and xycomponent spin correlations for the nearest neighbors like Eq. ( 9), one could obtain a variational ground state energy E var .One could optimize the variational energy for the XXZ model against the two variational parameters Jz and J⊥ .The procedure of this approach is slightly different from what has been described in the previous section.Here, one does not need the given order parameters that are obtained from the experiments or the numerical calculations of the XXZ model.

A. Spinon bands from mean-field theory
To study this precursory DSL, we numerically evaluate the spinon band structures of the mean-field Hamiltonian Eq. ( 4) with the supersolid ansatz.As the system has a U(1) degeneracy along the z-axis, without loss of generality, we choose a specific orientation such that ⟨S y i ⟩ = 0 (see Fig. 1).Now there remain three degrees of freedom to the ansatz, which we denote as ⟨S z 0 ⟩ (the z component of the blue sublattice), ⟨S z 2 ⟩ and ⟨S x 2 ⟩ (the z, x components of the green sublattice).We obtain the following mean-field Hamiltonian, where ⟨S µ i ⟩ is extracted from the experiments [20].More details can be found in Appendix A. The Yshape state is for the supersolid, where (M z 0 , M z 2 , M x 2 ) = (−0.54,0.27, 0.173)µ B /Co.To convert this M configuration into our spin order, we notice that the in-plane saturation moment is M x S ≈ 2.2µ B /Co, and the out of plane saturation moment is M z S ≈ 2.5µ B /Co [20].Therefore we obtain ⟨S i ⟩ according to = (−0.108,0.054, 0.0393).( 11) Here we do not plan to be quantitative.As it is a spinon mean-field theory without Gutzwiller projection, quantitatively reproducing experiments is impractical.But to be specific, we use the above order parameters in the calculation.
The spinon hopping can be deduced by comparing the Dirac velocity and/or the bandwidth of our INS calculation against the experimental ones [3], and we here choose t ≈ 0.0338 meV.The effective exchange J is obtained self-consistently.Namely, we vary J⊥ and Jz to make ⟨ Ŝx ⟩ and ⟨ Ŝz ⟩ to be consistent with the input experimental ones.Since ⟨S z 0 ⟩ = −2⟨S z 2 ⟩ always holds true (for both experiment and our evaluation), there are only two degrees of freedom to vary.The result is Jz = 4.394t, J⊥ = 4.102t.
The self-consistent spinon bands are obtained and shown in Fig. 2. The magnetic unit cell is six times larger than the lattice primitive cell due to the supersolid order and the π-flux background.The supersolid order triples the unit cell, and the π further doubles the unit cell.Thus, with the spin quantum number, there are 12 spinon bands in the magnetic Brillouin zone (MBZ) [the orange box in Fig. 2(c)], which is 1/6 of the lattice Brillouin zone.In the MBZ, four Dirac cones (2 valleys + 2 spins) are separated by a momentum around M. In Fig. 2(b), we further find the Dirac cones open up a small gap ∼ Jx ⟨S x ⟩.In the next subsection, we separate the Y-shape order into the Ising order along the z-axis and the planar order in xy-plane, and find that this small gap is opened by spin-up/down mixing due to the supersolid order.Since the time reversal symmetry is broken, we have made sure that the filled spinon bands have a vanishing Chern number in Appendix B, and thus there is no Chern-Simons field description or edge mode like the chiral spin liquid.

B. Effects of the supersolid order
In the absence of the spin supersolid order ⟨S z 0 ⟩ = ⟨S z 2 ⟩ = ⟨S x 2 ⟩ = 0, the mean-field Hamiltonian in the previous subsection gives rise to the gapless Dirac cones at the spinon Fermi level.There are two Dirac cones for each spin in the MBZ as shown in Fig. 3.The spinon bands are colored by ⟨ Ŝz (k)⟩ = ⟨ψ n,k | Ŝz |ψ n,k ⟩, where |ψ n,k ⟩ is the wavefunction for the n-th spinon band, Ŝz is the z-component of the spin operator and can be expressed as I 6×6 ⊗ 1 2 σ z in the basis of X k .To explore the properties of the Dirac spinons, in a usual treatment, one analyzes the symmetry properties of various mass gap terms for the Dirac cones and establishes the connection between the mass gaps with the physical spin observables.Here we already have the magnetic orders of the system in the beginning.We can then directly compute the evolution of the spinon spectrum in the presence of the spin supersolid.With the supersolid order, we have ⟨S z 0 ⟩ + 2⟨S z 2 ⟩ = 0 in Eq. ( 10).The supersolid order carries both the Ising antiferromagnetic order in z and the in-plane xy antiferromagnetic order, where the former is often referred to as the density wave order and the latter is referred to as the superfluid order in the hardcore boson language for the spins.Due to the presence of two such orders, it is a bit illuminating to separately study the impact on the gapless Dirac cones one by one.We begin with the Ising order by setting ⟨S x ⟩ = 0.As it is shown in Fig. 4, the degeneracy in Fig. 3 is lifted by Jz (⟨S z 0 ⟩ − ⟨S z 2 ⟩) ̸ = 0. Since the spin-z is still a good quantum number to label the spinon bands, the two cones from different spin sectors do not mix with each other, and the Dirac cones at the Fermi level remain gapless except for the breaking of the original four-fold degeneracy (see Fig. 4).In the presence of both Ising and in-plane orders for the spin supersolid, the J⊥ ⟨S x 2 ⟩ term does not commute with Ŝz , and this in-plane order immediately creates the gap for the Dirac spinons at the Fermi level by the spin mixing as shown in Fig. 5.

IV. SPIN DYNAMICS
To connect our results with the inelastic neutron scattering experiments, we further calculate the spinon continuum via the transverse dynamic spin structure factor, that is proportional to the inelastic neutron spectrum differential cross-section, as where |Ω⟩ (|n⟩) and E 0 (E n ) refer to the ground (excited) state and its energy, respectively.At the mean-field level, |Ω⟩ is to simply fill the spinon bands below the Fermi energy, and S + p = k f † k−p↑ f k↓ excites one spinon particlehole pair across the Fermi level.The summation over n includes all such excited spin-1 pairs.Since the neutron scattering events involve pairs of spinons, the energy and momentum transfer of one neutron are conserved by the total energy and total momentum of two spinons.This leads to a continuous spectrum in the INS measurements, that represents the fractionalized excitations [31,33,34] and differs from the magnon excitation with integer spins.
In Fig. 6, we show the continuum along the high symmetric momentum path.There are two additional features that arise here.First, the expected signals appear at M in addition to K at low energies, which agrees with the experimental results.Second, the model predicts the satellite peaks at the middle point from Γ and K. Fig. 6 shows the INS signals in the whole Brillouin zone with a fixed energy slice, including the ones around {Γ, M 1,2,3 } and {K, K ′ , K/2, K ′ /2}.These main features persist even with a weak external magnetic field as shown in Appendix C.
The position of these signals is a direct consequence of the coexistence of the DSL and the supersolid order.At low energies, the spinon particle-hole pairs can be created only across the Fermi energy among those gapped Dirac cones.Therefore, the allowed momentum transfer p's in Eq. ( 12) are about the separation of Dirac cones.As shown in Fig. 6, the distribution of the Dirac cones indeed matches the pattern of the INS signal in the momentum space.The spectral continuum at {K, K ′ , K/2, K ′ /2} is then a consequence due to the combination of the Dirac spinons and the supersolid order that has an ordering wavevector at K, K ′ .Since our theory is based on the free spinons and the Hilbert space constraint is not imposed on each site, the intra-Dirac-cone processes that contribute to the Γ point should not be favored by the antiferromagnetic spin interaction, and would in principle be suppressed.This can be shown with a random phase approximation to incorporate the spinon interactions [35].This is the caveat of the free spinon theory [36].
Although we have obtained the continuous excitation FIG. 7. The spin-wave excitation from linear spin-wave theory for the supersolid state.For the Y-shape configuration, the tilting angle θ from ẑ on the red and green sublattices (see Fig. 1) is determined by minimizing the classical exchange energy with cos θ = Jz/(Jz + J ⊥ ) and Jz = 1.73J ⊥ .Although this classical state is a ferrimagnet with a finite magnetization and differs slightly from the quantum supersolid, the spin-wave excitation should be a reasonable qualitative approximation of the spin-wave part of the excitations for the quantum supersolid.
from spinons, the fluctuation of the order parameters can still be important, and these are spin-wave excitations as magnons.These magnons can also contribute to the INS spectrum, but not usually in terms of continuous excitation.The magnons contribute to the INS spectrum as a well-defined quasiparticle-like dispersion.We calculate the spin-wave dispersion with the linear spin wave theory, and the results are depicted in Fig. 7.The gapless spectrum is due to the breaking of the U(1) symmetry.The spectrum at K is understood by the supersolid order with a wavevector K, and there is a visible gap at M point.In reality, the magnetic order is strongly renormalized by the quantum correction [20], and the bandwidth of the spin-wave excitation should be suppressed compared to the linear spin-wave calculation.Taking together with the spinon continuum, the actual INS signal should be a combination of both spinon continuum and the magnonic excitations [26,37].

V. DISCUSSION
There are two major questions that are addressed in this work.One is whether the continuous excitation in Na 2 BaCo(PO 4 ) 2 arises from the spinon pairs of any sort.The other question is whether Na 2 BaCo(PO 4 ) 2 has anything to do with the DSL.We provide a candidate answer to both questions by proposing a precursory DSL that combines the DSL and the supersolidity.It turns out that, our treatment is not really new.As we have recently learned, a similar approach was actually used by Ref. 38 to establish the supersolid spin order, instead of the spin dynamics.The progress in our work is to access the spin dynamics with this approach.
At the mean-field level we have obtained twospinon continuum and found concentrated signals at {M, K, K/2} points.These results provide a possible explanation for the bizarre INS results in Na 2 BaCo(PO 4 ) 2 .For the first question, phenomenologically two-magnon continuous excitations could also contribute to the INS continuum [39].This magnon continuum, however, usually occurs at higher energies instead of the low energy, and its spectral intensity is usually weak.The recent work in Ref. 3 actually calculated the two-magnon contribution and did not find an agreement with As for other candidate states, one natural candidate seems to be a Z 2 spin liquid with condensed bosonic spinons that generate the supersolid order.We find it difficult to reconcile the magnetic structure and the position of the continuous excitations of the Z 2 spin liquids that were classified with Schwinger bosons [40].Phenomenologically, however, it is possible that, if one directly introduces the supersolid on top of the π-flux Z 2 spin liquid in Ref. 40, one may establish a qualitatively consistent result with the experiment.We do not rule out this possibility here.Since the model for Na 2 BaCo(PO 4 ) 2 is the easy-axis XXZ model, we expect more advanced numerical techniques to be used for the spin dynamics.In fact, the J 1 -J 2 Heisenberg model on the triangular lattice was numerically found to stabilize the U(1) DSL with a very weak J 2 .The easy-axis XXZ model is more frustrated than the nearest-neighbor Heisenberg model, and is likely to be more close to the DSL.
Our free-spinon mean-field theory would certainly be renormalized once the Hilbert space constraint is implemented.This can be achieved by the renormalized meanfield theory or by the Gutzwiller projection [41].Except for the modulation and suppression of the spectral intensities, the qualitative results for the spinon continuum along with the proximate spin liquid state should persist, so we do not invoke such methods as we do not focus on the precise spectrum calculations in this study.With the supersolid, the precursory DSL develops a small gap, and this can potentially cause the confinement of the U(1) gauge theory.Since the U(1) DSL is regarded as the mother state for 2D magnets [2], one could imagine a weak spinon pairing that further stabilizes fractionalization [42].Such a state could be equivalent to a π-flux Z 2 spin liquid from the Schwinger boson classification [40].Nevertheless, as long as the fractionalization is present below the confinement length, one would expect continuous excitations at finite energies.Moreover, thermal fluctuation suppresses the supersolid above 0.15K, and our theory expects the system to behave like a gapless DSL.In fact, a finite thermal conductivity κ/T was observed at low temperatures above 0.15K [21] and is consistent with the expectation from DSL [43].
Beyond Na 2 BaCo(PO 4 ) 2 , other Co-based triangular lattice antiferromagnetic materials with the easy-axis XXZ model description could exhibit similar physics, i.e. the coexistence of the magnetic order and the spinon continuum (especially around the M points).In fact, a recent inelastic neutron scattering measurement in another Co-based triangular lattice antiferromagnet K 2 Co(SeO 3 ) 2 also found quite similar structures of continuous excitations, and this system is clearly in the easyaxis regime of the XXZ model with the supersolid order [44].Moreover, this measurement discovered the lowenergy satellite peaks at K/2 and K ′ /2 of the continuous excitations, and this is consistent with our theoretical expectation.The isostructural material Rb 2 Co(SeO 3 ) 2 is expected to behave quite similarly.
Ref. 2 has suggested searching for the precursory DSL among 2D frustrated lattices.One related such system that behaves like a DSL is the spin liquid candidate PrZnAl 11 O 19 where the Pr ions form a triangular lattice with spin-1/2 local moments [45,46].For these triangular lattice antiferromagnets with the non-Kramers doublets, the recent model only slightly deviates from the XXZ model by allowing one extra anisotropic interaction [47].Continuous excitations have actually been observed in PrZnAl 11 O 19 [48].In addition to the frustrated lattices, the bipartite lattices with frustrated interactions can still be promising candidates, as recently proposed [37,49,50].On more theoretical side, the fermionized vortex theory was used to describe the DSL on the triangular lattice [51] and may encounter some issues with imposing time reversal symmetry.Since the time reversal is broken by the supersolid, it might be interesting to revisit the fermionized vortex theory to explore the mass generation and supersolid in the DSL.An alternative approach is to directly fermionize the hardcore bosons that represent the spins [52].This approach is more direct, and can be more useful to obtain the variational wavefunction for the ground state.Numerically, it would be nice to directly access the DSL with the XXZ J 1 -J 2 model.With this ansatz, the spinon hopping part of the mean-field theory can be expressed as where and we have used the fact that 2M = b 2 .
Since the U(1) DSL ansatz gives rise to a π-flux background and doubling of the spinon unit cell, along with the presence of the supersolid ordering, the size of the magnetic Brillouin zone (MBZ) is 1/6 as large as the size of the lattice Brillouin zone in Fig. 8.In order to obtain the dispersion of the spinon excitations, we need to fold the lattice Brillouin zone (BZ) into MBZ and formally convert the sum over k ∈ BZ in Eqs.(A6) and (A7) into the summation over k ′ ∈ MBZ as where k ′ ∈ MBZ is implicitly assumed on the right hand side.Using the fact that 3K = b 1 − b 2 and 2M = b 2 , we can express the whole mean-field Hamiltonian in MBZ as T , and where the blank entries refer to 0, and Here I 2 and I 4 are 2×2 and 4×4 identity matrix, respectively.By diagonalizing Eq. (A9), we obtain 12 spinon bands in the MBZ as shown in the main text.

Appendix B: Vanishing Chern number of filled spinon bands
With the quantum supersolid order, the Dirac cones at the Fermi level acquire the mass gap, and the 6 filled spinon bands under the Fermi level are well separated from the upper 6 spinon bands.As the time-reversal symmetry is already broken by the supersolid order, it is then natural to know whether the filled spinon bands support a nontrivial Chern number.When one goes beyond the mean-field theory, the gapped spinon bands with a nontrivial Chern number could induce a Chern-Simon term in the U(1) gauge theory that is very much like the one in chiral spin liquids.
Since all these bands are now overlapping in energy, one needs to invoke the non-Abelian Berry connection to define the Chern number of the filled bands.The Chern number of the total 6 filled spinon bands is welldefined by an integral over the continuum magnetic Brillouin zone (MBZ) as where the non-Abelian Berry connection [53] A is a 6 × 6 matrix-valued one-form in the momentum kspace associated with the wavefunctions for the filled spinon bands ψ − k = (|ψ − 1k ⟩, ..., |ψ − 6k ⟩).Numerically we evaluate this Chern number using the Fukui method [54] in the M ×M -grid discrete momentum space k where the lattice field strength 3 ), and the link variable ) With the periodic boundary condition, the MBZ is a closed two-dimensional torus T 2 , guaranteeing the Chern number an integer.Then, with sufficiently fine grids, the Chern number on the discretized space C− equals to C − , and we find C− = 0 up to a grid choice M = 300.Thus, we conclude that the gap opening from the precursory DSL with the supersolid order is not topological and there is no Chern-Simons term.

Appendix C: Effect of external magnetic field
Here we further examine our model with the presence of a weak external magnetic field.For the in-plane magnetic field, the continuous U(1) symmetry will be absent, and the Goldstone mode will disappear.The magnetic excitation will be fully gapped.For out-of-plane field with B = B ẑ along the z direction, the continuous U(1) symmetry is still present in the model.Although the field slightly polarizes the spins, as long as the field is perturbative on the supersolid order, the Goldstone mode of the spin-wave excitation should persist in the weak field regime.For the spinon continuum, we do not expect a significant change as well, and the continuous excitations should persist in the weak field regime.In Fig. 9, we computed the INS results for the spinon continuum under the weak field.We have let the magnetic moments on the red and green sublattices to rotate while keeping their magnitude fixed.The net magnetic susceptibility is dM/dB = 1µ B T −1 /Co [20].We find that, the addition of B only slightly affects the spinon spectrum.One can understand this by reviewing the statements that were made in Fig. 4, where the addition of Jz splits up the spin-up and spin-down bands.Here the external field in the z direction does a similar thing, that is to further split them, while the gap size is not much affected since it is determined by J⊥ .As for the INS spectrum, the effect of external B z is almost negligible in our calculation.The predicted spectral continuum at the K and M points still exists.

Appendix D: Finite temperature behaviors
Since the spinon gap is induced by the supersolid magnetic order, when the thermal fluctuations suppress the magnetic order at low and finite temperature, the system would cross over to the Dirac spin liquid regime.Although this is the finite temperature paramagnetic phase, the physical properties are governed by the Dirac spin liquid.The expected finite-temperature phase diagram is depicted in Fig. 10, where the increased quantum fluctuation would further stabilize the spin liquid state at zero temperature.The similar picture was actually invoked for α-RuCl 3 [55] that also orders at low temperature.Thus, at the temperature above the ordering, the (gapless) Dirac spinons should give a heat capacity of T 2 and contribute to the thermal transport.As we have mentioned that, Ref. 21 already found the signature of itinerant excitations.Moreover, Ref. 56 further proprosed the gapless spin liquid from the finite temperature thermal Hall transport result.Although the thermal Hall result is not inconsistent with a gapless spin liquid, the generation of Berry curvature distribution for the spinons requires further understanding.while to raise some further discussion about the proximate states out of the quantum spin supersolid and some of the related properties.This may provide some further insights for the interested readers.The material Na 2 BaCo(PO 4 ) 2 is located in the fully frustrated regime of the XXZ model on the triangular lattice and thus has a sign problem for quantum Monte Carlo simulation.Most of the early analysis of the triangular lattice XXZ model FIG.10.The expected finite temperature behaviors, and the dashed line is for Na2BaCo(PO4)2.The system has a weak supersolid order at very low temperature.Thermal fluctuation suppresses the order and drives the system into the Dirac spin liquid regime.
was devoted to the less frustrated regime of the XXZ model, i.e.
J ⊥ < 0 where the quantum supersolid is well present before the knowledge on the fully frustrated side due to the absence of the sign problem [29,30].The intimate relation between the fully frustrated regime with J ⊥ > 0 and the less frustrated regime with J ⊥ < 0 was partially addressed in Ref. 23.It was shown that [23], in the perturbative regime J ⊥ ≪ J z , one can perform a unitary transformation to connect both regimes.This establishes the presence of the supersolid in the fully frustrated regime and argued to extend to the Heisenberg point based on the continuity argument.This unitary transformation that connects the fully frustrated and less frustrated regime is analogous to the one that was used in the pyrochlore spin ice context [57].
On the less frustrated side with J ⊥ < 0, due to the underlying U(1) symmetry of the XXZ model, a bosonvortex duality was adopted to explore the proximate states and phase transitions by viewing the spins as hardcore bosons [58].It was predicted from an emergent SU(2) symmetry that, the quantum supersolid is proximate to deconfined quantum criticality of the NCCP 1 universality class.Remarkably, this is identical to the Néel-VBS transition.From the recent duality argument, this criticality is dual to the N f = 2 fermionic quantum electrodynamics with Dirac fermions [59].Although this is on the less frustrated side and may not be directly re-lated to Na 2 BaCo(PO 4 ) 2 , it provides some insights about the connection between the quantum supersolid and the Dirac fermions, and may be useful for the full understanding of the nearby phases of the quantum supersolids in both regimes.

FIG. 1 .
FIG. 1.(a) The spin orientation of various magnetic phases for the easy-axis XXZ model on the triangular lattice.As the field B = B ẑ increases from zero, the system enters the 3sublattice supersolid (Y-shape), UUD, V-shape and the fully polarized states.The sublattices are plotted in red, green and blue.(b) The spinon hoppings tij are chosen with a gauge choice supporting a π flux.

FIG. 2 .
FIG. 2. (a) The spinon bands in MBZ with Jz = 4.394t, J⊥ = 4.102t, µ = 0.The weakly gapped Dirac cones are relevant for the low-energy INS signals.The bands below (above) the Fermi surface are in blue (red).The orange rectangle indicates the MBZ.(b) The spinon dispersion along the red dashed line in (a).(c) The hexagonal lattice Brillouin zone with high symmetry points.The orange rectangle is the MBZ.

FIG. 3 .
FIG.3.The band dispersions for Jz = 0 and J⊥ = 0.The energy is in the unit of t.The up-spin spinon bands in red is degenerate with the down-spin spinon bands in blue, and the up-spin and down-spin bands are degenerate so the color cannot distinguish them.The Dirac cones are explicitly shown.

FIG. 4 .
FIG. 4. The upper figure is the spinon band dispersions for Jz|⟨S z 0 ⟩| = 0.5t and J⊥ ⟨S x 2 ⟩ = 0. Energy is in the unit of t.The red (blue) Dirac cones formed by up-spin (down-spin) bands do not open the gap.The lower figure is the spinon dispersion near the Fermi level where the gapless feature persists.

FIG. 6 .
FIG. 6.(a) The INS spectrum along the high-symmetry path with the same parameters in Fig. 2. In the low-energy regime of the INS spectrum, there are two strong signals at Γ and M, and lighter ones at K and K/2.(b) The constant energy slice of (a) at E/t = 1 and is marked by the red arrow in (a).(c) The momentum dependence of S −+ (p, E) at E/t = 1.

ACKNOWLEDGMENTS
We particularly thank Liusuo Wu for sharing his data prior to the submission, and Cenke Xu and Jiawei Mei for discussion.This work is supported by NSFC with Grant No. 92065203, MOST of China with Grants No. 2021YFA1400300, and by the Collaborative Research Fund of Hong Kong with Grant No. C6009-20G and C7012-21G, by the Guangdong-Hong Kong Joint Laboratory of Quantum Matter, the NSFC/RGC JRS grant with Grant No. N HKU774/21, by the General Research Fund of Hong Kong with Grants No. 17310622 and No. 17303023, and by the Fundamental Research Funds for the Central Universities, Peking University.hopping [See Fig. 8].

FIG. 8 .
FIG. 8.The upper figure is the triangular lattice formed by the Co 2+ ions with the basis vectors a1 = a(1, 0) and a2 = a(1/2, √ 3/2).The hopping parameters tij of the spinons are indicated within the bond with a gauge choice supporting a π-flux background.In the lower figure, the black dashed lines show the lattice Brillouin zone and the orange solid lines show the magnetic Brillouin zone.The corresponding high symmetry points are also denoted.