To analyze the QSL phase of this model, we here adopt the well-known and widely-used Abrikosov fermion construction since it can be utilized to study both gapped and gapless phases, while the Schwinger boson formalism has the limitation to study gapped phases.^[27] In the Abrikosov fermion representation, the effective spin-1/2 operator S_i on site i is represented by where ${f}_{i,\alpha }^{\dagger }\,({f}_{i,\alpha })$ creates (annihilates) a spinon with the spin index α = ↑,↓ at site i, and σ is a vector of three Pauli matrices. The Hilbert space constraint $\sum _{\alpha }{f}_{i\alpha }^{\dagger }{f}_{i\alpha }=1$ on local fermion number is imposed to project out the unphysical states and faithfully reproduce the physical Hilbert space. Substituting Eq. (2) into the J₁–J₂ Hamiltonian Eq. (1), one obtains an interacting four-fermion system. Further performing a mean-field decoupling would reformulate the interacting fermionic system to a quadratic level.^[15] Specifically, by ignoring the pairing channel, the general quadratic spinon Hamiltonian with only spinon hopping sector is obtained as where we have maintained the SU(2) spin rotation symmetry of the original spin model and the local occupation constraint is relaxed such that only its average value satisfies $ \sum _{\alpha }\langle {f}_{i\alpha }^{\dagger }{f}_{i\alpha }\rangle =1$. The global chemical potential μ is introduced as a Lagrange multiplier to enforce such a constraint. Moreover, the mean-field parameters t_1,ij and t_2,ij represent the hopping amplitudes between NN and NNN sites, respectively. In the numerical studies, such as variational Monte Carlo approach, a similar mean-field Hamiltonian to Eq. (3) can also be exploited as a good starting point to construct the many-body variational wavefunction.

Here we stress again our purpose in the following is not to solve for the detailed ground state of a spin Hamiltonian as in Eq. (1). Instead, we assume that the system stabilizes a magnetically disordered QSL phase in the intermediate region 0.4 ≲ J₂ / J₁ ≲ 0.6, as suggested by a variety of numerical studies.^[22–26] Comparing with pursuing solving a spin model exactly, it might be a better, or at least as a supplementary strategy to start from the potential QSL states and then single out the nontrivial and robust experimental signatures that allow us to distinguish a QSL. We then start from a mean-field theory to proceed with our analysis since much can already be learned from a mean-field investigation. A full treatment of the original spin model requires the involvement of all quantum fluctuations of the parameters around the mean-field solution, but the robust and intrinsic experimental signatures could maintain even the fluctuations are included.

The spinons fulfill the projective symmetries of the square lattice, and the mean-field parameters t_1,ij and t_2,ij, also called mean-field ansatzs in the literatures, should be constrained by a systematic projective symmetry group (PSG) analysis,^[15] which results in a classification of all possible QSLs. It is shown^[15] that under open boundary condition of the square lattice, one can always fix the gauge such that the gauge transformations for translations satisfy where T₁, T₂ are two translations of the square lattice, and ${{\mathscr{G}}}_{T}(i)$ is the associated SU(2) gauge transformation. The spinons thus fulfill the gauge enriched^[15] translations ${\mathop{T}\limits^{\sim }}_{1}={{\mathscr{G}}}_{{T}_{1}}{T}_{1}$ and ${\mathop{T}\limits^{\sim }}_{2}={{\mathscr{G}}}_{{T}_{2}}{T}_{2}$. The parameter η_xy takes values ±1, then it generally divides all the possible QSLs on the lattice into two categories, i.e., zero-flux states correspond to η_xy = +1 and π-flux states correspond to η_xy = –1. In particular, the π-flux states harbor a fractionalized translation symmetry, which would result in a sharp signature in experiment. To illustrate this idea, we adopt the simplest two cases and the values of t_ij in Eq. (3) we have taken are shown in Fig. 1, where figure 1(a) corresponds to a zero-flux QSL, while figure 1(b) corresponds to a π-flux QSL. It can be easily verified that t_1,ij = t_1,ji = –t₁ on the (thick) red bonds in Fig. 1(b) involve a background π gauge flux penetrating per square plaquette, while this flux has no influence on the NNN hoppings. In fact, figures 1(a) and 1(b) correspond to a spinon Fermi surface QSL and a π-flux Dirac QSL, respectively.

Inelastic neutron scattering (INS) measurement represents the best experimental probe to directly detect the magnetic excitations, and the dynamical information of excitations is encoded into the dynamical spin structure factor where N is the total number of lattice sites and the summation runs over all the excited eigenstates |n〉 with ξ_n(q) being the energy of the n-th excited state with momentum q, while |G〉 stands for the spinon ground state with spinons filling the Fermi sea. In the numerical calculations, the delta function is taken to have a Lorentz broadening, δ(ω) = η/[π(ω² + η²)] with η = 0.1t₁. Moreover, since ${\boldsymbol S}_{q}^{+}= \sum _{k}{f}_{k+q,\uparrow }^{\dagger }{f}_{k,\downarrow }$, the summation in Eq. (5) should be over all possible spin-1 excited states that are characterized by one spinon particle–hole pair crossing the spinon Fermi energy with a total energy ω and total momentum q.^[9] In other words, the momentum-transfer q and energy-transfer ω(q) of the neutron should be shared between the spinon particle–hole pair where ω₁(k) [ω₂(k)] is the spinon (hole) excitation energy with momentum k and the minus sign comes from the nature of hole excitation. The above equations indicate the two-spinon spectrum continuum, which is often count as a manifestation of fractionalized excitations, is a general character of QSLs. In Figs. 2(a) and 2(b), we present the density plots of dynamical structure factor ${\mathscr{S}}(\boldsymbol q,\omega )$ along the high symmetry line Γ–M–Γ–X–M marked in Fig. 2(e) for zero-flux QSL and π-flux QSL, respectively. As depicted in Fig. 2(a), besides the obvious spectral continuum, a clear V-shape appears around the Γ point, which originates from the particle–hole pair excitations crossing near the Fermi surface.^[9] While for the π-flux Dirac QSL in Fig. 2(b), the obvious phenomenon becomes the low-energy cone features that originate from the inter- (large q) and intra-Dirac (small q) cone scatterings, which is just the reflection of a Dirac QSL with Dirac band touching, not the π-flux for the spinons of the QSL. However, these signatures and evidences seem not to be strong enough to confirm a QSL, as it has been explicitly shown that a simple spectrum continuum, including the V-shape feature, in the dynamical spin structure factor can also be explained in the scenarios of usual glassy and disorder-induced states.^[12–14] Therefore, additional signature of spectrum that is more unique to QSLs is expected to diagnose the fractionalized excitations.

For this purpose, next we consider the intrinsic fractionalized translation symmetry due to the π gauge flux of the π-flux QSL. According to Eq. (4), the translation symmetries of spinons should satisfy where ${\mathop{T}\limits^{\sim }}_{1}$ and ${\mathop{T}\limits^{\sim }}_{2}$ are the gauge enriched translations acting on the spinon degrees of freedom instead on the physical spins. For the zero-flux state with η_xy = 1, the translations commute as usual, while for the π-flux state with η_xy = –1, the anti-commutation of translations implies a fractionalization of symmetry. It was first realized that the crystal momentum fractionalization of internal spinon degrees of freedom has dramatic effects on the neutron spectrum.^[16,17] As a consequence, the periodicity of the upper excitation edge of the dynamic spin structure factor defined by is doubled.^[28] Therefore, for the zero-flux state, the upper two-spinon excitation edge should have the usual periodicity, while for the π-flux state, the upper two-spinon excitation edge exhibits an enhanced periodicity, that could serve as a sharp identification of fractionalized excitations beyond the simple spectrum continuum, since it is impossibly mimicked by any disorder-induced states. We illustrate the contour plots of the upper edge of the dynamic spin structure factor for the zero-flux QSL and π-flux QSL in Figs. 2(c) and 2(d), respectively. It is clear that the π-flux QSL exhibits a fractionalization pattern with an enhanced periodicity in Fig. 2(d), which is readily accessible to the INS measurements. The enhanced spectral periodicity with a folded Brillouin zone is the dynamical property rather than the static property, and can not be captured by the static spin structure factor.