Effect of interlayer processes on the superconducting state within t-J-U model: Full Gutzwiller wave-function solution and relation to experiment Michał Zegrodnik1,∗ and Józef Spałek2,1,† 1Academic Centre for Materials and Nanotechnology, AGH University of Science and Technology, Al. Mickiewicza 30, 30-059 Krakow, Poland 2Marian Smoluchowski Institute of Physics, Jagiellonian University, ul. Łojasiewicza 11, 30-348 Krakow, Poland (Dated: 13.09.2016) 7 1 The Gutzwiller wave function solution of the t-J-U model is considered for the bilayer high- 0 TC superconductor by using the so-called diagrammatic expansion method. The focus is on the 2 influence of the interlayer effects on the superconducting state. The chosen pairing symmetry is n a mixture of dx2−y2 symmetry within the layers and the so-called s± symmetry for the interlayer contribution. The analyzed interlayer terms reflect the interlayer electron hopping, the interlayer a exchangecoupling,andtheinterlayerpairhopping. Theobtainedresultsarecomparedwithselected J experimental data corresponding to the copper-based compound Bi-2212 with two Cu-O planes in 0 theunitcell. Forthesakeofcomparison, selected resultsforthecaseofthebilayerHubbardmodel 2 are also provided. This paper complements our recent results obtained for the single-plane high temperaturecuprates (cf. Ref. 26). ] n o PACSnumbers: 74.78.Na,84.71.Mn c - pr I. INTRODUCTION alyer pairing leads to deviations from the pure d-wave u character of the superconducting order parameter what s The crystal lattice of high-T cuprate superconduc- may be important in the context of some of the ARPES . C t tors is composed of blocks of Cu-O planes separated by experiments indicating that the gap symmetry is more a m charge-reservoir layers1. In the case of compounds with complex than a simple d-wave17–19. two (or more) Cu-O planes within a single block, the so- Another effect analyzed within the bilayer models is - d called bilayer splitting of the Fermi surface appears due the interlayer pair tunneling. It has been shown that n to the hybridization of electron states originating from the initially proposedinterlayerJosephsoncoupling aris- o theindividuallayers. Initially,itwaspredictedthatsuch ing as a second-order process in the interlayer electron c splitting should vanish in the nodal (k = k ) direction. hopping20 is insufficient to enhance significantly super- [ x y However, the ARPES analysis for Bi-22123 and YBCO4 conductivity due to a small value of the interlayer hop- 1 compounds has shown a nonzero splitting between the ping tz with respect to the intralayer value t (such cou- v hybridized bonding and antibonding bands in that di- pling would be ∝ t2/t). The pair tunneling term has z 2 rection. Possible explanation of appearance of such a also been introduced in a phenomenological manner to 4 splitting can be the nodal interlayer coupling5, mixing postulate the form of the in-phase gap function between 7 5 with the chain bonds6, or sensitivity of the quasiparti- two layers21. A different approach is to investigate the 0 cle spectrum in the nodal region to small effects such as interlayer pair hopping as originating from the matrix . the spin-orbit coupling7. It has also been argued that elements of the long-range Coulombic interaction12. 1 the main contribution to the nodal splitting comes from Here we analyze the influence of interalyer processes 0 7 the vertical hopping between the O 2pσ orbitals3. At such as the interlayer electron hopping, the interlayer 1 the same time, with the increasing number of layers in a pair hopping, and the interlayer exchange coupling, on : single Cu-O block the superconducting critical temper- the bulk superconducting state within the t-J-U model v ature is increased (for the number of layers n 6 3) in- ofthebilayerstructure. Onepfthepurposesofsuchanal- i X dicating that the interlayer processes influence the pair- ysis is to see to what extent the universal properties of r ing strength. Also, the optical Josephson plasma modes the single plane systems are preserved. To take into ac- a measurements2 show that the interlayer Josephson cou- countthestrongelectroniccorrelationsinthesystem,we pling strength can be correlated with the value of the usethediagrammaticexpansionGutzwillerwavefunction critical temperature, T . method (DE-GWF)22–25. The choice of the model has C Even though a vast majority of previous theoretical beendictatedby agoodquantitativeagreementbetween studies regarding the copper based materials is focused the theoretical results and principal experimental data on a single Cu-O plane, the bilayer Hubbard, t-J, and t- for the copper-basedmaterials, which has been reported J-U models have also been considered8–14. What distin- recently for the single layerversionof the same model26. guishes the latter models from the single-layer approach The gap symmetry selected in our analysis is a mixture is, among others, the appearance of the interlayer pair- of dx2−y2 symmetry within the layer and the so-called ing channels which can contribute to the total supercon- s± symmetry for the interlayer contribution9,27,28. The ducting gap. In this respect, a possible choice of pairing resulting gap function has the form ∆(k)=∆ (k)±∆ d ⊥ symmetries are discussed in Refs. 15 and 16. The inter- with the + (−) sign, corresponding to the bonding (an- 2 tibonding) band and ∆ (k) representing the intralayer twotermsoftheintralayerHamiltonian(2)representthe d d-wave gap. As reported in Refs. 27 and 28, such a Hubbardmodelwhich consistof the hopping and the in- state has been found as stable for the case of the t-J trasite Coulomb repulsion terms. The third term rep- model within the slave-boson analysis. In contradistinc- resents the antiferromagnetic exchange interaction. We tiontotheseresults,thecoexistenceoftheinterlayerand takeintoaccountthenearestandthenext-nearestneigh- intralayerpairingshasnotbeenfoundwithinofthevaria- bor hopping terms with t = −0.35 eV and t′ = 0.25|t|, tionalMonteCarlo(VMC) methodforthesamemodel9. respectively (unless stated otherwise). Moreover,the J ij Within our analysis of the bilayer t-J-U model, the ad- integralsareassumednonzeroonlyforthenearestneigh- mixtureoftheinterlayers± contributionleadstoasmall bors, for which J ≡ J ≡ 0.25|t|. The intraatomic ij lowering of the system energy with the interlayer gap Coulombrepulsion is set to U =14|t|, unless stated oth- magnitude being one order of magnitude smaller than erwise. The model can thus be considered as either the the intralayer one. As shown earlier25 for the case of a extendedHubbardort-J HamiltonianwithJ originating single-layer t-J model, our results are of the VMC qual- from p-d superexchange and a moderately high value of ity. U in these charge-transferMott-Hubbard systems26. For the sake of completeness, we present here also the The first two terms of the interlayer Hamiltonian Hˆ ⊥ results for the bilayer Hubbard model and analyze the represent the interlayer hopping and exchange, respec- influence of the differentiation between the values of the tively. The hopping parameters t⊥ are chosen so that in ij Hubbard U corresponding to the two layers. reciprocal space the single-particle part takes the form The structure of the paper is as follows. In Section II adnesdcrIiIbIewtheesbhaoswicecxopnlciceiptltybtehheinmdothdeelDHEa-mGiWltoFnimanethanodd H⊥t = ′′(cid:20)tbs+ t4z(coskx−cosky)2(cid:21)cˆ†klσcˆkl′σ , (4) kXll′σ asappliedtothebilayerstructure. InSectionIVwepro- vide the results ofour calculations andanalyze them. In where t introduces the k-independent bilayer splitting, bs particular,we compare our results with available experi- whereas the k-dependent part is taken in the usually mental data for the dispersion relation, the Fermi veloc- form6,29 and is responsible for a highly anisotropic char- ity, and the bilayer splitting, all in the nodal direction. acter of the splitting. In our analysis, we assume that The summary is provided in the last Section. Note also J⊥ is nonzero for i = j only. The remaining two terms ij that the results supplement our previous analysis for a of (3) represent the so-called on-site and off-site inter- single-layer situation26 that should be regarded as Part layer pair hoppings12, respectively. They originate from I of the series. In that paper the methodology behind the long-range Coulomb interaction. Follwing Ref. 12 selecting the t-J-U model is also discussed. we set U′ = −2U′′. In our analysis we take into ac- count only the nearest-neighbor interlayerpair hoppings (for |R − R | = 1). Nevertheless, the interlayer pair II. MODEL i j hopping generates four-site terms (two sites on the first layer and two on the second) and it is very difficult to We consider the t-J-U model for the bilayer high-TC take them into account within the DE-GWF method in superconductorrepresentedbythefollowingHamiltonian higherorders. Hence,theirinfluence isgoingtobe taken into account within the zeroth-order calculation scheme, Hˆ =Hˆ +Hˆ , (1) k ⊥ which is equivalent to the modified Gutzwiller approxi- mation (SGA)30–32. Such an approachis justified by the where fact the U′ and U′′ interaction parameters do not repre- Hˆ = ′t cˆ† cˆ +U nˆ nˆ + ′J Sˆ ·Sˆ , sent the dominant energies in the system (i.e., they are k ij ilσ jlσ il↑ il↓ ij il jl much smaller than U). Xijlσ Xil Xijl We have checked that the interlayer density-density (2) ′′ represents the intralayer part and Coulomb interaction term, Vll′nˆilnˆil′, has a very ill′ smallnegativeinfluenceontXhepairedstate. Thatiswhy ′′ ′′ Hˆ⊥ = t⊥ijcˆ†ilσcˆjl′σ + Ji⊥jSˆil·Sˆjl′ wedo notinclude this terminthe analysispresentedbe- iXjll′σ Xijll′ low. TheinfluenceoftheintralayerintersiteCoulombre- ′′ pulsionisalsoomittedhere;itsinfluenceontheSCstate +U′′ cˆ†il↑cˆ†jl↓cˆjl′↓cˆil′↑+cˆ†il↑cˆ†jl↓cˆil′↓cˆjl′↑ (3) has been discussed by us separately within the single- Xijll′ (cid:0) (cid:1) layer model in Ref. 33. ′′ +U′ cˆ†il↑cˆ†il↓cˆil′↓cˆil′↑, A methodological remark is in place here. Namely, in Xill′ principle the last term in (3) should also have a corre- ′ spondingtermin(2) ∼J cˆ† cˆ† cˆ cˆ . However, describes a number of possible interlayer dynamical pro- ijl il↑ il↓ jl↓ jl↑ X cesses. The primed (double primed) summation means it is of the order Jd2, where d2 = hnˆ nˆ i is the double i↑ i↓ that i 6= j (and l 6= l′), where i,j are indices of sites occupancy probability, and under the circumstance that within the layer and l,l′ enumerate the layers. The first J ≪U, as well as J is substantially smaller than |U′| it 3 canbesafelyneglected. Notealsothatthepresentformu- (l ,m )6=(j,l′) for all p, p′. The powers indices k (k ) p p 1 2 lationisbasedonreal-spacepairing,whereastheoriginal expresshowmanytimes the indices f onthe righthand p formulation was based on a direct Josephson tunneling of the equation have the value 1 (2) for a given term. of the Cooper pairs in k-space20,34,35. The coherentpair They fulfill the relation k +k = k. The maximal k, 1 2 tunneling across the planes has also been confirmed36. for which the terms in Eq. (9) are taken into account, represents the order of the expansion. By using the Wicks theorem in real space one can III. METHOD express all the averages in the noncorrelated state (h...i ≡ hΨ |...|Ψ i), which appear in the formulas for 0 0 0 To take into account the strong electronic correlations the ground state energy of the system in the form of di- we start with the Gutzwiller-type wave function in the agrams, with vertices attached to the lattice sites and form the edges corresponding to the so-called paramagnetic |ΨGi≡Pˆ|Ψ0i=Yil Pˆil|Ψ0i, (5) lScionijuelnls′,t≡Pliinjhelcˆls′†iσl↑focˆ≡r†jl′wh↓cˆih†i0liσ.cchˆIjnl|′∆σoiRu0,r|o=cral|tcRhueillas−tuiopRnesjrlc′wo|ne≤dtu3ac.kteinigntloinaecs-, When it comes to the SC lines we assume the d- where |Ψ i represents an uncorrelated state (chosen as 0 wave spin-singlet pairing with no intrasite pairing com- the superconducting uncorrelated state) and ponents, which corresponds to the strong correlation Pˆ ≡ λ |Γi hΓ|, (6) limit for the copper-based superconductors. The ap- il il ilil proach has been discussed extensively before24–26,37,38. XΓ For the interlayer pairing contribution we take the s± with λ being the variational parameters and |Γi rep- pairing symmetry9,27,28. il il resenting the states from the local basis The method described briefly above allows to express the expectation value of our Hamilotnian in the corre- |Γiil ∈{|∅iil,|↑iil,|↓iil,|↑↓iil}. (7) lated state The consecutive states representempty, singly, and dou- hHˆi = hΨG|Hˆ|ΨGi, (10) G bly occupiedlocalconfigurations,respectively. Tosignif- hΨ |Ψ i G G icantlysimplify the calculationsandimprovethe conver- gence, one customarily imposes the condition22,23 as a function of all the lines Pijll′σ, Sijll′ and the varia- tional parameters x , x . 1 2 Pˆ2 ≡1+x dˆHF , (8) Tocalculatethevaluesoftheparamagneticandsuper- il il il conducting lines, as well as the variational parameters, where x is a variationalparameter and dˆHF =nˆHFnˆHF, the grand-canonial potential F = hHˆiG−µGnG is min- il il il↑ il↓ imized, where µ and n are respectively the chemical nˆHF = nˆ − n , with n = hΨ |nˆ |Ψ i. By using G G ilσ ilσ l0 l0 0 ilσ 0 potential and the number of particles per lattice site de- Eqs. (6) and (8), one can express the parameters λ il termined in the correlated state. The minimization con- with the use of x . Furthermore, by assuming spatial il dition together with the normalization of |Ψ i, can be homogeneity within the layer we are left with the two 0 cast into the form of an effective Schrödinger equation variational parameters x (l = 1,2). For the case with l twoequivalentlayers,whichisanalyzedindetailhere,we Hˆ |Ψ i=E |Ψ i, (11) eff 0 eff 0 can additionally set x =x . However,at the end of the 1 2 paperwealsoconsiderthe casewithdifferentinteraction with the following effective Hamiltonian parameters for the two layers. The approach allows to express the expectation values, in the correlated state Hˆeff = teijffll′cˆ†ilσcˆjl′σ + ∆eijffll′cˆ†il↑cˆ†jl′↓+H.c. , |Ψ i, corresponding to each of the terms in the starting iXjll′σ Xijll′(cid:0) (cid:1) G (12) Hamiltonian(1). Forexample,forthecaseofthehopping where the effective hopping and the effective supercon- termthe correspondingexpectationvalue takesthe form ductinggapparametersaredefinedthroughtherelations hΨG|cˆ†ilσcˆjl′σ|ΨGi= ∂F ∂F ∞ k1! ′ xk11xk22hc˜†ilσc˜jl′σdˆHmF1f1...dˆHmFkfki0 , teijffll′ ≡ ∂Pijll′σ, ∆eijffll′ ≡ ∂Sijll′. (13) Xk=0 m1f1X...mkfk Such an approach has been commonly used in the (9) literature22,24,25,39 and as shown in Ref. 40, it is equiv- alent to the method based on the Lagrange multipliers where dˆH∅F ≡1, c˜(il†σ) ≡Pˆilcˆ(il†σ)Pˆil, and the index m corre- whichin turnguaranteesthe preservationof the statisti- sponds to lattice sites within the layer, whereas f labels cal consistency during the minimization procedure. the layers. The primmed summation on the right hand As one can see from (11) the uncorrelated state |Ψ0i side is restricted to (lp,mp)6=(lp′,mp′), (lp,mp)6=(i,l), is an eigenstate of Hˆeff, which has a BCS-like form with 4 the kinetic energy term and the pairing term leading to A. Results for the t-J-U model nonzero expectation values hcˆ† cˆ† i . Nevertheless, in il↑ jl′↓ 0 contrast to the BCS theory, here, the pairing is defined First, we analyze the influence of the interlayer elec- in real space and it is not due to the phonon-mediated tron hopping terms on the paired state, with no inter- interaction. Instead, it stems from the inter-electronic layer exchange interaction included (J⊥ =0). The mag- correlations taken into account in higher orders via dia- nitude of the k||-dependent part of the interlayer hop- grammaticexpansion. The terms for l′ 6=l in the second ping is governed by t (cf. Eq. (4)) and is responsible z summation of (12) corresponds to the interlayerpairing. for the anisotropy of the bilayer splitting. This term To determine |Ψ i we first transform the effective leads to a significant splitting in the antinodal direction 0 Hamiltonianto reciprocalspaceandthendiagnonalizeit (π,0) and no splitting in the nodal one (π,π). In some by using the Nambu representation. This allows for the papers referring to the bilayer cuprate superconductors, derivation of the self-consistent equations for the super- the nodal bilayer splitting is neglected, i.e., the situa- conducting and paramagnetic lines, as well as the chem- tion with tbs =0 is considered11,12,21. At the same time, ical potential. Additionally, the parameters {x} are de- considerations with the nearest-neighbor interlayer hop- l terminedvariationally. AfterthedeterminationofPijll′σ, ping included only (with tz = 0) have been also carried Sijll′, µ, x1, x2, we can calculate the correlated equiv- out8–10,13. In the latter situation the splitting between alents of the superconducting lines, hcˆ† cˆ† i , which thebondingandtheantibondingbandsisisotropicwhich il↑ jl′↓ G is rather loosely connected to the electronic structure of are regarded as the true gap parameters in the corre- the copper-based compounds. In Fig. 1 we show the lated state. In general, each non-correlated gap pa- rameter hcˆ† cˆ† i taken into account within the dia- influence of both the isotropic and anisotropic bilayer il↑ jl′↓ 0 splittings on the paired state. As one can see, the ef- grammatic procedure has its correlated correspondent hcˆ† cˆ† i . The symmetry of the correlatedgap reflects fect of increasing either tbs or tz is similar. Due to the il↑ jl′↓ G hoppingofelectronsbetweenthelayerstheintralayercor- the symmetry assumed for the non-correlated gap (in- related gap diminishes mainly in the overdoped regime, tralayer d-wave and interlayer s± wave). As it has been whereas the interlayer one increases. This is caused by pointed out in Refs. 24, 25, 33, the dominant contribu- thefactthatifthereisnointerlayerinteraction(J⊥ =0), tiontothesuperconductingphasecomesfromhcˆ† cˆ† i il↑ jl↓ G the hoppingmainlydilutes theintralayerpairbondsand with Ri−Rj =(a,0) for the d-wave intralayer pairing. thus weakens the pairing. Nevertheless, in the under- For the s± interlayer pairing the dominant contribution dopedregimethe gapparameter∆|| remains almostun- comes from hcˆ†il↑cˆ†jl′↓iG for Ri −Rj = (0,0) and l 6= l′. changed by the interlayer hopping.GAlso, the interlayer Those two intra- and inter-layer correlated gap parame- gapisapproximatelyoneorderofmagnitudesmallerthan ters are going to be denoted by ∆|| and ∆⊥ throughout the intralayer one. As one would expect, the case with G G the paper and lead to the following k-dependence of the t = t ≡ 0 leads to intralayer pairing only (∆⊥ = 0). z bs G overallgap In the zeroth order of the DE-GWF calculations, which are equivalent to the RMFT approach, one obtains no interlayer pairing even in the case of nonzero interlayer ∆ (k)=2∆||(cosk +cosk )±∆⊥, (14) G G x y G hopping. Next, we scrutinize the influence of the interlayer ex- where the ± sign corresponds to bonding and antibond- change coupling J⊥ on the SC gap. The results are pre- ing hybridized bands, respectively. sented in Fig. 2. It can be seen that the J⊥ term has a verysmallinfluenceontheintralayerpairingwhileitsin- terlayer correspondent is significantly increased (almost twice) when changing the value of J⊥ from 0.0 to 0.25. It should be noted that even for the case when J⊥ = J IV. RESULTS AND DISCUSSION (i.e., when the intra- and inter-layer exchange couplings areequal),theintralayergapparameterisstillaboutone Ouranalysisisperformedinthefourthorderofthedi- order of magnitude larger that the interlayer one. This agrammaticexpansion,ifnotstatedotherwise. Asshown effect should be related to the fact that the system is in Refs. 25 and 33, for the considered types of models, infinite in the (x,y) plane whereas it is limited in the satisfactory convergence is achieved in the 4-5 order of z-direction (two layers only). theexpansion. Wefocusmainlyonthet-J-U bilayercase In Fig. 3 we also show the influence of the interlayer analysis (Subsection A). As it has been already stated pair hopping processeson the superconducting gap. Fol- we set the intralayer model parameters to t=−0.35 eV, lowing Nishiguchi et al.12, we set U′ = −2U′′. As one t′ = 0.25|t|, U = 14|t|, J = 0.25|t|, unless stated other- could expect, the inclusion of the pair tunneling pro- wise. Inthefollowingweprovidethe valuesofthemodel cess enhances the intralayer SC gap and broadens the parameters in the units of |t|. For the sake of complete- SC stability regime. However, due to the mentioned ness we also present selected results for the case of the effect the maximum value of the interlayer gap is de- Hubbard model (Subsection B). creased. The positive effect on the interlayer correlated 5 (a) (b) (a) (b) 0.06 0.015 0.06 0.2 G 0.054 ap,G0.045 tbs=0.0 ap,G 0.01 tbs=-0.75 gap,0.045 G0.052 correlatedg0.00.1053 tbs=-0.75 correlatedg0.005 tbs=-0.25 correlated00.0.0135 GGGG ((((JJJJ====0000....0202)5)5)) (cG) 00..000186 0.2 0 0 0 0 0.15 0.3 0.45 0 0.15 0.3 0.45 0 0.15 0.3 0.45 0 0.15 0.3 0.45 doping, doping, doping, interlayercoupling,J (c) (d) FIG. 2. (a) Intra- and inter-layer correlated gap parameters 0.06 as a function of doping for two selected values of J⊥, for G tz=0.0 G tz =−0.3 and tbs =0. (b) and (c) the intra- and inter-layer p,0.045 p, 0.005 tZ=-0.45 correlated gap parameters, respectively, as a function of the edga 0.03 tz=-0.45 edga iδn=ter0la.2yewrheixchchiasncgleosceoutoploinpgtimfoarltdhoepsienlge.ctTedhevaelffueectofofdoJp⊥inigs correlat0.015 correlat0.0025 tZ=-0.15 mstaintoeraarendprtehsuersvtehde. principal single-plane features of the SC 0 0 0 0.15 0.3 0.45 0 0.15 0.3 0.45 (a) (b) doping, doping, 0.0045 0.08 U''=-0.9 FIG. 1. (a), (b) Dimensionless intra- (∆||) and inter-layer G G U''=0.0 (fsot∆er⊥Gpt)z-0c=.o2r5r0.e,lFaatonerddtbgtsbas=pvp0aartryhainmegientbteeerrtswl,abeyeoenrthg0aa.0psaaisnGfudzenr-co0t.i7ion5n,tohwfeidtwohphtionhlgee atedgap,00..0046 U''=0.0 atedgap,0.00.001053 U''=-0.9 dopingrange. (c),(d)theintra-andinter-layercorrelatedgap el el parameters, respectively, as a function of doping for tbs = 0 corr0.02 corr and tz varying between 0.0 and -0.45 with the step -0.15. 0 0 Similarlyasabove,theinterlayergapiszeroforthecasewith 0 0.15 0.3 0.45 0 0.15 0.3 0.45 tz =0. For nonzero tz both gaps vanish at the same doping. doping, doping, The results correspond to J⊥ =0. FIG. 3. Intralayer (a) and interlayer (b) correlated gaps as a function of doping for tz = −0.3, tbs = 0, J⊥ = 0 and for selected values of the pair hopping tunelling magnitude U′′ gapissignificanteventhoughthevaluesoftheU”param- varying between 0.0 and −0.9 with the step −0.3, assuming eterhavebeentakenasrelativelysmallincomparisonto U′ =−2U′′. The substantial value of U′′ 6−0.3|t| does not the Hubbard U. Note that the interlayer pair hopping reporoducetheexperimentalphasediagram with ∆G =0for term within this analysis originates from the long-range δ&0.4. Coulomb interaction12. From Figs. 1-3 one can see that the most important factor in changing effectively the gap is the Coulomb in- oftheinterlayerpairingmechanism. Suchachoiceisjus- teractions(U′,U′′)leadingtotheinterlayerpairhopping tified by the fact that the interlayer gap magnitude is processes. Theroleofthenumberoflayersintheelemen- one order of magnitude smaller than its intralayer core- tary cell in increasing the stability of the SC phase was spondant (cf. Fig. 1). In Fig. 4 (a) and (b) we plot the analyzed earlier within the mean-field approach41. calculatedtwosheetsofFermisurfaceandthedispersion In the remaining part of this subsection we show the relations close to the corresponding Fermi momenta for results for the case when the model parameters are set selectedvaluesofdopingclosetothe optimaldopingand to reproduce the selected principal characteristics mea- in the underdoped regime, respectively. The theoretical sured by ARPES experiments of Ref. 3. According to composed Fermi surface topology shown in Fig. 4 (a) is the experimental data obtained for Bi-2212, the bilayer roughlysimilartothatobtainedexperimentallyinRef. 3 splitting of the Fermi surface is highly anisotropic, with (cf. Fig. 1. inthatpaper)andthetheoreticallineswhich maximal value in the antinodal direction and small, but reflectthedispersionrelationsinFig. 4(b)arerelatively nonzero value in the nodal direction. Moreover, even close to the experimental points. However, it should be in the overdoped regime both the bonding and the anti- notedthatthe authorsofRef. 3provideonlythe critical bondingpartsoftheFermisurfaceappeartobehole-like. temperature of the corresponding samples which allows To reproduce the data within our model both t and t to determine the correspondingdoping regionsonly (un- z bs havetobetakenasnonzero(t =−0.1andt =−0.06). derdopedoroverdopedregion). Thespecificdopinglevels z bs Intheseconsiderationswehaveneglectedtheappearance ofthe samplesarenot providedexplicitly andthis is one 6 of the reasons why our theoretical results may not re- flectpreciselythosecomingfromtheexperiment. InFig. 4 (c) we plot the bilayer splitting in energy (top) and (a) =0.16 (b) (closetoopt.dop.) =0.125(underdopedreg.) in the reciprocal space (bottom) in the nodal direction, 3 0.01 antibonding both as a function of doping. The grayregionrepresents bonding 0 tinheRmefe.a3su.rTedhevaelxupeesrfimorednitffaellryenmtedaospuirnegdlveevreylswperaeksednotepd- ky2 [eV]-0.01 a(GWF) 1 -0.02 b(GWF) ing dependence of the nodal bilayer splitting value has a(exp.) been confirmed in Ref. 4 for YBCO. However, for this -0.03 b(exp.) 0 compound,thevalueof∆BS wasalmostfivetimeslarger 0 1 2 3 -0.01 0 0.01 0.02 than the one measured for Bi-2212. As one can see, our (c) kx (d) k-kF[A-1] theoreticalresultsareclosetotheexperimentalonesonly in the overdopedregime while in the underdopedregime V] 0.03 measuredvalues 2.2 tshheoucladlcbuelanteodtebdilathyeartsFpoliutrtniniegrseigtnaifil.c4a2nhtlayvdeercerpeoarsteesd. Iat [eBS0.015 GWF A] 2 0 V complete vanishing of the bilayer splitting below the op- e 1.8 timal doping for the YBCO samples. The inconsistence -1]A 0.016 measuredvalues v[F GWF(bonding) between the different sets of experimental data from the [S0.008 1.6 GWF(antibonding) B GWF measuredvalue one side, and the theoretical results from the other re- 0 1.4 quires a further analysis. 0 0.1 0.2 0.3 0 0.1 0.2 0.3 doping, doping, In Fig. 5 we show the effect of the interlayer coupling (e) (f) onthebilayersplittingandthenondalFermivelocity. As 0.4 G0.06 oadnlespeoecinnadcnernesaeteseo,snw.itHthhoewtJheve⊥eirnv,catrlhueeaesFiaennrgmdJie⊥vveeltnohcefiotbyrilirasayvteehrreysrpwslitetratoiknnlgyg -1k[]AF0.35 dgap, 0.04 iTwnhhteeicrhlsaayimsereaclssooiutupaallimtnigoosnsttthuaenkacehffsaepncltgaeecdseibnfyovrtFhthiaserefacbcotarorrerellay(ctevf.disiFgbialgep.. nodal 0.02.53 baonntibdoinngding correlate 0.020 0 0.1 0.2 0.3 0 0.1 0.2 0.3 0.4 2). doping, doping, Duetothenonzerosplittinginthenodaldirection,also two values of the nodalFermi velocity appear,whichare FIG. 4. (a) Fermi surface sheets close to the optimal dop- shown in Fig. 4 (d). The measured values reported in ing. (b) Theoretical (solid lines) and experimental (circles) dispersion relations in the underdoped regime. The exper- Ref. 3 are close to 2.0 eVÅ and they correspond to imental data have been taken from Ref. 3. (c) The nodal the bonding band. Also, other experimentalreports cor- bilayer splitting in energy (top) and in the reciprocal space respondingto the cupratesuperconductors,demonstrate (bottom), both as a function of doping. The shaded regions veryweakdopingdependenceofthenodalFermivelocity correspond roughly to the measured values reported in Ref. close to that value4,47,48. In this case, the weak doping 3. (d) Nodal Fermi velocity as a function of doping corre- dependenceisreproducedbythetheoryascanbeseenin sponding to the bonding and antibonding bands. The gray Fig. 4 (d). Moreover, the bonding and the antibonding (solid)linecorrespondstothemeasuredvaluetakenfromRef. Fermi velocities are very close. Near half filling both of 3 which refer to the bonding band. (e) Fermi momenta cor- them have almost the same value. It is caused by the respondingto thebondingand antibondingbands,both as a fact that the nodal bilayer splitting is decreasing as we function of doping. (f) Intralayer correlated gap magnitude as a function of doping. The presented fittings have been aregettingclosertothehalffilledsituation,whatisillus- obtained for the following values of the model parameters: tratedalsoinFig. 4(e). Forthesakeofcompleteness,in Fig. 4 (f) we also provide the correlated gap magnitude tu′n=its0.o3f6|,t|t)z.=−0.15, tbs =−0.06, U =13.5, J =0.25 (all in asa function ofdopingfor the same setofmodelparam- eters as the remaining results presented in this Figure. As one can see in Fig. 6, the effect of the splitting B. Results for the Hubbard model parametertbs issimilarasinthecaseofthet-J-U model. However,forlargevaluesoft ,adiscontinuityappearsin bs theδ-dependenceofboth∆|| and∆⊥ intheunderdoped Forthesakeofcompleteness,wehavealsoanalyzedthe G G regime (the curves for t = 0.75), which has not been selectedresultsforthecaseoftheHubbardmodel(J ≡0, bs J⊥ ≡ 0). For these calculations we set the interlayer observed earlier. hopping as non-dependent on k (t = 0). In the main We have also analyzed the influence of the dispropor- z partofthisSubsectionwealsoneglecttheinterlayerpair tionationbetween the Coulombrepulsionin the two lay- hopping term (U′ =U′′ =0). Nonetheless, the influence ers. In the following analysis U and U correspond to 1 2 of nonzero U′ and U′′ is shown briefly at the end. the Coulomb repulsion in the first and the second layer, 7 (a) (b) (a) (b) 0.02 0.045 2.1 []eVBS0.00.1053 J 0.25J 0.0 v[eV]AF1.29 baonntibdoinngd(inJg(J0.0)0.0) edgap,G0.00.4053 tbs=0.0 edgap,G0.01 tbs=-0.75 1.8 bonding(J 0.25) at at 0 1.7 antibonding(J 0.25) orrel 0.015 orrel 0 0.1do0p.2ing,0.3 0.4 0 0.1doping,0.2 0.3 c 0 tbs=-0.75 c 0 tbs=-0.25 0 0.1 0.2 0.3 0 0.1 0.2 0.3 doping, doping, FIG. 5. (a) The nodal bilayer splitting in energy as a func- tion of doping for different values of the interlayer exchange coupling varying between J⊥ = 0.0 and J⊥ = 0.25 with the FIG. 6. Intra- (a) and inter-layer (b)correlated gap vs. dop- step 0.05. (b) Nodal Fermi velocity as a function of doping ing for tz = 0 for tbs varying between 0.0 and 0.75 with the correspondingtothebondingandantibondingbandsfortwo step 0.25. For tbs =0 the interlayer gap is zero in the whole doping range. The data correspond to the Hubbard model values of thethe interlayerexchange coupling. with U =14 and (U′ =U′′ =0). respectively. Insuchasituationthe intralayercorrelated gaps, the average number of electrons per atomic site, (a) (b) andtheaveragedoubleoccupanciesinthetwolayerscan 0.06 0.06 have different values. All the mentioned quantities are G U2=14 G U2=9.8 p, 0.04 p, 0.04 presented in Figs. 7 and 8. It can be seen that the a a g g cboertrweelaentedUgaapnsdfoUr bointhcrelaayseerss. dFeocrreUase=as0thbeodthiffesurepnecre- corr. 0.02 GG corr. 0.02 1 2 2 0 0 conducting gaps vanish (cf. Fig. 8 (a)). In Fig. 8 (b) we (c) (d) 0.06 0.06 show that in the considered system electrons are pushed G U2=7 G U2=5.6 from the layer with stronger Coulomb repulsion to that p, 0.04 p, 0.04 a a with the weaker. Obviously, the double occupancy in g g the second layer increases with decreasing U2, as can orr. 0.02 orr. 0.02 c c be seen in Fig. 8 (c). For the case with U1 = 14 and 0 0 U = 0 the paired phase vanishes, ∆(1) = ∆(2) = 0 (e) 0.06 (f) 0.06 (c2f. Fig. 8(a)). However, after settiGng to inGterlayer G U2=4.2 G U2=2.8 pair hopping as nonzero, we can retain the stability of ap, 0.04 ap, 0.04 g g tthheispsiatiureadtiopnh,atsheeinpatirhiengsyosrtiegmina(tcifn.g fFroigm. t8he(dc)o)r.relIan- corr. 0.02 corr. 0.02 tion effects in the layer with nonzero Coulomb repulsion 0 0 0.1 0.2 0.3 0.1 0.2 0.3 istransferedtothesecondlayer(withnoCoulombrepul- doping, doping, sion) by the pair hopping process. Such situation brings into mind the systems composedof one superconducting FIG. 7. Correlated gaps for the bilayer Hubbard model in monolayer of correlated compound placed on a metallic thefirst(redsolid line)andsecond(bluesolid line)layersvs. layer in which the paired phase can be induced by the dopingfordifferentvaluesofCoulombrepulsioninthesecond proximity effect43. However, in such interpretation the layer, U2 (the explicit values are provided in the Figures). CooperpairtunnelingtermsincludedintheHamiltonian The value U1 corresponding to the Coulomb repulsion in the would have to originate from the Josephson coupling. first layer remains constant and equal to U1 = 14 while the interlayer electron hopping parameters are set to tbs = −0.1 andtz =0. In(a)wehaveU1 =U2 andbothcorrelatedgaps have the same values marked by blue line, while in (f) the V. SUMMARY value of U2 is so low that the correlated gap in the second layerispracticallyzero. Thedatacorrespondtonointerlayer We start with a general remark. The majority of the pair hoppingterm (U′ =U′′ =0). papers concerning the real space pairing in high-Tc su- perconductorsdealwithasingleCu-Oplaneasreflecting the intrinsic properties of the copper-based compounds. sic question still remains as to whether the model two- In our recent (preceding) paper26 we have shown that dimensionality, even if realistic, provides similar results indeed in such an approach one can rationalize, even in whenthe basicaspectsofthree-dimensionalityareincor- a quantitative manner, the principal ground-state prop- porated into the theoretical scheme. Our present work erties ofhigh-T superconductors. Nevertheless,the ba- provides an affirmative answer to this question, even C 8 (a) alsowithinouranalysisbothinterlayerhoppingtermsdo (b) not influence significantly the correlated gap amplitude, 0.03 =0.21 0.45 ∆||,inthe underdopedregimefornottoolargevaluesof G G0.015 G n |tbs| and |tz|. 0.3 n1 Themaininfluenceoftheinterlayerexchangecoupling G n 2 is seen inthe interlayersuperconducting gapmagnitude, 0 0 4 8 12 0 4 8 12 while the intralayer gap is only slightly reduced with in- Coulombrepulsion,U Coulombrepulsion,U 2 2 creasing J⊥ (cf. Fig. 2). Weak dependence of the in- (c) (d) 0.3 0.12 tralayerpairingonJ⊥ is alsomentioned inRef. 13. The nonzeroexchangeintegralJ⊥ explainsthenatureofanti- 2d0.15 dd1222 G0.08 GG fisertreommpatginngettiosmcoinnctluhdeeinthsualtaftoinrgnostnazteer4o5,a46n.dTsmhearlelfvoarelu,eist 0.04 of J⊥ one can rationalize the nature of AF state and at 0 0 thesametimeupholdthevalidityofthe principalresults 0 4 8 12 0 -0.2 -0.4 -0.6 concerning the superconducting phase obtained for the Coulombrepulsion,U interlayerpairhoppingU'' 2 single-plane systems. Among the considered interlayer processes the most FIG. 8. Intralayer correlated gap (a), the number of carriers significant influence on the SC phase has been reported per atomic site (b), and double occupancy probability (c) in the first (blue line) and second (red line) atomic layer, all fortheinterlayerpairhopping. Thecorrespondingterms as a function of Coulomb repulsion in the second layer. For havealsobeenincludedwithintheanalysisofthebilayer (a), (b), and (c) the interlayer pair hopping is set to zero Hubbard model in Ref. 12. Similarly as in our case, it is U′′ = U′ = 0. (d) Intralayer correlated gaps as functions necessary to include at least one of the off-site terms to of interlayer pair hopping for U2 = 0. We keep U′ = −2U′′ achieve a meaningful enhancement of the SC phase sta- fulfilled. For (a), (b), (c), and (d) the Coulomb repulsion in bility. Inoursituation,thisisadirectconsequenceofthe the first layer is set to U1 = 14 and the interlayer electron factthatweassumehcˆ† cˆ† i≡0,asaresultofthestrong hopping parameters are: tbs = −0.1 and tz = 0, while the il↑ il↓ selected value of doping is δ=0.21. on-siteCoulombrepulsion. Withintheanalysisshownin Ref. 12 the contribution to the pairing, resulting from the interplay between the three pair-hopping terms (one on-site and two off-site pair hoppings, cf. Eq. 3), seems though the detailed features are model-parameter de- nottobe straightforward. However,theiranalysisiscar- pendent. We believe that the direct relation between ried out within the spin-fluctuation picture for the case theory and experiment analyzed here and in our previ- of weak interactions, what means that their approach is ous work26, cancontribute significantly to the resolution very different from ours. of the fundamental question whether a purely electronic We have also analyzed the basic features which refer modelbasedonelectroniccorrelations(withoutaglue44) to the available data obtained by ARPES experiments canrationalizethesuperconductingpropertiesofthecop- (such as the Fermi surface two-sheet structure, the bi- per systems. layer splitting, the dispersion relations, the nodal Fermi Turning to our results here, we have analyzed the in- velocity). In this part of our analysis, the model param- fluence ofdifferent interlayerdynamicalprocessesonthe eters have been fitted and roughly reproduce the exper- superconducting state within the bilayer t-J-U model. imental data presented in Ref. 3. As one can see, the When it comes to the bilayer splitting, induced by the highly anisotropiccharacterof the bilayersplitting, with interlayer electron hopping, both the anisotropic (k||- a small value of the splitting in the nodal direction, has dependent term with the magnitude tuned by parame- been reproduced (cf. Fig. 4 (a)). The calculated value tertz)andtheisotropic(ofmagnitudetbs)contributions of the bilayer splitting is similar to that measuredin the havebeenanalyzed. Ouranalysisshowsthatbothterms experimentfor overdopedsamples (cf. Fig. 4(b)). How- have a similar influence on the supercodncuting gap pa- ever, in the underdoped regime significant discrepancies rameters (cf. Fig. 1). With the increasing interlayer between our theory and experiment appear (cf. Fig. 4 hoppingmagnitudesthe intralayerSC gapdecreasesand (c) and (e)). The mentioned report3 does not provide that of the interlayer part increases. For the case of the specific values of the doping levels, so the precise com- bilayerHubbardmodel,forlargevalueoftbs =−0.75the parison cannot by carried out. At this point, it should discontinuityappearsinthedopingdependenceoftheSC be noted that in disagreement with Ref. 3, Fournier et gap (cf. Fig. 6). As this effect has not been observed al.42havereportedavanishingbilayersplittingbelowthe experimentally, the tbs value must be essentially smaller optimal doping for the YBCO samples. Such a behavior than the value of |t|. is not reproduced within our approach. Finally, we also IthasbeenreportedrecentlyinRef. 13thattheresults show that the nodal Fermi velocities, corresponding to of the VMC calculations for the bilayer t-J-U model for the bonding andthe antibonding bands, are verysimilar d-wavesuperconductingphaseareinsensitivetot inthe and close to the value measured experimentally (≈ 2.0 ⊥ range |t | <0.5 for δ <0.2. As one can see from Fig. 1 eVÅ ) in the whole doping range (cf. Fig. 4 (d)). The ⊥ 9 presentanalysisforthebilayert-J-U modelcomplements can be induced by the Cooper pair tunneling processes our recent results26 for a single Cu-O planar supercon- (cf. Fig. 8 d). ductor. ThebilayerHubbardmodelconsideredheremightalso be ofinterestin referencetothe ultra-coldatomsystems The results presented for the case of the bilayer Hub- analyzed in recent years49,50. bardmodelwith the twodifferentvaluesofthe Coulomb repulsion for the two layers show that the dispropor- VI. ACKNOWLEDGEMENT tionation between U and U leads to different carrier 1 2 concentrations in the layers, what in turn has a nega- tive influence on the paired phase (cf. Fig. 8 a and b). We would like to acknowledge the financial sup- However,for the situation representinga metallic mono- port from the Grant MAESTRO, No. DEC- layerincontactwithastronglycorrelatedsuperconduct- 2012/04/A/ST3/00342fromtheNationalScienceCentre ing monolayer (U = 14 and U = 0), the paired phase (NCN) of Poland. 1 2 ∗ [email protected] 26 J.Spałek,M. Zegrodnik,andJ.Kaczmarczyk,Phys.Rev. † [email protected] B 95, 024506 (2017). 1 Shin-ichi Uchida, High Temperature Superconductivity, 27 M. Ubbens,P. A. Lee, Phys.Rev.B 50, 438 (1994). Springer Japan, Tokyo, 2015. 28 P. A. Lee, K. Kuboki, J. Phys. Chem. Solids 56, 1633 2 R.Kleiner,F.Steinmeyer,G.Kunkel,andP.Müller,Phys. (1995). Rev.Lett. 68, 2394 (1992). 29 O. K. Andersen, O. Jepsen, A. I. Lichtenstein, and I. I. 3 A. A.Kordyuket al., Phys.Rev.B 70, 214525 (2004). Mazin, Phys.Rev.B 49, 4145 (1994). 4 S.V.Borisenkoetal.,Phys.Rev.Lett.96,117004 (2006). 30 J. Jędrak and J. Spałek,Phys.Rev. B83, 104512 (2011). 5 D. Garcia-Aldea and S. Chakravarty, New J. Phys. 12, 31 J. Kaczmarczyk and J. Spałek, Phys. Rev. B 84, 125140 105005 (2010). (2011). 6 O. K. Andersen, A. I. Lichtenstein, O. Jepsen, and F. 32 M. Zegrodnik, J. Bünemann,and J. Spałek,NewJ. Phys. Paulsen, J. Phys. Chem. Solids 56, 1573 (1995). 16, 033001 (2014). 7 N.Harrison,B.J.Ramshaw,andA.ShekhterSci.Rep.5, 33 M. Abram, M. Zegrodnik, and J. Spałek, Arxiv: 10914 (2015). 1607.05399 (2016). 8 A. Medhi, S. Basu, and C. Kadolkar, Eur. Phys. J. B 72, 34 W.-Q. Chen, J. Y. Gan, T. M. Rice, and F. C. Zhang, 583 (2009). Europhys.Lett. 98, 57005 (2012). 9 T.A.MaierandD.J.ScalpinoPhys.Rev.B84,180513(R) 35 P. W. Anderson,Physics B 199& 200, 8 (1994). (2011). 36 J. D. Hettinger et al., Phys. Rev.Lett.74, 4726 (1995). 10 N. Lanata, P. Barone, and M. Fabrizio, Phys. Rev. B 80, 37 M. M. Wysokiński, J. Kaczmarczyk, and J. Spałek, Phys. 224524 (2009). Rev.B 92, 125135 (2015). 11 M.Mori,T.Tohyama,andS.Maekawa,J.Phys.Soc.Jpn. 38 M. M. Wysokiński, J. Kaczmarczyk, and J. Spałek, Phys. 75, 034708 (2006). Rev.B 94, 024517 (2016). 12 K. Nishiguchi, K. Kuroki, R. Arita, T. Oka, and H. Aoki, 39 Q.H. Wang, Z.D. Wang, Y. Chen, and F.C. Zhang, Phys. Phys. Rev.B 88, 014509 (2013). Rev.B 73, 092507 (2006). 13 K.-K.Voo, Phys. Lett. A 379, 1743 (2015). 40 J. Kaczmarczyk, Phil. Mag. 95, 563 (2014). 14 R. Rüger, L. F. Tocchio, R. Valenti, and C. Gros, New J. 41 K. Byczuk and J. Spałek,Phys. Rev.B 53, R518 (1996). Phys. 16, 033010 (2014). 42 D. Furnieret al., Nature Phys.6, 905 (2010). 15 A. Medhi, S. Basu, C. Y. Kadolkar, Phys. Rev. B 76, 43 O. Yuli, I. Asulin, O. Millo, D. Orgad, L. Iomin, and G. 235122 (2007). Koren, Phys.Rev.Lett. 101, 057005 (2008). 16 A. Medhi, S. Basu, C. Y. Kadolkar, Physica C 451, 13 44 P. W. Anderson,Science 317, 1705 (2007). (2007). 45 R.Coldea, S. M. Hayden,G. Aeppli,T. G. Perring, C. D. 17 G. Zhao, Phys.Rev.B 75, 140510(R) (2007). Frost,T.E.Mason,S.-W.Cheong,andZ.Fisk,Phys.Rev. 18 H. Dinget al., Phys.Rev.Lett. 74, 2784 (1995). Lett. 2001, 5377 (2001). 19 I. Vobornik et al., Physica C 317, 589 (1999). 46 S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. 20 S. Chakravarty, A. Sudbø, P. W. Anderson, and S. P. Rev.B 39, 2344 (1989). Strong, Science 261, 337 (1993). 47 X.J. Zhou et al.,Nature 423, 398 (2003). 21 C. Chen, A. Fujimori, C. Ting, and Y. Chen, arxiv. 48 A. A. Kordyuk,S. V. Borisenko, A. Koitzsch, J. Fink, M. 1211.3477 (2012). Knupfer,and Berger, Phys. Rev.B 71, 214513 (2005). 22 J. Bünemann, T. Schickling, and F. Gebhard, Europhys. 49 Y.Nishida, Phys.Rev.A 82, 011605(R) (2010). Lett. 98, 27006 (2012). 50 M. Klawunn, A. Pikovski, L. Santos, Phys. Rev. A 82, 23 F. Gebhard, Phys.Rev. B 41, 9452 (1990). 044701 (2010). 24 J. Kaczmarczyk, J. Spałek, T. Schickling, and J. Büne- mann, Phys.Rev.B 88, 115127 (2013). 25 J. Kaczmarczyk, J. Bünemann, and J. Spałek, New J. Phys. 16, 073018 (2014).