Instability of the $U$(1) spin liquid with a large spinon Fermi surface in the Heisenberg-ring exchange model on the triangular lattice

The Fermionic RVB state can be expanded in the Ising basis as

in which

denotes an Ising basis written in terms of the Fermion Fock state

is the corresponding wave function amplitude. Here and in the following, we will use $i_{k}$ and $j_{k}$ to denote the locations of the $k$-th up and the $k$-th down spin in the Ising basis $|\sigma_{1},....,\sigma_{N}\rangle$, rather than the two end points of the $k$-th valence bond. $\mathbf{A}$ is a $\frac{N}{2}\times\frac{N}{2}$ matrix with its matrix element given by $a(i_{k},j_{k^{\prime}})$ at its $k$-th row and $k^{\prime}$-th column.

In practice, we can treat the RVB amplitude $a(i_{k},j_{k})$ as variational parameter directly. The variational state constructed in this way will be referred to as the general RVB state(gRVB) in the following discussion. The number of variational parameters in the gRVB state increases very rapidly with the system size. As we will see in the next section, such an unfavorable feature of the gRVB state is compensated partly by the fact that the calculation of the energy gradient in the gRVB state is rather cheap.

The Fermionic RVB state can also be generated by Gutzwiller projection of mean field ground state of the following Bogoliubov-de Gennes Hamiltonian

Here the condition $\Delta_{i,j}=\Delta_{j,i}$ is imposed to enforce spin rotational symmetry of the variational ground state. In general, both $\chi_{i,j}$ and $\Delta_{i,j}$ can be chosen complex. In our work, $\chi_{i,j}$ and $\Delta_{i,j}$ will be chosen real and be restricted to the nearest neighboring bonds.They are otherwise free of any assumption.

$H_{MF}$ is usually referred to as a mean field ansatz or a variational ansatz of the Fermionic RVB state. To generate $|f-RVB\rangle$, we rewrite $H_{MF}$ in the following form

in which

Here $\bm{\chi}$ and $\bm{\Delta}$ are $N\times N$ matrices with $\chi_{i,j}$ and $\Delta_{i,j}$ as their matrix elements. $H_{MF}$ can be diagonalized by the following unitary transformation

in which $\bm{u}$ and $\bm{v}$ are $N\times N$ matrices. The diagonalized Hamiltonian takes the form of

in which $\bm{E}$ is a $N\times N$ diagonal matrix with positive definite diagonal matrix elements. When $\bm{\Delta}\neq 0$, the RVB amplitude $a(i_{k},j_{k})$ of the corresponding Fermionic RVB state is given by

When $\bm{\Delta}=0$, the RVB amplitude is given by

in which $\varphi_{m}(j_{k})$ is the eigenvector of the matrix $\bm{\chi}$ with eigenvalue $E_{m}$. Here $\mu$ denotes the chemical potential of the Fermionic spinon at half filling.

In practical variational calculation, it is often convenient to adopt a particle-hole transformation on the down-spin Fermions, which is given by

The mean field ground state of $H_{MF}$ is then constructed by filling up the lowest $N$ eigenvectors of $\mathbf{M}$. More specifically, we can expand the RVB state as

in which

is the reference state with all sites occupied by down spin Fermions(or the vacuum of $\tilde{f}$ operators). Here $i_{1},...,i_{N/2}$ denote the locations of the $N/2$ up spin Fermions. Note that they are also the locations of the $N/2$ holes of the down spin Fermion. The wave function in the particle-hole transformed picture then takes the form of

in which $\bm{\Phi}$ is a $N\times N$ matrix of the form

Here $\phi_{n}(i)$ denotes the $i$-th component of the $n$-th eigenvector of the matrix $\mathbf{M}$ with negative eigenvalue.

In this study, we use either the general Fermionic RVB state or RVB state generated from mean field ansatz to describe the ground state of the $J_{1}-J_{4}$ model. For the general RVB state, there will be $\frac{N(N-1)}{2}$ variational parameters to be optimized on a finite cluster with $N$ sites. For RVB state generated from mean field ansatz, we can choose either a $U(1)$ ansatz, in which $\bm{\Delta}=0$, or a $Z_{2}$ ansatz, in which both $\bm{\chi}$ and $\bm{\Delta}$ are nonzero. In the $U(1)$ ansatz, there are $3N$ variational parameters to be optimized, which are the hopping amplitude $\chi_{i,j}$ on the $3N$ nearest neighboring bonds. In the $Z_{2}$ ansatz, there are $6N+1$ variational parameters to be optimized, which are the hopping amplitude $\chi_{i,j}$ and the pairing amplitude $\Delta_{i,j}$ on the $3N$ nearest neighboring bonds and the chemical potential setting the Fermi level of the spinon. The gRVB state represents the most general form of a Fermionic RVB state and it may not be generated by any short-ranged mean field ansatz. Correspondingly, it contains a much larger number of variational parameters. The optimization of large number of variational parameters calls for very efficient optimization algorithm, which we will now introduce.

The steepest descent(SD) algorithm is the simplest optimization algorithm. It corresponds to setting the Hessian matrix proportional to the identity matrix. In the SD algorithm, the variational parameters are updated as follows

in which $\delta$ is the step length. The step length is usually adjusted by trial and error. A more intelligent way to adjust the step length is the following self-learning trick, in which the step length is updated according to the change in the direction of the gradient as follows

Here $\eta\in(0,1)$ is an acceleration factor, $\nabla E$ and $\nabla E^{\prime}$ are the gradients of the energy in the current and the previous step. Suitable choice of $\eta$ can accelerate significantly the optimization procedure at the initial stage, when the energy gradient is large.

The SD algorithm usually works very well at the initial stage of the optimization procedure. However, it loses efficiency when the optimization procedure encounters long and narrow valley with flat bottom in the energy landscape. In such a situation, there will be large fluctuation in the eigenvalues of the Hessian matrix. Approximating the Hessian matrix with an identity matrix is then clearly not a wise choice.

The stochastic reconfiguration algorithm is a widely used variational optimization algorithm. It mimics the effect of the Hessian matrix with a positive-definite Hermitian matrix $\mathbf{S}$ generated from the metric of the variational state in the variational spaceSeiji . More specifically, $\mathbf{S}$ is given by

It is easy to show that $\mathbf{S}$ is just the Hessian matrix of the following quantity with respect to the change of variational parameters

Here we assume that $|\Psi\rangle$ is fixed and $|\Psi^{\prime}\rangle$ is varying. $\Delta^{\mathrm{SR}}$ defined in this way can be interpreted as the distance between $|\Psi\rangle$ and $|\Psi^{\prime}\rangle$ in the Hilbert space.

In the SR algorithm, the variational parameters are updated as follows

in which $\delta$ is the step length. The self-learning acceleration trick is also applicable in the SR method.

The introduction of the $\mathbf{S}$ matrix in the SR method amounts to replace the naive distance in the Euclidean space of variational parameters with the distance in the Hilbert space. Such a regulation procedure can be very helpful when some variational parameters are nearly redundant. However, since $\mathbf{S}$ only depends on the variational state but not on the Hamiltonian, it can not approximate the effect of the true Hessian matrix correctly in certain situations. In practice, the SR method may still suffer slow convergence or even run away from true minimum. A better approximation of the effect of the Hessian matrix is needed.

The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method is a quasi-Newton method. It generates an iterative approximation for the (inverse) of the Hessian matrixNocedal from the gradient. The approximate (inverse)Hessian matrix is updated as follows

in which $k=0,1,2.....$

are the difference between successive variational parameters and energy gradient. $\bm{\alpha}_{k=0}$ is the initial guess of the variational parameters and $\nabla E_{k=0}$ is the energy gradient at the starting point. The Hessian matrix is initially set to be the identity matrix, namely $\mathbf{B}_{k=0}=\mathbf{I}$.

Using such an iterative approximation on the Hessian matrix, the variational parameters are updated as follows

Here $\delta$ is the step length. In principle the step length should be determined by a linear search in the direction of $\mathbf{B}_{k}\nabla E_{k}$. Such a linear search can be accomplished in principle in the variational Monte Carlo simulation by reweighting in the searching direction. However, to reduce computational expense, we choose a fixed step length by trial and error.

Another simpler method to go beyond the steepest descent method is the conjugate gradient(CG) methodNocedal . It corrects the searching direction as follows

in which $k=0,1,2,....$

or

Initially we set $\beta_{0}=0$.

The variational parameters are updated in the conjugate gradient method as follows

in which the step length $\delta$ should be determined by linear search in the direction of $\mathbf{d}_{k}$. In practice, to reduce computational expense we choose a fixed step length by trial and error.

We note that for both the BFGS and the CG method, the initial step of the optimization is just the SD update. As a result of the finite accuracy in the computation of the energy and the gradient in variational Monte Carlo simulation, we cut off the BFGS and the CG iteration at a finite depth. We find empirically that 10 steps of BFGS or CG iteration has the best balance between numerical efficiency and numerical stability, after which we restart the iteration by setting $k=0$.

Fig.1 compares the performance of the four algorithms in a typical situation. Here we optimize the variational energy of the $J_{1}-J_{4}$ model at $J_{4}=0.065$ with the $U(1)$ mean field ansatz. We start from the same initial guess and use the same normalized step length. We find that the BFGS method has the best numerical efficiency and stability among the four algorithms. The CG algorithm exhibits a similar numerical efficiency as the BFGS algorithm but with a significantly larger fluctuation. Both the steepest descent and the stochastic reconfiguration method fail to escape from the the local minimum around $E\approx-0.4976$ within 2000 optimization steps. In the following, we will mainly use the BFGS method. However, we note that the conjugate gradient method has the advantage that it does not need to store the approximate Hessian matrix, which is huge when the number of variational parameters is large.

Fig.4 plots the variational energies computed from the three kinds of RVB state as a function of $J_{4}$ for $0\leq J_{4}\leq 0.045$. The variational energies of all these three kinds of RVB state are very close to each other and are within statistical error perfect linear functions of $J_{4}$. Such a nearly zero curvature behavior in the variational energy indicates that the variational ground state is almost independent of $J_{4}$. In the case of the $U(1)$ RVB state, the optimized variational state reduces to the well known $\pi$-flux phase on the triangular lattice. This can be explicitly seen from the optimized variational parameters of the $U(1)$ ansatz, which indeed encloses a gauge flux of $\pi$ around every elementary rhombi of the triangular lattice. We find that the small difference in the optimized variational energies between the $Z_{2}$ and the $U(1)$ RVB state can be attributed to finite size effect. On a $24\times 24$ cluster, the optimized variational energy of the two states become indistinguishable within the statistical error of the variational Monte Carlo simulation. The equivalence between the optimized $U(1)$ and $Z_{2}$ RVB state in this low $J_{4}$ regime can also be seen from the static spin structure factor $S(\mathbf{q})$, which is defined as

Here we define the components of the wave vector as follows

in which $\mathbf{b}_{1,2}$ are the two reciprocal vectors of the triangular lattice. In Fig.5, we present the static spin structure factor of the optimized $U(1)$ and $Z_{2}$ RVB state at $J_{4}=0$. We find that $S(\mathbf{q})$ is within statistical error independent of $J_{4}$ in the whole $0\leq J_{4}\leq 0.045$ regime and features pronounced peaks at the wave vector corresponding to the 120 degree ordered phase, namely $\mathbf{q}=(\frac{2\pi}{3},\frac{4\pi}{3})$ and $\mathbf{q}=(\frac{4\pi}{3},\frac{2\pi}{3})$.

The energy difference between the general RVB state and the $U(1)$ or the $Z_{2}$ RVB state is more tricky. We find that the optimized gRVB state can not be generated by any short-ranged mean field ansatz. However, the static spin structure factor of the optimized gRVB state is essentially indistinguishable from that of the optimized $U(1)$ or $Z_{2}$ RVB state(see Fig5c and 5d). It is currently impossible to conduct the optimization of the gRVB state on a cluster significantly larger than the $12\times 12$ cluster. For example, on a $24\times 24$ cluster, the number of variational parameters in the gRVB state becomes $\frac{N(N-1)}{2}=165600$, which is too large to be reliably optimized.

In all, we find that the optimized RVB state in the whole $0\leq J_{4}\leq 0.045$ regime represents a state essentially equivalent to the $\pi$-flux phase on the triangular lattice, which is the closest counterpart of the 120 degree ordered phase in the space of Fermionic RVB state.

It is generally believed that in the large $J_{4}$ regime the SFS state is the most favorable variational ground state. Indeed, we find that the SFS state is locally extremely stable in this regime. The SFS state is generated by the following mean field ansatz

As a result of the translational and rotational symmetry of the mean field ansatz, the energy gradient must be also translational and rotational symmetric. The energy gradient in the SFS phase thus must vanish since the Gutzwiller projected state is invariant under a uniform rescaling of all variational parameters. In other words, the SFS phase is an exact saddle point in the space of Fermionic RVB state. To illustrate the local stability of the SFS state, we have performed variational optimization at $J_{4}=0.15$ starting from the following initial guess of the mean field ansatz

Here $r$ is a random number distributed uniformly in the range of $r\in(0,1)$, $\eta$ is a constant measuring the distance of the initial guess from the SFS ansatz. As is shown in Fig.6, even if we choose $\eta$ as large as $\eta=0.75$, the variational energy still converges to a value very close to that of the SFS state, which is $E\approx-0.576$ at $J_{4}=0.15$. The corresponding variational parameters also converge to that of the SFS ansatz.

However, the SFS state is only locally stable. To illustrate this point, we have performed variational optimization starting from fully random initial guess of the variational parameters. The optimization is done for $J_{4}=0.15$ on a $12\times 12$ cluster with periodic boundary condition. A preliminary trial shows that the optimized variational ground state exhibits an approximate $4\times 6$ modulation in its local spin correlation pattern. This is illustrated in Fig.7. We note that the modulation in the local spin correlation is very strong, ranging from almost pure spin singlet correlation to pure spin triplet correlation between nearest neighboring spins. The translational symmetry is seriously broken.

We then refine the optimization by assuming a $4\times 6$ periodicity in the variational parameters. We find that for both the $U(1)$ and the $Z_{2}$ mean field ansatz, the optimized variational energy converge to values much lower than that of the SFS state. Fig.8 and Fig.9 illustrate the convergence of the variational energies for randomly chosen initial guess of the $U(1)$ and the $Z_{2}$ mean field ansatz. For $U(1)$ mean field ansatz, we find that about $1/3$ of the variational optimization trials converge to the lowest variational energy of $E\approx-0.603$. For $Z_{2}$ mean field ansatz, about $1/2$ of the variational optimization trials converge to the lowest variational energy of $E\approx-0.606$. The optimized variational energy for both types of RVB states are about $5\%$ lower than the energy of the SFS state, which is $E\approx-0.576$ at $J_{4}=0.15$.

On the $12\times 12$ cluster studied here, the most general periodicity of the mean field ansatz that is compatible with the periodic boundary condition of the cluster is $C_{1}\times C_{2}$, in which $C_{1,2}=1,2,3,4,6,12$ denotes the period in the $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ direction. We have performed refined variational optimization assuming each of such periodicity. We find that the lowest variational energy is always obtained with the $4\times 6$ periodicity at $J_{4}=0.15$.

We find that the optimized variational energy as a function of $J_{4}$ exhibits a very small curvature in the large $J_{4}$ regime(see Fig.10). This implies that the variational ground state is almost independent of $J_{4}$ in this part of the phase diagram also. A naive linear extrapolation of the variational energy implies that the SFS state can not be the true ground state of the $J_{1}-J_{4}$ model for $J_{4}\leq 0.5$. To check if the SFS state can be stabilized at still larger $J_{4}$, we have performed variational optimization for $J_{4}=10$, a value that is too large to be realistic for real materials. We find that a symmetry breaking phase with the same $4\times 6$ periodicity in its spin correlation pattern is still significantly lower in energy than the SFS state. More specifically, we find that the optimized variational energy is $E\approx-16.09$ at $J_{4}=10$. This is again about $5\%$ lower than the energy of the SFS state, which is $E\approx-15.33$ at $J_{4}=10$. Fig.11 illustrates the local spin correlation pattern in the optimized state, which is similar to that at $J_{4}=0.15$. We thus conclude that the SFS state is never the best RVB state for the $J_{1}-J_{4}$ model and that the $4\times 6$ state is the best variational ground state in the whole range of $J_{4}\geq 0.09$. Such a state breaks both the translational and the point group symmetry of the model and exhibits $4\times 6$ modulation in its local spin correlation pattern.

The optimized variational ground state in the intermediate regime of $0.045\leq J_{4}\leq 0.09$ is characterized by a zigzag modulation in its local spin correlation pattern. As is illustrated in Fig.12, the zigzag pattern manifests itself most evidently in the antiferromagnetic correlated backbones running through the $2\mathbf{a}_{1}-\mathbf{a}_{2}$ or equivalent directions of the triangular lattice. The translational symmetry and the rotational symmetry are thus spontaneously broken. Looking more closely, we find that the zigzag pattern in the $0.045\leq J_{4}\leq 0.055$ regime differs from that in the $0.055\leq J_{4}\leq 0.09$ regime by an additional two-fold translational symmetry breaking. Fig.13 shows the optimized variational energy for the zigzag phase. The energy difference between the SFS state and the zigzag state is even more significant than that in the $4\times 6$ phase. A closer inspection also shows that the energy difference between the $U(1)$ and the $Z_{2}$ state is more significant in the zigzag regime than that in the $\pi$-flux and the large $J_{4}$ regime.

Instead of the zigzag phase found here, several translational invariant phases have been proposed in the intermediate $J_{4}$ regime in the literature. In Ref.[Grover, ], a nodal $d$-wave $Z_{2}$ state with nematic spinon dispersion is proposed to be the variational ground state in the intermediate regime. It is later found that another $Z_{2}$ state with $d+id$-wave spinon pairing and quadratic band touching(QBT) in the spinon dispersion is more favorable than the nodal d-wave state in a tiny range of $J_{4}$Xu . The QBT state hosts gapless spinon excitation with a finite density of state and a gapped gauge fluctuation spectrum as a result of the $d+id$-wave spinon pairing. Such an excitation characteristics is argued to be helpful to resolve the puzzle related to the anomalous gauge fluctuation correction to the specific heat in the SFS state. While we think that the anomalous gauge fluctuation correction in the SFS state is only a theoretical artifact of the conventional gauge theory argumentLi2 , it is nevertheless interesting to compare the variational energies of these novel states to our variational result. In Fig.14, we plot the variational energies of the major candidates of the variational ground state of the $J_{1}-J_{4}$ model in the intermediate $J_{4}$ regime. From this figure we see that the energy advantage of the QBT state over the nodal $d$-wave state is meaningless when compared to the huge energy difference between these states and the zigzag state we find from variational optimization. This establishes firmly the zigzag nature of the variational ground state in the intermediate regime.

A zigzag phase has also been reported in a similar parameter regime in a recent DMRG study of the $J_{1}-J_{2}-J_{4}$ modelMoore , in which $J_{2}$ denotes the exchange coupling between next nearest neighboring spins. Different from what we found here, the zigzag phase reported in Ref.[Moore, ] exhibits magnetic long range order. The Fermionic RVB state we adopted in this study is not expected to describe accurately such magnetic long range ordered phases. However, as we have seen in the case of the $\pi$-flux phase for $0\leq J_{4}\leq 0.045$, the Fermionic RVB state can nevertheless reproduce correctly the qualitative feature of the spin structure factor of the magnetic ordered phase. We find that this is also the case in the zigzag phase. In Fig.15 we plot the spin structure factor of the model at $J_{4}=0.06$ calculated from the optimized gRVB statenote1 . The spin structure factor is characterized by the prominent peaks at $\mathbf{q}=(\pm\frac{\pi}{2},\pi)$ and the weaker peak at $\mathbf{q}=(\pi,0)$. These are exactly the positions of the magnetic Bragg peaks in the classical zigzag phase(illustrated in Fig15c).

With these considerations in mind, it is better to interpret the $\pi$-flux phase in the small $J_{4}$ regime and the zigzag phase in the intermediate $J_{4}$ regime both as the closest approximation of the corresponding magnetic ordered phases, namely the 120 degree ordered phase and the zigzag ordered phase. Thus, the transition at $J_{4}=0.045$ should be better understood as a first order transition between two magnetic ordered phases, rather than the transition between a spin liquid phase and a valence bond solid phase. According to Ref.[Xu, ], the 120 degree ordered phase becomes degenerate with the nodal d-wave state at $J_{4}=0.0525$. Since the energy of the zigzag phase found here is lower than that of the nodal d-wave state at $J_{4}=0.0525$, we expect that the transition between the 120 degree ordered phase and the zigzag phase to occur at a smaller value of $J_{4}$ than $0.0525$. This is indeed the case in our calculation.

To provide further support to the variational results, we have performed exact diagonalization(ED) calculation on a $6\times 6$ cluster with periodic boundary condition. Here we focus on the identity representation. Using translational and point group symmetry, together with the spin rotational symmetry along the $z$-axis, we can reduce the number of symmetrized basis to 31,554,903. The Hamiltonian matrix is then diagonalized by the Lanczos method to obtain its lowest few eigenvalues. Fig.16 shows the ground state energy as a function of $J_{4}$ obtained from the ED calculation. Two features of the ground state energy curve are of particular interest to us. First, the curvature of the ground state energy is very small in both the small $J_{4}$ and the large $J_{4}$ regime, a trend in good agreement with the variational result presented above. Second, in the intermediate $J_{4}$ regime, level crossings between the ground state and the first excited state are observed, implying potential first order transition in the thermodynamic limit. Such transitions may be closely related to the appearance of the zigzag phase we obtained in the variational study. Clearly, to make more solid conclusion on this issue, more systematic ED calculation is needed.

We note that the finite size effect on a $6\times 6$ cluster can be very significant, especially when the considered state is gapless. For example, if we perform variational optimization on a $6\times 6$ cluster with periodic boundary condition, what we would obtain at $J_{4}=0.15$ is a symmetry broken state with a $\sqrt{3}\times\sqrt{3}$ modulation in its local spin correlation pattern( as is illustrated in Fig.17), rather the $4\times 6$ state we find on the $12\times 12$ cluster. The modulation of the local spin correlation in the $\sqrt{3}\times\sqrt{3}$ state is found to be much weaker than that in the $4\times 6$ state. Such a result is partially expected since the $4\times 6$ state is incompatible with the periodic boundary condition of the $6\times 6$ cluster. However, when we extend the optimized ansatz of such a $\sqrt{3}\times\sqrt{3}$ state to a $12\times 12$ cluster, we find that the variational energy is very bad. On the other hand, when we extend the optimized ansatz of the $4\times 6$ state on the $12\times 12$ cluster to a $24\times 24$ cluster, the variational energy is essentially unchanged. Fig.18 illustrates such contrasting behavior of the $\sqrt{3}\times\sqrt{3}$ state and the $4\times 6$ state. Thus, the dominance of the nearly uniform $\sqrt{3}\times\sqrt{3}$ state on the $6\times 6$ cluster is purely a finite size effect. This indicates that a sufficiently large cluster is needed to draw conclusion on the phase diagram of model.