Superconductivity, pair density wave, and Neel order in cuprates

TL;DR

This paper explores the coexistence of superconductivity, pair density waves, and Neel order in underdoped cuprates, revealing different favored states at half filling and finite doping through variational quantum Monte Carlo simulations.

Contribution

It introduces a comprehensive analysis of competing orders in cuprates, identifying the energetically favored states and their evolution with doping levels.

Findings

At half filling, the $pS_{\bf Q}+pT_0$ state is more stable.

At finite doping, the $dS_0+dT_{\bf Q}$ state becomes more favorable.

Results suggest proximity to a $\mathrm{Z_2}$ quantum spin liquid with a spinon Fermi surface.

Abstract

We investigate in underdoped cuprates possible coexistence of the superconducting (SC) order at zero momentum and pair density wave (PDW) at momentum ${\bf Q}=(\pi, \pi)$ in the presence of a Neel order. By symmetry, the $d$-wave uniform singlet pairing $dS_0$ can coexist with the $d$-wave triplet PDW $dT_{\bf Q}$, and the $p$-wave singlet PDW $pS_{\bf Q}$ can coexist with the $p$-wave uniform triplet $pT_0$. At half filling, we find the novel $pS_{\bf Q}+pT_0$ state is energetically more favorable than the $dS_0+dT_{\bf Q}$ state. At finite doping, however, the $dS_0+dT_{\bf Q}$ state is more favorable. In both types of states, the variational triplet parameters, $dT_{\bf Q}$ and $pT_0$, are of secondary significance. Our results point to a fully symmetric $\mathrm{Z_2}$ quantum spin liquid with spinon Fermi surface in proximity to the Neel order at zero doping, and to intertwined $d$-wave triplet PDW fluctuations and spin moment fluctuations along with the dominant $d$-wave singlet SC at finite doping. The results are obtained by variational quantum Monte Carlo simulations.

Equations

$$H=$$

$$H_v=$$

$$U_{ij}=-t_{1v}{\tau }_3+{\Delta }_{ij}{\tau }^{+}+{\Delta }_{ji}^{\ast }{\tau }^{-}.$$

$${\Delta }_{ij}=(b_x^2-b_y^2)(S_0+T_{\mathbf{Q}}{\eta }_i)+(b_x+i{b}_y)(S_{\mathbf{Q}}{\eta }_i+T_0),$$

$$x=\{\mu ,{\mu }_z,t_{2v},S_0,T_{\mathbf{Q}}\},d-\mathrm{wave}\mathrm{case},$$

$$x=\{\mu ,{\mu }_z,t_{2v},S_{\mathbf{Q}},T_0\},p-\mathrm{wave}\mathrm{case}.$$

$$|G⟩=\frac{P|{\psi }_0⟩}{\sqrt{⟨{\psi }_0|P|{\psi }_0⟩}},$$

$$E=\frac{1}{N}⟨G|H|G⟩,m=\frac{1}{N}\sum _i{\eta }_i⟨G|S_i^z|G⟩.$$

$$U_{ij}\to V_i{U}_{ij}{V}_j^{†},V_i=e^{i\mathbf{Q}\cdot {\mathbf{r}}_i{\tau }_3/2},$$

$$U_{i,i+\hat{x}}=-i{\tau }_0+S_{\mathbf{Q}}{\tau }_2,U_{i,i+\hat{y}}=-i{\tau }_0+S_{\mathbf{Q}}{\tau }_1,$$

$$W_{\mathcal{P}(1234)}=\mathcal{P}(U_1{U}_2{U}_3{U}_4),$$

$$\displaystyle\begin{array}[]{ccccc}g=&P_{x},&P_{y},&P_{d},&T.\\ V_{g}=&(-)^{x_{i}}i\tau_{1},&(-)^{y_{i}}i\tau_{2},&\frac{i}{\sqrt{2}}(\tau_{1}+\tau_{2}),&(-)^{x_{i}+y_{i}}.\end{array}$$

$$H_v={\psi }_{\mathbf{k}}^{†}[{\epsilon }_{\mathbf{k}}{\tau }_0+2S_{\mathbf{Q}}(\cos k_x{\tau }_2+\cos k_y{\tau }_1)]{\psi }_{\mathbf{k}}.$$

$$E_{\mathbf{k}}={\epsilon }_{\mathbf{k}}\pm 2|S_{\mathbf{Q}}|\sqrt{{\cos }^2{k}_x+{\cos }^2{k}_y}.$$

$$H_v=\sum _{ij}{\psi }_i^{†}{M}_{ij}{\psi }_j,$$

$$M{\phi }_n={ϵ}_n{\phi }_n,$$

$$U=(u_1,u_2,\cdots ,u_N),V=(v_1,v_2,\cdots ,v_N),$$

$$A=U{V}^{-1}.$$

$$|{\psi }_0⟩={\left(\sum _{ij}{c}_{i\uparrow }^{†}{A}_{ij}{c}_{j\downarrow }^{†}\right)}^{N_e/2}|0⟩,$$

$$d{x}^{\mu }|\mu ⟩=d{x}^a|a⟩,d{x}^a=⟨a|\mu ⟩d{x}^{\mu }\equiv e_{\mu }^{a}d{x}^{\mu }.$$

$$d{x}^{\mu }d{x}^{\nu }⟨\mu |\nu ⟩\equiv d{x}^{\mu }d{x}^{\nu }{g}_{\mu \nu }=d{x}^{a\ast }d{x}^b{\delta }_{ab}.$$

$${\partial }_{a}E=e_a^{\mu }{\partial }_{\mu }E,e_a^{\mu }{e}_{\nu }^a={\delta }_{\nu }^{\mu }.$$

$$d{x}^a=-{\partial }_{a}Edt.$$

$$e_{\mu }^{a}d{x}^{\mu }=-e_a^{\mu }{\partial }_{\mu }Edt.$$

$$g_{\sigma \mu }d{x}^{\mu }=-{\partial }_{\sigma }Edt,\to d{x}^{\mu }=-g^{\mu \nu }{\partial }_{\nu }Edt,$$

$$d{x}^{a\ast }=⟨\mu |a⟩d{x}^{\mu }.$$

$$d{x}^a=-\frac{\partial E}{\partial x^{a\ast }}dt.$$

$$⟨a|\mu ⟩d{x}^{\mu }=-\frac{\partial x^{\mu }}{\partial x^{a\ast }}{\partial }_{\mu }Edt.$$

$$g_{\nu \mu }d{x}^{\mu }=-\frac{\partial x^{a\ast }}{\partial x^{\nu }}\frac{\partial x^{\mu }}{\partial x^{a\ast }}{\partial }_{\mu }Edt=-{\partial }_{\nu }Edt,$$

$$⟨\nu |a⟩⟨a|\mu ⟩=⟨\nu |\mu ⟩=g_{\nu \mu },⟨\nu |a⟩=\frac{\partial x^{a\ast }}{\partial x^{\nu }}.$$

$$g^{\mu \nu }\to \sum _k^{{}^{\prime }}{\varphi }_k(\mu )\frac{1}{{\Lambda }_k}{\varphi }_k^{\ast }(\nu ),$$

Full text

I introduction

The investigation of the mechanism of high-temperature superconductivity in cuprates remains to be an exciting topic. One of the interesting proposals is Anderson’s resonating valence bond (RVB) state. Anderson87-RVB A suitable Hamiltonian describing such a system is the one-band $t$-$J$ model. ZhangRice88-tJ In this model, while the parent compound at half filling is automatically a Mott insulator for the charge degrees of freedom, the spin sector is much more intriguing. The RVB state may be viewed as a linear combination of configurations of the covering of spin singlets, a quantum spin liquid (QSL) with fractional spinon excitations. Chemical doping introduces mobile holes, leaving room for spin singlets to relocate and hence making the system a superconductor immediately.Lee06-review Initially, an $s$-wave RVB state is proposed in view of the experimental robustness of superconductivity against impurity scattering,swave but it is found later that the $d$-wave RVB state is energetically better,Gros88-VQMC ; dwave and can actually be robust against impurities in doped Mott insulators because of charge renormalization. RMFT

It should be pointed out, however, that in the undoped limit, there is a local charge-SU(2) symmetry following from the charge neutrality.su2 This symmetry relates various forms of RVB states. For example, the RVB state with a $\pi$-flux in the spinon hopping around a plaquette may be mapped to a state with $d$-wave pairing, etc. Such states are said to be gauge equivalent and describe the same spin liquid upon projection to the physical Hilbert space. The projective symmetry group (PSG) has been developed Wen02-PSG ; Wen2004 to classify all possible and physically distinct spin liquids. The theory is built on the assumption that the spins are in a quantum disordered state. In a pure two-dimensional (2D) model, this might be a reasonable assumption since the continuous spin-SU(2) symmetry cannot be broken spontaneously at finite temperatures. However, one can still ask whether the moments could order at zero temperature, or in the ground state. Including inter-layer coupling can extend such a ground state order into finite temperatures. In fact, the Neel order is observed experimentally in the parent undoped compounds of cuprates. When the spin ordering is taken into account, the notion of spin liquid, with fractional spinon excitations, seems to be less well-defined, as spin-$\frac{1}{2}$ spinons would have been confined to form spin-1 magnons. However, inelastic neutron scattering experiments show that the spin excitations away from the Neel vector are broadened significantly, well beyond the linear spin-wave description. Piazza15-neutron A recent interesting proposal is that even in the presence of Neel order, the spinons may be deconfined in a partial region of the Brillouin zone, pdeconf1 ; pdeconf2 ; Yu18-CPT although the same phenomenon could also be understood in terms of magnon-magnon scattering.magnon1 ; magnon2 ; Yu18-CPT It is therefore interesting to consider spinon states on the background of the Neel order.

The charge-SU(2) symmetry is broken at finite doping. The Neel order is weakened but persists at small finite doping. The effect of this order on superconductivity is also an intriguing topic. In this case, the spin-SU(2) symmetry is broken down to O(2). As a result, there is no longer sharp distinction between spin-singlet and spin-triplet Cooper pairings, and there is room for coexistence, Lu14-topo ; Gupta16-meanfield although the relative weight is not dictated by symmetry. There is a residual point group $C_{4}$ (with respect to a site), which leaves the spin moment invariant in the Neel state. (In 3D, inversion is also a symmetry.) This symmetry dictates what kinds of singlet and triplet could coexist. There are three irreducible representations for the $C_{4}$ group, namely, the 1D $A$ and $B$ representations, and the 2D $E$ representation. Since the completely symmetric $A$ representation is not a favorable pairing symmetry we will ignore it henceforth. The $B$ representation transforms as $d$-wave. The doubly degenerate $E$ representation transforms as $p$-wave. Therefore, the singlet and triplet should transform under $C_{4}$ identically either as $d$-wave or $p$-wave. These possibilities are illustrated in Fig.1. In the first case, the $d$-wave singlet Cooper pair at momentum ${\bf q}=0$ (a), can coexist with a $d$-wave triplet at momentum ${\bf Q}\equiv(\pi,\pi)$ (b). The latter is a pair density wave (PDW), namely, the Cooper pairing at the center-of-mass momentum of ${\bf Q}$. In the second case, the $p+ip^{\prime}$-wave singlet PDW at momentum ${\bf q}={\bf Q}$ (c) can coexist with the same $p+ip^{\prime}$-wave but triplet SC state (d). The four types of states are denoted in a self-explaining manner as $dS_{0}$, $dT_{\bf Q}$, $pS_{\bf Q}$ and $pT_{0}$. Interestingly, in the PDW states, an electron at momentum ${\bf k}_{1}={\bf k}+{\bf Q}$ pairs up with another one at ${\bf k}_{2}=-{\bf k}$. When ${\bf k}$ is on the so-called umklapp surface (US), so will be both ${\bf k}_{1}$ and ${\bf k}_{2}$, related by mirror symmetry. The scattering of such a Cooper pair across the US was argued to be the key mechanism that could generate not only the single-particle gap but also two-particle gap near the antinodes and hence the pseudogap. maurice On the other hand, the chiral $pT_{0}$ state was recently proposed Lu14-topo to be present on the background of the Neel order, in an effort to explain in the underdoped regime the opening of a mini gap in the quasiparticle excitations along the otherwise gapless (nodal) direction of the $dS_{0}$ state. Shen04-ARPES ; Tanaka06-ARPES ; Vishik12-ARPES ; Razzoli13-ARPES ; Peng13-ARPES Remarkably, if this were the case, the cuprate would be topologically nontrivial because of the chiral $p+ip^{\prime}$ pairing. ReadGreen00-topo ; QiZhang11-review

We are therefore motivated to investigate possible coexistence of superconducting (SC) order at zero momentum and PDW at momentum ${\bf Q}$ in the presence of a Neel order. We use the variational quantum Monte Carlo (VQMC) to treat the strong correlation effects. Our main results are as follows. At half filling, we find the novel $pS_{\bf Q}+pT_{0}$ state is energetically more favorable than the $dS_{0}+dT_{\bf Q}$ state (including the conventional $dS_{0}$ state). We find that on its own the $pS_{\bf Q}$ state is a fully symmetric $\mathrm{Z_{2}}$ QSL with spinon Fermi surfaces, but it is unstable toward Neel ordering. At finite low doping, however, the $dS_{0}+dT_{\bf Q}$ state is more favorable. In both types of states, the variational triplet parameters, $dT_{\bf Q}$ and $pT_{0}$, are of secondary significance since the ground state energy is degenerate with respect to varying levels of such variables (to some extent) upon optimization of the other parameters (including the chemical potential and hopping integrals). However, they do enhance the average Neel moment. Our results point to a novel QSL in proximity to the Neel order at zero doping, and to intertwined $d$-wave triplet PDW fluctuations and spin moment fluctuations along with the dominant $d$-wave singlet SC at finite doping.

In the rest of this paper, we specify the model and method in Sec.II, discuss the results in Sec.III, and we draw conclusions and make perspective remarks in Sec.IV.

II model and method

We begin with the $t$-$J$ model on the square lattice, described by the Hamiltonian

[TABLE]

Here $\sigma$ is the spin polarization, $N_{1,2}$ denote the first- and second-neighbor bonds with hopping integrals $t_{1}=t$ and $t_{2}=t^{\prime}$, respectively; $c_{i\sigma}$ is the electron annihilation operator, $J$ is the Heisenberg spin exchange, ${\bf S}_{i}$ is the local spin, and $n_{i}$ is the local density. The operator $P=\Pi_{i}(1-n_{i\uparrow}n_{i\downarrow})$ projects away any double occupancy. As typical parameters for cuprates, we take $(t,t^{\prime},J)=(0.4,-0.12,0.13)$ eV. Lee06-review At half filling, the charge degrees of freedom are frozen, and the model reduces to the Heisenberg model. Upon hole doping, the doped holes may move without causing double occupancy, leading to metallicity and superconductivity.

The Hamiltonian includes infinitely strong correlations, due to the fact that no double-occupancy is allowed. In this work, we tackle the problem by variational quantum Monte Carlo (VQMC), which takes care of the no-double occupancy condition exactly. This method has been used extensively previously for the same system, Gros88-VQMC ; Ogata99-VQMC ; Ogata02-VQMC ; Trivedi01-VQMC ; Trivedi04-VQMC ; Edegger07-review ; Tan08-twomode yielding considerable insights into the Neel state at half filling and the uniform $d$-wave SC state at finite doping. More recently, the VQMC has been extended to deal with essentially unlimited number of variational parameters. Sorella05-newton ; Umrigar05-newton ; Morita Here we will go beyond the uniform $d$-wave ansatz yet limit ourselves to a handful of motivated parameters as we now describe. For later convenience, we introduce the Nambu basis $\psi_{i}^{\dagger}=(c_{i\uparrow}^{\dagger},c_{i\downarrow})$. The variational Hamiltonian can be written as,

[TABLE]

Henceforth $\tau_{1,2,3}$ are Pauli matrices in the Nambu (or particle-hole) basis. The second-neighbor hopping $t_{2v}$, the chemical potential $\mu$ and the exchange field $\mu_{z}$ are all variational parameters, and $\eta_{i}=\pm 1$ is a staggered sign for A/B sublattice. Furthermore, on the $N_{1}$ bonds,

[TABLE]

Henceforth we fix $t_{1v}=1$ without loss of generality (since the only role of $H_{v}$ is to construct the trial wavefunction, see below), $\tau^{\pm}=(\tau_{1}\pm i\tau_{2})/2$, and $\Delta_{ij}$ is the real-space pairing function on a directed bond $\langle ij\rangle$: For ${\bf b}={\bf r}_{j}-{\bf r}_{i}=(b_{x},b_{y})$,

[TABLE]

where $S_{0}$ ($T_{\bf Q}$) is the singlet (triplet) part of the $d$-wave pairing, and $S_{\bf Q}$ ($T_{0}$) is the singlet (triplet) part of the chiral $p$-wave pairing. Note we assumed the triplets all have their $d$-vectors along $z$, the direction of the Neel moment. In this way, the total spin of a triplet Cooper pair is orthogonal to the Neel moment, a most favorable situation for triplets to develop on the Neel order induced by $\mu_{z}$. On the other hand, the momentum ${\bf Q}=(\pi,\pi)$ in the PDW state is obvious from the staggered sign $\eta_{i}$ over the lattice. The four cases of the pairing function $\Delta_{ij}$, when only one of the four coefficients in it is nonzero, are illustrated in Fig.1. To summarize, we consider the set of variational parameters

[TABLE]

Since the $d$-wave and $p$-wave states belong to different irreducible representations of the $C_{4}$ group, we consider them separately. We assume all parameters in $x$ are real, as this turns out to gain energy better.

The normalized trial ground state is constructed as

[TABLE]

where $|\psi_{0}\rangle$ is the ground state of the free variational Hamiltonian $H_{v}$ (which depends on the parameter set $x$), in the canonical ensemble with a definite total number of electrons, $N_{e}=N(1-\delta)$, on the lattice. Here $N$ is the number of sites and $\delta$ is the hole doping level. We use the standard Monte Carlo to calculate the average energy density $E$ and the Neel order $m$,

[TABLE]

We optimize the parameter set $x$ automatically by adapting to our case the method proposed previously.Sorella05-newton ; Umrigar05-newton ; Morita More technical details can be found in the Appendix. We typically consider (tilted) square lattices with $N=82$, $128$, and $200$, and the finite size effect is found to be insignificant by casual check up to $N=800$. To stabilize the wavefunction, we apply antiperiodic boundary condition for fermions along one or both directions, and we add a tiny local singlet pairing that works even without using the antiperiodic boundary condition.

III Results and discussions

III.1 Variational results at zero doping

In Fig.2 we show the variation of the energy (per site) versus the change $\Delta x$ of the parameter set $x$ during the automatic optimization of $x$ for the undoped system in (a) the $d$-wave state with $x=(\mu,\mu_{z},t_{2v},S_{0},T_{\bf Q})$, and (b) the $p$-wave state with $x=(\mu,\mu_{z},t_{2v},S_{\bf Q},T_{0})$. The process starts from an initial guess $x_{0}$, and reaches quickly the stationary point $x_{*}$, where the energy is minimized. The lingering around $x_{*}$ provides a measure of the statistical error in energy and in the optimized parameters. In both cases, the parameter $\mu_{z}$ is finite, leading to Neel order, see Fig.4(b) at zero doping. The best energy of the $d$-wave state in Fig.2(a) is about $E=-1.1645J$, or $\langle{\bf S}_{i}\cdot{\bf S}_{j}\rangle=-0.3322J$ on $N_{1}$ bonds, consistent with that reported in the literature.YokoyamaShiba88-0.321 ; Gros88-VQMC ; LeeFeng88-0.332 ; Lee06-review ; Edegger07-review Interestingly, the optimized energy is even lower in the $p$-wave state in Fig.2(b), $E=-1.1660J$. Translated as $\langle{\bf S}_{i}\cdot{\bf S}_{j}\rangle=-0.3330J$, the energy is so far the best using fermionic VQMC, and is fairly close to the results of bosonic VQMC ($-0.3344J$Liang-0.3344 ; Weng-0.3344 ), Green’s function QMC ($-0.3346J$ TrivediCepeley-0.3346 ) and stochastic series expansion QMC ($-0.3347J$Sandvik97-QMC ). Moreover, starting from a different initial guess $x_{0}$ we would obtain (not shown) a different optional $x_{*}$, but the energy is degenerate within statistical error. For example, in both cases of $d$-wave and $p$-wave, we obtain the same energy, respectively, by setting the triplet components to zero while optimizing the others. However, the energy is much poorer if we set the singlet components to zero instead. Specializing to the energetically more favorable $p$-wave case, we conclude the $p$-wave singlet PDW is relevant while the $p$-wave triplet is only marginal. This could be understood from the fact that locally the singlet on $N_{1}$ bonds gains energy from the anti-ferromagnetic spin exchange, while the triplet would be costly. The surprise is the singlet favors a center-of-mass momentum ${\bf Q}$, in contrast to the usual uniform $d$-wave ansatz (related to our parameter $S_{0}$). The latter was widely assumed in previous fermionic VQMC. To make sure that our result is not a finite-size artefact, we perform the optimization on lattices of various sizes, and the results are shown in Fig.2(c) as a scaling plot versus the linear size $L=\sqrt{N}$. The energy difference between the two types of states clearly survives in the limit of $1/L\rightarrow 0$, suggesting insignificant finite-size effect. Therefore, our results point to a novel type of ground states with $p$-wave singlet PDW. An interesting question is whether such a pairing could persist at finite doping. We will come back to this point in the next section.

III.2 A $\mathrm{Z_{2}}$ QSL with nested spinon Fermi surfaces and a pair of Dirac nodes

At zero doping, we find that if we switch off the Neel order (by setting $\mu_{z}=0$), the optimized energy is $E=-1.1466J$ for the $pS_{\bf Q}$ state, and $E=-1.1406J$ for the $dS_{0}$ state. The $pS_{\bf Q}$ state is still better. The energy difference is far beyond statistical error. Furthermore, we can set $\mu=t_{2v}=0$ without affecting the optimized energy. In this case the variational Hamiltonian is only composed of the $U_{ij}$-terms in Eq. (2). We take a closer look into such a $pS_{\bf Q}$ state to understand why it is better than $dS_{0}$, in terms of the PSG theory.Wen02-PSG ; Wen2004

There are staggered signs in $\Delta_{ij}$ or in $U_{ij}$ in the $pS_{\bf Q}$ state, see Fig.1(c). This can actually be gauged away, given the exact charge neutrality (at zero doping) and hence local charge-SU(2) gauge invariance in VQMC. After the gauge transformation

[TABLE]

where ${\bf r}_{i}=(x_{i},y_{i})$ is the coordinate vector, we obtain a uniform ansatz:

[TABLE]

with $U_{ji}=U_{ij}^{\dagger}$. The Wilson loop around an elementary plaquette counter-clockwise can be written as

[TABLE]

where ${\cal P}$ denotes a specific cyclic permutation for a given path, $U_{1,2}=i\tau_{0}+S_{\bf Q}\tau_{1,2}$, and $U_{3,4}=U_{1,2}^{\dagger}$. We find $W_{1234}=we^{i{\bf F}\cdot\tau/2}$, with $w$ a complex factor and ${\bf F}\propto-(S_{\bf Q},-S_{\bf Q},1)$. Since $[U_{1,2},W_{1234}]\neq 0$, the Wilson loops with respect to the same base point, obtainable by all cyclic permutations in Eq. (10), carry noncolinear fluxoid, hence the invariant gauge group (IGG) is $\mathrm{Z_{2}}$. This IGG dictates that the distortion to the ansatz Eq. (9) in the form of $U_{ij}\rightarrow U_{ij}e^{i{\bf a}_{ij}\cdot\tau/2}$ is gapped for all directions of the gauge field ${\bf a}_{ij}$. In comparison, there are massless U(1) gauge fluctuations in the $dS_{0}$ state.Wen02-PSG ; Wen2004 Note that while the energy from VQMC is strictly invariant under local gauge transformation of $U_{ij}$, it does change under the gauge distortion if ${\bf a}_{ij}$ cannot be gauged out. Therefore, the trial state $|G\rangle$ using a given ansatz $U_{ij}$ should work better when the gauge field is already massive, as in the case of Eq. (9), which is gauge equivalent to the (nonmagnetic) $pS_{\bf Q}$ state.

We further study the properties of our $\mathrm{Z_{2}}$ QSL. In the following we list the elements $g$ of the conventional symmetry groups, upon which $U_{ij}$ in Eq. (9) changes, but can be restored by a subsequent local gauge transformation $V_{g}$:

[TABLE]

In the above, $P_{x,y,d}$ are mirrors sending $x\rightarrow-x$, $y\rightarrow-y$ and $x\leftrightarrow y$, respectively, and $T$ is time reversal which acts on $U_{ij}$ as $i\tau_{2}U_{ij}^{*}(-i\tau_{2})=-U_{ij}$.Wen02-PSG ; Wen2004 Note the four-fold rotation is identical to $P_{x}P_{d}$, so the point group $C_{4v}$ is automatically covered. We see that our spin liquid enjoys full physical symmetries. The set $\{V_{g}\cdot g\}$ forms the PSG of the spin liquid, and can be labeled as Z2A$xy(12)n$. Wen02-PSG ; Wen2004 This PSG has not been realized in previous VQMC or self-consistent MF studies.

Using the ansatz in Eq. (9), the variational Hamiltonian becomes

[TABLE]

Here $\varepsilon_{\bf k}=2(\sin k_{x}+\sin k_{y})$. The quasiparticle dispersion is easily found to be

[TABLE]

Clearly, the spinon excitation is gapless, and in fact there are two Fermi pockets enclosing $\pm{\bf Q}/2$, see Fig.3. In addition, at the Fermi level there are two Dirac nodes at $\pm(\pi/2,-\pi/2)$. Importantly, the spinon Fermi surfaces and Dirac nodes are protected by the above PSG, and this serves as an indicator of the quantum order in such a gapless $\mathrm{Z_{2}}$ QSL. Wen02-PSG ; Wen2004

Finally, we discuss how the Neel order appears to lower the energy of the $pS_{\bf Q}$ state further. Even if we assume the spinons are free from the coupling to the massive gauge fields, residual interactions between spinons can induce an instability in the presence of perfect nesting between the spinon Fermi surfaces. This results in a finite exchange field $\mu_{z}$ staggered over A- and B-sublattices. The quasiparticle dispersion becomes, in the folded Brillouin zone, $E_{\bf k}^{\prime}=\pm\sqrt{E_{\bf k}^{2}+\mu_{z}^{2}}$. This eventually gaps out the spinons and breaks the spin rotation symmetry. We see an interesting example that the massiveness of gauge field fluctuations in a $\mathrm{Z_{2}}$ QSL is insufficient to guarantee its stability against magnetic ordering.

III.3 Finite doping

We have also performed systematic VQMC simulations at finite doping. In Fig.4(a) we show the energy versus hole doping in the various variational states. We find the energy in the $dS_{0}$ and $dS_{0}+dT_{\bf Q}$ states are degenerate within statistical error, and so are the $pS_{\bf Q}$ and $pS_{\bf Q}+pT_{0}$ states. In contrast, the energy becomes poorer (not shown) if we get rid of the singlet components in both cases. This further enforces our view that the relevant pairing states are all singlets: $dS_{0}$ and $pS_{\bf Q}$. Moreover, while in the undoped case we find the $pS_{\bf Q}$ (or $pS_{\bf Q}+pT_{0}$) state has lower energy, at finite doping, we find the $dS_{0}$ (or $dS_{0}+dT_{\bf Q}$) state is systematically more favorable, down to the lowest nonzero doping we accessed. In fact, by linear interpolation the transition between these two types of states would be at a tiny hole doping level. Fig.4(b) shows the average Neel moment versus the hole doping. We find while the energy may be degenerate in the $d$-wave case, or in the $p$-wave case, as shown in (a), the moment differs if the triplet component is included. For example, the Neel moment is larger in the $dS_{0}+dT_{\bf Q}$ state than the $dS_{0}$ state, and similarly for $pS_{\bf Q}+pT_{0}$ versus $pS_{\bf Q}$.

Taking the lower-energy $dS_{0}+dT_{\bf Q}$ state at finite doping, we believe the harmlessness of the triplet component $dT_{\bf Q}$ to the optimized energy points to soft triplet PDW fluctuations and spin fluctuations in the underdoped regime. The dominant $dS_{0}$ is just the $d$-wave singlet pairing, well perceived in doped cuprates. The secondary $dT_{\bf Q}$ state may become important at higher energy scales (e.g., above the superconducting transition temperature), where the umklapp scattering of such pairs may be the key process to generate the pseudogap near the antinodes in underdoped cuprates.maurice Unfortunately, this is already beyond the scope of VQMC for the ground state.

On the other hand, our results indicate the $pS_{\bf Q}+pT_{0}$ state is not favorable at finite doping. In particular, the chiral $pT_{0}$ state does not seem to help gain energy, at either zero or finite doping, with or without the Neel order. This appears to rule out the $pT_{0}$ state as a leading order parameter, which was proposed as an intriguing possibility on a tunable background of Neel order (and with tunable spin-exchange and Coulomb interaction on $N_{1}$-bonds).Lu14-topo Our VQMC treats the Neel order as one of the variational parameters for a given $t$-$J$ model. In another context, the $d$-wave pairing was shown by MFT to be stable, and the quasiparticle excitation becomes nodeless, on a strong Neel background.gm1 ; gm2 It is natural that the energetic ordering of the states in different symmetries depends on the starting model, and on the method by which the strong correlation effects are addressed.

IV conclusion

We performed VQMC study of the $t$-$J$ model for cuprates, putting on equal footing the singlet and triplet pairings at zero as well as at momentum ${\bf Q}=(\pi,\pi)$, along with the Neel moment. At zero doping we find the $p$-wave singlet PDW ($pS_{\bf Q}$) to be favorable on the background of the Neel moment. We also discussed in terms of charge-SU(2) gauge theory that if the Neel order is ignored, the $pS_{\bf Q}$ state is still favorable versus the $dS_{0}$ state, and describes a novel Z2A$xy(12)n$ QSL with spinon Fermi surfaces and Dirac nodes. However, this QSL is unstable toward Neel ordering. Above possibly a tiny doping level, we find the $dS_{0}+dT_{\bf Q}$ state is established. We also find the uniform $p$-wave triplet ($pT_{0}$) is not a leading order at zero or finite doping. Our results point to a novel $\mathrm{Z_{2}}$ QSL in proximity to the Neel order at zero doping, and to intertwined $d$-wave triplet PDW fluctuations and spin moment fluctuations along with the dominant $d$-wave singlet SC at finite doping. The PDW state $dT_{\bf Q}$ at finite doping may also provide the microscopic origin of the umklapp-scattering mechanism for the pseudogap maurice .

Two remarks are in order. First, we find our results at finite doping is qualitatively unchanged even if we set $t^{\prime}=0$. Second, at finite doping, the magnetic order may become more complicated, such as the stripe order or even incommensurate magnetic order.Keimer15-review Charge ordering is also possible. These orders reflect the fact that there are many intertwined or competing orders in strongly correlated systems such as the cuprates.Fradkin15-review These orders require a larger set of variational parameters to accommodate, and are left for further studies.

Acknowledgements.

The project was supported by the National Key Research and Development Program of China (under Grant Nos. 2016YFA0300401 and 2016YFA0300202), the National Natural Science Foundation of China (under Grant Nos.11574134, 11874205 and 11774306), and the Strategic Priority Research Program of Chinese Academy of Sciences (No. XDB28000000). The numerical simulations are performed in High-Performance Computing Center of Collaborative Innovation Center of Advanced Microstructures, Nanjing University.

V Appendix

The ground state $|\psi_{0}\rangle$ of $H_{v}$ can be obtained straightforwardly as follows. Formally, we rewrite $H_{v}$ in the Nambu space,

[TABLE]

where $M_{ij}$ is a matrix in the Nambu as well as real spaces. Suppose $\varphi_{n}$ is the BdG eigenstate of $M$,

[TABLE]

with eigen energy $\epsilon_{n}$. We write the eigenstate as $\varphi_{n}=(u_{n},v_{n})^{t}$, where $u_{n}$ is a column vector for the particle part, and $v_{n}$ a column vector for the hole part. There are in total $2N$ BdG states, and by particle-hole symmetry, the energies of these states appear pairwise in sign. We assume $\epsilon_{n}\leq 0$ for $n\leq N$, and construct the $N\times N$ matrices

[TABLE]

and subsequently,

[TABLE]

Then the ground state of $H_{v}$ can be written as

[TABLE]

where $N_{e}$ is the total number of electrons and $|0\rangle$ is the vacuum. The matrix $A$ depends on $x$ implicitly. In principle, $A$ can also be obtained first in the momentum space (with sublattice structure) and then transformed into the real space. But the above form is versatile. Finally, the normalized trial ground state for the $t$-$J$ model is given by Eq. (6).

We now describe how we minimize the energy $E$ by varying $x^{\mu}\in x$. Since we have a handful of variables to optimize, we try to perform the optimization automatically, rather than scanning over the parameter space de forte. The method has been developed in the literature. Umrigar05-newton ; Sorella05-newton ; Morita For self-completeness, here we review the method briefly, and we specify some subtleties to take care of in our case. Naively, the simplest way to update $x^{\mu}$ is the steepest descent method, $dx^{\mu}=-\partial_{\mu}Edt$, where $\partial_{\mu}E\equiv\partial E/\partial x^{\mu}$ and $dt$ is an artificial time step. However, this equation could cause instabilities in application. The reason is as follows. After $x$ has been updated as $x\to x+dx$, the wave function changes as $|G\rangle\rightarrow|G\rangle+dx^{\mu}|\mu\rangle$, with $|\mu\rangle\equiv\partial_{\mu}|G\rangle$. For brevity, summation over repeated indices is assumed henceforth, unless specified otherwise. Since the states $\{|\mu\rangle\}$ are not necessarily orthogonal, the effects of $\{dx^{\mu}\}$ are not independent, so that a small error in $dx^{\mu}$ may need a large $dx^{\nu}$ to compensate in a later stage. This is the root of instability. To make improvement, we form locally an orthonormal basis set $\{|a\rangle\}$ so that

[TABLE]

We observe that

[TABLE]

For real $dx^{a}$ and real metric tensor $g_{\mu\nu}=\langle\mu|\nu\rangle$, this is exactly the invariant line element (squared) on a curved manifold. We will consider the simpler real case for a moment, because of its appealing geometrical interpretation, leaving the complex case afterwards. We now project the energy gradient into the orthogonal frame,

[TABLE]

We then apply the steepest descent method in the orthogonal frame,

[TABLE]

We can now write this equation back to the $x^{\mu}$-frame,

[TABLE]

Multiplying both sides by $e^{a}_{\sigma}$ and contracting $a$, we obtain

[TABLE]

where $g^{\mu\nu}$ is the inverse metric tensor. This form is explicitly covariant, and could have been obtained by an educated guess.

The metric tensor would be real if the variational parameters couple to operators diagonal in real space, such as the Jastrow factors for Hubbard models. This is unfortunately not our case. We therefore need to consider more general cases of $dx^{a}$ and $g_{\mu\nu}$. By definition we have

[TABLE]

For a complex change $dx^{a}$, the steepest descent equation should be written as

[TABLE]

Going back to the $x^{\mu}$-frame, we obtain

[TABLE]

Multiplying both sides by $\langle\nu|a\rangle$ and summing over $a$, we obtain

[TABLE]

where we used the relations, by definition,

[TABLE]

Therefore we are led to the same result as Eq. (26). A possible problem is the resulting $dx^{\mu}$ may be complex. A simple solution to the problem is to disregard the imaginary part of $g_{\mu\nu}$ without jeopardizing the invariant line element in Eq. (22). The correctness and efficiency of doing so can only be judged by practice.

Moreover, even if the update is performed in the form of Eq. (26), instability could still appear if the metric tensor has small eigenvalues. For better stability in application, two further improvements are proposed in the literature. The first is to enlarge slightly the diagonal elements of the metric tensor, $g_{\mu\mu}\rightarrow(1+\epsilon)g_{\mu\mu}$ (no summation over $\mu$), with a positive small number $\epsilon$. The second is to disregard the dangerous eigenmodes of the metric tensor. Since $g_{\mu\nu}$ is Hermitian and positive definite, we can expand its inverse as

[TABLE]

where $\phi_{k}$ is an eigenvector of $g_{\mu\nu}$ with eigenvalue $\Lambda_{k}$, and the primed summation excludes small eigenvalues that would blow up statistical errors in $\partial_{\mu}E$ in Eq. (26).

We now specify how $\partial_{\mu}E$ and $g_{\mu\nu}$ are calculated. We recall that $|G\rangle$ can be expressed as

[TABLE]

Here $R=\{i_{\uparrow},j_{\downarrow}\}$ denotes a real-space configuration of the electrons without double-occupancy, and the coefficient $\alpha_{R}$ is the determinant of the sub-matrix $A_{i_{\uparrow}\in R,j_{\downarrow}\in R}$. In the $R$-basis, we define the matrix $h_{RR^{\prime}}=\langle R|H|R^{\prime}\rangle$, which is a sparse matrix. We take $|\alpha\rangle$ as the column vector composed of $\alpha_{R}$, or the wavefunction in the $R$-space, and define the average $\langle{\cal O}\rangle=\langle\alpha|{\cal O}|\alpha\rangle/\langle\alpha|\alpha\rangle$ for any (local or nonlocal) operator in $R$-space. We also define $l_{\mu}(R)=\partial_{\mu}\ln\alpha_{R}$ formally as a local operator in the $R$-space. (This operator would be ill-defined if the phase of $\alpha_{R}$ does not vary smoothly with the parameter $x$. For example, an arbitrary global phase could be assigned to a wavefunction without any physical significance. To avoid such an ambiguity, in practice it is advisable to construct the wavefunction as a continuous fiber over the parameter space.) In these formal terms, we have

[TABLE]

The last imaginary part drops out upon contraction with $dx^{\mu}dx^{\nu}$. The same is true in the imaginary part of $\langle l_{\mu}^{*}l_{\nu}\rangle$. So eventually we simply take the real part of $g_{\mu\nu}$, as we argued previously. All of the above averages can be calculated conveniently by Monte Carlo. We checked that the automatic update works perfectly for the $dS_{0}+dT_{\bf Q}$ state at any doping level, but meets instabilities for the $pS_{\bf Q}+pT_{0}$ state at finite doping. In the latter case, the wavefunction seems to be too sensitive at some isolated points of $x$, resulting in jumps of the energy followed by steady decay. In the worst case we perform scanning over each direction of the parameter space to obtain reliable results.