Spin Liquid versus Spin Orbit Coupling on the Triangular Lattice

In this paper, we explore the relationship between strong spin-orbit coupling and spin liquid physics. We study a very general model on the triangular lattice where spin-orbit coupling leads to the presence of highly anisotropic interactions. We use variational Monte Carlo to study both $U(1)$ quantum spin liquid states and ordered ones, via the Gutzwiller projected fermion construction. We thereby obtain the ground state phase diagram in this phase space. We furthermore consider effects beyond the Gutzwiller wavefunctions for the spinon Fermi surface quantum spin liquid, which are of particular importance when spin-orbit coupling is present.


Introduction
Quantum spin liquids (QSLs) are exotic phases of correlated electrons possessing highly entangled ground states, exotic fractionalized excitations, and typically, the absence of any magnetic order [1,2]. Historically, studies of QSLs focused on spin-rotationally invariant Heisenberg models, but in recent years, strongly anisotropic interactions arising from spin-orbit coupling have come under focus [3]. In the famous Kitaev honeycomb model, bond-dependent interactions lead to an exactly solvable model with a spin liquid ground state [4]. Remarkably, it was later shown that these directional interactions can be generated in real materials when spinorbit effects are present [5,6]. In turn, this has led to the recent discoveries of many candidate 'Kitaev' materials and has paved the way for the study of spin liquid physics in spin-orbital systems. One recent example of particular interest is the material YbMgGaO 4 [7,8,9,10,11]. This system very likely contains directional interactions of significant strength. Moreover, thermodynamic and inelastic neutron scattering measurements have been interpreted as supporting a QSL state with a Fermi surface of neutral spin-1/2 excitations, "spinons", in this material.
Spin-orbit generated interactions invariably lead to a strong breaking of spin-rotation symmetry. A consideration of this symmetry in spin liquids can then reveal new and unexpected physics. One striking feature is that the lowered symmetry allows for new distinct spin-liquid phases which do not exist in the rotationally invariant case [12,13]. There exists a systematic method of classifying these phases, given by the so-called projective symmetry group (PSG)

The Model
In many physical systems, the spin and orbital degrees of freedom of the localized electrons are highly entangled. In these cases, when the rotation symmetry is broken by the surrounding crystal structure, the spin-rotation symmetry is broken as well. Superexchange processes then lead to the generation of highly anisotropic terms in the effective spin Hamiltonian. In these strongly spin-orbit coupled systems, lattice symmetry transformations are accompanied by an equivalent transformation in spin space. Following Ref. [15], we consider the Hamiltonian where γ ij = 1, e 2πi/3 , e −2πi/3 for bonds ij along the a 1 , a 2 , a 3 directions, respectively (see Fig. 1a)). This is the most general nearest-neighbor Hamiltonian which is invariant under the symmetry generators of the system: the translations T 1,2 along the a 1,2 directions, the sixfold roto-reflection S 6 within the plane of the lattice, the twofold rotation C 2 around a bond in the a 3 direction, and time reversal Θ. (Note that the threefold rotation C 3 = S 2 6 and the inversion I = S 3 6 are both generated by the sixfold roto-reflection. Conversely, the sixfold roto-reflection S 6 = C 2 3 I is a combination of a 120 • rotation and an inversion.) The symmetry generators are all discrete and act simultaneously in real space and spin space. In particular, they transform the coordinates x 1 , x 2 of a general lattice point r = x 1 a 1 + x 2 a 2 as while they transform the spin components (S x , S y , S z ) as T 1,2 : (S x , S y , S z ) → (S x , S y , S z ),  Figure 1: a) The three lattice bonds a 1 , a 2 , and a 3 . The two commensurate magnetic orders we consider are b) stripe order and c) 120 • antiferromagnetic order.
Importantly, the Hamiltonian does not generically have a continuous spin-rotation symmetry because the XXZ terms H ± and H z break the SU (2) spin symmetry down to an in-plane U (1) spin symmetry, while the remaining terms H ±± and H ±z further break the U (1) spin symmetry down to discrete spin symmetries that are intertwined with appropriate lattice symmetries.
It is helpful to write the H ±± and H ±z terms in a slightly different form to further expose the symmetries: wheren ij is the unit vector pointing from site i to site j. The term H ±± has a 'clock' structure where spins would like to align along the 120 • bond directions, and the term H ±z also has a bond dependent structure that incorporates theẑ direction. There are several cursory reasons one may expect to find spin liquid physics in this model. For one, due to its strong frustration, the triangular lattice has a long and storied history as a spin liquid candidate [21,22,23,24,25]. Beyond that, the form of the anisotropic part of H is highly reminiscent of the interactions in the Kitaev honeycomb model [4], where the direction-dependent spin-spin interactions frustrate the coupling in a way which renders all magnetic orders energetically unfavorable.

Generalities of parton wavefunctions
The ground state wave function in a quantum spin liquid is completely symmetric under all the symmetries of the Hamiltonian. The PSG gives a systematic classification of the allowed spin liquid phases under such a set of symmetries [14]. In the process, it also gives a construction of a representative wave function for each phase. It is a surprising fact that, in many cases, the number of allowed spin liquid phases increases as the symmetry is reduced [13,12]. Spin liquids are fundamentally defined by their fractionalized quasiparticle excitations, whose behavior can be described phenomenologically by a mean-field Hamiltonian. The PSG classifies the fractionalized symmetry by identifying the allowed form of the mean-field Hamiltonians. In general, these excitations can realize the symmetries of the original Hamiltonian in a nontrivial manner.
One starts by writing the physical spin operator S i in terms of fermionic parton operators: The parton operators f iσ , f † iσ live in a larger Hilbert space than the spins S i . To remedy this, one must also include the strict gauge constraint on the allowed states: In this paper, we enforce Eq. (6) at the level of the wave function. This is accomplished by applying the Gutzwiller projection operator P to a state |ψ 0 in the fermionic space: The projected wave function |Ψ lives in the proper Hilbert space of spins and, with a suitable choice of |ψ 0 , is highly entangled in real space. Furthermore, with some minor improvements, such an ansatz can be made to give variational energies which are competitive with the most state of the art 2D DMRG calculations [23]. For the state |ψ 0 , we take a "mean field" wavefunction, which is the ground state of some quadratic fermion Hamiltonian. The parameters of that fiduciary Hamiltonian then become variational parameters in the ansatz. When the fermions are allowed to hop in the mean field Hamiltonian, the partons become deconfined in the corresponding spin liquid phase. In general, the quadratic mean-field Hamiltonian can be written as where α = 0, x, y, z. A local gauge transformation, such as f iσ → e iθ i σ z f iσ , changes H mf but leaves the physical spin operator S i unchanged. Since the physical wave function is unchanged, all mean-field Hamiltonians related by such local gauge transformations must be equivalent. The parton Hamiltonian H mf can therefore ostensibly break the symmetries of H as long as there exists a local gauge transformation which restores the symmetry. In this case, we say that the quasiparticle realizes the symmetry nontrivially. The role of the PSG is to determine the set of allowed mean-field Hamiltonians which cannot be connected to each other by such a gauge transformation. Importantly, H mf is always invariant under some global transformations Φ → Φ · W , where W ∈ G. The group G ⊇ Z 2 of such global transformations is known as the 'invariant gauge group' (IGG) and determines the form of the gauge group around which fluctuations of the gauge field may occur. In this work, we consider U (1) spin liquids with IGG = U (1). A more complete study would also include Z 2 spin liquids (IGG = Z 2 ). However, even restricting to nearest-neighbor couplings, there are at least 18 different Z 2 mean-field ansätze. To avoid this complexity, we neglect these candidate QSLs for the present work. This is at least consistent with recent work on several related triangular lattice spin systems, for which the U (1) spin liquids have proven to have competitive energies [21,23]. Furthermore, the spinon Fermi surface QSL suggested by several previous papers for YbMgGaO 4 falls into the U (1) class.

Six specific parton states
The PSG classification of U (1) QSLs for the space group of our model was done in Ref. [18]. There are 6 distinct nearest-neighbor mean-field Hamiltonians: The ground state of each mean-field Hamiltonian defines |ψ 0 for the corresponding type of QSL. We distinguish two versions for each mean-field Hamiltonian H (n) mf , which differ only in the way translation symmetry is realized. In the A states, translation acts in the usual way as t ij = −1 for all nearest-neighbor bonds ij . Conversely, in the B states, translation acts nontrivially; this is achieved by setting t ij = ±1 such that the unit cell is doubled and a π flux is thread through every other triangle. In the A1/B1/A2/B2 cases, there are no variational parameters (since the overall scale of the Hamiltonian leaves its ground state unchanged), while in the A3/B3 cases, there is a single variational parameter λ/t.
We note that, importantly, the spinon band structure determines the physical properties of the wave functions and that it is gapless in all 6 states. This is necessary for a U(1) spin liquid to be stable in two dimensions. We now discuss some aspects of these states.
The (A1) state has no mixing between the up and down spin states. In order to satisfy the constraint f † i f i = 1, the band structure then must contain a large Fermi surface. We refer to this state as the uniform Fermi surface (uFS) or spinon metal state. Notably, although the microscopic Hamiltonian H has only discrete symmetries, the mean-field Hamiltonian of this uFS state is spin-rotationally invariant. This accidental "emergent SU(2) symmetry" is surprisingly robust, and is not an accident of assuming a nearest-neighbor form for H mf . In fact, the PSG does not allow any spin-dependent terms [which would break SU(2) symmetry] in H mf , even for hoppings of arbitrary distance. The argument for this hinges on the fact that both time-reversal (Θ) and inversion (I) symmetries act trivially in this class. First, the operators which implement these symmetries both involve a complex conjugation, timereversal by definition and inversion due to a site-exchange which corresponds to a Hermitian conjugation. Then, since spin is even under inversion and odd under time reversal, it is odd under their combination, and so a spin-dependent term with any complex coefficient is forbidden in the presence of such a combined symmetry. The (B1) state also has no mixing of the spin states, but translations act nontrivially on the spinons. The unit cell is then doubled and the band structure contains two Dirac cones. We therefore refer to this state as the Dirac spin liquid state. The uFS and Dirac states are the two U (1) spin liquids that can also occur in rotationally invariant systems. Note, however, that spin-dependent quadratic terms are not generically prohibited in the case of the (B1) state and that they in fact appear at the level of third-nearest-neighbor hoppings.
The (A2) and (B2) states are called the 120 • clock spin liquid (ClSL) and the 120 • clock + π spin liquid (ClπSL), respectively. These ansätze do mix the spin flavors and orbital degrees of freedom by including bond dependent hoppings. The band structures in both cases contain protected Dirac cones at the Γ, M , and K points in the Brillouin zone.
The (A3) state, called the Rashba spin liquid (RSL), also has Dirac cones at the Γ, M , and K points when λ = 0 or t = 0, and a gap opens at the Γ point for intermediate values of λ/t. Finally, the (B3) state, called the Rashba + π spin liquid (RπSL), contains 4 bands and a small Fermi surface for intermediate values of λ/t.

Energetics of PSG Wave Functions
The PSG method gives us the full set of allowed free fermion wave functions that are invariant under the symmetries of the system once the gauge constraint, Eq. (6), is enforced. It tells us nothing, however, about the energies of these wave functions. The PSG gives us a starting ansatz, but is completely agnostic about which state may actually be the ground state.
One simple way to proceed is to work directly with the single particle wave functions by satisfying Eq. (6) on average: 1 However, such a mean field approach requires an infinite number of approximations, the resulting wave functions do not even live in the proper Hilbert space, and thus it cannot give reliable energy estimates. Instead, we carry out a variational analysis based on the fully projected wavefunctions in Eq. (7). We calculate the variational energy E s = Ψ s |H|Ψ s , where s indicates one of the six QSL ansätze.
The results are highly constrained by how the projective symmetries are implemented in the given mean-field Hamiltonian. In particular, the uniform Fermi surface and Dirac spin liquid states are completely SU (2) invariant, and therefore the expectation values of the J ±± and J z± terms vanish in these states. Similarly, while both the 'clock' and 'Rashba' Hamiltonians have some spin-orbit terms, only the Rashba Hamiltonians include spin-orbit terms both within and perpendicular to the xy plane. Consequently, the 'clock' wave functions also yield vanishing expectation values for the J ±z terms.
We performed a variational Monte Carlo simulation and measured the energies of each of our trial wave functions on finite size lattices for system sizes up to N = 32 × 32 sites. Each mean-field wave function, when projected, gives a different pattern of entangled spins, giving rise to different spin correlations. When λ = 0, none of the wave functions have any free parameters. Setting J ± = 1 and scaling to the thermodynamic limit, the corresponding energy densities are then given by A few observations are apparent. First, we see that the (ClπSL) and (RπSL) ansätze are never competitive energetically in our regimes of interest. While the Dirac state has the lowest energy at J ±± = 0, the clock and Rashba spin liquid states become energetically favorable for large positive and negative J ±± , respectively. The Rashba states (and only the Rashba states) have an energy which is modified by including λ = 0, which is beneficial only when J ±z = 0. In this case, we determine the optimal Rashba state for a given value of J ±z by numerically minimizing the energy with respect to λ/t. The results for a full comparison of energies are presented in Fig. 2a), which shows the state of lowest variational energy amongst the 6 QSLs for all values of J ±± and J ±z . (Note that the phase diagram is qualitatively similar for all values of J z .) Looking ahead, it has been suggested [11] that next-nearest-neighbor interactions may be important in stabilizing a spin liquid ground state for our model. We therefore also looked at the variational energies of our ansätze when XXZ-like next-nearest-neighbor interactions are added (see Eq. 13 in Sec. 3.4.2). In Fig. 2b), we plot the lowest energy states as a function of the next-nearest-neighbor coupling J 2 for J ±z = 0. Notice that the Fermi surface state only becomes competitive in energy for very large next-nearest-neighbor coupling.

Parton formulation of ordered states
The PSG wave functions can be used as a starting point on which magnetic order can be added. This is done by adding a site dependent magnetic field h i to the mean-field Hamiltonians, which define our trial states: Magnetic order can be induced in this way on top of any of the 6 QSL states. In practice, the lowest energies are found by using H Zeeman term in this case fully gaps the partons. Consequently, the usual Polyakov argument, which applies to an emergent U (1) gauge theory with fully gapped Dirac fermions in two dimensions, implies that monopole instantons proliferate and the Dirac spinons are confined. Thus, the projected wavefunction built from H mo describes a state adiabatically connected to a conventional magnetically ordered one. If | h i | → ∞, Eq. (11) describes classical magnetic order with | S i | = 1/2 on each site. If instead a finite field is used, the value of the magnetic moment can be greatly reduced. In general, the energy should be optimized with respect to the full set of Zeeman fields h i on all sites. In practice, such an optimization would have too many parameters. Instead, we guess an appropriate pattern for these fields, and then optimize |h|/t to give the lowest variational energy. For example, in the Heisenberg limit, we choose the field to have a constant magnitude but an orientation with a three-sublattice structure of total vector sum zero (the symmetry pattern of the 120 • state): where q, |h| and a phase φ are variational parameters. In this case, the optimal magnetic field of our simple ansatz gives a staggered magnetic moment | S i | ≈ 0.30, while the corresponding DMRG calculations for the triangular-lattice Heisenberg model find a staggered magnetic moment M ∼ 0.20 [24]. Including local correlations in our variational state, for example, by including Jastrow factors, will in general reduce the value of S further. It is interesting that our PSG analysis provides a general way to construct any ansatz satisfying the constraint of Eq. (6), even allowing us to construct energetically competitive magnetic states in addition to giving a general classification of all spin liquid states.

Extended model
Implementation of the above method shows that the nearest-neighbor Hamiltonian H nn in Eq. (1) is dominated by magnetic order. To find actual spin liquid physics, we therefore extend the model to include second-and third-neighbor interactions. Keeping the same relative XY anisotropy, we study the Hamiltonian: To avoid complications involving canted magnetic orders, we restrict our attention to the case of J ±z = 0. With this in mind, in this section, we undertake the somewhat ambitious goal of describing the entire four-dimensional phase diagram in terms of the free parameters J z , J ±± , J 2 , and J 3 , all relative to J ± = 1.
We first review what is already known about the ground state phase diagram of Eq. (13): • In the absence of second-and third-neighbor interactions (J 2 = J 3 = 0), the Luttinger-Tisza analysis of Ref. [15] tells us the magnetic order when S is treated as a classical vector. In that case, there is a phase transition from the 120 • staggered antiferromagnetic state [ordered at wavevector q 120 = ( 4π 3 , 0)] at small |J ±± | to a striped phase [ordered at wavevector q s = (0, 2π • There is also a great deal of literature on the quantum J 1 − J 2 model (J ±± = J 3 = 0), in the Heisenberg limit (J z = 2J ± ) [25,24]. In this case, growing evidence suggests that a spin liquid phase interpolates between the 120 • phase for small J 2 and the stripe phase at large J 2 .

VMC results
The advantage of using variational Monte Carlo with simple trial wave functions is that we are able to explore a huge phase space of our Hamiltonian. We consider several ansätze for magnetic order, taking the Zeeman field in the form of Eq. (12) with wavevector q v = (q, 0) or q v = (0, q), where q, |h|, and a phase φ are variational parameters, which allows for both commensurate and incommensurate ordering. In practice, we find that the energies of all our ansätze, except for the striped phase with q s = (0, 2π √ 3 ), are independent of φ, even when the U (1) symmetry is broken by H ±± . In the stripe phase, we find that the minimum energy is always obtained for φ = 0 when J ±± > 0, giving the ordering pattern seen in Fig. 1b), and for φ = π/2 when J ±± < 0, which rotates all spins by 90 • . In Fig. 3, we present our result for the full quantum J z − J ±± − J 2 phase diagram. Notice that our results agree well with the previously understood limits. When J 2 = 0, the system acts very similar to the classical case, with a transition between the 120 • and stripe orders around J ±± ≈ 0.20 + 0.05J z . When a second-neighbor interaction is added, we indeed see that a Dirac spin liquid appears between the 120 • and stripe phases. This phase is stable for small J ±± , but both large J 2 and J ±± favor stripe order, leading to the triangular shape of the spin liquid regime which we see in Fig. 3. It is also notable that the extent of the spin liquid phase shrinks dramatically when J z is lowered from the Heisenberg point. This is in agreement with the DMRG results on this model in Ref. [20]. We are also able to go beyond this model to look at the effect of the third-neighbor XXZ interaction J 3 . Since both the second-and third-neighbor sites are separated by two lattice bonds, a simple superexchange picture implies that such a term would be present in materials with J 3 ∼ 0.5J 2 . We will see that the effect of such a term is to enhance the size of the spin liquid regime.
First, we present the results in the classical limit. When J ±± = 0, the system has U (1) symmetry and we can solve for the classical magnetic order using the Luttinger-Tisza method since any coplanar magnetic order with a single ordering wave vector will satisfy the hard constraint that S i = 1/2 on every site. The result is that, in addition to the usual 120 • and stripe phases, J 3 favors two additional incommensurate magnetic phases, with ordering wave vectors at (q, 0) and (0, q). These phases can be thought of as the incommensurate versions of the 120 • and stripe phases, respectively. A third incommensurate order with wave vector (q, q) also appears classically, but we will ignore this as such a phase never appears in the quantum case. The full classical phase diagram is shown in Fig. 4a) and is independent of J z since the ordering is always in the xy plane.
Our VMC results on the quantum model agree remarkably well with the classical phase diagram, considering we have used completely different methods. Fig. 4b) shows the results for J z = 1.0. We see that the shapes of the magnetic phases are largely the same as in the classical case, but the intermediate region where the phases meet is occupied by a broad spin liquid regime.
In Fig. 5, we show the full J 2 − J 3 − J ±± phase diagram for J z = 1.0. In addition to the presence of incommensurate magnetic order, the major feature of the data is that the spin liquid regime is enhanced with respect to the J 3 = 0 case. The third-neighbor interaction provides further frustration and finds stripe order particularly unfavorable. The spin liquid phase therefore survives to a large value of J ±± when J 3 is included.
As mentioned previously, more accurate energies can be found by adding further variational parameters to the wave function, such as allowing for Jastrow factors [26,27] or performing a small number of Lánczos steps [28]. However, we find that supplementing the PSG wave functions in this way only gives small improvements in the energies, leading to very small shifts of the phase boundaries. In section IV, we look at how we can make qualitative changes to the spin liquid ansätze.
In summary, our variational Monte Carlo calculation allowed us to explore a huge parameter space of the Hamiltonian in Eq. (13) and to obtain quantitative results for the ground state in each parameter regime. When a second-neighbor interaction is added, the Dirac spin liquid appears as the ground state between the 120 • and stripe phases. This phase shrinks dramatically away from the Heisenberg limit, but is in fact enhanced when a small third-neighbor interaction is included.

Perturbative Correction to the Wave function
In this section, we take a more phenomenological approach to studying a quantum spin liquid in the presence of strong spin-orbit coupling. We propose modifications to the mean-field ansätze which can be implemented numerically in the variational wave functions.
The plain mean field ansätze are limited in the amount of complexity they can accommodate. The main issue with the VMC simulation in this context is that the two most energetically competitive states, the Fermi surface and the Dirac spin liquid ones, possess too much symmetry. Our trial wave functions have no coupling between the spin and orbital degrees of freedom, which is a feature one would expect to find in the Hamiltonian's true ground state. Furthermore, according to the PSG analysis, no fermion bilinear operators inducing such spin orbit coupling can be added to the uniform Fermi surface Hamiltonian, not even at the further-neighbor level.
Instead, we formulate a method to incorporate many-body effects which modify our wave functions. Inspired by the path integral formulation for an interacting quantum field theory, we consider the variational state where H = H ±± is defined in Eq. (1). This form is reminiscent of the Lánczos algorithm, where applications of large powers of an operator project a trial state into the ground state of the given operator. Indeed, if we let α → ∞, this operator projects into the ground state of H. Instead, however, we take a slightly different approach, and let α be a small perturbation on P|ψ 0 , treating |Ψ as a variational wave function.
There have been previous works combining the Lánczos algorithm with variational Monte Carlo [23,28]. This proceeds by applying a finite number of Lánczos steps and working with the wave function |Ψ (n) = (1 + n p=1 α p H p )|ψ 0 , where the series is truncated for some small n, and the α p are left as variational parameters. While this works well if the initial state is very close to the ground state of H, it is less effective as a phenomenological tool. The reason is that corrections at any finite order n necessarily scale to zero in the thermodynamic limit. When calculating the correction to an expectation value using |Ψ (n) , "disconnected" powers of the Hamiltonian are subtracted off in the numerator, but not in the denominator. The normalization factor in the denominator therefore necessarily grows faster than the numerator with system size. Additional powers of n are then needed to compensate for this fact, but a fully extensive correction is only found at n ∼ N .
Instead, we have found that the best way to work with the wave function in Eq. (14) numerically is to implement the correction perturbatively in α, but to all powers in n. To do this, we realize that the expectation value of any operator with respect to our improved wave function can be written as where · · · 0 is the expectation value taken with respect to the unperturbed wave function P|ψ 0 . It is now possible to expand Eq. (15) analogously to diagrammatic perturbation theory. For any Hermitian operator O, the expanded correction reads The subtle difference is that now, by including all powers of n, all terms in the denominator exactly cancel the higher order "disconnected" pieces in the numerator. In the VMC calculation, this is expressed by the fact that H ij H k ≈ H ij H k as |(ij) − (kl)| → ∞. This way, we are able to measure, in our numerical simulation, many-body corrections to the wave function which survive in the thermodynamic limit. In principle, applying the operator exp[−αH] to our unperturbed trial wave function could cause a phase transition, and we would no longer be working with a spin liquid state. For small α, however, we expect that the spin liquid ground state should be stable to such a perturbation. In the spinon metal, in a similar vein to Fermi liquid theory, we expect that these terms only give a correction to the self-energy of spinons near the Fermi surface [29].

Correction to the Energy
To begin, we measure the correction to the energy of the Dirac and uniform Fermi surface states, which arises from including the spin-orbit interaction in our variational wave function. We can directly measure the first and second order corrections numerically.
For any operator O, we write the n th order correction to the expectation value O from applying exp[−αH] as α n O (n) . Expanding Eq. (16) gives In our case, H = H ±± and H ±± 0 = H ±± H ± 0 = 0. Therefore, the spin-orbit part of the Hamiltonian is altered at order α, while the rotationally invariant part is corrected at order α 2 : In Fig. 6, we show the resulting scaling of H (1) ±± and H (2) ± to the thermodynamic limit. The result is that the spinon metal is more susceptible, compared to the Dirac state, to energetically beneficial corrections to H ±± and less susceptible to detrimental corrections to H ± and H z . Putting this together, we find that the optimal value of the variational parameter is α min ∼ J ±± /(J ± + J z ), which gives an energy correction ∆E ∼ −J 2 ±± /(J ± + J z ). More precisely, we find that the energy densities after the lowest-order corrections are given by This implies that that the Fermi surface state becomes energetically superior to the Dirac state between J ±± = 0.57 at J z = 0 and J ±± = 1.54 at J z = 2.0. One caveat, of course, is that these values of J ±± may fall outside the perturbative regime. Also, while smaller J z appears to be more favorable for the spinon Fermi surface, this is also the parameter regime which is more susceptible to magnetic order.

Correction to the Spin Structure Factor
Studying the improved variational wave function makes it clear that the spinon metal state in a spin-orbit coupled environment has several unique properties, despite the fact that the mean-field Hamiltonian retains its rotational invariance. Taking our analogy to Fermi liquid theory seriously, the spin-orbit interaction gives a momentum and spin dependent correction to the self energy. This appears as a momentum dependent correction to the structure factor, which we can again measure directly in our simulation. We differentiate between the various spin polarized contributions to the spin-spin correlation function: The first-order correction to the correlation function is The results are shown in Fig. 7. The corrections to the spin-polarized structure factor are direction-dependent broad peaks at the M points of the Brillouin zone which appear at order α ∼ J ±± /(J ± + J z ). Therefore, in a spinon metal with spin-orbit coupling, spin-spin correlations when measured with different spin polarizations are direction dependent. This type of measurement could prove to be an important test to show both the presence of spinorbit interactions and the absence of spontaneous symmetry breaking. Similar directional peaks can be seen in related models when spin-orbit terms are directly included in the ground state ansatz [30]. We note that these kinds of direction-dependent structure factors have already been measured experimentally by resonant elastic x-ray scattering in the honeycomb lattice iridate Na 2 IrO 3 [31].

General considerations
Thermal transport measurements can be a powerful tool for studying magnetic insulators. The idea is to set up a thermal gradient ∇T (which is analogous to an electric field) and then measure the heat current j th in response to it (which is analogous to an electric current). Any heat current in the insulator must be carried by the emergent quasiparticles, giving us a probe of the low energy excitations. The thermal conductivity, κ, can be defined within linear response as The spinon Fermi surface QSL is unusual due to the large number of gapless excitations. This leads to a predicted linear T term appearing in the diagonal component of κ, similar to what one would expect in a metal. The deconfined spinons carry heat in the same way physical electrons carry charge in an electrical conductor. A major difficulty is that many degrees of freedom, most notably phonons, can contribute to the diagonal thermal conductivity, making the measurement challenging. The thermal Hall conductivity, however, given by the off-diagonal component of κ, should not contain a phonon term. Furthermore, as explained in Ref. [32], it is very difficult to find an effect generated by magnons on the triangular lattice due to a cancellation of the contributions from neighboring edge sharing plaquettes. A large nonzero thermal Hall conductivity could therefore be a strong indicator of exotic physics. Indeed, in Ref. [32], the authors also predict that a spinon metal would display such an effect. However, the reasoning is very subtle, depending on a coupling of the orbital motion of the spinons to the external electromagnetic field through the interaction with the internal gauge field.
Here, we argue that there exists a distinct contribution to the thermal Hall conductivity in the spinon metal which is unique to spin-orbit coupled systems and relies only on a Zeeman coupling to the external electromagnetic field. For itinerant fermions with conserved charge, the presence of spin-orbit coupling can lead to a nontrivial Berry curvature which may induce an anomalous component of the charge Hall conductivity, in the absence of any Lorentz force. This mechanism of anomalous Hall conductivity was explored intensely for Rashba two-dimensional electron gases and in many other models. In the following, we adapt this idea to study the thermal conductivity of the Fermi surface QSL state.
The U(1) QSL states studied here have an emergent conserved charge, which is the fermion number associated with the emergent U(1) gauge symmetry. Consequently, at the parton level, we can define a current associated with this charge, and we may consider, formally, the emergent conductivity tensor σ qp µν defined with respect to the emergent current and a potential coupling to the associated emergent charge density. This is not the true electrical conductivity, since this is an insulator, and it is also not the thermal conductivity. Thus we proceed in two stages. First, we consider the anomalous emergent Hall conductivity of the spinons. Then, we relate it to the more easily measurable thermal Hall conductivity (in principle, the emergent conductivity should also be measureable, but it is not obvious how to do so).

Effective quasiparticle Hamiltonian
At the mean field level, the emergent Hall conductivity can be extracted as an integral over the Berry curvature of the occupied spinon bands. Within the simple PSG wave function, the spinon metal is spin-rotationally invariant and therefore has zero Berry curvature. On symmetry grounds, however, we expect that a Hall conductivity should microscopically arise. To estimate it, we consider the 'improved' wave function, and infer a self-energy correction which breaks spin-rotational symmetry and induces a non-zero Berry curvature.
The Berry gauge field (Berry connection) is defined for a single particle system as where |u k is defined as in the Bloch wave function. The anomalous Hall conductivity is then given by where the first (line) integral is taken around the Fermi surface ∂S, while the second (area) integral is taken over the area S spanned by it. This physical quantity is invariant under U (1) gauge transformations, as is immediately evident from its expression in terms of the Berry curvature B(k) = ∇ k × A(k).
To obtain the Berry curvature, we suppose that the system is described by an effective quasiparticle Hamiltonian including a self-energy correction Σ(k) and a Zeeman coupling to an external magnetic field B = hẑ: We determine the self-energy Σ(k) by requiring that the off-diagonal expectation value Π ↑↓ (k) ≡ f † k↑ f k↓ calculated using the improved wave function matches that calculated using the effective Hamiltonian H eff (k).
To proceed, we consider an improved wave function similar to that in Eq. (14): where now we takeH = H ±z . The reason for this change is that the previously-considered correction due to H ±± gives exactly zero contribution to Π ↑↓ because it conserves the total spin S z modulo 2. The analogous contribution due to H ±z , however, does contribute. We expect that the energetically optimal value of the variational parameter isα ∼ J ±z /J 0 , where J 0 is on the order of the other exchange couplings (J ± and J z ).
Using the same perturbative expansion as above, the first-order form of Π ↑↓ (k) becomes If we represent the spin-spin interaction in momentum space with a momentum-dependent form factorγ where we arrive at the last line after conserving spin and momentum in the zeroth-order expectation values as well as usingγ(−k) =γ(k) andγ(0) = 0. Performing similar manipulations on Π (1) L (k) and combining the two contributions gives where n kσ = f † kσ f kσ is a number operator and a = | a µ | is the lattice constant. Importantly, Λ is real and independent of both µ and ± due to the sixfold symmetry S 6 . Furthermore, in the limit of T |h|, the integrand is only non-zero in an annulus of thickness ∼ h/(aJ 0 ) around the Fermi surface of radius ∼ 1/a, and the integral can then be estimated as Λ ∼ h/J 0 .
The real and imaginary parts of Σ(k) are plotted in Fig. 8. Note that the complex phase of Σ(k) ∝γ(k) winds by 4π around the Γ point. we have attempted to consolidate several previous results on this topic. We began by looking at the U (1) PSG wave functions derived in Ref. [18]. Instead of working with these wave functions phenomenologically, we go beyond their simple mean-field analysis and find quantitative estimates of the energies of these ansätze using variational Monte Carlo. We also use VMC to give a complete picture of magnetic order in our model. Our results improve on the classical magnetic phase diagrams presented in Refs. [15,16]. In those works, a phase transition between the 120 • and stripe phases is found in the nearest-neighbor model, and it is conjectured that large spin fluctuations may lead to the presence of a nonmagnetic phase. In our work, by building on the PSG ansätze, we also find a phase transition between the two magnetic phases in a similar parameter regime. We further find that second-neighbor interactions are necessary to create a spin liquid ground state and we identify the Dirac spin liquid as the lowest energy state. This confirms and extends earlier studies of the isotropic Heisenberg model [23].
The only other calculation of the full quantum phase diagram in this model was given by the DMRG analysis in Ref. [20]. Our phase diagram agrees with the DMRG analysis when second-neighbor interactions are included. The XXZ anisotropy and J ±± interactions both work to limit the spin liquid phase to a very small region of parameter space. However, we go beyond this and also include a third-neighbor interaction, which we believe gives a more complete picture on the behavior of the spin liquid phase. We find that even a very small third-neighbor interaction can greatly stabilize the spin liquid regime.

Relevance to Materials
This model has recently attracted much attention for its potential relevance to the material YbMgGaO 4 . Experiments find enticing evidence for a spinon Fermi surface state from thermodynamic and inelastic neutron scattering measurements [9,11]. Our work addressed the theoretical basis for such physics.
Our results support the claim of Ref. [20] that YbMgGaO 4 likely falls outside of the spin liquid phase in the presence of only first-and second-neighbor interactions. We found, however, that a very small third-neighbor interaction can greatly increase the size of the spin liquid phase and may appear quite naturally in the material. However, using the simple PSG picture, we always find that the Dirac spin liquid is energetically favored over the spinon Fermi surface state.
While the above results do not support the spinon Fermi surface state, we did find some effects which could tilt the balance in its favor. We saw that the spin-orbit interactions favor the spinon Fermi surface over the Dirac spin liquid state when we include effects beyond the simple projected mean-field wave functions. This leaves open the possibility that the spinon metal could be energetically favorable, perhaps assisted by other factors such as disorder or a small ring-exchange interaction.
If we assume that a spinon metal state does exist, interesting features emerge due to spinorbit coupling. We showed how the spin-orbit interactions could explain the existence of broad peaks at the M points in the spin structure factor and also predicted that measurements of the spin-polarized structure factors would display specific polarization-dependent peaks reflecting the anisotropic interactions. We also propose that the spin-orbit coupled spinon metal state may have a rather large thermal Hall conductivity which could be a very clear signature of spin liquid physics in such a system.

Future directions and implications
Looking forward, we anticipate a number of implications for the results and techniques developed in this work. For our spin-orbit coupled triangular systems, we showed that the restrictions imposed on the standard Gutzwiller-projected free fermion states by the PSG are quite severe for several of the U (1) QSL states. Consequently, they are unable to adapt to strongly anisotropic interactions, and this may open the door to competition from Z 2 QSL states in the case of such anisotropic models. In turn, this would be of considerable interest as the Gutzwiller-based approach almost always favors U (1) states in Heisenberg models. The possibility of inducing fully gapped topological QSLs should be explored in the future by VMC techniques.
We argued that the thermal Hall effect should be a key signature of itinerant spinon excitations in spin-orbit coupled systems. While we obtained such an effect for the U (1) spinon Fermi surface state on the triangular lattice, it was in fact suppressed by the PSGmandated vanishing of effective spin-orbit coupling on the fermionic spinons at the free-particle level. Ultimately, this suppression owes itself to the presence of inversion symmetry, which, in conjunction with time-reversal symmetry, act on the spinons analogously to the way they do on real electrons. As is well known, the combination of inversion and time reversal in that context imply an exact two-fold Kramers degeneracy of the full electronic band structure, and a similar effect occurs here. When inversion is absent, for example, when an electric field is present normal to a two-dimensional electron gas, spin splitting occurs. The Rashba spin-orbit coupling induced by such a field is known to induce a large anomalous Hall effect in that context [33]. This strongly suggests that one should look for an enhanced thermal Hall effect in two-dimensional magnetic materials in which the magnetic layer has an asymmetric environment. This criteria, along with the requirement of large spin-orbit coupling, should assist in a search for this phenomenon.
Our methodology offers a consistent and quantitative method to compare QSLs and ordered phases for anisotropic magnetic Hamiltonians. This should have broad applicability to other materials such as the Kitaev compounds α-RuCl 3 , Na 2 IrO 3 , and Li 2 IrO 3 in all its structural variations, and to three-dimensional systems like rare earth pyrochlores and spinels. The ability of VMC-based methods to tackle large systems is a unique numerical advantage. We expect many insights from such studies in the future.