Interplay between Zhang-Rice singlets and high-spin states in a model for doped NiO$_2$ planes

Superconductivity found in doped NdNiO$_2$ is puzzling as two local symmetries of doped NiO$_2$ layers compete, with presumably far-reaching implications for the involved mechanism: a cuprate-like regime with Zhang-Rice singlets {\cblue is replaced by local triplet states at realistic values of charge-transfer energy, which would suggest a rather different superconductivity scenario from high-$T_c$ cuprates}. We address this competition by investigating Ni$_4$O$_8$ clusters with periodic boundary conditions in the parameter range relevant for the superconducting nickelates. With increasing value of charge-transfer energy we observe upon hole doping the expected crossover from the cuprate regime dominated by Zhang-Rice singlets to the local triplet states. We find that smaller charge-transfer energy $\Delta$ is able to drive this change of the ground state character when realistic values for nickel-oxygen repulsion $U_{dp}$ are taken into account. For large values of the charge-transfer energy, oxygen orbitals are less important than in superconducting cuprates as their spectral weight is found only at rather high excitation energies. However, a second Ni($3d$) orbital can easily become relevant, with either the $xy$ or the $3z^2-r^2$ orbitals contributing in addition to the $x^2-y^2$ orbital {\cblue to the formation of triplet states. In addition,} our result that $U_{dp}$ (acting between Ni and O) favors onsite triplets implies that correlation effects beyond purely onsite interactions should be taken into account when obtaining effective two-band models.


I. INTRODUCTION
Two-dimensional (2D) nickelates such as LaNiO 2 have been theoretically proposed long ago [1] as candidate materials for unconventional superconductors, but only recently has superconductivity been found in Sr-doped NdNiO 2 thin films [2]. This discovery could be classified as fulfilling the paradigm of high-T c superconductivity in systems similar to cuprates where both e g symmetries contribute at the Fermi surface [3]. Indeed, nickelate heterostructures were considered to be the most promising [4,5], but it took another decade until the superconductivity was found in quasi-2D films [2]. Shortly after the discovery of superconductivity in doped NdNiO 2 , the search for its mechanism began and it became clear that once again we have to do with unconventional superconducting (SC) materials. But the situation in the theory is somewhat similar to SC cuprates, where after their discovery in 1986 [6], the mechanism of high-T c superconductivity remains still puzzling [7,8].
Recently several groups calculated the electronic structure and tried to understand the essential differences to cuprate physics. Two remarkable differences between the CuO 2 planes of correlated insulator La 2 CuO 4 and the NiO 2 planes in NdNiO 2 became clearly evident. First, the nickelates are 'self-doped' due to the their rare-earth bands, so that an 'undoped' compound does not correspond to an undoped NiO 2 layer [9,10]. Regardless of the level of self-doping realized in superconducting NdNiO 2 , we focus here on understanding the undoped NiO 2 layers, which we regard as a kind of idealized 'parent compound'. The argument for studying a NiO 2 plane is a certain effective decomposition into a threedimensional band including rare-earth states and a more 2D * a.m.oles@fkf.mpi.de band of x 2 − y 2 character [11]. Second, the NiO 2 layers miss the apical oxygen ions present in the CuO 2 layers. However, the NiO 2 plane of the infinite 2D layer in NdNiO 2 is similar to the CuO 2 plane of CaCuO 2 , where the apical oxygens are missing as well. One might argue that the properties of these planes would be then similar but they in fact reflect two different parameter regimes of the charge-transfer model which describes them both.
In spite of great similarity between CuO 2 and NiO 2 planes, there are thus substantial differences in the electronic structure. While La 2 CuO 4 is a charge-transfer system, the NiO 2 planes in NdNiO 2 have large charge-transfer energy ∆ which indicates more a Mott-Hubbard system [12]. This is also confirmed in other works [9,[13][14][15][16][17][18] and thus Ni-O hybridization plays here a less important role than in cuprates. Arguments were also given that electronic structure in nickelate superconductors could be reduced to a one-band Hubbard model [19]. The undoped system has one hole in x 2 − y 2 Ni orbital which is the |b 1 symmetry state but Zhang-Rice (ZR) singlets (S = 0) |b 1 L b1 are much weaker in doped systems. These states compete with high-spin (S = 1) |a 1 b 1 states when doping increases [12].
But perhaps the most important difference between the above two classes of materials is that there is only one band which crosses the Fermi level in cuprates, while two bands cross the Fermi level in nickelates [3,17,20,21]. This suggests that indeed both e g orbitals could contribute at finite Sr doping to the properties of NdNiO 2 . In the impurity model one finds therefore a transition from the singlet to triplet local configuration at doped Ni 2+ ion [12], and we shall investigate here how this changes by going to a more extended Ni-O system than a single NiO 4 unit, i.e., beyond the impurity model. Furthermore, we show below that the transition to the regime where high-spin states contribute depends on the value of intersite Coulomb repulsion U dp .
Another question which we want to address here is the nature of high-spin states in doped materials. In a cuprate model the orbital which is close to the top x 2 − y 2 orbital is the second e g state of 3z 2 − r 2 symmetry (called also z 2 below), and S = 1 states could form by hole doping [22]. They would compete with ZR singlets [23]. Recent electronic structure results for (Nd,Sr)NiO 2 found using a combination of dynamical mean-field theory of correlated electrons and band-structure methods indicate a remarkable orbitalselective renormalization of the Ni 3d bands [24]. However, for a NiO 4 plaque in NdNiO 4 crystal fields are quite different and the orbital xy is the first orbital below x 2 − y 2 orbital occupied by one hole, while the z 2 orbital is at the bottom of the orbital states [25]. This sequence of the orbital states at a Ni ion was recently confirmed by quantum chemistry calculations [26]. It is thus challenging to ask to which 3d orbital doped holes will go in a doped material and whether this could have any physical consequences.
The paper is organized as follows. In Sec. II we explain the charge transfer model and its parameters. The methods to analyze finite systems, i.e., exact diagonalization and the variational cluster approximation (VCA) are introduced in Sec. III. The numerical results are presented and discussed in Sec. IV. The paper is summarized in Sec. V. In the Appendix we analyze briefly the consequences of the partial filling of xy orbital on Ni ions and conclude that the general conclusions concerning the possibility of the participation of high-spin states in the ground state are not affected.

II. THE CHARGE-TRANSFER MODEL
We introduce the multiband d − p Hamiltonian for a NiO 2 plane where we consider a 2D 2 × 2 cluster shown in Fig. 1 (with periodic boundary conditions) which includes four orbitals per NiO 2 unit cell: two e g orbitals {3z 2 − r 2 , x 2 − y 2 } at each Ni + ion and one bonding 2p σ orbital (either 2p x or 2p y ) at each oxygen ion in the 2D plane (12 is the total number of ions in the periodic Ni 4 O 8 cluster), Here the first two terms in the Hamiltonian (1) stand for the kinetic energy: H dp includes the d − p hybridization ∝ t pd and H pp includes the interoxygen p − p hopping ∝ t pp , whered † mα,σ (p † jν,σ ) is the creation operator of an electron at nickel site m (oxygen site i) in an orbital α. At Ni ions α ∈ {z,z}, where z andz stands for 3z 2 − r 2 and x 2 − y 2 orbitals, while at O ions ν ∈ {x, y} stands for p x and p y orbital, with up or down spin, σ ∈ {↑, ↓}. The model includes two 3d orbital states of e g symmetry at Ni + ions, and one 2p σ bonding oxygen orbital state at oxygen ions in the (a, b) plane, either 2p x for the Ni-O bond mi a or 2p y for the bond mi b. We study an isolated NiO 2 plane in NdNiO 2 , so do not consider 2p z orbitals at apical oxygen positions.
In the following we will use shorthand notation, and instead of {x 2 − y 2 , 3z 2 − r 2 } we shall write {(z), (z)} -this emphasizes the fact that z axis is chosen as the quantization axis for this e g orbital basis, while "()" brackets are used here to distinguish two Ni(3d) orbitals from O(2p σ ) orbitals, labeled as {x, y}. The elements of the matrices {t mα;jν } and {t iµ;jν } are assumed to be non-zero only for nearest neighbor nickeloxygen d − p pairs, and for nearest neighbor oxygen-oxygen p − p pairs. The next nearest neighbor hopping is neglected. (The nonzero {t mα;jν } and {t iµ;jν } elements are standard and are listed, e.g. in the Appendix of Ref. [27]). The one-particle (level) energies are included in H diag , where the onsite elements of kinetic energy are bare level energies and local crystal fields. The cluster geometry and precise forms of different terms are standard; these terms were introduced in the previous realizations of the three-band d − p model for CuO 2 planes [28] and for other transition metal oxides [27,29]. The diagonal part H diag depends on electron number operators. It takes into account the effects of local crystal fields and the difference of reference orbital energies (here we employ the hole notation, ∆ > 0), between d and p orbitals (for bare orbital energies) where ε d is the average energy of all 3d orbitals, i.e., the reference energy before they split in the crystal field due to the surrounding oxygens. We fix this reference energy for all 3d orbitals to zero, ε d = 0, and use only ∆ = −ε p and thus we write, The first sum is restricted to oxygen sites {i}, while the second one runs over nickel sites {m}. The crystal-field splitting strength vector (∆ z 2 ) describes the splitting between two e g orbitals; here we assume ∆ z 2 = 0 to find the lowest value of ∆ for which the ground state of the NiO 2 plane changes. A finite value of ∆ z 2 = 2.0 eV is assumed in the Appendix.  [33] and depend on three Racah parameters [34]. For systems where only one kind of orbitals are involved, two Kanamori parameters suffice: (i) intraorbital Coulomb interaction U and (ii) Hund's exchange J H . Here for a pair of e g electrons in different orbitals one finds J H = 4B + C (a similar expression can be also written for t 2g electrons), Interorbital Coulomb interactions ∝n iαniβ are expressed in terms of orbital electron density operators for a pair α < β,n iα = σniασ = σd † iασdiασ . Orbital spin oper-atorsˆ S iα ≡ {Ŝ x iα ,Ŝ y iα ,Ŝ z iα } appear in the Hund's exchange term, −2J Hˆ S iα ·ˆ S iβ . In a Mott insulator, charge fluctuations are quenched and electrons localize due to large energy of the lowest multiplet state, (U − 3J H ) t, associated with highspin charge excitation and determining the fundamental Mott gap. For the present Ni + ions, Hund's exchange J H stabilizes high-spin excited states of two holes per site, with spin S = 1.
Local interactions at oxygen ions, H p int , are again rotationally invariant in the orbital space [33] and depend on two Kanamori parameters: (i) intraorbital Coulomb interaction U p and (ii) Hund's exchange J p between a pair of 2p electrons in different orbitals, Interorbital Coulomb interactions ∝n p iµn p iν are expressed in terms of orbital electron density operators for a pair µ < ν,n p iµ = σn p iµσ = σp † iµσpiµσ . Orbital spin operators,ˆ S p iµ ≡ {Ŝ px iµ ,Ŝ py iµ ,Ŝ pz iµ }, appear in the Hund's exchange term, −2J pˆ S iµ ·ˆ S iν . In a Mott insulator, charge fluctuations are quenched and electrons localize due to large energy of the fundamental Mott gap, (U p −3J p H ) t, associated with highspin charge excitation of p 2 in case of p 1 (or p 4 in case of p 5 ) in the ground state configuration. Hund's exchange J p H stabilizes then high-spin S = 1 states at oxygen sites with two holes.
The parameters of the charge-transfer model Eq. (1) used in the present paper are given in Table I. They are compared with the parameters used in Ref. [12]. The only difference is somewhat smaller (pdσ) = 1.3 eV hybridization element which we consider to be more realistic, see Table I.

A. Exact diagonalization
We start with the impurity model and investigate the ground state of (NiO 4 ) 6− cluster [12] which corresponds to the (NiO 3− 2 ) unit in the NiO 2 plane doped by one hole. This cluster has a small basis of orbital states {(z), (z), x, y} and is diagonalized exactly for N ↑ = N ↓ = 1 holes, i.e., one ↑-spin and one ↓-spin hole. Next we calculate the ground states of the model (1) for the (NiO 2 ) 4 cluster with periodic boundary conditions using the Lanczos algorithm [35] for the orbital basis {(z), (z), x, y} in each unit cell. We start with the undoped system, i.e.. 4 holes in this cluster, in the N ↑ = N ↓ = 2-hole configuration. Doped systems contain N ↑ = 3, N ↓ = 2 and N ↑ = N ↓ = 3 holes, respectively.

B. Variational Cluster Approach
We complement our exact diagonalization study by using the variational cluster approach (VCA) [36], which embeds the four-unit-cell cluster into a larger system. This gives us some some effective access to larger lattices, e.g. momentumresolved one-particle spectral densities or density of states. We use (Lanczos) exact diagonalization to obtain the Green's function of a cluster G Cl , consisting of four Ni and six oxygen ions. The cluster self-energy Σ Cl is then extracted from G Cl and inserted into the one-particle Green's function G of the thermodynamic limit. The approximation consists of replacing the (unknown) thermodynamic-limit self-energy by that of the cluster, i.e., by setting Σ ≡ Σ Cl .
According to the self-energy functional theory [37], the optimal cluster self-energy is the one optimizing the thermodynamic grand potential, The parameters { τ } that can be varied to find optimized Σ τ are any one-particle parameters used when solving the small cluster. We focus here on an overall fictitious chemical potential µ to ensure thermodynamic consistency [38] as well as a staggered antiferromagnetic field h favoring Néel order. Neither has been seen to significantly affect results reported here.
Static observables like spin or orbital densities as well as dynamic quantities like the one-particle density of states of k-resolved one-particle spectra are then obtained from the approximated Green's function G. The kinetic energy is thus fully included while interaction effects are truncated to those accessible to the directly solved cluster. A slight difficulty lies in the Coulomb repulsion U pd between Ni and O ions, as this can connect orbitals inside the cluster with others outside. In contrast to inter-cluster hopping, such inter-cluster interactions cannot rigorously be included in the VCA, but have to be truncated or approximated. We have here allowed periodic boundary conditions for these interactions. While this would likely be too crude an approximation if intersite interactions are close to driving an ordered phase (e.g. a charge-density wave) [39][40][41], it has been shown good enough to stabilize uniform charges in the study of a strongly related three-band model for cuprates [42]. Since the oxygen ions have lower occupation in our present work, the approximation can here be expected to be even less critical.

A. Impurity model
First we analyze the weight distribution in the ground state of the doped (NiO 2 ) 2− cluster with 2 holes per unit cell in Ni(3d 10 ) and O(2p 6 ) states, i.e., for the filling by N ↑ = N ↓ = 1 hole. In Fig. 2(a) we see a rapid change from the singlet (S = 0) to triplet (S = 1) state at ∆ c 7.7 eV. This value is somewhat lower than ∆ c = 8.1 eV given in Ref. [12] as we use here a somewhat smaller hybridization, see Table I. But overall we find a very good qualitative agreement with a distinct transition from the ZR singlet state dominated by the b 1 L b1 configuration to the a 1 b 1 state, where both holes are mainly at Ni ion in the high-spin configuration. We observe that this transition is governed by the interplay between the d − p hybridization and Hund's exchange at Ni ion, J H .  Table I: (a) our calculation reproduces qualitatively the transition between two ground state symmetries as in Fig.  4 in Ref. [12]; the transition is seen at somewhat lower ∆ 7.7 eV; (b) somewhat enhanced triplet symmetry for finite U dp = 1 eV.
The ground state for ∆ < ∆ c consists of one hole occupying the |(z) and the oxygen state with the same symmetry constructed out of four bonding {p † iσ |0 } states, i.e., and the ZR singlet at site m is, The weight of the ZR singlet is shown in Fig. 2 and labeled as b 1 L b1 . At ∆ > ∆ c the above weight becomes negligible and the largest weight is found instead for the |ψ m = |a 1 b 1 state which is antisymmetric for spins in {(z), (z)} orbitals and labeled as a 1 b 1 .
Increasing intersite Coulomb repulsion to U dp = 1 eV moves the above transition to a value of ∆ c roughly lowered by U dp , in agreement with the expectation from the mean field approximation. Otherwise the weights of different configurations are almost the same, so one finds that this intersite charge-charge repulsion influences the energies of the two competing states and not their internal structure.

B. Model of NiO2 plane
We next consider a cluster of 2 × 2 unit cells (each containing one Ni and two O ions) with periodic boundary conditions. The undoped system corresponds then to one hole per Ni, i.e., to two holes with spin up and down, N ↑ = N ↓ = 2 holes within a Ni 4 O 8 cluster. Orbital-resolved densities are shown in Fig. 3 and go from an almost even distribution between x 2 − y 2 and 2p orbitals at small ∆ to nearly complete localization on x 2 − y 2 at large ∆. Since there is on average just one hole per one NiO 2 unit cell, there are only very few configurations where two holes may interact by U dp and this parameter does almost not change the ground state at U dp = 0, cf. Figs. 3(a) and 3(b). Consider first the undoped system, i.e., the one with one hole per NiO 2 unit cell. The hole is then predominantly in the orbital of x 2 − y 2 symmetry. The weights of these configurations increase with increasing charge transfer energy ∆ when the holes hybridize more weakly with the surrounding oxygens. The many-hole wave function may be written here as where K σ and L σ are the sets of Ni ions which are occupied by σ-spin electrons, and a (z)(z)(z)(z) and b (z)(z)(z)(z) are the coefficients of two competing states. For large charge-transfer energy ∆ > 6 eV the holes concentrate within the |x 2 − y 2 state at each Ni site in the undoped system. The hole density for the second e g orbital almost vanishes and the high-spin states play no role.
Going to a doped system, there are many more configurations in the real space and a very complex wave function similar to Eq. (11) will arise. Clearly, specifying the hole configurations in real space is unpractical. Therefore we shall focus attention on the leading two-hole configurations, such as these considered for an isolated NiO 4 cluster, see Here the first line is the ZR singlet and the second line stands for the |a 1 b 1 state considered before in Sec. IV A, and {a m , b m } are the probability amplitudes of each of the elemental configurations. Of course, the ground state is more complex and includes also other configurations. Taking |a m | 2 and |b m | 2 as the respective probabilities, we obtain the main local configurations with their weights {b 1 L b1 , a 1 b 1 } displayed in Fig. 4. The finite weights of other configurations displayed as well in Fig. 4 are obtained in a similar way numerically and demonstrate that delocalization of holes over oxygens in the regime of small ∆. Analyzing these weight distributions shows that the character of the wave function changes when the charge transfer energy increases from small to large ∆. Figure 4 refers to the lowest doping achievable on a cluster with four unit cells, namely one additional hole, i.e., x = 1 4 doped hole per Ni site, or N ↑ = 3 and N ↓ = 2. In contrast to the single Ni ion results shown in Fig. 2, we no longer see a sharp jump, but rather a gradual transition. This is due to  Fig. 1 as obtained for increasing ∆, N ↑ = 3, and N ↓ = 2 for: (a) U dp = 0 and (b) U dp = 1.0 eV. The crossing point between the two ground state components with the largest weights, b1L b 1 and a1b1, defines the crossover point between the two types of the ground state. It is moved to a smaller value of ∆ when U dp = 1.0 eV. the transition from a rotationally invariant impurity to a translational invariant lattice, where states of, e.g. b 1 and a 1 symmetry, are allowed to hybridize. In agreement with this wave function, we observe that the most important two configurations for the ground state are: |b 1 L b1 and |a 1 b 1 states. They stand for the ZR singlet state and for the high-spin S = 1 state stabilized by Hund's exchange when ∆ > ∆ c . As observed  Fig. 1 with one or two holes at Ni sites, as obtained for increasing ∆, N ↑ = N ↓ = 3, and for: (a) U dp = 0, and (b) U dp = 1.0 eV. for the impurity model in Fig. 2 and as expected, chargetransfer energies ∆ 7.0 eV lead to a large weight in the ZR-singlet state, while the triplet state localized at Ni dominates at large ∆. Again, non-local Coulomb repulsion U dp shifts this transition precisely in the range of ∆ ≈ 6 − 8 eV expected to apply to nickelates, with U dp > 0 favoring the local triplet state already at lower crystal fields.
As already apparent in the impurity model Fig. 2, other states in addition to the 'pure' ZR-singlet and onsite S = 1 states, more wave functions contribute, especially around the When two holes are added to the Ni 4 O 8 cluster, i.e., N ↑ = N ↓ = 3, one finds locally once again the low-spin ground state with ligand holes |b 1 L b1 for small values of ∆ < 5 eV, see Fig. 5. This regime corresponds to doped cuprates. For even smaller values of ∆ one finds stronger delocalization of holes onto the oxygen orbitals seen by the component of the |d 10 L 2 state. But for ∆ 7 eV as in nickelates [12], this state plays no role and the ground state is qualitatively different-it has the largest weight for the highspin state at Ni sites, |a 1 b 1 . This trend is enhanced by the realistic finite value of U dp = 1.0 eV, cf. Figs. 5(a) and 5(b).
The evolution of the ground state in the doped systems may be qualitatively characterized by the transition from the ground state dominated by the low-spin |b 1 L b1 local states to the more localized at Ni sites high-spin |a 1 b 1 states, as shown in Fig. 6. The delocalization over oxygen orbitals occurs easier in the low doping regime x = 1 4 , while for higher dop-  Table I. The LHB at negative energies is predominantly occupied by holes in x 2 − y 2 orbitals. In the UHB the occupation of z 2 states is enhanced while that of oxygen states is suppressed with increasing U dp .
ing of x = 1 2 the high-spin components dominate already for ∆ < 6.0 eV, taking the realistic finite value of U dp = 1.0 eV.

C. Discussion: comparison between CuO2 and NiO2 plane
To illustrate the nature of the electronic states, we consider the occupied and empty hole states now. As the (NiO 2 ) 3− plane is negatively charged, it is more convenient to use here the hole notation, i.e., the occupied states at low energy in the lower Hubbard band (LHB) are the states occupied by holes. In the undoped systems this means the filling of one hole per each NiO 2 unit. First we consider such a system and show that the undoped NiO 2 plane is insulating, see Fig. 7. Whether or not this corresponds to the real situation in NdNiO 2 is an open question-we suggest that this system is in a poor metallic state after considerable self-doping. It increases the hole concentration in NiO 2 planes beyond one hole per NiO 2 unit cell when Nd ions have a smaller positive change than Nd 3+ . We emphasize that due to the large value of ∆ [12], nickelates are in the Mott-Hubbard regime of the Zaanen-Sawatzky-Allen diagram of correlated insulators [43]. This is also illustrated by two Hubbard bands: the LHB and the upper Hubbard band (UHB), separated by large Mott-Hubbard gap in the undoped system, see Fig. 7. The LHB and the UHB have mainly the states at Ni sites, with little oxygen admixture. The ground state contains predominantly the holes within the x 2 − y 2 orbitals in the LHB, while the z 2 states have large weight for the lower edge of the UHB. This weight is enhanced for U dp = 1.0 eV, cf. Figs. 7(a) and 7(b). Moreover, the character of the lowest UHB states changes with U dp , going from a somewhat ZR-singlet-like band (albeit with z 2 admixture and reduced 2p content) to an almost pure z 2 band.
To illustrate this point, let us first come back the the ZR-singlet states arising in an analogous model applied to cuprates [42]. With a lower value of ∆ = 3.0 eV [42], the CuO 2 plane is in the charge-transfer regime. Hybridization between x 2 − y 2 orbitals with oxygen 2p states is here naturally stronger, so that oxygen content of occupied states rises. In the electronic structure around the Fermi energy, shown in Fig. 8, one find the lower part of the UHB to be dominated by x 2 − y 2 states. Nevertheless, substantial oxygen-p character is present, especially around (π/2, π/2) at the very lowest edges of the spectrum.
When U dp is increased from 0 to 1 eV, the gap increases and orbital contributions to the bands are affected, but overall band shapes and their order remain the same, see Fig. 8(b). For instance, one finds slight changes in the orbital character of the lowest unoccupied states: in addition to x 2 − y 2 and 2p orbitals, even some z 2 weight is now present, especially around (0, π) and (π, 0). As can be seen in Fig. 4, hole doping only induces a very small triplet component for these parameters. The z 2 weight seen in the spectra is thus more easily understood as arising from delocalization of the ligand hole onto neighboring z 2 orbitals. When U dp becomes relevant, it pushes more of the ligand holes into z 2 states. Despite this modification of orbital makeup, lowest states of the UHB are naturally explained in terms of a single ZR-singlet band for both U dp = 0 and 1 eV, while a z 2 band with a minimum at (0, 0) comes at higher energy.
Let us now come back to the nickelate regime with a larger ∆ = 7.0 eV, for which k-resolved spectra are shown in Fig. 9; they correspond to the density of states and the Hubbard subbands of Fig. 7. At first sight, one immediately sees the stronger z 2 character near the lower edge of the UHB, both without or with U dp . The main oxygen bands, on the other hand, are pushed to even higher energies outside the the energy window depicted. The usual ZR-singlet picture does thus not apply, as one can also infer from the sizable triplet component seen for ∆ = 7.0 eV in Fig. 4.
For U dp = 0, however, the density of states in Fig. 7(a) clearly also shows some 2p and x 2 − y 2 weight in the first hole-doping states. The spectrum Fig. 9(a) suggests the corresponding lowest unoccupied band to be related to the ZRsinglet band of the cuprate scenario Fig. 8(b). The shape is similar, with more pronounced minima at (π, 0) and (0, π), but the orbital makeup differs from cuprates, with much re-FIG. 8. The orbital-resolved and k-resolved spectra for the cuprate parameters with ∆ = 3.0 eV along the high symmetry lines in the 2D Brillouin zone for two values of U dp : (a) U dp = 0 and (b) U dp = 1.0 eV; the other parameter as in Table I. Blue and red shading stands for x 2 − y 2 and z 2 orbital states at Ni ions, while green shading corresponds to oxygen states. duced 2p and much increased z 2 character. A description of the lowest states in terms of a ZR-singlet-like band appears still plausible, even though the different orbital makeup and the resulting substantial triplet admixture suggest that effective interaction parameters can be quite different from those describing hole doping in cuprates.
At U dp = 1 eV, finally, both the density of states in Fig. 7(b) and the k-resolved spectrum in Fig. 9(b) show the lowest unoccupied states to be of almost pure z 2 character. While some 2p and x 2 − y 2 character is still present, the system is now more easily understood as a 'classical' Mott-Hubbard insulator, where the gap separates Hubbard subbands built mainly by the orbital states at Ni ions. Accordingly, the onsite triplet dominates over the ZR-singlet in Fig. 4.

V. SUMMARY
In summary, we have investigated the ground state of NiO 2 planes in nickelate superconductors at increasing hole doping. First of all, we find that the local high-spin states |a 1 b 1 are important for ∆ 7.0 eV, as expected for nickelate superconductors. However, there is a crossover transition between the two regimes, dominated by low-spin and high-spin states, FIG. 9. The orbital and k-resolved spectra for the nickelate parameters (Table I) along the high symmetry lines in the 2D Brillouin zone for: (a) U dp = 0.0 eV and (b) U dp = 1.0 eV; the other parameter as in Table I. Blue and red shading stands for x 2 − y 2 and z 2 states at Ni ions, while green shading corresponds to oxygen states. when the charge transfer energy increases.
Second, we would like to point out that U dp is an important parameter to model the situation in doped NiO 2 planes. The Coulomb repulsion between Ni and O sites influences the stability of low-spin ZR states as this interaction is to a large extent missing in the competing high-spin states. The value of ∆ 7.0 eV is in the critical regime where the weight distribution and the nature of the ground state depends in a subtle way on the value of ∆.
Finally, as our third conclusion we wish to point out the increasing admixture of high-spin states at Ni ions with increasing charge-transfer energy ∆ in the ground state of doped nickelates. We have shown that the z 2 orbitals are the next to be occupied by doped holes in case any extra splitting between the z 2 and the x 2 − y 2 orbitals is absent. In this way the holes concentrate at Ni sites and the population of the S = 1 states is increased. In nickelate films the orbitals which are preferentially occupied are instead xy, but they lead to the same global result that high-spin states arise locally. In the Appendix we present arguments how the upper Hubbard band may change when the xy orbitals are filled instead by holes rather than the z 2 ones.
In the absence of electrostatic crystal field, the hybridization with 2p orbitals indeed slightly favors the z 2 orbitals over the t 2g states, but the difference is not large. Electro-static crystal fields due to the planar geometry and missing apical oxygens can thus easily tip the balance. However, the emerging scenario remains equivalent: At ∆ ≈ 7.0 eV and U dp = 0, hole doping still leads to large singlet weight, while triplet states (now involving the xy orbitals) dominate for U dp = 1 eV.
In general, we find that hole doping distributes over the entire system, in contrast to the onset of local S = 1 states [22]. This issue is however still open and should be investigated further. At the moment we can conclude that the contribution of high-spin states to the electronic structure arises in any case, independently of which Ni orbitals are occupied by a second hole in doped systems.

ACKNOWLEDGMENTS
We thank Andres Greco, Krzysztof Rościszewski, and George A. Sawatzky for very insightful discussions. We thank Ali Alavi for informing us about the content of Ref. [26] prior to publication. This research was supported in part by the National Science There is some discussion concerning the order of the d orbitals, apart from the commonly accepted view that holes in the 'undoped' reference system are preferentially found in x 2 − y 2 orbitals. Due to the planar cluster geometry (see Fig. 1), one can argue that the charge repulsion coming from the oxygen ions favors the electron occupation of the 3z 2 − r 2 orbital and that the xy orbital should come below the x 2 − y 2 orbital, as indeed suggested by the crystal-field splittings [25]. Such arguments were also given for instance in Refs. [44,45] and are supported by quantum chemistry calculations of hole levels [26].
On the other hand, one can argue that effective crystal-field splitting is often determined by hybridization rather than electrostatic forces. In this picture, putting a hole into the 3z 2 −r 2 orbital is favorable because it can then delocalize better than the hole which resides in the xy orbital. Such effects resulting in the level order opposite to the one expected from electrostatic charge have, e.g. been observed in iridates with their extended 5d wave functions [46]. In a (b) direction, the corresponding hopping parameters are t pd /2 √ 3/2 0.866 eV for the hopping between the 3z 2 − r 2 orbital and the p x (p y ) orbital vs. (pdπ) 0.75 eV for hopping between xy and p y (p y ). Based on the analysis of Ref. [47], where hybridization effects are found to dominate over point-charge effects, it has been argued [12] that the 3z 2 −r 2 orbital should be considered preferentially.
A closer look at the results reported in Ref. [47] indicates, however, that the hybridization effects are not very large and in fact the strongest in the 'charge-transfer' limit of a smaller crystal field separating oxygen and transition metal ion. When the hole becomes more localized onto the transition metal ion, impact of hybridization with oxygen is reduced, so that the energy difference between xy and 3z 2 − r 2 shrinks. We thus investigate the aspect of level ordering in this appendix.
To do so, we extend the model to include three 3d orbitals {x 2 − y 2 , 3z 2 − r 2 , xy} per each Ni ion as well as two (p x and p y ) at each oxygen, and study it using the VCA. We then have to truncate the Hilbert space of the four-unit-cell cluster, but wish to point out that a three-unit-cell cluster with the full Hilbert space gave equivalent results. Since the x 2 −y 2 orbital was nearly always found to be half filled, we restrict possible states to these with at least two holes in these orbitals in the cluster shown in Fig. 1. In other words, we allow at most two of the four x 2 − y 2 orbitals to be completely filled.
Parameters referring to the orbitals included in the model with a smaller basis set of the main text, as summarized in Tab. I, are kept here. In addition, Hund's exchange coupling J p H = 0.8 eV on oxygen ions as well as corresponding interorbital U p = U p − 2J p H was introduced. Additional hoppings were chosen following Ref. [47]: t pp = −0.35 eV between nearest neighbor parallel oxygen orbitals and t xy = t pd /2 = 0.65 eV connecting the xy orbital to oxygen. For the signs resulting from orbital phases, see Ref. [47].
Without any static point-charge fields, the xy orbital is indeed found to be above the 3z 2 − r 2 orbital (in the hole notation), so that one would encounter the latter first when hole doping the compound. Below we shall follow here the first principle calculations which give the xy orbitals as the energetically closest ones to x 2 − y 2 [25]. If this is not assumed, for U dp = 0, holes enter first a band composed of a mixture of x 2 − y 2 , 3z 2 − r 2 , and 2p orbitals, see Fig. 10(a). 3z 2 − r 2 states are mixed into this band, as seen before in Fig. 9(a), while xy states come at slightly higher energies and do not mix noticeably with the lowest hole band.
However, the splitting between the 3d states, 3z 2 − r 2 and xy, is very small indeed and moderate explicit crystal fields with the sum going over all Ni ions and n m,z 2 referring to the hole density in the 3z 2 − r 2 orbital, could easily reverse their order. Such a situation is shown in Fig. 10(b) for a moderately large ∆ z 2 = 2.0 eV, which is in line with recent quantumchemistry results [26]. Nevertheless, conclusions from the main text remain valid: for U dp = 0, lowest hole-doping states are composed on x 2 − y 2 and 2p states. In fact, their lower 3z 2 −r 2 -population and nearly absent xy character makes this band rather more similar to the ZR-singlet band known from cuprates.
Once U dp = 1 eV, Mott-Hubbard gap increases (Fig. 11) and doped holes hardly go into the x 2 − y 2 orbital, regardless of whether ∆ z 2 = 0 or ∆ z 2 = 2 eV. For ∆ z 2 = 0, the lowest hole-doping states are now of almost exclusively 3z 2 − r 2 character, see Fig. 11(a). Once ∆ z 2 = 2.0 eV, the crystal field FIG. 10. k-resolved one-particle spectra obtained for the 7-orbital model which includes three {x 2 − y 2 , 3z 2 − r 2 , xy} orbitals per Ni site and two {px, py} orbitals at oxygen sites in NiO2 unit cell of Fig. 1. In (a) there is no explicit crystal field between x 2 − y 2 and 3z 2 − r 2 orbitals (5), so that the xy orbital lies slightly above the lowest 3z 2 − r 2 states. In (b) the crystal field ∆ z 2 = 2.0 eV, see Eq. (5), rises the 3z 2 − r 2 orbital, but the basic scenario of holes entering a ZR-singlet band remains intact. The parameters as in Tab. I, and ∆ = 7.0 eV, U pd = 0. Yellow shading refers to the xy, red to 3z 2 − r 2 , and blue to x 2 − y 2 orbitals. pushes the 3z 2 − r 2 states towards higher energies and holes enter a nearly pure xy band instead, see Fig. 11(b).
We thus conclude that the xy orbital is likely to be relevant to doped NiO 2 planes, but that this does not affect our results: for U pd = 0, a band with robust x 2 − y 2 character orbital hosts the lowest hole-doping states regardless of which orbital comes next, see Fig. 10. For U pd = 1.0 eV, conversely, holes enter preferentially either xy or 3z 2 − r 2 orbital, depending on their relative crystal field, see Fig. 11, and high-spin S = 1 states form locally.
Since the xy orbital bonds to the orbitals orthogonal to those coupling with both the x 2 − y 2 and z 2 orbitals, x 2 − y 2 states see less oxygen-mediated hybridization with xy that with z 2 states. The ZR-singlet-like band thus remains more clearly separated from the xy band than from the 3z 2 − r 2 states. In the full system, however, the 'other' orbitals in addition to x 2 − y 2 have been shown to hybridize strongly with rare-earth states and also with each other [48]. This clear separation of xy and x 2 − y 2 orbitals may thus be an artifact of our model, while the stronger mixing of x 2 − y 2 -states with a second band -as discussed in the main text-may be more FIG. 11. As in Fig. 10, but for U pd = 1.0 eV. In (a) without crystal field, the lowest hole-doped states are now 3z 2 − r 2 , while they have xy character when ∆ z 2 = 2.0 eV (5) in (b). The color convention is the same as in Fig. 10. realistic [11].