Projected BCS Theory for the Unification of Antiferromagnetism and Strongly Correlated Superconductivity

The intimate connection between antiferromagnetism and superconductivity is at the core of high-temperature superconductivity. Here, we put forward the projected BCS theory for the unification of antiferromagnetism at half filling and strongly correlated superconductivity at moderate doping. Specifically, it is shown that the projected BCS theory provides excellent trial states for the exact ground states of the $t$-$J$ model in the square lattice, generating the unified phase diagram as a continuous function of hole concentration. Precisely capturing antiferromagnetism at half filling, which is ultimately a consequence of the strong correlation between Cooper pairs, the projected BCS theory is able to produce better trial states for strongly correlated superconductivity at moderate doping than the resonating valence bond state. Finally, we discuss various ramifications of the projected BCS theory.

The intimate connection between antiferromagnetism and superconductivity is at the core of hightemperature superconductivity. Here, we put forward the projected BCS theory for the unification of antiferromagnetism at half filling and strongly correlated superconductivity at moderate doping. Specifically, it is shown that the projected BCS theory provides excellent trial states for the exact ground states of the t-J model in the square lattice, generating the unified phase diagram as a continuous function of hole concentration. Precisely capturing antiferromagnetism at half filling, which is ultimately a consequence of the strong correlation between Cooper pairs, the projected BCS theory is able to produce better trial states for strongly correlated superconductivity at moderate doping than the resonating valence bond state. Finally, we discuss various ramifications of the projected BCS theory.
Despite the large variety in material properties, there is a certain list of common features robust across hightemperature superconductors. One of the most salient features in such a list is the proximity between antiferromagnetism and superconductivity. An important question is exactly how these two phenomena are connected together.
The spin singlet is the smallest unit of antiferromagnetism. Since the BCS state is the condensate of spinsinglet electron pairs, i.e., Cooper pairs, it is rather natural that antiferromagnetism, at least with short-range order, is closely connected with superconductivity. The resonating valence bond (RVB) state is the trial state constructed under this very rationale. Specifically, the RVB state is the projected BCS state with all the components containing doubly occupied sites projected out [1][2][3][4][5][6][7][8]. Known as the Gutzwiller projection when fully applied, this projection can be also partially imposed [9][10][11].
A problem is that the RVB state reduces to the spin liquid with fluctuating spin singlets rather than the longrange-ordered antiferromagnet at half filling. Quantum antiferromagnets in two dimensions can, in principle, have various types of spin liquids and valence bond solids [12]. However, does this mean that antiferromagnetism with true long-range order is necessarily incompatible with superconductivity?
Here, we show that it is possible to construct a unified theory of antiferromagnetism (AF) and strongly correlated superconductivity (SCSC). A gist of the unified theory is that the Gutzwiller projection, not commuting with the BCS Hamiltonian, should be treated better. Specifically, in this theory, the Gutzwiller projection is applied onto the BCS Hamiltonian itself: the BCS Hamiltonian is directly diagonalized in the Gutzwillerprojected Hilbert space [13,14]. For convenience, let us call this theory the projected BCS theory.
To provide the validity of the projected BCS theory, we perform exact diagonalization of the projected BCS Hamiltonian and compare the so-obtained exact ground state with that of the t-J model. Actually, the projected BCS theory has been previously analyzed by using the similar numerical technique [13,14]. While providing promising results, however, the previous analysis has suffered from finite-size effects. Here, we overcome this problem by devising the proper overlap in the presence of particle number fluctuations. As a result, it is shown that the projected BCS theory provides excellent trial states for the exact ground states of the t-J model in the square lattice for a wide range of hole concentration. In particular, we provide a rigorous proof for the equivalence between the exact ground states of the t-J model and the projected BCS theory at half filling. More importantly, we obtain the unified phase diagram of the t-J model as a continuous function of hole concentration, capturing both AF at half filling and SCSC at moderate doping. The Fermi liquid state is also naturally captured at sufficiently large doping. This phase diagram provides concrete evidence for the split between the pairing amplitude and the real superconducting order parameter at low doping, which can in principle explain the pseudogap phenomenon.
It is emphasized that, precisely capturing AF at half filling, the projected BCS theory is able to produce better trial states for SCSC at moderate doping than the RVB state, demonstrating the merit of the unification of AF and SCSC, as advocated by the SO(5) theory [15]. Finally, the projected BCS theory can be used to reveal the relevance of the s-wave pairing symmetry as well as the geometrical frustration to high-temperature superconductivity.

Results
High-temperature superconductivity and the fractional quantum Hall effect.
There is a remarkable similarity between the RVB state for arXiv:1912.13064v1 [cond-mat.str-el] 30 Dec 2019 high-temperature superconductivity and the composite fermion (CF) state [16] for the fractional quantum Hall effect (FQHE). Simply put, the RVB state is the BCS state projected onto the Hilbert space with no double occupancy: where P G is the Gutzwiller projection operator, and ψ BCS is the BCS wave function describing the electronelectron paired state composed of Cooper pairs. Similarly, the CF state is the electron-vortex bound state projected onto the lowest Landau level (LLL): where P LLL is the LLL projection operator, and ψ Jastrow is the Jastrow-correlated wave function describing the electron-vortex bound state composed of CFs [16]. The electron-electron pairing in the RVB state is a consequence of the intricate interplay between the no-doubleoccupancy constraint and the residual hopping effect. Similarly, the electron-vortex binding in the CF state is a consequence of that between the LLL constraint and the residual electron-electron interaction. One of the most important and convenient properties of these two states is that both Cooper pairs and CFs themselves can be treated as weakly interacting. The 'plain vanilla' version of the RVB theory assuming essentially non-interacting Cooper pairs can provide a semiquantitative understanding of high-temperature superconductivity [8]. Similarly, the FQHE can be understood as the integer quantum Hall effect of weakly interacting CFs [16]. The trial wave function based on weakly interacting CFs has been incredibly successful, basically explaining all known properties of the FQHE [17].
The similarity between high-temperature superconductivity and the FQHE goes even deeper in that there are fundamental difficulties in the field theoretical implementation of the required projection operators, P G and P LLL , respectively. Specifically, various gauge theories, often in terms of slave particle methods, have been widely used to implement P G [18]. Similarly, the Chern-Simons gauge theory has provided a reasonable low-energy description of the FQHE (and also the CF sea state) by implementing P LLL in the mean-field and one-loop levels [19][20][21]. Despite these successes, however, it is necessary to treat the projection operators better to capture the effects of strong correlation more accurately.
There are several important examples in the FQHE, where the correlation between CFs cannot be ignored. The first example is the Wigner crystal occurring at low filling factors in the LLL [22], where CFs, rather than electrons, form the Wigner crystal [23,24]. The second example is the famous 5/2 state, the even-denominator FQHE state occurring in the half-filled second Landau level [25]. Let us elaborate on the second example below.
It is generally believed that the 5/2 state is a paired quantum Hall state [26][27][28][29][30]. Despite this general belief, however, there are various crucial questions to be answered. Among such questions are the proof of pairing other than the gap, the pairing symmetry, the validity of the non-Abelian braiding statistics, and so on. Considering the effectiveness of the trial state method in the FQHE, the better trial state can be constructed, the more conclusive answers can be obtained.
One of the most promising trial states for the 5/2 state is the Moore-Read (MR) Pfaffian state [26,27] along with its particle-hole (PH) conjugate, called the anti-Pfaffian state, and their PH-symmetrized combination [31,32]. While the MR Pfaffian state can be written explicitly in the planar geometry, it is generally more convenient to construct it by directly diagonalizing the model Hamiltonian, whose exact ground state is known to be the MR Pfaffian state. Such a model Hamiltonian is the threebody repulsive interaction, often denoted as H 3 , imposing an energy penalty to closely-packed three-electron clusters. Consequently, checking the validity of the MR Pfaffian state amounts to computing the overlap between the exact ground states of the Coulomb interaction and H 3 [31,32].
By the same token, to capture the long-range antiferromagnetic order at half filling, one needs to better take into account the strong correlation between Cooper pairs. To this end, we consider the projected BCS Hamiltonian.
Projected BCS Hamiltonian. The projected BCS Hamiltonian can be written as follows: where where t is the hopping parameter, ∆ ij is the pairing amplitude, and µ is the chemical potential. For the d-wave pairing symmetry, ∆ ij = ±∆ for the x and y directions, respectively. For the s-wave pairing symmetry, ∆ ij = ∆ for both directions. The particle number operator n i is given by n i = σ c † iσ c iσ . In the meantime, the t-J Hamiltonian can be written as follows: where H t is the same as in Eq. (4), and where J is the spin exchange coupling constant with J > 0 [33][34][35]. Note that, unless mentioned otherwise, all the energy scales are measured in units of t throughout this work. Similar to what is done for the 5/2 state, we would like to compute the overlap between the exact ground states of the t-J model and the projected BCS theory. Before performing actual computations, however, it is important to check the symmetries of the Hamiltonians, which can crucially affect the overlap. To this end, we perform group theoretical analyses, whose details are provided in Methods.
More importantly, there is a certain limit, where the value of the overlap is exactly known. That is, the overlap should be exactly unity at half filling.
Equivalence at half filling. The equivalence between the t-J model and the projected BCS theory at half filling can be proved analytically. With the details of the proof presented in Methods, here, we provide an overall summary.
We begin by decomposing the Schrödinger equation of the projected BCS theory in terms of a system of simultaneous equations connecting between different number sectors. Then, we take the limit of large µ to maximize the electron number under the no-double-occupancy constraint, which leads to half filling. In this limit, we carefully expand the weighting amplitude of the energy eigenstate in each particle number sector in orders of 1/µ. Finally, by keeping only the most dominant terms, it is shown that the projected BCS theory becomes entirely equivalent to the Heisenberg model, i.e., the t-J model at half filling.
Several aspects of the proof are worthwhile to mention. First, the equivalence is not just for the ground state, but rather the entire energy eigenstates at half filling. Second, the proof works regardless of the pairing symmetry or the lattice structure. As shown by actual computations, the overlap is exactly unity at half filling not only for the d-wave pairing symmetry, but also for the s-wave pairing symmetry in the square lattice. The same is true for the triangular lattice. Third, the equivalence at half filling does not necessarily guarantee a high overlap upon doping. The overlap at finite doping depends crucially on the pairing symmetry and the lattice structure. It is shown that the overlap can remain high with finite ∆ for the d-wave pairing symmetry in the square lattice, but not for the s-wave pairing symmetry. The overlap is generally quite low at moderate doping in the triangular lattice.
To actually compute the overlap between the exact ground states of the t-J model and the projected BCS theory, now, let us consider the proper overlap in the presence of particle number fluctuations.
Proper overlap in the presence of particle number fluctuations. The exact ground state of the projected BCS theory can be expanded in terms of its com-ponents in various particle number sectors: where h is the number of holes, related with those of sites, N s , and electrons, N e , via h = N s − N e . The RVB state can be also similarly expanded. Meanwhile, any exact ground state of the t-J model has a fixed number of particles. This means that each particle number sector has its own exact ground state of the t-J model, forming a set of {|φ t-J h }. Now, the question is how to compute the overlap between |ψ PBCS and {|φ t-J h }. A traditional method is to project |ψ PBCS onto a specific particle number sector and compute its overlap with |φ t-J h in that sector, i.e., Let us call this the number-projected overlap squared (NPOS). As shown previously, the NPOS can provide a reasonable measure of the overlap except that there is a certain degree of arbitrariness in choosing the particle number. In other words, the average particle number is already implicitly determined by the chemical potential µ in H PBCS . Therefore, it is redundant to choose both particle number sector and chemical potential. We fix this problem by weighting the NPOS over all possible particle number sectors: which we call the number-weighted overlap squared (NWOS). Note that the weighting factor can be computed via |λ h | 2 = ψ PBCS |P h |ψ PBCS , where P h is the number projection operator. It is instructive to perform an explicit analysis of the 2 × 2 system to demonstrate how the NWOS is actually implemented. This analysis is also useful to illuminate the equivalence between the t-J model and the projected BCS theory at half filling. See Methods for details. Below, we present the numerical analysis of the 2 × 2 and larger systems via full-fledged exact diagonalization.
Unified phase diagram. Figure 1 shows the phase diagram of the t-J model in the square lattice determined via the NWOS between the exact ground states of the t-J model and the projected BCS theory with the d-wave pairing symmetry. Note that the NWOS is computed as a function of ∆ and µ, both of which are varied for all possible values. Also, note that N s is appropriately chosen for the proper tessellation of the square lattice. See Methods for details.
There are several important features to be noted. First, as expected from the equivalence at half filling, the NWOS is unity regardless of ∆ in the limit of large Phase diagram of the t-J model in the square lattice determined via the number-weighted overlap squared (NWOS) between the exact ground states of the t-J model and the projected BCS theory with the d-wave pairing symmetry. The spin exchange coupling constant of the t-J model is set to be J/t = 0.5. Denoted in color with bright yellow being unity and black being zero, the NWOS is computed as a function of pairing amplitude ∆ and chemical potential µ. As one can see, the NWOS is exactly unity regardless of ∆ at sufficiently large µ, proving that the long-range antiferromegnetic order is precisely captured by the projected BCS theory. Note that the NWOS also becomes unity as µ approaches −∞, i.e., the vacuum, where both spin exchange and pairing terms play negligible roles. The number of sites is varied as Ns  µ. This confirms that AF is precisely captured by the projected BCS theory.
Second, the NWOS becomes also unity in the limit of negatively large µ, corresponding to the vacuum. Here, the exact ground states of the t-J model and the pro-jected BCS theory are trivially identical since both spin exchange and pairing terms play negligible roles. Near this limit, the ground state is the usual Fermi liquid state.
Finally, clearly being non-zero, the optimal ∆ producing the maximum NWOS forms a rather well-defined curve as a function of µ, especially for larger systems at N s = 16 and 20. Note that, at N s = 18, there is a large black dome around µ = 0, masking the otherwise high NWOS region. The high NWOS can be obtained in this region if the first excited state of the projected BCS theory, which is nearly degenerate with its exact ground state, is used to compute the NWOS. In conclusion, there is a well-defined curve of the non-zero pairing amplitude as a function of doping. Actually, it is physically meaningful to express the phase diagram as a function of hole concentration x rather than µ. To this end, we convert µ to x by inverting the relationship x = 1 − i n i /N s with n i implicitly depending on µ. Figure 2 shows the phase diagram of t-J model in the square lattice as a function of ∆ and x. It is now possible to determine the optimal pairing amplitude producing the maximum NWOS, ∆ max , as a function of x. As one can see from Fig. 2, ∆ max forms a well-defined curve clearly lifted from ∆ = 0, especially at N s = 16 and 20.
Note that the spin exchange coupling constant of the t-J model is fixed to be J/t = 0.5 in Figs. 1 and 2. See Methods to find how the phase diagram changes with variation of J/t. Now, one may wonder if the non-zero value of ∆ actually means superconductivity. To answer this question, we compute the superconducting order parameter by using the exact ground states of the projected BCS theory with ∆ max determined above as a function of x.
Superconducting order parameter. One of the most advantageous features of the projected BCS theory is that its exact ground states have intrinsic particle number fluctuations. Thus, it is possible to directly compute the superconducting order parameter, c i↑ c j↓ , instead of taking the large-distance limit of the off-diagonal longrange order (ODLRO), c i↑ c j↓ c † k↓ c † l↑ , with i (k) and j (l) being the nearest neighbors, but i and k taken to be largely separated.
Specifically, we compute the superconducting order parameter as follows. First, we obtain ∆ max as a function of x as shown in Fig. 3 (a). Then, we compute c i↑ c j↓ by using the exact ground states of the projected BCS theory with ∆ max as an input parameter. Figure 3 (b) shows the resulting superconducting order parameter. It is important to note that the superconducting order parameter vanishes as x decreases in contrast to ∆ max , which approaches a finite value. This split between the pairing amplitude and the real superconducting order parameter is due to the strong correlation between Cooper pairs.
To elucidate the role of the strong correlation more clearly, we compare c i↑ c j↓ with the "bare" superconducting order parameter c i↑ c j↓ 0 , which would be the superconducting order parameter if there were no Gutzwiller projection: where E 0,k = ξ 2 0,k + ∆ 2 0,k with ∆ 0,k = 2∆ max (cos k x − cos k y ) and ξ 0,k = −2t(cos k x + cos k y ) − µ 0 . Above, the translational invariance is assumed. Note that µ 0 is determined from the condition that the "bare" hole concentration is the same as x: 1 Ns k ξ 0,k /E 0,k = x. As shown in Fig. 3 (b), the bare superconducting order parameter does not vanish at low doping similar to ∆ max .
While important to identify superconductivity, the superconducting order parameter is not a direct physical observable by itself. For a physical observable, the expectation value of the pairing energy, ∆ max c i↑ c j↓ , is  shown in Fig. 3 (c). Similar to the superconducting order parameter, the expectation value of the pairing energy vanishes at low doping while its bare counterpart, ∆ max c i↑ c j↓ 0 , does not. It is interesting to mention that a similar split between the pairing amplitude and the real superconducting order parameter was previously observed [5,7].
Comparison with the RVB state. To put the results of the projected BCS theory into prospective, we perform the similar NWOS analysis for the RVB state.
To this end, it is necessary to construct the RVB state precisely tailor-made for each finite system. Specifically, we would like to find the amplitude of the RVB state for each given basis state with spin up and down electrons located in {r ↑ , r ↓ }, respectively: whereg ij =g(r i − r j ) is the Fourier transform of g k = ∆ k /(ξ k + E k ) with E k = ξ 2 k + ∆ 2 k , ∆ k = 2∆(cos k x − cos k y )+δ, and ξ k = −2t(cos k x +cos k y )−µ. The indices  Figure 4 shows the phase diagram of the t-J model in the square lattice determined via the NWOS between the exact ground state of the t-J model and the RVB state with the d-wave pairing symmetry. Due to its inability to capture the long-range antiferromagnetic order, the RVB state has zero NWOS with the exact ground state of the t-J model at sufficiently large µ for N s = 4, 10, 16, 20. For N s = 8 and 18, the RVB state is not even available beyond certain large µ, where the Gutzwiller projection completely annihilate the BCS state, making the RVB state a null state.
Similar to what is done for the projected BCS theory, the RVB phase diagram can be also plotted as a function of x in Fig. 5. Again, the RVB state has zero NWOS with the exact ground state of the t-J model at half filling. Away from half filling, overall, the RVB phase diagram shows irregular behaviors depending on N s with broad vacant regions. Even when the RVB phase diagram behaves reasonably well at N s = 10 and 20, the optimal pairing amplitude of the RVB state exhibits a rather different behavior from that of the projected BCS theory at low doping. That is, the optimal pairing amplitude of the RVB state vanishes as x approaches zero.
Before discussing the different behaviors of the optimal pairing amplitude, however, it is important to check which among the two states, the RVB state and the exact ground state of the projected BCS theory, has the higher overlap with the exact ground state of the t-J model as a function of x. Figure 6 shows the maximum NWOS, O 2 max , of the projected BCS theory in comparison with that of the RVB state. As one can see, O 2 max of the projected BCS theory is pinned to be exactly unity at half filling, while that of the RVB state approaches a seemingly random value as x → 0 and becomes en-tirely zero strictly at x = 0. More importantly, O 2 max of the projected BCS theory is reduced somewhat upon initial doping but bounces back to a high value as x increases. Meanwhile, O 2 max of the RVB state is generally lower than that for the projected BCS theory except for a narrow range of x at low doping. It is interesting to speculate if the RVB state can provide a good trial state at this range of doping, where the spin glass phase can occur.
In summary, the projected BCS theory provides better trial states for the ground state of the t-J model than the RVB state, capturing both AF at half filling and SCSC at moderate doping.
S-wave pairing symmetry. As mentioned previously, the equivalence at half filling does not necessarily guarantee a high overlap upon doping, while so for the d-wave pairing symmetry in the square lattice. We would like to investigate what happens for the s-wave pairing symmetry in the square lattice. Figure 7 shows that, in this case, the NWOS is always maximized along ∆ = 0 for all finite µ, meaning that the s-wave pairing cannot be formed at finite doping. Note that the NWOS is unity regardless of ∆ in both limits of µ = ∞ and −∞, corresponding to half filling and vacuum, respectively. According to the proof of equivalence, the projected BCS theory is entirely equivalent to the Heisenberg model at half filling regardless of the pairing symmetry. The vacuum states are trivially identical.
It is important to note that the maximum NWOS of the s-wave pairing symmetry along ∆ = 0 is bound to be lower than that of the d-wave pairing symmetry occurring at the optimal curve of ∆ in Fig. 1, considering both pairing symmetries generate the same Hamiltonian at ∆ = 0. In conclusion, the d-wave pairing symmetry is preferred over the s-wave counterpart for the t-J model in the square lattice.
Geometrical frustration. Now, we investigate how the overlap depends on the lattice structure, especially with geometrical frustration. Specifically, we perform the NWOS analysis for the projected BCS theory in the triangular lattice.
Similar to the square lattice, we consider both s-or d-wave pairing symmetries. While simply a constant for the s-wave pairing, the pairing amplitude is given as ∆ ij = ∆e 2iθij for the d-wave pairing with θ ij being the angle between the x axis and the nearest-neighborconnecting vector r ij = r i − r j . Figure 8 shows that, at moderate doping, the projected BCS theory provides poor trial states for the t-J model in the triangular lattice regardless of the pairing symmetry.
In conclusion, geometrical frustration, at least, in the form of the triangular lattice is detrimental to the formation of superconductivity in contrast to the original rationale behind the RVB state. This result is consistent with the experimental observation that high-temperature superconductivity occurs almost always in tandem with  unfrustrated AF at half filling.

Discussion
In this work, it is shown that the projected BCS theory with the d-wave pairing symmetry can provide excellent trial states for the exact ground states of the t-J model in the square lattice for a wide range of hole concentration. In particular, the unified phase diagram, capturing both AF at half filling and SCSC at moderate doping, is obtained by computing the overlap between the exact ground states of the t-J model and the projected BCS theory with the d-wave pairing symmetry. A main breakthrough in this work is to devise the overlap properly taking into account particle number fluctuations.
By utilizing the natural presence of particle number fluctuations, the superconducting order parameter can be computed directly from the exact ground state of the projected BCS theory. The split between the optimal pairing amplitude and the superconducting order parameter can in principle provide an explanation for the pseudogap phenomenon at low doping.
For future work, it would be interesting to investigate various spectral properties of the projected BCS theory by computing both normal and anomalous Green's functions in the presence of intrinsic particle number fluctuations. This investigation could shed light on the issue of strange metal behaviors.
Methods Symmetries. Group theoretical analyses can facilitate exact diagonalization. Specifically, the basis states expanding the Hilbert space can be conveniently categorized into appropriate symmetry sectors, blockdiagonalizing the Hamiltonian.
In this work, we are interested in the sub-Hilbert space with both total momentum and z-component of the total spin being zero. In this situation, it is sufficient to focus on the point group symmetries.
Our goal is to understand the point group symmetries of the t-J and projected BCS Hamiltonians. Note that, being a local operator, the Gutzwiller projection commutes with all the point group symmetries. Therefore, one does not need to consider the Gutzwiller projection for the discussion of the point group symmetries.
Let us first discuss the point group symmetries in the square lattice. The point groups representing the symmetries of the square lattice are given by where e is the identity operator, C n the rotation operator about the z axis by the angle of 2π/n, σ x/y the reflection operator with respect to the x/y axis, and σ +/− d the reflection operator with respect to the y = ±x line. With spin degrees of freedom, the point groups should be enlarged to include the magnetic (or color) group [36]. In other words, one needs to consider the effects of the spin-flip operator f . Specifically, the point groups of the t-J model in the square lattice can be determined by the fact that both H t and H J commute with all the symmetry operators in C 4v as well as the spin-flip operator f . In other words, the t-J model has C II 4v ≡ {e, f } C 4v as its point groups in the square lattice.
Meanwhile, the point groups of the projected BCS Hamiltonian depend on the pairing symmetry. To understand this, it is important to note that H ∆ does not commute with some of the symmetry operators in C II Consequently, in the square lattice, the point groups of the projected BCS Hamiltonian with the s-wave pairing symmetry are simply given by C I 4v ≡ C 4v , excluding the spin-flip operator. In the case of the d-wave pairing symmetry, the spin-flip operator can be combined with some of the symmetry operators in C 4v . As a result, it can be shown that the point groups of the projected BCS Hamiltonian with the d-wave pairing symmetry are given by Note that the point groups of the projected BCS Hamiltonian (C I 4v and C III 4v for the s-and d-wave pairing symmetries, respectively) are only the subgroups of those of the t-J Hamiltonian (C II 4v ). In general, if two Hamiltonians have different symmetries, the overlap between their ground states would be suppressed.
Fortunately, the effects of both f and C 4 on H ∆ can be absorbed into the sign change of the pairing amplitude ∆. The sign, more generally, the phase of the pairing amplitude only affects the relative phase of the superconducting ground state between different particle number sectors, while keeping the component of the state in each sector invariant. Consequently, if properly defined to take into account particle number fluctuations, the overlap would not be automatically suppressed. See the main text for details.
Finally, let us consider the point group symmetries in the triangular lattice. The point groups representing the symmetries of the triangular lattice are C 6v . Since f still works as a symmetry element for the t-J model, the point groups of the t-Jmodel are given by C II 6v ≡ {e, f } C 6v . Similar to the square lattice, the projected BCS Hamiltonian with the s-wave pairing symmetry has C I 6v ≡ C 6v in the triangular lattice due to the lack of the spin-flip symmetry. Even worse, the projected BCS Hamiltonian with the d-wave pairing symmetry has only C 2 as its point groups in the triangular lattice.
Proof of equivalence at half filling. The energy eigenstates of the projected BCS theory can be expanded in terms of their component in each particle number sector: where h is the hole number, related with the site number N s and the electron number N e via h = N s −N e . λ h is the weighting amplitude in each hole number sector. Note that h is always an even number since we are interested in the paired state. For convenience, let us rewrite the Schrödinger equation for the projected BCS theory as follows: where H t,c,a = P G H t,c,a P G with H c and H a being the creation and annihilation parts of H ∆ , respectively. The particle number operator is defined as N = i n i . Then, the Schrödinger equation can be written component by component: (i) for h = 0, and (iii) for h = N s , Half filling can be obtained in the limit of infinite chemical potential, µ → ∞. In this limit, to satisfy Eq. (14), E = −N s µ + O(1/µ α ) with α ≥ 0. Note that λ h cannot diverge since they are the coefficients of a normalized wave function. Meanwhile, according to Eq. (16), O(λ Ns−2 ) = O(µλ Ns ). Under the assumption that λ h = O(1/µ ν h ), this means that ν Ns = ν Ns−2 + 1.
In this scaling, one can rewrite the two component-bycomponent equations for h = 0 and 2 by keeping only the most dominant terms: which can be combined to generate the Schrödinger equation solely at h = 0 as follows: Under the no-double occupancy constraint imposed by the Gutzwiller projection, one can show that where the second line is obtained by examining the matrix elements for all possible spin configurations in the (i, j) sites, i.e., | ↑↑ , | ↑↓ , | ↓↑ , | ↓↓ . Consequently, the final Schrödinger equation at half filling becomes as follows: which is nothing but the Hamiltonian of the Heisenberg model.
It is interesting to note that the equivalence between the projected BCS theory and the Heisenberg model at half filling has been previously argued by using a different analytical method [14], where the no-double-occupancy constraint is relaxed so that the projected BCS Hamiltonian is replaced by the BCS Hamiltonian with finite on-site repulsive interaction U . Eventually, the equivalence is obtained in the limit of large U . Note that the equivalence has been also argued by using a path-integral approach [37].
Explicit analysis of the 2 × 2 system. The entire Hilbert space of the 2 × 2 system can be expanded by the following seven basis states: (i) three states in the half-filled sector, (ii) three states in the two-hole sector, and finally (iii) one state in the vacuum sector, Note that we are interested in the Hilbert space with zero momentum (i.e., translationally invariant) and zero z-component of the total spin (i.e., spin-flip invariant). In terms of the ordered basis set {|e 1 , |e 2 , |e 3 }, the t-J Hamiltonian can be written in the half-filled sector as follows: generating the following ground state: where N 0 = 2/3 is the normalization constant. Note that |e 3 couples with neither |e 1 nor |e 2 since they have different point group symmetries. Similarly, in terms of the ordered basis set {|e 4 , |e 5 , |e 6 }, the t-J Hamiltonian can be written in the two-hole sector as follows: generating the following ground state: where α = √ J 2 + 32t 2 − J, β = 4 √ 2t, and N 2 = 1/ |α| 2 + |β| 2 . Similar to half filling, |e 6 couples with neither |e 4 nor |e 5 .
Meanwhile, the projected BCS Hamiltonian with the d-wave pairing symmetry can be written as the following block diagonal form in terms of the ordered basis set {|e 1 , |e 2 , |e 4 , |e 5 , |e 3 , |e 6 , |e 7 }: The characteristic polynomial of each block matrix is given as follows: where the subscript −/+ indicates the parity associated with the π/2 rotation followed by the spin flip. After a careful examination, one can show that the ground state occurs in the negative parity sector if µ gets sufficiently large. In this situation, it is convenient to observe that there is an easy solution of C − (ε), whose eigenvalue is −4µ with the corresponding eigenstate given by |e 1 + √ 2|e 2 . The other three eigenstates in the negative parity sector including the ground state should be orthogonal to this eigenstate. Therefore, the ground state at sufficiently large µ takes the following form: where N = 1/ 3/2 + |f | 2 + |g| 2 is the normalization constant, and f and g are some appropriate functions of ∆/t and µ/t, whose details are not provided here except that |f |, |g| 1 for sufficiently large µ. It is important to note that |ψ PBCS has the component of |e 1 − 1 √ 2 |e 2 in the half-filled sector, which coincides exactly with the ground state of t-J model in the same sector, |φ t-J 0 , as expected from the proof of equivalence at half filling.
Actually, it is interesting to perform the similar analysis for the projected BCS theory with the s-wave pairing symmetry, whose ground state also coincides exactly with the ground state of t-J model at half filling for basically the same reason mentioned above. Thus, there is a certain degree of robustness in the equivalence between the projected BCS theory and the Heisenberg model at half filling.
Now, an important question is if |ψ PBCS can capture the ground state of t-J model accurately as µ is reduced. That is, how accurately f and g in Eq. (31) can capture α and β in Eq. (28) as a function of µ? This can be checked via the NWOS: As one can see from Fig. 1 at N s = 4, the NWOS remains high upon doping. Note that there is a phase transition along the certain critical line of µ and ∆, beyond which the ground state of the projected BCS theory is no longer given by Eq. (31). While providing various successful small-system checks, the 2 × 2 system is way too small to determine if the optimal pairing amplitude can become non-zero at finite doping. To this end, one has to analyze larger systems by using full-fledged exact diagonalization.
Tessellation of the lattice. Let us begin with the tessellation of the square lattice, which can be tiled with tilted square unit cells. The tilted square unit cells are considered to observe the rotational symmetry of the lattice.
The distance between any two adjacent tilted square unit cells can be written as d sq = √ n 2 + m 2 (in units of lattice constant) for any non-negative integers n and m without the loss of generality, which is also the side length of the titled square unit cell. This means that the number of sites included in the tilted square unit cell is N s = d 2 sq = n 2 + m 2 . Since we are interested in the systems with even number of particles with zero z-component of the total spin, both m and n should be either even or odd integers; N s = 2, 4, 8,10,16,18,20,26,32,34,36, · · · .
For illustration, some of the tilted square unit cells are depicted in Fig. 9. Note that, if either n or m is zero (e.g., N s = 4 and 16), or n = m (e.g., N s = 8 and 18), the tessellation with tilted square unit cells has additional reflection symmetries; one for the horizontal axis, one for the vertical axis, and two for the diagonal axes. That is, the point groups become C 4v for N s = 4, 8, 16, 18. Due to the exponential growth of the Hilbert space, we are able to perform exact diagonalization only up to N s = 20 in this work. Figure 10 shows the number of basis states in the common logarithm scale, log N b , as a function of hole concentration x = 1 − N e /N s . Note that, here, the number of basis states is computed within the restricted Hilbert space, where both total momentum and z-component of the total spin are zero without using any other point group symmetries. As one can see, the total number of basis states can be more than 10 millions for N s = 20. Roughly speaking, the total number  of states increases by one order of magnitude as the number of sites increases by two. This means that the next available system at N s = 26 would have more than 10 billion basis states, which are beyond the current computing capacity. Now, let us switch gears to the tessellation of the triangular lattice. Similar to the square lattice, the distance between any two adjacent tilted hexagonal unit cells can be written as d tr = √ n 2 + nm + m 2 (in units of lattice constant) for any non-negative integers n and m without the loss of generality. After a moment of deliberation, one can show that the number of sites included in the tilted hexagonal unit cell is also given by N s = d 2 tr = n 2 + nm + m 2 . Both m and n should be even integers for N s to be an integer number; N s = 4, 12, 16, 28, 36, · · · . Meanwhile, the 120-degree spin order can exist only when m = n (mod 3). Combined all together, this means that the only accessible system in this work via exact diagonalization is that with N s = 12.
Phase diagram of the t-J model with variation of J/t. Figure 11 shows how the phase diagram of the t-J model changes with variation of J/t. As one can see from Fig. 11 (a), the optimal pairing amplitude producing the maximum NWOS is obtained along ∆ = 0 for small J/t, meaning that there is no pairing correlation at this parameter regime. The optimal pairing amplitude begins to be lifted from ∆ = 0 beyond a critical value of J/t roughly larger than 0.3. As J/t increases further, the optimal pairing amplitude also gets higher, expanding the regime of superconductivity.
Too large values of J/t may not be physically meaningful in the single-band model, where the t-J model is obtained as the large-U expansion of the Hubbard model with J/t being proportional to t/U [33,34]. It is, however, worthwhile to note that the spin exchange coupling can be also generated by the d-p hybridization mechanism in the three-band model, where the value of J/t can go beyond the t/U scaling [35] Construction of the RVB state. We begin with the usual form of the BCS state, which can be given in the momentum space as follows: where u 2 k = 1 2 (1 + ξ k /E k ) and v 2 k = 1 2 (1 − ξ k /E k ) with ξ k = k − µ and E k = ξ 2 k + ∆ 2 k .g ij =g(r i − r j ) is the Fourier transform of g k = v k /u k , and N is the normalization constant. Note that the second line is obtained under the condition that u k = 0 for all k. We will address what happens if this condition is violated later.
The Gutzwiller projection can be conveniently implemented if one works with the real-space basis states, which are already in the projected Hilbert space with no double occupancy. For example, consider that the spin up and down electrons are located in a given configuration set of {r ↑ , r ↓ }. Then, the amplitude of the RVB state for this configuration set can be given as follows: where i and j run though the coordinates of all spin up and down electrons, respectively [3]. Now, let us come back to what happens if u k = 0 at some momenta. In this situation, which can occur for the d-wave pairing symmetry, the amplitude of the RVB state takes a rather complicated form. Specifically, one has to take into account the contributions from u k = 0 separately from those from otherwise. For example, let us say that u k = 0 at k = k. Then, the plane wave states at k and −k are "preoccupied" by the spin up and down electrons, respectively, rather than forming Cooper pairs.
Consequently, the amplitude of the RVB state should be modified as follows: × det (e ik l ·r m↑ ) det (e −ik l ·r n↓ ), (35) where k l is the l-th preoccupied momentum, and r m↑ and r n↓ are the coordinates of the m-th spin-up and n-th spin-down electrons, respectively, belonging to the preoccupied plane wave states. The sum is taken over all possible combinations of choosing {r ↑ , r ↓ } out of {r ↑ , r ↓ }. p is the permutation parity occurring when all creation operators are rearranged in a predetermined convention.
Unfortunately, computing Eq. (35) turns out to be quite time-consuming since there are simply too many different combinations. To reduce the computing time, we employ a trick of adding a very small constant to the pairing amplitude, i.e., ∆ k → ∆ k + δ, which eliminates the preoccupied momenta, leaving only the single combination, where all electrons are Cooper-paired. At the end of computation, we take the limit of vanishing δ. In Figs. 4 and 5, we take δ/t = 0.0001.