Universal Chemical Formula Dependence of $Ab$ $Initio$ Low-Energy Effective Hamiltonian in Single-Layer Carrier Doped Cuprate Superconductors -- Study by Hierarchical Dependence Extraction Algorithm

We explore the possibility to control the superconducting (SC) transition temperature at optimal hole doping $T_{c}^{\rm opt}$ in cuprates by tuning the chemical formula (CF). $T_{c}^{\rm opt}$ can be theoretically predicted from the parameters of the \textit{ab initio} low-energy effective Hamiltonian (LEH) with one antibonding (AB) Cu$3d_{x^2-y^2}$/O$2p_{\sigma}$ orbital per Cu atom in the CuO$_2$ plane, notably the nearest neighbor hopping amplitude $|t_1|$ and the ratio $u=U/|t_1|$, where $U$ is the onsite effective Coulomb repulsion. However, the CF dependence of $|t_1|$ and $u$ is a highly nontrivial question. In this paper, we propose the universal dependence of $|t_1|$ and $u$ on the CF and structural features in hole doped cuprates with a single CuO$_2$ layer sandwiched between block layers. To do so, we perform extensive \textit{ab initio} calculations of $|t_1|$ and $u$ and analyze the results by employing a machine learning method called Hierarchical Dependence Extraction (HDE). The main results are the following: (a) $|t_1|$ has a main-order dependence on the radii $R_{\rm X}$ and $R_{\rm A}$ of the apical anion X and cation A in the block layer. ($|t_1|$ increases when $R_{\rm X}$ or $R_{\rm A}$ decreases.) (b) $u$ has a main-order dependence on the negative ionic charge $Z_{\rm X}$ of X and the hole doping $\delta$ of the AB orbital. ($u$ decreases when $|Z_{\rm X}|$ increases or $\delta$ increases.) We elucidate and discuss the microscopic mechanism of (a,b). We demonstrate the predictive power of the HDE by showing the consistency between (a,b) and results from previous works. The present results provide a basis for optimizing SC properties in cuprates and possibly akin materials. Also, the HDE method offers a general platform to identify dependencies between physical quantities.


I. INTRODUCTION
One of the grand challenges in condensed matter physics is the design of superconductors (SCs) with high transition temperature T c .The diverse distribution of T T opt c ≃ 0.16|t 1 |F SC (1) was proposed [14], in which the dimensionless SC order parameter F SC mainly depends on u.The u dependence of F SC is reminded in Appendix A and is summarized below.F SC is zero for u ≲ 6.5, and increases sharply with increasing u ≃ 6.5 − 8.0 (weak-coupling regime), reaching a maximum at u opt ≃ 8.0 − 8.5 (optimal regime); then, F SC decreases with increasing u ≳ 9.0 (strong-coupling regime).For a given material, |t 1 | and u can be calculated by using the multiscale ab initio scheme for correlated electrons (MACE) [15][16][17][18], which allowed to establish Eq. (1).Thus, a key point for materials design of higher-T opt c cuprates is to elucidate the universal chemical formula (CF) dependence of |t 1 | and u.In previous works on cuprates such as Hg1201, Bi 2 Sr 2 CuO 6 (Bi2201), Bi 2 Sr 2 CaCu 2 O 8 (Bi2212) and CaCuO 2 [19] as well as Hg1223 [20], the nontrivial dependence of |t 1 | and u on the interatomic distances and the CF has been partly clarified.However, the more general CF dependence of |t 1 | and u is required to obtain a thorough understanding of the CF dependence of T opt c .
The goal of this paper is twofold.First, we propose a machine learning procedure that is tailor-made to extract the nonlinear dependencies of a given quantity y on other quantities x i from the main-order to the higher-order.This procedure is denoted as Hierarchical Dependence Extraction (HDE).Second, we propose the universal CF dependence of |t 1 | and u in single-layer cuprates, by performing explicit ab initio calculations of the AB LEH for a training set that is representative of single-layer cuprates (including copper oxides, oxychlorides and oxyfluorides), and applying the HDE to analyze the results and construct expressions of |t 1 | and u.We generalize the existing MACE procedure to obtain the crystal parameters as a function of the chemical variables (the radii and charges of the cations and anions in the block layer), and in fine the AB LEH as a function of the chemical variables.The combination of the generalized MACE (gMACE) procedure with the analysis of the results by the HDE is denoted as gMACE+HDE.We demonstrate the predictive power of the HDE by showing the consistency between the universal CF dependencies of |t 1 | and u obtained by employing the gMACE+HDE and previous results on Hg1201, Bi2201, Bi2212, CaCuO 2 and Hg1223 [19,20].This paper is organized as follows.Section II gives an overview of the main-order dependence (MOD) of |t 1 | and u on the chemical variables.Section III describes the HDE and the gMACE methodologies employed in this paper, and how the HDE is applied to analyze the results of the gMACE calculation in the gMACE+HDE.Section IV details the results on the MOD of |t 1 | and u on the chemical variables, and proposes the microscopic mechanism underlying to this MOD.Section V discusses the results from the perspective of Eq. ( 1), and proposes guidelines to optimize the value of T opt c in future design of single-layer cuprates for which the gMACE calculation is performed.Section VI is the conclusion.Appendix A reminds the u dependence of F SC from Ref. [14].Appendix B details the choice of the training set of single-layer cuprates.Appendix C gives details on the HDE.Appendix D gives details on the gMACE procedure and the values of the intermediary quantities obtained in the gMACE calculation.Appendix E gives the analysis of the competition between variables in the MOD.Appendix F gives details on the hole doping dependence of the screening for each compound in the training set, and possible implications on the SC properties.Appendix G gives complements on the density of states near the Fermi level in hole-doped oxychlorides.(I) |t 1 | mainly depends on the crystal ionic radii R X and R A of the apical anion X and cation A in the block layer that separates two CuO 2 layers.(See Fig. 1 for an illustration of X and A.) The MOD up to the second order (MOD2) is 1   |t 1 | MOD2 = 0.534 − 0.000327 R 9.99 X + 7.99R 7.76   A .
Qualitatively, |t 1 | increases when R X or R A decreases (see Fig. 1).The microscopic mechanism is summarized as follows: Reducing R X and/or R A reduces the chemical pressure that pushes atoms apart from each other inside the crystal.This reduces the cell parameter a and thus the distance between Cu atoms in the CuO 2 plane, which increases the overlap between AB orbitals located on neighboring Cu sites, and thus |t 1 |.This result is consistent with that in [20], in which |t 1 | increases when applying physical uniaxial pressure P a along a direction.(The application of P a reduces a.) 1 In the MOD2s given in Eqs. ( 2), ( 3), ( 4) and ( 5 4)] and the screening ratio R = U/v [Eq.(5)].R A and Z A (R X and Z X ) are the crystal ionic radius and ionic charge of the A cation (apical X anion).n AB is the average number of electrons per AB orbital.
(II) u mainly depends on the negative ionic charge Z X of the apical anion and the average number of electrons n AB in the AB orbital.(We have n AB = 1 − δ , where δ is the hole doping.)The MOD2 is Qualitatively, u increases when |Z X | decreases or n AB increases (see Fig. 1).To understand the origin of the MOD of u in Eq. (3), we decompose u = vR/|t 1 |, where v is the onsite bare Coulomb interaction and R = U/v is the screening ratio, and we examine the MODs of v and R in (III) and (IV) below.(As shown below, the MOD of u is mainly determined by the MOD of R.) (III) v mainly depends on R A and the positive ionic charge Z A of the cation.The MOD2 is v MOD2 = 18.874 − 2.787R 0.92 A 1 + 87.90Z −9.12A .(4) Qualitatively, v increases when R A decreases or Z A increases (see Fig. 1).The microscopic mechanism is summarized as follows: Reducing R A or increasing Z A modifies the crystal electric field (namely, the Madelung potential created by cations and anions in the crystal) felt by the Cu3d x 2 −y 2 and O2p σ electrons.This stabilizes the in-plane O2p σ orbitals with respect to the Cu3d x 2 −y 2 orbitals, which increases the Cu3d x 2 −y 2 /O2p σ charge transfer energy ∆E xp .(Details are given later in Sec.IV.)This reduces the Cu3d x 2 −y 2 /O2p σ hybridization, which increases the Cu3d x 2 −y 2 atomic character of the AB orbital.This increases the localization of the AB orbital and thus v.This result is consistent with [20].
(IV) The CF dependence of R is more complex than that of |t 1 | and v, but we identify a rough MOD2 of R on Z X and n AB , which is R MOD2 = 4.224 − 3.906 |Z X | 0.01 + 0.00078n −9.99 AB .(5) Qualitatively, R increases when (i) |Z X | decreases or (ii) n AB increases (see Fig. 1).The microscopic mechanism of (i) and (ii) is briefly summarized below.(Details are given later in Sec.IV.) (i) Decreasing |Z X | reduces the negative charge of the apical anion.This reduces the negative Madelung potential (MP − ) created by the apical anion and felt by the electrons in the nearby CuO 2 plane.This reduces the energy of the electrons in the CuO 2 plane, and also reduces the Fermi energy.As a consequence, the empty states become higher in energy relative to the Fermi level.This reduces the screening from empty states, and thus, increases R.
(ii) The decrease in R with decreasing n AB (increasing δ ) is consistent with [19], in which the increase in δ causes the rapid decrease in R and thus u.(In [19], calculations were made at fixed |Z X | = 2 and varying n AB = 1.0, 0.9 and 0.8, which corresponds to δ = 0.0, 0.1 and 0.2.)The rapid decrease in u eventually suppresses F SC [14] so that the system ends up in the metallic state, in agreement with the experimental ground state in the overdoped region.
Remarkably, the MOD2 of R [Eq.(5)] is very similar to the MOD2 of u [Eq.( 3)].Thus, the MOD2 of u = vR/|t 1 | is dominated by the MOD2 of R. Consistently, in the ab initio result, the diverse distribution of u ≃ 7 − 10.5 originates mainly (albeit not exclusively) from the diverse distribution of R. Indeed, the relative variation between the minimum and maximum ab initio values is 27% for v ≃ 13.5 − 16.5 eV, 30% for |t 1 | ≃ 0.40 − 0.55 eV, 48% for R ≃ 0.22 − 0.34, and 50% for u ≃ 7 − 10.5.(The variation in u is reproduced by that in R.) Still, the relative variation in |t 1 | and v is nonnegligible compared to that in R, so that the CF dependencies of |t 1 | and v also contribute to the CF dependence of u beyond the MOD2 in Eq. (3).In Sec.IV, we will decompose u = vR/|t 1 |, and we will discuss in detail the CF dependence of |t 1 |, v and R.

A. Framework of HDE
Here, we summarize the HDE methodology to construct a descriptor for a given physical quantity y as a function of other given quantities in the variable space V = {x i , i = 1..N V }, where i is the variable index.(We make complete abstraction of the physical meaning of these variables.)We restrict the presentation to the essence of the procedure; complementary discussions and computational details are given in Appendix C, and the motivation of the HDE procedure is discussed in detail in Appendix C 1.
General problem -Starting from the variable space V, we define the candidate descriptor space The elements x(p) of C V are called "candidate descriptors", and are functions of the variables x i in V; their analytic expression depends on variational parameters which are encoded in the vector p. [The expression of x(p) and definition of variational parameters are given later.]The general problem consists in finding p opt such that x opt = x(p opt ) is the best candidate descriptor for y among the elements of C V , i.e.
where f is a fitness function that describes how well y is described by x(p opt ).
Fitness function -We choose the fitness function where ρ is the (for k 1 ̸ = 0), which is prominently used in the HDE procedure as discussed later.The value of f is insensitive to the values of k 0 and k 1 in Eq. ( 9), and has an intuitive interpretation irrespective of the values of k 0 and k 1 .(See Appendix C 2 for details.) Expression of the candidate descriptor -The descriptor is constructed iteratively by adding factors that contain dependencies of y on x i from the lowest order to the highest order.
The expressions of the candidate descriptors that are considered in this paper at generation g = 1 and g ≥ 2 are, respectively, where is the best candidate descriptor at generation g − 1, i g is one of the indices between 1 and N V , and the parametric operator ⋆ (ζ ,β ,α) (dubbed hereafter as the wildcard operator) is defined as The variational parameter vector at g is .N, we obtain p (N) = (p 1 , p 2 , ..., p N ) that corresponds to p in Eq. (7).At g, the values of variational parameters that maximize the fitness function are denoted as ).The wildcard operator in Eq. ( 12) is versatile and can represent any algebraic operation depending on the values of (ζ , β , α), as discussed in detail in Appendix C 3. In the practical procedure, (ζ g , β g , α g ) are optimized together with i g to maximize f [y, x (g) ], so that the character of the wildcard operator is automatically adjusted to describe y as accurately as possible.Also, note that the wildcard operator can mimick the identity operator for specific values of (ζ , β , α) (see Appendix C 3), so that x opt (g−1) is included in the candidate descriptor space at g.This guarantees that Finally, note that even though y has an affine dependence on x(p) when f [y, x(p)] = 1, x(p) has a nonlinear dependence on x i in the general case.Thus, nonlinearity in the dependence of y on x i may be described by the above formalism.
Procedure to construct the descriptor for y as a function of x i -To obtain x opt (N) , we employ the following procedure, which is denoted as 10) that maximizes f[y, x (1) ], and we obtain x opt (1) .Then, we increment g, we determine p 11) that maximizes f[y, x (g) ], and we obtain x opt (g) .We iterate up to g = N. (The computational details of the optimization are given in Appendix C 4.) We obtain From this descriptor, we deduce the estimated expression of y as a function of x i , as in which the coefficients k 0 and k 1 are calculated by an affine regression of y on x opt (N) .
Completeness of the dependence of y on x i -The procedure P[y, V] allows to probe the completeness of the dependence of y on x i .When we perform P[y, V], we assume the following conjecture: D[y, V] The dependence of y is entirely contained in V = {x i }. [This implies that there exists p such that f [y, x(p)] = 1.]After performing P[y, V], the validity of D[y, V] can be checked by examining the value of in which If not, there are two possibilities: either D[y, V] is incorrect, or D[y, V] is correct but the procedure P[y, V] is insufficient to capture the whole dependence of y on the variables x i in V.In the scope of this paper, we assume D[y, V] is incorrect.
In case D[y, V] is incorrect, it is possible to improve the descriptor by replacing V with a superset of V, as mentioned in the next paragraph.In case D[y, V] is correct, we obtain x opt (N) on which y has an affine dependence, and we obtain the expression y N of y as a function of x i in Eq. (15).[The correctness of D[y, V] is the guarantee that the affine interpolation in Eq. ( 15) is accurate.]Main-order dependence of y on x i -The expression of x opt (N) [Eq.( 14)] has a hierarchical structure which reveals the hierarchy in the dependencies of y on x i : The lower (higher) values of g correspond to the lower-order (higher-order) dependencies of y on x i , and the variable x i g contains the g thorder dependence of y.Note that, when incrementing g, we allow only up to one variable to be introduced in Eq. ( 14): This allows to select the variable which corresponds to the g th -order dependence.Also, note that several variables may be in close competition in the g th order dependence, as discussed later.
In this paper, we mainly discuss the MODs that are contained in g ≲ 2 − 3; the full list of dependencies for g from 1 to N is given in Sec.S1 of Supplemental Material (SM) [21].For g ≲ 2 − 3, we define the MOD of y on x i up to the g th order [MODg], as y MODg = y g , which is Eq. ( 15) with N = g.Because y MODg depends on no more than g different variables in V, it is possible to represent graphically y MODg as a function of x i for g ≤ 2 by using e.g. a color map, as done in Fig. 1.
The accuracy of the MODg is quantified in f (g) [y, V].Note that, in the general case, the MODg of y on x i is but a rough description of y: Typically, f (g) [y, V] ≃ 0.65 − 0.95 for g ≲ 3.
The accurate description requires to take into account the higher-order dependencies of y on x i beyond g ≲ 3 as well.Nonetheless, the MODg contains the principal mechanism of the dependence of y; also, for g ≤ 2, we have f (g) [y, V] ≃ 1 in some particular cases, as seen later.
Competition between the x i in the dependence of y -At the generation g, we determine the variable x i that corresponds to x i opt g , as mentioned before.We also examine the competition between variables as detailed below.We define the maximal fitness of the variable x i at the generation g = 1 and g ≥ 2 as, respectively, which is the maximal value of the fitness function that is obtained at g by enforcing x i opt g = x i .Also, we define the score of the variable x i at the generation g as We have s (g) [y, x i ] = 1 if x i corresponds to x opt i g that is found in the optimization, and s (g) [y, x i ] = 0 if x i has the lowest f(g) [y, x i ] among the variables in V. Complementary discussions are given in Appendix E. In the following, the analysis of the competition between variables is performed in Appendix E and briefly mentioned in the main text.
Practical HDE procedure -In practice, given y and V = {x i }, the HDE is denoted as HDE[y, V] and is employed as follows.We perform P[y, V].The output quantities are f ∞ [y, V] and p opt g ) for g from 1 to N.Then, we check the validity of D[y, V] by examining the value of f ∞ [y, V].If D[y, V] is incorrect, we may attempt to replace V by a superset of V and relaunch the procedure.
If D[y, V] is correct, we perform the below restricted procedure, denoted as rHDE[y, V].We attempt to simplify the dependence of y by eliminating the variables on which y has a higher-order dependence, in decreasing order.Namely, we take the optimized variable indices i (Each time we increment j, we remove the variable that corresponds to the highest-order dependence.)We check whether is the minimal subset of V that describes y entirely.

B. Framework of gMACE
Next, we summarize the ab initio MACE scheme and its generalization to the gMACE in the present paper.[See Fig. 2 for an illustration.]We restrict the presentation to the essence of the procedure; details are given in Appendix D.
Chemical formula dependence of crystal parameters -Obtaining the CF dependence of the AB LEH requires to extend the standard MACE scheme [15][16][17][18][19][20][22][23][24].The latter allows to calculate the AB LEH by starting from a given CF together with the crystal symmetry and crystal parameter (CP) values.(The CP values are usually taken from experiment.)To calculate the AB LEH as a function of the CF, the missing step is to calculate the CP as a function of the CF.We add this missing step to MACE by introducing the CF variables (the ionic radii and charges of cations and anions in the crystal) and calculating the CP as a function of the CF by performing the structural optimization, instead of relying on the experimental CP values.[Even though the experimental CP may be more accurate than the optimized CP, the experimental CP is not always available, and the structural optimization allows to obtain the systematic CF dependence of the CP.]We only assume the symmetry of the primitive cell during the structural optimization (see Appendix D 1 for details).
Calculation of the AB LEH -After obtaining the CP for a given CF by employing the structural optimization, the AB LEH is calculated by following the MACE procedure as employed in [20], which is summarized below in the successive steps (i-v).This procedure combines the generalized gradient approximation (GGA) [25] and the constrained random phase approximation (cRPA) [26,27], and is denoted as GGA+cRPA.
(i) Starting from the CF and the CP values, we first perform a DFT calculation.This allows to obtain the DFT electronic structure at the GGA level.
(ii) From the DFT electronic structure, we compute the maximally localized Wannier orbitals (MLWOs) [28,29] that span the medium-energy (M) space.The M space consists in the 17 bands with Cu3d-like, in-plane O2p-like and apical X2p-like character near the Fermi level (see Fig. 2).The band with the highest energy in the M space contains most of the AB character, and the bands outside the M space form the high-energy (H) space.The MLWO centered on the atom l and with m-like orbital character is denoted as (l, m), where the atom index l, l ′ takes the values Cu, O, O ′ , X (O and O' denote the two in-plane O atoms in the unit cell, and the atomic positions are given in Appendix D 1).The orbital index m, m ′ takes the values x, z, xy and yz for the Cu3d x 2 −y 2 , Cu3d 3z 2 −r 2 , Cu3d xy and Cu3d yz orbitals (the Cu3d zx and Cu3d yz orbitals are equivalent), p, p π and p z for the in-plane O2p σ , O2p π and O2p z orbitals, and p x , p z for the apical X2p x and X2p z orbitals (the X2p y and X2p x orbitals are equivalent).Note that O2p σ , O'2p σ , O2p π and O'2p π correspond respectively to O2p x , O'2p y , O2p y and O'2p x .
(iii) From the M space, we compute the AB MLWO by using a procedure that is detailed in Appendix D. The AB MLWO centered on the Cu atom located in the unit cell at R is denoted as w R .We also obtain the AB band that corresponds to the dispersion of the AB MLWO.Then, the other 16 bands in the M space are disentangled [30] from the AB band.(See Fig. 2 for an illustration of the AB band and disentangled M band dispersions.)(iv) From the AB MLWO, we compute in which Ω is the unit cell, R 1 = [100], and h GGA is the oneparticle part at the GGA level.We also compute the onsite bare Coulomb interaction as in which v(r, r ′ ) is the bare Coulomb interaction.(v) From the AB band, the disentangled M bands and the H bands, we compute the cRPA screening ratio R (and also u = U/|t 1 |) as follows.We compute the cRPA effective interaction W H (r, r ′ ) at zero frequency, whose expression is found in Eqs.(D1) and (D4) in Appendix D 1.We deduce the onsite effective Coulomb interaction U by replacing v(r, r ′ ) with W H (r, r ′ ) in Eq. (22).We deduce R = U /v and u = U/|t 1 |.

C. Training set of cuprates and gMACE+HDE procedure
Below, we define the training set of cuprates that is considered in this paper.(For each CF in the training set, we apply the gMACE procedure described in Sec.III B.) Then, we describe how the HDE is applied to analyze the results of the gMACE calculation and extract the CF dependence of the AB LEH parameters.
Training set -The training set of cuprates is defined below and illustrated in Fig. 3. (Detailed discussions on the choice of the training set are given in Appendix B.) The training set includes N tr = 36 CFs, including both experimentally confirmed and hypothetical SC compounds.The general CF is A'A 2 CuO 2 X 2 and the block layer consists in A'A 2 X 2 .For the undoped compound, we have A' = Hg, ∅, A = Ba, Sr, Ca, La, Y, Sc, and X = O, F, Cl.For the doped compound, we use the same procedure as in [19,20]: We use the virtual crystal approximation [31] to substitute part of the A or A' chemical element by the chemical element whose atomic number is that of A or A' minus one.We consider hole doping δ = 0.0, 0.1 and 0.2 (that is, up to 20% hole doping).This range includes the experimental range in which the SC state is observed.We have The ionic charges are related to δ and n AB as follows: Note that all CFs in the training set do not need to correspond to experimentally confirmed SC cuprates.For a given CF, the gMACE result can be used as part of the data that is analyzed by the HDE procedure to infer the systematic CF dependence of |t 1 | and u, by making complete abstraction of whether or not the corresponding crystal structure can be stabilized in experiment.gMACE+HDE procedure -We apply the HDE to express the AB LEH parameters as a function of the CF variables.In addition, to gain further insights on the underlying microscopic mechanism, we apply the HDE to express the intermediate quantities within the gMACE as functions of each other.Namely, we consider N s = 4 levels of variables whose definition is guided by the hierarchical structure of gMACE as illustrated in Fig. 2. At each step s, we define the variable space V s = {x s i , i = 1..N V s }, and we use the HDE to express x s i as a function of the variables in V s ′ (with s ′ < s).For each s, the variables in V s are chosen as follows.(These variables are illustrated in Fig. 2.) R A ′ is the crystal ionic radius of the second cation A' = Hg 1−x Au x in the block layer (R A ′ is set to 0 if A' = ∅).Values of the variables in V 1 that are considered in this paper are given in Appendix D 3. R A , R X and R A ′ are expressed in Å throughout this paper.DFT band structure (s = 3): We consider

Crystal parameters (s
These variables consist in characteristic energies within the M space; they are defined below, and their choice is further discussed in Appendix D 2. First, we include the absolute value of the onsite energies ε l m of all MLWOs (l, m) in the M space.(Note that ε l m < 0 because the onsite energy is below the Fermi level.)Second, we include the nonzero hopping amplitudes |t l,l ′ m,m ′ | between the MLWOs (l, m) and (l ′ , m ′ ) in the unit cell, where (l, m) = (Cu, x), (O, p).We do not take into account hoppings in the unit cell that do not involve (Cu, x) Here, we first detail (I).Then, we detail the items (1-5) in Fig. 5.
Second, the dependence of Third, we discuss details of the MOD2 of |t 1 | on V 1 , which was given in Eq. ( 2).The MOD2 is not sufficient to describe |t 1 | entirely, but the main-order dependence of Indeed, we have f and is represented in Fig. 6(d,e).Note that |t xp | dominates over ∆E xp in the MOD2 (see the score analysis in Appendix E), which is also visible in Fig. 6(e): The color map has a horizontal-like pattern.24), and this is consistent with previous works.In the case of Hg1223, we have 12 of [20].(The latter result was obtained by modifying a artificially without modifying other CPs.)The result in Eq. ( 24) is more general, because Eq. ( 24) is established for N tr = 36 CFs rather than one CF, and it accounts for the CF dependence of all CPs.

FIG. 6. Details on the dependence of |t
The MOD2 in Eq. ( 24) can be interpreted as follows.|t 1 | represents the hopping amplitude between neighboring AB FIG. 7. Details on the dependence of 25)], the MOD2 of ∆E xp on V 2 [Eq.( 28)], the MOD2 of d z X on V 1 [Eq.( 29)], the MOD2 of a on V 1 [Eq.( 27 and is represented in Fig. 7(b,g).The score analysis shows that the dominance of a in the dependence of The a dependence of |t xp | is consistent with results on Hg1223 [20].In particular, the exponent −3.75 in Eq. ( 25) is very close to that obtained for Hg1223, in which |t xp | ∝ a −3.86 (see [20], Fig. 12).This suggests |t xp | scales as a −3.75 universally and irrespective of the crystalline environment outside the CuO 2 plane.This is intuitive because the Cu3d x 2 −y 2 and O2p σ orbitals extend mainly in the CuO 2 plane as illustrated in Fig. 2.
In addition, |t 1 | increases when a decreases according to Eqs. ( 25) and (24).This is consistent with [20], in which the pressure-induced decrease in a is the main cause of the pressure-induced increase in |t 1 |.For completeness, we per- and is represented in Fig. 7(f,k).Qualitatively, |t 1 | MOD1 increases when a decreases, which is consistent with [20] and also with Eqs. ( 25) and (24).For completeness, note that there is a quantitative difference between Eq. ( 26) and [20]: In the latter, we have |t 1 | ∝ a −2.88 (see [20], Fig. 12).The difference may be explained as follows.

973, and
which is represented in Fig. 7(e,j).The score analysis confirms that R A and R X correspond respectively to x i opt 1 and x i opt 2 (see Appendix E).The MOD2 in Eq. ( 27) is interpreted as follows.Qualitatively, a MOD2 increases when R A or R X increases, which is consistent with the hard sphere picture illustrated in Fig. 8: The interatomic distances and thus the cell parameter a increase when the ionic radii are larger.Thus, the hard-sphere FIG. 8. Illustration of the hard sphere picture for the crystal of the single-layer cuprate.In this picture, cations and apical anions are assumed to be rigid spheres that touch each other, and we have picture reproduces the qualitative dependence of a on the ionic radii.Furthermore, the MOD2 in Eq. ( 27) is qualitatively consistent with experiment, in which a increases with increasing R A .For instance, for X = Cl, we have a = 3.87 Å for A = Ca [32,33], a = 3.97 − 3.98 Å for A = Sr [34,35], and a = 4.10 Å for A = Ba [10].(These values are in correct agreement with the values a = 3.88 Å for A = Ca, a = 4.00 Å for A = Sr and a = 4.15 Å for A = Ba obtained in this paper by employing the structural optimization.)And, we have R Ca = 1.14 Å, R Sr = 1.32 Å and R Ba = 1.49Å [36].(The values of R A are given in Appendix D 3, and the R A dependence of experimental a is also emphasized in [10].)Also, for X = O, we have a ≃ 3.78 Å in La 2 CuO 4 [37] and a ≃ 3.88 Å in HgBa 2 CuO 4 [4].(These values are in correct agreement with the values a ≃ 3.82 in La 2 CuO 4 and a ≃ 3.92 Å in HgBa 2 CuO 4 obtained in this paper by employing the structural optimization.)And, we have R La = 1.17 Å and R Ba = 1.49Å [36].
(4) Dependence of ∆E xp on V 2 -∆E xp is not determined entirely by V 2 , but the MOD2 reveals the main-order mechanism that controls the value of is not completely correct.As for the MOD2, we have (28) and is represented in Fig. 7(c,h).Even though d z X corresponds to x i opt 1 , the dependence of ∆E xp on d z X is equally important to that on a (see the score analysis in Appendix E).Consistently, in Fig. 7(c), the color map has a diagonal-like pattern.
Qualitatively, ∆E xp increases when d z X decreases or a decreases in Eq. (28).The increase in ∆E xp with decreasing d z X can be understood as follows.When d z X decreases, the distance between the apical X and the CuO 2 plane decreases, so that the MP − created by the apical X anion and felt by the Cu and in-plane O is stronger.This increases the energy of both Cu3d x 2 −y 2 and O2p σ electrons.The MP − felt by Cu is the strongest, because the Cu is closer to the apical X compared to the in-plane O. [See Fig. 9(b) for an illustration.]Thus, the energy of Cu3d x 2 −y 2 electrons increases more than that of in-plane O2p σ electrons.As a consequence, ∆E xp increases.
The increase in ∆E xp with decreasing a is consistent with [20], and the mechanism is reminded here.When a decreases, the distance between the in-plane O and Cu is reduced, so that the MP − created by the in-plane O anions and felt by the Cu is stronger.[See Fig. 9(c) for an illustration.]This increases the energy of Cu3d x 2 −y 2 electrons with respect to that of in-plane O2p σ electrons, which increases ∆E xp .
979, and is represented in Fig. 7(d,i).The score analysis confirms that R A and R X correspond respectively to x i opt 1 and Qualitatively, d z XMOD2 increases when R A or R X increases in Eq. ( 29).This can also be understood by considering the hard sphere picture (see the right panel in Fig. 8).In the ab initio result, we always have

the values of d z
A and d z X are given in Appendix D 3), so that the apical X is farther from the CuO 2 plane compared to the A cation.In the hard sphere picture, increasing the ionic radius R A of the A cation pushes the apical X even farther from the CuO 2 plane, which increases d z X .The same mechanism occurs when R X increases.
Dependence of ∆E xp on V 1 -On (4,5), for completeness, we discuss the dependence of

982, and the MOD3 is
The score analysis confirms that R A and R X correspond respectively to x i opt 1 and x i opt 2 , and n AB corresponds to x i opt 3 but is slightly in competition with Z A (see Appendix E).
Qualitatively, ∆E xpMOD3 increases when R A decreases, R X decreases or n AB increases in Eq. (30).This is consistent with Eqs. ( 28), ( 29) and (27).Also, ∆E xp decreases when n AB decreases (the hole doping δ increases).This is consistent with [19] ] between A and in-plane O, which is smaller than the distance d A,Cu = a 2 /2 + (d z A ) 2 between A and Cu.Thus, we have

B. Chemical formula dependence of v
Here, we first detail (III).Then, we detail the item (6) in Fig. 5. [The items (2-5) have already been discussed in the previous section].4).The score analysis shows that R A and Z A correspond to x i opt According to both Eq. ( 31) and Eq. ( 32), v increases when |t xp | decreases or ∆E xp increases, which is consistent with [20].
C. Chemical formula dependence of R Here, we first detail (IV).Then, we discuss the items (7-10) in Fig. 5. (5)] reveals the main-order mechanism of the dependence of R. Namely, R has a very rough MOD2 on |Z X | and n AB , which is consistent with the below discussion.The score analysis shows that Z X and n AB correspond respectively to x i opt Qualitatively, R increases when (i) |Z X | decreases or (ii) n AB increases in Eq. ( 5) (see also Fig. 1).The microscopic mechanism of (i) and (ii) is detailed below.change substantially.However, because the Fermi level is reduced as discussed above, the empty states become higher in energy relative to the Fermi level (so that ε e increases).This reduces the screening from empty states, and thus, increases R.
(ii) The decrease in R with decreasing n AB (increasing δ ) is discussed in [19], and the microscopic mechanism is reminded here.When δ increases, the hole doping of O sites increases.This reduces the negative charge of in-plane O anions.This reduces the MP − created by the in-plane O and felt by the nearby Cu. [See Fig. 9(c) for an illustration.]This reduces the energy of Cu3d electrons, and also reduces the Fermi energy.(Indeed, the AB band at the Fermi level has Cu3d x 2 −y 2 character.)On the other hand, the O2p electrons are less affected.However, because the Fermi level is reduced, the occupied O2p states become closer to the Fermi level (so that |ε o | decreases).This increases the screening from occupied states, and thus, reduces R.
For completeness, the hole doping dependence of R is discussed in detail in Appendix F, which is summarized here.Even though R increases when |Z X | decreases as discussed above, the decrease in |Z X | also accelerates the decrease in R with decreasing n AB , which is not captured by Eq. ( 5).This is why the three color points with the lowest R in Fig. 11(b) deviate from the linear interpolation, which is a major cause of the relatively low value of f (2) [R, V 1 ] = 0.642.The MOD2 in Eq. ( 5) may be combined with the results in Appendix F to obtain a more accurate picture of the dependence of R on |Z X | and n AB .
Thus, from A' = ∅ to A' = Hg 1−δ Au δ , |ε X p x | is reduced by 1.62 eV for a given value of |ε X p z |.On the other hand, for A' = ∅, the value of k 0 is universal irrespective of A and X.This suggests that the decrease in |ε X p x | from A' = ∅ to A' = Hg 1−δ Au δ does not depend on A or X, but rather on the presence of the A' atom and the subsequent change in crystal structure and crystal electric field.A possible explanation is the following: For A' = Hg 1−δ Au δ , there is a A' atom close to the apical X (see Fig. 2), and the apical X2p orbital overlaps with the A'5d orbitals.This may cause the apical Xp x MLWO to catch antibonding X2p/A'5d character 3 , which may destabilize the apical Xp x MLWO (i.e. increase its onsite energy and thus reduce |ε X p x |).This is consistent with the isosurface of the apical Xp x MLWO in Fig. 2 for A' = Hg 1−δ Au δ : We see that the apical Xp x MLWO has a slight A'yz/zx character near the A' atom.
Second, |ε X p x | increases when R X decreases.This is consistent with [20]: When a decreases, the occupied bands in the M space (including apical X bands) become farther from the Fermi level [see [20], Fig. 3(h-k)].This is because the MP − created by the in-plane O and felt by the Cu increases [as illustrated in Fig. 9(c)], so that the energy of the Cu3d x 2 −y 2 increases (and this shifts the Fermi level upwards), whereas the apical X orbitals are less affected.And, the decrease in R X contributes to decrease a, as discussed previously.35) is understood as follows.First, |ε o | decreases when n AB decreases (i.e. the hole doping increases).This is consistent with [19], and the cause is interpreted as a rigid shift of the Fermi level upon hole doping.Because the partially filled AB band is the only band at the Fermi level, increasing the hole doping reduces the number of electrons in the AB band, which shifts the Fermi level downwards.As a result, the energy of occupied bands relative to the Fermi level increases.These include the highest occupied band outside the AB band, whose energy is ε o .Thus, ε o < 0 increases, so that |ε o | decreases.
Second, |ε o | increases when R A decreases.The microscopic mechanism is discussed in detail in Appendix G.
(10) Dependence of ε e on V 1 -The MOD of ε e on |Z X | and n AB in Eq. ( 36) is understood as follows.First, ε e increases when |Z X | decreases.The mechanism was summarized in Sec.II [see (IV), item (i)], and is detailed here.Apical X anions with the negative charge Z X emit a MP − that increases the energy of surrounding electrons, in particular those in the M bands, because the Cu and in-plane O atoms are in the vicinity of apical X. Reducing |Z X | reduces the MP − from apical X felt by the Cu and in-plane O. [See Fig. 9(b) for an illustration.]This reduces the energy of M bands.The Fermi energy is also reduced, because it is determined by the partially filled AB band which is in the M space.Thus, the empty bands become farther from the Fermi level.This is consistent with results on Hg1223 in [20] [see e.g.Appendix E1].
Second, ε e increases when n AB decreases.This is because (i) the Fermi level is shifted downwards when n AB decreases because of the hole doping of the AB band [as discussed in (9)]: As a consequence, the empty bands become farther from the Fermi level, and thus ε e increases.In addition, (ii) when n AB decreases, the M bands are stabilized, which further shifts the M bands and thus the Fermi level downwards.This is because the hole doping of in-plane O increases when n AB decreases as mentioned before: This reduces the negative charge of the in-plane O ions, and thus the Madelung potential created by the in-plane O ions.[See Fig. 9(c) for an illustration.]This stabilizes the M bands, which increases the energy gap between M bands and empty bands.

V. DISCUSSION
Here, we discuss the universality and accuracy of the obtained expressions of |t 1 |, v and R. Also, we propose prescriptions to optimize T opt c in future design of SC cuprates and akin materials.

A. Universality and accuracy of the expressions of AB LEH parameters
The CF dependencies of |t 1 |, v and R obtained in this paper offer reliable guidelines for design of SC cuprates and akin compounds, if we assume that (A) the training set considered in this paper is representative of the diversity in single-layer cuprates, and (B) the CF dependence of the AB LEH parameters |t 1 | and u is correctly captured (at least qualitatively) by the GGA+cRPA version of MACE employed in this paper.
(A) is supported in Appendix B; below, we support (B).The detailed discussion on (B) is necessary, because the GGA+cRPA is the simplest level of the MACE scheme; more sophisticated versions of MACE have been employed in previous works, such as the constrained GW (cGW ) supplemented with self-interaction correction (SIC) at the cGW -SIC level [17] and level renormalization feedback (LRFB) at the cGW -SIC+LRFB level [18].Eq. ( 1) was determined by solving AB LEHs at the cGW -SIC+LRFB level [19].
The GGA+cRPA is expected to capture correctly the qualitative materials dependence of |t 1 | and u [19,20] while avoiding the complexity of the cGW -SIC+LRFB calculation.For instance, u is smaller in Bi2201 compared to Bi2212 at the cGW -SIC+LRFB level, and this result is reproduced qualitatively by the GGA+cRPA [19].In addition, in Hg1223, the qualitative pressure dependence of |t 1 | and u is captured by the GGA+cRPA [20].
Still, it should be noted that the accurate prediction of the materials dependence of T opt c cannot be done by considering the materials dependent |t 1 | and u at the GGA+cRPA level (the values in Fig. 4).Indeed, the GGA+cRPA has a limitation at the quantitative level.For instance, in Bi2201 and Bi2212, u is underestimated at the GGA+cRPA level, and the difference between the values of u (|t 1 |) at the GGA+cRPA and cGW -SIC+LRFB levels is around 10% (5%) [19].The uncertainty on u may cause a significant uncertainty on F SC and thus T opt c , because T opt c strongly depends on u via F SC in Eq. ( 1), especially when u is located in the weak-coupling region u ≃ 6.5 − 8.0 [14].Thus, the quantitative improvement from the GGA+cRPA level to the cGW -SIC+LRFB level is required to tackle the accurate prediction of T opt c ≃ 0.16|t 1 |F SC [14] from the values of |t 1 | and u.
Nonetheless, we restrict to the simplest GGA+cRPA level in the present paper, because (i) the cGW -SIC+LRFB calculation is complex and computationally expensive, and (ii) the GGA+cRPA is expected to capture correctly the materials dependence of the AB LEH at least qualitatively [19,20] as discussed above.Quantitative prediction of the CF dependence of T opt c by using Eq. ( 1) requires the CF dependence of the cGW -SIC+LRFB result, which is left for future studies.
Also, note that some of the results obtained in this paper are independent of the restriction to the GGA+cRPA level, and are expected to remain valid for the AB LEH at the cGW -SIC+LRFB level.This is the case of the items (3,5,6,7) in Fig. 5.For instance, the dependencies of R and v on V 3 in Eqs. ( 33) and (31) make complete abstraction of the level of the electronic structure from which the AB LEH is calculated (GGA in the case of GGA+cRPA, or GW +LRFB [18] in the case of cGW -SIC+LRFB).Thus, the dependencies of R and v on V 3 in Eqs. ( 33) and ( 31) are expected to be rather universal.(The detailed dependencies of and ε e at the GW +LRFB level on the CF are left for future studies.)However, on |t 1 |, the dependence on V 3 in Eq. ( 24) does not take into account the removal of exchange-correlation double counting [16] that is done at the cGW -SIC+LRFB level.This corrects |t 1 | by a term which is materials dependent [19,20].(The detailed materials dependence of this term is left for future studies.)should be satisfied.To satisfy Eq. ( 39), the values of |t 1 |, v and R may be tuned rather independently if we consider their distinct dependencies on the CF.For instance, Z A appears in the MOD of v [Eq.( 4)] but not in that of |t 1 | [Eq.( 2)], so that tuning Z A allows to tune v without affecting |t 1 | substantially.
Note that the dependence of |t 1 | on the CF in Eq. ( 2) predicts an upper bound T opt c,max ≃ 140 K for T opt c in single-layer cuprates at ambient pressure.Indeed, according to Eq. ( 2), there is a maximal value of |t 1 |, which is |t 1 | max ≃ 0.53 eV (|t 1 | max ≃ 0.58 eV if we consider g = 15 instead of g = 2).And, in the universal u dependence of F SC [14], the maximal value of F SC at u = u opt is F SC,max ≃ 0.13.Thus, for The upper bound |t 1 | max corresponds to R X = R A = 0, which cannot be reached in experiment.However, reducing R X or R A may allow to increase |t 1 | and make |t 1 | as close as possible to |t 1 | max .Reducing R X or R A can be done by e.g.replacing A or X by an isovalent atom that has a smaller ionic radius.This contributes to reduce the interatomic distances between atoms in the crystal and thus a according to Eq. ( 27), by mimicking the effects of physical P a at the chemical level.
Under pressure, |t 1 | may be further increased due to the application of physical P a , and may reach or exceed the upper bound |t 1 | max at ambient pressure. 4Note that the P a -induced increase in |t 1 | eventually causes a rapid decrease in u and thus F SC by making u fall into the weak-coupling region [20]; this may be countered by tuning Z A and Z X to increase v and R, so that the criterion in Eq. ( 39) is still satisfied at high P a .
Finally, for completeness, we discuss the particular case of X = F, Cl in Appendix F. Indeed, for X = F, Cl, the decrease in u with decreasing n AB is sharper compared to X = O, which is an effect beyond the MOD2 in Eq. ( 3).The origin of this sharper decrease is the sharper decrease in R for X = F, Cl compared to X = O, as mentioned before and detailed in Appendix F. Possible implications on SC properties are discussed in Appendix F. The qualitative insights obtained in this paper may also be useful for design of other SC materials whose crystal structure is similar to that of cuprates, namely, other strongly correlated electron materials whose low-energy physics may be described by the AB LEH on the two-dimensional square lattice.These include nickelates [38] and the recently proposed Ag-and Pd-based compounds [24,39].
More generally, the gMACE+HDE procedure employed in this paper may be used for design of strongly correlated electron materials including those which do not have SC properties.Even more generally, the HDE procedure offers a general platform to extract dependencies between any physical quantities y and x i , regardless their physical meaning.These quantities may be either calculated within a theoretical framework or measured in an experiment.multi-layer cuprates [14,19,20].In particular, the values of u correspond to the range u ≃ 6.5 − 10.5 in which F SC is nonzero (see [14], Fig. 10), and thus is nonzero according to Eq. ( 1).This suggests the training set is a good platform to study the microscopic mechanism of the materials dependence of |t 1 | and u and thus T  The essence of the HDE was presented in Sec.III A. Here, as a complement, we discuss in more detail the motivation of the HDE procedure.
The complete elucidation of the dependence of y on V = {x i } requires several conditions: (a) Probe the completeness of the dependence of y on V; (b) Extract the hierarchy in the dependencies of y on x i ; (c) Capture the nonlinear dependence of y on x i ; (d) Obtain an explicit expression of y as a function of x i .In addition, it is desirable to (e) keep the procedure as simple as possible, and (f) preserve the sparsity in the approached expression of y (namely, the approached expression of y should depend on as few x i as possible).
On (a), probing the completeness of the dependence of y on V allows to clarify whether or not it is possible to construct a perfect descriptor for y as a function of x i .This provides useful information prior to (b,c).
On (b), clarifying the hierarchy in the dependencies of y on x i is necessary to pinpoint the MOD of y, namely, the variables x i on which y depends the most.This allows to construct an approximation of y as a function of these x i by combining (b) and (d), which consists in y MODg which is defined and discussed in the main text.The MODg of y contains the principal microscopic mechanism of the dependence of y, and thus, provides useful clues in the context of materials design.
On (c), y has a nonlinear dependence on x i in the general case.For instance, in previous works on cuprates [19,20], it has been shown that the AB LEH parameters have a nonlinear dependence on band structure variables and crystal structure variables.
On (d), the explicit expression of y as a function of x i is necessary to understand whether y increases or decreases with increasing x i .[Such information is not captured by (b) alone.]As mentioned above, the combination of (b) and (d) allows to obtain y MODg .
On (f), the sparsity allows to facilitate the physical interpretation of the approached expression of y.Indeed, if an accurate approximation of y can be constructed from only a few x i , the distinct contributions of the x i to the dependence of y may be analyzed and discussed more easily.
Already existing procedures such as linear regression, polynomial regression, and also symbolic regression and sparse regression fulfill part of the above conditions (a-e), but not all of them.This is discussed below.
Linear regression fulfills (d) and (e), but not other points.Although it is the simplest approach (e) to obtain an explicit expression of y as a function of the x i (d), the dependence of y on x i is nonlinear in the general case as discussed before, so that (c) is not fulfilled.
Polynomial regression and symbolic regression fulfill (c) in addition to (d) and (e), but not other points.Polynomial regression has good approximation properties provided that y = h({x i }) is a continuous function of the x i : In that case, y can be approached uniformly by a polynomial of the x i according to the Stone-Weierstrass theorem.Symbolic regression allows to go beyond the polynomial regression without introducing assumptions on the expression of y as a function of the x i .However, these techniques do not allow to fulfill (a), (b) and (f) in the general case.
Sparse regression allows to perform polynomial and symbolic regression by fulfilling (f).Regularization techniques allow to penalize expressions of y that depend on many x i , allowing to construct an expression of y as a function of a reduced number of x i (f).However, they do not offer a clear framework to fulfill (a) or (b).
The HDE is designed to fulfill all conditions (a-f).First, the HDE allows to probe the completeness of the dependence (a) by examining the value of f ∞ [y, V] in Eq. ( 16).Second, the HDE allows to extract the hierarchy in dependencies of y on x i (b) thanks to the recurrent expression of the candidate descriptor x (g) in Eq. (11).In the HDE[y, V], when we increment g, we successively add terms to x opt in Eq. ( 14) that correspond to the higher-order dependencies of y on {x i } (the order increases with g).The competition between dependencies can also be analyzed by calculating the score of each variable in Eq. (20).Third, the wildcard operator in Eq. ( 12) encompasses polynomials of variables, so that the nonlinear dependence of y may be captured at an acceptable level of accuracy (c).(Still, note that its approximation properties are not perfect, as discussed below.)Fourth, the HDE allows to obtain an explicit expression of y as a function of x i (d).Fifth, the HDE procedure is simple and deterministic (e), and can be applied even if the size of the training set is relatively small (we should have N t ≥ 3).Sixth, the HDE procedure allows to preserve the sparsity by introducing no more than one variable x i in Eq. ( 14) when g is incremented.
Note that, in the HDE[y, V], the ratio usually decreases with increasing g, but the amplitude of the ratio does not necessarily decrease with increasing g.Namely, when incrementing g, the effect on the dependence of y on x i may be corrective (i.e.f (g) [y, V] increases by a small amount so that ∆ f (g) is small), but the effect on the amplitude of x opt may be significant: The amplitude of the additional term (which is encoded in |∆x opt (g) |) is not necessarily small.
Note that the approximation properties of the HDE are not perfect due to the compromise between approximability and simplicity in the present framework.Namely, we enforce sparsity by allowing only one variable x i in V to be introduced in x opt (g) when g is incremented, and the expression of ∆x opt (g) is kept as simple as possible to facilitate its interpretation.On one hand, this choice does not allow to represent all polynomials of x i , which requires to introduce several variables at once in x opt (g) when g is incremented.On the other hand, this allows to obtain the hierarchy between the dependencies of y on x i in a simple manner.In the scope of this paper, the consistency between results provided by the HDE and previous works suggests the HDE in its present form is reliable.Possible extensions of the HDE that involve more sophisticated expressions of ∆x opt (g) are left for future studies.
Finally, on (f), it should be noted that sparsity is distinct from dimensionality reduction.Dimensionality reduction techniques such as principal component analysis (PCA) [42] have been used in e.g.recent studies on iron-based SC materials [43].The PCA allows to reduce the size of the variable space V = {x i , i = 1..N V }, by extracting principal components that are linear combinations of the x i .(Details can be found in e.g.Ref. [44].)The m th principal component is denoted as The principal components are sufficient to reproduce the diverse distribution of values of x i in the training set, so that an approached expression of y may be constructed as a function of a few z (m) .However, in Eq. (C3), the weights v (m) i may be nonzero for many variables x i in the general case.Thus, if we try to express y as a function of {z (m) , m} instead of {x i , i = 1..N V }, the expression of y will depend on a few z (m) but may depend on many x i , so that the sparsity is not preserved.To interpret the expression of y, we discuss the dependence of y on the x i rather than on the z (m) , because the variables x i have a physical meaning (see Sec. III C).

Comments on the choice of fitness function
Here, we discuss in further detail the choice of the fitness function in Eq. ( 8).Namely, we remind the definition of the Pearson correlation coefficient that is used in Eq. ( 8), and we discuss the invariance property of the fitness function [Eq.( 9)].
Pearson correlation coefficient -The Pearson correlation coefficient ρ(y, x) between two variables y and x is defined as follows.The variables y and x are represented by two data samples {y[ j], j = 1..N t } and {x[ j], j = 1..N t }, where N t is the size of the training set.The sample mean of x, sample covariance of y and x and sample variance of x are defined as, respectively: and we calculate ρ(y, x) as The value of ρ(y, x) is a real number between −1 and +1.
If ρ(y, x) = 1 [ρ(y, x) = −1], then x and y are perfectly correlated (anticorrelated), and there exist k 0 and k 1 such that the equation Invariance property of the fitness function -The invariance property [Eq.( 9)] of the fitness function comes from the property of the Pearson correlation coefficient.[In Eq. (C8), we have The invariance property simplifies the search of the best candidate descriptor: The fitness function [Eq.( 8)] encodes only the linear dependence of y on x(p), which is sufficient to judge how well y is described by x(p). 5From the viewpoint of the fitness function, the candidate descriptor x(p) is equivalent to any other candidate descriptor in the equivalence class and we do not need to distinguish between elements of E[x(p)] during the optimization in Eq. (7).The values of k 0 and k 1 are determined after x opt has been determined, by performing the affine regression in Eq. ( 15).In addition, the value of the fitness function has a rather universal meaning (at least at a fixed value of N t ), which facilitates the judgment of the quality of the candidate descriptor x(p).This is a consequence of the invariance property [Eq.( 9)].For a given y and x(p), the value of f [y, x(p)] allows to deduce how well y is described by x(p), irrespective of the scale or order of magnitude of y and x(p).The description is perfect for f [y, x(p)] = 1.0.Empirically (but rather universally at least at N t = 36), we observe in Sec.IV that the description is almost perfect for f [y, x(p)] ≳ 0.98.[In this paper, we assume that x(p) ] ≲ 0.97, the description is good, but corrective higher-order dependencies of y on x i may be missing.If f [y, x(p)] ≃ 0.6 − 0.9, the description is rough or very rough, but x(p) captures correctly the MOD.Typically, for the MOD2s discussed in Sec.II and Sec.IV, the values of f (2) [y, x opt (2) ] are above 0.6, and most of the time above 0.8−0.9.There are also particular cases in which f (2) [y, x opt (2) ] ≳ 0.98, so that the MOD2 is sufficient to describe y entirely.
(i) The wildcard operator can represent basic algebraic operations.First, multiplication and division are accounted for by: in which the arbitrarily large yet finite factor ζ is traced out by using the invariance property of the fitness function.Furthermore, Eq. (C11) encompasses products and ratios between x and exponents of x ′ .Second, addition and subtraction are accounted for by: for ζ = ±1.Furthermore, Eq. (C12) encompasses any linear combination of x and x ′ if |ζ | ̸ = 1.Third, exponents and polynomials of any variable x are accounted for by: (ii) The wildcard operator can represent the reset and identity operators.The reset operator is defined as and allows to replace x by x ′ .If x ′ = x, Eq. (C15) becomes the identity operator.[The possible values of (ζ opt , β opt , α opt ) which represent the identity operator are (ζ → +∞, β = 1, α = 1) and also (ζ > 0, β , α = 0).]The fact that the wildcard operator is able to mimic the identity operator guarantees that the best candidate descriptor from the generation g is included in the candidate descriptor space for the generation g + 1.In particular, this implies Eq. ( 13), which is a desired property of the fitness function: Given the best candidate descriptor at g, the iteration from g to g + 1 adds the g + 1 th order dependence, and taking into account the g + 1 th order dependence necessarily improves the descriptor.
(iii) The wildcard operator can represent a perturbation of x by x ′ .For small |ζ | and α = 1, we have in Eq. ( 12), so that x ⋆ (ζ ,β ,α=1) x ′ yields x corrected by a small term that depends on x ′ .

Computational details of the HDE procedure
Here, we give computational details of the P[y, V] procedure that is employed to obtain x opt (N) in Eq. ( 14).First, we use a finite number N of generations.In this paper, we use N = 15 in Sec.IV C, and we check that the fitness function varies slowly with increasing N at N ≃ 15, so that it is reasonable to stop at N = 15.
Second, at fixed g, we determine p opt that maximizes f [y, x(p)] in Eq. ( 7) as follows.We scan C V by computing f [y, x(p)] for each x(p) in C V .Note that this is computationally expensive as discussed below; however, this allows to avoid falling into any local extremum of f [y, x(p)].Another possibility to reduce the computational cost would be to prepare an initial guess for p and then optimize p by using e.g. a gradient descent algorithm.Such extensions and the detailed study of the dependence of the result on the initial guess for p are left for future studies.
At fixed g, the variational parameters are p g = (α g , β g , ζ g , i g ) as mentioned in Sec.III A. The number of indices i g is given by the total number of variables N V in V = {x i }.As for β , we consider only N β = 2 values, which are β = 0, 1; this choice preserves the polyvalence of the wildcard operator as illustrated in Fig. 14.As for α and ζ , we use a finite grid of values for the optimization of the fitness function.The numbers of values of α and ζ in the grid are denoted as N α and N ζ , respectively.Below, we detail the choice of this grid.
The computational cost increases rapidly with N α and N ζ : At fixed g, the number of candidate descriptors to be evaluated is In practice, we reduce the computational cost by employing a multi-step optimization, as detailed in (i,ii) below.
(i) First, we optimize (α g , ζ g ) on a coarse grid together with (β g , i g ).The coarse grid is the following:  (ii) Then, we keep the values of β opt g and i opt g that were determined in (i), and we refine the optimization of (α g , ζ g ) on a finer grid, which is constructed iteratively as follows.We introduce the iteration index j, starting from j = 2.The fine grid is built around α .We iterate up to j = 3.

Numerical aspects and pitfalls
Here, we discuss a few limitations and numerical pitfalls in the HDE procedure.First, note that if β >0, the values of x [or x opt (g−1) in Eq. ( 11)] must be nonzero when calculating Namely, if the variable x is represented by the sample {x[ j], j = 1..N t } introduced in Appendix C 2, we should have x[ j] ̸ = 0 for all j.This implies x i [ j] ̸ = 0 for all x i in V, because we have x opt (1) = x α 1 i 1 : If x i 1 has zero values, then x opt (1) will have zero values if α 1 > 0, or x opt (1) will diverge if α 1 < 0. In our calculations, some of the variables x i in V have zero values: For instance, in V 1 , the variable R A ′ has zero values if A'= ∅.Thus, we introduce an infinitesimal offset δ off > 0 in all variables.(Namely, we replace x i [ j] by x i [ j] + δ off .)We choose the value δ off = 10 −4 which is small enough to be negligible with respect to the difference between values of x i Third, the maximal value of f[y, x (g) ] may be obtained for several candidate descriptors with the same i g but different α g .These are in the same equivalence class [Eq.(C9)], but this is nontrivial.This is discussed in Sec.S4 of SM [21].In this case, we choose to apply the following convention: We choose α such that |α| is minimal and α > 0. (Then, if g = 2 and if there are several values of ζ for which the fitness function has the same value, we choose ζ such that |ζ | is minimal and ζ > 0.) In practice, if we use this convention, the optimization of α on the coarse grid then on the fine grid yields α = |α| min = 0.01.
A concrete example is |Z X |, whose value is either 1 and 2. In this case, f [R, |Z X | α ] = 0.503 for all values of α ̸ = 0.This is why we have e.g.α 1 = 0.01 in Eq. (5).Note that this choice should not affect the physical interpretation of the MOD of R on |Z X | in Eq. (5).
and if we choose another convention, e.g.we choose α such that |α| is minimal and α < 0, then we obtain Finally, note that the complexity of the expression y N of y in Eq. ( 15) increases with increasing g.However, incrementing g up to a large value of N allows to probe the completeness of the dependence of y on V, as discussed in Sec.III A. The principal microscopic mechanism of the dependence of y on V is usually contained in the MODg for g ≲ 3.
Appendix D: Details on the gMACE procedure 1. Computational details of the gMACE procedure Here, we give computational details of the ab initio calculations in the gMACE procedure.We also detail the procedure that is employed to construct the AB MLWO.a. Structural optimization and DFT calculation -We use Quantum ESPRESSO [45,46] and optimized normconserving Vanderbilt pseudopotentials [47,48] with the GGA-PBE functional [25].To model hole doping, we use the virtual crystal approximation [31] as done in [19,20] and as mentioned in the main text: The pseudopotential of A or A' cation is interpolated with that of the chemical element whose atomic number is that of A or A' minus one.We use a planewave cutoff of 100 Ry for the wavefunctions.The full Brillouin zone is sampled by using a k-point grid of size 8 × 8 × 8 in the structural optimization, 12 × 12 × 12 in the self-consistent DFT calculation, and 6 × 6 × 6 in the non-selfconsistent DFT calculation (8 × 8 × 4 if A'=Hg 1−δ Au δ ).We use a Fermi-Dirac smearing of 0.002 Ry. b.Crystal parameters and symmetry -The primitive vectors of the Bravais lattice and the positions of atoms are given in Table II.During the structural optimization, we optimize the values of a, c, d z A and d z X altogether.Then, we use the optimized values in the self-consistent and non-self-consistent DFT calculations.
For completeness, note that the high-symmetry structure that is considered in Table II may be lowered in experiment.Namely, atoms may undergo displacements around their ideal positions listed in Table II, which creates incommensurate modulations in the crystal structure.The origin of these displacements is the nonideal Goldschmidt tolerance factor of the perovskite-like structure of the cuprate (i.e. the ratio between the radii of atoms in the crystal is not ideal), which causes a geometric mismatch between the block layer and the CuO 2 plane.This happens e.g. in the case of Bi2201 [49,50].Also, in the case of Sr 2 CuO 2 F 2 , the F atoms are distorted, and the doping may introduce excess F at different positions [8].
It is possible to account for the structural distortion in a simplified manner, by considering a distortion that is restricted to the unit cell and ignores the incommensurate character of the distortion, as done in [19] in the case of Bi2201.(See Appendix C of [19].)However, in Bi2201, the effect of the distortion |t 1 | and u is small [19]: If we compare the AB LEH obtained from the crystal structure with and without distortion, |t 1 | does not vary and U varies by no more than 3%.In the present paper, for simplicity, we do not consider the distortion and always assume the ideal atomic positions in Table II.x y z c. Wannier orbital -Here, we detail the procedure to construct the AB orbital.We use the RESPACK code [19,51].The initial guess consists in an atomic Cu3d x 2 −y 2 orbital centered on the Cu atom in the unit cell.In the outer window, the spillage functional [29] is minimized to obtain the AB subspace and band dispersion, then the spread functional [28] is minimized to obtain the AB MLWO.In previous works [17][18][19][20], the outer window that is used to construct the AB orbital consists in the M space, from which the N B lowest bands are excluded to avoid catching the Cu3d x 2 −y 2 /O2p σ bonding (B) character.However, the characteristics of the AB MLWO and values of AB LEH parameters depend slightly on N B : For instance, the AB LEH parameters may vary by a few percent if N B varies by 2 or 3 [17,19].Although this small dependence of the AB MLWO on N B does not change the physics of the AB LEH, it may prevent the very accurate comparison between AB LEHs that are derived from different CFs.
To estimate accurately the materials dependence of the AB LEH, we employ a modified procedure to construct the AB orbital, in which we remove the B character without excluding the B bands from the outer window.This procedure is based on the antibonding-bonding (ABB) transformation [24], and is described below.We use the whole M space as the outer window, and we consider three initial guesses for the MLWOs: The Cu3d x 2 −y 2 atomic orbital centered on Cu, and the O2p σ (O'2p σ ) orbital centered on O (O').Then, we minimize the spillage functional [29] to obtain the band dispersion that corresponds to the three MLWOs with Cu3d x 2 −y 2 and O/O'2p σ character.This band dispersion consists in the partly filled AB band plus the two fully occupied B bands; see e.g.Fig. 8 in [19] for an illustration.Then, we discard the B band dispersion and we keep only the AB band, which contains the AB subspace and has been determined without introducing the dependence on N B .Finally, we obtain the AB MLWO from the AB band by considering the AB band as the outer window, reinitializing the Cu3d x 2 −y 2 initial guess and minimizing the spread functional [28].(Note that the ABB transformation allows to derive a three-orbital Hamiltonian that includes the AB orbital and two B orbitals; here, contrary to the ABB transformation, we discard the B subspace.)d.Polarization We use the RESPACK code [19,51].The cRPA polarization at zero frequency is expressed as [51]: and in the above equations, η is an infinitesimal positive number, q is a wavevector in the Brillouin zone, G, G ′ are reciprocal lattice vectors, nk is the Kohn-Sham one-particle state with energy ε nk and wavefunction ψ nk , T nk = 1 if nk belongs to the AB band and T nk = 0 else.The cRPA effective interaction is deduced as in which v is the bare Coulomb interaction.We use a planewave cutoff of 8 Ry for the calculation of the cRPA polarization and dielectric function.

Choice of the variables in V 3
As a complement to Sec.III C, the choice of the variables in V 3 is discussed below.The DFT band structure is entirely described by the Kohn-Sham energies ε nk and orbitals ψ nk .However, it is too complex to consider the whole set {ε nk , ψ nk } as the variable space V 3 .Instead, we consider variables that capture the essential characteristic energies in the band structure.To choose the variables, we use as a guide the M space.Prioritizing the M space for the choice of variables is natural, for two reasons.First, the M space determines the characteristics of the AB Wannier orbital including |t 1 | and v [20], because the M space contains the Cu3d x 2 −y 2like and O2p σ -like bands and the AB band.Second, the M space plays a prominent role on the cRPA screening and thus on R = U/v.Indeed, the screening increases when the charge transfer energies between occupied and empty bands decrease [see Eqs.(D1) and (D2)], and the charge transfer energies are the smallest for the occupied M bands.

Values of intermediate quantities within the gMACE procedure
Here, we give the values of intermediate quantities in the variables spaces V s with s = 1..4.On V 1 , the values of R A , R X , R A ′ , |Z X | and Z A that are considered in this paper are represented in Fig. 15.(The numerical values are given in Sec.S1 of SM [21].)The values of R A , R X and R A ′ are the crystal ionic radii values from [36].(For completeness, we mention that ions are assumed to be 6-coordinate in the calculation of these values.)For hole-doped compounds, the partial ion substitution is accounted for as follows.If A'= ∅, the ion A becomes A 2−δ Ãδ where the atomic number of Ã is that of A minus one; otherwise, the ion A' is Hg 1−δ Au δ .In this case, we interpolate linearly the values of ionic radii; namely, we consider   (36)], s (2) [ε e , n AB ] = 1 and n AB is not in competition with other variables.This corresponds to the scenario (S1).Thus, |Z X | and n AB correspond respectively to  3) and (5), we discuss the n AB dependence of u and R when the CF at δ = 0 is fixed.The values of u and R for all compounds in the training set are shown in Fig. 18(a).For X = F, Cl and especially (X, A) = (F, Ba), (Cl, Ba) and (Cl, Sr), we observe a sharp decrease in u between 10% and 20% hole doping (i.e. between n AB = 0.9 and n AB = 0.8).This decrease in u is caused by the sharp decrease in R as seen in Fig. 18(a,b).Although the decrease in R with increasing δ (or decreasing n AB ) is consistent with the MOD2 in Eq. ( 5), the sharper decrease in R for X = F, Cl compared to X = O is not captured by Eq. ( 5).For completeness, we discuss the origin of the sharper decrease in R for X = F, Cl in detail here.
The decrease in the δ dependence of R may be quantified by considering the quadratic interpolation where R 0 , R  16), and also (ii) the lower value of the M space bandwidth W M for X = F, Cl.On (ii), we see in Fig. 18(d) that the three compounds with (X, A) = (F, Ba), (Cl, Ba) and (Cl, Sr) have the lowest W M at δ = 0.2.The above discussed dependence of |R 2 | on W M (ii) is interpreted as follows.If W M is lower, then the charge transfer energies between the occupied states in the M space and the empty states in the M space (namely, the empty part of the AB band) are smaller.Thus, the intra-M space cRPA screening will be stronger.A rough scaling of the intra-M space screening is 1/W M [see Eq. (D1)], so that the smaller W M , the more R will decrease when W M further decreases.[And, for X = F, Cl, W M decreases with increasing δ as seen in Fig. 18(d).]On the other hand, in X = F, Cl, the screening channel between the occupied states and empty states outside M is relatively weak compared to X = O due to the larger ε e (see Fig. 16).This suggests the intra-M space screening dominates over the other screening channels for X = F, Cl. (Note that, besides the decrease in W M , the decrease in |ε o | with increasing δ also participates in increasing the intra-M space screening.) The higher |R 2 | in X = F, Cl compared to X = O implies a nontrivial point on the origin of the SC in oxychlorides and oxyfluorides.[Confirmation of this point requires to improve the derivation of the AB LEH at the cGW -SIC+LRFB level, which is left for future studies.]Even though u and R are higher in the undoped compound for X = F, Cl compared to X = O [according to Eq. ( 3) and ( 5)], u may become lower in the hole doped compound at optimal hole doping for X = F, Cl compared to X = O due to the higher |R 2 |.In particular, u may fall into the weak-coupling regime (u ≃ 6.5 − 8.0) at δ ≃ 0.2 for X = F, Cl, which may correspond to the opti- ≃ 27 K [12], and this is realized at δ = 0.18, which is close to δ = 0.2 at which the sharp decrease in R happens.If u is in the weak-coupling regime at δ = 0.18, then the lower T c compared to other cuprates is caused by the decrease in F SC with decreasing u in the weakcoupling regime.First, we discuss (a).Because ε o is the highest energy of the occupied bands outside the AB band, the value of |ε o | mainly depends on the orbitals that form the density of states near the Fermi level, excluding the Cu3d x 2 −y 2 and O2p σ orbitals.Thus, we examine the character of the occupied bands near the Fermi level, for the three above CFs.To do so, we represent the partial density of states (pDOS) for several types of ML-WOs in Fig. 19.Near the Fermi level, the Cu3d x 2 −y 2 /O2p σ pDOS corresponds to the AB subspace, and the character of the highest occupied bands outside the AB band is given by the pDOS other than Cu3d x 2 −y 2 /O2p σ .Outside the AB band, the density of states near the Fermi level is dominated by the O2p z pDOS and Cu3d yz/zx pDOS (see Fig. 19).Thus, the value of |ε o | mainly depends on the O2p z and Cu3d yz/zx ML-WOs.In particular, |ε o | increases if the O2p z and Cu3d yz/zx orbitals are deeper in energy.Now, we discuss (b), i.e. why the O2p z and Cu3d yz/zx orbitals are deeper in energy when R A decreases.If R A is reduced, then a is reduced [see Eq. ( 27)], and the in-plane O is closer to the A cation.Thus, the MP + from the A cation that is felt by the in-plane O is stronger.[See Fig. 9(a) for an illustration.]Thus, the O2p MLWOs in the M space are stabilized, i.e. their onsite energy is reduced.This shifts the O2p pDOS downwards, i.e. farther from the Fermi level.With respect to the O atom in the unit cell, the positions of the four nearest A cations (in Cartesian coordinates) are ±a/2y ± d z A z (y and z are unitary vectors along the y and z directions that are considered in Table II): The O2p σ orbital extends along x direction and thus avoids the A cations, whereas the O2p z orbital extends along z direction, so that the O2p z electrons are closer to the A cation compared to the O2p σ electrons.This explains why the O2p z electrons are prominently affected by the MP + from the A cation.Similarly, the Cu3d yz/zx orbitals extend along the z direction, and may be prominently affected by the MP + from the A cation.The above discussion is supported by the values of the onsite energies of the MLWOs: From (i) to (iii), ε O p decreases by 0.97 eV, whereas ε O p z decreases by 1.16 eV, so that the O2p z orbital is indeed more stabilized than the O2p σ orbital.Also, ε Cu x decreases by 0.26 eV, whereas ε Cu yz/zx decreases by 0.48 eV, so that the Cu3d yz/zx orbital is indeed more stabilized than the Cu3d x 2 −y 2 orbital.
II. OVERVIEW: MAIN-ORDER DEPENDENCE OF AB LEH PARAMETERS ON CHEMICAL FORMULA Here, we give an overview of the MOD of AB LEH parameters on the CF that is obtained by applying the gMACE+HDE.(Details on the results are given later in Sec.IV, and prescriptions to optimize T opt c based on these results are proposed in Sec.V.) The MODs of |t 1 | and u are summarized below in (I,II) and illustrated in Fig. 1.

FIG. 1 .
FIG. 1. Upper panel: Simplified representation of the CuO 2 plane and the surrounding crystalline environment.We show the square lattice formed by the Cu atoms in the CuO 2 plane (the in-plane O atoms are not shown), and the isosurface of the AB orbital centered on one of the Cu atoms (yellow is positive, blue is negative).We also show the A cations and apical X anions near the CuO 2 plane.Lower panel: Representation of the MODs of the AB LEH parameters on the CF obtained from gMACE+HDE.We show the MOD of |t 1 | [Eq.(2)] and u = U/|t 1 | [Eq.(3)] together with the onsite bare Coulomb interaction v [Eq.(4)] and the screening ratio R = U/v [Eq.(5)].R A and Z A (R X and Z X ) are the crystal ionic radius and ionic charge of the A cation (apical X anion).n AB is the average number of electrons per AB orbital.

FIG. 2 .
FIG. 2. Summary and illustration of the gMACE scheme and variables that are considered in the gMACE+HDE procedure.On s = 2, the primitive cell vectors and atomic positions are given in Appendix D 1. On s = 3, we show the band structure in the case of HgSr 2 CuO 4 [the high-symmetry points are Γ = (0, 0, 0), D= (1/2, 0, 0) and X= (1/2, 1/2, 0) in coordinates of the reciprocal lattice].The dashed (solid) black bands are those outside (inside) the M space.The red band is the AB band, and the cyan bands are the 16 other M bands after disentanglement from the AB band.We also show the principal characteristic energies that are discussed in the main text: Energy ε e (ε o ) of the lowest empty band (highest occupied disentangled M band) outside the AB band, and onsite energies ε Cu x , ε O p , ε X p x and ε X p z of the maximally localized Wannier orbitals (MLWOs) whose character is Cu3d x 2 −y 2 (abbreviated as Cux), O2p σ (Op), X2p x (Xp x ) and X2p z (Xp z ).(The Cux/Op charge transfer energy is ∆E xp = ε Cu x − ε O p ).We also show the isosurfaces of the Cux, Op, Xp x and Xp z MLWOs (yellow is positive, blue is negative).On s = 4, we show the isosurface of the AB MLWO centered on the Cu atom, and a schematic illustration of the square lattice formed by the Cu atoms in the CuO 2 plane.
where a, c and c ⊥ are the cell parameters, and d z A (d z X ) is the distance between the CuO 2 plane and the A cation (apical X anion).The coordinates of primitive vectors and atoms in the unit cell are given and discussed in Appendix D 1.All variables in V 2 are expressed in Å throughout this paper.
FIG. 4. Ab initio values of |t 1 |, u = U/|t 1 |, v and R = U/v obtained by employing the gMACE for all CFs in the training set.For each color point, the corresponding CF is shown in Fig. 3.

3 and x i opt 4 ,
FIG. 5. Summary of the MODs between quantities obtained within the gMACE+HDE.Upper panel: MOD2s of |t 1 |, v and R on the CF variables in V 1 .[The items (I), (III) and (IV) have been summarized in Sec.II, and are discussed in detail in Sec.IV.] Lower panel: MODs between intermediate quantities within the gMACE.[The items (1-10) are discussed in detail in Sec.IV.]

1 (
884 [see Fig. 6(a)], and the dependence of |t 1 | on x opt (2) is shown in Fig. 6(b).Note that the dependence of |t 1 | on R A is equally important to that on R X in Eq. (2), even though R X corresponds to x i opt see the score analysis in Appendix E).

( 2 )
FIG. 7. Details on the dependence of |txp | on V 2 , ∆E xp on V 2 , a on V 1 , d z X on V 1 , and also |t 1 | on V 1 and ∆E xp on V 1 for completeness.Panel (a): Values of f (g) [y, V] in the HDE[|t xp |,V 2 ], the HDE[∆E xp ,V 2 ], the HDE[a,V 1 ], the HDE[d z X ,V 1 ], the HDE[|t 1 |,{a}], and the HDE[∆E xp ,V 1 ].Panels (b-f): Representation of the MOD1 of |t xp | on V 2 [Eq.(25)], the MOD2 of ∆E xp on V 2 [Eq.(28)], the MOD2 of d z X on V 1 [Eq.(29)], the MOD2 of a on V 1 [Eq.(27)], and the MOD1 of |t 1 | on {a} [Eq.(26)].Panels (g-k): Dependence of |t xp | and |t 1 | on their respective x opt (1) , and ∆E xp , d z X and a on their respective x opt (2) .[Thepanels (g) to (k) correspond to the panels (b) to (f), respectively.]The CF that corresponds to each color point is shown in Fig.3, and the solid black line shows the linear interpolation.

FIG. 9 .
FIG. 9. Schematic illustration of the positive Madelung potential (MP + ) and negative Madelung potential (MP − ) created by the cations and anions within the crystal.Panel (a): MP + created by the cation A and felt by the nearest neighbor apical X (green arrows), in-plane O (red arrows), and Cu (blue arrows).Panel (b): MP − created by the apical X anion and felt by the nearest neighbor Cu (blue arrow), in-plane O (red arrows), and cations (green arrows).Panel (c): MP − created by the in-plane O anion and felt by the nearest neighbor Cu (blue arrows) and other in-plane O (red arrows).In the panel (a), the MP + created by A and felt at a distance d of A scales as V A (d) = Z A /d.The distance d A,X =

1 and x i opt 2 unambiguously
(see Appendix E).

FIG. 12 .
FIG. 12. Dependence of y = |ε X p x |, |ε o | and ε e on V 1 .Panel (a): Values of f (g) [y, V] in the HDE[y,V 1 ].Panels (b-d): Representation of the MOD2 of y on V 1 in Eqs.(34), (35),(36).Panels (e-g): Dependence of y on x opt (2) which corresponds to the panels (b-d), respectively.The CF that corresponds to each color point is shown in Fig.3, and the solid black line shows the linear interpolation.
|ε X p z | (in Appendix D 3), we see that |ε X p x | ≃ k 0 + k 1 |ε X p z |has an affine dependence on |ε X p z |, but the coefficients k 0 are different for A' = ∅ and A' = Hg 1−δ Au δ whereas the coefficients k 1 are nearly identical.Namely, the affine regression yields

( 9 )
Dependence of |ε o | on V 1 -The MOD of |ε o | on n AB and R A in Eq. (

B
. Prescriptions to optimize T opt c By considering the MODs given in Sec.II, we propose prescriptions to optimize T opt c in future design of SC cuprates.The overall strategy is to maximize |t 1 | while keeping u as close as possible to its optimal value u opt ≃ 8.0−8.5.Indeed, T opt c increases with both |t 1 | and F SC in Eq. (1), and the universal u dependence of F SC has a maximum at u opt (see Sec. I).To realize u = U/|t 1 | = vR/|t 1 | as close as possible to u opt , the criterion vR/|t 1 | ≃ 8.5 (39)
VI. CONCLUSION We have proposed the universal CF dependence of the AB LEH parameters |t 1 | and u in single-layer cuprates, by proposing the HDE procedure and applying it to analyze the results of ab initio calculations of the AB LEH for various singlelayer cuprates.The results and especially the MODs given in Sec.II provide insights to optimize T opt c in future design of SC cuprates as proposed in Sec.V.
opt c .Note that for the CFs in the training set, the M space has a similar structure, which facilitates the comparison between the variables in V 3 .Namely, the number of bands in the M space is N M = 17 for all compounds in the training set.(See the band structures in Sec.S2 of SM [21].)In other single-layer compounds such as Bi2201 with T opt c ≃ 10 − 40 K [1, 41], we have N M = 23 due to the presence of 6 additional bands from the BiO block layer.However, we not consider cuprates with a Bi2201-like CF and crystal structure, in order to simplify the comparison between compounds.This is not expected to weaken the generality of the result, because the diverse distribution in experimental T opt c ≃ 24 − 94 K is already reproduced by the CFs in the training set, as mentioned above.Possible extensions of the training set to other cuprates including multi-layer cuprates and Bi2201-like single-layer cuprates are left for future studies.
Appendix C: Details on the HDE procedure 1. Motivation of the HDE
C22) so that the physical interpretation of the MOD1 (R increases when |Z X | decreases) remains the same irrespective of the selected convention.The relatively large values of |k 0 | and |k 1 | in the above MOD1 and in the MOD2 [Eq.(5)] with respect to the ab initio values of R ≃ 0.22 − 0.34 are a consequence of the small value of |α|.Still, the range of ab initio values of R is correctly reproduced by the MOD2, as seen in Fig. 1.
TABLE II.Primitivevectors a, b, c of the Bravais lattice and atomic positions in Cartesian coordinates.We consider c ⊥ = 0 if A' = Hg 1−δ Au δ and c ⊥ = a/2 if A' = ∅.The atomic positions are entirely determined by a, c, c ⊥ , d z A and d z X .Note that there are two O atoms, two A atoms and two X atoms in the unit cell.The first and second O atoms in the unit cell are denoted as O and O', respectively.
Appendix F: Hole doping dependence of screeningHere, as a complement to Sec.IV C and the MOD2 of u and R in Eqs. ( 1 and R 2 are determined entirely by the ab initio values of R at δ = 0.0, 0.1 and 0.2: R 0 is R at δ = 0 in Fig. 18(b), and the values of R 1 and R 2 are given in Fig. 18(c).The decrease in R is encoded in R 2 < 0 for all compounds, and the sharp decrease in R for (X, A) = (F, Ba), (Cl, Ba) and (Cl, Sr) is reflected in the high value of |R 2 | compared to the other compounds.Other compounds with X = F, Cl also have higher values of |R 2 | compared to X = O.Namely, we have |R 2 | ≤ 0.55 for X = O, |R 2 | ≥ 1.14 for X = F, Cl, and |R 2 | ≥ 3.55 for (X, A) = (F, Ba), (Cl, Ba) and (Cl, Sr).Possible causes of the higher value of |R 2 | are (i) the sharp decrease in |ε o | from δ = 0.1 to δ = 0.2 for X = F, Cl (see Fig.

FIG. 18 .
FIG. 18. Panel (a): Values of u and R taken from Fig. 4. Panel (b): Values of R as a function of hole doping δ = 1 − n AB .The dashed curves show the quadratic interpolation of the δ dependence of R by using Eq.(F1).Panel (c): Values of the coefficients R 1 and R 2 in Eq. (F1).Panel (d): Values of W M as a function of R 2 .For each color point, the corresponding CF is shown in Fig. 3.In the panel (c), the colors of the crosses correspond to the values of A, X and A' in Fig.3.
957, so that a describes |t 1 | reasonably, but not perfectly.This is because |t 1 | has not only a dominant dependence on |t xp | but also a small dependence on ∆E xp , and ∆E xp is not described entirely by a as seen later in (4).The MOD1 of |t 1 | on a is dependence of |t 1 | by fixing the other CP values; the present result is more general because it accounts for the materials dependence of other CP values via the structural optimization.
[20] | depends on both |t xp | and ∆E xp [Eq.(24)], and ∆E xp depends on the crystalline environment outside the CuO 2 plane contrary to |t xp |. [For instance, as seen later in (5), the MOD2 of ∆E xp depends not only on a but also on d zX .]Theresultin[20]captures the a

TABLE I .
List of experimentally confirmed SC materials whose chemical formula is included in the training set.T exp cis the experimental value of T c from the reference given in the third column.