Superconducting mechanism for the cuprate Ba 2 CuO 3 + δ based on a multiorbital Lieb lattice model

For the recently discovered cuprate superconductor Ba 2 CuO 3 + δ , we propose a lattice structure which resembles the model considered by Lieb to represent the vastly oxygen-deﬁcient material. We ﬁrst investigate the stability of the Lieb-lattice structure and then construct a multiorbital Hubbard model based on ﬁrst-principles calculation. By applying the ﬂuctuation-exchange approximation to the model and solving the linearized Eliashberg equation, we show that s -wave and d -wave pairings closely compete with each other and, more interestingly, that the intraorbital and interorbital pairings coexist. We further show that if the energy of the d 3 z 2 − r 2 band is raised to make it “incipient” with the lower edge of the band close to the Fermi level within a realistic band ﬁlling regime, s ± -wave superconductivity is strongly enhanced. We reveal an intriguing relation between the Lieb model and the two-orbital model for the usual K 2 NiF 4 structure where a close competition between s - and d -wave pairings is known to occur. The enhanced superconductivity in the present model is further shown to be related to an enhancement found previously in the bilayer Hubbard model with an incipient band.


I. INTRODUCTION
More than 30 years have passed since the discovery of the high-T c cuprates, but a full understanding of their physics remains one of the most challenging problems in the condensed matter physics [1].However, one strong consensus has been reached: The CuO 2 planes play an essential role in the occurrence of superconductivity.Namely, the cuprates have a layered perovskite crystal structure, where a copper atom is surrounded by oxygens, typically with an octahedral coordination.Since the octahedron is elongated in the c-axis direction, the crystal field splitting makes the 3d x 2 −y 2 orbital have the highest energy among the 3d orbitals.Hence, the d 9 electron configuration results in a situation where the electronic structure can be regarded as basically a single-band system.Indeed, some of the present authors have shown that there is a systematic material dependence, in which T c is basically increased as the one-band character (3d x 2 −y 2 ) becomes stronger, i.e., when the energy of the 3d 3z 2 −r 2 orbital is lowered below that of 3d x 2 −y 2 , which is realized for higher apical oxygen heights [2][3][4][5].
The recent experimental discovery by Li et al. [6] of another type of cuprate superconductor, Ba 2 CuO 3+δ , is remarkable in this context.The material, having a layered structure, exhibits T c = 73 K, which is much higher than that of "214" La 2−x Sr x CuO 4 [7] with T c 40 K, but more interestingly, a dramatic feature, among others not seen in conventional cuprates, is that a large amount of oxygen deficiencies exist in the CuO 2 planes [8].Details of the sample preparation is reported to be as follows: Ba 2 CuO 3+δ samples are synthesized in a tetragonal symmetry at a much higher pressure (18 GPa) than usual and at a temperature of 1000 • C in a polycrystalline form.This is in contrast with a lower pressure synthesis in which only an orthorhombic phase is synthesized.This implies that the tetragonal phase, even if metastable, is stabilized with the high-pressure synthesis.The excess oxygens O δ are also added in the Cu-O planes by the extremely high pressure synthesis with δ 0.2.This immediately raises a puzzle regarding the origin of the high T c because the CuO 2 planes should simply be disrupted at this level of O deficiency from the conventional Cu-O plane.Another notable feature in Ba 2 CuO 3+δ is that the combination of the oxygen content of 3 + δ 3.2 and the +2 valence of Ba should make the electron configuration significantly deviate from d 9 ; namely, an unprecedented amount of holes (as large as ≈ 40%) exist.This sharply contrasts with the conventional wisdom for the cuprates that superconductivity is optimized around 15% hole doping [1].Yet another curious feature is that each CuO octahedron is compressed rather than elongated along the c axis with the apical oxygen height smaller than the inplane Cu-O distance, so that the Cu 3d 3z 2 −r 2 orbital should  [45] software.Cu site 4 is not displayed for clarity.be higher approaching that of 3d x 2 −y 2 , and so a multiband, multiorbital situation is expected.These features are all in strong contradiction with the high-T c condition for the conventional cuprates, which suggests that an alternative pairing mechanism may be at work in this new material.Indeed, a number of theoretical studies have proposed various pairing mechanisms based on various lattice structures [9][10][11][12][13][14].The experimental finding of Ba 2 CuO 3+δ may also shed a light on the previous finding of Sr 2 CuO 3+δ [15][16][17], which also possesses a large amount of oxygen deficiencies and a T c as high as 90 K but a much lower superconducting fraction than in the Ba 2 CuO 3+δ .
Given this background, a theoretical challenge is that how we can construct a model and fathom the structure of the gap function for the material, which has hugely oxygendepleted CuO 2 planes.Assuming that the deficiencies are ordered, some candidates for the crystal structure have been proposed.Liu et al. propose a chain-type structure, which actually exists in Sr 2 CuO 3 [10,[18][19][20].Li et al. predict a ladder-type lattice based on an automated structure inversion method [12].Le et al. propose a structure where a matrix of Ba 2 CuO 4 with CuO 2 planes is embedded in the background of Ba 2 CuO 3 [11].More recently, another type of lattice dubbed as the brick-wall model has been proposed [14].
Thus, the lattice structures considered so far (other than the conventional K 2 NiF 4 type) have one-dimensional natures in some sense or other, but an experiment [6] suggests the material has tetragonal symmetry.This has motivated us to propose here another structure, depicted in Figs.1(a) and 1(c), as a candidate for the undoped Ba 2 CuO 3 ("213" composition), where by doping we mean adding excess oxygens.We call the proposed structure the "Lieb-lattice type," since it resembles the model considered by Lieb [21] if we focus on the Cu sites 1, 2, and 3 in Fig. 1 and ignore Cu site 4, which is shown to be electronically irrelevant.The model considered by Lieb possesses a flat band in the band structure, and, in the context of magnetism, it is theoretically proven that ferromagnetism occurs at half-filling when the on-site repulsive interaction U is turned on.A superconducting mechanism exploiting the flat band of the Lieb lattice has also been proposed [22].Lieb originally considered a class of models with different numbers of sublattice sites, and superconductivity in such a model in a quasi-1D structure has also been studied with the density-matrix renormalization group [23].Here, however, we shall see that the model derived in the present study is actually distinct from the original (single-orbital) Lieb model, since the present material inherently has a multiorbital nature, as we shall show.
We start with an investigation of the stability of the Lieb lattice in terms of the total energy and phonon calculations for the lattice structure of Ba 2 CuO 3 and calculate its electronic band structure.The obtained band structure is then used to construct multiorbital models, for which we apply the fluctuation exchange (FLEX) approximation [24,25] to study the superconductivity.We show that s-wave and d-wave pairings closely compete with each other, where we find a peculiar case of coexisting intraorbital and interorbital pairings.We further show that superconductivity is strongly enhanced if we increase the energy of the d 3z 2 −r 2 band (from its original position obtained by first-principles calculation for Ba 2 CuO 3 ) to make it "incipient" [23,[26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44], where the lower band edge comes close to the Fermi level within a realistic band filling regime.In an even wider scope, we reveal that the Lieb model has an intimate relation with the two-orbital model of the K 2 NiF 4 structure where a close competition between s-wave and d-wave pairings is known to occur [9].We finally point out a relation between the enhanced superconductivity in the present models and an enhancement found previously in the bilayer Hubbard model with an incipient band.

II. FORMULATION
We consider the Lieb-lattice-type structure for Ba 2 CuO 3 , where the in-plane oxygen deficiencies are ordered as shown in Fig. 1(a).Details of the Ba atom position and the unit cell of this structure are depicted in Appendix A. We perform structural optimization using the Vienna Ab Initio Simulation Package (VASP) [46,47].Here, we adopt the generalized gradient approximation formulated by Perdew, Burke, and Ernzerhof for the exchange-correlation energy functional [48] and the projector augmented wave method [49] without the inclusion of the spin-orbit coupling, and take an 8 × 8 × 8 k-mesh with a plane-wave cutoff energy of E cut = 650 eV.We also examine the dynamical stability of the Lieb lattice by performing phonon calculation.We employ the finite displacement method as implemented in the PHONOPY software [50] in combination with VASP.We took a 2 × 2 × 2 supercell and a 3 × 3 × 3 k-mesh.Other conditions such as the energy cutoff are the same as those adopted in the structural optimization, which is always the case for phonon calculations in general.
For the optimized lattice structure, we obtain the electronic band structure taking a 6 × 6 × 6 k-mesh with a planewave cutoff energy of E cut = 550 eV.From the electronic band structure, we extract the maximally localized Wannier functions [51,52] using the WANNIER90 code [53].Here we disregard very small hopping parameters to simplify the multiorbital Lieb lattice models (see Appendix B for details).
In order to take account of the electron correlation effects beyond the first principles band calculation, we further introduce the on-site multiorbital interactions with a Hamiltonian, Here, i denotes the sites, μ, ν indicate the orbitals, σ represents the spins, c † iμσ creates an electron, and n iμσ = c † iμσ c iμσ .Interactions are U , the intraorbital repulsion; U , the interorbital repulsion; J, Hund's coupling; and J , the pair hopping.We do not consider electron-phonon interactions, since our aim is to investigate an electronic mechanism of superconductivity.To analyze the many-body effect, here we adopt the FLEX approximation.
In the FLEX approximation, renormalized Green's function is determined self-consistently from the Dyson equation, where the self-energy is calculated by taking the bubble and ladder diagrams that consist of the irreducible susceptibility, which is calculated from the renormalized Green's function G.Here q = (q, ω) stands for the wave vector q and the Matsubara frequency ω, l i denotes the orbitals, T is the temperature, and N is the number of k points.In order to avoid double counting of the effect of the electron interaction already considered in the first principles calculation, we subtract the ω = 0 component of the self-energy Re (k, 0) during the self-consistent loop following Ref.[54].Hence, the Fermi surface of the models remains unchanged even after the correlation effects are taken into account by the FLEX calculation.It should be noted that the double counting is not rigorously avoided since (k, ω = 0) in FLEX is not the same as that in the DFT calculation.
Using the obtained Green's function along with the spin ( χs ) and charge ( χc ) susceptibilities, which are matrices for multiorbital systems with the interaction matrices given as we obtain the effective spin-singlet pairing interaction, which is plugged into the linearized Eliashberg equation, where the gap function μν (k) is also a matrix.G in Eq. ( 8) is the renormalized Green's function obtained from the FLEX calculation.Through G, the mass renormalization and finitelifetime effects are taken into account.The maximum eigenvalue λ of this equation reaches unity at T = T c , so that λ calculated at a fixed temperature can be a measure of T c .Throughout the present study, we calculate λ at T = 0.01 eV.We refer to the eigenfunction of the linearized Eliashberg equation as the gap function.Note that since Eq. ( 8) is a linearized equation, the absolute value of the gap function does not have any physical meaning, and only its relative magnitude among different orbital components and its symmetry are relevant.Both the Green's functions and the gap functions are obtained first in the orbital representation, which can be transformed into the band representation with a unitary transformation.Green's functions will be presented by taking its absolute value.Also, Green's functions and the gap functions will be presented for the lowest Fermionic Matsubara frequency iπ k B T , and the effective pairing interactions αββ α will be presented at the lowest bosonic Matsubara frequency 0.
Assuming a rigid band obtained for the Lieb-lattice-type model, we vary the band filling in a regime that contains a case corresponding to Ba 2 CuO 3+δ with a realistic δ ≈ 0.2.In the calculation, we take 2048 Matsubara frequencies and a 16 × 16 × 2 k-point mesh.We have checked that calculation taking 32 × 32 × 2 k-point mesh gives essentially the same results.

A. Stability of the Lieb structure for Ba 2 CuO 3
We start with the stability of the Lieb-lattice-type structure for Ba 2 CuO 3 .The obtained total energy of the optimized lattice structure is E tot (Lieb) = −33.95eV.For comparison, we have also performed structural optimization for the chain-type structure of Ba 2 CuO 3 , shown in Fig. 1(b) (see Appendix A (0,0,π) (0,0,0) (0,0,0) (π,π,0) (π,0,0) for the actual structure), whose total energy is estimated as E tot (chain) = −33.99eV.Thus, the total energies of the two structures turn out to be quite close to each other; considering the accuracy of the first principles calculation, a difference of 40 meV can be reversed, e.g., by the effects of the correlation and/or excess oxygens not taken into account here [55].Then, given the fact that the chain-type structure is known to be realized in an existing material Sr 2 CuO 3 [18][19][20] but inconsistent with the tetragonal symmetry, the Lieb-latticetype structure may be considered as a realistic candidate for the lattice structure of Ba 2 CuO 3+δ .We further calculate the phonon dispersion as presented in Fig. 2. We find that no imaginary phonon modes are present for this lattice structure, which suggests dynamical stability of the present Lieb-lattice structure.
As for the apical oxygen position determined by structural optimization, its distance measured from the in-plane Cu site turns out to depend on the site: 2.18 Å above Cu site 1, 1.96 Å above Cu sites 2 and 3, and 1.82 Å above Cu site 4. Namely, the sites with smaller oxygen coordination numbers have lower apical oxygen heights.The average value is 1.98 Å, which is substantially smaller than the value (2.42 Å) for La 2 CuO 4 .We may note that this is qualitatively consistent with the experimental value of 1.86 Å [6], if we consider the fact that the excess oxygens in the actual material will increase in-plane holes that should attract the apical oxygens.In fact, for a cuprate La 2−x Sr x CuO 4 , it is actually observed experimentally in Ref. [56] that hole doping by substituting La with Sr (which has nearly the same ion radius as La) induces reduction of the apical oxygen height from 2.42 Å (x = 0) to 2.30 Å (x = 0.2).

B. Electronic band structure of the Lieb lattice
From the optimized crystal structure, we obtain the electronic band structure as presented in Fig. 3, where we also display the weight of the Cu d x 2 −y 2 and Cu d 3z 2 −r 2 orbital characters.The bands originating from the orbitals at Cu site 4, which has no neighboring oxygens in the plane, have energies very low, so that site 4 is irrelevant in the model building.We then extract tight-binding models downfolded in terms of the maximally localized Wannier orbitals.We consider the two e g orbitals centered at each of the Cu sites 1, 2, and 3, which results in a six-orbital model, where the ].An important point here is that the low apical oxygen height makes the upper (lower) bonding band mainly originated from the d 3z 2 −r 2 (d x 2 −y 2 ) at site 1, because the energy levels of these orbitals are inverted from those in the conventional cuprates, as suggested in previous studies [6,9,11].The top two bands are the antibonding bands from the hybridization between site 1 and sites 2 and 3. Here, the orbitals at sites 2 and 3, hybridized with site 1 orbitals, have energy higher than site 1 d x 2 −y 2 and d 3z 2 −r 2 orbitals because the apical oxygens at sites 2 and 3 are closer to Cu.
In the actual material Ba 2 CuO 3+δ , the Fermi level should be shifted downward to intersect the middle two bands because the oxygen content is larger than in Ba 2 CuO 3 .Therefore, we further construct a two-orbital model that extracts the two middle bands, as shown in Fig. 4(b), by considering d x 2 −y 2 and d 3z 2 −r 2 Wannier functions centered at Cu site 1.In this model, the orbitals of Cu sites 2 and 3, as well as the oxygen 2p orbitals, are implicitly taken into account through the Wannier functions.In this two-orbital model, we shall sometimes refer to the d x 2 −y 2 and d 3z 2 −r 2 Wannier orbitals as orbitals 1 and 2, respectively.Note that although we cannot rule out the possibility of a high-spin state, in the following analysis we assume that the ground state is in a low-spin state.
We vary the band filling in the following FLEX calculation assuming a rigid band in these models, where the correspondence between the oxygen content and the band filling is as The parameter values adopted are U = 2.0 eV, J = J = U/10, and U = U − 2J.In the two-orbital model, we use U = 2 eV, which is somewhat smaller than that of the conventional cuprates [62][63][64], to take into account the wide spread of the Wannier functions across neighboring Cu atoms.
follows.In Ba 2 CuO 3 , the nominal Cu valence is +2, so the electron configuration is 3d 9 on average, that is, (three e g electrons) × (four sites) = (12 e g electrons) per unit cell.Since the e g orbitals at Cu site 4 are fully occupied by electrons, there are 8 e g electrons per unit cell in the six-orbital model.Namely, the band filling n, defined as the number of electrons per spin per unit cell, is n = 4. Similarly, n = 3 corresponds to Ba 2 CuO 3.25 .The band filling of the six-orbital model subtracted by two gives that of the two-orbital model because the bottom two occupied bands are ignored in the latter; namely, n = 1 and n = 2 in the two-orbital model correspond to Ba 2 CuO 3.25 and Ba 2 CuO 3 , respectively.An interesting point here is that the oxygen content of O 3.25 , which is close to the actual experimental situation of ≈O 3.2 [6] and implies a large amount of hole doping (≈50%) in the usual sense of the term, in fact corresponds to half-filling in the present two-orbital model of the Lieb-type lattice.We shall indeed see that the electron correlation effect are maximized and thus superconductivity is optimized around this band filling.
For later reference, we have also obtained, via structural optimization, the electronic band structure for the K 2 NiF 4structured Ba 2 CuO 4 by VASP taking a 8 × 8 × 8 k-mesh with a plane-wave cutoff energy of E cut = 550 eV [Figs.3(d) and 3(e)].There, we construct a two-orbital model by extracting the e g orbitals centered at the Cu sites [Fig.4(d)].We shall discuss the relation between the Lieb lattice-type and K 2 NiF 4 -type structures in the Discussions section.

C. Superconductivity
We now move on to the FLEX calculation for superconductivity.We start with the two-orbital model for the Lieb lattice.Figure 5 plots the eigenvalue of the Eliashberg equation λ at T = 0.01 eV for the s-and d-wave pairing symmetries.We can see that the two pairing symmetries give somewhat close values of λ, where the s wave slightly dominates within the parameter regime studied.In both symmetries, λ is maximized at n = 1.This is because electron correlation effects are maximized around half-filling, and, as stressed above, we can notice this band filling corresponds to the oxygen content of O 3.25 , close to the actual material.
Before we go any further, we have to carefully examine the gap functions and the definition of the pairing symmetry.Figure 6 displays the gap functions in both orbital and band representations.Note above all that the gap function in the present multiorbital system is a matrix.A curious finding in Fig. 6 is that the pairing symmetry is inverted between diagonal and off-diagonal matrix elements, i.e., s-wave diagonal elements are accompanied by d-wave off-diagonal ones [Figs.6(a (c) and 6(d)] representations.We can trace the curious phenomenon back to the hybridization between the d x 2 −y 2 and d 3z 2 −r 2 orbitals, where the hybridization has sign structure in a d x 2 −y 2 symmetry.In Fig. 5, we have abbreviated the diagonal s-wave with off-diagonal d wave as s wave, and the diagonal d with off-diagonal s as d wave.We will adopt this abbreviation hereafter.
If we focus on the orbital-diagonal elements of the gap function, the d wave has the cos(k x ) − cos(k y ) form usually encountered in similar analysis of the conventional cuprates.In real space, this corresponds to a nearest-neighbor pairing, whose wave function changes its sign upon 90-deg rotation.On the other hand, the s-wave gap roughly has the form cos(k x ) + cos(k y ), which implies that this is basically an extended s-wave with a pair residing on nearest-neighbor sites.We can compare the band representation of the s-wave gap function with Green's function |G(k)| in Figs.6(e) and 6(f).The ridges in |G(k)| represent the Fermi surface, and we have two pieces for the Fermi surface in this two-orbital model as seen from |G 11 | and |G 22 |.We can then realize that the gap has a sign inversion across the two Fermi surfaces; namely, we have here the so-called s±-wave gap function, as depicted in the left panel of Fig. 7(b), which is reminiscent of the ironbased superconductors [28,57] as far as the Fermi surfaces are concerned.
If we now turn to the orbital-off-diagonal elements of the gap function, the symmetries (s or d) are exchanged from the diagonal elements as we have noted, but we can also notice in Fig. 6 that the amplitudes of the off-diagonal elements are comparable to those of the diagonal ones, which indicates that the interorbital pairing has significant contributions.Now, interorbital spin-triplet pairing has been studied in the multiorbital Hubbard model that has degenerate orbitals [58][59][60][61], but in the present case the gap function has the symmetry αβ (k) = βα (−k) that signifies a singlet pairing.We can intuitively grasp, in real space, the coexistence of the interand intraorbital pairings in Fig. 7(a), where both of intraorbital and interorbital pairs reside on nearest neighbors, as we have explained.Since pairing between electrons having large energy difference is unlikely, this kind of pairing is peculiar to systems where the two orbitals are close in energy.
We now grasp how the coexisting inter-and intraorbital pairings arise.We find, in the typical parameter regime considered, that the diagonal components of Green's function (at the lowest Matsubara frequency) are larger than the offdiagonal ones, with their real part much larger than the imaginary part, and that they satisfy Re where m, l = orbital 1 (d x 2 −y 2 ) or 2 (d 3z 2 −r 2 ).We can thus roughly extract the contributions of these components in the linearized Eliashberg equation as for the intraorbital [Eqs.(9) and (12)] and interorbital [Eqs.(10) and (11)] pair scattering channels with the Feynman diagrams for the pairing interaction αββ α as depicted in Figs.8(i) and 8(j) and Figs.9(i) and 9(j).See Table I.
To identify which interaction parameters in the Hamiltonian dominate these pairing interactions, we calculate the eigenvalue λ against the band filling for various choices of (U, U , J, J ) as shown in Fig. 10.Here we permit breaking the orbital rotational symmetry (U = U − 2J) in order to extract the effect of each interaction.The pairing interaction αββ α for the varied interactions is depicted in Figs. 8 and 9.We find in Fig. 10 that increasing the interorbital interactions U , J, and J , enhances λ.On the other hand, Fig. 10(c A salient feature here is that αββ α s all have peaks around q = (π, π ).From Eqs. ( 9)-( 12 (j).Note that the pairing interactions are plotted over the range 0 q x , q y 2π to display the peak structure around (π, π ) clearly.
that the value of d-wave λ (≈0.1, not shown) is smaller than that of s-wave λ for the original value of t 12 .Also, in the absence of t 12 , the off-diagonal component, 12 (j).Note that the pairing interactions are plotted over the range 0 q x , q y 2π to display the peak structure around (π, π ) clearly.arrive at a mechanism in which the effect of the interorbital hybridization can enhance the superconductivity through the coexisting intra-and interorbital pairings.
Finally, let us turn to the six-orbital model in its FLEX results.We show the band-filling dependence of λ in Fig. 11, and the gap functions and Green's functions in Fig. 12, for s-and d-wave pairings.We find that the results are similar to those obtained for the two-orbital model.Namely, s wave and d wave closely compete with each other, and the intra-and interorbital pairing components coexist.Further understanding of the relation between the two-and six-orbital models is given in the Discussions section below.

D. Dependence on the level offset between the two orbitals
So far we have adopted the tight-binding parameter values (the hoppings and the on-site energies) estimated by firstprinciples band calculation and Wannierization.The obtained maximum value of λ is about 0.3 at T = 0.01 eV, which is not large enough to explain the T c ≈ 73 K experimentally found in Ba 2 CuO 3+δ [6].If we stick to these parameter values, some additional pairing mechanisms (e.g., electron-phonon interaction) which boost the T c would be required.We note, however, that there are some ambiguities in the parameter values in the Hamiltonian, particularly the energy level offset between d x 2 −y 2 and d 3z 2 −r 2 orbitals.First, in the actual samples used in the experiment with excess oxygens, the holes doped into the planes may attract the apical oxygen more strongly than theoretically estimated here, as we have mentioned above.Indeed, the theoretical average apical height, 1.98 Å, estimated here for Ba 2 CuO 3 , is larger than the experimental value (1.86 Å).The excess oxygens themselves would repel the apical oxygens, but this effect should be insignificant at Cu site 1, where no additional oxygens can be coordinated.The lowered apical oxygens at site 1 would further push up the d 3z 2 −r 2 energy level.Second, the level offset may be affected by correlation effects that are not taken into account in the first-principles calculation.For instance, Ref. [65] studied a nickelate superlattice system where the d x 2 −y 2 and d 3z 2 −r 2 levels are inverted, and found that the correlation effect taken into account with the dynamical mean-field theory pushes the d 3z 2 −r 2 band just above the Fermi level to make its Fermi surface vanishing.
With these considerations, let us probe how the eigenvalue λ against the band filling changes when the level offset is varied, in order to seek a possibility for a further enhancement of superconductivity.In Fig. 13, the eigenvalue of the Eliashberg equation of the two-orbital model is plotted against the band filling for various values of E , where the d 3z 2 −r 2 level is changed.The result does reveal intriguing features.If we first focus on the region around quarter-filling n = 0.5 (one electron per two orbitals, corresponding to an oxygen content of 3 + δ = 3.375), the λ is enhanced as E increases.Here we find that the dominating pairing symmetry changes from s wave to d wave.This is because the system approaches a half-filled single-band system for higher d 3z 2 −r 2 level, favoring d x 2 −y 2 -wave pairing [66].This resembles the situation in the conventional cuprates with nearly d 9 electron configuration (three electrons per two e g orbitals), where a sufficiently low d 3z 2 −r 2 level results in an effective single-band system comprising the d x 2 −y 2 orbital and hence favors d-wave We note that the FLEX calculation did not converge for some cases around n = 1 due to large spin fluctuations, and hence the data points for those cases are missing.
pairing [2][3][4][5]: The difference is that the d 3z 2 −r 2 is moved away from the main band upward in the present case or downward in the latter.However, the eigenvalue λ here is not so high as in the typical single-band cuprates such as HgBa 2 CuO 4 because the band width is narrower.
If we now turn to the region around n = 1, E smaller than the original value somewhat enhances the s-wave λ.We find that a smaller E increases 1111 and 1122 (not shown), and we speculate that the near degeneracy of the two orbitals favors the interorbital pairing through the enhancement of 1122 (the intraorbital scattering of interorbital pairs).The enhancement of λ upon reducing E further confirms our statement in Sec.III C that the interorbital hybridization, which induces interorbital pairing, is favorable for superconductivity.
What is even more interesting and realistic is the case of larger E in the n = 1 regime (which does correspond to 3 + δ = 3.25, close to the experimental situation [6]).There, λ corresponding to the s wave in Fig. 13 is strongly enhanced when E is increased to some extent from its original value.In this parameter regime, we find that the s wave strongly dominates over the d wave.The maximum value of the s-wave λ is as large as 0.6, which is close to the value obtained for HgBa 2 CuO 4 , a superconductor with T c 100 K [67].
In fact, for E where λ is optimized, the bottom of the d 3z 2 −r 2 band is close to the Fermi level (as indicated in inset of Fig. 13).Recently, such a band lying just above or below the Fermi level is referred to as an "incipient band" and has received attention, especially in the context of the iron-based superconductors [27][28][29][30][31][32][33][34][35]38], where hole bands lying just below the Fermi level are observed in some materials [30,31,[68][69][70][71][72][73][74][75].In a wider context, the possibility of the occurrence or strong enhancement of superconductivity due to an incipient band has long been proposed for multiband Hubbard models on various types of lattices [23,26,36,37,[39][40][41][42][43][44][76][77][78][79][80][81][82].A salient feature in these cases is that the gap function typically exhibits nodeless "s±-wave" symmetry.Indeed, the gap function obtained for the present two-orbital model exactly has a nodeless s±-wave symmetry, as displayed in Fig. 14.Another prominent feature in the gap function is that the interorbital pairing (off-diagonal element in the orbital representation) is suppressed compared to the case with the original value of E (Fig. 6).This may seem to contradict with what we have concluded previously, namely, that the interorbital pairing is favorable for the superconductivity.We shall further discuss this issue in Sec.IV B. Also, if we look in Fig. 15 at the effective pairing interactions for the incipient band case, we can find that 1221 , which describes the interorbital scattering of intraorbital pairs, is strongly enhanced.We shall come back to the relation between the present model and those in the previous studies also in Sec.IV.

A. Relation with the K 2 NiF 4 -type structure
In the present Lieb-lattice model, s± and d x 2 −y 2 waves are found to compete with each other.A close competition between s± and d x 2 −y 2 has been found in a previous theoretical study, where the K 2 NiF 4 -type structure with a reduced apical oxygen height was adopted [9,83].It is an intriguing problem how these are possibly related.
In Fig. 4(d), we have shown the first-principles band structure of Ba 2 CuO 4 in the K 2 NiF 4 -type structure.We can notice that the band structure of the two-orbital model for the K 2 NiF 4 -type structure is similar to that of Ba 2 CuO 3 in the Lieb-lattice structure [Fig.4(b)], except for the band width.This is in fact understandable because, in the Lieb lattice, site-2 and site-3 orbitals extend toward site 1 [Fig.4(c), lower right panel], so that they can be regarded as playing a role of the oxygen 2p σ orbitals in the K 2 NiF 4 structure; in this structure, the oxygen 2p orbitals have an energy somewhat lower than the Cu 3d orbitals, while in the Lieb lattice, the site-2 and site-3 orbitals have somewhat higher energies than the site-1 orbitals, presumably because the apical oxygen heights at sites 2 and 3 are lower than at site 1.The bandwidth of the latter is narrower than the former because the electron hoppings between Cu site-1 orbital and the Cu site-2 and site-3 orbitals are smaller than those between Cu 3d and O 2p orbitals.We can therefore state that, starting from the conventional CuO 2 plane and removing the oxygens to form a Lieb lattice, we unexpectedly encounter an analog of the CuO 2 plane, on a smaller energy scale.Since the band structures are similar, so are the FLEX results.In Figs.16 and 17, we show the FLEX result for the two-orbital model of the K 2 NiF 4 -type structure.We again end up with a close competition between s±-wave and d x 2 −y 2 -wave pairings, which is qualitatively consistent with the previous random-phase approximation study [9].We note that the eigenvalue λ for the K 2 NiF 4 -type structure in Fig. 16 is larger than that for the Lieb-lattice-type structure in Fig. 5 because of the wider bandwidth, which might seem more consistent with the experiment [6] from the viewpoint of the high T c , but this, of course, is not the case because the K 2 NiF 4 -type structure does not take into account the large amount of oxygen vacancies observed experimentally.
From the above consideration, we can make an interesting observation that the relation between the two-orbital and sixorbital models for the Lieb lattice is analogous to the relation between the single-band Hubbard model and the three-band d-p model in the conventional cuprates.

B. Relation with the bilayer model
We have found that s±-wave superconductivity is strongly enhanced in the Lieb-lattice model when the d 3z 2 −r 2 band is raised so that it becomes incipient.Such a strong enhancement of s±-wave superconductivity reminds us of the bilayer Hubbard model on a square lattice [39,42,43,[76][77][78][79][80][84][85][86][87][88][89][90][91][92], where enhanced superconductivity is found when one of the bands is incipient [39,42,43,[76][77][78][79][80].However, the bilayer  Hubbard model is a single-orbital (one orbital per site) system with two sites per unit cell.The enhancement of superconductivity in such a system is mediated by spin fluctuations originating from the on-site U , which gives rise to the pairing interaction between the bonding and antibonding bands, both of which have equal weight of the two orbitals in a unit cell.Similar enhancement mechanism of superconductivity in multiorbital systems has been discussed in the context of the iron-based superconductors, where portions of the electron and hole Fermi surfaces having the same orbital character interact via the effective interaction ( llll in the present notation) enhanced mainly by the intraorbital U [93].In the present Lieb-lattice model and also in the two-orbital model of the K 2 NiF 4 -type structure, the situation is distinct in that the orbital character is quite different between the two bands.The lower band has a strong d x 2 −y 2 character, while the upper band is dominated by d 3z 2 −r 2 character, so that here the interorbital interactions (U , J, J ) should be the key.The present view that interorbital interactions play an important role can be reinforced from the pairing interactions presented in Fig. 15, where the interorbital pair scattering vertices ( 1221 and 1212 ) are large, and also from the gap functions in Fig. 14, where the sign of the gap function is reversed between d x 2 −y 2 and d 3z 2 −r 2 orbitals.Figure 14 also shows that the orbital and band representations resemble each other, which implies the two bands have different orbital characters.
Can we identify the reason why superconductivity is enhanced even when the incipient band has an orbital character different from that of the main band?Let us propose a succinct way to understand this.As shown in Appendix C, the single-orbital (with one orbital per site) Hubbard model on a bilayer square lattice [Fig.18(a)] with an on-site interaction U can be transformed into a two-orbital Hubbard model on a (monolayer) square lattice with all the on-site intra-and interorbital interactions being U/2 [94].In this transformation, the bonding and antibonding orbitals in the bilayer system, comprising the two sites connected by the vertical interlayer hopping t ⊥ , translate to the d x 2 −y 2 and d 3z 2 −r 2 orbitals in the monolayer system, with E playing a role of 2t ⊥ .In the bilayer Hubbard model, whose noninteracting band structure is depicted in Fig. 18(b), superconductivity is found to be strongly enhanced when one of the bands becomes (nearly) incipient upon increasing t ⊥ [39,42,43,[76][77][78][79].
From the viewpoint of the transformation introduced here, the result for the two-orbital model of the Lieb lattice thus corresponds to that of the bilayer model in that superconductivity is strongly enhanced when the upper band becomes incipient upon increasing E .The gap function of the Lieb two-orbital model also resembles that obtained for the bilayer model.We show in Fig. 18(c  the same unit cell (connected by t ⊥ ), while in the two-orbital model a mixing of intra-and interunit cell pairings takes place.
The interorbital repulsion U is usually known to enhance charge or orbital fluctuations [95], which generally compete with spin fluctuations in mediating Cooper pairing because the charge and orbital (spin) fluctuations give rise to attractive (repulsive) pairing interactions.An intriguing point to note in the present case is that U plays a crucial role in enhancing spin fluctuations, as can be captured from the analogy with the bilayer Hubbard model where spin fluctuations solely dominate.
The two-orbital to bilayer transformation appears to be approximately valid only when E is not too small; namely, s wave dominates over d wave when E is small in the present two-orbital model, while d wave is dominant in the bilayer model with small t ⊥ [77].This is presumably because the cos(k x ) − cos(k y ) form of the hybridization cannot be transformed into the bilayer square lattice form of the Hamiltonian (with the off-diagonal elements in Eq. (C5) in Appendix C having tetragonal symmetry), so that the transformation loses its validity as E becomes smaller than the interorbital hopping between d x 2 −y 2 and d 3z 2 −r 2 orbitals.When E is large, on the other hand, the interorbital hybridization loses its significance, so that the transformation becomes more valid.However, when E is small, the effect of the hybridization will be prominent, so that d wave gives way to s wave, in contrast to the bilayer model case, because, as mentioned in Sec.III C, the interorbital hopping makes the s-wave pairing more favorable.
We mentioned in Secs.III C and III D that the interorbital hybridization is favorable for superconductivity.However, we find that this is not the case when E is large as in the incipient band situation; there, if we turn off the interorbital hopping t 12 , the interorbital component of the gap function 12 vanishes, but λ is enhanced; namely, the better correspondence to the bilayer model is more favorable for superconductivity.Since the interorbital and the incipient-band enhanced pairing mechanisms are essentially different, whether the presence of interorbital pairing is favorable for superconductivity or not depends on the magnitude of the level offset between the two orbitals.
We note that apart from the problem of Ba 2 CuO 3+δ , the occurrence of s±-wave superconductivity in the cuprates was also proposed for a model [96] of highly overdoped CuO 2 monolayer grown on Bi 2 Sr 2 CaCu 2 O 8+δ [97].In this model, the d 3z 2 −r 2 band lies below the d x 2 −y 2 band as in the conventional cuprates, but due to the large amount of holes, the Fermi level not only intersects the d x 2 −y 2 band but also intersects the top of d 3z 2 −r 2 band.This resembles the incipient-band situation of the present two-orbital model, if we make an electron-hole transformation.

V. CONCLUSIONS
In the present study, we have proposed the Lieb lattice as a candidate for the lattice structure of the newly discovered superconductor Ba 2 CuO 3+δ .We have shown from the total energy that the proposed lattice structure is almost as stable as the chain-type structure that is known to exist.The dynamical stability of the proposed structure is also shown through a phonon calculation.
Applying a FLEX approximation to the relevant twoorbital and six-orbital models derived from the first-principles band calculation, we find that coexistence of intra-and interorbital pairings arises due to the relatively small energy level offset between the d x 2 −y 2 and d 3z 2 −r 2 orbitals.As for the pairing symmetry, s-wave and d-wave pairings closely compete with each other, with the former dominating.Superconductivity is optimized around the band filling corresponding to the oxygen content of 3 + δ = 3.25, which is close to that of the actual material.
The maximum eigenvalue of the Eliashberg equation is not large enough to explain the observed T c .While a cooperation with other pairing such as phonons may be necessary to fully understand the experiment, we have proposed an alternative scenario for explaining the observed T c by varying the level offset between the two orbitals, which is motivated from the consideration that the level offset may be larger in the actual material Ba 2 CuO 3+δ than its first-principles estimation for Ba 2 CuO 3 .We have indeed found that s±wave superconductivity is strongly enhanced when the d 3z 2 −r 2 band is raised in energy so that it becomes nearly incipient around the band filling corresponding to the oxygen content of 3 + δ = 3.25.In this situation, in contrast to the case with smaller level offset, the amplitude of the interorbital pairing is small, while the interorbital pair scattering plays an essential role.
We have then noted that both the band structure and the FLEX results resemble those of the two-orbital model for the K 2 NiF 4 -type structure.We have traced its origin back to the fact that the Cu orbitals at sites 2 and 3 in the Lieb lattice play the role of the oxygen 2p σ orbitals in the K 2 NiF 4 structure, so that the electronic structure of the former is analogous to that of the latter.
From this observation, we have further pointed out a relation between the two-orbital model for the Lieb lattice and the Hubbard model on the bilayer square lattice.When one of the bands is incipient, the two models exhibit similar results regarding the enhancement of superconductivity and the nodeless form of the gap function.The resemblance suggests that the transformation between the two-orbital model and the bilayer model, which is shown to be rigorous when the intraand interorbital interactions are equal, is valid to some extent even when the interactions are not equal.
In the present study, we have focused on the 2-1-3 composition and varied the band filling assuming a rigid band.It will be an interesting and important future problem to explicitly investigate the effect of the "+δ" excess oxygens.

APPENDIX A: THE CRYSTAL STRUCTURES OF THE LIEB-LATTICE-TYPE AND THE CHAIN-TYPE Ba 2 CuO 3
In the main text, we have proposed a Lieb-lattice-type structure for Ba 2 CuO 3 and discussed the stability for comparison with the chain-type structure.For their stacking geometry, we show the detail in Fig. 19 as the side views of these structures that we actually use in the first-principles calculation.For the Lieb lattice, layers are stacked in such a way that the in-plane components of the translation vector between adjacent layers are always ≈( 14 , 1 4 ) in units of the lattice constants, as depicted by a side view [Fig.19(a)] and also by a top view (Fig. 20 stacking pattern breaks tetragonal symmetry as can be seen from Fig. 20; in fact, one way to strictly preserve this symmetry is to have four layers in a unit cell, where the in-plane components of the translation vector between neighboring layers are ≈( 14 , 1 4 ), ( 1 4 , − 1 4 ), (− 1 4 , − 1 4 ), and (− 1 4 , 1 4 ), but that would result in a very large unit cell.We adopt the structure depicted in Figs.19 and 20 to reduce the size of the unit cell and hence the calculation cost.In practice, however, we find that the structure we adopt approximately preserves tetragonal symmetry in that the tight-binding parameters of the obtained models possess tetragonal symmetry within the accuracy ≈10 −4 eV.As explained in Appendix B below, these small parameters are disregarded in the FLEX  In order to simplify the multiorbital models, we disregard small hopping parameters |t| < 1.0 × 10 −2 eV in the sixorbital model or |t| < 5.0 × 10 −2 eV in the two-orbital model.

APPENDIX C: THE RELATION BETWEEN THE TWO-ORBITAL MODEL AND THE BILAYER MODEL
Let us explain here the relation between the single-orbital (one orbital per site) bilayer model and the two-orbital model [94].We label the two sites in a unit cell of the bilayer lattice as i = 1, 2, and the orbitals in the two-orbital model as a, b = α, β.The on-site interaction part of the Hamiltonian of the bilayer Hubbard model is Namely, we end up with a two-orbital model where the on-site intra-and interorbital interactions all have the same strength, U/2.We can also show how the kinetic energy part of the Hamiltonian is transformed.Let c kiσ , c † kiσ be the Fourier transform of c miσ , c † miσ .The kinetic energy part is then given in momentum space as For instance, for the bilayer model on a square lattice with only the in-plane nearest-neighbor hopping t and the vertical interplane hopping t ⊥ [Fig.18 We can thus see that the term t ⊥ contained in ε corresponds to the energy offset E between the two orbitals.

2 FIG. 1 .
FIG. 1. Cu-O plane in (a) the Lieb-lattice-type structure and (b) the chain-type structure.The apical oxygen positions (not displayed) are all occupied in both cases.Ba atoms reside at the same sites as in the K 2 NiF 4 -type structure.(c) A bird's-eye view of the Lieb structure with VESTA [45] software.Cu site 4 is not displayed for clarity.

FIG. 5 .
FIG.5.Eigenvalue λ of the two-orbital model for s-wave (intraorbital s and interorbital d) and d-wave (intraorbital d and interorbital s) pairings plotted against the band filling.The parameter values adopted are U = 2.0 eV, J = J = U/10, and U = U − 2J.In the two-orbital model, we use U = 2 eV, which is somewhat smaller than that of the conventional cuprates[62][63][64], to take into account the wide spread of the Wannier functions across neighboring Cu atoms.
) and 6(c)], while d-wave diagonal elements are accompanied by s-wave off-diagonal ones [Figs.6(b) and 6(d)].This occurs both in orbital [Figs.6(a) and 6(b)] and band [Figs.6 FIG. 6.For the two-orbital model, (a) the intraorbital s wave with interorbital d-wave gap functions.The portions of the gap functions that change sign across the wave vector (π, π ) are indicated by yellow arrows (see the text).(b) The intraorbital d wave with interorbital s-wave gap functions are displayed in the orbital representation.Panels (c) and (d) represent them in the band representation, respectively.Note that since Eq.(8) is a linearized equation, the absolute value of the gap function does not have any physical meaning, and only its relative magnitude among different orbital components and its symmetry are relevant.(e) The absolute value of Green's function in the orbital representation, while (f) shows them in the band representation.The parameter values adopted are n = 1.0,U = 2.0 eV, J = J = U/10, and U = U − 2J.

FIG. 7 .
FIG. 7. (a) Schematics of the coexisting intra-and interorbital nearest-neighbor pairings.In the upper figure, we depict the energy level of the two orbitals centered at site 1 of the Lieb-type lattice and the electrons occupying those orbitals and forming nearestneighbor pairs.Note that "nearest neighbor" here refers to the nearest-neighbor unit cells rather than the sites.(b) The Fermi surface at n = 1.0 for the original band structure of the two-orbital model (left) and that of the band structure in the absence of the interorbital hopping (right).The signs in the gap functions are indicated for the s± wave (left) and d wave (right).
) shows increasing the intraorbital U initially enhances λ, which is rounded off for larger U .If we compare this with Figs.9(a)-(c) and 9(e)-(g), we reveal that increasing U , J, and J enhances 1221 and 1212 , which in turn enhances λ.On the other hand, increasing U enhances 1111 [Figs.8(a) and 8(d)] and 1122 [Figs.8(e) and 8(h)], but suppresses 1221 [Figs.9(a) and 9(d)] and 1212 [Figs.9(e) and 9(h)].The increase of 1111 enhances the intraorbital pairing while the increase of 1122 enhances the interorbital pairings, but the suppression of 1221 degrades intraorbital pairings while the suppression of 1212 degrades interorbital pairings, which is probably the origin of the nonmonotonic behavior against the U variation.

FIG. 8 .
FIG. 8.The interaction dependence of the effective intraorbital interactions, 1111 [(a)-(d)], and interorbital interactions, 1122 [(e)-(h)], of intraorbital pairs at the lowest Matsubara frequency.The bottom panels depict the Feynman diagram of 1111 (i) and 1122(j).Note that the pairing interactions are plotted over the range 0 q x , q y 2π to display the peak structure around (π, π ) clearly.

FIG. 9 .
FIG. 9.The interaction dependence of the effective interorbital interactions of intraorbital pairs, 1221 [(a)-(d)], and of interorbital pairs, 1212 [(e)-(h)], at the lowest Matsubara frequency.The bottom panels depict the Feynman diagram of 1221 (i) and 1212(j).Note that the pairing interactions are plotted over the range 0 q x , q y 2π to display the peak structure around (π, π ) clearly.

n
(band filling) J = J' = 0.1 eV J = J' = 0.2 eV J = J' = 0.J' = U 10. s-wave eigenvalue λ plotted against the band filling in the two-orbital model for the interaction values varied in different ways.(a) U dependence with fixed U = 2.0 eV and J = J = 0.2 eV, (b) J, J dependence with fixed U = 2.0 eV and U = 1.6 eV, and (c) U dependence with U = U − 2J and J = J = U/10.

FIG. 11 .
FIG. 11.Eigenvalue λ in the six-orbital model for s-wave and d-wave pairings plotted against the band filling.The interaction parameters are U = 2.5 eV, J = J = U/10, and U = U − 2J.In the six-orbital model, since the Wannier orbitals are localized to each Cu atom, we use U = 2.5 eV, which is close to values evaluated for the conventional cuprates [62-64].

FIG. 12 .
FIG. 12.The s-wave gap function in the six-orbital model in the orbital representation (a) or the band representation (b), along with Green's function orbital representation (c) or band representation (d).The parameter values adopted are n = 3.0, U = 2.5 eV, J = J = U/10, and U = U − 2J.

FIG. 13 .
FIG. 13.The largest eigenvalue λ of the two-orbital model plotted against the band filling for various values of the level offset, E ≡ E d 3z 2 −r 2 − E d x 2 −y 2 .The original value is E = 0.47 eV.The pairing symmetry is s wave, except for those symbols marked with a dashed square where the symmetry is d wave.The interaction parameters are U = 2.0 eV, J = J = U/10, and U = U − 2J.The inset depicts the bare band structure with E increased by +0.5 eV.The horizontal black line represents the Fermi level for n = 0.9.We note that the FLEX calculation did not converge for some cases around n = 1 due to large spin fluctuations, and hence the data points for those cases are missing.

FIG. 14 .
FIG. 14.The s-wave gap function for the two-orbital model in the incipient-band case in the orbital representation (a) and the band representation (b), along with Green's function in the orbital representation (c) and the band representation (d).The interaction parameters are U = 2.0 eV, J = J = U/10, U = U − 2J, n = 0.9, and E = E original + 0.5 eV.

4 FIG. 16 .
FIG. 16.Eigenvalue λ in the two-orbital model of the K 2 NiF 4type structure for s-wave and d-wave pairings plotted against the band filling.The interaction parameters are U = 3.0 eV, J = J = U/10, and U = U − 2J.

FIG. 17 .
FIG. 17.The s-wave and d-wave gap functions in the two-orbital model of the K 2 NiF 4 -type structure in the orbital representation (a) or band representation (b), along with Green's function in the orbital representation (c) or band representation (d).The parameter values are n = 1.0,U = 3.0 eV, J = J = U/10, and U = U − 2J.
) the gap function of the bilayer Hubbard model in the band representation.(Note that the band representation in the bilayer model corresponds to the orbital representation in the two-orbital model.)The parameter values are determined from those for the Lieb two-orbital model withE increased by +0.5 eV using the transformation given in Appendix C. The nodeless s±-wave gap function indeed resembles that of the two-orbital model.By "nodeless," we mean that the gap within each band does not change sign, not only on the Fermi surface but over the entire Brillouin zone.Also, the enhancement of the interaction 1221 for the interorbital scattering of intraorbital pairs (Fig.15) corresponds to the dominant pairing interaction in the bilayer model that induces interband scattering of intraband pairs.All these resemblances between the two-orbital model and the bilayer model suggest that the transformation is approximately valid even when U = U = J.On the other hand, if we look more closely, the gap function of the bilayer model is nearly constant within each band [Fig.18(c)],whereas that of the present model exhibits momentum dependence, which roughly has a cos(k x ) + cos(k y ) + const.form.Namely, in the bilayer model, the pairing in real space occurs basically within (c) s-wave 1=lower band, 2=upper band) [band rep.]

FIG. 18 .
FIG. 18.(a) The bilayer square lattice, (b) band structure of this model, and (c) the s-wave gap functions of the bilayer Hubbard model in the band representation.The parameter values are n = 0.9, t = −0.12,t = 0.01, t ⊥ = 0.5, and U = 3.

3 FIG. 19 .
FIG. 19.Side views of (a) the Lieb-lattice-type and (b) chaintype structures, depicted with VESTA software.In panel (a), we indicate the translation vectors a 1 , a 2 , and a 3 .Parallelograms formed by blue lines delineate the unit cells.

FIG. 20
FIG. 20.A top view of the stacked Cu-O layers of the Lieblattice-type structure.The translation vectors are also indicated as in the previous figure.
APPENDIX B: THE TIGHT-BINDING PARAMETER VALUES FOR THE LIEB-LATTICE MODELSHere we give the tight-binding parameters of the six-orbital [Fig.4(a)] and two-orbital [Fig.4(b)]Lieb-lattice models obtained by the first-principles calculation with WANNIER90.

TABLE I .
The pairing interactions.