Effect of interlayer processes on the superconducting state within t-J-U model: Full Gutzwiller wave-function solution and relation to experiment

TL;DR

This paper investigates how interlayer effects influence the superconducting state in bilayer high-Tc cuprates using the full Gutzwiller wave-function approach within the t-J-U model, comparing results with experimental data and simpler models.

Contribution

It provides a detailed analysis of interlayer processes on superconductivity using the full Gutzwiller wave-function in the t-J-U model, extending previous single-layer studies.

Findings

Interlayer effects significantly modify the superconducting pairing symmetry.

Results align with experimental data for Bi-2212.

Comparison with bilayer Hubbard model highlights the importance of interlayer interactions.

Abstract

The Gutzwiller wave function solution of the $t$-$J$-$U$ model is considered for the bilayer high-T$_C$ superconductor by using the so-called diagrammatic expansion method. The focus is on the influence of the interlayer effects on the superconducting state. The chosen pairing symmetry is a mixture of $d_{x^2-y^2}$ symmetry within the layers and the so-called $s^{\pm}$ symmetry for the interlayer contribution. The analyzed interlayer terms reflect the interlayer electron hopping, the interlayer exchange coupling, and the interlayer pair hopping. The obtained results are compared with selected experimental data corresponding to the copper-based compound Bi-2212 with two Cu-O planes in the unit cell. For the sake of comparison, selected results for the case of the bilayer Hubbard model are also provided. This paper complements our recent results obtained for the single-plane high temperature cuprates [cf. J. Spa{\l}ek, M. Zegrodnik, and J. Kaczmarczyk, Phys. Rev. B {\bf 95}, 024506 (2017)]

Equations

$$\hat{\mathcal{H}}={\hat{\mathcal{H}}}_{\parallel }+{\hat{\mathcal{H}}}_{\perp },$$

$${\hat{\mathcal{H}}}_{\parallel }=\text{sideset}{}^{\prime }{\sum }_{ijl\sigma }{t}_{ij}{\hat{c}}_{il\sigma }^{†}{\hat{c}}_{jl\sigma }+U\sum _{il}{\hat{n}}_{il\uparrow }{\hat{n}}_{il\downarrow }+\text{sideset}{}^{\prime }{\sum }_{ijl}{J}_{ij}{\hat{\mathbf{S}}}_{il}\cdot {\hat{\mathbf{S}}}_{jl},$$

$$\begin{split}\mathcal{\hat{H}_{\perp}}=&\sideset{}{{}^{\prime\prime}}{\sum}_{ijll^{\prime}\sigma}t^{\perp}_{ij}\hat{c}^{\dagger}_{il\sigma}\hat{c}_{jl^{\prime}\sigma}+\sideset{}{{}^{\prime\prime}}{\sum}_{ijll^{\prime}}J^{\perp}_{ij}\hat{\mathbf{S}}_{il}\cdot\hat{\mathbf{S}}_{jl^{\prime}}\\ +&U^{\prime\prime}\sideset{}{{}^{\prime\prime}}{\sum}_{ijll^{\prime}}\big{(}\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl\downarrow}\hat{c}_{jl^{\prime}\downarrow}\hat{c}_{il^{\prime}\uparrow}+\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl\downarrow}\hat{c}_{il^{\prime}\downarrow}\hat{c}_{jl^{\prime}\uparrow}\big{)}\\ +&U^{\prime}\sideset{}{{}^{\prime\prime}}{\sum}_{ill^{\prime}}\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{il\downarrow}\hat{c}_{il^{\prime}\downarrow}\hat{c}_{il^{\prime}\uparrow},\end{split}$$

$$\mathcal{H}^{t}_{\perp}=\sideset{}{{}^{\prime\prime}}{\sum}_{\mathbf{k}ll^{\prime}\sigma}\bigg{[}t_{\textrm{bs}}+\frac{t_{z}}{4}(\cos k_{x}-\cos k_{y})^{2}\bigg{]}\;\hat{c}^{\dagger}_{\mathbf{k}l\sigma}\hat{c}_{\mathbf{k}l^{\prime}\sigma}\;,$$

$$|{\Psi }_G⟩\equiv \hat{P}|{\Psi }_0⟩=\prod _{il}{\hat{P}}_{il}|{\Psi }_0⟩,$$

$${\hat{P}}_{il}\equiv \sum _{\Gamma }{\lambda }_{il}|\Gamma {⟩}_{ilil}⟨\Gamma |,$$

$$|\Gamma {⟩}_{il}\in \{|\varnothing {⟩}_{il},|\uparrow {⟩}_{il},|\downarrow {⟩}_{il},|\uparrow \downarrow {⟩}_{il}\}.$$

$${\hat{P}}_{il}^2\equiv 1+x_{il}{\hat{d}}_{il}^{\text{HF}},$$

$$\begin{array}{rl}⟨{\Psi }_G|& {\hat{c}}_{il\sigma }^{†}{\hat{c}}_{j{l}^{\prime }\sigma }|{\Psi }_G⟩=\\ & \sum _{k=0}^{\infty }\frac{1}{k!}\text{sideset}{}^{\prime }{\sum }_{m_1{f}_1...m_k{f}_k}{x}_1^{k_1}{x}_2^{k_2}⟨{\tilde{c}}_{il\sigma }^{†}{\tilde{c}}_{j{l}^{\prime }\sigma }{\hat{d}}_{m_1{f}_1}^{\text{HF}}...{\hat{d}}_{m_k{f}_k}^{\text{HF}}{⟩}_0,\end{array}$$

$$⟨\hat{H}{⟩}_G=\frac{⟨{\Psi }_G|\hat{H}|{\Psi }_G⟩}{⟨{\Psi }_G|{\Psi }_G⟩},$$

$${\hat{H}}_{\text{eff}}|{\Psi }_0⟩=E_{\text{eff}}|{\Psi }_0⟩,$$

$$\hat{\mathcal{H}}_{\textrm{eff}}=\sum_{ijll^{\prime}\sigma}t^{\textrm{eff}}_{ijll^{\prime}}\hat{c}^{\dagger}_{il\sigma}\hat{c}_{jl^{\prime}\sigma}+\sum_{ijll^{\prime}}\big{(}\Delta^{\textrm{eff}}_{ijll^{\prime}}\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl^{\prime}\downarrow}+H.c.\big{)},$$

$$t_{ijl{l}^{\prime }}^{\text{eff}}\equiv \frac{\partial \mathcal{F}}{\partial P_{ijl{l}^{\prime }\sigma }},{\Delta }_{ijl{l}^{\prime }}^{\text{eff}}\equiv \frac{\partial \mathcal{F}}{\partial S_{ijl{l}^{\prime }}}.$$

$${\Delta }_G(\mathbf{k})=2{\Delta }_G^{||}(\cos k_x+\cos k_y)\pm {\Delta }_G^{\perp },$$

Full text

I Introduction

The crystal lattice of high-$T_{C}$ cuprate superconductors is composed of blocks of Cu-O planes separated by charge-reservoir layersUchida2015 . In the case of compounds with two (or more) Cu-O planes within a single block, the so-called bilayer splitting of the Fermi surface appears due to the hybridization of electron states originating from the individual layers. Initially, it was predicted that such splitting should vanish in the nodal ($k_{x}=k_{y}$) direction. However, the ARPES analysis for Bi-2212Kordyuk2004 and YBCOBorisenko2006 compounds has shown a nonzero splitting between the hybridized bonding and antibonding bands in that direction. Possible explanation of appearance of such a splitting can be the nodal interlayer coupling Garcia2010 , mixing with the chain bonds Andersen1995 , or sensitivity of the quasiparticle spectrum in the nodal region to small effects such as the spin-orbit coupling Harrison2015 . It has also been argued that the main contribution to the nodal splitting comes from the vertical hopping between the $O$ $2p_{\sigma}$ orbitalsKordyuk2004 . At the same time, with the increasing number of layers in a single Cu-O block the superconducting critical temperature is increased (for the number of layers $n\leqslant 3$) indicating that the interlayer processes influence the pairing strength. Also, the optical Josephson plasma modes measurementsKleiner1992 show that the interlayer Josephson coupling strength can be correlated with the value of the critical temperature, $T_{C}$.

Even though a vast majority of previous theoretical studies regarding the copper based materials is focused on a single Cu-O plane, the bilayer Hubbard, $t$-$J$, and $t$-$J$-$U$ models have also been consideredMedhi2009 ; Maier2011 ; Lanata2009 ; Mori2006 ; Nishiguchi2013 ; Voo2015 ; Ruger . What distinguishes the latter models from the single-layer approach is, among others, the appearance of the interlayer pairing channels which can contribute to the total superconducting gap. In this respect, a possible choice of pairing symmetries are discussed in Refs. Medhi2007, ; Medhi2007_2, . The interalyer pairing leads to deviations from the pure $d$-$wave$ character of the superconducting order parameter what may be important in the context of some of the ARPES experiments indicating that the gap symmetry is more complex than a simple $d$-$wave$ Zhao2007 ; Ding1995 ; Vobornik1999 .

Another effect analyzed within the bilayer models is the interlayer pair tunneling. It has been shown that the initially proposed interlayer Josephson coupling arising as a second-order process in the interlayer electron hoppingChakravarty1993 is insufficient to enhance significantly superconductivity due to a small value of the interlayer hopping $t_{z}$ with respect to the intralayer value $t$ (such coupling would be $\propto t_{z}^{2}/t$). The pair tunneling term has also been introduced in a phenomenological manner to postulate the form of the in-phase gap function between two layersChen2012 . A different approach is to investigate the interlayer pair hopping as originating from the matrix elements of the long-range Coulombic interactionNishiguchi2013 .

Here we analyze the influence of interalyer processes such as the interlayer electron hopping, the interlayer pair hopping, and the interlayer exchange coupling, on the bulk superconducting state within the $t$-$J$-$U$ model of the bilayer structure. One pf the purposes of such analysis is to see to what extent the universal properties of the single plane systems are preserved. To take into account the strong electronic correlations in the system, we use the diagrammatic expansion Gutzwiller wave function method (DE-GWF)Bunemann2012 ; Gebhard1990 ; Kaczmarczyk2013 ; Kaczmarczyk2014 . The choice of the model has been dictated by a good quantitative agreement between the theoretical results and principal experimental data for the copper-based materials, which has been reported recently for the single layer version of the same modelSpalek2016 . The gap symmetry selected in our analysis is a mixture of $d_{x^{2}-y^{2}}$ symmetry within the layer and the so-called $s^{\pm}$ symmetry for the interlayer contributionMaier2011 ; Ubbens1994 ; Lee1995 . The resulting gap function has the form $\Delta(\mathbf{k})=\Delta_{d}(\mathbf{k})\pm\Delta_{\perp}$ with the $+$ ($-$) sign, corresponding to the bonding (antibonding) band and $\Delta_{d}(\mathbf{k})$ representing the intralayer $d$-$wave$ gap. As reported in Refs. Lee1995, ; Ubbens1994, , such a state has been found as stable for the case of the $t$-$J$ model within the slave-boson analysis. In contradistinction to these results, the coexistence of the interlayer and intralayer pairings has not been found within of the variational Monte Carlo (VMC) method for the same modelMaier2011 . Within our analysis of the bilayer $t$-$J$-$U$ model, the admixture of the interlayer $s^{\pm}$ contribution leads to a small lowering of the system energy with the interlayer gap magnitude being one order of magnitude smaller than the intralayer one. As shown earlierKaczmarczyk2014 for the case of a single-layer $t$-$J$ model, our results are of the VMC quality.

For the sake of completeness, we present here also the results for the bilayer Hubbard model and analyze the influence of the differentiation between the values of the Hubbard $U$ corresponding to the two layers.

The structure of the paper is as follows. In Section II and III we show explicitly the model Hamiltonian and describe the basic concept behind the DE-GWF method as applied to the bilayer structure. In Section IV we provide the results of our calculations and analyze them. In particular, we compare our results with available experimental data for the dispersion relation, the Fermi velocity, and the bilayer splitting, all in the nodal direction. The summary is provided in the last Section. Note also that the results supplement our previous analysis for a single-layer situationSpalek2016 that should be regarded as Part I of the series. In that paper the methodology behind selecting the $t$-$J$-$U$ model is also discussed.

II Model

We consider the $t$-$J$-$U$ model for the bilayer high-TC superconductor represented by the following Hamiltonian

[TABLE]

where

[TABLE]

represents the intralayer part and

[TABLE]

describes a number of possible interlayer dynamical processes. The primed (double primed) summation means that $i\neq j$ (and $l\neq l^{\prime}$), where $i,j$ are indices of sites within the layer and $l,l^{\prime}$ enumerate the layers. The first two terms of the intralayer Hamiltonian (2) represent the Hubbard model which consist of the hopping and the intrasite Coulomb repulsion terms. The third term represents the antiferromagnetic exchange interaction. We take into account the nearest and the next-nearest neighbor hopping terms with $t=-0.35$ eV and $t^{\prime}=0.25|t|$, respectively (unless stated otherwise). Moreover, the $J_{ij}$ integrals are assumed nonzero only for the nearest neighbors, for which $J_{ij}\equiv J\equiv 0.25|t|$. The intraatomic Coulomb repulsion is set to $U=14|t|$, unless stated otherwise. The model can thus be considered as either the extended Hubbard or $t$-$J$ Hamiltonian with $J$ originating from $p$-$d$ superexchange and a moderately high value of $U$ in these charge-transfer Mott-Hubbard systemsSpalek2016 .

The first two terms of the interlayer Hamiltonian $\hat{\mathcal{H}}_{\perp}$ represent the interlayer hopping and exchange, respectively. The hopping parameters $t^{\perp}_{ij}$ are chosen so that in reciprocal space the single-particle part takes the form

[TABLE]

where $t_{\textrm{bs}}$ introduces the $\mathbf{k}$-independent bilayer splitting, whereas the $\mathbf{k}$-dependent part is taken in the usually form Andersen1995 ; Andersen1994 and is responsible for a highly anisotropic character of the splitting. In our analysis, we assume that $J^{\perp}_{ij}$ is nonzero for $i=j$ only. The remaining two terms of (3) represent the so-called on-site and off-site interlayer pair hoppings Nishiguchi2013 , respectively. They originate from the long-range Coulomb interaction. Follwing Ref. Nishiguchi2013, we set $U^{\prime}=-2U^{\prime\prime}$. In our analysis we take into account only the nearest-neighbor interlayer pair hoppings (for $|\mathbf{R}_{i}-\mathbf{R}_{j}|=1$). Nevertheless, the interlayer pair hopping generates four-site terms (two sites on the first layer and two on the second) and it is very difficult to take them into account within the DE-GWF method in higher orders. Hence, their influence is going to be taken into account within the zeroth-order calculation scheme, which is equivalent to the modified Gutzwiller approximation (SGA)Jedrak2011 ; Kaczmarczyk2011 ; Zegrodnik2014 . Such an approach is justified by the fact the $U^{\prime}$ and $U^{\prime\prime}$ interaction parameters do not represent the dominant energies in the system (i.e., they are much smaller than $U$).

We have checked that the interlayer density-density Coulomb interaction term, $\sideset{}{{}^{\prime\prime}}{\sum}_{ill^{\prime}}V_{ll^{\prime}}\hat{n}_{il}\hat{n}_{il^{\prime}}$, has a very small negative influence on the paired state. That is why we do not include this term in the analysis presented below. The influence of the intralayer intersite Coulomb repulsion is also omitted here; its influence on the SC state has been discussed by us separately within the single-layer model in Ref. Abram2016, .

A methodological remark is in place here. Namely, in principle the last term in (3) should also have a corresponding term in (2) $\sim J\sideset{}{{}^{\prime}}{\sum}_{ijl}\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{il\downarrow}\hat{c}_{jl\downarrow}\hat{c}_{jl\uparrow}$. However, it is of the order $Jd^{2}$, where $d^{2}=\langle\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}\rangle$ is the double occupancy probability, and under the circumstance that $J\ll U$, as well as $J$ is substantially smaller than $|U^{\prime}|$ it can be safely neglected. Note also that the present formulation is based on real-space pairing, whereas the original formulation was based on a direct Josephson tunneling of the Cooper pairs in $\mathbf{k}$-spaceChakravarty1993 ; Anderson1994 ; Chen2012_2 . The coherent pair tunneling across the planes has also been confirmedHettinger1995 .

III Method

To take into account the strong electronic correlations we start with the Gutzwiller-type wave function in the form

[TABLE]

where $|\Psi_{0}\rangle$ represents an uncorrelated state (chosen as the superconducting uncorrelated state) and

[TABLE]

with $\lambda_{il}$ being the variational parameters and $|\Gamma\rangle_{il}$ representing the states from the local basis

[TABLE]

The consecutive states represent empty, singly, and doubly occupied local configurations, respectively. To significantly simplify the calculations and improve the convergence, one customarily imposes the conditionBunemann2012 ; Gebhard1990

[TABLE]

where $x_{il}$ is a variational parameter and $\hat{d}^{\textrm{HF}}_{il}=\hat{n}_{il\uparrow}^{\textrm{HF}}\hat{n}_{il\downarrow}^{\textrm{HF}}$, $\hat{n}_{il\sigma}^{\textrm{HF}}=\hat{n}_{il\sigma}-n_{l0}$, with $n_{l0}=\langle\Psi_{0}|\hat{n}_{il\sigma}|\Psi_{0}\rangle$. By using Eqs. (6) and (8), one can express the parameters $\lambda_{il}$ with the use of $x_{il}$. Furthermore, by assuming spatial homogeneity within the layer we are left with the two variational parameters $x_{l}$ ($l=1,2$). For the case with two equivalent layers, which is analyzed in detail here, we can additionally set $x_{1}=x_{2}$. However, at the end of the paper we also consider the case with different interaction parameters for the two layers. The approach allows to express the expectation values, in the correlated state $|\Psi_{G}\rangle$, corresponding to each of the terms in the starting Hamiltonian (1). For example, for the case of the hopping term the corresponding expectation value takes the form

[TABLE]

where $\hat{d}^{\textrm{HF}}_{\varnothing}\equiv 1$, $\tilde{c}^{(\dagger)}_{il\sigma}\equiv\hat{P}_{il}\hat{c}^{(\dagger)}_{il\sigma}\hat{P}_{il}$, and the index $m$ corresponds to lattice sites within the layer, whereas $f$ labels the layers. The primmed summation on the right hand side is restricted to $(l_{p},m_{p})\neq(l_{p^{\prime}},m_{p^{\prime}})$, $(l_{p},m_{p})\neq(i,l)$, $(l_{p},m_{p})\neq(j,l^{\prime})$ for all $p$, $p^{\prime}$. The powers indices $k_{1}$ ($k_{2}$) express how many times the indices $f_{p}$ on the right hand of the equation have the value $1$ ($2$) for a given term. They fulfill the relation $k_{1}+k_{2}=k$. The maximal $k$, 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 express all the averages in the noncorrelated state ($\langle...\rangle_{0}\equiv\langle\Psi_{0}|...|\Psi_{0}\rangle$), which appear in the formulas for the ground state energy of the system in the form of diagrams, with vertices attached to the lattice sites and the edges corresponding to the so-called paramagnetic lines, $P_{ijll^{\prime}\sigma}\equiv\langle\hat{c}^{\dagger}_{il\sigma}\hat{c}_{jl^{\prime}\sigma}\rangle_{0}$, or the superconducting lines, $S_{ijll^{\prime}}\equiv\langle\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl^{\prime}\downarrow}\rangle_{0}$. In our calculations we take into account lines for which $|\Delta\mathbf{R}|=|\mathbf{R}_{il}-\mathbf{R}_{jl^{\prime}}|\leq 3$.

When it comes to the SC lines we assume the $d$-$wave$ spin-singlet pairing with no intrasite pairing components, which corresponds to the strong correlation limit for the copper-based superconductors. The approach has been discussed extensively before Wysokinski2015 ; Wysokinski2016 ; Kaczmarczyk2013 ; Kaczmarczyk2014 ; Spalek2016 . For the interlayer pairing contribution we take the $s^{\pm}$ pairing symmetryMaier2011 ; Ubbens1994 ; Lee1995 .

The method described briefly above allows to express the expectation value of our Hamilotnian in the correlated state

[TABLE]

as a function of all the lines $P_{ijll^{\prime}\sigma}$, $S_{ijll^{\prime}}$ and the variational parameters $x_{1}$, $x_{2}$.

To calculate the values of the paramagnetic and superconducting lines, as well as the variational parameters, the grand-canonial potential $\mathcal{F}=\langle\hat{H}\rangle_{G}-\mu_{G}n_{G}$ is minimized, where $\mu_{G}$ and $n_{G}$ are respectively the chemical potential and the number of particles per lattice site determined in the correlated state. The minimization condition together with the normalization of $|\Psi_{0}\rangle$, can be cast into the form of an effective Schrödinger equation

[TABLE]

with the following effective Hamiltonian

[TABLE]

where the effective hopping and the effective superconducting gap parameters are defined through the relations

[TABLE]

Such an approach has been commonly used in the literature Kaczmarczyk2013 ; Kaczmarczyk2014 ; Bunemann2012 ; Wang2006 and as shown in Ref. Kaczmarczyk2014phmag, , it is equivalent to the method based on the Lagrange multipliers which in turn guarantees the preservation of the statistical consistency during the minimization procedure.

As one can see from (11) the uncorrelated state $|\Psi_{0}\rangle$ is an eigenstate of $\hat{H}_{\textrm{eff}}$, which has a BCS-like form with the kinetic energy term and the pairing term leading to nonzero expectation values $\langle\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl^{\prime}\downarrow}\rangle_{0}$. Nevertheless, in contrast to the BCS theory, here, the pairing is defined in real space and it is not due to the phonon-mediated interaction. Instead, it stems from the inter-electronic correlations taken into account in higher orders via diagrammatic expansion. The terms for $l^{\prime}\neq l$ in the second summation of (12) corresponds to the interlayer pairing.

To determine $|\Psi_{0}\rangle$ we first transform the effective Hamiltonian to reciprocal space and then diagnonalize it by using the Nambu representation. This allows for the derivation of the self-consistent equations for the superconducting and paramagnetic lines, as well as the chemical potential. Additionally, the parameters $\{x_{l}\}$ are determined variationally. After the determination of $P_{ijll^{\prime}\sigma}$, $S_{ijll^{\prime}}$, $\mu$, $x_{1}$, $x_{2}$, we can calculate the correlated equivalents of the superconducting lines, $\langle\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl^{\prime}\downarrow}\rangle_{G}$, which are regarded as the true gap parameters in the correlated state. In general, each non-correlated gap parameter $\langle\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl^{\prime}\downarrow}\rangle_{0}$ taken into account within the diagrammatic procedure has its correlated correspondent $\langle\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl^{\prime}\downarrow}\rangle_{G}$. The symmetry of the correlated gap reflects the symmetry assumed for the non-correlated gap (intralayer $d$-$wave$ and interlayer $s^{\pm}$ wave). As it has been pointed out in Refs. Kaczmarczyk2013, , Kaczmarczyk2014, , Abram2016, , the dominant contribution to the superconducting phase comes from $\langle\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl\downarrow}\rangle_{G}$ with $\mathbf{R}_{i}-\mathbf{R}_{j}=(a,0)$ for the $d$-$wave$ intralayer pairing. For the $s^{\pm}$ interlayer pairing the dominant contribution comes from $\langle\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{jl^{\prime}\downarrow}\rangle_{G}$ for $\mathbf{R}_{i}-\mathbf{R}_{j}=(0,0)$ and $l\neq l^{\prime}$. Those two intra- and inter-layer correlated gap parameters are going to be denoted by $\Delta^{||}_{G}$ and $\Delta^{\perp}_{G}$ throughout the paper and lead to the following $\mathbf{k}$-dependence of the overall gap

[TABLE]

where the $\pm$ sign corresponds to bonding and antibonding hybridized bands, respectively.

IV Results and discussion

Our analysis is performed in the fourth order of the diagrammatic expansion, if not stated otherwise. As shown in Refs. Kaczmarczyk2014, and Abram2016, , for the considered types of models, satisfactory convergence is achieved in the 4-5 order of the expansion. We focus mainly on the $t$-$J$-$U$ bilayer case analysis (Subsection A). As it has been already stated we set the intralayer model parameters to $t=-0.35$ eV, $t^{\prime}=0.25|t|$, $U=14|t|$, $J=0.25|t|$, unless stated otherwise. In the following we provide the values of the model parameters in the units of $|t|$. For the sake of completeness we also present selected results for the case of the Hubbard model (Subsection B).

IV.1 Results for the $t$-$J$-$U$ model

First, we analyze the influence of the interlayer electron hopping terms on the paired state, with no interlayer exchange interaction included ($J^{\perp}=0$). The magnitude of the $\mathbf{k}^{||}$-dependent part of the interlayer hopping is governed by $t_{z}$ (cf. Eq. (4)) and is responsible for the anisotropy of the bilayer splitting. This term leads to a significant splitting in the antinodal direction ($\pi$,[math]) and no splitting in the nodal one ($\pi$,$\pi$). In some papers referring to the bilayer cuprate superconductors, the nodal bilayer splitting is neglected, i.e., the situation with $t_{bs}=0$ is considered Mori2006 ; Nishiguchi2013 ; Chen2012 . At the same time, considerations with the nearest-neighbor interlayer hopping included only (with $t_{z}=0$) have been also carried out Maier2011 ; Lanata2009 ; Medhi2009 ; Voo2015 . In the latter situation the splitting between the bonding and the antibonding bands is isotropic which is rather loosely connected to the electronic structure of the copper-based compounds. In Fig. 1 we show the influence of both the isotropic and anisotropic bilayer splittings on the paired state. As one can see, the effect of increasing either $t_{\textrm{bs}}$ or $t_{z}$ is similar. Due to the hopping of electrons between the layers the intralayer correlated gap diminishes mainly in the overdoped regime, whereas the interlayer one increases. This is caused by the fact that if there is no interlayer interaction ($J^{\perp}=0$), the hopping mainly dilutes the intralayer pair bonds and thus weakens the pairing. Nevertheless, in the underdoped regime the gap parameter $\Delta^{||}_{G}$ remains almost unchanged by the interlayer hopping. Also, the interlayer gap is approximately one order of magnitude smaller than the intralayer one. As one would expect, the case with $t_{z}=t_{\textrm{bs}}\equiv 0$ leads to intralayer pairing only ($\Delta^{\perp}_{G}=0$). 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 hopping.

Next, we scrutinize the influence of the interlayer exchange coupling $J^{\perp}$ on the SC gap. The results are presented in Fig. 2. It can be seen that the $J^{\perp}$ term has a very small influence on the intralayer pairing while its interlayer correspondent is significantly increased (almost twice) when changing the value of $J^{\perp}$ from 0.0 to 0.25. It should be noted that even for the case when $J^{\perp}=J$ (i.e., when the intra- and inter-layer exchange couplings are equal), the intralayer gap parameter is still about one order of magnitude larger that the interlayer one. This effect should be related to the fact that the system is infinite in the ($x$,$y$) plane whereas it is limited in the $z$-direction (two layers only).

In Fig. 3 we also show the influence of the interlayer pair hopping processes on the superconducting gap. Following Nishiguchi et al. Nishiguchi2013 , we set $U^{\prime}=-2U^{\prime\prime}$. As one could expect, the inclusion of the pair tunneling process enhances the intralayer SC gap and broadens the SC stability regime. However, due to the mentioned effect the maximum value of the interlayer gap is decreased. The positive effect on the interlayer correlated gap is significant even though the values of the $U"$ parameter have been taken as relatively small in comparison to the Hubbard $U$. Note that the interlayer pair hopping term within this analysis originates from the long-range Coulomb interactionNishiguchi2013 .

From Figs. 1-3 one can see that the most important factor in changing effectively the gap is the Coulomb interactions ($U^{\prime}$, $U^{\prime\prime}$) leading to the interlayer pair hopping processes. The role of the number of layers in the elementary cell in increasing the stability of the SC phase was analyzed earlier within the mean-field approachByczuk1995 .

In the remaining part of this subsection we show the results for the case when the model parameters are set to reproduce the selected principal characteristics measured by ARPES experiments of Ref. Kordyuk2004, . According to the experimental data obtained for Bi-2212, the bilayer splitting of the Fermi surface is highly anisotropic, with maximal value in the antinodal direction and small, but nonzero value in the nodal direction. Moreover, even in the overdoped regime both the bonding and the antibonding parts of the Fermi surface appear to be hole-like. To reproduce the data within our model both $t_{z}$ and $t_{\textrm{bs}}$ have to be taken as nonzero ($t_{z}=-0.1$ and $t_{\textrm{bs}}=-0.06$). In these considerations we have neglected the appearance of the interlayer pairing mechanism. Such a choice is justified by the fact that the interlayer gap magnitude is one order of magnitude smaller than its intralayer corespondant (cf. Fig. 1). In Fig. 4 (a) and (b) we plot the calculated two sheets of Fermi surface and the dispersion relations close to the corresponding Fermi momenta for selected values of doping close to the optimal doping and in the underdoped regime, respectively. The theoretical composed Fermi surface topology shown in Fig. 4 (a) is roughly similar to that obtained experimentally in Ref. Kordyuk2004, (cf. Fig. 1. in that paper) and the theoretical lines which reflect the dispersion relations in Fig. 4 (b) are relatively close to the experimental points. However, it should be noted that the authors of Ref. Kordyuk2004, provide only the critical temperature of the corresponding samples which allows to determine the corresponding doping regions only (underdoped or overdoped region). The specific doping levels of the samples are not provided explicitly and this is one of the reasons why our theoretical results may not reflect precisely those coming from the experiment. In Fig. 4 (c) we plot the bilayer splitting in energy (top) and in the reciprocal space (bottom) in the nodal direction, both as a function of doping. The gray region represents the measured values for different doping levels presented in Ref. Kordyuk2004, . The experimentally measured very weak doping dependence of the nodal bilayer splitting value has been confirmed in Ref. Borisenko2006, for YBCO. However, for this compound, the value of $\Delta_{\mathrm{BS}}$ was almost five times larger than the one measured for Bi-2212. As one can see, our theoretical results are close to the experimental ones only in the overdoped regime while in the underdoped regime the calculated bilayer splitting significantly decreases. It should be noted that Fournier et al. Furnier2010 have reported a complete vanishing of the bilayer splitting below the optimal doping for the YBCO samples. The inconsistence between the different sets of experimental data from the one side, and the theoretical results from the other requires a further analysis.

In Fig. 5 we show the effect of the interlayer coupling on the bilayer splitting and the nondal Fermi velocity. As one can see, with the increasing $J^{\perp}$ the bilayer splitting also increases. However, the Fermi velocity is very weakly dependent on the $J^{\perp}$ value and even for rather strong interlayer couplings the changes in $v_{F}$ are barely visible. The same situation takes place for the correlated gap which is also almost unaffected by this factor (cf. Fig. 2).

Due to the nonzero splitting in the nodal direction, also two values of the nodal Fermi velocity appear, which are shown in Fig. 4 (d). The measured values reported in Ref. Kordyuk2004, are close to $2.0$ eVÅ  and they correspond to the bonding band. Also, other experimental reports corresponding to the cuprate superconductors, demonstrate very weak doping dependence of the nodal Fermi velocity close to that value Zhou2003 ; Kordyuk2005 ; Borisenko2006 . In this case, the weak doping dependence is reproduced by the theory as can be seen in Fig. 4 (d). Moreover, the bonding and the antibonding Fermi velocities are very close. Near half filling both of them have almost the same value. It is caused by the fact that the nodal bilayer splitting is decreasing as we are getting closer to the half filled situation, what is illustrated also in Fig. 4 (e). For the sake of completeness, in Fig. 4 (f) we also provide the correlated gap magnitude as a function of doping for the same set of model parameters as the remaining results presented in this Figure.

IV.2 Results for the Hubbard model

For the sake of completeness, we have also analyzed the selected results for the case of the Hubbard model ($J\equiv 0$, $J^{\perp}\equiv 0$). For these calculations we set the interlayer hopping as non-dependent on $\mathbf{k}$ ($t_{z}=0$). In the main part of this Subsection we also neglect the interlayer pair hopping term ($U^{\prime}=U^{\prime\prime}=0$). Nonetheless, the influence of nonzero $U^{\prime}$ and $U^{\prime\prime}$ is shown briefly at the end.

As one can see in Fig. 6, the effect of the splitting parameter $t_{\textrm{bs}}$ is similar as in the case of the $t$-$J$-$U$ model. However, for large values of $t_{\textrm{bs}}$, a discontinuity appears in the $\delta$-dependence of both $\Delta^{||}_{G}$ and $\Delta^{\perp}_{G}$ in the underdoped regime (the curves for $t_{\textrm{bs}}=0.75$), which has not been observed earlier.

We have also analyzed the influence of the disproportionation between the Coulomb repulsion in the two layers. In the following analysis $U_{1}$ and $U_{2}$ correspond to the Coulomb repulsion in the first and the second layer, respectively. In such a situation the intralayer correlated gaps, the average number of electrons per atomic site, and the average double occupancies in the two layers can have different values. All the mentioned quantities are presented in Figs. 7 and 8. It can be seen that the correlated gaps for both layers decrease as the difference between $U_{1}$ and $U_{2}$ increases. For $U_{2}=0$ both superconducting gaps vanish (cf. Fig. 8 (a)). In Fig. 8 (b) we show that in the considered system electrons are pushed from the layer with stronger Coulomb repulsion to that with the weaker. Obviously, the double occupancy in the second layer increases with decreasing $U_{2}$, as can be seen in Fig. 8 (c). For the case with $U_{1}=14$ and $U_{2}=0$ the paired phase vanishes, $\Delta_{G}^{(1)}=\Delta_{G}^{(2)}=0$ (cf. Fig. 8(a)). However, after setting to interlayer pair hopping as nonzero, we can retain the stability of the paired phase in the system (cf. Fig. 8 (d)). In this situation, the pairing originating from the correlation effects in the layer with nonzero Coulomb repulsion is transfered to the second layer (with no Coulomb repulsion) by the pair hopping process. Such situation brings into mind the systems composed of one superconducting monolayer of correlated compound placed on a metallic layer in which the paired phase can be induced by the proximity effect Yuli2008 . However, in such interpretation the Cooper pair tunneling terms included in the Hamiltonian would have to originate from the Josephson coupling.

V Summary

We start with a general remark. The majority of the papers concerning the real space pairing in high-Tc superconductors deal with a single Cu-O plane as reflecting the intrinsic properties of the copper-based compounds. In our recent (preceding) paper Spalek2016 we have shown that indeed in such an approach one can rationalize, even in a quantitative manner, the principal ground-state properties of high-TC superconductors. Nevertheless, the basic question still remains as to whether the model two-dimensionality, even if realistic, provides similar results when the basic aspects of three-dimensionality are incorporated into the theoretical scheme. Our present work provides an affirmative answer to this question, even though the detailed features are model-parameter dependent. We believe that the direct relation between theory and experiment analyzed here and in our previous workSpalek2016 , can contribute significantly to the resolution of the fundamental question whether a purely electronic model based on electronic correlations (without a glueAnderson2007 ) can rationalize the superconducting properties of the copper systems.

Turning to our results here, we have analyzed the influence of different interlayer dynamical processes on the superconducting state within the bilayer $t$-$J$-$U$ model. When it comes to the bilayer splitting, induced by the interlayer electron hopping, both the anisotropic ($\mathbf{k}^{||}$-dependent term with the magnitude tuned by parameter $t_{z}$) and the isotropic (of magnitude $t_{\textrm{bs}}$) contributions have been analyzed. Our analysis shows that both terms have a similar influence on the supercodncuting gap parameters (cf. Fig. 1). With the increasing interlayer hopping magnitudes the intralayer SC gap decreases and that of the interlayer part increases. For the case of the bilayer Hubbard model, for large value of $t_{\textrm{bs}}=-0.75$ the discontinuity appears in the doping dependence of the SC gap (cf. Fig. 6). As this effect has not been observed experimentally, the $t_{\mathrm{bs}}$ value must be essentially smaller than the value of $|t|$.

It has been reported recently in Ref. Voo2015, that the results of the VMC calculations for the bilayer $t$-$J$-$U$ model for $d$-$wave$ superconducting phase are insensitive to $t_{\perp}$ in the range $|t_{\perp}|<0.5$ for $\delta<0.2$. As one can see from Fig. 1 also within our analysis both interlayer hopping terms do not influence significantly the correlated gap amplitude, $\Delta_{G}^{||}$, in the underdoped regime for not too large values of $|t_{\mathrm{bs}}|$ and $|t_{z}|$.

The main influence of the interlayer exchange coupling is seen in the interlayer superconducting gap magnitude, while the intralayer gap is only slightly reduced with increasing $J^{\perp}$ (cf. Fig. 2). Weak dependence of the intralayer pairing on $J^{\perp}$ is also mentioned in Ref. Voo2015, . The nonzero exchange integral $J^{\perp}$ explains the nature of antiferromagnetism in the insulating stateColdea2001 ; Chakravarty1989 . Therefore, it is tempting to conclude that for nonzero and small values of $J^{\perp}$ one can rationalize the nature of AF state and at the same time uphold the validity of the principal results concerning the superconducting phase obtained for the single-plane systems.

Among the considered interlayer processes the most significant influence on the SC phase has been reported for the interlayer pair hopping. The corresponding terms have also been included within the analysis of the bilayer Hubbard model in Ref. Nishiguchi2013, . Similarly as in our case, it is necessary to include at least one of the off-site terms to achieve a meaningful enhancement of the SC phase stability. In our situation, this is a direct consequence of the fact that we assume $\langle\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{il\downarrow}\rangle\equiv 0$, as a result of the strong on-site Coulomb repulsion. Within the analysis shown in Ref. Nishiguchi2013, 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 not to be straightforward. However, their analysis is carried out within the spin-fluctuation picture for the case of weak interactions, what means that their approach is very different from ours.

We have also analyzed the basic features which refer to the available data obtained by ARPES experiments (such as the Fermi surface two-sheet structure, the bilayer splitting, the dispersion relations, the nodal Fermi velocity). In this part of our analysis, the model parameters have been fitted and roughly reproduce the experimental data presented in Ref. Kordyuk2004, . As one can see, the highly anisotropic character of the bilayer splitting, with a small value of the splitting in the nodal direction, has been reproduced (cf. Fig. 4 (a)). The calculated value of the bilayer splitting is similar to that measured in the experiment for overdoped samples (cf. Fig. 4 (b)). However, in the underdoped regime significant discrepancies between our theory and experiment appear (cf. Fig. 4 (c) and (e)). The mentioned reportKordyuk2004 does not provide specific values of the doping levels, so the precise comparison cannot by carried out. At this point, it should be noted that in disagreement with Ref. Kordyuk2004, , Fournier et al. Furnier2010 have reported a vanishing bilayer splitting below the optimal doping for the YBCO samples. Such a behavior is not reproduced within our approach. Finally, we also show that the nodal Fermi velocities, corresponding to the bonding and the antibonding bands, are very similar and close to the value measured experimentally ($\approx 2.0$ eVÅ ) in the whole doping range (cf. Fig. 4 (d)). The present analysis for the bilayer $t$-$J$-$U$ model complements our recent resultsSpalek2016 for a single Cu-O planar superconductor.

The results presented for the case of the bilayer Hubbard model with the two different values of the Coulomb repulsion for the two layers show that the disproportionation between $U_{1}$ and $U_{2}$ leads to different carrier concentrations in the layers, what in turn has a negative influence on the paired phase (cf. Fig. 8 a and b). However, for the situation representing a metallic monolayer in contact with a strongly correlated superconducting monolayer ($U_{1}=14$ and $U_{2}=0$), the paired phase can be induced by the Cooper pair tunneling processes (cf. Fig. 8 d).

The bilayer Hubbard model considered here might also be of interest in reference to the ultra-cold atom systems analyzed in recent yearsNishida2010 ; Klawunn2010 .

VI Acknowledgement

We would like to acknowledge the financial support from the Grant MAESTRO, No. DEC-2012/04/A/ST3/00342 from the National Science Centre (NCN) of Poland.