Superconductivity, antiferromagnetism and phase separation in the two-dimensional Hubbard model: A dual-fermion approach
aa r X i v : . [ c ond - m a t . s t r- e l ] D ec Superconductivity, antiferromagnetism and phase separation inthe two-dimensional Hubbard model: A dual-fermion approach
Junya Otsuki , Hartmut Hafermann , and Alexander I. Lichtenstein Department of Physics, Tohoku University, Sendai 980-8578, Japan Institut de Physique Th´eorique (IPhT), CEA, CNRS, 91191 Gif-sur-Yvette, France Institute of Theoretical Physics, University of Hamburg, 20355 Hamburg, Germany (Dated: September 17, 2018)The dual-fermion approach offers a way to perform diagrammatic expansion around the dynamicalmean-field theory. Using this formalism, the influence of antiferromagnetic fluctuations on the self-energy is taken into account through ladder-type diagrams in the particle-hole channel. The resultingphase diagram for the (quasi-)two-dimensional Hubbard model exhibits antiferromagnetism and d-wave superconductivity. Furthermore, a uniform charge instability, i.e., phase separation, is obtainedin the low doping regime around the Mott insulator. We also examine spin/charge density wavefluctuations including d-wave symmetry. The model exhibits a tendency towards an unconventionalcharge density-wave, but no divergence of the susceptibility is found.
PACS numbers: 71.10.-w, 71.10.Fd, 75.10.-b
I. INTRODUCTION
Magnetism and superconductivity appear nearby intypical phase diagrams of transition-metal and heavy-fermion compounds. Magnetism is related to the Mottinsulating state and heavy-fermion formation, which canbe described in terms of local correlations. On the otherhand, unconventional superconductivity requires spatialcorrelations to be taken into account. For a comprehen-sive understanding, therefore, one needs a unified treat-ment of local correlations and spatial fluctuations, whichhas been a theoretical challenge in the field of stronglycorrelated electron systems.A long-standing problem, which may be related tomagnetism and superconductivity, is the pseudo-gapstate in the low-doped regime of cuprates [1]. One ofthe candidates for its origin is a hidden order, i.e. thestaggered flux state or the d-density wave (d-DW) [2].There are some experiments which indicate broken time-reversal symmetry in the pseudo-gap regime [3, 4]. The-oretically, the mean-field approximation based on theslave-boson representation yields a d-DW in the t - J model [5–7]. However, no clear evidence for the tran-sition has been found in the Hubbard model [8–11].Another feature, which possibly emerges near the Mottinsulator, is a uniform charge instability, i.e., phase sep-aration between two states with different electron den-sity. It was pointed out for the t - J model on the basisof energy arguments [12], and was indeed demonstratednumerically in the one-dimensional system [13] and ininfinite dimensions [14]. In contrast to the d-DW, thephase separation has been observed also in the Hubbardmodel by means of various numerical methods [15–23],while quantum Monte Carlo investigations reported noevidence of phase separation [24, 25].The unresolved problems described above motivate usto investigate the two-dimensional Hubbard model as aprototypical model of strongly correlated electron sys-tems, and to develop new theories which could clarify these issues. The dynamical mean-field theory (DMFT)provides a description of the Mott transition [26] and itscluster extensions provide a route to the d-wave super-conductivity (d-SC) in the doped regime [27, 28]. Thed-SC has indeed been obtained in several numerical cal-culations [18, 29–31]. We note that cluster DMFT par-ticularly accounts for short-range correlations in additionto the local ones.A different kind of extension of single-site DMFT hasbeen worked on, which, in contrast to cluster exten-sions, aims at incorporating long-range correlations [32–38]. The common idea of these approaches is to intro-duce an additional step of solving the lattice problemin a certain way after the DMFT equations are solved.The various formulations differ (i) physically, in the setsof diagrams which are summed beyond DMFT and (ii)technically, how double counting of correlation effects isavoided, that may arise when two different methods arecombined.Rubtsov et al. introduced an auxiliary fermion whichmediates itinerancy of electrons [39, 40]. With this dualfermion, a perturbation expansion around the DMFT hasbeen made possible without the double-counting prob-lem; the zeroth-order approximation in this theory cor-responds to DMFT, and spatial correlations are system-atically incorporated by summing up a series of diagrams.In particular, ladder diagrams similar to those in the fluc-tuation exchange approximation (FLEX) [41–44] yielddescriptions of collective modes (long-range fluctuations).Indeed, it has been shown that inclusion of the ladderdiagrams in the dual-fermion approach leads to param-agnon excitations that exhibit antiferromagnetic (AFM)fluctuations in the paramagnetic state [45, 46]. At thesame time, the ladder approximation yields suppressionof the AFM phase transition in two dimensions [45, 46]and the expected critical exponents in case that phasetransitions are found [47], demonstrating that long-rangefluctuations essential for the critical behavior are appro-priately included. Therefore, the dual-fermion approach T U AFM ( λ sp = 0.90)(0.94)(0.98)AFM (DMFT)4 t / U Figure 1: (Color online) A phase diagram at half filling, δ = 0. with ladder-type diagrams provides a combined descrip-tion of strong local correlations and long-range correla-tions.Although first results of the ladder approximation havebeen presented in 2009 [45, 48], its exemplary results fordoped Mott insulators have been limited because of sometechnical difficulties arising from strong AFM fluctua-tions. In this paper, we overcome these limitations andpresent systematic results for the doped regime of thetwo-dimensional Hubbard model. We address possiblephase transitions of the d-DW and the phase separationin the doped Mott insulator as well as the d-SC. Our re-sults reveal further characteristics of the ladder approxi-mation.The rest of this paper is organized as follows. Inthe next section, we first present phase diagrams ob-tained in this investigation to give an overview of ourresults. Afterwards, the dual-fermion formalism and theself-energy equation are presented in Section III. Suc-ceeding Sections IV–VII present detailed numerical re-sults and related formulas for the AFM susceptibility,superconductivity, phase separation, and unconventionaldensity waves. The paper is closed with discussions inSection VIII. II. OVERVIEW
Prior to presenting formalism and detailed numericalresults, we first give an overview of our results obtained inthis paper. We investigate the two-dimensional Hubbardmodel: H = X k σ ǫ k c † k σ c k σ + U X r n r ↑ n r ↓ , (1)with ǫ k = − t (cos k x + cos k y ). The number operator n r σ is defined by n r σ = N − P kq c † k σ c k + q σ e i q · r , where T δ U = 8 AFM ( λ sp = 0.90)(0.94)(0.98)AFM (DMFT)d-SCPS Figure 2: (Color online) Phase diagrams under doping δ =1 − n for U = 8. N denotes the number of lattice sites. We take t = 1 asthe unit of energy.In two-dimensional systems, the AFM transition is for-bidden at T > χ ∼ e cβ of the suscep-tibility at low temperatures [50, 51]. Our approximationindeed shows no AFM transition within calculated tem-peratures. To quantify the AFM fluctuations, we definea “phase boundary” by the points where the fluctuationsexceed a certain criterion (see Section IV for details). Wemay regard this line as a phase boundary in quasi-two di-mensions. The phase diagram at half filling obtained inthis way is shown in Fig. 1. We plot three phase bound-aries corresponding to different criteria. In DMFT, thereexists a real phase transition, which is plotted for com-parison.According to a cluster DMFT calculation with a para-magnetic bath [52], the Mott transition takes place at U ≃ T ≃ . δ = 1 − n for U = 8. The d-SC is ob-tained in the region T . .
05 and δ . .
18. The su-perconducting transition temperature T c monotonicallyincreases approaching half filling ( δ = 0). This behav-ior is reminiscent of the FLEX [43, 44] and differs fromthat in cluster DMFT, where the d-SC phase exhibits amaximum at finite doping [29, 31]. We consider that themonotonic behavior of T c in our results is due to insuf-ficient treatment of short-range spin fluctuations, whichwill be discussed in Sec. VIII.In the low-doping regime above T c , we found a phaseseparation. The line T PS in Fig. 2 shows the spinodal line,where the uniform charge susceptibility diverges. Thephase separation extends up to δ ≃ .
15. At
T < T PS ,the solution is thermodynamically unstable because of ∂n/∂µ <
0. Thermodynamic stability is acquired by in-homogeneous coexistence of regions with different dopinglevels: Mott insulating regions with δ = 0 and metallicregions with larger doping δ = 0. The phase boundaryfor the d-SC has been computed with the homogeneoussolution, which, in fact, is thermodynamically unstablein the region δ . .
15. Therefore, pure d-SC only realizesin the region 0 . . δ . .
18, while it may not occur for δ . .
15. We have also examined the possibility of a d-DW. We find that the d-DW dominates over DWs withother symmetries, but the corresponding susceptibilityshows no divergence.
III. DUAL-LADDER APPROXIMATIONA. Dual action
In the dual-fermion approach, the lattice model issolved in two steps. First, an effective impurity model,which is the same as in DMFT, is solved with the aid ofsome numerical methods. The local correlations, whichare essential for formation of the Mott gap, are fully takeninto account at this stage. In the next step, an interactinglattice model is constructed by quantities evaluated in thefirst step, and is solved by a diagrammatic perturbationtheory. This way, spatial fluctuations are incorporatedin addition to the local correlations in DMFT. In the fol-lowing, we first give a brief summary of the dual-fermionapproach [39, 40, 46].It is convenient to work in the path-integral rep-resentation. The partition function Z is written interms of Grassmann variables c r ( τ ) and c ∗ r ( τ ): Z = R Q r D [ c ∗ r c r ] e −S [ c ∗ ,c ] . The action S is given by S [ c ∗ , c ] = X r S imp [ c ∗ r , c r ] + X ω k σ ( ǫ k − ∆ ω ) c ∗ ω k σ c ω k σ , (2)where the Fourier transform of c ( τ ) is defined by c ω = β − / R β c ( τ ) e iωτ with ω being the fermionic Matsubarafrequency. The first term S imp describes the effectiveimpurity model of DMFT [26]: S imp [ c ∗ r , c r ] = − X ωσ ( iω + µ − ∆ ω ) c ∗ ω r σ c ω r σ + U Z dτ n r ↑ ( τ ) n r ↓ ( τ ) . (3)The hybridization function ∆ ω is actually canceled outin Eq. (2), but an approximate solution may depend on∆ ω . A condition for determining ∆ ω will be discussedlater.In order to construct a lattice model for which thesolution of S imp is the starting point, Rubtsov et al. introduced an auxiliary fermion which “decouples” thekinetic-energy term [39, 40] [the second term in Eq. (2)].This fermion is termed dual fermion and represented by f . The dual fermions locally hybridize with the electronsand mediate the electron itinerancy. The point is thatthe transformed action written with c and f variables hasonly local terms concerning c variables. Therefore, onecan integrate out c variables at each site independently.This process corresponds to solving the effective impu-rity problem expressed by S imp . The local hybridizationbetween c and f introduces effective interaction termsamong the f variables, which are local in space but non-local in the time domain.The resulting partition function thus consists only ofthe dual variables f . Hence, our task now is to solvethe dual system described by ˜ Z = R Q r D [ f ∗ r f r ] e − ˜ S [ f ∗ ,f ] .The action ˜ S is given by˜ S [ f ∗ , f ] = − X ω k σ ( ˜ G ω k ) − f ∗ ω k σ f ω k σ + ˜ V [ f ∗ , f ] , (4)with the bare dual Green’s function ˜ G ω k defined by˜ G ω k = ( g − ω + ∆ ω − ǫ k ) − − g ω . (5)Here g ω = −h c ω r σ c ∗ ω r σ i imp is the impurity Green’s func-tion, with h· · · i imp being a thermal average with respectto the action S imp . The first term in Eq. (5) correspondsto the lattice Green’s function in DMFT. Subtractingthe second term excludes double counting of local cor-relations. The term ˜ V denotes local interactions, whichinclude many-body interactions as well as a two-bodyterm. The point of the transformed action ˜ S is that thebare propagator ˜ G and the interaction ˜ V fully includelocal correlations. Hence, the (undressed) dual fermions f may be regarded as particles which involve all the localinteraction processes. Residual interactions between thedressed particles are described by ˜ V .Once the dual Green’s function ˜ G ω k = −h f ω k σ f ∗ ω k σ i ˜ S is evaluated, it is readily transformed to the Green’s func-tion G ω k of the original electrons by means of the exactrelation G − ω k = ( g ω + g ω ˜Σ ω k g ω ) − + ∆ ω − ǫ k . (6)Here, we introduced the dual self-energy ˜Σ ω k =( ˜ G ω k ) − − ˜ G − ω k . It is clear from this expression that˜Σ ω k = 0 leads to the DMFT formula for the latticeGreen’s function.The formalism presented above is still exact. In thefollowing, two approximations will be made. Firstly, weretain only two-body interactions in ˜ V . Secondly, we per-form a perturbation expansion with respect to ˜ V to sumup a certain set of diagrams for ˜Σ ω k . These approxima-tions rely on the idea that the DMFT is a good startingpoint for Mott insulators, and hence the spatial corre-lations may be dealt with perturbatively. We can alsoendorse this treatment by arguments based on the 1 /d expansion, which will be discussed later. B. Interaction vertex for dual fermions
We retain only two-body interactions in ˜ V neglectingterms involving more than three particles. Thus, ˜ V reads˜ V = − X kk ′ q X σ σ σ σ γ σ σ σ σ ωω ′ ; ν f ∗ kσ f ∗ k ′ + q,σ f k ′ σ f k + q,σ , (7)where k = ( ω, k ) and q = ( ν, q ) with ν being the bosonicMatsubara frequency. The interaction coefficient γ cor-responds to the vertex evaluated in the effective impuritysystem. It is defined through h c c c ∗ c ∗ i imp = g g ( δ δ − δ δ ) + T g g γ g g . (8)Here, a simplified notation is used such as 1 ≡ ( ω , σ ).Using energy conservation, we parameterize the fre-quency dependence of γ as γ σ σ σ σ ωω ′ ; ν ≡ γ ( ω,σ ) , ( ω ′ + ν,σ ) , ( ω ′ ,σ ) , ( ω + ν,σ ) . (9)The antisymmetric nature of γ leads to the relation γ σ σ σ σ ωω ′ ; ν = − γ σ σ σ σ ω,ω + ν ; ω ′ − ω . The interaction ˜ V is rep-resented by the diagram in Fig. 3(a).We consider the spin dependence of γ . Using the Paulimatrix σ ξ ( ξ = 0 , x, y, z , including the unit matrix σ ),we transform γ as γ σ σ σ σ = 12 X ξξ ′ γ ξξ ′ σ ξσ σ σ ξ ′ σ σ . (10)Without magnetic field, γ ξξ ′ is diagonal andthere are only two independent components: γ = diag( γ ch , γ sp , γ sp , γ sp ). By inverting Eq. (10),we obtain γ ch ≡ γ = 12 X σσ ′ γ σσ ′ σ ′ σ , (11) γ sp ≡ γ zz = 12 X σσ ′ σσ ′ γ σσ ′ σ ′ σ , (12)which corresponds to interactions in charge andlongitudinal-spin channels, respectively. The transverse-spin channel γ ⊥ ≡ γ ↑↓↑↓ = ( γ xx − iγ xy ) is equivalent to γ sp , since we are considering the paramagnetic state. C. Self-energy
We evaluate the dual self-energy ˜Σ ω k taking the two-body interaction ˜ V in Eq. (7) into account. In principle,one can apply any numerical method as well as approx-imations for this purpose. Here, we use a perturbationtheory and sum up certain diagrams based on physicalconsiderations and a 1 /d analysis. = + + (a) (b)(c) Figure 3: Diagrammatic representations of (a) bare interac-tion ˜ V for dual fermions, (b) the dual self-energy ˜Σ in theladder approximation, and (c) the equation for the renormal-ized vertex Γ. We first discuss the choice of diagrams from a physicalpoint of view. Near the Mott insulator, the importantingredient are spin fluctuations, which can be taken intoaccount by ladder-type diagrams. Indeed, the ladder ap-proximation gives a magnetic spectrum which exhibitslow-energy spin excitations (so-called paramagnons), asa consequence of strong AFM fluctuations [45, 46]. Inthe following, we present the self-energy formula in theladder approximation, stressing on the SU(2) symmetryfor the spin indices.We first evaluate the renormalized vertex Γ collectingsuccessive particle-hole excitations. Since the propagator˜ G is independent of the spin component, the vertex Γ α foreach channel α = ch , sp independently obeys the Bethe-Salpeter equation:Γ αωω ′ ; ν q = γ αωω ′ ; ν + T X ω ′′ γ αωω ′′ ; ν ˜ χ ω ′′ ; ν q Γ αω ′′ ω ′ ; ν q , (13)where ˜ χ is defined by˜ χ ω ; ν q = − N X k ˜ G ω k ˜ G ω + ν, k + q . (14)Figure 3(c) shows the diagrammatic representation of theabove equation. The self-energy is evaluated with therenormalized vertex Γ as follows:˜Σ ω k = − TN X ω ′ k ′ γ ch ωω ′ ;0 ˜ G ω ′ k ′ + T N X ν q ˜ G ω + ν, k + q ( V ch + 3 V sp ) ωω ; ν q , (15)where the effective interaction V α is defined by V αωω ′ ; ν q = T X ω ′′ γ αωω ′′ ; ν ˜ χ ω ′′ ; ν q (cid:2) αω ′′ ω ′ ; ν q − γ αω ′′ ω ′ ; ν (cid:3) . (16)The corresponding diagram for the dual self-energy isshown in Fig. 3(b). (a) (b) (c) Figure 4: Self-energy diagrams: (a) Example of a leadingorder diagram in terms of 1 /d but not included in the presentapproximation for physical reasons (see text), (b) an exampleof the next-leading contributions in the 1 /d expansion, and(c) zero contribution by the self-consistency condition. D. A perspective from a /d expansion The momentum dependence of the self-energy disap-pears in the limit of d = ∞ dimensions. This means thatDMFT provides the exact solution of fermionic modelsin d = ∞ [26, 55]. Since the dual-fermion approach offersan expansion around DMFT, it is reasonable to classifydiagrams in terms of 1 /d . In this view, we reconsider theself-energy diagram presented above.We consider the large- d limit with the scaling t ∝ / √ d [26, 55]. The local Green’s function G r =0 of theoriginal electrons is of zeroth order in 1 /d , while its dualcounterpart vanishes, ˜ G r =0 = 0, by the self-consistencycondition given later (we omit the ω index for simplicity).Hence, the dual Green’s function has only intersite com-ponents ˜ G r =0 which scale as ˜ G r =0 ∼ G r =0 ∼ O (1 / √ d ).The second-order diagram for ˜Σ ω k (the second term inFig. 3(b) with the renormalized vertex replaced by thebare interaction) has a contribution of order O (1 / √ d ).The ladder diagrams summed up in Fig. 3(c) are of thesame order, because the factor 1 /d arising from the twopropagators is canceled by the lattice summation. There-fore, all diagrams included in the second-term of Fig. 3(b)provide the leading contributions of order 1 / √ d . Ac-tually, ladder diagrams in the particle-particle channelalso have contributions of the same order (e.g., the di-agram in Fig. 4(a)). However, we may neglect themsince particle-particle fluctuations only have a minor ef-fect in the doped Mott insulator. Diagrams containinghigher-order vertices only appear at second-to-leading or-der [e.g., Fig. 4(b)], since the local diagram like Fig. 4(c)vanishes as explained later. It means that the three-bodyand higher-order interactions do not enter to leading or-der of the 1 /d expansion. Indeed, it has been numericallyconfirmed that the ladder-type diagrams dominate overdiagrams built from the three-particle vertex [45]. Inconclusion, the dual-ladder self-energy in Eq. (15) con-stitutes the leading correction to the DMFT around the d = ∞ limit. E. Self-consistency condition
So far, the hybridization function ∆ ω is arbitrary.We discuss here how to determine ∆ ω . The conditionis that the scheme should reduce to DMFT if no self-energy corrections are taken into account: ˜Σ ω k = 0. Thefollowing self-consistency condition fulfills this require-ment [39, 40, 46]: X k ˜ G ω k = 0 . (17)It is clear from Eq. (5) that when ˜Σ ω k = 0, this condi-tion leads to the DMFT self-consistency condition. Fur-thermore, it eliminates the contribution from the Hartreediagram (the first term in Fig. 3(b)). Similarly, all dia-grams which have a propagator connecting the same site(local loops) give no contribution [e.g., Fig. 4(c)]. F. Technical details
We solve the effective impurity problem using thehybridization-expansion solver (CT-HYB) [56] of thecontinuous-time quantum Monte Carlo method [57, 58].The vertex γ ωω ′ ; ν as well as g ω are computed. We ap-plied an efficient implementation for the vertex calcula-tion [59]. The vertex γ ωω ′ ; ν is computed in a small en-ergy window, while the energy cutoff for g ω can be takensufficiently large (10 –10 Matsubara frequency points)in the CT-HYB algorithm. To be concrete, we restrictthe frequencies of γ ωω ′ ; ν to | ω | , | ω ′ | ≤ (2 n c + 1) πT and | ν | ≤ m c πT . Typically, we take n c = m c = 20 for T & .
1, and up to n c = m c = 60 for lower temperatures.Such a small cutoff compared to the one for g ω is possi-ble because the frequency summation of γ ωω ′ ; ν is alwaystaken with ˜ G ω k , which decays faster than the ordinaryGreen’s function, ˜ G ω k ∼ − ǫ k /ω . We note that the neg-ative bosonic frequencies, ν <
0, need not be computedsince we have the relation γ αωω ′ ; ν = ( γ α − ω − ω ′ ; − ν ) ∗ .The quantities g ω and γ ωω ′ ; ν are plugged into thedual-lattice calculations. The momentum summation(convolution) in Eq. (15) is evaluated in the real space.Here we can use FFT to reduce O ( N ) calculation into O ( N log N ). On the other hand, we simply add the fre-quency ν in Eq. (15), since the frequency summation doesnot simplify considerably in the imaginary-time domaindue to the full frequency dependence of γ ωω ′ ; ν . The lat-tice size N is fixed at N = 32 ×
32 (excepting Fig. 6).This size is sufficiently large for our purpose of reveal-ing possible phase transitions. A larger system size isnecessary to observe the critical behavior, which will bediscussed in the next section. We compute ˜Σ ω k and ˜ G ω k iteratively until they are converged. To get convergence,we mix new and old data of ˜Σ ω k . The weight of thenew data ranges from 0.5 down to 0.02 by checking thetendency toward convergence.After ˜Σ ω k is obtained, we update the bath ∆ ω and goback to the impurity problem. We use the formula [40]∆ new ω = ∆ old ω + ξ ˜ G ω, r =0 / [ g ω ( g ω + ˜ G ω, r =0 )]. Here, ξ is themixing parameter and typically we take ξ = 0 . ω is unstable. In thiscase, we need an elaborate treatment of ˜Σ ω k to avoid thenumerical instability (see Appendix A for details). IV. ANTIFERROMAGNETIC SUSCEPTIBILITY
We first present numerical results for the AFM sus-ceptibilities at half filling, n = 1. We shall show datafor lower temperatures than in the previous calculationof the ladder approximation [45], and discuss the pointof our calculations to achieve convergence in the criticalregime. The results may also be regarded as a benchmarkof our calculations.The spin and charge susceptibilities of the original elec-trons, χ αν q , are connected to the reducible vertex part ofthe dual fermions by an exact relation [46, 60, 61]. Inthe ladder approximation, we use Γ α in Eq. (13) for thereducible vertex to obtain explicit expression for χ αν q as χ αν q = χ ν q + T X ωω ′ X ω ; ν q Γ αωω ′ ; ν q X ω ′ ; ν q , (18)where χ ν q = − TN X ω k G ω k G ω + ν, k + q , (19) X ω ; ν q = − N X k G ω k G ω + ν, k + q R ω k R ω + ν, k + q , (20)with R ω k = g − ω (∆ ω − ǫ k ) − . Equation (18) reduces tothe DMFT formula if the DMFT limit, ˜Σ ω k = 0, is taken.Figure 5(a) shows the temperature dependence of theinverse of the static AFM susceptibility χ sp ν =0 , Q , where Q = ( π, π ) is the nesting vector. The DMFT result isalso plotted for comparison. From these data, we findthat the susceptibility does not diverge in the ladder ap-proximation, while the DMFT susceptibility diverges andobeys the Curie-Weiss law.In the regime where fluctuations are strong, or moreprecisely, at T . T DMFTN with T DMFTN being the DMFTN´eel temperature, the method presented in Appendix Ais essential to achieve convergence. Here we only mentionthat this method does not change the equations, and issimply a way of obtaining a converged solution. A spu-rious divergence of χ sp , which may arise during the iter-ation, is removed. In this procedure, the main quantitywe need to check is the dimensionless matrix ˆ A definedby ( ˆ A ) ωω ′ = γ sp ωω ′ ;0 ˜ χ ω ′ ;0 , Q . The condition for the diver-gence of Γ sp in Eq. (13), and hence of χ sp is λ sp = 1where λ sp denotes the largest eigenvalue of ˆ A [62]. Thetemperature dependence of λ sp is shown in Fig. 5(b). Itturns out that λ sp approaches 1 with decreasing T in theladder approximation. T (a) 1 / χ sp ν =0, Q U = 48 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0.1 1 T (b) λ sp U = 48 Figure 5: (Color online) Temperature dependence of (a) thestatic AFM susceptibility χ sp ν =0 , Q and (b) the largest eigen-value λ sp of the matrix ˆ A at half-filling n = 1. The closedsymbols and open symbols show results in the present ap-proximation and within DMFT, respectively. -3 -2 -1
0 1 2 3 4 5 6 7 8 9 – λ s p T ∼ exp( – ∆ / T ) ∆ = 1.9 ∆ = 1.7 U = 4 ( N = 32 × N = 64 × N = 128 × U = 8 ( N = 32 × N = 64 × N = 128 × Figure 6: (Color online) A scaling plot: 1 − λ sp as a function of1 /T . Results for different system sizes, N = 32 ×
32, 64 × × − λ sp ∝ exp( − ∆ /T ). In the critical regime, the susceptibility diverges ex-ponentially toward T = 0: χ ∼ e β ∆ [50, 51]. It followsthat λ sp approaches 1 according to 1 − λ sp ∝ e − β ∆ . Inorder to check this behavior, we plot 1 − λ sp as a func-tion of 1 /T in Fig. 6. Results for larger system sizes, N = 64 ×
64 and 128 × − λ sp . − .It indicates that the slow decays for N = 32 ×
32 and64 ×
64 observed at 1 /T & − λ sp ∝ e − β ∆ indicatedby the solid lines. We thus conclude that our approxima-tion correctly reproduces the N´eel temperature of T N = 0required from the Mermin-Wagner theorem. V. SUPERCONDUCTIVITYA. Formulas for pairing susceptibilities
In this section, we discuss the superconductivity in thedoped regime. We first derive a formula for the pairingsusceptibility of the dual fermions. The susceptibilitiesof the dual fermions can be transformed to those of theoriginal electrons [46, 60, 61]. Actually, numerical trans-formations cannot be performed in the case of unconven-tional (momentum-dependent) order parameters becausethe susceptibility matrix is too large to store in memory[see Eq. (22)]. However, since the diverging point is com-mon to both susceptibilities, we can determine the tran-sition temperature from the dual-fermion susceptibilitywithout transforming to the electron susceptibility.We consider Cooper pairs with opposite spin directionsof the constituent electrons. With a form factor φ k whichdepends on both k and ω , the order parameter Φ is ex-pressed as Φ = P k φ k h f k ↑ f − k ↓ i ˜ S . The static susceptibil-ity for this pairing is defined by P kk ′ φ k ˜ P kk ′ φ ∗ k ′ where˜ P kk ′ = h f k ↑ f − k ↓ f ∗− k ′ ↓ f ∗ k ′ ↑ i ˜ S . (21)The Bethe-Salpeter equation for this Green’s function iswritten as˜ P kk ′ = ˜ P k δ kk ′ − TN X k ′′ ˜ P k Γ pp kk ′′ ˜ P k ′′ k ′ , (22)where ˜ P k = ˜ G k ˜ G − k . (23)For the irreducible vertex part Γ pp , we take account ofeffective interactions mediated by the spin and chargefluctuations. Hence, Γ pp is given in terms of the renor-malized vertex in Eq. (13) as [46, 48]Γ pp kk ′ = − Γ ↑↓↓↑ ω, − ω ′ ; ω ′ − ω, k ′ − k + Γ ↑↓↑↓ ω,ω ′ ; − ω − ω ′ , − k − k ′ + γ ↑↓↓↑ ω, − ω ′ ; ω ′ − ω . (24)The first term in Eq. (24) incorporates the charge andlongitudinal spin fluctuations, and the second term thetransverse spin fluctuations. The third term subtractstheir double counting. A diagrammatic representationfor Γ pp is shown in Fig. 7.Without magnetic field, the pairing susceptibility isclassified according to the total spin of the pair. For thispurpose, we replace the pair operator by its symmetrizedor anti-symmetrized form: f k ↑ f − k ↓ → √ f k ↑ f − k ↓ ∓ f k ↓ f − k ↑ ) . (25)Here, − corresponds to the spin singlet and + to thespin triplet. The corresponding pairing susceptibility isexpressed as ˜ P ± k,k ′ = ˜ P k,k ′ ± ˜ P k, − k ′ . (26) + − = Figure 7: The pairing interaction (the irreducible vertex forthe pairing susceptibility) Γ pp in the ladder approximation.The box with stripes stands for the renormalized vertex Γ inFig 3(c). Hence, the inversion of the fermionic frequency and mo-mentum, k = ( ω, k ) → − k = ( − ω, − k ), transforms ˜ P ± kk ′ as P ± k,k ′ = ± P ± k, − k ′ = ± P ±− k,k ′ = P ±− k, − k ′ . From Eq. (22),we obtain the equation for ˜ P ± kk ′ ,˜ P ± kk ′ = ˜ P k ( δ k,k ′ ± δ k, − k ′ ) − TN X k ′′ ˜ P k Γ pp ± kk ′′ ˜ P ± k ′′ k ′ , (27)where the (anti-)symmetrized vertex Γ pp ± kk ′ is defined byΓ pp ± kk ′ = (Γ pp k,k ′ ± Γ pp k, − k ′ ) /
2. Their explicit expressionsreadΓ pp+ kk ′ = 14 h (3Γ sp − Γ ch ) ω, − ω ′ ; ω ′ − ω, k ′ − k − γ sp ω, − ω ′ ; ω ′ − ω i + ( ω ′ → − ω ′ ) , (28)Γ pp − kk ′ = 14 h − (Γ sp + Γ ch ) ω, − ω ′ ; ω ′ − ω, k ′ − k + 2 γ sp ω, − ω ′ ; ω ′ − ω i − ( ω ′ → − ω ′ ) , (29)where ( ω ′ → − ω ′ ) is symbolic for the terms appearingbefore it with ω ′ replaced by − ω ′ .The dimension of the matrices is too large to solveEq. (27) numerically. We instead deal with an eigen-value problem to determine the transition temperatureand to extract the dominant pairing fluctuations. Nearthe transition temperature, we may neglect the first termin Eq. (27) to obtain the linear equationˆ K ± φ = λ SC φ, ( ˆ K ± ) kk ′ = − TN ˜ P k Γ pp ± kk ′ . (30)We can demonstrate from the explicit form of Γ pp ± kk ′ thatthe eigenvalues λ SC are purely real. The condition for thedivergence of the susceptibility is λ SCmax = 1 with λ SCmax being the largest eigenvalue. The corresponding eigen-function φ k gives the form factor of the order parameter. B. Numerical results
We evaluated the largest eigenvalues λ SCmax of Eq. (30)by a kind of power method. In this calculation, we en-forced a particular spatial symmetry to pick up an eigen-function belonging to a certain irreducible representation(see Appendix B for details). In this way, we computed spin singletA spin triplet -1-0.5 0 0.5 1 k y / π –0+-1-0.5 0 0.5 1 k y / π –0+-1-0.5 0 0.5 1 k y / π –0+-1 -0.5 0 0.5 1 k x / π -1-0.5 0 0.5 1 k y / π –0+-1-0.5 0 0.5 1 k y / π –0+-1 -0.5 0 0.5 1 k x / π -1-0.5 0 0.5 1 k y / π –0+ -1-0.5 0 0.5 1 k y / π –0+-1-0.5 0 0.5 1 k y / π –0+-1-0.5 0 0.5 1 k y / π –0+-1 -0.5 0 0.5 1 k x / π -1-0.5 0 0.5 1 k y / π –0+-1-0.5 0 0.5 1 k y / π –0+-1 -0.5 0 0.5 1 k x / π -1-0.5 0 0.5 1 k y / π –0+ A B B E u ( x )E u ( y ) Figure 8: (Color online) Momentum dependence of the eigen-functions φ ω , k of Eq. (30) for U = 8, δ = 0 .
14 and T = 0 . ω = πT . Either the even-frequency part φ even ω , k or theodd-frequency part φ odd ω , k is plotted depending on which isallowed by the Pauli principle.
10 types of pairings (2 spin symmetries × φ k be-come real. Finally, we define even- and odd-frequencyparts, φ even ω k = φ ω k + φ − ω k and φ odd ω k = φ ω k − φ − ω k , tosee the frequency dependence. We have confirmed thateither φ even or φ odd vanishes to fulfill the Pauli principle,e.g., φ odd = 0 for the spin-singlet with symmetry A .We first show eigenfunctions φ ω k obtained in the waydescribed above. Figure 8 shows the momentum de-pendence of φ ω k with the lowest Matsubara frequency, ω = πT . The main feature is that some functions haveonly minimal nodes required from the symmetry and therest have additional nodes. In the A symmetry, for ex- λ S C T singlet A1gA2gB1gB2gEutriplet A1gA2gB1gB2gEu Figure 9: (Color online) Temperature dependence of theeigenvalues λ SC of Eq. (30) for U = 8 and δ = 0 . ample, there is no node for the triplet, while a line nodeexists on the Fermi level for the singlet (i.e., extendeds-wave symmetry, cos k x + cos k y ).Which type of superconductivity actually occurs is ex-amined from the temperature dependence of λ SC . It canbe seen from Fig. 9 that λ SC for the spin-singlet B (d x − y ) symmetry crosses 1 as expected. The transitiontemperature T c is estimated to be T c ≃ .
030 for theseparameters. The doping dependence of T c is plotted inthe phase diagram in Fig. 2. VI. PHASE SEPARATION
Our next interest lies in the paramagnetic state above T c and near the Mott insulator. In this regime, we foundan instability of the uniform charge fluctuations. Fig-ure 10 shows the temperature dependence of the chem-ical potential µ for several values of doping δ = 1 − n for U = 8. The decrease of µ below T ≃ T = 0 .
1, somelines for different doping levels intersect. It means that µ is a non-monotonic function of δ at low temperaturesas shown in the inset of Fig. 10. This behavior indicatesa phase separation as explained below.At T = 0 . δ and δ . Actually,the Mott insulator with δ = 0 is also a solution in thiscase. Hence, there are three solutions ( δ = 0 < δ < δ ),two of which ( δ and δ ) are thermodynamically stableand one ( δ ) is unstable. In order to make the averagedoping ¯ δ at 0 < ¯ δ < δ , the system becomes spatiallyinhomogeneous between the Mott insulator with δ = 0and the metallic state with δ = δ .We define the temperature T PS for the phase separa-tion by the point where two lines intersect in Fig. 10.It corresponds to the so-called spinodal point where the µ T δ = 00.10.2 0.8 1.2 1.6 2 2.4 0 0.1 0.2 0.3 δ T = 0.2510.1580.100 Figure 10: (Color online) Temperature dependence of thechemical potential µ . The doping δ is varied from 0 to 0.2in 0.02 steps. The inset shows µ as a function of δ for fixed T . uniform charge susceptibility diverges. The result for T PS is plotted in the phase diagram of Fig. 2. Below T PS , thehomogeneous solution is thermodynamically unstable.Before concluding the paragraph, we comment on atechnical issue related to these observations. The uni-form charge susceptibility can be computed in two differ-ent ways: Either from the chemical potential as discussedabove or from the correlation function as presented inSection IV for the spin channel. In the present approx-imation, the two results are not consistent. Indeed, wefound no divergence of the uniform charge susceptibilitycomputed from the correlation function. In this case, theone computed from the chemical potential is more reli-able in the sense that derivative of the self-energy with re-spect to µ is taken strictly, while the correlation functionincorporates only a part of the corresponding diagrams.To improve consistency, i.e., to obtain divergence in thecorrelation function, we need more elaborated treatmentof the irreducible vertex to satisfy the Ward identity [63]. VII. UNCONVENTIONAL DENSITY WAVES
In the previous section, we discussed the phase sepa-ration taking place near the Mott insulator. In this sec-tion, we examine the possibility of another phase, whichhas been discussed extensively, namely, the staggered fluxstate or the d-DW state [2, 5, 8–10]. To make our formu-lation general, we consider both spin and charge channels( α = sp , ch), arbitrary wave vectors q , and arbitrary spa-tial symmetry. The order parameter Ψ αη q is defined byΨ αη q = X ω k σ σ ασσ ψ η k h f ∗ ω k σ f ω k + q σ i ˜ S , (31)where σ ch = σ and σ sp = σ z . The index η labelsdifferent form factors ψ η k . The d-DW corresponds to Figure 11: The local-current configuration in the d-DW state. Ψ ch , d Q with the form factor ψ d k = i (cos k x − cos k y ), whilethe ordinary DW is given by ψ s k = 1. In the real-space representation, the d-DW exhibits a local current i h f ∗ r f r + x i − i h f ∗ r + x f r i which aligns as in Fig. 11. Fol-lowing the same reasoning as for the pairing correla-tions, we consider susceptibilities of the dual fermions.The susceptibility corresponding to Ψ αη q is given by P kk ′ ψ ηk ˜ χ αkk ′ ; q ( ψ ηk ′ ) ∗ with˜ χ αkk ′ ; q = 12 X σσ ′ σ ασσ σ ασ ′ σ ′ h f ∗ kσ f k + qσ f ∗ k ′ + qσ ′ f k ′ σ ′ i ˜ S . (32)The susceptibility formula in Eq. (18) does not give riseto the unconventional DW, since the irreducible vertex islocal in this formula. Higher-order processes need to betaken into account.In order to see which processes enhance unconventionalDWs, it is instructive to analyze influences of intersiteinteractions H int in the mean-field approximation or therandom-phase approximation (RPA) [64, 65]. We con-sider the nearest-neighbor repulsion V and the AFM ex-change interaction J : H int = 12 X h i,j i X σσ ′ h V c † iσ c iσ c † jσ ′ c jσ ′ + Jc † iσ c iσ ′ c † jσ ′ c jσ i . (33)The staggered susceptibility χ αη Q in RPA is given by χ αη Q = [( χ η Q ) − − I αη Q ] − , (34)where χ η Q = − ( T /N ) P k | ψ η k | G ω k G ω k + Q . The sign of I αη Q determines whether the fluctuations are enhanced( I αη Q >
0) or suppressed ( I αη Q < I αη Q from U , V and J are summarized in Table I [64].There are two types of diagrams in the RPA: The bubble-type (upper table row) and the ladder-type (lower). Itturns out that the ladder-type diagram for V and J maycause an unconventional CDW and/or SDW [66]. Forthe Hubbard model without V and J , the local repulsion U gives rise to AFM fluctuations, which effectively havean effect similar to that of the J -term. Hence, there isa chance that the unconventional CDW is induced byprocesses beyond RPA. This spin-fluctuation mediatedinteraction is taken into account, for example, by thediagram in Fig. 12(a). This process is indeed included inthe irreducible vertex derived by a functional derivativeof the Luttinger-Ward functional in the FLEX [41].0 Table I: Effect of different diagrams and interaction terms onthe staggered susceptibilities in the RPA. The signs + and − indicate enhancement and suppression, respectively. U V J +SDW − CDW +CDW +SDW+CDW+SDW +uSDW+uCDW +uCDW (a) (b)
Figure 12: (a) An exemplary susceptibility diagram whichcontributes to the unconventional CDW. (b) A diagram forthe dual-fermion irreducible vertex Γ ′ in Eqs. (36) and (37).The box with stripes stands for the renormalized vertex Γ inFig 3(c). In the dual-fermion approach, the effective interactionmediated by spin fluctuations can be constructed usingthe renormalized vertex Γ in Eq. (13) and Fig. 3(c). Wenote that Γ contains both bubble-type and ladder-typediagrams and more if it is written with U , since the inter-action vertex γ is fully antisymmetrized. Linearizing theBethe-Salpeter equation for the static susceptibility as inthe case of superconductivity, we obtain the eigenvalueequationˆ L α q ψ = λ DW ψ, ( ˆ L α q ) kk ′ = − TN G ω k G ω, k + q Γ ′ αkk ′ . (35)The irreducible vertex part Γ ′ is given in terms of Γ asΓ ′ ch kk ′ = −
12 (3Γ sp + Γ ch ) ω,ω ; ω ′ − ω, k ′ − k , (36)Γ ′ sp kk ′ = 12 (Γ sp − Γ ch ) ω,ω ; ω ′ − ω, k ′ − k . (37)Figure 12(b) shows a diagram corresponding to the ver-tex Γ ′ . The matrix ˆ L in Eq. (35) is non-hermitian sothat the eigenvalues λ DW are complex numbers in gen-eral. By numerical calculations, we found that λ DW con-sists of purely real numbers as well as complex numbers.We have computed eigenvalues which has the largest realpart by means of the Arnoldi method [67]. As in the caseof superconductivity, we obtained 5 eigenvalues with dif-ferent spatial symmetries for each channel α = sp, ch.Figure 13 shows thus obtained eigenvalues for q = Q as a function of temperature. The largest fluctuation isthe ordinary antiferromagnetism (sp-A ) as expected. λ D W T sp A1gA2gB1gB2gEuch A1gA2gB1gB2gEu Figure 13: (Color online) Temperature dependence of theeigenvalues λ DW of Eq. (35) with q = Q . The parametersare U = 8 and δ = 0 . The cusp around T = 0 . T side and acomplex number on the low- T side. On the other hand,at low temperatures, there is a significant enhancementof the leading eigenvalue corresponding to the eigenfunc-tion with B (d x − y ) symmetry in the charge channel.This fluctuation corresponds to the staggered flux stateor the d-DW state. This eigenvalue may exceed the onein the AFM spin channel at lower temperatures. An ex-trapolation of λ DW to lower temperature indicates thatthe transition to the d-DW state takes place at T ∼ . VIII. DISCUSSIONS AND SUMMARY
We have applied the dual-fermion approach to thetwo-dimensional Hubbard model. The AFM fluctuationshave been taken into account by the ladder diagrams,which constitute the leading correction to the DMFTin terms of a 1 /d expansion. Practically, there was aconvergence issue which prevents calculations below themean-field critical temperature. By solving this technicalproblem, we are able to extend the applicability of theapproach to significantly lower temperatures.Possible phase transitions have been investigated fromsusceptibilities in both the particle-hole and particle-particle channels. For both calculations, we used aphysically equivalent irreducible vertex representing spin-fluctuation mediated interactions. Thus, we comparedfluctuations of spin/charge DW and singlet/triplet su-perconductivity for all spatial symmetries including d-SC and d-DW. The conclusion obtained is summarized1in the phase diagrams in Section II. The leading insta-bility under doping is the d-SC as expected. The d-DWfluctuations also show a clear tendency toward divergenceat low temperatures. However, the estimated transitiontemperature for the d-DW is below the T c of d-SC, forall parameters considered. Our result hence supports theabsence of a d-DW, which was yielded by several approx-imations [8–11].At temperatures above T c and low doping, we ob-served phase separation between the Mott insulator with δ = 0 and metallic region with δ = 0. The existenceof the phase separation agrees with other numerical cal-culations [17–19, 22, 23], but conflicts with QMC re-sults [24, 25]. Provided that the instability is not anartifact, a possible reason for the discrepancy is the sys-tem size: We employ N = 32 ×
32 = 1024 lattice sites,while the QMC calculations were performed for smallersize,
N < δ ≃ .
15 for U = 8,and therefore a pure d-SC occurs only in the limited re-gion 0 . . δ . .
18. This estimation is in quantitativeagreement with the cluster DMFT [18].The dual-fermion approach is complemental to thecluster DMFT among theories based on DMFT. Theladder approximation in the dual fermion, in particu-lar, aims at incorporating long-range fluctuations, whilethe cluster DMFT incorporates only short-range corre-lations. Hence, it would be informative to summarizeconsistency and inconsistency between those results toclarify characteristics of two complemental approaches.The instabilities reported from the cluster DMFT areconsistent with ours: The phase separation as well asd-SC take place under doping [17–19] and the d-DW ispredominated by the d-SC [9]. What can be reproducedby the dual-fermion approach but not by cluster DMFTis the critical behavior of the susceptibilities [47], sincea feedback of low-energy two-particle excitations to theself-energy is essential for it. Our results for the AFMsusceptibility exhibit a strong departure from the Curie-Weiss law, and are consistent with a critical temperatureof T N = 0 expected from the Mermin-Wagner theorem.This aspect will be considered in more detail elsewhere.The short-range correlations, on the other hand, playan important role near the Mott insulator. Its influencemay arise in the doping dependence of T c . In our dual-fermion calculations, it turned out that T c computed byneglecting the phase separation, namely, computed withthe thermodynamically unstable solution, show no down-turn as approaching the Mott insulator from finite dop-ing. In contrast, the cluster DMFT yields the dome shapeof the d-SC phase [29, 31].Further development beyond the dual-fermion ap-proach has recently been attempted [68]. The so-calleddual boson theory introduces a bosonic counterpart of thedual fermion for the purpose of treating intersite inter-actions beyond mean-field theory [69] and collective ex-citations [63, 70]. Furthermore, we expect that the dualboson in the spin channel yields formation of a intersite singlet, resulting in a reduction of T c near the Mott insu-lator. This effect may be brought about by the couplingbetween spins and the vector bosonic field, which can betreated exactly by the recently developed algorithm [71]based on the CT-QMC method. Acknowledgments
We acknowledge useful discussions with Y. Kuramoto,H. Yokoyama, H. Tsunetsugu, N. Tsuji, and M. Kitatani.A part of the computations was performed in the ISSPSupercomputer Center, the University of Tokyo. Two ofus (JO and HH) acknowledge hospitality of the ISSP dur-ing the NHSCP2014 workshop. This work was supportedby JSPS KAKENHI Grant Number 26800172. HH ac-knowledges support from the FP7/ERC, under GrantAgreement No. 278472-MottMetals.
Appendix A: Stabilization of self-energy calculation
In two-dimensional systems, the AFM susceptibilitydiverges exponentially for T approaching zero [50, 51].Since the ladder approximation explicitly takes the AFMfluctuations into account in the self-energy, the criticallylarge fluctuations complicate the convergence of the self-energy iterations. In this Appendix, we present how torelieve this difficulty to get better convergence.To show our idea of how to avoid the instability inthe critical regime, we begin with the FLEX equations,which encounter the same problem in a simpler form. InFLEX, the susceptibility χ ( q ) is given by χ ( q ) = χ ( q )1 − U χ ( q ) . (A1)Here, χ ( q ) is computed with the dressed Green’s func-tion G ( k ), which is determined self-consistently. Eventhough χ ( q ) is positive and finite in the converged solu-tion, the right-hand side may diverge or become negativeduring the iteration if a trial G ( k ) is not sufficiently closeto the solution. The self-energy evaluated from this sus-ceptibility is divergent. This is the source of instabilityof the iteration.To avoid this instability, we manipulate χ ( q ) so that χ ( q ) does not diverge. Specifically, we replace χ ( q ) with χ ′ ( q ) = min[ χ ( q ) , (1 − η ) /U ] , (A2)where η is a small constant. If η is too small, say η = 10 − , the instability may not be taken away. Empiri-cally, the FLEX iteration is stable down to η ≃ − . Wenote that this replacement is done only to avoid the in-stability and to approach the solution. If the trial G ( k ) issufficiently close to the actual solution, this replacementis no longer necessary. The converged solution thereforeis well defined. In order to assure this, we should checkwhether the condition U χ < − η is satisfied after the2iteration is converged. If this is not the case, it meansthat the actual solution of the equation is in the region U χ > − η .We apply the above trick to the dual self-energy inEq. (15). Introducing a matrix notation for the fermionicfrequencies, (ˆΓ ν q ) ωω ′ ≡ Γ ωω ′ ; ν q , ( ˆ χ ν q ) ωω ′ ≡ ˜ χ ω ; ν q δ ωω ′ and ( ˆ V αν q ) ωω ′ ≡ V αωω ′ ; ν q , and omitting indices α , ν and q for simplicity, the renormalized vertex Γ in Eq. (13) andthe effective interaction V in Eq. (16) are rewritten insimpler forms ˆΓ = ˆ γ + T ˆ γ ˆ χ ˆΓ , (A3)ˆ V = T ˆ γ ˆ χ [2ˆΓ − ˆ γ ] . (A4)We diagonalize the dimensionless matrix ( T ˆ γ ˆ χ ) accord-ing to ˆ U − ( T ˆ γ ˆ χ ) ˆ U = λ. (A5)We note that the eigenvalues λ i are complex in general.Using the diagonal matrix λ , the vertex ˆΓ and ˆ V areexpressed as ˆΓ = ˆ U (1 − λ ) − ˆ U − ˆ γ, (A6)ˆ V = ˆ U λ (1 − λ ) − (1 + λ ) ˆ U − ˆ γ. (A7)It is clear from these expressions that the iteration be-comes unstable once one of the eigenvalues λ i exceeds 1(or more precisely, Re λ i > λ i withRe λ ′ i = min(Re λ i , − η ) , (A8)during the iterations. We remind that one needs to verifythat the condition Re λ i < − η is fulfilled after conver-gence is reached. Appendix B: Spatial symmetry of pairingcorrelations
In this Appendix, we present how to obtain eigenvec-tors with specific symmetry in Eqs. (30) and (35). Ina power method and related algorithms, one computesa matrix-vector product ˆ Kφ (old) to obtain a new vector φ (new) . At this point, we restrict φ (new) to a subspacedefined by the projection operator P : φ (new) = P ˆ Kφ (old) . (B1) Thus, eigenvectors in the subspace are selectively com-puted.We consider an explicit expression for the projectionoperator P . In the square lattice, engenvectors φ be-long to one of the irreducible representations D in thepoint group D . The symmetry property of D is summa-rized in the character table in Table II. There are 5 irre-ducible representations, D = A , A , B , B , E u , and5 types of symmetry operations, C = E, C , C , C ′ , C ′′ . Table II: The character table for the point group D [72]. Forthe two-dimensional representation E u , the diagonal elementof the representation matrix is shown. The operation C de-notes π/ z axis, and C ′ and C ′′ denote π rotations around the x axis and the line x = y , respectively. E C C C ′ C ′′ A − − − − − − u ( x ) 1 0 − u ( y ) 1 0 − − The value σ = +1 or − C , i.e., C φ = σφ , while σ = 0 means thatthe operation changes the basis [e.g., π/ C transforms E u ( x ) to E u ( y )].We can project an arbitrary vector φ onto the ir-reducible representation D by enforcing the symmetryproperty given in Table II. Hence, the projection opera-tor P ( D ) may be decomposed into a product P ( D ) = Y C Q ( C , σ ( C , D )) . (B2)The operator Q ( C , σ ) singles out vectors which have theeigenvalue σ of the operation C . Since the eigenvalue σ is either +1 or −