Superconductivity studied by solving ab initio low-energy effective Hamiltonians for carrier doped CaCuO$_2$, Bi$_2$Sr$_2$CuO$_6$, Bi$_2$Sr$_2$CaCu$_2$O$_8$, and HgBa$_2$CuO$_4$

We numerically analyze superconductivity (SC) in the cuprate superconductors by using ab initio effective Hamiltonians consisting of the antibonding combination of Cu $3d_{x^2-y^2}$ and O $2p_{\sigma}$ orbitals. We perform variational Monte Carlo calculations for the four carrier doped cuprates with diverse experimental optimal SC critical temperature $T_{c}^{\rm opt}$: CaCuO$_2$ ($T_{c}^{\rm opt} \sim 110$ K), Bi$_2$Sr$_2$CuO$_6$ ($T_{c}^{\rm opt} \sim 10$-$40$ K), Bi$_2$Sr$_2$CaCu$_2$O$_8$ ($T_{c}^{\rm opt} \sim 85$-$100$ K), and HgBa$_2$CuO$_4$ ($T_{c}^{\rm opt} \sim 90$ K). Materials and hole doping concentration ($\delta$) dependencies of the SC order parameter $F_{\rm SC}$ and the competition with spin/charge order show essential and quantitative agreements with the available experiments in the following points: (1) The ground state is commonly the SC state, which is severely competing with the charge/spin stripe and antiferromagnetic states. (2) $F_{\rm SC}$ shows amplitude consistent with the superfluid density measured in the muon spin resonance and its dome structure found in $\delta$ dependence shows consistency with that of the SC gap in the tunneling and photoemission measurements. We further find insights into the universal SC mechanism: (I) $F_{\rm SC}$ increases with the ratio $U/|t_1|$, indicating that $U/|t_1|$ is the principal component controlling the SC. Here, $U$ and $t_1$ are the onsite Coulomb repulsion and the nearest neighbor hopping, respectively, in the Hamiltonians. (II) A universal scaling $T_{c}^{\rm opt}\sim 0.16 \lvert t_1 \rvert F_{\rm SC}$ holds. (III) SC is enhanced and optimized if $U$ is increased beyond the real available materials. It is further enhanced by decreasing the offsite interaction. The present findings provide useful clues for the design of new SC materials with even higher $T_{c}^{\rm opt}$.


I. INTRODUCTION
The mechanism and origin of the large superconducting (SC) gap, high superfluid density, and high critical temperatures T c observed in high-T c superconductors, such as copper oxides, remain a central challenge in condensed matter physics.In these copper oxides, the d-wave SC state is severely competing with other orders, such as spin and charge stripes or antiferromagnetic (AFM) states, and the observed T c widely ranges from above 130 K to below 10 K. Understanding and reproducing these diverse phenomena without relying on adjustable parameters is hence desirable, especially when clarifying their origin.When ab initio calculations are able to reproduce systematic materials dependence quantitatively by solely relying on their crystal structures, it provides us with valuable insight into the universal mechanism behind and into the principal components for the enhancement of SC beyond existing materials.
Many studies have suggested severe competitions of the SC with charge/spin stripe and AFM states theoretically based on simplified Hubbard-like or t-J Hamiltonians as models of the cuprate superconductors [1][2][3][4][5][6][7][8].A positive correlation between U and T c or SC tendency was also pointed out by taking U as an adjustable parameter in the Hubbard type Hamiltonians [9].However, ab initio studies without adjustable parameters are few and it is not clear whether the diversity of the materials dependence can be accounted for.There still exists a limited number of ab initio studies: The phase diagram including the SC phase was reproduced by solving the ab initio Hamiltonian for a particular case of Hg compound [10].
Ab initio Hamiltonians were derived by Nilsson et al. for several cuprate compounds, which reported an empirical observation without solving the Hamiltonians that the experimental optimal T c is generally higher for larger U/|t| in their parameters [11,12].The relation of the charge transfer energy to T c was also pointed out [13,14].Aside from the cuprates, there exist some ab initio studies on strong-coupling superconducting materials such as the iron-based superconductors [15], fullerene [16], and nickelates [17] to discuss the superconducting properties.However, to our knowledge, there exist no systematic studies on the SC properties by solving solely ab initio Hamiltonians without adjustable parameters with the help of accurate many-body solvers to reveal the origin of the diverse materials dependence, especially for the challenging cuprates.Unless reproducing the materials dependent properties quantitatively, the universal mechanism would also not be able to be identified confidently either.
In this paper we show properties of typical cuprate superconductors calculated by solving the ab initio Hamiltonians of four families of materials, namely carrier doped CaCuO 2 , HgBa 2 CuO 4 (abbreviated as Hg1201 hereafter), Bi 2 Sr 2 CuO 6 (Bi2201 hereafter), and Bi 2 Sr 2 CaCu 2 O 8 (Bi2212 hereafter) [12], by applying a state-of-the-art quantum many-body solver based on the variational Monte Carlo (VMC) algorithm [18,19], including the combination with neural network [20,21] elaborated from earlier variational algorithms [22][23][24][25].It is experimentally known that the optimum critical temperature T opt c is realized at around δ = 0.1-0.15[26,27] for doped CaCuO 2 (T opt c ∼ 110K) and Hg1201 (T opt c ∼ 90K), and at around δ = 0.15-0.25 [28][29][30] for doped Bi2212 (T opt c ∼ 85-100K) and Bi2201 (T opt c ∼ 10-40K).We elucidate that similarity and diversity among the four families, especially the amplitude of the SC order parameters and T opt c in the experiments are accounted for by using the present ab initio results, which provides us with insights into the materials dependence and the universal mechanism: We emphasize that our ab initio analyses contain essential differences from most of the Hubbard models studies.One important difference is the presence of the realistic off-site interactions.We will clarify that this crucially stabilizes the charge uniform SC state without clear stripe long-range order.
The dominance of the SC for all the four families is successfully demonstrated.In addition, the δ dependence of F SC has a dome structure with the peak at δ ∼ 0.05-0.1 consistently with the experimental indications.On the other hand, the dome peak of T c appears at larger δ > 0.1 in the experiment.This shift from the dome peak of F SC is understood from the decreasing renormalization factor with decreasing doping, which does not affect F SC but T c .
Although T opt c has a variety among these four families, we show universally that (I) the higher SC order parameter F SC at the optimal doping basically results from a larger ratio U/|t 1 |, where U is the on-site repulsive Coulomb interaction and t 1 is the nearest neighbor hopping in our ab initio parameters of single-band effective Hamiltonian for the AB orbital of strongly hybridized Cu 3d x 2 −y 2 and O 2p σ orbitals.Furthermore, we show that (II) T opt c is determined by the scaling T opt c ∼ 0.16|t 1 |F SC .The δ dependence of the local energy suggests a universal superconducting mechanism: Though the bare in-teraction is strongly repulsive, the Mottness converts it to the strong effective attraction required for the Cooper pairing.
Despite monotonic increase of F SC with U/|t 1 | within the existing materials, we further show that (III) larger U/|t 1 | beyond the ab initio values makes the peak of F SC followed by the reduction.We find that F SC can be roughly 30% more enhanced than the ab initio case when U is 20% increased beyond the ab initio value by preserving the transfer and other off-site interaction parameters.We also show that F SC is further enhanced to as much as the double of the existing material by the additional reduction of the off-site Coulomb interaction.These searches beyond the ab initio parameters for the existing materials offer a guide for future experimental materials design.This paper is organized as follows.Section II presents the methods and computational details: First, the effective Hamiltonians studied in this paper are summarized in Sec.II A. Then, we give the numerical methods in Sec.II B and define the physical quantities in Sec.II C. We present in Sec.III the results for each of the four families of compounds.Based on the obtained ab initio results, in Sec.IV, we further explore the direction to enhance and optimize the SC order parameter by controlling the effective interaction parameters beyond the ab initio values, to gain insights into the future materials design.In Sec.V, we discuss our analyses.Summary and conclusion are given in Sec.VI.

A. Ab initio effective Hamiltonian
Within this paper, we solve the ab initio single-band effective Hamiltonians for CaCuO 2 , Hg1201, Bi2201, and Bi2212 as derived in Ref. [12].It should be noted that the single band is constructed from the AB orbital of strongly hybridized Cu 3d x 2 −y 2 and O 2p σ orbitals and not from the atomic single orbital of Cu 3d x 2 −y 2 .This is justified by very large hybridization gap of the AB and bonding (B) or nonbonding (NB) orbitals.See Appendix D of Ref. [12] and the last paragraph of Sec.V in this paper.Here the transfer and interaction parameters are derived at values close to the experimental optimum hole concentration (at 10% doping for CaCuO 2 and Hg1201 and at 20% for the two Bi compounds).This choice is appropriate in this paper, because properties at optimum hole concentration are the central subject.The effective Hamiltonians have the form with Here i, j are the unit cell indices of the maximally localized Wannier function [31,32] and c † iσ (c iσ ) is the corresponding creation (annihilation) operator of an electron of spin σ at the site i.The number operator is given by n i = σ n iσ and n iσ = c † iσ c iσ .Hopping amplitudes t ij in the kinetic energy H kin depend on the relative coordinate vector r i − r j by assuming the translational invariance of the crystal structure.The direct effective Coulomb interaction given by H U is scaled by the on-site interaction U , and off-site interaction H V is the sum over the combination of the site i and j, which is proportional to V ij .Leading values for all of the four materials are listed in Table I.For longer ranged part of t ij and V ij , see Sec.S1A in Supplemental Material (SM) [33].Note that the Hamiltonian parameters for Hg1201 in Ref. [12], which we employ, are improved from Ref. [34].It results in slightly different physical quantities on the quantitative level in the present solution in comparison to Ref. [10].
From four different materials, we learn that the realistic range of available ab initio Hamiltonian parameters is estimated to be 6 (see also Table IV).In this paper, we investigate whether the diversity of the SC properties can be understood within this range of parameters.We note that the effective Hamiltonian parameters in Eq. ( 1) are restricted to a single CuO 2 layer.However, in the case of multilayer cuprates CaCuO 2 and Bi2212, the distance between CuO 2 layers is comparable to the cell parameter along x direction, so that interlayer coupling parameters (given in Table VI in Appendix for Bi2212) also exist in the effective Hamiltonian [12], and its amplitude V l n ≲ 0.6 eV is comparable to that of the intralayer off-site interaction V n ≲ 0.9 eV.This interlayer coupling is ignored in Eq. ( 1), but potentially plays a role in the SC properties.This role is actually examined in Sec.III C 1 in the case of Bi2212, which ensures that the SC order parameter F ∞ SC (defined later in Sec.II C) and physical quantities are essentially not affected by the interlayer coupling within the present case of CaCuO 2 and Bi2212 as we discuss in Sec III C 1. Thus, we restrict to Eq. ( 1) even in the case of CaCuO 2 and Bi2212.We employ this "single-layer approximation" for all the four materials throughout this paper unless otherwise noted.

B. Numerical Methods
We solve the Hamiltonian in Eq. ( 1) by applying the many-variable variational Monte Carlo (mVMC) method [18,19] with a trial wave function of the form |ψ⟩ = P G P J P dh |ϕ pair ⟩.Here we consider the Gutzwiller factor P G = exp(−g i n i↑ n i↓ ) [35], the Jastrow correlation factor P J = exp( i<j α ij n i n j ) [36,37], the doublon-holon correlation factor [38], and a generalized pairing wave function of the form m , and f iσ,jσ ′ .We also supplement the mVMC method with the restricted Boltzmann machine (RBM) [20] and the firstorder Lanczos step to improve the wave function and also to take the zero limit of the variance extrapolation to improve the estimate by following the variance extrapolation method [39][40][41] and by using the simple mVMC, mVMC +Lanczos, and mVMC+RBM results (see Appendices A, B, and Ref. [42] for the detailed procedure.
Competing states with spin and/or charge order or strong fluctuations can be studied by imposing the meanfield order at the initial trial wave function [19].The correlated metallic state without any symmetry breaking can also be studied by using the ground-state wave function of the noninteracting system as an initial state.These initial states are then relaxed to lower the energy by optimizing the variational parameters.If the competing states exist, the optimization leads to multiple locally stable solutions.The true ground state is determined by comparing the total energy after taking the variance extrapolation described in Appendix B and if possible after careful size extrapolation to see the thermodynamic limit.
The computational details are the following.For all numerical solutions of finite-size lattices subsequently presented in this paper, we assumed the antiperiodic-periodic boundary condition on a N = L × L square lattice of length L, where N is the number of sites on the single-layer system by ignoring the interlayer coupling except for Bi2212.For Bi2212, we examine the two-layer system to examine the effect of interlayer coupling, because a unit cell of Bi2212 contains two layers and the interlayer coupling could play roles in the SC.Within a layer, hoppings and interactions were taken into account up to t 9 and V 9 , i.e., up to the 2D distance R = (3,3) in the unit of the Cu-Cu distance, while contributions smaller than 0.001 eV were ignored.Unless explicitly mentioned this is applied throughout the whole paper.We concentrate solely on hole doped cases, where the hole doping is given via δ = 1 − N e /N and N e is the number of electrons in the system.

C. Physical Quantities
The physical quantities discussed in this paper are defined as follows: The total energy per site E/N = ⟨H⟩ /N is calculated after the variance extrapolation as it is summarized in Appendix B. To see whether the state has spin and charge order, we compute the spin structure factor and the charge structure factor where is the spin-1/2 operator and σ σσ ′ is the Pauli matrix (= (σ x σσ ′ , σ y σσ ′ , σ z σσ ′ )).The long-range order is determined whether S s (q)/N or S c (q)/N remains after taking the limit N → ∞.The SC long-range order is measured by the d-wave SC correlation function Here ∆ † d (r i ) describes the order parameter of the form where the d x 2 −y 2 -wave symmetry is included via the form factor f d (r) = δ ry,0 (δ rx,1 + δ rx,−1 ) − δ rx,0 (δ ry,1 + δ ry,−1 ).( 6) Then, we deduce P d (r) over the long-range part, as where r = (r x , r y ) includes all sites within the square (−L/2, L/2] 2 and M is the number of lattice points satisfying √ 2L/4 < |r| < √ 2L/2.We define the SC order parameter in the thermodynamic limit, F ∞ SC by

III. AB INITIO RESULTS
In this section we present our calculated ab initio results on the four families of compounds.We first analyze the results for doped CaCuO 2 in detail and then compare it with the result of doped Hg1201.Doped Bi2201 and Bi2212 suffer from the experimental uncertainty of the apical oxygen position due to the supermodulation of the crystal structure.Since it generates a variance of the Hamiltonian parameters if one assumes the translational symmetry of the Hamiltonian parameters, we show the properties by indicating this range.All of the four materials show dominance of the d-wave SC in the doping concentration dependence and the calculated results reproduce the experimental materials dependence of the strength of SC, which makes it possible to extract the universal properties and systematic trends as well.

A. Doped CaCuO2
Here, we present the results of our calculation for the doped CaCuO 2 in the following order: (1) δ dependence of SC properties; in particular, P ∞ d and the d-wave SC order parameter F ∞ SC .(See the definitions in Sec.II C.) (2) Competition between the SC state and other competing states including stripe states.

Properties of superconducting phase
First, we discuss the r dependence of the pairing correlation P d (r) for L × L lattice and its long-ranged part Pd (L): Fig. 1 shows P d (r) and Pd (L) for several choices of square-lattices with the linear size L from 24 to 36 and hole doping δ from 0.028 (2.8 %) to 0.247 (24.7 %).For each value of δ, we observe that Pd does not significantly depend on L, suggesting the existence of a strong SC long-range order in the thermodynamic limit in this ground-state candidate.This is indeed confirmed by a size extrapolation, i.e., plot of Pd (L) as a function of 1/L to estimate P ∞ d in the limit L → ∞ via linear regression, as shown in the insets of Fig. 1(a)-(h).The linear 1/L scaling was employed in Ref. [10] and is expected to work because of Dirac-type linear dispersion for the In each case, we show P d (r) for the square lattice size L = 24, 30, 36.For the same distance r we employ the largest value of correlation.We perform the same procedure in later figures.Inset of each panel: Size extrapolation of Pd (L) to the thermodynamic limit P ∞ d , whose numerical value is listed in Table II.The gray line shows the linear approximation.Statistical errors originating from the Monte Carlo sampling are smaller than the symbol size.
quasiparticle excitation of the d-wave superconductor at the nodal points.Here, we have shown the data calculated from the transfer and interaction parameters in the Hamiltonian fixed at 10% hole doping for simplicity as we addressed above.However we can test its robustness by taking δ dependent transfer and interaction parameters.The result is shown in Appendix C and the difference is small.
After taking the size extrapolation to the thermodynamic limit, we show the δ dependence of the order parameter F ∞ SC calculated from Eq.( 8) and P ∞ d = lim L→∞ Pd (L) in Fig. 2 and the numerical values in Table II.This shows a rapid increase of F ∞ SC from 0 at δ = 0 up to δ ∼ 0.05 as a function of δ followed by a plateau around 0.05 ≤ δ ≤ 0.1 and monotonic decrease with further increasing δ above around 0.1 in the thermodynamic limit.
The dome structure ubiquitously observed for T c in the cuprates is qualitatively similar to the δ dependence in F ∞ SC , but the peak for F ∞ SC is located at somewhat lower δ ∼ 0.05 than the case of experimental T c , where the optimum δ is observed to be δ ∼ 0.12 [43].We will discuss this discrepancy in Sec.V.However, the monotonic decrease of F ∞ SC with increasing δ for δ ≥ 0.1 is consistent with the universal trend of δ dependence of the SC gap identified from the angle-resolved photoemission spectra (ARPES) and the scanning tunneling microscope (STM) of the cuprates in general [44,45], though the SC gap in the experimental estimate contains an ambiguity associated with the contribution from the pseudogap.In the mean-field picture, the SC gap is the product of the order parameter F ∞ SC and the effective attractive interaction.If we consider the experimentally observed maximum gap ∼ 50 meV [44,45] and F ∞ SC ∼ 0.13, the characteristic scale of attractive interaction is as large as ∼ 0.4 eV.This imposes a constraint on theories for the SC mechanism.The sharp increase of F ∞ SC between δ = 0 and 0.05 and the subsequent dome structure are similar to the earlier study by a VMC method and a cluster dynamical mean field study [9,46] for the Hubbard model, where rapidly increasing F SC from 0 at δ = 0 already reaches F SC ∼ 0.1 at δ = 0.03 in the present notation.In case of the Hubbard model, however, it was argued that the ground state is actually not SC but stripe-ordered states [4][5][6].

Competition of SC, stripe, and AFM states
Now, we analyze the competition between the SC and other states.The energies of the SC state and other states at L = 24 are given in Fig. 3.We see that the SC state has the lowest energy in the region from δ = 0.05 to δ = 0.25, indicating that the SC phase is dominant in the ground state of doped CaCuO 2 .The SC ground state, however, is severely competing with the C4S8, and C3S3like stripe states, and AFM-type state within the energy difference of 5-10 meV.Here, CmSn denotes the charge and spin ordered stripe state with the periodicity of m lattice spacing for the charge modulation and the period n for the spin order.Spin and charge real-space patterns and structure factors are explicitly illustrated in Secs.S2-S4 of SM [33] for the C4S8, C3S3-like, and AFM-like excited states.In the region studied here, we find only C4S8 and C3S3 as candidates of the competing stripe order, which has similarity to an earlier study of the simple Hubbard model with the next-nearest-neighbor transfer in the range 0.2 < |t 2 /t 1 | < 0.3 and 0.05 < δ < 0.25 [6] and an ab initio study [10].Note that C4S8 at δ = 1/8 and C3S3 at δ = 1/6 have a particular commensurability energy gain.In Fig. 3 (b), a dip exists in the AFM states at δ = 0.167.At the moment, the origin of the dip is not clear.that the AFM states up to δ = 0.16 and the C4S8 stripe state at finite doping around δ ∼ 0.12 do have the longrange order as one can see in Figs. 14 (a) and 15, respectively.However, although the initial trial states are ordered mean-field states, long-ranged order seems absent or is very weak after the optimization in VMC calculations for the C3S3 stripe states and seems replaced by well-developed short-ranged correlations at δ ̸ = 0.This is the reason why we add "-like" for the case of C3S3.The AFM and stripe states are in any case excited states of the SC ground states at δ ≥ 0.05.
On the comparison between SC and stripe or AFM states, similar severe competitions were reported in Hubbard models.However, the present ab initio results have a crucial difference, where the charge uniform SC phase is the dominant ground state, while the ground states of Hubbard models irrespective of the presence or absence of t 2 mostly have the stripe long-range order.The reason for this difference originates from the presence of realistic off-site interaction in the ab initio case as was already pointed out in Ref. [10].We discuss this point in the comparison to the Hubbard models in Appendix E

B. Doped Hg1201 Compared with CaCuO2
Now, we present our results in the case of Hg1201, and compare it to CaCuO 2 .In this comparison, we find (i) the positive correlation between The pairing correlation for Hg1201 Hamiltonian and the size extrapolation are shown in Fig. 4 for δ = 0.146, indicating the existence of the SC long-range order.The size of the order parameter is F ∞ SC ∼ 0.09 as compared with ∼ 0.116 for doped CaCuO 2 at δ ∼ 0.12, respectively, which are both close to each optimal concentration.The difference in F ∞ SC between Ca and Hg compounds can be compared with the difference in U/|t 1 | = 8.10 for CaCuO 2 and 7.35 for Hg1201.As we discuss later, ) with 1 ≤ i ≤ 9 in general, but V i /t 1 and t i /t 1 dependencies are weaker as compared to dependence on U/|t 1 | in the realistic parameter range.See Appendix F for the example of |t 2 /t 1 | dependence.See also Fig. 9 for the V 1 /t 1 dependence.In both cases, the change in F SC at the optimum doping is at most 10% in the realistic parameter range.In fact, in the comparison of Bi2212, Bi2201, CaCuO 2 , and Hg1201, |t 2 /t 1 | is ∼ 0.30, 0.27, 0.25 and 0.20, respectively, which does not have systematic correlation with T opt c .Appendix F shows tiny anticorrelation of |t 2 /t 1 | and F SC , but is practically negligible at the optimal doping.After careful examination of other parameters as well, the difference of F ∞ SC in these four compounds studied is concluded to be ascribed to the difference in U/|t 1 |.
The materials dependent F ∞ SC may also be compared with the difference in T opt c ∼ 110 and 90 K for CaCuO 2 and Hg1201, respectively, because T c may be proportional to the order parameter F ∞ SC .Since T c has the dimension of energy and should also be scaled by the overall characteristic energy scale t 1 , T c may be proportional to |t 1 |F ∞ SC .In fact, the ratio of T opt c /|t 1 |F ∞ SC as a non-dimensional quantity is ∼ 0.16(1) at the optimal doping δ ∼ 0.12 for CaCuO 2 and ∼ 0.16(2) at the optimal point δ ∼ 0.15 for Hg1201 as we show in Table III, supporting the hypothesis that T opt c is universally given from the relation at the optimal doping.
In the Uemura plot [47], it was observed from the muon spin resonance (µSR) measurement that T opt c is proportional to the ratio between the superfluid density n s , here interpreted as F ∞ SC / √ 2, and the effective mass m * .Since the mass enhancement from the bare band mass m 0 , namely m * /m 0 at the optimal hole density may be similar in the cuprates, T opt c is indeed expected to be roughly proportional to |t 1 |F ∞ SC according to the Uemura plot, because the inverse band mass is essentially determined by the dominant transfer t 1 .In addition, n s estimated from the relaxation rate σ ∼ 2 µs −1 from the µSR measurement for the cuprates with T c ∼ 80 − 100 K corresponds to n s m 0 /m * ∼ 4 × 10 21 cm −3 .If we cut out the volume including one Cu atom with the c axis length ∼ 6 Å as in [47] irrespective of the unit cell volume to compare with F ∞ SC defined as the value per Cu atom, this corresponds to F ∞ SC ∼ 0.10 by considering the definition F ∞ SC = √ 2n s and m * /m 0 ∼ 5 assumed in Ref. [47].Then it is also quantitatively consistent with the present result of F ∞ SC ∼ 0.10 at the optimum doping.These indicate that our results indeed capture the realistic situation more or less quantitatively.We will show in Sec.V that F ∞ SC also agrees with the estimate from the angle resolved photoemission spectra (ARPES) in case of Bi2212 and Bi2201.

C. Doped Bi2201 and Bi2212
In this subsection, we discuss ab initio results for Bi2212 and Bi2201 and compare them each other.Unfortunately, in these two compounds, an uncertainty exists in the experimental crystal parameters that causes the uncertainty in the effective Hamiltonian parameters as well.Especially, the distance d z Oap between an apical oxygen and the nearest Cu atom is not fully precisely determined and the available experimental data have considerable variations [48][49][50][51][52][53].This uncertainty is also re- for doped CaCuO2 at δ = 0.125, and Hg1201 at δ = 0.146 (these values of δ are chosen in accordance with those closest to experimental optimal values, δ ∼ 0.12 and 0.15, respectively, and to allow the numerical size extrapolation easier).Note that the ratio T opt c /(|t1|F ∞ SC ) is given as a nondimensional quantity by using 1 eV = 1.16 × 10 4 K.The parentheses in the last digit indicate the error bar.lated to the structural distortion and long-period modulation of the CuO 2 plane arising from the effect from the BiO layer [54,55] as we discuss in Sec.V. Recent ab initio studies have clarified that this uncertainty leads to a possible variety of effective Hamiltonian parameters, especially owing to the variation of the apical oxygen position [12].
In principle, the structural optimization in ab initio calculations is desired to predict the stable atomic position.However, such an optimization in strongly correlated electron systems is at the moment not necessarily accurate enough and we leave this task for future studies.Instead, in this paper, we admit a range of Hamiltonian parameters and discuss the consequence.
As analyzed in Ref. [12], the apical oxygen position sensitively affects the effective Hamiltonian parameters, primarily the value of U .For Bi2212, the value U ∼ 4.2 eV in Table I I, namely U = 4.2 eV for Bi2212 and U = 4.4 eV for Bi2201 and then discuss in Sec.III C 3 the possible range of SC properties originating from this uncertainty later.

Bi2212
We begin with the results for Bi2212.Figure 5 shows P d (r) and Pd (L) at δ = 0.167 for L from 16 to 36 by switching off the interlayer transfers and interactions.Namely we first show the results obtained by solving the single-layer Hamiltonian despite the actual two-layer unit cell of Bi2212.The case of δ dependent Hamiltonian is also shown in Appendix C, where the difference is small.Similarly to CaCuO 2 , we identify a SC ground state with large Pd and it does not change significantly by increasing L, where the size extrapolation gives P ∞ d ∼ 0.0151 as shown in the inset of Fig. 5, which corresponds to a SC order parameter F ∞ SC ∼ 0.12.This relatively strong value of all four considered compounds.Again the enhanced F ∞ SC originates from the larger U/|t 1 | in accordance with the observation in the comparison of Hg1201 and CaCuO 2 in Sec.III B. This large F ∞ SC is also consistent with the high T c (up to ∼ 100 K) [30].We will discuss more intricate aspect in Sec.V.
The competition with other phases is seen in δ dependence of the total energy shown in Fig. 19 in Appendix I. Similarly to CaCuO 2 , the SC state is the ground state in most of the doping concentration, while it is severely competing with spin and charge ordered states.Now, we extend the calculation by switching on the interlayer terms and solve the two-layer Hamiltonian obtained in Ref. [12] to examine the effects of interlayer coupling.For the calculations we take two identical layers (in terms of intralayer parameters for t ij , V i , and U where i and j are intralayer combination), coupled by the interlayer terms listed in Table VI (Appendix H).The interlayer contributions are restricted to the leading interlayer hopping term of size t l 0 = −0.098eV and nearest-and next-nearest-neighbor interlayer interaction of size V l 0 = 0.643 eV and V l 1 = 0.463 eV.See Ref. [12] for more details.A comparison between the single-and two-layer cases of P d (r) for L = 16 and δ = 0.167 is shown in Fig. 6.We see that the long-range average of P d (r) in the two-layer case is close to that of P d (r) in the single layer case, which demonstrates that P d (r) is not significantly affected by the interlayer Coulomb interaction and hopping parameters.Indeed, for the single layer we found a long-range average of the SC correlation function of P single d = 0.0225, while for the two-layer compound the average is P two d = 0.0213.The corresponding values of the SC order parameter are F single,∞ SC = 0.150 and F two,∞ SC = 0.146, which differ by only ∼ 2.8%.The essentially same behavior between the single-and twolayer cases may not depend on the system size in accordance with the result of a two-layer Hubbard model at the optimum doping [56].This similarity may be attributed to (i) the relatively small leading interlayer hopping parameter of t l 0 = −0.098eV and also (ii) the robustness of the SC solution against the interlayer Coulomb interaction parameter, because the pairing occurs essentially only within a layer.

Bi2201
In the case of Bi2201, we again solve the single-layer Hamiltonian.The competition with other phases are seen in δ dependence of the total energy shown in Fig. 20 in Appendix I. Figure 7 shows P d (r) and its extrapolation to the thermodynamic limit, which suggests the stable longranged SC order.Again, the obtained value F ∞ SC ∼ 0.10 is consistent with the rule that larger U/|t 1 | leads to larger F ∞ SC because U/|t 1 | is the second largest among the four materials in the estimate shown in Table I.The smaller U/|t 1 | relative to Bi2212 leads to weaker SC.However, on a more quantitative aspect, we need to be careful about the uncertainty of the Hamiltonian parameter.We will discuss this issue below.

Effect of structural uncertainty on possible variation of SC properties
Since we have the uncertainty of the Hamiltonian parameters particularly for the interaction as we discussed above, we here monitor the effects of modifying the effective interactions for Bi2212 and Bi2201, which well represent the effect of variant apical oxygen position as shown in Appendix C of Ref. [12].Namely, Table I with preserved transfer parameters fixed at each ab initio value, together with interaction scaling represents most of the effect of the apical oxygen shift and we scale the Hamiltonian (1) such that where α = ξ = 1 represents the ab initio case given in Table I, α scales the on-site Coulomb interaction term H U , and ξ scales the remaining off-site interactions V i .Since the apical oxygen shift alters the interaction parameters in the way α ∼ ξ [12], we examine the dependence on α = ξ below.Figures 8 (a I and when we admit the uncertainty range of the interaction parameters for Bi2212 and Bi2201.The scaling Eq.( 9) proposed for CaCuO 2 and Hg1201 is also valid in the Bi compounds and the experimental T opt c is within the inferred range.We conclude that its materials dependence is well captured for the four studied materials (see also Fig. 11 in Sec.VI).
We realize that the fragility and diversity of T opt c experimentally observed in the range 10 < T opt c < 40 K for Bi2201 is accounted for by the range of actual apical oxygen position.This range may be caused by the type of dopant atoms, impurities, and the spatial inhomogeneity caused by the supermodulation, which may depend on samples and the amplitude of the modulation.In fact, it was observed that the d z Oap periodically varies as much as 6% in accordance with the supermodulation for Bi2212 and a comparable modulation may exist for Bi2201 as well, which can be the origin of the experimental uncertainty [54,55].The basic origin of this diversity and relatively low T c among the four families of compounds is attributed to relatively small U/|t 1 | in the lower uncertainty range, at which the SC order becomes sensitively damaged by a slight decrease of U/|t 1 |.We discuss more general aspects of the interaction dependence in Sec.IV.Even when we admit the uncertainty range, the general trend about the weaker SC of Bi2201 than those of Bi2212 is well explained by this ab initio result.It can also safely be addressed that Bi2212 has one of the strongest SC and T opt c among the four families comparably to CaCuO 2 .The effects of apical oxygen position on Hamiltonian parameters are discussed in Sec.V and in Appendix L.

IV. RESULTS BEYOND AB INITIO: INTERACTION DEPENDENCE OF SUPERCONDUCTING ORDER
We now study SC properties beyond the ab initio results.Ab initio results in the previous section successfully reproduce the experimental trend and have revealed that U/|t 1 | is the principally important Hamiltonian parameter to control the SC order parameter.Therefore, it is intriguing to examine the optimum Hamiltonian parameters to maximize the order parameter and hence the optimum T c beyond the existing materials for the purpose of materials design to seek for higher T c superconductors.We present U/|t 1 | dependence of F SC as well as dependence on off-site interaction when tuning the interaction parameters artificially away from the ab initio value while keeping transfer parameters fixed at ab initio values.We here take an example of CaCuO 2 Hamiltonian at δ = 0.167 and monitor the effect of α and ξ dependencies defined in Eq. (10).
We examine three types of scaling to go beyond the ab initio Hamiltonian: (i) Scale only the on-site Coulomb interaction by α with fixed ξ = 1.0 (H(α, 1.0)), (ii) Scale the full interaction part equally by using α = ξ (namely, H(α, α)), and (iii) Fix α = 1.2 and scale the off-site Coulomb interaction uniformly via ξ by employing H(1.2, ξ), by considering the fact that (ii) shows the maximum SC order at α = ξ = 1.2.For the cases (i) and (ii) we chose scaling values α ranging from 0.6 up to 4.0 while for the third case the range from 0.0 to 2.0 is chosen.Since the size dependence is not appreciable, we study L = 24 lattice.The SC order parameter F SC are shown in Fig. 9 (a) at δ = 0.167 hole doping.
The results show that the SC order parameter can be enhanced with the amount of around 30% from the ab initio value when the interaction parameter is tuned to α ∼ 1.2 for (i) and around 20% at α = ξ ∼ 1.2 for (ii), which may allow T opt c as much as ∼ 130-140 K, when compared to the ab initio results for CaCuO 2 .
In the tuning (iii), we find that the order parameter further increases up to F ∞ SC ∼ 0.22 by decreasing ξ, which is twice as large as the ab initio case.However, we keep in mind that on-site and off-site interactions cannot independently be controlled in the usual experimental conditions.The present result offers a guide to enhance the SC in designing artificial structure and metamaterials including surface and interface, where quicker screening of the off-site interaction is desirable by keeping the on-site interaction at the optimal value (in this case α = 1.2).Another limitation to be considered is the competition with the stripe and AFM order.As far as we restrict the on-site interaction within the ab initio range, the SC energy is always lower than that of the stripe state as one sees in Fig. 9 (b).Here, we show the competition with the C4S8 because it is established that the most severe competitor is the C4S8 state.However, the decreasing of energy difference, such as at (α = 1.2, ξ = 0), where the difference is ≤ 2 meV, may result in the thermal destruction of the SC order.For ξ ≥ 0.2 the SC is still a stable ground state, while at small ξ, F SC becomes nearly twice of the ab initio value for the doped CaCuO 2 .
The order parameter decreases when U/|t 1 | is too large beyond the realistic range of the cuprates studied in this paper as we see in Fig. 9 (a).This reduction was already pointed out on the level of the Hubbard model [9].The reduction at large U/|t 1 | is studied in Appendix J, which shows nontrivial power-law dependence of F SC on U/|t 1 |.The reduction itself may be easily understood on a qualitative level from the suppression of charge fluctuation with increasing U/|t 1 |, which also suppresses the quantum entanglement caused by the suppression of both spin singlet fluctuation and dynamical exciton generation as was reported in the literature [9,57].It was pointed out that the enhanced quantum entanglement can be achieved by the fractionalization of electrons [57,58], which may be maximized at the optimum U/|t 1 |.

V. DISCUSSION
By assuming the same ratio T c /(|t 1 |F ∞ SC ) ∼ 0.16 with CaCuO 2 and Hg1201, we can infer the range of T c arising from the uncertainty of the apical oxygen position and resultant uncertainty of U/|t 1 | listed in Table IV.The range of inferred T c listed in Table IV for Bi2212 is consistent with the experimentally observed range of T opt c within the error bars.This suggests that the sample dependence of T opt c may be accounted for by the sample dependence of the apical oxygen position.
We note that the order parameter F ∞ SC ∼ 0.10 obtained here for Bi2212 and Bi2201 also shows consistency with the result obtained by using the machine learning of the angle-resolved photoemission spectroscopy data for Bi2212 and Bi2201, respectively, at the optimum dop- ing [59], which gave ⟨c k↑ c −k↓ ⟩ ∼ 0.065 at the antinodal point for Bi2212 and the momentum averaged value ∼ 0.063 for Bi2201 at the optimum doping.These are translated commonly to In the case of Bi2201, the estimated 0.16|t 1 |F ∞ SC is also listed in Table IV.The comparison with the sample dependence of experimental T c suggests that the true apical oxygen position is distributed near the lower bound d z Oap ∼ 2.45 Å, if it is spatially uniform.Alternatively if the supermodulation exists, the lower bound of d z Oap in the modulation may be close to 2.45 Å, because it governs the SC order as the bottleneck.It is desired to test this inference by precise and simultaneous measurements of the relation between T c and d z Oap for Bi2212 and Bi2201.The order parameter F ∞ SC increases with decreasing hole doping for δ > 0.05 as we find in Fig. 2, which follows the same trend as the SC gap as we discussed, but is slightly different from the dome structure known for T c in the cuprates, where the peak of the dome is located at higher δ.Complete and quantitative understanding of this different trend is not the scope of this paper and is left for future studies.However, the origin of this difference can be inferred to be attributed to the increasing damping and incoherence of electrons in the underdoped region toward the metal-insulator transition as was analyzed before [44,45], which is represented as the enhanced self-energy of the normal electrons toward the Mott insulator.In the experimental conditions, the atomic doping/substitution introduces atomically spatial inhomogeneity, which is ignored in the present study and may cause the localization of carriers, which could also be the origin of slower increase of T c upon doping.
In Appendix K, we show a qualitative difference of the momentum distribution between the optimal and underdoped hole concentrations, which suggests a signature of the increased damping at lower carrier concentration within the present ab initio study.We also discuss in Appendix K the subtlety and complexity in the underdoped region due to the pseudogap formation, involved in several quantities, which is not the scope of this paper.
The δ dependence of the energy decomposed to the kinetic, on-site and off-site interaction energies are analyzed in Appendix G.It should be noted that the on-site interaction energy E U = ⟨H U ⟩ has a convex curvature as a function of δ, which contributes to the effective attraction of electrons despite the original strong repulsion U .In fact, the U/|t 1 | dependence of the local attraction is qualitatively consistent with F SC in Fig. 9. See Appendix G for more details.The local effective attraction may cause the Cooper pairing as well as the stripe formation.It was pointed out that the convex curve of the local energy generates bistable excitations, one in the underdoped side and the other in the overdoped side, inducing the electron fractionalization and the enhanced quantum entanglement through the quantum tunneling of the two excitations [57,58].This line of further research is an important future subject.
The importance of the apical oxygen position has been pointed out from various viewpoints [60][61][62][63].In this paper, we have elucidated the crucial role of controlling U in general in the single-band description, which quantitatively explains the variation of T opt c and its uncertainty in the Bi compounds.In addition, the modulation of the SC gap with the modulation of the apical oxygen position in Bi2212 has indicated that the longer d z Oap induces the smaller SC gap [55].This is consistent with the trend of F ∞ SC , which decreases when α is increased in the realistic range as shown in Fig. 8.This indicates that Bi2212 is located already slightly above the peak in the α dependence of F ∞ SC , which corresponds to α = 1.2 for doped CaCuO 2 shown in Fig. 9.The present observation is also in accordance with the effect observed by laser irradiation aiming at the displacement of the apical oxygen position [64].The control of d z Oap if possible in a spatially uniform fashion may help to optimize the SC in which the disturbance and pair breaking by the randomness caused by the inhomogeneous supermodulation in the case of the Bi compounds could be avoided.
The effects of the displacement of the apical oxygen position on the Hamiltonian parameters were discussed in Ref. [12].In Appendix L, we readdress this issue in relation to earlier work.
We have mainly focused on the quantities at optimal doping to clearly extract the diversity of the materials dependence, where the experimental subtlety due to the effects from extrinsic randomness as well as the complexity of physical quantities arising from the pseudogap formation in the underdoped region is irrelevant.The behavior of suppressed T c , the SC carrier density n s , and the coherent spectral weight arising from the pseudogap formation in the underdoped region are left for future studies.See Appendix K for more details.
In general, the d-wave SC correlation has a dip at (1,1) distance.The origin is speculated as follows: The singlet pair between electrons at (0,0) and (1,0) sites dynamically interfere with the singlet pair between (1,1) and (1,0), because these two singlets are incompatible.The same is true for the singlets, which share (0,1) site.This double interfere makes the SC correlation of the pair between (0,0) and (1,1) smaller very generally.However, at the moment, we do not know the origin of particularly large dip around δ = 1/8, which is left for future study.
There exist studies proposing the possibility of the pair-density-wave (PDW) states [65].However, in the present ab initio Hamiltonians, the PDW correlation remains small as is shown in Fig. S7 of SM [33], where peak is absent at nonzero momentum.If the stripe longranged order coexists in the SC ground state, the PDW order must be trivially accompanied.However, as one sees in Figs.14(b) and S6, the stripe correlation exists but remains small and is scaled to zero in the thermodynamic limit.
In this paper, we have employed the single-band Hamiltonian for the AB orbital, because the hybridization gap between AB and B or NB orbitals is so large (∼ 8 − 9 eV) that the B and NB bands are more or less completely filled and inactive (see Fig. 10(b) of Ref. [34] and Appendix D of Ref. [12]).When one starts from the three bands constructed from Cu 3d x 2 −y 2 and O 2p σ atomic orbitals, the analysis would be more complicated.See Appendix M for this issue.

VI. SUMMARY AND CONCLUSIONS
We have studied the superconductivity in the ab initio Hamiltonians for CaCuO 2 , Hg1201, Bi2201, and Bi2212 derived by using the experimental crystal structure in Ref. [12] without adjustable parameters .The dominance of SC order against severely competing stripe states and AFM state in a wide range of hole concentration is shown in the solutions for the ground state of all four materials obtained from the variational Monte Carlo calculations, which agrees with the experimental results.
The SC order parameter F ∞ SC at the optimal doping shows consistency with the superfluid density measured in the µSR and the machine learning analysis of the ARPES data for Bi2212 and Bi2201.F ∞ SC decreases with increasing doping for the doping concentration δ > 0.05, showing a similarity to the SC gap reported in the STM and ARPES measurements.On the other hand, F ∞ SC quickly decreases to zero toward δ = 0 for δ < 0.05 forming a dome structure which has a similarity to experimental T c , but the dome peak appears at slightly lower δ for the calculated F ∞ SC .This may be attributed to the reduced renormalization factor suggested by the broadened momentum distribution.
From the comparison of the four materials, we have revealed that U/|t 1 | is a crucial parameter to control the strength of the SC order; larger U/|t 1 | materials show larger SC order parameter F ∞ SC in the realistic materials.This explains that T c and the SC gap at the optimum doping are larger for CaCuO 2 than Hg1201, where T opt on-site and off-site interaction independently, further optimization of the SC order parameter as much as the factor 2 larger value beyond the available compounds synthesized so far without falling into other competing states can be achieved as the theoretical maximum value in the present mechanism.By increasing |t 1 | as well as the whole parameter values uniformly, T opt c should obviously increase accordingly.These offer a clue for the materials design in the future.Here we discuss the stability of the AFM, C3S3, and C4S8 states in the thermodynamic limit for CaCuO 2 .To do so, each state is stabilized on different lattice sizes and the spin and charge structure factors (S s (q) and S c (q)) are calculated.We plot S s (q) as a function of 1/L or 1/L 2 depending on the cases of the presence or absence, respectively, of the AFM order by following the convention and perform a linear extrapolation.The result is shown in Fig. 14, where in Fig. 14(a) the peak height of the spin structure factor follows a linear trend with a nonzero offset, indicating a stable long-ranged AFM order.In Fig. 14(b), the peak of the spin structure is scaled to zero in the thermodynamic limit, which indicates the absence of the long-ranged AFM order in the SC state at this doping.We did not find the coexistence of the SC and AFM order at other doping, either.The size scaling of the charge-stripe states is shown in Fig. 15.Although the plot of the size dependence is not sufficient due to very demanding computational cost for larger sizes, the trend suggests that only the C4S8 state shows a clear long-ranged order of spin and charge at around δ = 0.125 in the thermodynamic limit.The C3S3 state seems to collapse to a paramagnetic state at the chosen doping of δ = 0.207, at which C3S3 has relative stability.Hence we use the labeling "C3S3-like" instead of C3S3 in the main text.We do not go into details of the size dependence for the stripe-like states, because they are in any case excited states.Here we discuss the crucial difference of the ab initio results from the extensively studied Hubbard models without the off-site interaction.In the Hubbard model studies, irrespective of the presence of the next-nearest transfer t 2 or absence of it, in the hole doped region, a broad consensus seems to be formed, where the ground state has stripe type long-ranged charge order [3,5,6,[67][68][69][70][71][72].The charge and spin stripe states were also suggested to coexist with weak SC order in some cases [5,72], but other studies did not find the SC order [6,68,69,71].
In the Hubbard model studies, the numerical methods have a variety including the present VMC, density matrix renormalization group, constrained path quantum Monte Carlo, tensor network and density matrix embedding theory, which have their own advantages and disadvantages and they are complementary.When they can be compared with reliable solutions, the above Hubbard model calculations were benchmarked, which have shown comparable accuracy when compared between each state-ofthe-art updated version.Even in quantum spin models, the level of accuracy of the above methods is roughly similar [21].We also refer to recent thorough comparisons [73].
On the other hand, when realistic off-site interaction is included, the ground state is reported to be charge homogeneous superconducting state [10].The off-site interaction substantially suppresses the SC order but the charge and spin stripe states are more damaged by the frustration introduced by the off-site interaction.The role of off-site interaction for the stabilization of the SC state relative to the stripe state was directly demonstrated in Ref. [10] by the comparison of ab initio result and that of the Hamiltonian obtained by switching off only the off-site interactions.We confirmed the similar behavior for the present Hamiltonians.The absence of developed stripe correlations in the SC ground state is seen in Fig. 14 herein and in Fig. S6 of SM [33].
Appendix F: t2 dependence Here we show t 2 dependence of F SC by starting from the ab initio Hamiltonian for hole doped CaCuO 2 with other Hamiltonian parameters fixed, where the effect of t 2 is monitored beyond the ab initio value primarily within the realistic range of 16 and 17 for L = 24 lattice.F SC slightly decreases with increasing t 2 , which is qualitatively consistent with a different approach [63].However, the variation of F SC is at most 10% near the optimum doping in the realistic range.Furthermore, F SC has essentially no t 2 dependence in the range 0.0 ≤ |t 2 /t 1 | ≤ 0.2.On the other hand, the period of the stripe order is known to sensitively depend on t 2 [5,6,68,69,71,72] in the ground states of Hubbard models and it may alter the superconductivity if it coexists, while the present charge uniform  superconducting ground state is quite different.V.
The result shows that the total energy is concave as a function of δ of course, which is required from the thermodynamic stability, while only E U exhibits convex behavior with c 2 < 0. Because the effective particle interaction is given by the δ 2 term, we find that the local quantity E U is the origin of the attraction while the non-local energies contribute to c 2 > 0 in the total energy ensuring the thermodynamic stability.This supports the idea that the local strong correlation (repulsion) called Mottness turns to originate the effective attraction of the electrons, which is the underlying mechanism of both of the Cooper pairing and charge segregation such as the stripes.The attraction is understood from the following Mottness: The local energy is retained high in the Mott insulator because of U .However, if the carrier is doped, this is released by acquiring the itinerancy in a nonlinear fashion as a function of δ which yields c 2 < 0 and the attraction.In accordance with the α dependence of F SC , c 2 shows a similar behavior as −c 2 = 1.47, 2.24, 2.13, 1.85 and 1.53 at α = 1.0, 1.1, 1.2, 1.3 and 1.4, respectively with a peak around α ∼ 1.1-1.2.Bi2212 The effective interlayer Hamiltonian parameters for Bi2212 are shown in Table VI, where V l n is the interlayer Coulomb interaction, and t l n the interlayer hoppings.As defined in Ref. [12], the notation n = 0 represents the interaction or hopping between interlayer nearest-neighbor Cu atoms (located one above the other), and n ≥ 1 represents the interaction or hopping between a Cu atom and its (n + 1)th interlayer nearest-neighbor (located above or below its intralayer nth nearest-neighbor).F SC and the SC gap ∆ SC against the naive expectation may be another example.Although F SC and the SC gap ∆ SC grow on top of the pseudogap as we revealed in the case of F SC , n s seems to be severely suppressed by the pseudogap around the antinodal point.This trend is indeed seen in the comparison of the muon penetration depth and the SC gap measurement [45,[74][75][76], which causes difficulty in the comparison of our calculated result and experimental indications in the underdoped region.Such a complexity is expected to be small at the optimum doping region, while the prominent materials dependence is seen most prominently at the optimum doping.This is the reason why we focus on the materials dependence at the optimum doping.
In Table VII we show the Hamiltonian parameters of Bi2201 when the apical oxygen is artificially shifted.As is already addressed in Ref. [12], U decreases with decreasing d z Oap because of increased screening from electrons at the apical oxygen, which is consistent with the claim in Ref. [61].In Ref. [60], |t 2 /t 1 | increases with increasing d z Oap , for instance, in the comparison of Hg1201 and Bi2201.However, it shows opposite trend in Table VII and in the results in Ref. [12].References [60,61] focused on the effect of orbitals whose energies are above the Fermi level such as Cu 4s and apical O 2p z orbitals.

The increasing d z
Oap makes those levels lower because farther distance to the negative CuO 2 layer charge for the apical O 2p z orbital and farther distance to the negatively charged apical oxygen for the Cu4s orbital makes the Madelung potential lower.This lowering induces a larger hybridization with the AB orbital located around the Fermi level constructed from Cu 3d x 2 −y 2 and in-plane O 2p σ orbital (our target band), which we call in-plane CuO 2 AB orbital.This increased hybridization especially enhances t 2 .However, this is not the whole story.The effective next-nearest-neighbor hopping t 2 in the singleband Hamiltonian is also altered by the effect from the bands below the Fermi level such as d z orbital.Increasing d z Oap makes the lowering of the d z level and hence causes the decrease in the hybridization with in-plane CuO 2 AB orbital which cancels the effect of the increased hybridization of the orbital above the Fermi level as was pointed out in Ref. [63].More precisely, the apical O 2p z and Cu 3d z orbitals are strongly hybridized and they form bonding and antibonding bands and we need to consider all of these contributions, which are taken into account quantitatively in our calculations in the derivation of the effective Hamiltonian.As a consequence, the present estimate for t 2 /t 1 has large difference from that in the estimates of Ref. [60] in some cases.For instance, in the case of Hg1201, |t 2 /t 1 | ∼ 0.36 for Hg1201 in Ref. [60], while ∼ 0.20 in the present study.Although Ref. [60] is based on complex approximations to estimate t 2 /t 1 only by taking into account the contribution from the band above the Fermi level, recent standard way employs the maximally localized Wannier orbitals and its Hamiltonian matrix elements for the estimate of hopping, by considering all the bands contribution near the Fermi level.This is much simpler, more straightforward and transparent for the estimate of the lattice fermion Hamiltonian parameters, which are used in Ref. [12] as the basis of our VMC calculations.In fact, recent estimates of |t 2 /t 1 | for Hg1201 are 0.20, and 0.23 in Refs.[77,78] (in the revised manuscript), respectively, which are consistent with the present 0.20.
As we already mentioned in Appendix F, larger |t 2 /t 1 | slightly but quantitatively suppresses F SC in the realistic parameter range, which is opposite to the prediction in Ref. [60] but is consistent with Ref. [63].Furthermore, and most importantly, too small dependence of optimal F SC on |t 2 /t 1 | clarified in this paper makes the role of |t 2 /t 1 | on F SC highly questionable.We find that the effect of the apical oxygen position on the superconductivity is primarily to control U .

Appendix M: Comparison to approach using multiband Hamiltonian
There exists recent work based on the atomic orbitals containing Cu 3d x 2 −y 2 and O 2p σ orbitals [79][80][81].Of course, multiorbital Hamiltonians should give essentially a similar answer if the derivation and the solving procedure are appropriate.On the other hand, the Hamiltonian becomes more complex with larger number of paramters, as it should be.
However, the AB band and NB or B band are well separated with the hybridization gap (band center separation is ∼ 9 eV and the direct gap is 5-6eV for the cases we studied in this paper).See Appendix D of Ref. [12] for detailed analyses.In this circumstance, we can safely start from the picture of single-band Hamiltonian derived from the AB band of strongly hybridized Cu 3d x 2 −y 2 and O2p σ orbitals only, because the B orbitals are more or less completely filled and inactive.See Appendix D of Ref. [12].See also Fig. 10(b) of Ref. [34], where the completely filled B bands are verified for the Hg-based cuprate through all the relevant hole densities and this is universal in the curate superconductors.In the single AB band description, the B degrees of freedom are downfolded and give the renormalization to the AB orbital description.Since the AB-B hybridization gap is large, the perturbative downfolding procedure to renormalize and eliminate B and NB degrees of freedom works well as a good approximation [12].This is based on the multiscale ab initio scheme for correlated electrons (MACE) with refined GW approximation supplemented by the level renormalization feed back [34].Except for the AB band, all the bands are either well below or above the Fermi level so that they can be perturbatively taken by the partial trace summation to give renormalizations to the AB degrees of freedom.
Of course, one can start from the three-band Hamiltonian where the charge transfer gap and covalency are relevant parameters.However, it should end up with this AB/B description after the basis transformation if one focuses on the low-energy physics in the realistic situation of the cuprates.The effect of the parameters of the charge transfer energy and the d-p covalency were taken into account in our downfolding procedure from the three-band to a single-band AB Hamiltonian in Ref. [12].For instance, larger charge transfer gap results in poorer screening and larger correlation (U ) as is confirmed in Ref. [12] [see the comparison of Table IV with Tables I and II in Ref. [12]].Therefore, the three-band parameters are encoded in U/|t 1 | and other parematers in the AB single-band description indirectly and systematically in a complex manner.
There exist several recent analyses based on multiband Hamiltonian containing Cu 3d x 2 −y 2 and O2p σ atomic orbitals for the cuprates or Hubbard-type models [79][80][81].Three-band Hamiltonian constructed from the Cu 3d x 2 −y 2 and O2p σ atomic orbitals derived and listed in Table IV of Ref. [12] shows rough consistency with the proposal by Ref. [79], in which the authors claim stronger superconducting order for smaller charge transfer gap ∆E xp and larger d-p transfer t xp in the notation of Ref. [12].Naively, one would expect that smaller ∆E xp makes stronger screening on the AB band and decreases effective U in the single-band picture, while larger t xp directly leads to larger t 1 .Both result in smaller U/|t 1 |, which appears to contradict the statement claimed in the present paper that larger U/|t 1 | leads to larger superconducting order parameter in the realistic parameter region.
However, one needs to be careful about the parameter region employed in Ref. [79].When one sees U d (the direct on-site repulsion between atomic d orbital) dependence of the superconducting order parameter in Fig. 2 of Ref. [79], one clearly finds that the superconducting order parameter decreases from U d = 10 to 18 or 14 in the energy scale of t pp , transfer between neighboring O 2p orbital.When one compares the parameters with those in Table IV of Ref. [12], and compares with Tables I and II of Ref. [12], one notices that U d = 10 in Ref. [79] already corresponds to the region around the optimum of U/|t 1 | in the present single-band description and further increase of U d drives the system into the strong coupling region with larger U , where the superconducting order decreases with increasing U/|t 1 | as one can see in our result in Fig. 9(a).This perfectly explains the U d dependence between U d = 10 and 18 in Ref. [79] as well as the t pd = t xp and p level, ϵ p dependences, because larger t pd = t xp and smaller ϵ p = ∆E xp both make smaller U/|t 1 | as was mentioned above and leads to larger superconducting order parameter in the strong-coupling region as is clarified in Fig. 9(a).It clarifies that the region studied in Ref. [79] is outside of the real materials dependence of the cuprates studied here.We further need more detailed studies of the correspondence between multiorbital and the present single-band description in real materials systematically, which is left for future studies.

C3S3
Using the same analysis as above the C3S3 state with charge period three and spin period three can be identified.The spin-and charge structure factors, depicted in Fig. 26 (a) or (b), show peaks at q CDW = (2π/3, 0) and (4π/3, 0) for the charge structure factor, and q SDW = (2π/3, π) and (4π/3, π) for the spin structure factor.In (c) the corresponding momentum distribution is also shown.Fig. 27 (a) or (b) shows the real space configuration of spin and charge with the period of three in r x direction.Compared to the C4S8 state, the C3S3 state realized here has significantly smaller peaks at q CDW in S c (q).As shown in Appendix C of the main text, this may already  For the AF state at δ = 0.167 hole doping on a L = 24 square lattice for CaCuO 2 , the spin-and charge structure factor, and momentum distribution is shown in Fig. 28.In (a) the strong peak at q = (π, π) in S s (q) indicates a strong antiferromagnetic correlations, corresponding to the well known checkerboard spin configuration pattern.For the charge structure factor on the other hand (shown in (b)), there are no sharp peaks.Over the full range, the function is smooth, resulting in a charge homogeneous state.The momentum distribution n(k) is shown in (c), suggesting non Fermi liquid behaviour.

S4. Uncertainty of the apical oxygen position in Bi2212 and Bi2201
SC correlations P d (r) as a function of distance r for several choices of the scaling factor α = ξ with the Hamiltonian H(α, α) are presented in Fig. 29 for Bi2212 and Bi2201.

FIG. 1 .
FIG.1.SC correlation function P d (r) for CaCuO2 at different hole dopings: (a) δ = 0.028, (b) δ = 0.045, (c) δ = 0.087, (d) δ = 0.101, (e) δ = 0.125 (here L = 28 instead of L = 30), (f) δ = 0.167, (g) δ = 0.208, (h) δ = 0.247.In each case, we show P d (r) for the square lattice size L = 24, 30, 36.For the same distance r we employ the largest value of correlation.We perform the same procedure in later figures.Inset of each panel: Size extrapolation of Pd (L) to the thermodynamic limit P ∞ d , whose numerical value is listed in TableII.The gray line shows the linear approximation.Statistical errors originating from the Monte Carlo sampling are smaller than the symbol size.

FIG. 2 .
FIG. 2. The SC order parameter FSC as a function of δ for doped CaCuO2.The gray filled circles show the values of FSC(L) at L = 24 square lattice calculated from Pd (L) shown in Fig. 1 by using FSC(L) = Pd (L) .The red squares are the size extrapolated values F ∞ SC calculated from P ∞ d .Inset shows the corresponding Pd (L) at L = 24 and P ∞ d .
Figures 14 and 15  in Appendix D show size dependence of spin and charge structure factors.Although the demanding computation cost does not allow larger system calculation, the available size dependence supports

FIG. 3 .
FIG. 3. (a) Variance extrapolated energies of doped CaCuO2 for various ground state candidates, SC, charge and spin stripe states C3S3 and C4S8, and AFM state as a function of hole doping δ on a L = 24 square lattice.All energies are subtracted by the function F(δ) = −12.76470• δ + 6.44626 for better visibility.(b) Energy difference ∆E for the variance extrapolated data from (a).

28 L = 24 FIG. 4 .
FIG.4.SC correlation function P d (r) for Hg1201 at δ = 0.146 which is close to the experimental optimal doping, at the square lattice size L = 24, 28, 30, 32 and 36.Error bars are smaller than the symbol size.Inset: Size extrapolation of P d to the thermodynamic limit L → ∞.Because of relatively scattered data we employ the average of the two biggest lattice sizes (L = 32, 36) for the size extrapolation L → ∞.In fact, the value at the largest sizes is consistent with the systematic δ dependence observed near δ = 0.146 after the size extrapolation (not shown).The error bar at 1/L = 0 is estimated as the biggest difference to F ∞ SC from the given data.

16 FIG. 5 .
FIG. 5. SC correlation function P d (r) for Bi2212 at δ = 0.167 for the square lattice size L = 16, 24, 30, 36.Inset: Size extrapolation of Pd to the thermodynamic limit.The gray line shows the linear extrapolation to 1/L → 0 by using the values for L = 16, 24, 30, 36.Error bars are smaller than the symbol size.

FIG. 6 .
FIG.6.SC correlation function P d (r) for Bi2212 at δ = 0.167 in the cases of the single-and the explicit two-layer calculations.In the single-layer case a 16 × 16 square lattice was considered, which corresponds to a 16 × 16 × 2 lattice in the two-layer case.Error bars are smaller than the symbol size.

24 FIG. 7 .
FIG.7.SC correlation function P d (r) for Bi2201 at δ = 0.167 at the square lattice sizes L = 24, 30, 36.Inset: Size extrapolation of Pd to the thermodynamic limit.The gray line shows the linear extrapolation to 1/L → 0 by using the values for L = 24, 30, 36.Error bars are smaller than the symbol size.
FIG. 9. (a) FSC over scaling α or ξ (effective scaling of U/|t1| and Vi/U ) of the SC state on the L = 24 square lattice and δ = 0.167 hole doping.The inset is an enlarged plot around the peak 0.9 < α, ξ < 1.5.For further details see the main text.(b) Variance extrapolated energy difference ∆E between the SC and C4S8 as a function of ξ in the case H(α = 1.2, ξ).

cFIG. 10 FIG. 11 .
FIG. 10.F ∞SC as a function of U/|t1| for the four cuprate compounds at δ = 0.167 plotted from the list in Tables III and IV.

FIG. 15 .
FIG. 15.Size dependence of the charge and spin stripe state after the optimization of variational parameters for different lattice sizes and fillings.For C3S3 (red crosses) δ = 0.207 and L = 12, 24 and for C4S8 (blue dot) δ = 0.125 and L = 16, 24.(a) 10 3 × Ss(qSDW)/N vs 1/L.(b) 10 3 × Sc(qCDW)/N vs. 1/L 2 .qSDW and qCDW are the wave numbers at the peak of Ss and Sc, respectively, for the spin density wave (SDW) and charge density wave (CDW) modulations, respectively.The dashed lines are only to guide the eyes.

Appendix E :
Comparison to Hubbard model studies

Appendix G:
Analysis of δ dependence of energy In Fig. 18 the total energy per site E tot and the onsite Coulomb part E U (see panel (a), (c), respectively) are shown for doped CaCuO 2 .Each energy contribution is subtracted by a linear function F (δ) = b 0 + b 1 δ for better visibility, where b 0 and b 1 are listed in Table V (see gray lines in (a), (c)) and are shown in (b), (d).The subtracted energies are fitted by a quadratic function P(δ) = c 0 +c 1 δ +c 2 δ 2 to examine the curvature.Explicit values of the parameters are given in Table

FIG. 18 .
FIG.18.Energy per site as a function of δ in the SC state on a L = 24 square lattice.Total energy E total /N (a), and E total /N − F (δ) (b) as well as on-site Coulomb energy EU /N (c), and EU /N − F (δ) (d) are plotted.Here, a δ-linear function F (δ) defined below is subtracted in (b) and (d) for better visibility.Note that F (δ) is different depending on the type of the energy as seen in TableV.The gray line in the left column indicates F (δ) = b0 + b1δ, which is subtracted in (b) and (d).Right column: Energies after subtraction of F (δ) are fitted via a quadratic function P(δ) = c0 + c1δ + c2δ 2 .Fitting parameters are given in TableV.

FIG. 20 .
FIG. 20.Variance extrapolated energies of Bi2201 for various ground-state candidates (SC, charge-spin stripe (C3S3 and C4S8), and AFM state) as a function of hole doping δ on a L = 24 square lattice.(a) Total energies per site subtracted by F(δ).All energies are subtracted by the function F(δ) = −15.38797• δ + 7.79753.(b) Energy difference ∆E for the variance extrapolated data from (a).

TABLE I .
[12]nitio single-band effective Hamiltonian for CaCuO2, Hg1201, Bi2212, and Bi2201 taken from Ref.[12].U is the on-site interaction.The nth-neighbor hopping amplitude and Coulomb interaction are denoted as tn and Vn, respectively.Interlayer hoppings and interactions are neglected here.All values are given in eV.

TABLE II .
Size extrapolated SC correlation function P ∞ d and order parameter FSC for doped CaCuO2 for several choices of doping δ.The fitting error is of the order of ∼ 10 −4 for P ∞

TABLE III .
Comparison of long-ranged SC correlation P ∞ d and the order parameter F ∞ SC with U/|t1| as well as comparison between |t1|F ∞ SC and T opt c is intermediate and the uncertainty range is between 4.0 and 4.7 eV for U by considering that d z Oap may range from 2.25 Å to 2.45 Å.On the other hand, the value U ∼ 4.4 eV for Bi2201 in Table I is the upper bound and the uncertainty ranges from U ∼ 4.4 to 3.5 eV by considering that d z Oap may range from d z Oap = 2.6 to 2.45 Å.We first present in Secs.III C 1 and III C 2 the results obtained from the parameters shown in Table Table IV summarizes the size-extrapolated P ∞ SC for the two Bi compounds when we use the Hamiltonian parameters listed in Table

TABLE IV .
Comparison of long-ranged SC correlation P ∞ d and the order parameter F ∞ SC with U/|t1| for Bi2212 and Bi2201 at δ = 0.167, which is chosen in accordance with the experimental optimal values.The range of values represents the uncertainty range arising from the uncertainty of the apical oxygen position.This estimate helps the inference of the correct apical oxygen position (see text).

TABLE VII .
Ab initio single-band effective Hamiltonian for Bi2201 when the apical oxygen is shifted.The energy unit is eV.

TABLE VIII .
[12]nitio effective Hamiltonian for CaCuO2, Hg1201, Bi2212, and Bi2201.Hoppings and screened Coulomb interactions derived for the single-band effective Hamiltonians are taken from[12].The nth neighbor hopping amplitude and Coulomb interaction are denoted as tn and Vn.Interlayer hoppings and interactions are neglected here.All values are given in eV., that in the thermodynamic limit, this stripe state will fall into a paramagnetic state.S3.AF state for CaCuO 2 away from half filling indicate