Electromagnetic signatures of chiral quantum spin liquid

Quantum spin liquid (QSL) has become an exciting topic in interacting spin systems that do not order magnetically down to the lowest experimentally accessible temperature; however, conclusive experimental evidence remains lacking. Motivated by the recent surge of theoretical and experimental interest in a half-filled Hubbard model on the triangular lattice, where chiral QSL can be stabilized, we investigate the electromagnetic signature of the chiral QSL to aid experimental detection. We systematically studied the electrical charge and orbital electrical current associated with a spinon excitation in the chiral QSL based on parton mean-field theory and unbiased density-matrix renormalization group calculations. We then calculated both longitudinal and transverse optical conductivities below the Mott gap. We also conduct quantum field theory analysis to unravel the connection between spinon excitation and emergent and physical gauge fields. Our results show that the chiral QSL phase has a clear electromagnetic response even in a Mott insulator regime, which can facilitate the experimental detection of this long-sought-after phase.


I. INTRODUCTION
Quantum spin liquid (QSL) states are interacting quantum spin systems that do not order magnetically down to zero temperature.This absence of magnetic order leads to a quantum-disordered ground state with characteristic long-range quantum entanglement, fractionalized excitations, and its associated emergent gauge fields.Consequently, it has been challenging to understand and characterize QSL since its inception [1].The experimental detection of QSL states becomes even more difficult due to the lack of a conventional order parameter [2].However, recent developments in both the theoretical and experimental fronts have led to a continuous surge of interest in analyzing and detecting this illusive state of matter [3].Examples range from the discovery of various iridates/ruthenates compounds as candidate materials to realize proximate Kitaev physics [4,5] to the observation of topological spin liquids in the Rydberg atom quantum simulator [6] and quantum processor [7].The appearance of QSL requires suppressing the magnetic orders, and therefore frustrated magnets are the playground for hunting for QSL.In this regard, the triangular lattice Hubbard model (TLHM) has always remained a centerpiece of attention.
In the large U limit of the TLHM at half-filling, the effective low energy Hamiltonian is an antiferromagnetic Heisenberg model which stabilizes the conventional 120 ○ (Néel) order [8][9][10].However, it is widely believed that the ground state of the TLHM drifts toward a QSL state when the correlations become weaker but remain above the Mott transition [11].Recently, various densitymatrix renormalization group studies (DMRG) [12][13][14][15], and matrix product state (MPS) [16] analyses on TLHMs have predicted the evidence for a Kalmeyer-Laughlin type chiral quantum spin liquid (cQSL) phases [17], see Fig. 1 for a schematic phase diagram.
The TLHM can be realized in certain materials.Previous experimental work has shown characteristic evidence for a QSL phase in certain organic Mott insulators [18][19][20][21].Although, the controversy over the gapped [20] or gapless [22] nature of the underlying excitations still remains.In another triangular lattice material, YbMgGaO 4 , the gapless character is well supported by the nuclear magnetic resonance [23] and muon spin rotation [24] experiments, as well as evidence of a spinon Fermi surface revealed by neutron scattering studies [25].Therefore, it is necessary to look for some smoking-gun signatures that can decipher the true nature of the QSL phase.
Motivated by the identification of the cQSL phase in TLHM and its potential relevance in several compounds, here, we systematically analyze its electromagnetic responses.Despite being a Mott insulator, there is a remnant electromagnetic response due to the virtual hopping of electrons [26,27].Assuming a cQSL phase, which spontaneously breaks time-reversal symmetry (TRS), we analyze the corresponding effective spin model [13] within FIG. 1.A schematic phase diagram for the triangular lattice Hubbard model at half-filling with a metallic phase at small U , followed by a putative cQSL phase with non-vanishing chiral order parameter χ = ⟨Si ⋅ (Sj × S k ⟩ at an intermediate coupling regime U1 ∼ 9t, and U2 ∼ 11t [12,15], and a magnetic ordered Néel state at strong coupling.Note that χ = 0 in the other two phases.arXiv:2304.08635v2[cond-mat.str-el]24 May 2023 the parton mean-field spinon description and obtain the associated orbital magnetization and electrical polarization.We also performed unbiased DMRG calculations on the half-filled TLHM at an intermediate coupling U (U 1 < U < U 2 ).Our numerical analysis further supports the mean-field results for the electromagnetic responses.To have a universal picture, we additionally employ the quantum field theory description to elucidate explicitly the relationship among the emergent and the physical gauge fields and low energy spinon excitations in the cQSL.
To relate our theoretical framework to experiments, we compute the transverse optical conductivity (within the spinon description), which is associated with the magneto-optical Faraday rotation (MOFE) where l is the thickness in the direction of light propagation with frequency Ω, n is the index of refraction, and σ ′ xy (Ω) is the real part of the optical conductivity in 3D.Our electromagnetic response functions, including the orbital magnetization profile and the structure of Θ F , provide a clear experimental signature of the cQSL.
For completeness, we also analyze the behavior of the dynamic spin-structure factor and lay out the possible connection with the relevant experiments.
The rest of this paper is organized as follows: in Sec.II, we provide the spinon description of the cQSL in TLHM.In Sec.III and Sec.III A, we provide details of the derivation for electrical polarization and orbital magnetization.The DMRG calculations supporting our mean-field calculations are given in Sec.III B. Sec.III C provides a picture based on quantum field theory.In Sec.III D, and Sec.III E, we compute the dynamic spin structure factor and transverse optical conductivity with the electrical polarization, respectively.Finally, we discuss the implications of our results and conclude in Sec.IV.

II. MODEL
We start with the TLHM at half-filling with the corresponding Hamiltonian written as where c † iσ creates an electron at site i with spin σ, and U is the strength of the onsite Coulomb repulsion.In the strong coupling limit (U ≫ t), the charge degrees of freedom are gapped out, and the relevant microscopic model can be analyzed in terms of an effective spin model.Within a second-order perturbation expansion in t/U , the corresponding spin Hamiltonian reads ∑ ⟨ij⟩ S i ⋅ S j , where J (2) = 4t 2 /U is the antiferromagnetic Heisenberg coupling.However, in the intermediate coupling regime, i.e., U ≳ t, the above second-order perturbation does not completely capture the lowenergy dynamics, and we need to include higher-order spin corrections.Such a procedure leads to further neighbor spin exchange terms, including ring exchange-like interactions [28].Therefore, although a Néel order is preferred at larger U , incorporating subleading order correction modifies the overall magnetic order at an intermediate U .Previous theoretical works [26,[29][30][31] have reported the existence of two critical coupling strengths U 1 ∼ 9t, and U 2 ∼ 11t.The current consensus is that TLHM hosts a putative QSL phase in the intermediate regime between U 1 and U 2 , eventually becoming a Néel ordered state at a larger U > U 2 .
Motivated by these previous studies and recent developments in the DMRG results [12,15], we adopt a phenomenological chiral spin liquid model to describe its concomitant features.The effective Hamiltonian, which hosts cQSL as a ground state, is written as where the associated exchange couplings are written in terms of the parameters of the original low-energy spin model.Here ⟨ij⟩ denotes the nearest sites and ⟪ijk⟫ denotes three sites in a unit triangle.It was argued in Refs.[13,26] that the four-spin ring exchange term (see SM [32] for details) is responsible for the appearance of the chiral term in Eq. ( 3).
Here, we focus on the model as in Eq. ( 3) and analyze it within a mean-field description.We utilize the standard parton decomposition of the spins as S i = 1 2 f † iα σ αβ f iβ , where f † iα creates a neutral spinon excitation with spin α at site i, and σ denotes the vector of Pauli matrices (the repeated indices are assumed to be summed over).This fractionalization leads to an enlargement of the Hilbert space.Therefore, one needs to implement a local constraint (f † iα f iα = 1) to project to the physical Hilbert space.Plugging this back into Eq.( 3) and assuming a nonzero mean-field decomposition as m ij = ⟨f † iα f jα ⟩, we obtain a noninteracting spinon Hamiltonian as (see supplementary material (SM) [32] for details) where the primed summation corresponds to all the permutations between the three neighboring sites i, j, k.
Here, we adopted a mean-field decomposition only in the particle-hole channel, although a more general decomposition with both particle-particle and particle-hole channel may provide a qualitatively better description of the emergent spinon spectrum [33,34].
Assuming the translational invariance, we simplify the mean-field order parameter m ij = m 0 e iϕij , where m 0 is the amplitude, and ϕ ij 's are bond-dependent phases.Subsequently, we capture the physics of the Hamiltonian in Eq. ( 4) with a simplified model as Focusing on a three-site cluster, the hopping amplitude t, and the phases ψ ij 's are related to the parameters in Eq. (4) as However, the phases ψ ij 's and the hopping t remain undetermined.To further progress, we utilize Lieb's theorem [35], which states that a fermion hopping on a bipartite lattice realizes its ground state with π-flux per bipartite plaquettes.Since the triangular lattice is monopartite, we consider a decorated lattice comprised of doubled unit cells [see Fig. 2(a)] with the hopping amplitudes between different neighboring sites such that the total flux within the rhombus-shaped unit cell is π.In such a construction, we can do further simplification and solve Eq. (6a), and Eq.(6b) to show that [32] with the constraint, the total flux within a triangle is π/2.Note that m 0 still remains undetermined.A particular choice of ψ ij is shown in Fig. 2(a) to realize the staggered flux configurations between the up and the down triangles, where θ = 0 corresponds to π/2 flux within a triangle.TRS is preserved for θ = π/2.Diagonalizing the Hamiltonian in Eq. ( 4) obtains the corresponding spinon band structure.The uniform flux phase (θ = 0) leads to a gapped spinon spectrum, as shown in Fig. 2(b).Note that the spectrum becomes gapless for the staggered flux configuration with θ = π/2, and remains gapped for any other choice of θ.The spinon spectrum is doubly degenerate for the spin-up and spindown components.The gapped bands acquire a nonzero Chern number in the uniform flux configuration.Using the link variable formulation [36], we obtain the total Chern number distribution for the bands as C = {2, −2} in the cQSL phase.Therefore, it is expected to host chiral spinon edge modes and exhibit quantized Hall thermal conductivity at low temperatures [37].

III. ANALYSIS AND RESULTS
Now we discuss the main results of this paper by focusing on the electromagnetic signatures in the cQSL phase.Despite a charge-neutral Mott insulator, the virtual hopping of electrons leads to a nonvanishing expectation value of the charge fluctuations and circulating loop currents in the cQSL phase [27].In fact, such features are expected in spin liquid systems [38][39][40].The relevant operators for the charge fluctuations and loop currents in the TLHM read [27] δ ρi,jk = e where ⟨ijk⟩ denotes an elementary triangle in the lattice, e is the electronic charge, rij is the unit vector along the bond ⟨ij⟩.The forms of δ ρi,jk and Îij,k are uniquely determined by the transformation of these quantities with respect to the following symmetry operations: SU(2) spin rotation, TRS, and inversion operation.We now compute the expectation values of the above operators in the spinon ground state.In this regard, we construct the real space spinon Hamiltonian on a finite system of linear size L = 30 and obtain the eigenvalues of the corresponding eigenfunctions of the 2L 2 × 2L 2 Hamiltonian [32].At first, we rewrite the above operators in Eq. (8a), and Eq.(8b) in the spinon degrees of freedom using the same mean-field decomposition as in Sec.II.For explicit numerical analysis, we need to fix the mean-field parameters.For subsequent analysis in this section, we work in units of t = 1.This leads to a solution of m 0 in terms of J, and Jχ from Eq. (7) as where ρ 0 = e 8m0t 3 U 3 , and are parameters that depend on the amplitude of the mean-field.Note that we added the contributions of the spin degrees of freedom in Eq. ( 9) [32] because of the degenerate spin bands and hence skipped the spin indices.
To obtain the total charge fluctuation and the loop current for a particular site or a bond, we need to add the contributions of all the shared triangles [27,40].Utilizing the mean-field expressions in Eq. ( 9) for the relevant operators, we calculate their expectation values in the spinon ground state [32] for a finite system, as mentioned before.The numerical estimates converge beyond the linear size L ∼ 20.In the periodic boundary conditions (PBC), each isolated triangle leads to identical estimates for the charge fluctuation and loop current expectation values.Consequently, there are neither charge redistributions nor circulating loop currents in the cQSL ground state.However, we obtain novel localized charge profiles and loop currents around the system's edges in a finite system i.e. with open boundary conditions (OBC).The corresponding results are shown in Fig. 3(a).The arrows around the edge signify the magnitude and direction of the localized currents.All values are in units of 2|I 0 |.The magnitude of the loop currents is slightly larger (∼ 0.7145) around the corners [C 1 in Fig. 3(a)] which are formed by either an up or down triangle, whereas they are smaller (∼ 0.6338) around corners which are composed of both an up and a down triangle [C 2 in Fig. 3(a)].Note that the loop currents quickly saturate (∼ 0.6764) as we move away from the corners along the edges and are consistent with the inversion and C 6 rotation symmetries.
Similarly, a finite charge fluctuation redistributes localized charges around the system's edges, as shown by blue and red circles.In this case, all numbers are shown in the unit of 2ρ 0 .Like the loop currents, the charge profile quickly saturates away from the corners.The maximum charge fluctuations (+0.0217/ − 0.0185) happens around the corner C 1 , whereas the minimum fluctuation (+0.0201/−0.0163)occurs around the corners C 2 .
The key feature is that the smaller the number of shared triangles for a particular site or a bond, the more the corresponding charge fluctuations or localized currents are, respectively.Most interestingly, the charge separation around edges leads to the formation of a unique dipole moment distribution that can be observed experimentally.

A. Case of a localized spinon
The cQSL supports spinons as its low-energy excitation.At the sample edge, there exists a gapless chiral spinon edge mode due to the non-trivial topology of the spinon bands.However, the spinon excitations are gapped inside the bulk.In a clean system with translational invariance in bulk, there are no charge fluctuations or loop currents in bulk [see Fig. 3(a)].Here we focus on an isolated/localized spinon excitation in bulk and discuss its associated electromagnetic responses.
In a clean cQSL, the lower spinon bands with spin up and down are fully occupied.To create a spinon hole, we demand that a specific spin in the spin Hamiltonian does not participate in the fractionalization into spinons.In the mean-field description, this can be achieved by setting the chemical potential for spinons at the pinning site [see Fig. 3(b,c)] to be high so that spinons will not occupy the defect site within the low-energy dynamics.This creates a localized spinon hole at the pinning site.Now, we consider a system as before with the defect formed by a large chemical potential at the pinning site as shown in Fig. 3(b,c), and impose periodic boundary conditions (PBC).Performing a similar analysis as in Sec.III [32], we notice a redistribution of the charge profile around the localized spinon hole, and a build-up of localized circulating loop current [see Fig. 3(b,c)].As before, all the numbers for charge and current are in units of 2ρ 0 and 2|I 0 |, respectively.We notice that the circulating loop current around the spinon hole site has the opposite chirality compared to the loop current flowing along the edge [see Fig. 3(a)] in the clean system with OBC.On the other hand, dipole moments formed by the charge redistribution are anti-aligned with the edge dipole moments in the clean system.In the latter case, we only focus on the nearest-neighbor location around the pinning site.Note that the charge profile quickly vanishes away from the pinning center.

B. DMRG calculations
To validate the above mean-field calculations, we next study the Hamiltonian Eq. ( 2) by using an unbiased DMRG method.Our DMRG calculations focus mainly on the 4-leg cylinder, retaining up to D = 4000 U(1) states.We summarize our DMRG results in Fig. 3(d,e) with U = 10t, i.e., deep in the cQSL regime.Here we show the left half of the cylinder for simplicity.We identify that the persistent electric current exists only close to the boundary, manifested by the nontrivial topology and spontaneous TRS breaking of the cQSL phase.The local electrical current quickly reduces from the boundary to the bulk.In the deep bulk, the net current is vanishingly small.
To create a spinon hole, we can add a local magnetic field H loc = V i (n i↑ − n i↓ ) to the Hamiltonian Eq. ( 2).(In practice, we add two local magnetic fields and ensure they are separated far away.One pinning point is shown as the green dot in Fig. 3(e), and the other is in the other half of the cylinder that is not shown here.)The local magnetic field pins the spin locally and forbids it from fractionalizing into delocalized spinons, therefore creating a spinon hole.Around the spinon hole, nonzero electric currents emerge in bulk.Importantly, around the pinned spinon hole, we identify the formation of a loop current (as indicated by the dashed arrow).It is also clear that the electrical charge distribution deviates from the average filling 1 required for the Mott insulator, as shown in Fig. 3(c,e).The general picture of this loop current and charge distribution associated with a spinon hole agrees with the prediction of the mean-field calculations in Sec.III A. Because of the finite size effect in the narrow direction in DMRG calculations, the current and charge distribution does not respect C 6 rotation symmetry along the spinon.
Numerical estimates: The DMRG results allow us to estimate the magnitude of the mean-field order parameter m 0 .Firstly, we provide a rough estimate of J, Jχ in Eq. (3) based on Refs.[12,13,32].Inserting characteristic values such as t = 1 eV, U = 10 eV, and χ ∼ −0.35 [13], we obtain J ∼ 0.37 eV and Jχ ∼ 0.15 eV.Here, χ is the nonzero chiral order parameter as defined in Fig. 1.Since the eigenfunctions of the Hamiltonian in Eq. ( 5) do not depend on the magnitude of t, we can compare the loop current magnitudes around the edge of the system obtained by DMRG with our mean-field analysis.Our estimates provide a mean-field amplitude m 0 ∼ 0.1.Utilizing this in Eq. ( 7), we obtain an order of magnitude for our phenomenological hopping parameter t ∼ 0.02 eV.Plugging in the magnitude (obtained by DMRG) of the enclosed loop current around our localized spinons, we estimate an emergent orbital magnetization ∼ 0.01 µ B , where µ B is the Bohr magneton.

C. Quantum field theory description
The orbital electrical current associated with a spinon can also be understood from the quantum field theory perspective, which sheds further light on the origin of the orbital electrical current.One hallmark of the QSL is the fractionalization of spins and the appearance of an emer-gent gauge field.Understanding the coupling between the emergent gauge field and the physical electromagnetic fields is crucial for the electromagnetic detection of the QSL.In terms of the parton description, the electron operator can be written as c σ = bf σ , where b is a boson operator that carries the electron charge e, and f σ is a fermionic spinon operator that carries the spin- b boson is gapped with g > 0 in the cQSL which is a Mott insulator.However, there is still a diamagnetic response in A − a due to the local current loop in the presence of a magnetic field, similar to Landau diamagnetism in metal, albeit the current loops are strongly localized.Since the b boson is gapped, we can integrate it out to obtain an effective Lagrangian as where the first term on the right-hand side is the Chern-Simon term obtained by integrating out f σ that fills topological Chern bands with a Chern number C (C = 1 in our model).Here χ b accounts for the diamagnetic susceptibility due to the gapped boson b, χ B is the susceptibility of the background [43].It is clear from the Chern-Simon term that a spinon carries π/C flux of a [44].The physical magnetic field associated with the emergent magnetic field, which can be seen from Eq. ( 11) by minimizing L with respect to B ≡ ∇ × A, is: Hence a spinon excitation induces an orbital electrical current with a total flux of B equal to χ b π/(χ b + χ B )C.

D. Dynamic spin-structure factor
In the previous sections, we established that spinon excitations in the cQSL phase carry orbital electrical loop currents and charges.Now, we proceed to investigate the electromagnetic response of a cQSL in terms of optical conductivity and Faraday rotation.Before considering the optical conductivity, which involves higher-order spinon correlation functions, we consider the standard dynamic spin-structure factor (DSSF) in the framework of spinon description.DSSF is an essential physical quantity that is routinely used as an experimental tool to probe the nature of the magnetic ground state and is defined as where N s denotes the number of sites, and q, ω denotes the probe momentum and frequency, respectively.With the two-sublattice structure as illustrated in Fig. 2(a), we first rewrite Eq. ( 12) in terms of spinon operators.The above expression simplifies upon utilizing the spectral representation with the weighted summation over the sub-lattice resolved spin-structure factors.The latter is written as [32] S ηζ (q, ω) = where {A, B} denotes the two sublattice degrees of freedom, and E n denotes the eigen energy of the n -th excited state.Note that we added the contributions from the degenerate spin up and down bands, and consequently skipped the indices as before.Rewriting in the diagonal basis and summing the sublattice degrees of freedom, we obtain the DSSF in our phenomenological cQSL.In Fig. 4(b), we show the DSSF profile.Note that we've adopted a normalization where the absolute maximum is set to unity for convenience.The excited state |n⟩ contains one pair of spinon hole and spinon excitation, or spinon exciton, as evident from Eq. ( 13).
We notice that apart from a relatively strong peak centered in a narrow region around the edge of the BZ at the K point, there are almost no sharp features within the BZ.The broad continuum in the BZ reflects the ab-sence of any long-range magnetic order, i.e., there are no well-defined magnon excitations at a given momentum q with energy ω.The relatively broad/diffused bands (illustrated by the white halos) correspond to a two-spinon continuum.At q = 0, the DSSF corresponds to the vertical spinon exciton, as is evident from Eq. ( 13).In this case, the wave function overlap between the wave function of the spinon hole in the occupied band and the spinon in the unoccupied band is zero at the same momentum and subsequently leads to a vanishing weight distribution around Γ point as seen in Fig. 4(b).To illustrate this, we also plot the scattering density of states g(ω, q) = ∑ k δ(ω − ε k+q − q) in Fig. 4(a), where there is a finite spectral weight around the Γ point.The absence of spectral weight around the Γ point is common to the cQSL phase in other lattices, viz.kagome [45,46].In re- ality, the fluctuations of the emergent gauge field around the mean-field saddle point mediate the attraction between the spinon hole and the spinon, which has been neglected in the present discussion.However, even in this case, the spectral weight around Γ point will vanish due to the zero overlap of the eigenfunctions [45].

E. Optical conductivity and Faraday rotation
Finally, we focus on the main result of our work by showing that optical responses below the Mott gap can be used to probe the emergent cQSL state in the TLHM [38,47,48].The longitudinal and the transverse optical conductivity in this regime become nonvanishing because of the finite electronic polarization.Following the work by Bulaevskii et al. [27], we obtain the corresponding expression for a three-site problem as where a is the lattice constant, and t, U are the parameters defined as before in Eq. (2).The above two expressions are particularly relevant as we deal with a triangle lattice.However, note that within a lattice framework, we need to add the contributions of all the triangles surrounding a particular site i to obtain the total polarization P. The latter naturally couples to an external electric field as −P ⋅ E(t).Consequently, the associated optical conductivity within the linear response theory reads [48][49][50] where |ψ 0 ⟩, and |ψ n ⟩ are the ground and excited states, respectively, H |ψ n ⟩ = E n |ψ n ⟩ ∀n ∈ {0, 1, 2, . ..},V is the volume, and ω n = E n − E 0 , where E 0 is the energy of the ground state.Note that the above expression is valid in the frequency regime much less than the energy scale (U ) associated with the charge gap in the Hubbard model, i.e., ̵ hω ≪ U .Additionally, broken TRS in the chiral phase immediately implies non-vanishing off-diagonal components (a ≠ b).This leads to a finite MOFE signal proportional to the real part of the transverse optical conductivity defined in Eq. ( 1).
We proceed as before in Sec.III D by rewriting the polarization operator in terms of spinon degrees of freedom.Readers are referred to Ref. [32] for the details of the calculations.However, in stark contrast to the DSSF analysis, here we need to consider the correlation functions involving eight spinon operators as is evident from Eq. ( 15) [32].We perform numerical integration in Mathematica with a quasi-Monte Carlo routine and obtain the transverse and longitudinal optical conductivity as a function of the frequency as shown in Fig. 5.Both the real (σ ′ ) and imaginary (σ ′′ ) part of the quantities are shown in panel (a) and panel (b), respectively.Similar to Sec.III D, we adopted a normalization in which the absolute maximum of the quantities is set to unity.
We notice that σ ′ xy (ω) changes sign at a frequency ω 0 ∼ 9 t that is almost twice the spinon gap around the BZ edge at the M point. of current experiments [51].

IV. DISCUSSION AND CONCLUDING REMARKS
This paper provides extensive mean-field analysis for the electromagnetic response of a cQSL phase.We started from a phenomenological cQSL Hamiltonian as in Eq. ( 3) and analyzed the spectrum of fractionalized excitations in terms of spinon mead field theory.Despite being deep inside the Mott insulator regime, where the charge degrees of freedom are gapped, we obtain a nonvanishing electrical loop current distribution and charge fluctuations associated with a localized spinon excitation.Additionally, we performed unbiased DMRG calculations in the triangular lattice Hubbard model at the intermediate coupling regime, where the cQSL is stabilized.The DMRG results confirm the physical picture of the parton mean-field results, where both approaches provide similar structures of the loop currents and charge redistributions in the cQSL phase, as illustrated in Fig. 3.The DMRG calculations further allow us to estimate the magnitude of the electrical charge and orbital current associated with a spinon excitation.Assuming a typical value of t = 1 eV and U = 10 eV, we estimate the electrical current and charge around the localized spinons to be around 17 µA, and ±0.1% of e, respectively.In addition, we performed quantum field theory analysis to unravel the connection between the spinon excitation and emergent and physical gauge fields, which clearly shows that a flux of the physical magnetic field dresses a localized spinon.
The electromagnetic characteristics of spinon excitations immediately imply a nonvanishing optical response in the cQSL.We compute the optical response functions by focusing on the optical conductivity.The nonvanishing transverse optical conductivity σ ′ xy (ω) below the Mott gap can be considered a smoking gun signature of the underlying chiral nature of the QSL.Since a finite σ ′ xy (ω) signifies a non-zero Faraday rotation angle Θ F , our predictions can be directly tested by suitable optical techniques such as MOFE or Kerr effect.Since σ ′ xy (ω) changes sign as the frequency increases, an experimental signature of cQSL would be to see if, as a function of incoming photon frequency, the Faraday angle changes sign or not.For completeness and as an intermediate step, we also analyze the dynamic spin-structure factor of the cQSL as illustrated in Fig. 4(b).The absence of sharp features signifies no well-defined magnon excitations in the QSL.
In cQSL, each unit triangle carries an orbital current.However, this orbital current cancels in the bond shared by two neighboring triangles for a translationally invariant system.This cancelation is not perfect in the presence of impurities or near edges, leaving finite orbital magnetization localized around impurities.Therefore, the orbital magnetization localized around impurities already serves as a signature of time-reversal symmetry breaking in QSL.This defect-induced orbital magnetization can be distinguished from spinons, which are dynamical excitations (despite being gapped) of cQSL.Depending on the protocol to tune the system into the cQSL, spinons can be created at different system locations, and the protocol can control their density.In contrast, the orbital magnetization localized around impurities does not depend on the protocol.
Compared to our previous theoretical work on Kitaev materials [40], here, TRS is spontaneously broken due to considerable charge fluctuations in a Hubbard model at intermediate coupling strength.As noted in our quantitative estimates for the loop current or associated charge polarization, the latter translates into a larger electromagnetic response.Note that the associated gauge structure for the cQSL in the TLHM is U(1), whereas the Kitaev spin liquid has a Z 2 gauge structure.
In summary, we show that spinon excitations in cQSL carry an electrical charge and orbital current, despite the system being a Mott insulator.Such an electromagnetic response can be detected experimentally using the MOFE or Kerr effect.Therefore, our work provides a clear electromagnetic signature of the cQSL, which helps determine the nature of nonmagnetic states observed in certain materials realizing the triangular lattice Mott insulator.

V. ACKNOWLEDGEMENT I. HUBBARD MODEL TO CHIRAL SPIN MODEL
In this section, we outline the steps for phenomenologically obtaining the chiral spin model starting from a Hubbard model on a triangular lattice at half-filling.The corresponding Hamiltonian is written as where ⟨ij⟩ corresponds to the nearest-neighbor tight-binding model on a triangular lattice, σ corresponds to the spin degrees of freedom, and U is the strength of the local Hubbard repulsion.In this case, the low-energy effective spin Hamiltonian in the strong-coupling limit (U ≫ t) can be obtained through Schrieffer-Wolff transformation (SWT) [S1-S3] as eff + H eff , H where the exchange couplings are given by [S3] = − In a previous theoretical work [S4], it was shown that the ring exchange term leads to an induced chirality in the low-energy spin dynamics of the TLHM.It is worth mentioning that such a flux phase was previously pointed out in triangular lattice material κ-(ET) 2 Cu 2 (CN) 3 by Motrunich within the mean-field description in Ref. [S5].Motivated by these studies, we consider the following chiral spin liquid model as where J, and Jχ are provided in Ref. [S4] as R , Jχ = 3J ) , where χ = ⟨S i ⋅ (S j × S k ⟩ is the non-vanishing chiral order parameter in the emergent chiral QSL state.Self-consistent density-matrix renormalization analysis in Ref. [S4, S6] shows that it is non-vanishing in a wide region between U c1 ∼ 9U and U c2 ∼ 10.75U .Therefore, the parameters of the phenomenological Hamiltonian in the main text are directly related to the parameters of the original Hubbard model.

II. PHENOMENOLOGICAL DESCRIPTION
We assume that the underlying fractionalized excitations are spinons, and correspondingly rewrite the spin degrees of freedom as [S7] Here, f † iα creates a spinon at site i with spin α, and σ denotes the vector of the Pauli matrices.We utilize the product relation for the Pauli matrices as Plugging this back into Eq.(3) of the main text, we first obtain the Heisenberg part as where n i = f † iα f iα and we assume the summation over the repeated indices unless explicitly mentioned.We note that the spinon description manifestly enlarges the physical Hilbert space.To remain in the physical Hilbert space, we utilize the half-filling constraint per site as ∑ α f † iα f iα = 1.However, in this work, we assume this constraint to be loosely applicable in an average way in the spirit of mean field theory.Next, we perform only particle-hole mean-field decomposition to rewrite the Hamiltonian in Eq. (S7) as where we have ignored the last term in Eq. (S7), and introduced a mean-field ansatz as m ij = ⟨f † iα f jα ⟩.Similarly, we rewrite the chiral term as [S8] The total Hamiltonian is then H pheno = H heisen + H chiral .Without going into a self-consistent mean-field analysis, we assume a particular form of the mean fields and benchmark our analysis with our unbiased DMRG calculations.
Assuming translational invariance, we choose m ij = m 0 e iϕij , where m 0 is the amplitude of the order parameter and ϕ ij are the bond-dependent phases.Now ignoring the amplitude and phase fluctuations, we can write the Hamiltonian as Combining the hopping phases ϕ ij 's, the above Hamiltonian can be written in a compact form as ., where t is the spinon hopping amplitude and is related to our phenomenological parameters as where ψ ij + ψ jk + ψ ki = Φ 0 with Φ 0 being total flux enclosed within a single triangular plaquette.At this point, all ϕ ij /ψ ij remains undetermined.Now we utilize Lieb's theorem [S9] to determine the phases.According to the theorem, a fermion hopping on a bipartite lattice realizes the ground state with π-flux square plaquettes.Consequently, we consider a doubled unit cell such that one up and one down triangle jointly form the rhombus-like bipartite unit cell as shown in Fig. 1(a), and impose a π-flux in the doubled unit cell.This still leaves us with various choices for the bonddependent phases ψ ij , ψ jk , and ψ ki forming the triangular plaquette ⟨ijk⟩.A generic choice of the bond-dependent phases is shown in Fig. 1(a), where θ is some arbitrary angle specifying whether both the up and down triangles have the same flux π/2 or some staggered flux configurations as π/2 ± θ (both adding to π in the rhombus-shaped unit cell).

A. Topological spinon bands
We now move on to compute the spinon bands within the phenomenological flux phases in the triangle lattice.First of all, the primitive and the reciprocal lattice vectors of the original triangular lattice are given by where a is the lattice constant.The corresponding nearest-neighbor vectors as given by We also show the corresponding Brillouin zone (BZ) in Fig. S1(a) with the high-symmetry points as Since the previous flux configuration doubles the unit cell as a 1 → 2a 1 , the tight-binding Hamiltonian in the sub-lattice basis (see Fig. 1(a) in the main text) is written as (note that we consider the uniform flux configuration with θ = 0) Translating into the momentum space, we obtain The above Hamiltonian can be written in a compact form with Pauli matrices as H eff = −2t ∑ k d k ⋅ σ where The gapped spinon spectrum is obtained by diagonalizing the Hamiltonian in Eq. (S16).The dispersion is given by (see Fig. 1(b) in the main text).The spectrum remains gapped for any other choice of θ, except at θ = π 2 when the gap closes as depicted in Fig. 1(c) in the main text.We computed the Chern number in the gapped phase using link variable method [S10], and find the Chern numbers for the bands to be ±2 (upon adding the spin degenerate bands) [see Fig. S1(b)].

B. Orbital electrical current and charge fluctuation in the gapped phase
In this section, we first provide the steps leading to an emergent non-vanishing loop electrical current distribution in the CSL phase.The current operator in the single-band Hubbard model reads as [S11] where rij is the unit vector connecting two sites i, j, and the localized current flows within a triangular loop.Rewriting in terms of the spinons, we obtain Note that in the last line, we removed the spin-label.Since the spinon bands are degenerate in the spin degrees of freedom, we have added the contributions from both spin channels.In a similar spirit, we can re-express the charge fluctuation operator in spinon language as [S11] where we again added the spin degeneracy in the last line of the above equation.
As we are interested in the real space loop current and chare fluctuation profile, we consider a real-space calculation to evaluate the loop current expectation values.The Hamiltonian in Eq. (S16) is therefore written on a 2D triangle lattice of size L×L unit cells.Since each unit cell contains two sub-lattice sites (A & B), there are 2L 2 spinon operators f † iσ in the system.The lattice vectors are chosen as a 1 = (2, 0), and a 2 = (1/2, √ 3/2).The spinon operators at a site i is written as In terms of the spinon vector f † , we can write the the Hamiltonian in Eq. (S16) as H eff = f † H f , where H is written as a 2L 2 × 2L 2 matrix.On the diagonal basis, we can rewrite the spinon operators as where U is the diagonalizing matrix for Hamiltonian H, and ζ u/d (m, n) are the diagonal spinon operators corresponding to the spectrum as shown in Fig. S2.In panel (a), we show the band dispersion with a topological gap ∆ in the PBC, while in panel (b), the spectrum in the case of OBC is shown with an edge mode inside the bulk gap ∆.The spectrum for a localized spinon at site r 0 inside the bulk is shown in panel (c).To realize the latter scenario, we impose a large onsite chemical potential V 0 at the site r 0 within the unit cell.Due to the large energy, this specific site will not host any spinons and can be thought of as a localized spinon hole.The physical situation might be some empty defect sites, or some magnetic impurity sitting inside the bulk of the system.

Analysis of the expectation values
Once we know all the eigenenergy and the eigenstates of the Hamiltonian, it is straightforward to obtain the average of the loop current operator in the ground state.The target quantity of our interest is ⟨f † α (R mn )f β (R ′ m ′ n ′ )⟩ for {α, β} ∈ {A, B}, where R mn denotes the position of the site at R mn = mâ 1 + nâ 2 .The analysis goes as follows where ⟨ζ † l ζ l ′ ⟩ = δ ll ′ , and the 2L 2 × 2L 2 matrices W RR ′ are defined as follows Here, B RR ′ are L 2 ×L 2 matrices corresponding to the non-zero connections allowed by the orientations of the triangles.The corresponding results for the charge fluctuation and loop current distribution are shown in Fig. 3(a-c) in the main text.

III. DYNAMIC SPIN STRUCTURE FACTOR
In this section, we provide the details of the analysis of the dynamical spin structure factor (DSSF).The latter is defined as where i, j corresponds to the position of the unit cell containing two sub-lattice sites.Note that the unit-cell i has two sub-lattice sites labeled by A, and B. Consequently, we can rewrite the above equation as dte iωt [ ⟨S i,A (t) ⋅ S j,A (0)⟩ + ⟨S i,B (t) ⋅ S j,B (0)⟩ + e iq⋅a1 ⟨S i,A (t) ⋅ S j,B (0)⟩ + e −iq⋅a1 ⟨S i,B (t) ⋅ S j,A (0)⟩ ], where we have explicitly written down the DSSF in sub-lattice resolved coordinates.It is straightforward to show next that a typical sub-lattice resolved term is given by ({η, ζ} ∈ {A, B}) S ηζ (q, ω) = ∑ where in the last line, we have summed over the degenerate spin degrees of freedom and omitted the spin indices, and |n⟩ corresponds to an excited eigenmode with energy E n .We perform the numerical integration in Mathematica and the corresponding plots are shown in Fig. 4 in the main text.We approximate the delta function δ(x) as δ(x) = 1 π γ 2 x 2 +γ 2 , and considered γ = 0.1 for numerical purposes.

IV. OPTICAL CONDUCTIVITY IN THE CSL PHASE
Finally, in this section, we provide the details of the analysis for the transverse and longitudinal optical conductivity in the CSL phase.Since the parent compound is a Mott insulator, we do not have any mobile charges; however, the charge fluctuations in the insulating phase will lead to electrical polarization which couples to the external electric field and lead to finite optical conductivity.The corresponding electrical susceptibility is given by where P = P x x + P y ŷ is the total polarization in the system.For a triangular plaquette, the corresponding expression can be obtained from the charge fluctuation operators [S11].Here, we consider a single site at i 0 embedded in the lattice as shown in Fig. S3.The polarization for each triangle is thereafter written as Now, we can add all these contributions to obtain the total polarization per site i 0 as We further utlize the Pauli matrix identities and rewrite the above expression in the spinon language as where δ αβγδ = 2δ αδ δ βγ − δ αβ δ γδ .Consequently, we can obtain the total polarization which is the sum over two sub-lattice polarizations as )

FIG. 2 .
FIG. 2. (a) Phenomenological spinon model on a triangular lattice with bond-dependent hoppings and a two-sublattice unit cell illustrated within the orange-dashed box.The hopping phases allow π-flux within each rhombus-shaped bipartite plaquette (see the main text for more discussion).For θ = 0, both the up and a down triangle forming the rhombus acquire uniform π/2 fluxes, where the flux configuration is staggered for any nonzero θ.The spinon spectrum for the uniform (gapped, θ = 0) and the staggered (gapless, θ = π/2) flux configuration are shown in panels (b) and (c), respectively.

FIG. 3 .
FIG. 3. (a) An illustration of the localized loop current and charge distributions around the edge of a finite system of linear size L = 30 (in open boundary condition) within the mean-field spinon description of the spin model in Eq. (3).For illustrative purposes, we do not show the explicit distribution of the loop currents within the bulk.Note that the loop current and charge fluctuation quickly vanish after a few lattice spacings inside the bulk.The loop currents (b) and charge distribution (c) around the localized spinon hole site were obtained in periodic boundary conditions with the same system size.The red and blue colors signify the opposite signs of charge redistribution.The numbers are presented in the unit of 2I0 and 2ρ0, respectively (see the main text).(d) and (e), Plots of the local electric currents on the triangular Hubbard model for (d) without and (e) with a spinon hole located at the position labeled by the green color obtained by the DMRG calculations.The loop current emerges around the local magnetic field.The red arrows represent the direction of the loop current.The numbers around the bonds label the absolute value of the current in the unit of et/ ̵ h.(f) Plot of the charge redistribution around the local spinon hole on a finite-size system (illustrated by the green region) obtained by the DMRG calculations.For clarity, the numbers are in the unit of 10 −4 e.

1 2 .
In cQSL, f σ fermions form Chern bands as was shown in Sec.II.The fractionalization dictates that the charged boson is coupled to both the physical gauge field A and an emergent gauge field a as b → b exp[i(A − a)], while the spinon is coupled only to the emergent gauge field as, f σ → f σ exp(ia).The effective low-energy Lagrangian for the b boson has the standard Ginzburg-Landau form (we use the unit ̵ h = e = c = 1)[40][41][42]

FIG. 4 .
FIG.4.The normalized scattering density of states g(q, ω) (a) (see definition in the main text) and the dynamic spin structure factor S(q, ω) (b) along the high symmetry points of the triangular lattice Brillouin zone (BZ).The analysis is performed neglecting any spinon interactions which are present when the dynamics of the emergent gauge field are considered explicitly.

FIG. 5 .
FIG. 5.The real (a) and the imaginary (b) part of the normalized transverse (dashed line) and longitudinal (solid line) optical conductivity as a function of the frequency of the incident light.

FIG. S1 .
FIG. S1.(a) The Brillouin zone of a triangular lattice with the high-symmetry points marked by the filled circles.(b) The spectrum of gapped spinons within our mean-field approximation with the staggered phase θ = 0.The topological Chern numbers are labeled on the two bands.Each spinon band is doubly degenerate in the spin degrees of freedom.

)
FIG. S2.(a) The eigenvalue distribution of the real-space Hamiltonian (PBC) as a function of the system size.The topological protection manifests into a gap in the eigenvalue distribution.(b) The eigenvalue distribution of the real-space Hamiltonian as a function of the system size in open boundary conditions with the edge modes in the gap.(c) The eigenvalue distribution of the real-space Hamiltonian as a function of the system size for the triangular lattice in periodic boundary conditions.We set a local chemical potential V0 at site r0 for the A-sublattice.For very large V0, there are two additional states which are shown by the symbols around zero energy and the other one at higher energy dictated by V0.