Introduction.-A hallmark of strongly correlated systems is the emergence of novel degrees of freedom at low energies from strong correlations. A prototype case is fractionalized excitations-fundamentally different from excitations in weakly interacting limit-such as spinons in herbertsmithite ZnCu 3 ðOHÞ 6 Cl 2 [1] and YbMgGaO 4 [2][3][4]. A particularly intriguing possibility arises in quantum spin liquids because their emergent fermionic excitations can form a Fermi surface in momentum space, rendering the properties of these insulators akin to those of conventional metals. The two-dimensional (2D) triangular lattice-based organic compounds EtMe 3 Sb½PdðdmitÞ 2 2 and κ-ðETÞ 2 Cu 2 ðCNÞ 3 [5][6][7][8] are among the most famous candidate materials believed to host such a critical spin liquid (CSL) with an emergent spinon Fermi surface (SFS) [9,10]. A four-spin ring exchange is needed to describe these materials [11][12][13]. An outstanding challenge is how to demonstrate and reveal the presence of fractionalized fermionic excitations, particularly with regards to the SFS.
On the theory side, it has been proposed to study the emergent SFSs in CSLs through the singular peaks in the spin structure factor (SSF)-those that arise from realspace power-law decaying spin correlations-which can be related to the locations of the SFS [9]. Using this procedure, recent density matrix renormalization group (DMRG) results reported the possible SFS of the spin-1=2 model on a triangular lattice with a four-spin ring exchange [12,13] and in the Kitaev model on a honeycomb lattice [14]. However, it is still difficult to reconstruct the actual shape of the SFS through the DMRG results of the SSFs based on small system sizes.
An alternative quantity to describe long-range entangled states is the entanglement entropy (EE) [15], such as the von Neumann EE and the Renyi EE (REE), which are obtained from reduced density matrix of a subsystem by tracing out the degrees of freedom outside this subsystem. The EE plays an important role in several fields, ranging from quantum information to condensed matter physics [16], and has been measured experimentally [17]. It is believed that the EE of the ground states in most local Hamiltonians satisfies the "EE area law" [18]: when a system is divided into subsystems, the EE is proportional to the area of the boundary between the two subsystems at the leading order.
Violations of the EE area law do exist in various cases. In one dimension, they are found in several quantum critical systems [19]. In higher dimensions, these violations are associated with the presence of a SFS in momentum space. The most well-known examples are the ground states of free fermions with Fermi surfaces [20,21], where the violation is logarithmic, i.e., the EE is proportional to the surface area multiplied by a factor that grows logarithmically with the subsystem size. Intriguingly, the EE in these noninteracting systems takes the Widom formula [20][21][22], where the coefficient of the leading term in the dependence of EE on the subsystem size captures the geometric information of the Fermi surface and that of the subsystem. For gapless electronic systems, calculations perturbative in the interactions [23] show that such a violation retains the same form as that of a free Fermi gas. Recently, it has been suggested that the EE associated with the composite Fermi liquid phase of the half-filled Landau level (ν ¼ 1=2) is also described by the Widom formula [24]. By contrast, for frustrated strongly correlated electrons, as in Hubbard models, or spin systems, as in Heisenberg models, all with a possible emergent SFS, the EE has not been much explored [25,26].
Our present work goes beyond previous efforts and provides a generic procedure to reconstruct the geometry of the emergent SFS. We present the first variational Monte Carlo (VMC) study of the EE to test the conjectured Widom formula for strongly correlated systems. Employing a widely discussed example of a CSL with an emergent SFS, we introduce a direct probe of emergent fractionalized excitations using the real-space EE together with examining the singularity of the SSF. Remarkably, we show that the leading order of the EE has the form of the Widom formula multiplied by a previously unknown factor of 2. This numerical factor captures the presence of two free gapless modes associated with two flavors of spinons. From this formula, we provide the basis [27] for a systematic methodology to explicitly reconstruct the emergent SFS geometry. We remark that using the SSF or EE individually only allows you to test the existence or not of fractionalized excitations (i.e., a "yes" or "no" answer), but a combined methodology is necessary to recover the full shape of the SFS. We also remark that we employ VMC only for simplicity: our methodology can be used if other techniques are employed, such as the quantum Monte Carlo technique or DMRG. With the only caveat that it is advisable to employ several trial states to search for self-consistency to remove the bias uncertainty intrinsic of variational procedures, our procedure is quite generic.
Entanglement entropy and Widom formula.-We will first provide robust numerical evidence for the validity of the Widom formula in a CSL with emergent SFS. A typical ground-state wave function (WF) to represent the possible CSL on a triangular lattice [5][6][7][8] is the Gutzwiller projected Slater determinant:jψi ¼ P G jψ 0 i, where the Gutzwiller
forbids double occupation on each site, and jψ 0 i is the ground state of the meanfield Hamiltonian on the triangular lattice
The Gutzwiller projector P G is crucial to avoid a trivial Fermi surface of real electrons, while still allowing a possible SFS. This variational WF is known to be accurate for the quasi-1D J 1 -J 2 spin-1=2 chain with four-spin exchanges [28], providing a reasonable starting point for our effort. We begin by considering an isotropic system with a total number of sites N s ¼ L × L [t 0 ¼ t in Fig. 1(a)].
Based on the Widom formula, the REE associated with a subsystem consisting of L A × αL A sites along the a 1 ≡ ð1; 0Þ and a 2 ≡ ð1=2; ffiffi ffi 3 p =2Þ directions (lattice constant a ≡ 1), as illustrated in the bottom left portion of Fig. 1(a), can be concisely expressed as follows (derivations in Supplemental Material [29])
where _ ¼ means the leading logarithmic contribution in REE, α represents the ratio between the linear length of the subsystem (L A ) and that of the whole system (L), i.e., α ¼ L A =L, and c eff is effectively the number of free gapless modes in the low-energy limit. Additionally, A sf refers to the cross section of the SFS, which is determined by the span in the momenta between right or left moving patches (k n FR=L ) of the SFS along any particular observation 1), we plot S 2 =½ð1 þ αÞL A vs lnðL A Þ. The slopes of the lines give the prefactor of the leading term of the REE. We fix the whole system size and choose the subsystem size to be L A × αL A with α ¼ 1 and 1=2. The data on L ¼ 18 with α ¼ 1 (red squares) are consistent with Ref. [26]. We can clearly observe the proportionality
direction n. This is illustrated in the top right portion of Fig. 1(a), where the emergent SFS is expected to be circular.
We have carried out the VMC simulations on the triangular lattice with the whole system size fixed to be L × L, with L up to 20. We calculated the REE associated with a subsystem of L A × αL A sites, where both L A and αL A are less than or equal to L=2. The resulting REE vs L A is plotted in Fig. 1(b), which shows that S 2 =½ð1 þ αÞL A vs ln L A has the same slope for different choices of α within error bars. The proportionality S 2 ∼ ð1 þ αÞL A ln L A provides direct evidence that the REE of the CSL studied here satisfies the Widom formula Eq. ( 1). The slope in Fig. 1
In order to pin down the explicit formula for the REE of a CSL, additional information is needed to determine the values of c eff and A n sf separately, as addressed next.
Spin structure factor.-Using the VMC described earlier, we calculated the SSF D q ≡ ð1=N s Þ× P i;j hS i • S j ie iq•ðr i -r j Þ with the spin operator S i ¼ P σ;σ 0 1 2 c † iσ σ σ;σ 0 c iσ 0 . It is known that for an arbitrary observation direction n, D q should show singular peaks at q ¼ 0 and k n FR -k n FL , which are associated with forward and backward scattering processes. The information of D q can be used to determine the cross section of the emergent SFS whose surface unit vector is perpendicular to n, i.e., A n sf ¼ jk n FR -k n FL j [28]. In the isotropic case, A n sf ¼ A sf is independent of the direction.
In Fig. 2 we show the numerical data for the SSF on a triangular lattice with 30 × 30 sites. Figure 2(a) gives a 3D side view of the SSF in the Brillouin zone (BZ), denoted by the black hexagon, where we can see a sharp singular point at q ¼ 0 and weaker singular lines on the surface whose locations are theoretically suggested to be q ¼ k n FR -k n FL . Figure 2(b) shows the 3D top view of D q . In the present finite-size calculations, the singular lines on the 3D D q surface are more clearly revealed near the BZ boundary, while the weaker singular lines inside the BZ are masked by the sharper singular point at q ¼ 0. From Fig. 2(b), we can determine the location of the full singular lines by fitting k n FR -k n FL [27], which allows us to extract the (average) cross sections of the emergent SFS to be 5.24 AE 0.05. When this value for the cross section is combined with the slopes of the normalized REE vs lnðL A Þ shown in Fig. 1(b), we obtain c eff ≃ 2.01 AE 0.02. This value indicates the presence of two free gapless modes for each "independent" 1D patch in the low-energy limit [30], so it should be universal for all shapes of convex critical Fermi surfaces. If we introduce anisotropy into the system, c eff should remain the same.
Visualizing emergent spinon Fermi surface.-The explicit formula for the EE obtained above can be used to reveal the emergent SFS directly. For an isotropic system, since the shape of an emergent SFS is circular, and its diameter can be extracted once the REE is calculated. To address a more general case, we focus on a triangular lattice system with anisotropy. Specifically, we consider a Gutzwiller-projected WF with hopping amplitudes t along each ladder (AEa 1 directions) and t 0 along the zigzag directions [AEa 2 and AEða 1 -a 2 Þ] that couple different ladders as shown at the bottom right of Fig. 1(a). For an illustration, we use t 0 =t ¼ 0.7 to obtain the REE associated with the subsystems. FIG. 2. The 3D plot of SSF for the isotropic case. (a) Side view of the SSF within the BZ for the isotropic case (black hexagon represents the BZ), employing a triangular 30 × 30 cluster. There is a sharp and singular peak at q ¼ 0, which corresponds to the uniform real-space power-law decaying behavior. The much weaker singular lines near the boundary of the BZ correspond to the oscillating real-space behavior caused by the presence of the SFS. (b) Top view of the SSF. The blue lines are fitting results that are seen to match well with the singular lines for the SFS in the SSF on a triangular lattice consisting of 30 × 30 sites. The details of the fitting method are in Supplemental Material [29].
Because of numerical and computational time limitations, below we choose three subsystem geometries to obtain the REE and thereby construct the anisotropic SFS. Specifically, we calculate the REE for a subsystem with L A × αL A sites, where we consider ratios α ¼ 1=2; 1; 2. The REE results in these systems are shown in Fig. 3(a). Setting c eff ¼ 2 for the present anisotropic system the formula for REE becomes
where A a 1 or 2 represents the cross sections of the SFS projected onto the a 1 or 2 axis. We can write down three equations, corresponding to α ¼ 1=2, 1, and 2, respectively: (i)
We can choose any two out of the three equations to obtain the values of A a 1 or 2 . Since there are three choices, we can obtain three numerical approximations for A a 1 or 2 , which we average over to reduce the statistical error. We find that A a 1 ≃ 1.53π and A a 2 ≃ 1.89π, based on which the shape of the emergent SFS can be constructed as illustrated in Fig. 3(b). The green and the blue lines represent the cross sections of A a 1 and A a 2 . Since there is an inversion center for the emergent surface in momentum space [27], once A a 2 is known, we can draw its inverted partner, denoted as Ãa 2 (brown line) in Fig. 3(b). The dashed lines are perpendicular to A a 1=2 and Ãa 2 , respectively. Connecting all the intersections of the dashed lines results in the red hexagonal shape, which provides the leading-order approximation to the shape of the emergent SFS. In principle, we can improve the accuracy of the shape if we perform more (time consuming) REE calculations using different subsystem geometries [27].
For comparison, we also show the shape of the SFS in Fig. 3(b) (light-gray ellipse) obtained by extracting k n FRk n FL from the SSF. The exact numerical results for the SSF are shown in Supplemental Material [29]. The emergent SFS reconstructed from the REE results is quite consistent with the light-gray ellipse in Fig. 3(b), which provides additional support for our procedure. With (costly) additional values of α our results will be even closer to the ellipse. We remark that in strongly correlated systems, where analytical methods are difficult to use and numerical simulations only can be performed on small clusters, it may be difficult (or sometimes impossible) to determine the locations of k n FR -k n FL and thus the here proposed EE probe becomes the only practical procedure, exhibiting its unique value. From this overarching perspective, the present work builds up a foundation for using the EE to probe emergent SFSs in general cases.
Conclusion and outlook.-In this work, we examined the entanglement properties of a CSL with an emergent SFS. Numerically, we have proved the validity of a generalized Widom formula Eq. (1) [29] for this type of strongly correlated systems. Based on this formula, we provide a general procedure to reveal and construct the shape and size of emergent SFSs, by examining the singularity of the SSF and the real-space EE. This is an advance over previous efforts that relied on the singular peaks in the SSFs to locate the SFS by DMRG, because using only the latter the whole shape of the SFS cannot be obtained. In addition, we have obtained the universal factor c eff ¼ 2 that describes two free gapless modes in a CSL employing robust numeral calculations, without "guessing" this value in advance.
The current work can be straightforwardly generalized to CSLs of higher-spin (S ≥ 1) systems. Of particular interest is the 6H-B phase of S ¼ 1 Ba 3 NiSb 2 O 9 [31,32] that was recently suggested to realize a CSL with three flavors of fermionic spinons, forming a large SFS [33]. From our perspective, it is always possible to write down a Gutzwillerprojected WF of three flavors of fermions to represent the S ¼ 1 CSL. Based on the results presented here, we conjecture that the leading EE in this case also satisfies the Widom formula, but with c eff ¼ 3. Finally, our work points to new prospects for deepening the understanding of correlated systems such as heavy-fermion materials, in which the nature of quantum spins and Fermi surface plays a crucial role [34]. Examining the quantum entanglement properties promises a conceptually new way of elucidating their quantum phases and criticality.