Effective Bi-Layer Model Hamiltonian and Density-Matrix Renormalization Group Study for the High-T_c Superconductivity in La₃Ni₂O₇ under High Pressure

The successfully synthesized infinite-layer nickelate superconductors^[1–7] have provided another platform for research of the microscopic origin of unconventional superconductivity. Similar to high-T_c cuprates, the infinite-layer nickelates with nominal Ni⁺ (3d⁹) cations contain 3d_x²–y² orbital degrees of freedom on a quasi-two-dimensional Ni–O square lattice.^[8,9] However, the electronic states of oxygen 2p orbitals are far below the Fermi level and have a much reduced 3d–2p mixing due to the larger separation of their on-site energies, and their low-energy electronic structures are more likely to fall into the Mott–Hubbard than the charge-transfer regime with a super-exchange energy being at least an order of magnitude smaller than in cuprates.^[10] As a consequence, the parent infinite-layer nickelates may be modeled as a self-doped Mott insulator with two types of charge carriers, and the low-temperature upturn of the electric resistivity^[1,11] arises from the magnetic spin scattering between low-density conduction electrons from the rare earths and localized Ni-3d_x²–y² magnetic moments.^[12,13] Upon Sr doping, the superconducting T_c is just around 9–15 K in the Nd_0.8Sr_0.2NiO₂ thin films,^[1] and the maximum T_c of 31 K has been achieved in Pr_0.82Sr_0.18NiO₂ films under high pressure.^[14] So far, no bulk crystals can be synthesized and show superconductivity.

Recently, it has been reported that the superconductivity with possible T_c ≈ 80 K is observed in the single crystal of La₃Ni₂O₇ with pressure between 14.0 and 43.5 GPa using high-pressure resistance and mutual inductive magnetic susceptibility measurements.^[15] Density functional theory (DFT) calculations for the two nearest intra-layer Ni cations in a bilayer Ruddlesden–Popper (RP) phase^[15,16] have suggested that both the 3d_x²–y² and 3d_z² orbitals of Ni cations strongly mix with oxygen 2p orbitals. The 3d_z² orbitals via the apical oxygen usually have a large inter-layer coupling due to the quantum confinement of the NiO₂ bilayer in the structure, and the resulting energy splitting of Ni cations can dramatically change the distribution of the averaged valence state of Ni^2.5+. The numerical results further indicated that the superconductivity emerges coincidently with the metallization of the σ-bonding bands under the Fermi level, consisting of the 3d_z² orbitals with the apical oxygen connecting Ni–O bilayers.^[17,18] These distinct features are important clues for the high-T_c superconductivity in this RP double-layered perovskite nickelates, which are different from the infinite-layer nickelate superconductors.

In this work, based on the DFT electronic structure,^[15,16] we propose an effective bi-layer model Hamiltonian including both 3d_z² and 3d_x²–y² orbital electrons of the nickel cations, which are different from the single orbital bi-layer Hubbard model.^[19–27] The main feature of the model is that the 3d_z² electrons form inter-layer σ-bonding and anti-bonding bands via the apical oxygen anions between two layers, while the 3d_x²–y² electrons hybridize with the 3d_z² electrons within each NiO₂ plane. Due to their special spatial symmetries of two e_g orbitals, the intra-layer hopping of the 3d_z² orbital electrons and the inter-layer hopping of the 3d_x²–y² orbital electrons are very small, and can be neglected. The chemical potential difference of these two orbital electrons ensures that the 3d_z² orbitals are close to half-filling and the 3d_x²–y² orbitals are near quarter-filling. The strong on-site Hubbard repulsion gives rise to an inter-layer antiferromagnetic super-exchange of the 3d_z² orbital electrons J. Applying pressure can increase the coupling strength J and self-dope additional holes on the 3d_z² orbitals with the same amount of electrons doped on the 3d_x²–y² orbitals. By performing numerical density-matrix renormalization group (DMRG) calculations on a minimum one-dimensional setup and focusing on the large J and small doping of 3d_z² orbital limit, we observe a charge-density wave (CDW) for both 3d_z² and 3d_x²–y² electrons but with different wavelengths. We also find instability of superconductivity for both orbitals from the equal-time spin singlet pair–pair correlations. We attribute the pairing on 3d_z² orbital to the formation of inter-layer singlet pairs, and the pairing on 3d_x²–y² orbitals from the hybridization of the two orbitals.

Effective Model Hamiltonian. Let us focus on the bi-layer RP bulk single crystals of La₃Ni₂O₇ [see Figs. 1(a) and 1(b)]. A simple electron count gives Ni^2.5+, i.e., 3d^7.5 state for both Ni cations, and the previous experiments indicated that La₃Ni₂O₇ is a paramagnetic metal. Ni^2.5+ is usually believed to be given by mixed valence states of Ni²⁺ (3d⁸) and Ni³⁺ (3d⁷), corresponding to the half-filled of both 3d_z² and 3d_x²–y² orbitals and singly occupied 3d_z² with empty 3d_x²–y² orbitals, respectively. With a bilayer RP phase, two 3d_z² orbitals via apical oxygen anions usually have a large inter-layer coupling due to the quantum confinement of the NiO₂ bilayer in the structure, which causes a large energy splitting of Ni cations, see Fig. 1(c). Then the distribution of the averaged valence state of Ni^2.5+ can be dramatically changed.

Zoom In Zoom Out Reset image size

We first consider an isolated Ni–O–Ni three-site system, which includes only the half-filled 3d_z² orbitals [see Fig. 1(c)],

$$\begin{eqnarray}\begin{array}{ll}H= & v\displaystyle \sum _{\sigma }({d}_{2\sigma }^{\dagger }{p}_{\sigma }-{d}_{1\sigma }^{\dagger }{p}_{\sigma }+{\rm{h}}.{\rm{c}}.)+{\epsilon }_{p}{n}_{p}\\ & +{\epsilon }_{d}({n}_{d1}+{n}_{d2})+U{n}_{d1\uparrow }{n}_{d1\downarrow }+U{n}_{d2\uparrow }{n}_{d2\downarrow },\end{array}\end{eqnarray} \tag{ 1 }$$

where v is the hybridization between the 3d_z² orbital of nickel and the p_z orbital of apical oxygen, _d (_p ) is the energy for the 3d_z² (p_z ) orbital, and U is the on-site Coulomb repulsion for the 3d_z² orbital of nickel. Notice that the signs for hybridization are different for the two 3d_z² orbitals. When U ≫ v,_p ,_d , the lower energy configurations are shown in (i) and (ii) in Fig. 1(c), where the oxygen p_z orbital is doubly occupied and the 3d_z² orbital is half-filled. The state of these two 3d_z² orbitals occupied by parallel electrons has a higher energy. By considering the virtual transition to higher energy states, we can derive an anti-ferromagnetic super-exchanges between these two 3d_z² orbitals as

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}=J{{\boldsymbol{S}}}_{1}\cdot {{\boldsymbol{S}}}_{2},\end{eqnarray} \tag{ 2 }$$

with

$$\begin{eqnarray}&&J=\displaystyle \frac{4{v}^{4}}{{(U+\varDelta )}^{2}}\left(\displaystyle \frac{1}{U}+\displaystyle \frac{1}{U+\varDelta }\right),\,\,\varDelta ={\epsilon }_{d}-{\epsilon }_{p}.\end{eqnarray} \tag{ 3 }$$

Then a simplified bi-layer model Hamiltonian of a square planar coordinated Ni cations consisting of the 3d_z² and 3d_x²–y² orbitals characterizes the effective low-energy physics,

$$\begin{eqnarray}\begin{array}{ll}H= & -{t}_{{x}^{2}-{y}^{2}}\displaystyle \sum _{\langle ij\rangle ,\sigma ,a=1,2}({c}_{a,i,\sigma }^{\dagger }{c}_{a,j,\sigma }+{\rm{h}}.{\rm{c}}.)\\ & -{\mu }_{{x}^{2}-{y}^{2}}\displaystyle \sum _{i,a=1,2}{n}_{a,i}^{c}+U\displaystyle \sum _{i,a=1,2}{n}_{a,i\uparrow }^{c}{n}_{a,i\downarrow }^{c}\\ & -{t}_{{x}^{2}-{y}^{2},{z}^{2}}\displaystyle \sum _{i,\sigma ,a=1,2}({d}_{a,i,\sigma }^{\dagger }{\mathop{c}\limits^{\sim }}_{a,i,\sigma }+{\rm{h}}.{\rm{c}}.)\\ & -{t}_{{z}^{2}}\displaystyle \sum _{i,\sigma }({d}_{1,i,\sigma }^{\dagger }{d}_{2,i,\sigma }+{\rm{h}}.{\rm{c}}.)\\ & +J\displaystyle \sum _{i}{{\boldsymbol{S}}}_{1,i}^{d}\cdot {{\boldsymbol{S}}}_{2,i}^{d}-{\mu }_{{z}^{2}}\displaystyle \sum _{i,a}{n}_{a,i}^{d},\end{array}\end{eqnarray} \tag{ 4 }$$

with ${c}_{a,i,\sigma }^{\dagger }$ (${d}_{a,i,\sigma }^{\dagger }$) creates a 3d_x²–y² (3d_z² ) electron at the i-th site for the layer a = 1, 2, ${\mathop{c}\limits^{\sim }}_{a,i,\sigma }=({c}_{a,i+x,\sigma }+{c}_{a,i-x,\sigma }-{c}_{a,i+y,\sigma }-{c}_{a,i-y,\sigma })/2$, t_z² is the hopping for the 3d_z² electron between the two layers and is set as the energy unit, t_x²–y² is the hopping for the 3d_x²–y² electron in each layer, t_{x²–y²,z²} is the intra-layer hybridization of the 3d_x²–y² electron and 3d_z² electrons between the nearest neighbors and the signs for the vertical and horizontal directions are opposite,^[15] and μ_x²–y² and μ_z² are the chemical potentials for the 3d_x²–y² and 3d_z² orbitals of nickels, respectively. Here μ_z² should be much smaller than μ_x²–y² to ensure that the 3d_z² orbital is near half-filling, while the 3d_x²–y² orbital is close to quarter-filling. In this model, we have ignored the intra-layer (inter-layer) hopping for the 3d_z² (3d_x²–y² ) orbitals and double occupancy of the 3d_z² orbital is not allowed, i.e., the local constraint ${n}_{a,i}^{d}={n}_{a,i,\uparrow }^{d}+{n}_{a,i,\downarrow }^{d}\lt 2$ has been imposed. Actually this effective bi-layer model Hamiltonian is different from the single orbital bi-layer Hubbard model,^[19–27] and the distinct low-energy physics can be expected.

When the intra-layer hybridization between two orbitals t_{x²–y²,z²} is absent, the 3d_z² orbitals are decoupled with the 3d_x²–y² orbitals. Because the intra-layer hopping of the 3d_z² orbitals can be neglected, the ground state of the two-site half-filled 3d_z² orbitals is an isolated singlet. The filling of 3d_x²–y² is close to n^c = 1/2 per lattice site, so the ground state of the interacting 3d_x²–y² electrons behave as a paramagnetic metal. With a finite hybridization t_{x²–y²,z²} , the 3d_z² electrons can hop on the lattice, and the isolated pairs of the 3d_z² electrons can gain coherence and the system may display superconductivity in the large inter-layer coupling limit. This picture shares a similarity to the superconductivity theory of the metallization of σ-bonding band.^[17,18]

Numerical Results of DMRG Study. In order to explore the low energy physics of the effective model Eq. (4), we employ the DMRG method^[28,29] to numerically solve a minimum one-dimensional setup which captures the double-layer structure with the length L = 32 as shown in Fig. 2. According to the DFT results,^[15,16] we choose the hopping parameter t_z² as the energy unit, and t_x²–y² =0.8, t_{x²–y²,z²} = 0.4. For the simplicity, we ignore the Hubbard repulsion for the 3d_x²–y² orbitals because they are far away from half-filling. Here we mainly focus on the large J limit, so the local antiferromagnetic spin coupling is set as J = 0.5. Under the averaged filling n = 3/4, we have a scan of Δ_μ = μ_x²–y² –μ_z² to ensure that the 3d_z² orbitals are 1/16 hole doping away from half-filling and the 3d_x²–y² orbital electrons have 9/16 electron filling. With these chosen parameters, the 3d_z² orbital electrons sit closely to the Mott-insulating limit, while the 3d_x²–y² orbital electrons are in the large doping limit, resembling the orbital selective Mott physics.^[30] In our DMRG calculations, the maximal number of the states is kept up to m = 18000 with a truncation error < 5 × 10⁻⁶. In the following, we will plot the numerical results with large kept states m to indicate the convergence of DMRG calculations.

Zoom In Zoom Out Reset image size

In Fig. 3(a), we give the local charge density distributions in real space of both the 3d_z² and 3d_x²–y² electrons along one row, and the numerical results for other rows have the same values. The obtained results with finite m display nice convergence. A CDW ordering with periodicity of 5 is clearly seen for the 3d_z² orbital electrons, while there seems to be a CDW ordering with wavelength 3 for the 3d_x²–y² electrons, but the modulation is not as regular as that in the 3d_z² orbitals.

Zoom In Zoom Out Reset image size

In the DMRG method, the calculation of the spin-spin correlations is more demanding. Thus we can also apply a pinning magnetic field with the strength h_m = 0.5 at one site in the left edge, which allows us to probe the magnetic structure by calculating the local spin density.^[31] In Figs. 3(b) and 3(c), the local spin densities of both the 3d_x²–y² and 3d_z² orbitals are displayed. Here the numerical results are shown for one row and the absolute values of the different rows are almost the same. We can see that the spin density of both orbital electrons is short-range (disordered), and the system exhibits a paramagnetic behavior.

In order to consider the superconductivity instability, we calculated the equal-time spin singlet superconducting pair–pair correlation function between bond i (formed by site (i, 1) and (i, 2)) and bond j (formed by site (j, 1) and (j, 2)), which is defined as $D(i,j)=\langle {\hat{\varDelta }}_{i}^{\dagger }{\hat{\varDelta }}_{j}\rangle$, where

$$\begin{eqnarray}&&{\hat{\varDelta }}_{i}^{\dagger }=({\hat{c}}_{(i,1),\uparrow }^{\dagger }{\hat{c}}_{(i,2),\downarrow }^{\dagger }-{\hat{c}}_{(i,1),\downarrow }^{\dagger }{\hat{c}}_{(i,2),\uparrow }^{\dagger })/\sqrt{2}.\end{eqnarray} \tag{ 5 }$$

We set the vertical 3d_z² bond as the reference bond in the calculation of the pair–pair correlations. In Fig. 4(a), we show the numerical results for the 3d_z² bonds, and find that the envelope of the pair–pair correlation displays an algebraically decay |D(r)| ∼ r^−K_sc with exponent K_sc = 1.5(1). This pair–pair correlation is always positive and oscillated in real space, consistent with the strong local singlet pairings of the 3d_z² electrons. On the other hand, for the 3d_x²–y² electrons, the envelope of this pair–pair correlation exhibits an algebraic decay with exponent K_sc = 0.74(7) displayed in Fig. 4(b). Moreover, this pair–pair correlation oscillates on the lattice with a special sign structure + – –.

Zoom In Zoom Out Reset image size

Since both the fitted exponents K_sc of the envelopes for 3d_z² and 3d_x²–y² orbital electrons are smaller than 2, the SC correlations are strong enough so that a quasi-long-range order emerges, implying the divergence of the static pair–pair SC susceptibility characterized as

$$\begin{eqnarray}&&{\chi }_{{\rm{sc}}}(T)\sim {T}^{-(2-{K}_{{\rm{sc}}})},\,\,{\rm{when}}\,T\to 0.\end{eqnarray} \tag{ 6 }$$

Moreover, the pair–pair correlations for 3d_x²–y² orbital electrons exhibit oscillation with sign changes, so the emerged superconductivity may be regarded as a pair-density wave ordered state. We note that the intra-layer pairing for the 3d_x²–y² orbitals is stronger than the inter-layer pairing for the 3d_z² orbitals. Though the inter-layer super-exchange favors the formation of singlets for the 3d_z² orbital, these isolated singlets need to hybridize with the itinerant 3d_x²–y² orbital electrons to gain coherence. Therefore, it is reasonable that the intra-layer pairing is stronger than the inter-layer pairing, even though the latter is argued to be the origin of electron pairing in the system.

In summary, we have proposed an effective bi-layer model Hamiltonian to describe the low energy physics of high-T_c superconductivity La₃Ni₂O₇ under high pressure.^[15] In this effective model, we have argued that the σ-bonding band formed from the 3d_z² orbitals via the apical oxygen can be metalized due to the hybridization with the itinerant 3d_x²–y² orbitals, displaying unconventional high-T_c superconductivity. We have also performed DMRG study on a minimum one-dimensional setup with the length L = 32. The DMRG results show instability of CDW-modulated superconductivity. The dominant spin singlet pair–pair correlation is from the 3d_x²–y² orbitals, displaying a pair-density-wave quasi-long-range order. Though this numerical calculation for a minimum setup, we would like to argue that the obtained properties can be used to justify the validity of the proposed effective bi-layer model Hamiltonian for the understanding the high-T_c superconductivity in La₃Ni₂O₇ under high pressure.

Note Added. During the preparation of this work, several theoretical studies^[32–37] appeared on arXiv and the electronic structures and possible pairing instabilities of the high-T_c superconductivity of La₃Ni₂O₇ under high pressure are independently discussed.

Acknowledgments