Magnetic exchange coupling in cuprate-analog $d^9$ nickelates

Motivated by the recent discovery of superconductivity in doped NdNiO$_2$, we study the magnetic exchange interaction $J$ in layered $d^9$ nickelates from first principles. The mother compounds of the high-$T_{\rm c}$ cuprates belong to the charge-transfer regime in the Zaanen-Sawatzky-Allen diagram and have $J$ larger than 100 meV. While this feature makes the cuprates very different from other transition metal oxides, it is of great interest whether layered $d^9$ nickelates can also have such a large $J$. However, one complexity is that NdNiO$_2$ is not a Mott insulator due to carrier doping from the block layer. To compare the cuprates and $d^9$ nickelates on an equal basis, we study RbCa$_2$NiO$_3$ and $A_{2}$NiO$_{2}$Br$_2$ ($A$: a cation with the valence of $2.5+$), which are recently designed theoretically by block-layer engineering. These nickelates are free from the self-doping effect and belong to the Mott-Hubbard regime. We show that these nickelates share a common thread with the high-$T_{\rm c}$ cuprates in that they also have a significant exchange interaction $J$ as large as about 100 meV.


I. INTRODUCTION
The discovery of superconductivity in doped nickel oxides Nd 0.8 Sr 0.2 NiO 2 [1,2] has attracted intensive interests both in experiment [3][4][5][6][7][8][9][10][11][12][13] and theory [3,, because the nickelate might be an analog of the well known high-T c superconductor, cuprates. Recently, the doping dependence has been explored both theoretically [38] and experimentally [8,9], and the presence of the superconducting dome has been confirmed [8,9]. The maximum superconducting transition temperature T c is about 15 K, not very high compared to that of the high-T c cuprates. However, because the Bardeen-Cooper-Schrieffer (BCS) phonon mechanism cannot explain the observed T c [17], the superconducting mechanism is most likely unconventional, in which the electron correlations play an important role. A recent observation of d-wave like superconducting gap also supports this scenario [12]. Here, a natural question arises: is there any possibility to realize T c as high as the cuprates in nickelates?
In the cuprates, the superconductivity emerges by doping carriers into the antiferromagnetic Mott insulator. In the theoretical Zaanen-Sawatzky-Allen diagram [57], the cuprates belong to the charge-transfer type in which the energy difference between the oxygen 2p and cupper 3d orbitals (∆ dp ) controls the metal-insulator transition. Since ∆ dp in the cuprates is small among transition metal oxides, the magnetic exchange coupling J takes one of the largest values known (∼130 meV) [58]. This large value of J is certainly a characteristic feature of the cuprates, which makes the cuprates very different from other transition metal oxides.
On the other hand, in the case of the nickelate NdNiO 2 , ∆ dp is larger than that of the cuprates [59]. Thus, naively, we expect smaller J for nickelates. Indeed, a recent experimental estimate using the Raman spec- * yusuke.nomura@riken.jp troscopy gives J = 25 meV [6]. However, it should be noted that the origin of small J in NdNiO 2 may be ascribed to another notable difference from the cuprates, namely, NdNiO 2 is not a Mott-insulator due to the selfdoping effect. In NdNiO 2 , orbitals in the Nd layer form extra Fermi pockets on top of the large Fermi surface formed by the Ni 3d x 2 −y 2 orbital, and the Ni 3d x 2 −y 2 orbital is hole-doped, i.e., the filling of the Ni 3d orbitals deviates from d 9 [59]. The self-doping naturally explains the absence of Mott-insulating behavior in NdNiO 2 . Although it has been shown that the Ni 3d x 2 −y 2 orbital forms a two-dimensional strongly-correlated system [17], J at the d 9 configuration with half-filled d x 2 −y 2 orbital is masked by the self-doping. The experimental estimate should be understood as the J value including the effect of the self-doping, not the J value at the ideal d 9 configuration. One of the reasons for the controversy in theory about the size of J [21][22][23][24][25][26][27][28][29][30] is ascribed to the ambiguity in calculating J (whether we calculate J at d 9 filling or J including the self-doping effect). In any case, it is a non-trivial problem whether we can justify the mapping onto a simple spin model to understand the property of NdNiO 2 . This fact makes NdNiO 2 an imperfect analog of the cuprates.
Recently, there is a proposal to design cuprate-analog nickelates without the complication of the self-doping [18] [60]. Since NdNiO 2 is a layered material, one can systematically propose nickelate family materials by changing the composition of the "block-layer" [61] between NiO 2 layers. Proposed dynamically stable nickelates have smaller Fermi pockets of the block-layer orbitals than NdNiO 2 . In some materials, the self-doping is completely suppressed, and the ideal d 9 system with halffilled 3d x 2 −y 2 orbital is realized. An ab initio estimate of Hubbard U using the constrained random-phase approximation (cRPA) [62] shows that the correlation strength U/t (t: nearest-neighbor hopping) is comparable to that of cuprates [18]. Therefore, once such nickelates are synthesized, the mother compounds will be a Mott insulator similarly to the cuprates, and the effective model becomes the Heisenberg model, which gets rid of the ambiguity in calculating J.
In this paper, we study the strength of J in the two ideal d 9 nickelates, which are free from the self-doping (see Sec II for the details of the materials). We estimate the J value by the following three methods [63]. First, we start from a single-orbital Hubbard model derived in Ref. 18 and then evaluate J by the expansion in terms of t/U . Second, we perform an energy mapping between the classical Heisenberg model and the total energy of different magnetic configurations calculated by the LDA+U (LDA: local density approximation) method. Third, we employ a scheme based on the so-called local force theorem. Hereafter, we simply call these three methods "strong-coupling expansion", "energy mapping method", and "local force approach", respectively. We show that the three independent estimates show reasonable agreement and conclude that the d 9 nickelates have sizeable J (about 100 meV) comparable to that of the cuprates. Therefore, the proposed d 9 nickelates provide an interesting playground to explore the cuprate-analog high-T c superconductivity.
The paper is organized as follows. In Sec. II, we introduce two ideal d 9 nickelates, RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A: a cation with the valence of 2.5+) and discuss the advantage over NdNiO 2 . In Sec. III, we explain the three methods employed in the present study, and we show the results in Sec. IV. Section V is devoted to the summary.

II. MATERIALS: d 9 NICKELATES
In Ref. [18], various layered nickelates have been systematically proposed. They are classified into "1213", "1214", "H 2 ", and "G" families, depending on the composition and the type of the block-layer [61]. Among the four families, the compounds without the self-doping exist in the 1213 and G families. We here take RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A: a cation with the valence of 2.5+) for a representative of the ideal d 9 nickelates belonging to 1213 and G families, respectively (see Figs. 1(a) and (c) for the crystal structure). In the following, we employ Ba 0.5 La 0.5 as A. The phonon calculations have shown that both RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ) are dynamically stable [18]. We take the crystal structure optimized in Ref. [18], and perform density-functional theory (DFT) calculations. Figs. 1(b) and 1(d) show the paramagnetic DFT band structure for RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ), respectively. As is shown in Ref. [18], only the Ni 3d x 2 −y 2 orbital crosses the Fermi level. As far as the topology of the band structure is concerned, these systems are more similar to the cuprates than NdNiO 2 .
The advantages of studying these nickelates rather than NdNiO 2 are as follows. First, it is still controversial whether the role of Nd-layer (block-layer) orbitals is essential or not. If the hybridization between Ni 3d and Nd-layer orbitals is substantial, the Nd-layer orbitals are not only a charge reservoir, but they might give Kondolike physics [2,3,24,33,43]. In the cases of the d 9 nickelates, RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ), the block-layer orbitals do not show up at the Fermi level, and this controversy can be avoided. We can also exclude the possible role of the 4f orbitals with localized moments proposed in Refs. [44,45].
Another controversial issue for NdNiO 2 is to which orbitals the doped holes go (d 9 L vs. d 8 , where L denotes a hole in a ligand oxygen). In the case of the cuprates (charge-transfer insulator), the holes are doped into the oxygen 2p orbitals. On the other hand, the nickelates have larger ∆ dp and are classified as Mott-Hubbard type [3,6,10,17,21,23]. Because there are nonzero hybridization between Ni 3d x 2 −y 2 and O 2p orbitals, some of the holes should be doped into oxygen 2p orbitals [28,37]. However, the amount should be smaller than that of the cuprates.
When the system is Mott-Hubbard type and the holes mainly reside in the Ni 3d orbitals, another issue arises: which model is more appropriate, the single-orbital or multi-orbital model? In other words, whether the doped d 8 configuration favors high-spin state or low-spin state. If the crystal field splitting between Ni 3d x 2 −y 2 and the other 3d orbitals is much larger than the Hund's coupling, holes stay within the Ni 3d x 2 −y 2 orbital, and the singleorbital model is justified. On this issue, several studies insist that Ni 3d multi-orbital nature cannot be ignored [21,31,32,35,36,[39][40][41][42]. To resolve this issue, we certainly need more experimental evidences. In the cases of RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ), compared to NdNiO 2 , the Ni 3d x 2 −y 2 orbital is more isolated in energy space from the other 3d orbitals [see Figs. 1(b) and 1(d)]: In the case of NdNiO 2 , due to the dispersion along the k z direction, the position of the Ni d 3z 2 −r 2 band becomes close to the Fermi level on the k z = π/c plane; however, such k z dependence is much weaker in RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ). Considering also the above-mentioned absence of the complication from the self-doping, in this study, we adopt the single-orbital Hubbard model as a minimal model for RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ). In the absence of the carrier doping, we can further map onto the spin model with the exchange coupling J.  → (π/a, π/a, 0) → (0, 0, 0) → (0, 0, π/c) → (π/a, 0, π/c) → (π/a, π/a, π/c) → (0, 0, π/c) (The symbols are different because the primitive cell of RbCa2NiO3 and A2NiO2Br2 is tetragonal and bace-centered tetragonal, respectively).

A. Strong-coupling expansion
When the single-orbital Hubbard model is a good description, the magnetic interactions can be obtained by strong-coupling perturbation expansion. The superexchange interaction J s (with t 4 -order correction term) and cyclic ring-exchange interaction J c are given by J s = 4t 2 /U − 24t 4 /U 3 and J c = 80t 4 /U 3 , respectively [64][65][66]. If we effectively take into account the effect of the ringexchange interaction in the nearest-neighbor interaction J, the J value becomes with S = 1/2.

B. Energy mapping method
Within the LDA+U [67][68][69][70], we perform the magnetic calculations. Here, U is introduced into the Ni 3d orbital subspace. We employ 2 × 2 × 1 supercell consisting of four conventional cells. We simulate two different magnetic solutions: one is Néel type [(π/a, π/a, 0) antiferromagnetic order] and the other is stripe type [(π/a, 0, 0) antiferromagnetic order]. We calculate the energy difference ∆E between the two antiferromagnetic solutions. When we assume the two-dimensional classical spin-1/2 Heisenberg model up to next-nearest-neighbor magnetic interaction J , ∆E per formula unit is given by ∆E = J/2−J J/2. We estimate J with this equation.

C. Local force approach
Based on the Néel-type solutions of the LDA+U calculations, we estimate J and J using the local force theorem [63]. Here, we employ the so-called Lichtenstein formula, which is recently developed in the low-energy Hamiltonian with the Wannier orbitals [71][72][73], given by, where i, j represent the atomic sites, and ω n = (2n+1)πT denotes the Matsubara frequency. Here, we set P = 0 (1) when the spins at i and j sites are aligned parallel (antiparallel) to each other.
The energy comparison between the Néel-and stripetype antiferromagnetic solutions is performed using 9 × 9×7 and 9×9×3 k-mesh for RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ), respectively. We treat Ba 0.5 La 0.5 by the virtual crystal approximation. The energy cutoff is set to be 100 Ry for the Kohn-Sham wave functions, and 400 Ry for the electron charge density.
For the estimate of J based on the local force approach, we first construct the maximally localized Wannier functions [79,80] for the Néel-type antiferromagnetic band structure using RESPACK [81,82]. For RbCa 2 NiO 3 , we use 5 × 5 × 5 k-mesh for the construction of Wannier orbitals. We put Ni d, O p, Ca d, and interstitial-s (located at the interstitial positions surrounded by Ni + , Ca 2+ , and Rb + cations) projections. The interstitial orbitals are stabilized because they feel attractions from the surrounded cations [17]. Then, we obtain 104 orbital (per spin) tight-binding Hamiltonian. For A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ), we employ 5 × 5 × 3 k-mesh for constructing Wannier orbitals. We derive 232 orbital (per spin) tight-binding Hamiltonian using the projections of Ni d, O p, Br p, A d, and interstitial-s (located at the interstitial positions surrounded by Ni + , A 2.5+ , and Br − ions) orbitals.
In the calculation of Eq. (2), we employ 16 × 16 × 16 kmesh and set the inverse temperature β = 200 eV −1 for both cases. We use the intermediate representation basis for the Matsubara frequency summation [83][84][85], and set the cutoff parameter Λ = 10 5 , which is sufficiently larger than W β where W is the band width.

IV. J IN d 9 NICKELATES
In the previous study [18], the effective single-orbital Hamiltonians for RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ) are constructed using maximally-localized Wannier functions [79,80] and cRPA [62]. The derived nearest-neighbor hopping and Hubbard parameters are t = −0.352 eV, U = 3.347 eV for RbCa 2 NiO 3 , and t = −0.337 eV, U = 3.586 eV for A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ). Then, the strong-coupling expansion described in Sec. III A gives J = 122 meV and J = 109 meV for RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ), respectively (see Appendix A for the estimate from threeorbital d-p model). Figures 2(a) and 2(b) show the band structure calculated by the LDA+U method for the Néel-type antiferromagnetic state. While the Hubbard U in the singleorbital Hubbard model is the Coulomb repulsion between the Wannier orbitals made from the Ni 3d x 2 −y 2 orbital with O 2p tails, the U interaction in the LDA+U calculation is the Coulomb repulsion between the Ni 3d orbitals. To make the difference clearer, we call U in the LDA+U calculation U . In Figs. 2(a) and 2(b), we have used U = 3 eV.
In contrast to the case of the LDA+U calculation for NdNiO 2 , where the system stays metal even in the presence of antiferromagnetic order [14,23,25], both systems become insulating. The top of the valence band has mainly Ni 3d character, in agreement with the classification into the Mott-Hubbard type insulator. We see that both systems are insulating even at smaller U (= 1 eV). For all the U region we studied (1-5 eV), there exists well defined spin-1/2 Ni spin moment. The results suggest that, if these d 9 nickelates are synthesized, they become antiferromagnetic Mott insulator as in the cuprates.
Figure 2(c) shows the energy difference ∆E per formula unit between the Néel-and stripe-type antiferromagnetic solutions. ∆E decreases as U increases, which is a natural behavior given that ∆E is governed by J and the origin of J is the superexchange interaction.
In Figs. 2(a) and 2(b), the band dispersions obtained by the Wannier tight-binding Hamiltonian, which are used in the local-force approach, are also shown. The Wannier bands well reproduce the LDA+U magnetic band dispersions.
From ∆E in Fig. 2(c), we perform the order estimate of J by the energy mapping method with assuming J /J = 0.05 (Sec. III B) [86]. Then J is given by J = ∆E/0.45. We also estimate J using the local force approach (Sec. III C).
These results on top of the J value estimated by the strong-coupling expansion (see above) are summarized in Figs Although the energy mapping method and local force approach are based on the same LDA+U calculations, we see that there is a discrepancy between the two results at small U values (although the difference is no more than two times). It should be noted that the former method sees the global change of the energy between the completely different magnetic patterns, whereas the latter approach only sees the local landscape around the Néel-type solutions. For larger U , the agreement between these two results becomes better as is expected: The system can be mapped to the classical spin model with a constant J regardless of the assumed magnetic structure in the local force approach.
Overall, all the three estimates of J lie within 60-140 meV, and we conclude that the d 9 nickelates have a sizable J of the order of 100 meV. The size of J is not far smaller than the value of the cuprates (∼130 meV [58]). Therefore, although these d 9 nickelates belong to the Mott-Hubbard type in the Zaanen-Sawatzky-Allen diagram [57], they share an interesting common feature with the high T c cuprates in that the exchange coupling is significantly large.

V. SUMMARY
One of the remarkable features of the high T c cuprates is the large exchange coupling J, whose size is as large as 130 meV [58]. In the present study, we have evaluated the size of J for d 9 nickelates from first principles. While the cuprates having small ∆ dp belong to the charge-transfer type in the Zaanen-Sawatzky-Allen diagram [57], nickelates with larger ∆ dp belong to the Mott-Hubbard type. To answer how large J can be ex-pected in the Mott-Hubbard insulating d 9 nickelates, we studied RbCa 2 NiO 3 and A 2 NiO 2 Br 2 (A = Ba 0.5 La 0.5 ), which was recently proposed theoretically and shown to be free from the self-doping in Ref. 18. By means of the strong-coupling expansion, energy mapping method, and local force approach, we found that J for these nickelates is as large as 100 meV, which is not far smaller than that of the cuprates. This result suggests that the d 9 nickelates and cuprates share a notable common feature in the Mott insulating phase, although the former and latter belong to the Mott-Hubbard and charge-transfer regime, respectively.
Finally, we note that the proposed d 9 nickelates might give rare examples of realizing the square-lattice Hubbard model with sizeable J in real materials. Recent numerical studies show that the phase diagram of the doped Hubbard model is under severe competition between the stripe state with charge/spin modulation and d-wave superconductivity [87][88][89][90]. Therefore, once synthesized, the d 9 nickelates will serve as a valuable testbed system to understand the superconductivity in the Hubbard-like model. They are also an important reference to understand the superconducting mechanism in the cuprates, because they would tell us whether the charge-transfer nature in the cuprates is essential in the high-T c superconductivity or not. In the main text, we estimate J by the strong-coupling expansion starting from the single-band Hubbard model. Here, we show that the J value is also on the order of 100 meV even when we perform the strong-coupling expansion based on the so-called d-p model consisting of Ni 3d x 2 −y 2 and two O 2p orbitals. In the strong-coupling expansion of the d-p model for the filling of one hole per unit cell, J is given by J = 4t 4 dp ∆ 2 dp U dd + 4t 4 dp ∆ 2 dp (∆ dp + U pp /2) , where t dp is the hopping between Ni 3d x 2 −y 2 and O 2p orbitals, U dd and U pp are the onsite Coulomb repulsion for Ni 3d x 2 −y 2 and O 2p orbitals, respectively. Using the RESPACK [81,82] based on the cRPA method [62,91] combined with the maximally-localized Wannier functions [79,80], we constructed three-orbital d-p model from first principles. We consider the double counting effect in the Hartree term and ∆ dp is given by ∆ dp = ∆ DFT dp + U dd n DFT d /2 − U pp n DFT p /2, where the superscript DFT stands for the DFT value, and n is the hole occupation.