{'fraction_non_alphanumeric': 0.04223031093243512, 'fraction_numerical': 0.009634016044620705, 'mean_word_length': 4.100868280066506, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Grand canonical Gutzwiller approximation for magnetic inhomogeneous systems

The Gutzwiller approximation (GA) for Gutzwiller-projected grand canonical wave functions with fugacity factors is investigated in detail. Our systems in general contain inhomogeneity and local magnetic moments. In deriving renormalization formulae, we also derive or estimate terms of higher powers of intersite contractions neglected in the conventional GA. We examine several different constraints, i.e., local/global spin-dependent/independent particle-number conservation. Out of the four, the local spin-dependent constraint seems the most promising at present. An improved GA derived from it agrees with the variational Monte Carlo method better than the conventional GA does. The corrections to the conventional GA can be interpreted as two-site correlation including the phase difference of configurations. Furthermore, projected quasi-particle excited states are orthogonal to each other within the GA. Using these states, spectral weights are calculated. We show that asymmetry between electron addition and removal spectra can appear by taking into account the higher powers of the intersite contractions in the case of the d-wave superconductors and the Fermi sea; the addition is smaller than the removal. However, the asymmetry is quite weak especially near the Fermi level. In contrast, projected s-wave superconductors can have the opposite asymmetry (addition larger than removal) especially near the Fermi level. In addition, formulae from the other three constraints are also derived, which may be useful depending on purposes.


I. INTRODUCTION
This paper concerns calculation of expectation values using projected wave functions in inhomogeneous systems. In order to study electronic systems with repulsive on-site interactions, Gutzwiller proposed projected wave functions 1 of the form P G |Ψ 0 with the Gutzwiller projection operator, to prohibit electron double occupancy on each site. Here, n iσ = c † iσ c iσ with c † iσ (c iσ ) being the creation (annihilation) operator of site i and spin σ.
Expectation values of operators by this projected wave function can be evaluated by the variational Monte Carlo method (VMC) numerically exactly within statistical errors. However, the VMC requires lots of computational effort for some issues. In addition, it needs one run for each parameter set, whereas an analytical method can generate more general formulae that often provide us some hint to understand the system. Thus, instead of the VMC, an analytical approximation called the Gutzwiller approximation (GA) is used on occasions, i.e., with |Ψ N ≡ P G |Ψ N 0 , where |Ψ N 0 have a fixed particle number N . The factor g O is the Gutzwiller renormalization factor for the operatorÔ. If one chooses a noninteracting or mean-field approximated wave function as |Ψ N 0 , the expectation value in the r.h.s. of Eq. (2) can be easily evaluated. The renormalization factor for the hopping term denoted by g t is smaller than unity because it is more difficult to hop in the presence of the strong on-site Coulomb repulsion between electrons; that for the exchange interaction denoted by g s is larger than unity because each site is more often singly occupied to avoid the other electrons. The GA was first introduced for the Hubbard model by Gutzwiller 2,3 , then reformulated by Ogawa et al. 4 A clear description of the method has been given by Vollhardt 5 . It was also applied to a mean-field theory for the t-J model by Zhang et al. 6 Improvements of the GA by taking more intersite correlations have been made by several authors. 7,8,9 The GA usually produces qualitatively correct results although it is reported that there are also qualitative differences in some cases 10 .
The original formulation of the GA implicitly assumes that a wave function before the projection has a fixed particle number N (in the following, we call it the "canonical scheme"). If the particle number of a wave function has fluctuation (the"grand canonical scheme"), then the Gutzwiller projection reduces the particle number (see Appendix A). Such reduction of the particle number may arouse a question whether the GA as Eq. (2) is valid because this equation seems to say that the wave functions before and after the projection have similar properties except for the double occupancy; are they similar if they have different particle numbers? To avoid such an unclear path, Anderson and Ong 11 , and Edegger et al. 12 formulated a grand canonical GA by taking the canonical scheme as a guide. Namely, one can force the projection not to change the average particle number, by gluing to P G a fugacity factor that compensates the particlenumber reduction. To our knowledge, the fugacity factor was first seen in a preliminary form in the paper by Yokoyama and Shiba 13 to relate the canonical and the grand canonical VMC. Gebhard 14 introduced positionand spin-dependent fugacity factors for calculational con-venience of the 1/d expansion whose d → ∞ limit corresponds to the GA. They also appear in the construction of the gossamer superconductivity by Laughlin 15 . Then, Wang et al. 16 used position-dependent but spinindependent fugacity factors for inhomogeneous systems.
The fugacity factors allow us freedom to choose a relation between the particle numbers before and after the projection, and the renormalization depends on this choice. Recently, Ko et al. 17 pointed out that two contradictory formulae of the Gutzwiller renormalization factors in the literature actually come from two different choices of the fugacity factors. That is, (i) the fugacity factors are determined so that the projection conserves the local particle density of each spin direction at each site, or (ii) so that the projection conserves the total particle number for each spin direction (this is the usual canonical-scheme constraint). Mainly for the square lattice antiferromagnet, they used the canonical scheme, and introduced additional position-and spin-dependent fugacity factors, then calculated each renormalization factor as a ratio of probabilities for the physical process.
In this paper, we examine in detail several different choices of fugacity factors that impose local/global spin-dependent/independent particle-number conservation. We adopt the grand canonical scheme, and derive general formulae. Some of our formulae are different from those by the canonical derivation. Furthermore, corrections to the conventional GA are also estimated or derived by taking intersite correlations into account. The structure of the paper is as follows: Secs. II and III are devoted for the case (i), and Sec. IV for (ii). First in Sec. II, we derive renormalization of the hopping and the pairing amplitude, the local spin moments and the exchange interaction from the local spin-dependent constraint. We test the formulae of the hopping amplitude by comparing with the VMC. Physical interpretations are given for newly derived terms. Subsequently in Sec. III, we also check orthogonality and excitation energies of projected Bogoljubov quasiparticle states, and discuss asymmetry between positive and negative bias spectra. Next, in Sec. IV, formulae from the global spin-dependent constraint are derived. The formulation there includes cases where the particle numbers before and after the projection are unequal. In addition, grand canonical GAs with local/global spin-independent constraints are briefly discussed in Sec. V.
In our impression, the grand canonical scheme simplifies calculation in many cases because it is free from complicated configuration counting. Furthermore, systematic improvement is straightforward by including terms from larger clusters in the linked-cluster expansion. 1 The formulation we use is similar to the 1/d expansion by Metzner and Vollhardt 18 , and Gebhard 14 . The lowestorder theory in the uniform non-superconducting limit of our formulation for the case (i) is equivalent to d → ∞ limit of the 1/d expansion. However, in inhomogeneous systems and in the presence of the second and the third neighbor hopping, it is not clear if 1/d is a good ex-pansion parameter. In addition, considering future improvements of the theory, it may be difficult to define terms of very high order in 1/d. Therefore, we naively use the linked-cluster expansion as Gutzwiller's original formulation 1 , then expand it in a power series of intersite contractions and neglect high order terms. Furthermore, we do not adhere to making derived formulae into the form of Eq. (2).
Throughout this paper, we use the following notation: A wave function before a projection is denoted by |Ψ 0 and it does not have a definite particle number and may have some inhomogeneity in general. Then, the wave function after the projection is represented by |Ψ = P |Ψ 0 , where P is a generalized projector that includes fugacity factors defined later. The expectation values of an arbitrary operatorÔ by these wave functions are denoted by Furthermore, In addition, S i denotes the spin operator at site i.

II. LOCAL CONSTRAINT
The Gutzwiller projection changes electron-density distribution in inhomogeneous systems in general. However, by introducing fugacity factors, one can force desired electron-density distribution. We prefer to start from the grand canonical GA with a local constraint for each spin direction, namely, for any i and σ. Note that this local constraint is different from the canonical scheme constraint that conserves the total particle number. However, this "local canonical" constraint simplifies the resultant formulae as shown in the following. For example, some of low order corrections to the GA vanish automatically. Furthermore, with this constraint, projected Bogoljubov quasiparticle states are approximately orthogonal to each other, and excitation energies are approximatively obtained by diagonalizing a renormalized Hamiltonian (shown in Sec. III). In general, S x i 0 and S y i 0 can be finite. Such cases will be discussed only in Sec. II D, and otherwise S x i 0 = S y i 0 = 0 and c † iσ c jσ 0 = 0 are assumed. Furthermore, although we have d-wave superconductors in mind, there may be deviation from d-wave in inhomogeneous magnetic systems, and c † i↑ c † i↓ 0 (on-site pairing before the projection) can be non-zero. We discuss effect of c † i↑ c † i↓ 0 = 0 in Sec. III G, and otherwise assume c † i↑ c † i↓ 0 = 0 for any i. We also do not consider triplet pairing of the form c † iσ c † jσ 0 and set it to zero for any i, j, σ. The generalization to c † iσ c † jσ 0 = 0 is straightforward.

A. Condition for fugacity factors
The projected wave function is defined as |Ψ = P |Ψ 0 with P ≡ i P i , where The local up and down particle numbers are controlled by λ 1 2n iσ iσ , and the fugacity factors λ iσ will be determined later to satisfy Eq. (6). In order to derive their explicit forms, let us calculate the density of σ-spin electron at site i, In principle, by applying the Wick theorem, these expectation values can be exactly evaluated. In practice, however, such calculation is quite difficult to carry out because too many terms appear by the Wick decomposition.
To approximate it, remember that intersite contractions, n ijσ and ∆ ij , are much smaller than on-site contractions, n iσ . An approximation to take the leading order with respect to the intersite contractions corresponds to the GA. Here, we take only on-site contractions. Then, l = i terms cancel out between the numerator and the denominator, namely, Therefore, the condition to determine λ iσ is given by λ iσ (1 − n iσ )/Ξ i = 1. By solving the simultaneous equations for up and down spins, we obtain The corrections to P 2 0 and n iσ can be calculated by taking into account intersite contractions between site i and other sites l = i. Let us calculate terms proportional to |n il↑ | 2 . Such terms appear by the Wick decomposition of n i↑ P 2 0 . We take on-site contractions for the sites other than i, l, and thus we only need to consider n i↑ P 2 l 0 . The operators inn i↑ = c † i↑ c i↑ are contracted with those in c † l↑ c l↑ or c l↑ c † l↑ in P 2 l . Then, the operators for the down spin are replaced by n l↓ or 1 − n l↓ . Namely, such contribution is written as In other words, the terms proportional to |n il↑ | 2 vanish when λ iσ is set as Eq. (11). Similarly, terms proportional to ∆ 2 il also vanish. Therefore, with Eq. (11), we have n iσ = n iσ + O(n 4 ijσ ) + O(∆ 4 ij ). Estimated corrections to λ iσ are also of the order of n 4 ijσ or ∆ 4 ij .

B. Hopping and pairing amplitude
For the hopping term, similar calculation can be carried out. Namely, for i = j, where we took on-site contractions except one intersite contraction (that is necessary) in the numerator. Then, the Gutzwiller renormalization factor is given by 17 The next order in fact involves one more site other than i and j, but the second and the third order of the intersite contractions for such contribution vanish when λ iσ is set as Eq. (11). Therefore, the third order term involves only sites i and j. Namely, taking more contractions between i and j in Eq. (12), .
The formula for c † i↓ c j↓ is obtained by replacing as ↑⇔↓ and ∆ ij ⇒ −∆ ji . The n ij↑ |n ij↓ | 2 term in Eq. (15) is from repulsive correlation between down-spin holes due to the Pauli principle: all of the four configurations in Fig. 1 contribute to c † i↑ c j↑ 0 , but only (a) does to c † i↑ c j↑ . Then, taking into account repulsion between down-spin holes, (a) has less weight than the estimate by the conventional GA that neglects this correlation.
On the other hand, the n ij↓ ∆ * ij ∆ ji term is from superconducting correlation; negative for ∆ ij = ∆ ji (singlet), and positive for ∆ ij = −∆ ji (triplet). This term seems related to the phase difference between the four configurations in Fig. 2, which appear in |Ψ 0 (before the projection). Our rough explanation in the case of n ij↑ ≃ n ij↓ is as follows: Suppose ζ a , ζ b , ζ c , ζ d are coefficients of the configuration (a,b,c,d) in |Ψ 0 , and assume they are real numbers. Then, ζ a ζ d contributes to n ij↑ , and −ζ b ζ c contributes to n ij↓ . Remember that the conventional GA can be derived by taking the ratio of the probability of configurations 6,17 ; it implicitly assumes that ζ a ζ d and −ζ b ζ c have the same sign. Turning on the superconducting correlation, configurations (a) and (b) as well as (c) and (d) start to correlate. Then, their contribution to the singlet order parameter before the pro- The magnitude of this quantity, however, is small if ζ a ζ d and −ζ b ζ c have the same sign. Therefore, to strengthen the singlet superconducting correlation, all of ζ a ζ d , ζ b ζ c should be small. Accordingly, the weight of Fig. 1(a) should be smaller than the estimate by the conventional GA.  We test this formula for a simplest case, i.e., the standard uniform non-magnetic BCS d-wave superconductor, The conventional GA as Eq. (14), the generalized GA as Eq. (15), and the VMC are compared in Fig. 3 for the nearest-neighbor hopping. Here, the generalized GA is done using 200×200 sites, and practically the finite-size effects are negligible; errors only come from neglect of the higher order of the intersite contractions. µ is adjusted to satisfy each hole concentration for each point. The VMC is carried out using 30×30 sites with x-antiperiodic y-periodic boundary condition. The hopping amplitude is averaged over every bond, and the statistical errors are negligible in the scale of this figure. For comparison with the GAs, µ is also adjusted to equalize the doping before the projection with that after. At small ∆ v , the generalized GA agrees with the VMC very well. As ∆ v increases, the deviation becomes larger. This is possibly because O(∆ 4 ij ) term neglected in Eq. (9) may start to make an important contribution. In inhomogeneous systems, there may be deviation from the d-wave. Since it is rather difficult to force the local constraint for the VMC in inhomogeneous systems, let us test the formula using a simpler non-d-wave, namely, a uniform p-wave superconductor by redefining in Eq. (16). The nearest-neighbor hopping amplitude in the x-direction is plotted in Fig. 4. The generalized GA shows a good overall agreement with the VMC. It especially reproduces characteristic peak at ∆ v ∼ 2t caused by the n ij↓ ∆ * ij ∆ ji term in contrast to the conventional GA.
The superconducting order parameters c † i↑ c † j↓ can be calculated similarly to the hopping term, i.e., . . The ∆ * ij |∆ ji | 2 term contains the attractive correlation between holes of i ↓ and j ↑; if ∆ ji is finite, i ↓ and j ↑ tend to be simultaneously occupied or unoccupied, and it is less likely that only one of them is occupied. Accordingly, this effect increases weight of the configurations in Fig. 5(a) and (d), and appears as the positive correction in Eq. (18). The ∆ * ij n ij↑ n * ij↓ term represents roundabout correlation between i ↑ and j ↓ through i ↓ and j ↑ as depicted in Note that Eqs. (15,18) is mainly aimed at |i − j| = 1. For next-nearest neighbors, O(n 4 ij ) of |i − j| = 1 may be comparable to O(n 2 i ′ j ′ ) of |i ′ −j ′ | = 2 and the former may be dominant especially in high dimensions. In general, as i and j separate from each other, the approximation by Eqs. (15,18) may lose accuracy.

C. Spin moment and exchange interaction
By definition, the local spin-z component at each site is not renormalized, i.e., For the exchange interaction term S i · S j , we take up to the second order of intersite contractions. Using symbols ↑, ↓, and +1, −1, interchangeably, it is written as with .
Here, Eq. (20) seems different from what is derived as a ratio of probabilities for the physical process using the canonical scheme with the fugacity factors by Ko et al. 17 We speculate that it possibly does not take into account all of the contractions above.
To compare with the result by 1/d expansion by Gebhard 14 , set ∆ ij = ∆ ji = 0 and consider antiferromagnets. By setting n iσ = n jσ in Eqs. (20,21), these equations are reduced to which are equivalent to the formula by the 1/d expansion 14 . However, when ∆ ij = 0, renormalization of S z i S z j is not reduced to such a simple form, and we need the original formula, Eq. (20). Note that Eqs. (23,24) can be used also for non-superconducting ferromagnets. Namely, the local constraint leads to the conclusion that antiferromagnets and ferromagnets are renormalized similarly. This is in distinct contrast to results of the GA with global constraint as will be discussed in Sec. IV F.

D. Systems with nonzero spin-xy components
This choice of fugacity factors encounters difficulties when S x i 0 or S y i 0 is finite. Let us redo the derivation including S ± i ≡ S ± i 0 : This is solved to give For a spin moment, S z i = S z i 0 , and This renormalization factor for S ± i is larger than unity because it is not bound by the local constraint. Since xy component is renormalized differently from z component, approximation depends on humans' choice of z-axis. This asymmetry is probably related to what is discussed by Ko et al. 17 The most reasonable choice of z-axis we think is making it parallel to S i 0 at each site. Then, S ± i = 0 for any i. It is equivalent to formulating a GA with constraints n i↑ +n i↓ = n i and S i = S i 0 . However, such a GA may yield very complicated renormalization factors for intersite terms. One way to avoid such a complexity is to use spin-independent constraint as shown in Sec. V A.

III. LOCAL CONSTRAINT: EXCITED STATES
The GA with the position-and spin-dependent constraint discussed in the previous section has an advantage in constructing plausible excited states which are approximately orthogonal to each other as shown below.
For shorthand notation, we use where N L is the number of lattice sites. Then, the subscript of this new operator runs from 1 to 2N L , and we represent it by single Greek symbols as c ρ . Furthermore, we definen A. Bogoljubov de Gennes (BdG) equation As a preparation, let us begin with deriving a BdG equation by minimizing the Gutzwiller-approximated energy following the procedure by Wang et al. 16 In the following, we work more on general properties of a BdG equation with the Gutzwiller projection, and do not use any Hamiltonian explicitly. However, what we have in mind is inhomogeneous t-J-type models, where the t ij term with i = j may represent local impurity potentials. The zero-temperature grand potential Ω ≡ H tJ − µ i,σn iσ can be approximated by the GA, and represented by a function of n ρζ , namely, We do not show the explicit form of Ω GA because it can be derived straightforwardly by using the formulae in the previous section. In the derivation, one can choose the level of the approximation: If one takes only the leading order of the intersite contractions, formulae in the nonmagnetic case are equivalent to those derived by Wang et al. 16 and Li et al. 19 If an improved Gutzwiller approximation such as Eq. (15) is used, a more accurate solution can be obtained in principle, although it may be more difficult to find self-consistent solutions. The chemical potential µ is determined to adjust the particle number N to satisfy N = −∂Ω GA /∂µ. The other variables are functional of Ψ 0 and determined by minimizing Ω GA , Assuming Ψ 0 |Ψ 0 = 1, then By combining Eqs. (35,36), where Then, Ω GA takes an extremum when |Ψ 0 is an eigenstate of H BdG , namely, The main differences from usual BdG Hamiltonian are the local renormalization factors in front of t ij and J ij , and the effective local chemical potential terms − i µ iσniσ with µ iσ = −∂Ω GA /∂n iσ − µ which come from n iσ -dependence of the renormalization factors. Local modulations of t ij and J ij tend to be enhanced by the local renormalization factors, and impurity potentials tend to be screened by the local chemical potentials. 20

B. Quasi-particles
We rewrite H BdG in a matrix form, The 2N L × 2N L matrix H ρζ can be diagonalized using a unitary matrix U , namely, Then, using the Hamiltonian is diagonalized as, The ground state of this effective Hamiltonian is |Ψ 0 ≡ n γ n |0 . Suppose the ground state is well approximated by P |Ψ 0 . Naively, one may assume that excited states are constructed by P γ † n |Ψ 0 . This form of excited states was first introduced for uniform systems by Zhang et al. 6 For fugacity factors in P , we use those in the ground state even for the excited states. It probably correspond to assuming that the quasi-particles γ n are not very localized and that the change of the particle distribution is negligible.

C. Orthogonality of the excited states
The orthogonality of these excited states can be checked by expanding γ n using Eq. (43). For example, for E n > 0, E m > 0, Here, we have to mind a discrepancy between creation and annihilation operators; Then, as the leading-order theory, we take on-site contractions except one intersite contraction. Thanks to Eq. (11), renormalization factors are reduced to a simple form, i.e., and we obtain simple results, does not satisfy this relation, but we are lucky enough to use the off-site pairing assumption, c † iσ c † iσ 0 = 0, then c † iσ P 2 c † iσ 0 ≈ 0, and can exclude such an exception. Our projected superconducting state includes the Fermi sea as a special case, and these relations look different at a sight from those derived for the Fermi sea by Fukushima et al. 21 In fact, however, these are identical if one remembers that P contains fugacity factors.
Using Eq. (49-52), one can transform back from c ρ to γ n to yield That is, the excited states are orthogonal to each other within the GA.

D. Excitation energy
Let |0 and |n denote the normalized ground and excited states, Neglecting the second order of the difference in n ρζ , Therefore, the excitation energies are approximately the same as eigenenergies of the effective Hamiltonian.

E. Density of states
To calculate the local density of states, we need matrix elements, | n|c † iσ |0 | 2 and | n|c iσ |0 | 2 . First of all, using Eq. (53) with n = m, the normalization of the excited states can be replaced as Then, we expand γ n in |n using Eq. (43) and use simple relations similar to Eqs. (49-52), namely, for any i, j, σ, τ . These formulae are true also for i = j, more explicitly, These relations are exact if exact λ iσ are used.
Since the indices of the renormalization factor are only from those of the operator between two P 's, we can transform back to γ n to yield where ρ = (i, σ) as Eq. (31). The common renormalization factor g t0 iiσ tells us that the positive and negative bias spectra are symmetric. This symmetric density of states is also obtained by the canonical scheme GA 21,22 . We go one-step further about this point in the next subsection.
F. Electron addition-removal asymmetry caused by higher order terms The conventional BCS theory tells us that the quasiparticle excitation spectra are symmetric between positive and negative bias. However, local density of states of high-T c superconductors measured by the STM is highly asymmetric and there is an argument that attributes this asymmetry to strong electron correlation. 11 Namely, electron addition may be more difficult than electron removal because the injected electron may be repelled by the other electrons due to their strong Coulomb repulsion. It is controversial whether the projected quasi-particle states have symmetric spectra or not. The GA gives symmetric spectra 21,22 if only quasi-particle excitation is considered (incoherent excitations may cause asymmetry 22 ). In contrast, the spectra calculated by the VMC show asymmetry. 10 To discuss this point, here we calculate corrections to the results in the former sections. When these corrections are taken into account, the orthogonal relation, Eq. (53), may not be satisfied any more. Therefore, in the following, we assume that the systems are almost uniform; in the uniform limit the wave number is a good quantum number due to the translational symmetry, and thus excited states are orthogonal. The next order corrections contain only site i and j similarly to those in the hopping term. We put general formulae in Appendix B, and here only show a special case of n i↑ = n i↓ = n i /2, n ij↑ = n ji↑ = n ij↓ = n ji↓ ≡ n ij , ∆ ij = ∆ ji = ∆ * ij = ∆ * ji . Then, with Eqs. (57-60) are rewritten as Since A ij ≥ 0 and α i ≥ 0, the corrections are positive for the electron removal, and negative for the addition. For more careful analysis, we also need to check the normalization. Including the corrections, Eqs. (49-52) are rewritten as These corrections to the normalization are all negative, and they do not seem to cancel the asymmetry in Eqs. (67-70). Therefore, these results suggest that the higher-order GA exhibits asymmetric spectra whose electron addition spectra are smaller than the removal. This asymmetry is consistent with the VMC calculations for excitation spectra by Chou et al. 10 , and for spectral weights by Bieri and Ivanov 23 and Yang et al. 24 For more explicit comparison, we calculate the spectral weights, and show them in Fig. 6 for both the conventional and the generalized GA. Here, we include t ′ and t ′′ in addition to Eq. (16) for a better correspondence to the high-T c superconductors. In the case of the standard BCS theory, Z + = |u k | 2 , Z − = |v k | 2 . Then, for each k-point below the Fermi level, one can find a corresponding point k ′ above the Fermi level such that E k ′ = E k , u k ′ = v k , v k ′ = u k . Then, summation of the contribution from these two points to the spectra is unity for both addition and removal spectra 25 because |u k | 2 + |u k ′ | 2 = |v k | 2 + |v k ′ | 2 = |u k | 2 + |v k | 2 = 1. Accordingly the excitation spectra are symmetric. The results of the conventional GA are Z + = g t |u k | 2 , Z − = g t |v k | 2 ; namely, the spectra are just renormalized by g t , and are symmetric as the standard BCS theory. In contrast, by including the corrections to them, Z − decreases and Z + increases, which can cause the asymmetry in the spectra. These Z ± are consistent with the VMC results 23,24 . Note that A ij is finite even for ∆ ij = 0, i.e., the Fermi sea also has the asymmetry, which is also consistent 24 . Similarly to the hopping term, Eqs. (67-74) are more accurate for small |i − j|. Hence, the Fourier transformed results may include errors from the summation over large |i − j|. It will be checked in the future studies by including higher order terms. Since this asymmetry appears as a deviation from the conventional GA, it is rather small (especially near the Fermi level), and does not look like what is seen in the STM experiment.

G. Opposite asymmetry in projected s-wave superconductors
We speculate that the origin of the asymmetry may not be so simple as the intuition that electron addition may be more difficult than removal because electrons repel each other. Here, we show a counterexample against this simple scenario. That is to say, projected s-wave superconductors can have the opposite asymmetry; the electron addition spectra are larger than the removal. Such projected s-wave superconductors may be realized if the pairing interaction is isotropic because d-wave does not gain energy from diagonal J ij . Even if J ij is finite only for nearest neighbors, the mean-field approximation in very overdoped systems converges to extended s-wave solutions. To be more precise, this opposite asymmetry is related to finite on-site pairing before the projection, c † i↑ c † i↓ 0 = 0 and not really related to the symmetry of the gap. Then, even for d-wave, inhomogeneity causes deviation from the d-wave, and c † i↑ c † i↓ 0 can be nonzero in general. Therefore, strongly disordered d-wave superconductors could have similar properties.
To take ∆ ii = 0 into account, we have to redo the derivation from the beginning. Then, Ξ i and λ i should be replaced by In fact, this generalization makes analytical treatment very difficult, and in the following we take only the leading order of the intersite contractions. Accordingly, its ∆ ii → 0 limit corresponds to the conventional GA (not the generalized GA in Sec. III F).
The most important matrix elements are Here, Eq. (80) is exact because of P G c i↑ P G = c i↑ P G and P G c † i↑ c † i↓ = 0. The renormalization factor in Eq. (79) is obviously larger than that in Eq. (80). Hence, these matrix elements suggest that the asymmetry is the opposite to that of d-wave. We have also calculated other matrix elements, and after the Fourier transform, Z ± for a uniform system are obtained as plotted in Fig. 7. In comparing Figs. 6 and 7, we should keep in mind that a better approximation is used for Fig. 6, and the asymmetry as in Fig. 6 does not appear in Fig. 7. Nevertheless, the characteristic asymmetry near the Fermi surface in Fig. 7 is very strong, and may remain even in calculation with a better precision. Most likely, Eqs. (79,80) mainly contribute to the asymmetry because the vicinity of the Fermi level changes most dramatically.

H. Physical consideration for the asymmetries
For the projected s-wave superconductors, the physical origin of the asymmetry may be understood as follows. We have been using terms "addition" and "removal" but these are in fact named from the ground state's view. If one takes the complex conjugate, this addition (removal) matrix elements can be regarded as removal from (addition to) an excited state. Let us adopt the excited states' view for a while. In the s-wave BCS superconducting state (before the projection), a Cooper pair may be formed more or less on-site, which is a resonance of the doubly occupied state and the empty state. When c iσ is operated to this wave function, it chooses the the "originally doubly occupied" state. Then, in P c iσ |Ψ 0 , the opposite spin state, iσ, is occupied with high probability. Accordingly, it is and easy to remove iσ electron. In contrast, for P c † iσ |Ψ 0 , it is impossible to add iσ electron. Finally, let us turn back to the ground state's view and review the arguments above. Then, the removal is difficult, but the addition is easy.
The asymmetry in the d-wave superconductors and the Fermi sea needs more consideration because it appears by higher order correlations. Fig. 8(a) shows a configuration in the ket |Ψ 0 contributing to c j↑ P c † i↑ P 0 . The first term in Eq. (68) represents direct correlation between c j↑ and c † i↑ . The second term comes from the repulsive correlation between down holes, which reduces weight of this configuration. On the other hand, for c † i↑ P c j↑ P 0 , both configurations (a) and (b) in |Ψ 0 contribute. However, when electron density is high, (b) is dominant because empty sites are rare. Then, correlation between down holes increases the weight of (b). Since this effect appears only at high density, the second term in Eq. (67) accompanies the factor α i .

IV. GLOBAL CONSTRAINT (QUASI-CANONICAL GA)
In the former sections, we have used the local constraint. However, if one needs the GA as an approximate method of the VMC, it may be preferred to require the usual canonical constraint, i.e., the total particle number constraint for each of up and down spins, where N after σ is the total number of σ electrons after the projection. Although one takes N after σ = i n iσ in the usual canonical GA, the particle numbers before and after the projection can be different in general. In fact, in the VMC, they are different; the particle number projection P N after σ is usually applied together with P G , and the chemical potential of |Ψ 0 is more like a variational parameter and does not control the particle number after the projection. In the following we use notation, with N L the total number of sites. Our purpose here is formulating a grand canonical GA that gives results of the canonical scheme, by imposing Eq. (81). If the total spin moment is nonzero, it must be reasonable to choose the spin-z axis parallel to the global moment so that i S x i 0 = i S y i 0 = 0. In that case, local xy components, S x i 0 , S y i 0 , may be finite in general. Then, similarly to the local-constraint formulation in Sec. II D, xy components of the local spin moments are renormalized differently from their z components. The canonical scheme condition for the total spin moment restricts the spin-z renormalization factor to the vicinity of unity, but the other directions are free from it. We expect that this spin-rotational asymmetric renormalization is a property from the canonical condition and should exist even in exact calculation. Furthermore, if the total spin moment is zero and local moments point to various directions, we have no idea how to choose the z-axis. Here, to avoid such complexity, we assume S x i 0 = S y i 0 = 0, and as in the former sections

A. Condition for fugacity factors
To control the total particle numbers, we need a factor in the form λ 1 2 P in iσ σ , namely, the fugacity factors λ σ do not have the site index. Accordingly, the projected wave function is defined as |Ψ = P |Ψ 0 with The formula for n iσ has the same form as Eq. (9), and only λ σ is different, i.e., Note that Ξ i is still site-dependent because it contains local electron densities. By inserting it into Eq. (81), we obtain In inhomogeneous systems, λ σ is solved numerically from Eq. (84) in general. An important point of this uniform fugacity approach is that the local electron density is also renormalized as in Eq. (83), and n iσ = n iσ 0 in general. When λ σ is solved and inserted into Eq. (83), the corrections to n iσ are of the second order of intersite contractions as will be explicitly shown in Sec. IV D.
The local spin-z component is renormalized as where symbols ↑, ↓ and +1, −1 are interchangeably used.

B. Hopping term
The Gutzwiller renormalization factor of the hopping term is given by Next order corrections to this formula involves another site, which may make important contribution for second or third neighbor hopping in some systems. Using a iσ ≡ (1 − λσ)(1 − n iσ ) + λ σ n iσ , it is written as (87) The corrections to λ iσ and Ξ i affect only from third order and not relevant to the equation above. Note that a iσ goes to zero in the uniform limit with n after σ = n σ . In the uniform systems, by omitting irrelevant site indices and using R σ ≡ n after σ /n σ , we obtain C. Exchange term with zero total spin-z component and λ ↑ = λ ↓ The general formulae for the exchange interaction term are too lengthy to present here. Our main interest is in systems with zero total spin-z component and λ ↑ = λ ↓ which includes such as non-magnetic systems, antiferromagnets, and stripes. Hence, for simplicity, we restrict ourselves to this case in the following except for Sec. IV F that treats ferromagnetic systems. The generalization to nonzero total spin-z component is straightforward but one has to work with more complexities.
When λ ↑ = λ ↓ ≡ λ, Eq. (85) is reduced to For the exchange interaction term S i · S j , we take up to the second order of intersite contractions. Assuming that m i is small, namely, where g s ij ≡ λ 2 (Ξ i Ξ j ) −1 . This is the result of the conventional GA. However, note that, if m i is of the order of unity, terms such as n 2 ijσ m i , ∆ 2 ij m i , have about the same order of contribution as n 2 ijσ , ∆ 2 ij , and formulae above should be modified as derived below. D. Beyond the "conventional" GA When m i ∼ n iσ , terms neglected in the previous derivation may grow, and we need to redo the derivation from the beginning. It is known that expectation values by Gutzwiller-projected states can be written in the form of a linked cluster expansion. 1 The terms we need here includes contribution from clusters one-site larger than those in the previous derivation.
We relegate detailed derivation to Appendix C, and only show final results here. The renormalization of the particle densities is given by where a lσ = (1 − λ)(1 − n lσ ) + λn lσ . The formula for n iσ is obtained by replacing as ↑⇔↓ and ∆ il ⇒ −∆ li . Then, the new equation to determine λ is given by i n iσ = N after . The solution can be written as λ ≈ λ (0) + λ (2) , where λ (0) is λ determined by Eq. (84), and λ (2) is the correction to it represented by Here, every λ is replaced by λ (0) in the r.h.s. Using this new λ, we can calculate spin terms, Here, λ in the second order terms can be replaced by λ (0) . Note that the contribution that involves a third site l in S z i cancels that in S z i S z j by the subtraction of S z i S z j . There is no correction of this order for S x i S x j + S y i S y j in Eq. (92).
Although several authors 7,8,9 formulated improved canonical GAs by taking nearest-neighbor correlations similarly to ours, our result is different from any of them even if we neglect the second-and the third-neighbor terms. The origin of this discrepancy is not clear at present.

E. Antiferromagnets
As an explicit example, we show the formulae for the square lattice antiferromagnet. For periodic systems, we can restrict the summation over the site index i to only inside of the unit cell. In the presence of the antiferromagnetic moments, n ij between second-or third-neighbor pairs may be comparable to or larger than that of the nearest neighbor pairs. To take into account these terms, we define n ijσ = χ, χ ′ , χ ′′ , ∆ ij = ∆ ji = ∆, ∆ ′ , ∆ ′′ , for the nearest, the second, the third neighbor pairs, and assume these are real numbers. In addition, n i↑ = n A and n i↓ = n B for A-sublattice, n i↑ = n A and n i↓ = n A for B-sublattice, and then m = (n A − n B )/2. By omitting irrelevant site indices, where Eq.(103) is for a nearest neighbor pair. In general, ∆ ij = ∆ ji may occur; in that case these equations need modification with a little more complexities.

F. Ferromagnets
Here, we show show a remarkable difference from the local constraint. That is, ferromagnets are renormalized very differently from antiferromagnets in contrast to results by the local constraint in Sec. II C. For ferromagnetic wave functions without superconductivity, we can set ∆ ij = 0, and N after σ = i n iσ 0 . Then, in the uniform cases, this GA with the global constraint is equivalent to the one with the local-constraint. Therefore, by setting m i = m in formulae in Sec. II, we obtain those for the ferromagnets, Many of the formulae derived with the global constraint in this section contain a σ , but it goes to zero in this ferromagnetic limit. It is consistent with no appearance of a σ in the local-constraint formulation.
In fact, the renormalization for the spin-z component represented by Eq. (104) is different from the one derived by Zhang et al. 6 using a probability argument of the canonical GA, namely, S z i S z j = g s S z i S z j 0 with g s defined by Eq. (106). However, we speculate that our result is more reasonable because spin moment term m 2 is not renormalized; the canonical constraint prevents the spin-z component from growing larger, in contrast to the antiferromagnetic moments, which are not bound by the canonical constraint. It may be clearer if we take the limit of small m for Eq. (104), and compare with Eq. (91) for antiferromagnets.

G. Effect of N after = N before
Projections reduce the Hilbert space. Hence, many wave functions may be equivalent to each other after the projection even if they are different before the projection. Here, we demonstrate it explicitly by the particle num-ber projection. Let us start from uniform non-magnetic cases. Define two BCS states, whereN ≡ iσn iσ . Under the particle number projection P N , Namely, wave functions P N |Ψ 0 and P N |Ψ ′ 0 are equivalent whereas |Ψ 0 and |Ψ ′ 0 are nonequivalent. At a sight, the quasi-particle excited states of these two BCS states look different because γ † kσ does not commute withλN /2 . However, c † kσ and c k −σ in γ † kσ in fact yield the same state, and thus P N γ † kσ |Ψ 0 and P N γ ′ † kσ |Ψ ′ 0 are equivalent, where γ ′ † kσ is a quasi-particle operator for |Ψ ′ 0 . Therefore, even if the average particle number of |Ψ 0 is not N after , one can make that of |Ψ ′ 0 equal to N after by choosingλ to satisfy then, using |Ψ ′ 0 , the GA can be applied with the local constraint n iσ = n iσ 0 . Accordingly, we can use convenient properties derived in Secs. II and III.
Such a transformation to relate the global constraint to the local one may be possible also for inhomogeneous systems, but there seems to be a problem. The particle numbers can be controlled by fugacity factors iσλn iσ /2 iσ as Eq. (109). One can chooseλ iσ in |Ψ ′ 0 to satisfy n iσ = n iσ 0 . Then, c † iσ is replaced byλ 1 2 iσ c † iσ as well as c iσ is replaced byλ − 1 2 iσ c iσ (in annihilating the electron, the fugacity factor caused by the creation should be canceled). Accordingly, |Ψ ′ 0 is a rather strange wave function because the quasi-particles may not satisfy the fermion commutation relation. Then one may need to redefine a proper quasi-particle set for |Ψ ′ 0 . Furthermore, since the quasi-particle operators do not commute with the fugacity's operator, definition of the excited states depends on their order, in contrast to the uniform cases.
We have originally speculated that the difference between the particle numbers before and after the projection may cause such asymmetry of the spectra as discussed in the latter part of Sec. III. Let us look again at Eqs. (61,62). Note that the electron removal matrix element is proportional to n iσ (density of the iσ electrons), while the addition is to 1−n i (density of the empty sites). Nevertheless, it is compensated by the fugacity factor and the renormalization factors for the removal and the addition are the same. Then, one may speculate that some asymmetry may appear if we destroy this balance by changing the fugacity factors. However, according to our analysis here, the particle number difference does not have much effect, and it does not cause any asymmetry at least in the uniform systems.

V. SPIN-INDEPENDENT CONSTRAINT
At present, our main interest is in the GAs with spindependent constraints in the former sections because they seem to be more convenient to investigate antiferromagnets, stripe state, impurity systems, and so on. However, in systems with a more complicated spin configuration, the GAs with spin-independent constraints may be useful. Therefore, here we work on it, but only take the leading order with respect to the intersite contractions.

A. Local constraint
A grand canonical GA with a spin-independent constraint, was introduced by Wang et al. 16 In non-magnetic cases (n i↑ = n i↓ ), this is identical to the GA with the spindependent constraint in Secs. II and III. However, in magnetic cases, the results of these two GAs are different. Accordingly, the ferromagnetic homogeneous limit 17 with the spin-independent constraint is not equivalent to that of the canonical GA, but is reduced to the GA for charge-canonical spin-grand-canonical functions explained in the next subsection. Since the formulae for g t and g s have been already derived in Refs. 16,17, here we derive them in a slightly more general form by assuming that S x i and S y i are finite. By replacing λ iσ by λ i in the derivation in Sec. II D and using S ± i ≡ S ± i 0 , we obtain In this case, the derivation of g s is rather simple because S i is nonzero only when it is operated to states where site i is singly occupied. 16 Note that S i S i 0 is automatically satisfied, and there is no complexity appeared in Sec. II D. However, in this formulation, projected quasiparticle excited states are not orthogonal for magnetic systems in general.
B. Global constraint (charge-canonical spin-grand-canonical GA) It is possible that a wave function before the projection is an eigenstate of the total particle number, but not any eigenstate of the total spin. One can formulate a GA also for such systems. In this case, the condition to impose is the canonical condition for the total particle number, iσ n iσ = N after , namely, Here, λ does not depend on σ. The renormalization formulae for S ± i = 0 are the same as those in Sec. IV if λ σ is replaced by λ. The generalization to S ± i = 0 is straightforward.
The results for antiferromagnets are equivalent to those in Sec. IV. The limit to uniform ferromagnets without superconductivity can be taken by setting ∆ ij = 0, N after σ = i n iσ 0 and dropping site indices, 17 , g s σ = n n − 2n ↑ n ↓ 2 .
(122) These are different from our quasi-canonical derivation in Eqs. (104,105). In fact, in this spin-grand-canonical formulation, g s for ferromagnets is the same as that for antiferromagnets. This discrepancy is explained as follows: If the wave function before the projection is an eigenstate of the total spin-z, then the renormalization is represented by Eqs. (104,105). If not, by Eq. (122). The Gutzwiller projection tends to magnify spin moments as explained in Appendix A. However, in the canonical scheme, only xy components of the spins are allowed to be enhanced because of the canonical constraint. On the other hand, the spin-grand-canonical case is free from this constraint, and spins are renormalized more isotropically.

VI. SUMMARY AND DISCUSSION
We have derived various formulae using the grand canonical Gutzwiller approximation with several different constraints imposed by the fugacity factors for inhomogeneous magnetic systems. The formulation with the local particle number conservation yields more simple formulae. On the other hand, the global particle number constraint is more convenient in comparing with the VMC.
In Secs. III, we have discussed the asymmetry of the density of states. Although the incoherent spectra are not taken into account in this paper, we speculate that they appear at much higher energy scale than the coherence peaks. The conventional BCS theory tells us that the quasi-particle excitation spectra are symmetric between positive and negative bias. In contrast, with the Gutzwiller projection, some asymmetry appears. One may think that electron addition is always more difficult than electron removal if repulsion between electrons is strong. However, we doubt if such simple intuition works, and speculate that the asymmetry depends on the Hamiltonian. As a counterexample, we have shown that the projected s-wave superconductor may have the opposite asymmetry. Namely, even with the strong repulsion, the addition spectra can be larger than the removal. We could be able to consider in this way. Let us take two (normalized) Gutzwiller-projected wave functions, |ψ and |φ . Suppose they are the ground states of Hamiltonians, H ψ and H φ , respectively. Furthermore, we assume that |ψ and |φ are also excited states of the other Hamiltonian, H φ and H ψ , respectively. Then, φ|c † iσ |ψ is an electron addition matrix element to the ground state of H ψ . However, if one takes its complex conjugate, it is an electron removal matrix element to the ground state of H φ . Note that | ψ|c iσ |φ | = | φ|c † iσ |ψ |. That is, if an electron addition matrix element for a Hamiltonian is small, an electron removal matrix element for a different Hamiltonian is also small. Therefore, the asymmetry is most likely determined not only by the projection, but also by the Hamiltonian.
Our ultimate purpose is to find good variational wave functions for systems with strong on-site Coulomb repulsion. Once the fugacity factors are introduced, one projected wave function can be related to a number of different unprojected wave functions, each of which is accompanied with fugacity factors of each definition. Therefore, one should probably choose a definition of fugacity factors that matches their purpose. We speculate that the projected optimized solution is similar in any choice of the fugacity factors if the calculation is done accurately enough.
There may be slight disagreement with the results by Ko et al. That is, their comparison between the VMC and the GA seems to say that the GA with the positionand spin-dependent constraint has larger errors than that with the global constraint. However, according to our estimation, the formulation with the position-and spindependent constraint has much smaller errors.
To improve the approximation, one can use techniques of the series expansion method, such as the finite cluster method for calculating higher order terms. One can also use the Padé approximation for extrapolation if necessary, but maybe it is enough just to neglect small terms without extrapolation. Since this linked-cluster expansion can be done analytically, there is a possibility to minimize the energy analytically in contrast to the VMC.
Namely, in the canonical scheme, P G does not change the particle number. On the other hand in the grand canonical scheme, a wave function |Ψ 0 is not an eigenstate ofN , and the particle number is distributed with the distribution function, where P N is an operator which projects onto terms with particle number N . Suppose ρ (0) N is sharply peaked at a mean value N , and fluctuation around it is of O( 1/N ). Then, |Ψ 0 can be regarded as a wave function almost identical with |Ψ N 0 in the thermodynamic limit. In contrast in the projected case, the average particle number of P G |Ψ 0 is in fact different from that of P G |Ψ N 0 . When N is large, an electron has more chance to meet another electron on a certain site. In other words, P G excludes more states with large N , and thus the peak position of the N distribution is shifted to a smaller value. Such distribution change was explicitly estimated by Edegger et al. 12 as summarized below. The distribution function after the projection, can be related to that before by ρ N = g N ρ 0 N with where C is a constant independent of N coming from the normalization of the wave functions. The GA can estimate g N by the ratio of the relative sizes of the projected and unprojected Hilbert spaces as, where N L is the number of lattice sites and N ↑ (N ↓ ) is the number of up (down) spins. We here discuss renormalization of spins using Eq. (A3). Since [P G , S i ] = 0, if a wave function before the projection is an eigenstate of ( i S i ) 2 and/or i S z i with eigenvalues S(S + 1), M , respectively, then these quantities are not changed by P G . If it is not such an eigenstate, P G may change the distribution of S, M . If N is fixed, by changing M in Eq. (A3), one can see that P G excludes more states with small M . As a result, the most "probable" M increases. In fact, Eq. (A3) correctly reproduce the limit of fully polarized states (N σ = 0), which are obviously not affected by P G . By rotating spin axes, and repeating this argument for the x, y directions, we conclude that P G excludes more states with small S 2 . Physically, this can be explained as follows: To make small S 2 , electrons have to cancel their spin moments, and then they have high chance to meet each other on a certain site. When S 2 is large, the spin of an electron tend to orient the same direction as those of the other, then electrons prefer to stay at different sites, and not affected by the projection so much.
The Gutzwiller projection makes more singly occupied sites, and local spin moments also tend to be magnified (this is probably related to increase of S 2 ). In the canonical scheme, magnitude of uniform (ferromagnetic) moments are restricted by the canonical constraint, whereas non-uniform (e.g. antiferromagnetic, sinusoidal) moments are free from the canonical constraint and can be enhanced. On the other hand, in the spin-grandcanonical scheme, the total spin-z component does not have such restriction, and ferromagnetic moments can be also enhanced (shown in Sec. V B).

APPENDIX B: ELECTRON ADDITION/REMOVAL MATRIX ELEMENTS
The general expressions of Eqs. (67-74) are written as follows usingn iσ ≡ 1 − n iσ , α iσ ≡ n iσ (1 − n iσ ) −1 and A ij ≡ n ji↑ n * ji↓ + ∆ * ji ∆ ij (i = j): c † j↑ P c i↑ P 0 P 2 0 ≈ g t0 ii↑ n ji↑ + n ji↓ α j↑Ãij n i↓nj↓ , c † j↓ P c † i↑ P 0 c j↓ P c i↑ P 0 c † j↑ P 2 c i↑ 0 P 2 0 ≈ n ji↑ − n ji↓ α i↑ α j↑Ãij n j↓ni↓ , (B5) APPENDIX C: HIGHER ORDER TERMS FOR THE GLOBAL CONSTRAINT By taking up to the second order of the intersite contractions, P 2 0 is represented by Here, Ξ lm contains all the terms of the second order of intersite contractions in P 2 0 . The division by l Ξ l cancels single site contribution and simplifies the expression. For calculating n i↑ , we need By taking the ratio between Eqs. (C1) and (C2), and neglecting fourth order terms, contribution from disconnected clusters disappears. Then, we obtain renormalization of particle densities as Eq. (93).
To determine λ (2) , we use the equation for 1 − n i↑ + n i↓ , namely, By replacing λ by λ (0) + λ (2) , zeroth order term cancels between the l.h.s. and the first term of the r.h.s. The later also contains λ (2) ; in fact, λ (2) in the other terms can be negligible because they are multiplied to other intersite contractions. Regarding λ (2) as the same order as n 2 ijσ and ∆ 2 ij , and neglecting high order terms, Eq. (95) is obtained.