$\textit{Ab initio}$ derivation of effective Hamiltonian for La$_{2}$CuO$_{4}$/La$_{1.55}$Sr$_{0.45}$CuO$_{4}$ heterostructure

TL;DR

This paper develops a method to derive effective low-energy Hamiltonians for interfaces in strongly correlated electron systems and applies it to cuprate heterostructures, revealing significant parameter changes due to lattice relaxation.

Contribution

The paper extends the multi-scale ab initio scheme to nonperiodic systems and demonstrates its application to cuprate interfaces, highlighting the impact of ionic relaxation on Hamiltonian parameters.

Findings

Parameters differ significantly from bulk La2CuO4, especially Coulomb and level offset.

Lattice relaxation causes CuO6 octahedra distortion and gradual layer dependence of parameters.

Relaxation dramatically alters the level offset and interface occupation numbers.

Abstract

We formulate a method of deriving effective low-energy Hamiltonian for nonperiodic systems such as interfaces for strongly correlated electron systems by extending conventional multi-scale $\textit{ab initio}$ scheme for correlated electrons (MACE). We apply the formalism to copper-oxide high $T_{\rm c}$ superconductors in an example of the interface between overdoped La$_{2-x}$Sr$_{x}$CuO$_{4}$ and Mott insulating La$_{2}$CuO$_{4}$ recently realized experimentally. We show that the parameters of the $E_g$ Hamiltonian derived for the La$_{2}$CuO$_{4}$/La$_{1.55}$Sr$_{0.45}$CuO$_{4}$ superlattice differ considerably from those for the bulk La$_{2}$CuO$_{4}$, particularly significant in the partially-screened Coulomb parameters and the level offset between the $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals, $\Delta E$. In addition, we investigate the effect of the lattice relaxation on the $E_g$ Hamiltonian by carefully comparing the parameters derived before and after the structure optimization. We find that the CuO$_{6}$ octahedra distort after the relaxation as a consequence of the Madelung potential difference between the insulator and metal sides, by which the layer dependence of the hopping and Coulomb parameters becomes more gradual than the unrelaxed case. Furthermore, the structure relaxation dramatically changes the $\Delta E$ value and the occupation number at the interface. This study not only evidences the importance of the ionic relaxation around interfaces but also provides a set of layer-dependent parameters of the $E_g$ Hamiltonian, which is expected to provide further insight into the interfacial superconductivity when solved with low-energy solvers.

Figures

Tables

Table 1. Table 1: Comparison of calculated (calc.) and experimental (expt.) structural parameters of La 2-x Sr x CuO 4 . h O subscript ℎ O h_{\mathrm{O}} is the distance between a Cu atom and its apical oxygen site. The calculations are performed for four different structures. We assign the labels for the undoped cases ( x = 0 𝑥 0 x=0 ) as N1 and N2 and heavily doped cases ( x = 0.45 𝑥 0.45 x=0.45 ) as D1 and D2. In N1 and D1, the lattice constants and internal atomic coordinates are fully relaxed, whereas in N2 and D2 the planer lattice constant is fixed to match that of the substrate LaSrAlO 4 (see the main text for detail). For the tetragonal phase, a 𝑎 a value in the table are 2 ​ a tetra 2 subscript 𝑎 tetra \sqrt{2}a_{\mathrm{tetra}} where a tetra subscript 𝑎 tetra a_{\mathrm{tetra}} is the lattice constant of the tetragonal unit cell, and b 𝑏 b value is left blank.

Label	$x$	$a$ (Å)	$b$ (Å)	$c$ (Å)	$h_O$ (Å)
Calc. (fully relaxed)
N1	0.00	5.398		13.176	2.450
	0.15	5.383		13.247	2.420
	0.30	5.381		13.237	2.390
D1	0.45	5.388		13.186	2.346
Calc. (relax $c$ with $a_{\mathrm{tetra}}=a_{\mathrm{LSAO}}$)
N2	0.00	5.310		13.354	2.481
D2	0.45	5.310		13.377	2.393
Expt. (Ref. Radaelli et al., 1994, 10 K)
	0.00	5.335	5.415	13.117	2.420
	0.15	5.325	5.349	13.197	2.414
	0.30	5.312		13.228	2.390

Table 2. Table 2: Hopping and Coulomb parameters of the bulk E g subscript 𝐸 𝑔 E_{g} Hamiltonians calculated with different doping level x 𝑥 x and structural parameters. The indices m , n 𝑚 𝑛 m,n label the MLWFs; 0 and 1 correspond to d z 2 subscript 𝑑 superscript 𝑧 2 d_{z^{2}} and d x 2 − y 2 subscript 𝑑 superscript 𝑥 2 superscript 𝑦 2 d_{x^{2}-y^{2}} , respectively. The meanings of other parameters are the followings. Δ ​ E Δ 𝐸 \Delta E : The on-site potential difference between the d x 2 − y 2 subscript 𝑑 superscript 𝑥 2 superscript 𝑦 2 d_{x^{2}-y^{2}} and d z 2 subscript 𝑑 superscript 𝑧 2 d_{z^{2}} orbitals; t , t ′ , t ′′ 𝑡 superscript 𝑡 ′ superscript 𝑡 ′′ t,t^{\prime},t^{\prime\prime} : nearest, next-nearest, and second next-nearest hopping parameters on CuO 2 plane, respectively; v , U 𝑣 𝑈 v,U : on-site Coulomb interactions; v n , V n subscript 𝑣 𝑛 subscript 𝑉 𝑛 v_{n},V_{n} : nearest-neighbor Coulomb interactions; J v , J subscript 𝐽 𝑣 𝐽 J_{v},J : exchange interactions; | U / t | 𝑈 𝑡 |U/t| : scaled correlation strength for the d x 2 − y 2 subscript 𝑑 superscript 𝑥 2 superscript 𝑦 2 d_{x^{2}-y^{2}} orbital.

		N1 ($x=0$)	N2 ($x=0$)	D1 ($x=0.45$)	D2 ($x=0.45$)
	$(m,n)$
$\Delta E$		$0.778$	$1.034$	$0.499$	$0.690$
$t$	$(0,0)$	$-0.063$	$-0.059$	$-0.089$	$-0.083$
	$(0,1)$	$0.170$	$0.169$	$0.204$	$0.200$
	$(1,1)$	$-0.463$	$-0.490$	$-0.478$	$-0.500$
$t^{\prime }$	$(0,0)$	$-0.010$	$-0.009$	$-0.015$	$-0.014$
	$(0,1)$	$0.000$	$0.000$	$0.000$	$0.000$
	$(1,1)$	$0.092$	$0.095$	$0.096$	$0.101$
$t^{\prime \prime }$	$(0,0)$	$-0.008$	$-0.007$	$-0.013$	$-0.011$
	$(0,1)$	$0.027$	$0.027$	$0.033$	$0.033$
	$(1,1)$	$-0.074$	$-0.078$	$-0.087$	$-0.091$
Spread of the MLWF (Å2)	$(0,0)$	4.00	3.99	4.72	4.85
	$(1,1)$	3.26	3.19	4.03	4.05
Bare Coulomb parameters
$v$	$(0,0)$	13.85	13.83	11.84	11.66
	$(0,1)$	12.55	12.53	10.74	10.53
	$(1,1)$	14.89	14.90	13.42	13.30
$v_n$	$(0,0)$	3.33	3.37	3.27	3.28
	$(0,1)$	3.63	3.68	3.56	3.59
	$(1,1)$	4.10	4.17	4.10	4.16
$J_v$	$(0,1)$	0.60	0.59	0.51	0.48
Partially screened Coulomb parameters
$U$	$(0,0)$	3.94	3.91	3.02	2.89
	$(0,1)$	2.63	2.60	1.93	1.80
	$(1,1)$	3.92	3.89	3.07	2.92
$V_n$	$(0,0)$	0.63	0.64	0.61	0.59
	$(0,1)$	0.72	0.72	0.69	0.66
	$(1,1)$	0.85	0.86	0.83	0.81
$J$	$(0,1)$	0.50	0.49	0.40	0.38
$\left|U/t\right|$	$(1,1)$	8.47	7.94	6.42	5.84

Table 3. Table S1: Hopping and Coulomb parameters of the bulk E g subscript 𝐸 𝑔 E_{g} Hamiltonians calculated with different doping level x 𝑥 x and structural parameters. Here, the MLWFs are constructed with the wider outer window including 14 valence bands per CuO 2 layer. The indices m , n 𝑚 𝑛 m,n label the MLWFs; 0 and 1 correspond to d z 2 subscript 𝑑 superscript 𝑧 2 d_{z^{2}} and d x 2 − y 2 subscript 𝑑 superscript 𝑥 2 superscript 𝑦 2 d_{x^{2}-y^{2}} , respectively. The meanings of other parameters are the followings. Δ ​ E Δ 𝐸 \Delta E : The on-site potential difference between the d x 2 − y 2 subscript 𝑑 superscript 𝑥 2 superscript 𝑦 2 d_{x^{2}-y^{2}} and d z 2 subscript 𝑑 superscript 𝑧 2 d_{z^{2}} orbitals; t , t ′ , t ′′ 𝑡 superscript 𝑡 ′ superscript 𝑡 ′′ t,t^{\prime},t^{\prime\prime} : nearest, next-nearest, and second next-nearest hopping parameters, respectively; v , U 𝑣 𝑈 v,U : on-site Coulomb interactions; v n , V n subscript 𝑣 𝑛 subscript 𝑉 𝑛 v_{n},V_{n} : nearest-neighbor Coulomb interactions; J v , J subscript 𝐽 𝑣 𝐽 J_{v},J : exchange interactions; | U / t | 𝑈 𝑡 |U/t| : scaled correlation strength for the d x 2 − y 2 subscript 𝑑 superscript 𝑥 2 superscript 𝑦 2 d_{x^{2}-y^{2}} orbital.

		N1 ($x=0$)	N2 ($x=0$)	D1 ($x=0.45$)	D2 ($x=0.45$)
	$(m,n)$
$\Delta E$		$1.257$	$1.645$	$0.950$	$1.302$
$t$	$(0,0)$	$-0.008$	$-0.010$	$-0.010$	$-0.010$
	$(0,1)$	$0.214$	$0.210$	$0.271$	$0.275$
	$(1,1)$	$-0.471$	$-0.496$	$-0.500$	$-0.522$
$t^{\prime }$	$(0,0)$	$-0.039$	$-0.034$	$-0.053$	$-0.051$
	$(0,1)$	$0.000$	$0.000$	$0.000$	$0.000$
	$(1,1)$	$0.088$	$0.089$	$0.093$	$0.094$
$t^{\prime \prime }$	$(0,0)$	$-0.009$	$-0.004$	$-0.020$	$-0.016$
	$(0,1)$	$0.034$	$0.036$	$0.041$	$0.045$
	$(1,1)$	$-0.072$	$-0.077$	$-0.089$	$-0.094$
Spread of the MLWF (Å2)	$(0,0)$	1.37	1.11	2.09	1.69
	$(1,1)$	3.69	3.58	4.04	4.17
Bare Coulomb parameters
$v$	$(0,0)$	21.94	22.51	19.78	20.96
	$(0,1)$	15.77	15.97	14.46	14.72
	$(1,1)$	14.53	14.56	13.56	13.35
$v_n$	$(0,0)$	3.65	3.73	3.60	3.69
	$(0,1)$	3.86	3.93	3.85	3.92
	$(1,1)$	4.04	4.10	4.08	4.12
$J_v$	$(0,1)$	0.79	0.79	0.73	0.73
Partially screened Coulomb parameters
$U$	$(0,0)$	6.67	6.85	5.21	5.42
	$(0,1)$	3.34	3.35	2.56	2.48
	$(1,1)$	3.82	3.79	3.11	2.94
$V_n$	$(0,0)$	0.73	0.75	0.69	0.69
	$(0,1)$	0.78	0.79	0.76	0.75
	$(1,1)$	0.84	0.84	0.83	0.80
$J$	$(0,1)$	0.69	0.68	0.61	0.61
$\left|U/t\right|$	$(1,1)$	8.11	7.64	6.22	5.63

Equations

$$\mathcal{H}$$

$$\displaystyle+\frac{1}{2}\sum_{\sigma\rho}\sum_{\bm{R}\bm{R}^{\prime}}\sum_{mn}\big{\{}U_{m\bm{R}n\bm{R}^{\prime}}a_{m\bm{R}}^{\sigma\dagger}a_{n\bm{R}^{\prime}}^{\rho\dagger}a_{n\bm{R}^{\prime}}^{\rho}a_{m\bm{R}}^{\sigma}$$

$$\displaystyle+J_{m\bm{R}n\bm{R}^{\prime}}(a_{m\bm{R}}^{\sigma\dagger}a_{n\bm{R}^{\prime}}^{\rho\dagger}a_{m\bm{R}}^{\rho}a_{n\bm{R}^{\prime}}^{\sigma}+a_{m\bm{R}}^{\sigma\dagger}a_{m\bm{R}}^{\rho\dagger}a_{n\bm{R}^{\prime}}^{\rho}a_{n\bm{R}^{\prime}}^{\sigma})\big{\}},$$

$$(1-{\delta }_{mn}{\delta }_{\mathbf{R}{\mathbf{R}}^{\prime }})t_{m\mathbf{R}n{\mathbf{R}}^{\prime }}-{\delta }_{mn}{\delta }_{\mathbf{R}{\mathbf{R}}^{\prime }}{\mu }_m$$

$$=⟨{\varphi }_{m\mathbf{R}}⟩{\mathcal{H}}_{\mathrm{KS}}{\varphi }_{n{\mathbf{R}}^{\prime }},$$

$$U_{m\mathbf{R}n{\mathbf{R}}^{\prime }}=⟨{\varphi }_{m\mathbf{R}}{\varphi }_{n{\mathbf{R}}^{\prime }}⟩W_{\mathrm{eff}}{\varphi }_{m\mathbf{R}}{\varphi }_{n{\mathbf{R}}^{\prime }},$$

$$J_{m\mathbf{R}n{\mathbf{R}}^{\prime }}=⟨{\varphi }_{m\mathbf{R}}{\varphi }_{n{\mathbf{R}}^{\prime }}⟩W_{\mathrm{eff}}{\varphi }_{n{\mathbf{R}}^{\prime }}{\varphi }_{m\mathbf{R}},$$

$$W_{\mathrm{eff}}=(1-v{P}_r{)}^{-1}v$$

$$P_{\mathbf{G}{\mathbf{G}}^{\prime }}^d(\mathbf{q},\omega )$$

$$\times \left[\frac{1}{\omega +G_{\alpha \beta }(\mathbf{k},\mathbf{q})}-\frac{1}{\omega -G_{\alpha \beta }(\mathbf{k},\mathbf{q})}\right],$$

$${\rho }_{\alpha \beta \mathbf{k}\mathbf{q}}(\mathbf{G})$$

$$F_{\alpha \beta }^d(\mathbf{k},\mathbf{q})$$

$$\times w_{\alpha \mathbf{k}+\mathbf{q}}{w}_{\beta \mathbf{k}},$$

$$G_{\alpha \beta }(\mathbf{k},\mathbf{q})$$

$$w_{\alpha \mathbf{k}}=\sum _m|U_{\alpha m}^{(\mathbf{k})}{|}^2,$$

$${\varphi }_{m\mathbf{R}}(\mathbf{r})$$

$$=\frac{1}{\sqrt{N}}\sum _{\mathbf{k}}{e}^{-\mathrm{i}\mathbf{k}\cdot \mathbf{R}}\sum _{\alpha }{U}_{\alpha m}^{(\mathbf{k})}{\psi }_{\alpha \mathbf{k}}(\mathbf{r}),$$

$$U_{m0n\mathbf{R}}=\frac{4\pi }{\Omega }\sum _{\mathbf{q}\mathbf{G}{\mathbf{G}}^{\prime }}{e}^{-\mathrm{i}\mathbf{q}\cdot \mathbf{R}}{\rho }_{m\mathbf{q}}(\mathbf{G}){ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-1}(\mathbf{q},\omega =0){\rho }_{n\mathbf{q}}^{\ast }({\mathbf{G}}^{\prime }),$$

$$J_{m0n\mathbf{R}}=\frac{4\pi }{\Omega }\sum _{\mathbf{q}\mathbf{G}{\mathbf{G}}^{\prime }}{\rho }_{mn\mathbf{R}\mathbf{q}}(\mathbf{G}){ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-1}(\mathbf{q},\omega =0){\rho }_{mn\mathbf{R}\mathbf{q}}^{\ast }({\mathbf{G}}^{\prime }),$$

$${\rho }_{mn\mathbf{R}\mathbf{q}}(\mathbf{G})=\frac{1}{N|\mathbf{q}+\mathbf{G}|}\sum _{\mathbf{k}}{e}^{-\mathrm{i}\mathbf{k}\cdot \mathbf{R}}⟨{\psi }_{m\mathbf{k}+\mathbf{q}}^{(w)}⟩e^{\mathrm{i}(\mathbf{q}+\mathbf{G})\cdot \mathbf{r}}{\psi }_{n\mathbf{k}}^{(w)}.$$

$${ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q},\omega )={\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}-[v(\mathbf{q}+\mathbf{G}){]}^{\frac{1}{2}}{P}_{\mathbf{G}{\mathbf{G}}^{\prime }}^r(\mathbf{q},\omega )[v(\mathbf{q}+{\mathbf{G}}^{\prime }){]}^{\frac{1}{2}},$$

$$\begin{gathered} {\tilde{\rho }}_{mn\mathbf{R}\mathbf{q}}(\mathbf{G})\approx \frac{1}{N|\mathbf{q}+\mathbf{G}{|}_{rep.}} \\ \times \sum _{\mathbf{k}}{e}^{-i\mathbf{k}\cdot \mathbf{R}}⟨{\psi }_{m\mathbf{k}+\mathbf{q}}^{(w)}⟩e^{\mathrm{i}(\mathbf{q}+\mathbf{G})\cdot \mathbf{r}}{\psi }_{n\mathbf{k}}^{(w)}, \end{gathered}$$

$$\frac{1}{|\mathbf{q}{|}_{rep.}}={\left[\frac{1}{{\Omega }_{S_{\mathbf{q}}}}{\int }_{S_{\mathbf{q}}}\frac{1}{|\mathbf{k}{|}^2}d\mathbf{k}\right]}^{\frac{1}{2}}.$$

Full text

I Introduction

Interface is at a frontier of condensed matter research and properties and functions not attainable in bulk crystals are the subjects of recent extensive studies. Among all, superconductivity is one of the hottest topics, where interface atomic layers often show properties superior to the bulk in terms of the critical temperature and its stability. Examples are interfaces of copper-oxide superconductors Zhou et al. (2010); Logvenov et al. (2009); Wu et al. (2013); Zheng et al. (2014) and the iron-based superconductors such as FeSe grown on the substrate such as SrTiO3 Ge et al. (2014); He et al. (2013); Qing-Yan et al. (2012).

Among all, recent experimental realization of the pinning of the critical temperature for the interface between overdoped La2-xSrxCuO4 and Mott insulating La2CuO4 Wu et al. (2013) has inspired several theoretical studies Misawa et al. (2016). Experimentally, $T_{\rm c}$ is pinned at 40 K, which is the highest critical temperature of the bulk even when the doping concentration is varied in a wide range of $0.2<x<0.5$ in the overdoped side La2-xSrxCuO4, indicating stable self-optimization of the superconductivity by the interface structure.

However, mechanisms of such fascinating phenomena at interfaces are difficult to identify experimentally in general, because interfaces give only tiny (negligible) contribution to thermodynamic quantities. In addition, the surface sensitive probes such as photoemission and scanning tunnel microscope spectroscopies are not suitable in contrast to surfaces. Therefore, even the lattice constant is hardly determined.

Given this situation, the role of first principles studies, which can predict lattice parameters and atomic positions, becomes more important. Furthermore, in the case of interfaces in strongly correlated electron systems, we need to take into account the effect of electron correlations properly. To study the correlation effect from first principles, derivation of low-energy effective Hamiltonian describing the degrees of freedom near the Fermi level is useful Imada and Miyake (2010); Kotliar et al. (2006). However, the lack of the periodicity makes calculations challenging and so far there exist only few applications to interface systems Hirayama et al. (2012).

To derive effective low-energy Hamiltonians for interfaces in strongly correlated electron systems from first prinpiples, we need to extend the formalism developed for the bulk systems Imada and Miyake (2010). In the formalism for the bulk periodic systems, the low-energy effective Hamiltonians are derived without any adjustable parameters based on the multi-scale ab initio scheme for correlated electrons (MACE) after taking the partial trace summation for the degrees of freedom out of the target low-energy space. The partial trace summation is taken by following either the constrained random phase approximation (cRPA) or the constrained GW (cGW) approximation Hirayama et al. (2013, 2017). All of these formalisms at the moment almost always use the experimental lattice structure and parameters of materials, although the lattice relaxation and optimization can be used without relying on the experimental values if the materials are not available and one wishes to design. However, for the interface, even for the experimentally available systems, the lattice parameters are in many cases not available. Therefore, for the interface calculation, we need first to relax and optimize the lattice structure and predict the precise lattice parameters.

In this paper, we propose a formalism by extending MACE to develop a scheme suitable for nonperiodic systems such as interfaces first by implementing the lattice relaxation near the interface. This procedure is next combined with the conventional MACE treatment. To show its performance, we take an example of the interface between overdoped La2-xSrxCuO4 and Mott insulating La2CuO4. To gain insight beyond the previous work Misawa et al. (2016), we examine the effect of lattice relaxation, which was not studied before.

We derive the two-band effective Hamiltonian consisting mainly of the antibonding band formed from Cu $3d_{x^{2}-y^{2}}$ and O $2p_{\sigma}$ orbitals and the Cu $3d_{3z^{2}-r^{2}}$ ($3d_{z^{2}}$) band for La compounds by extending the bulk studies Hirayama et al. (2018); Jang et al. (2016). The reason why we do not employ the one-band Hamiltonian is that the above two bands are severely hybridized in the case of La2CuO4 Matt et al. (2018).

The structure of the paper is the following: In Sec. II, we present the method. In Secs. III and IV, we show the results of lattice relaxation and derived two-band Hamiltonian parameters for bulk and interface systems, respectively. Secs. V and VI are devoted to discussion and conclusion of the paper.

II Method

We derive ab initio low-energy effective Hamiltonians for La2CuO4 and its heterostructures by employing the cRPA scheme Aryasetiawan et al. (2004) and the maximally localized Wannier function (MLWF) method Marzari et al. (2012) on top of density functional theory (DFT) calculations.

We first choose the target low-energy subspace ($d$ subspace) to construct the effective Hamiltonian for $d$-subspace electrons. In the case of La2CuO4, we adopt the two-band Hamiltonian comprising the Cu $3d_{x^{2}-y^{2}}$-like and $3d_{z^{2}}$-like orbitals near the Fermi level. Unlike the case of HgBa2CuO4 where the one-band Hamiltonian should be a good minimum model, the two $3d$ orbitals are strongly entangled in La2CuO4. Therefore, the minimum effective Hamiltonian of La2CuO4 should include at least two orbitals. Note that, in what follows, the “Cu $3d_{x^{2}-y^{2}}$-like orbitals” are actually the antibonding orbital of the strongly hybridizing copper $3d_{x^{2}-y^{2}}$ and oxygen $2p_{\sigma}$ orbital.

The form of the effective two-orbital Hamiltonian is

[TABLE]

where $a_{n\bm{R}}^{\sigma\dagger}$ ($a_{n\bm{R}}^{\sigma}$) is a creation (annihilation) operator of an electron with spin $\sigma$ in the $n$th Wannier orbital located in the unit cell at position $\bm{R}$, $t_{m\bm{R}n\bm{R}^{\prime}}$ is the hopping parameter, $\mu_{m}$ is the orbital-dependent on-site potential, $n_{m\bm{R}}^{\sigma}=a_{m\bm{R}}^{\sigma\dagger}a_{m\bm{R}}^{\sigma}$, and $U_{m\bm{R}n\bm{R}^{\prime}}$ and $J_{m\bm{R}n\bm{R}^{\prime}}$ are effective Coulomb and exchange interactions, respectively. Let $\mathcal{H}_{\mathrm{KS}}$ denote the Kohn-Sham (KS) Hamiltonian of DFT calculation and $W_{\rm eff}$ be the partially-screened Coulomb interaction, the hopping and Coulomb parameters are expressed as follows:

[TABLE]

with $\ket{\phi_{m\bm{R}}}=a_{m\bm{R}}^{\dagger}\ket{0}$ and $\delta_{ij}$ being the Kronecker delta.

In the calculation of the effective interaction parameters in Eqs. (3) and (4), we exclude the screening contribution associated with the $d$-$d$ polarization processes within the $E_{g}$ manifold of the $3d$ electron subspace (namely, $d$ subspace) and use the partially screened interactions $W_{\rm eff}$: The $d$-subspace screening contribution is considered when we analyze the low-energy Hamiltonians, therefore, we need to exclude it in deriving $U_{m\bm{R}n\bm{R}^{\prime}}$ and $J_{m\bm{R}n\bm{R}^{\prime}}$ to avoid the double counting of the screening Aryasetiawan et al. (2004). Then $W_{\rm eff}$ is given by

[TABLE]

where $P_{r}$ is given by $P_{r}=P-P_{d}$ with the full polarization $P$ and $d$-subspace polarization $P_{d}$, and $v$ is the bare Coulomb interaction.

When the $d$ subspace is not isolated from the high-energy subspace ($r$ subspace), which is usually the case in cuprates, it is necessary to handle the entanglement to construct the $d$-subspace polarization $P_{d}$ Miyake et al. (2009); Şaşıoğlu et al. (2011). In this study, we employ the simple approach of Ref. Şaşıoğlu et al., 2011, where the matrix element of $P_{d}$ in the plane-wave basis is given as

[TABLE]

Here, $\epsilon_{\alpha\bm{k}}$ and $\ket{\psi_{\alpha\bm{k}}}$ are the KS eigenvalue and wavefunction of the $\alpha$th band at the momentum $\bm{k}$, respectively, $\bm{G}$ is the reciprocal lattice vector, $\epsilon_{\mathrm{F}}$ is the Fermi level, $\theta(x)$ is the Heaviside step function, and $\delta$ is a small negative value. The term $w_{\alpha\bm{k}}$ in Eq. (8) is the weight of the target Wannier orbital defined as

[TABLE]

with $U_{\alpha m}^{(\bm{k})}$ being a unitary matrix that transforms $\psi_{\alpha\bm{k}}(\bm{r})$ into $\phi_{m\bm{R}}(\bm{r})$ as

[TABLE]

where $\psi_{m\bm{k}}^{(w)}$ is the Wannier-gauge Bloch wavefunction. Hence, $w_{\alpha\bm{k}}$ can be obtained straightforwardly from the results of the MLWF calculation. If the KS wavefunction $\ket{\psi_{\alpha\bm{k}}}$ is strictly represented only within the subspace of the $d$ electrons, we obtain $w_{\alpha\bm{k}}=1$. Therefore, the term $w_{\alpha\bm{k}+\bm{q}}w_{\beta\bm{k}}$ in $F^{d}_{\alpha\beta}(\bm{k},\bm{q})$ becomes exactly one when both the virtual state $\ket{\psi_{\alpha\bm{k}+\bm{q}}}$ and the occupied state $\ket{\psi_{\beta\bm{k}}}$ belong to the $d$ subspace, and the $d$-$d$ transition in $P^{r}_{\bm{G}\bm{G}^{\prime}}(\bm{q},\omega)=P_{\bm{G}\bm{G}^{\prime}}(\bm{q},\omega)-P^{d}_{\bm{G}\bm{G}^{\prime}}(\bm{q},\omega)$ is duly excluded.

Once we obtain $P_{d}$ by Eq. (II), we can calculate the partially screened Coulomb parameters in Eqs. (3) and (4) as

[TABLE]

where $\Omega$ is the crystal cell volume, $\epsilon^{-1}(\bm{q},\omega)$ is the inverse of the symmetric dielectric matrix, and $\rho_{m\bm{q}}(\bm{G})=\rho_{mm\bm{0}\bm{q}}(\bm{G})$ with $\rho_{mn\bm{R}\bm{q}}(\bm{G})$ being defined as follows:

[TABLE]

The symmetric dielectric matrix is defined as

[TABLE]

where the bare Coulomb interaction $v(\bm{q})$ in reciprocal space is given by $v(\bm{q})=4\pi/\Omega|\bm{q}|^{2}$.

III $E_{g}$ Hamiltonian of bulk systems

To clarify specific properties at the LCO/LSCO interface, it is essential to first understand properties of the bulk systems including the doping-level dependence of the low-energy Hamiltonians. In this section, we carefully compare the $E_{g}$ Hamiltonians of the non-doped La2CuO4 and the overdoped La1.55Sr0.45CuO4.

All of the DFT calculations in this work were performed by using Quantum ESPRESSO Giannozzi et al. (2009), which implements the plane-wave pseudopotential method. We employed the Pewdew–Burke–Ernzerhof (PBE) exchange-correlation potential Perdew et al. (1996) and the optimized norm-conserving Vanderbilt (ONCV) pseudopotentials Hamann (2013) from the SG15 table Schlipf and Gygi (2015). The kinetic-energy cutoff was set to 100 Ry. We employed the tetragonal structure of La2-xSrxCuO4 (space group: $I4/mmm$), where the Sr doping was modeled by the virtual crystal approximation (VCA).

When constructing the low-energy Hamiltonians of the bulk systems, we employed the conventional unit cell (see Fig. 1) because it is more convenient than the primitive cell to apply consistent Wannierization parameters, particularly the frozen window, between the bulk and heterostructure systems. In the conventional unit cell calculations, we employed the 8$\times$8$\times$4 $\bm{k}$ points for the Brillouin zone (BZ) integration with the smearing width of 0.02 Ry.

III.1 Structural properties

First, we fully relaxed the tetragonal structure of La2-xSrxCuO4 at different doping levels $x$ until both the force convergence criteria $|F|<10^{-6}$ Ry$\cdot$Bohr*-1* and the stress convergence criteria $|\sigma|<0.5$ kbar are satisfied. The results are compared with the available experimental data in Table 1. According to the paper of Radaelli et al. Radaelli et al. (1994), the tetragonal phase of La2-xSrxCuO4 is stable at 0 K only when $x>0.21$, below which the orthorhombic phase becomes the most stable at low temperatures. Nonetheless, it should be reasonable to use the tetragonal structure for the derivation of the effective Hamiltonian because no discontinuity in the $T_{\mathrm{c}}$ has been observed at the tetragonal-to-orthorhombic phase transition Radaelli et al. (1994). The in-plane lattice constants of the orthorhombic phase are approximately equal to $\sqrt{2}a_{\mathrm{tetra}}$. Therefore, we show $\sqrt{2}a_{\mathrm{tetra}}$ instead of $a_{\mathrm{tetra}}$ in Table 1 for the purpose of comparison.

It is observed in the table that our DFT results based on the VCA reasonably well reproduce the experimental cell parameters within 1% error as well as the trend of the doping-level dependence, thus validating the reliability of the VCA. In Table 1, we also show the optimized values of the out-of-plane lattice constant $c$ with $a$ and $b$ being fixed to those of LaSrAlO4 (LSAO), which was used as the substrate for growing the LCO/LSCO bilayer thin film Zhou et al. (2010). Since the in-plane lattice constants of LSAO is slightly smaller than those of LCO and LSCO, the LSAO substrate induces the compressive strain along the $a$ and $b$ axes, leading to a slightly larger $c$ value due to positive Poisson ratio.

III.2 Construction of MLWFs

Second, we construct the MLWFs of La2CuO4 and La1.55Sr0.45CuO4 to see the doping-level and strain dependence of the hopping parameters and the Coulomb interaction. Since the target $E_{g}$ orbitals are not isolated, the resulting MLWFs are rather sensitive to the chosen range of the energy window. In this study, we set the outer window by band index. Here, the outer window specifies the Hilbert subspace, within which the MLWFs are constructed. We included 8 valence bands per CuO2 layer in the outer window (see Fig. 2), which was the narrowest window setting to match the MLWF band structures with the KS ones around the Fermi level both for the non-doped and doped systems. Also, the frozen window Marzari et al. (2012) was used to perfectly reproduce the original KS band of the $d_{x^{2}-y^{2}}$-like orbitals at the Fermi energy. The resulting two MLWFs are rather extended as shown in Fig. 3.

It is possible to derive more localized Wannier orbitals by including more valence bands in the outer window. To see the influence of the outer window range on the effective Hamiltonian, we also show the calculated parameters when we included 14 valence bands per CuO2 layer as in Ref. Hirayama et al. (2018) in the supplemental material (SM). This condition still excludes the bonding and non-bonding states formed by the Cu $d_{x^{2}-y^{2}}$ and in-plane O $p_{\sigma}$ orbitals from the energy window but newly includes the bonding state formed by the Cu $d_{z^{2}}$ and apical-oxygen $p_{z}$ orbitals. Therefore, the resulting $d_{z^{2}}$-like MLWF becomes more localized than that of Fig. 3(a), while the changes of the $d_{x^{2}-y^{2}}$ orbital shape and parameters are small (see SM). As the two-band ($E_{g}$) description of the effective Hamiltonian, we believe that the choice of 8-band outer window is more appropriate, because the $d_{z^{2}}$ orbital is strongly hybridized with apical oxygen $p_{z}$ orbital. This strong hybridization is ignored if we employ the 14-band outer window.

III.3 cRPA calculation

In the present cRPA calculations for bulk, we considered the particle-hole excitations within 150 bands (75 bands/f.u.) in calculating the polarization, which corresponds to include 84 unoccupied bands (42 bands/f.u.) up to $\sim$ 21 eV above the Fermi level. The kinetic-energy cutoff for the polarization was set to 20 Ry, which was sufficiently large for the symmetric dielectric matrix to reach the large $|\bm{G}|$ limit of $\epsilon_{\bm{G}\bm{G}^{\prime}}\approx\delta_{\bm{G}\bm{G}^{\prime}}$. The same computational conditions were employed for the cRPA calculation of interfaces. The number of unoccupied bands was selected to make the cRPA calculation of the complex LCO/LSCO heterostructure feasible, but the screened Coulomb parameters, particularly the on-site Coulomb interaction, will be reduced further with increasing the number of bands. Fortunately, our setting still gives reasonably converged values of $U$. For example, in the N1 case, we obtained 3.92 eV for the on-site Coulomb interaction of the $d_{x^{2}-y^{2}}$ orbital. This value reduced to the almost converged result of 3.67 eV which were obtained with 500 bands. Therefore, the presented results for $U$ are within $\sim$ 7% error (overestimate) from the converged values. Our $U$ value is in reasonable agreement with those of the previous cRPA Werner et al. (2015); Jang et al. (2016) and cGW Hirayama et al. (2018) studies as described in Appendix A.

The values of effective Coulomb interactions computed by cRPA are determined by the shape (spatial spread) of MLWFs and the strength of screening. The doping and strain affect both of them. To see only the effect of the change in the shape of WLWFs, we also compute the Wannier matrix elements of the bare Coulomb interaction by replacing $W_{\rm eff}$ with $v$ in Eqs. (3) and (4).

III.4 Doping and strain dependence

Figure 2 compares the calculated KS and MLWF band structures of LCO and LSCO along the high-symmetry lines of the BZ. Here, the fully optimized structures (N1 and D1 in Table 1) are used. While the valence-band structures of the two systems are similar, a simple rigid-band picture seems insufficient for explaining the difference. To see the difference more clearly, we constructed atomic-like 17 MLWFs of LCO and LSCO from the isolated 17 valence bands consisting of the whole Cu $3d$ and O $2p$ manifold, and compared their onsite energy levels relative to the Fermi energy. We then observed that upon hole doping the energy levels of O $2p$ orbitals increased significantly, particularly for the apical oxygen, whereas those of Cu $3d$ orbitals changed only slightly. Hence, the hole doping by the chemical substitution reduces the energy-level difference of O $2p$ and Cu $3d$ orbitals. The observed shift of the O $2p$ energy levels can be attributed to the change of the electrostatic potential induced by the substitution, which mainly affects the O $2p$ orbitals of oxygen located near the La sites (see Fig. 1).

The hole doping also affects the hopping and Coulomb parameters as can be seen in Table S1. The significant change occurs for the partially screened on-site Coulomb parameters, whose reduction amounts to $\sim$ 21–31%. This reduction can be attributed to two different factors. One is the change of the spread of the MLWFs. Since the hole doping reduces the energy-level difference between the O $2p$ and Cu $3d$ orbitals as mentioned above, the $2p$-$3d$ hybridization becomes stronger, leading to more extended $E_{g}$ Wannier functions. The change of the MLWFs resulted in the 10–16% reduction of the bare on-site Coulomb interactions. The other factor is the change of the screening strength. Upon hole doping, the energy levels of the O $2p$ orbitals become closer to the Fermi energy. This enhances the screening channel in $P_{r}$ and explains further reduction in $U$. As for the intraorbital hopping parameters, we see a slight increase upon doping. As a result of the decrease of the Coulomb interaction and the slight increase of the hopping parameters, the strength of correlation, which can be measured by the ratio $|U/t|$ for $d_{x^{2}-y^{2}}$ orbital, changes significantly. $|U/t|$ changes from $\sim$ 8.5 (7.9) to $\sim$ 6.4 (5.8) for the fully-relaxed structure (relaxed structure with the constraint of $a=a_{\mathrm{LSAO}}$).

The level offset between the $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals $\Delta E=\mu_{x^{2}-y^{2}}-\mu_{z^{2}}$, which has been pointed out to be a key parameter to explain the material dependence of $T_{c}$ Sakakibara et al. (2010, 2012a, 2012b), is also affected by hole doping. It decreases by 0.279 eV and 0.344 eV for the fully-relaxed structure and the relaxed structure with the constraint of $a=a_{\mathrm{LSAO}}$, respectively.

Next, we discuss the strain dependence. To see the effect of the compressive stress along the $ab$-plane induced by the LSAO substrate, we also calculated the hopping and Coulomb parameters with the N2 and D2 structures of Table 1. As shown in Table S1, the compressive stress changes the Coulomb parameters only slightly but significantly increases the level offset $\Delta E$ by $\sim$ 0.19–0.26 eV. This tendency agrees with the previous theoretical result Sakakibara et al. (2012b).

We see that the doping and strain affect the parameters in the effective Hamiltonians. However, the changes of parameters upon doping and/or strain are usually neglected in previous studies. It is intriguing to study the superconducting amplitude and its competition with the charge inhomogeneity such as stripes by solving the present Hamiltonian using highly accurate low-energy solvers.

IV Eg Hamiltonian of heterostructures

To derive the $E_{g}$ Hamiltonian of the LCO/LSCO interface, we employ the superlattice (SL) structures schematically shown in Fig. 4. The heterostructure was constructed by stacking tetragonal unit cells of LCO and LSCO along $c$ axis. For the LCO and LSCO structure units, we employed the N2 and D2 conventional cells in Table S1, respectively. Since each unit cell of LCO contains 7 atoms, the structural model of a $(M,N)$ superlattice contains 7$\times(M+N)$ atoms, where $M$ and $N$ are the numbers of unit cells in the LCO and LSCO regions, respectively. Second, we performed DFT calculations of the $(M,N)$ superlattice and optimized the lattice constant along the $c$ axis as well as the internal coordinates, while the $a$ value was fixed to $a_{\mathrm{LSAO}}$. Third, we constructed MLWFs of all Cu $E_{g}$ orbitals, corresponding to $2\times(M+N)$ total orbitals, and calculated the hopping, bare Coulomb, and partially-screened Coulomb parameters. To see the convergence of the parameters with respect to the number of layers, we calculated hopping parameters for $(M,N)=$ (4,4), (6,6), and (8,8). The Coulomb parameters were calculated only for $(M,N)=$ (4,4) and (6,6) due to the computational limitations. In the structural optimization, we employed the 8$\times$8$\times$1 $\bm{k}$ points. In the subsequent cRPA calculations, the 8$\times$8$\times$2 $\bm{k}$-point mesh was employed in order to use the tetrahedron method for an accurate treatment of the summation over $\bm{k}$ in Eq. (II). In the calculations of SLs, the Brillouin zone becomes highly anisotropic. In this case, we found that the original definition of Eq. (14) needs to be modified to perform the stable cRPA calculation as is discussed in detail in Appendix B. Our modified $\tilde{\rho}_{mn\bm{R}\bm{q}}(\bm{G})$ gives physically correct $R$ dependence of the Coulomb parameters as shown in Fig. 11.

In this study we do not study the effect of interlayer atomic diffusion, which makes the structure of the interface slightly obscured Logvenov et al. (2009).

IV.1 Effect of structural optimization

The formation of an LCO/LSCO interface introduces an abrupt change of the electrostatic potential near the interface, which is energetically unfavorable. Therefore, the internal coordinates can deviate from the bulk values, which is likely to influence the DFT band structures and the associated $E_{g}$ Hamiltonian.

Figure 5 shows the layer dependence of the apical oxygen height obtained after performing the structural optimization. Since the mirror plane symmetry of the CuO2 planes is lost due to the interface, the $h_{\mathrm{O}}$ values of the apical oxygen atoms above and below a CuO2 plane are different from each other. To distinguish these two, we use $h_{\mathrm{O}}^{\pm}$ as defined in Fig. 4. After the relaxation, the $h_{\mathrm{O}}^{+}$ value, which is the distance to the apical oxygen on the LCO side, tends to increase from the bulk value, whereas $h_{\mathrm{O}}^{-}$ shows the opposite tendency. The difference $|h_{\mathrm{O}}^{+}-h_{\mathrm{O}}^{-}|$ becomes significant at the interface (layers $\pm 1$) and sharply decreases as going away from the interface. These structural changes were commonly observed in the studied SLs with different sizes.

All of these behaviors can be understood qualitatively from the Mardelung potential. The ionic charge of La and Sr near the interface generates electric field from the LCO side to the LSCO side at the interface, because the CuO2 planes closest to the interface, one in the LCO side and the other in the LSCO side, are separated by a LaO plane in the LCO side and a (La,Sr)O plane in the LSCO side. Namely, because of the charge imbalance between the LaO and (La,Sr)O planes, it generates an electric field in the direction from the LSCO to the LCO sides, which causes the upward shift of the oxygen and downward shift of the copper atoms in the configuration illustrated in Fig. 6. This electric field also induces the electron transfer from the LCO side to the LSCO side and the electron itinerancy (interlayer electron hopping) makes its transfer range vertical to the interface wider. It also makes a wider range of atomic-position shift in a self-consistent fashion. Since the electric field is strongest at the interface, the distortion is of course largest at the interface and becomes smaller at points far from the interface.

The inversion symmetry breaking due to the structure relaxation induces a stretch (contraction) of the octahedron along $c$ axis in the LSCO (LCO) side (Fig. 6), whose effect is discussed in detail in Sec. V. In addition, in plane Cu and O are not aligned in plane any more because of anti-phase distortion between Cu and O ions. The deviation of the O-Cu-O angle from 180∘ amounts to $\sim$4*∘* at the interface.

The structural relaxation considerably affects the electronic structure of the heterostructure as evidenced in Fig. 7. For example, the orbital energies of the $d_{x^{2}-y^{2}}$-like orbital along the line M–X change very sharply at the interface before the structure optimization, in accord with the previous DFT calculation Misawa et al. (2016). After the optimization, the orbital-energy shift occurs and the energy change becomes more gradual. The shift of the energy-level is noticeable in the layers $\pm 1$, particularly around the point X, but it is far smaller in the other layers ($\pm 2$, $\pm 3$). This behavior is consistent with the rather strong deformation of the O-Cu-O angle in the layers $\pm 1$ and its rapid recovery in the layers $\pm 2$ observed in Fig. 5. Therefore, these results indicate that the structure optimization influences the hopping parameters mainly at the layers $\pm 1$, which is investigated in the subsequent section.

IV.2 Layer-dependent hopping and Coulomb parameters

Figure 8 shows the layer dependence of the level offset $\Delta E$, the occupation number $n_{m}=\braket{a_{m\bm{0}}^{\dagger}a_{m\bm{0}}}$, the dominant part of the hopping parameters and the spread of the MLWFs calculated for the $(4,4)$, $(6,6)$ and $(8,8)$ SLs. To see the effect of the structure optimization, we compare the results before (open symbols) and after (filled symbols) the optimization in the figure.

We observe that the structural optimization affects all of the parameters especially at the layers $\pm 1$, leading to a more gradual layer dependence. The change of the slope is particularly significant in the layer dependence of the level offset $\Delta E$ [Fig. 8(a)], hole concentration given by $3-n_{z^{2}}-n_{x^{2}-y^{2}}$ [Figs. 8 (b) and (c)], and the interband hopping [Fig. 8 (g)]. For example, the hole concentration at the layer $+1$ changes from 0.01 to 0.13 by the structure optimization. The latter is close to an optimal doping level $x=0.15$ of the bulk La2-xSrxCuO4 at which a maximum $T_{\mathrm{c}}$ has been observed.

The layer dependence of the onsite Coulomb parameters calculated for the $(4,4)$ and $(6,6)$ SLs are shown in Fig. 9. As in the case of bulk systems discussed in Sec. III, the effective Coulomb interactions tend to be larger in the LCO side, which can be attributed to the smaller Wannier spread and the weaker screening. The structure relaxation makes the layer dependence of interaction parameters gradual. As a result, at the layer $+1$, the $U$ value of the $d_{x^{2}-y^{2}}$ orbital is 3.24 eV, which is about 15% smaller than the bulk LCO value.

We also calculated the layer dependence of the off-site Coulomb $V_{n}$ and exchange $J_{x^{2}-y^{2}\bm{0};z^{2}\bm{0}}$ interactions (see Fig. S1 of the SM). We observe that the layer dependence of the off-site Coulomb parameter $V_{n}$ is weak; $V_{x^{2}-y^{2}}$ and $V_{z^{2}}$ values change in the range of 0.84–0.89 eV and 0.65–0.67 eV, respectively. By contrast, the layer dependence of the exchange parameter $J_{x^{2}-y^{2}\bm{0};z^{2}\bm{0}}$ is rather significant, which changes from 0.47 eV at the layer $-3$ to 0.39 eV at the layer $+3$. These tendencies are consistent with the doping-level dependence of $V_{n}$ and $J$ observed in the bulk systems (Table S1).

Figure 10 shows the layer dependence of the scaled correlation strength $|U/t|$ calculated for the $d_{x^{2}-y^{2}}$ orbital. Near the LCO/LSCO interface, the $|U/t|$ value decreases from the bulk LCO value of 7.94 to 6.68 at the layer $+1$, which amounts to a 15% reduction. Since the relative stability of the superconducting phase over other competing phases changes rather sensitively with the $|U/t|$ value, considering the 15% reduction of $|U/t|$ would be necessary to explain the unique properties of the superconductivity at interfaces quantitatively.

Finally, we discuss the convergence of parameters to bulk values. We see that the calculated Coulomb parameters and the Wannier spread of the SLs did not reach the bulk LCO and LSCO values even at the layers farthest from the interface, indicating rather strong sensitivity of these parameters to the Madelung potential difference induced by the interface. This issue is expected to be resolved by using a much larger SL, which was not pursued in this study owing to the computational limitations. Notwithstanding, since the two different sizes of the SL calculation shows more or less the same behaviors, the Coulomb parameters near the interface are likely to be already converged and therefore reliable enough for studying the superconductivity at the LCO/LSCO interface by low-energy solvers.

V Discussion

Near the interface, $\Delta E$, transfers, hole concentrations and the interactions all show substantially more gradual change with moving from the LSCO side to the LCO side than the unrelaxed case. We here discuss that all the above characteristic features are explained by a basic principle “the nature relaxes to avoid discontinuous changes”. This principle is manifested concretely in the real material in the following way.

Before the structural relaxation, the electrostatic potential of the LaO layer changes abruptly, and it generates rather strong electric field at the interface. However, the shift of the atomic position and the electronic charge redistribution described in Sec. IVA leads to more gradual transition of the electron concentration between the LCO and LSCO sides. Since the negatively-charged apex oxygen approaches (moves away from) the CuO2 layer on the LCO (LSCO) side on average as is illustrated in Fig. 6, the electronic level at the CuO2 layer is raised (lowered) in the LCO (LSCO) side in comparison to the unrelaxed lattice. It enhances the electron transfer from the LCO side to the LSCO side originally caused by the electric field at the interfacial charge imbalance (see Sec. IVA) as is seen in Fig. 8(b), (c). Since the onsite level of $d_{x^{2}-y^{2}}$-like orbital is higher than that of $d_{z^{2}}$-like orbital, the electron density decreases mainly from the $d_{x^{2}-y^{2}}$-like orbital in the LCO side and increases mainly in the the $d_{z^{2}}$-like orbital in the LSCO side as shown in Figs. 8(b) and (c).

The level offset between the $d_{x^{2}-y^{2}}$-like and $d_{z^{2}}$-like orbitals, $\Delta E$, is determined mainly by the ligand field of the six oxygen ions surrounding Cu. The contraction (stretch) of the octahedron for the LCO (LSCO) side leads to the decrease (increase) of the average distance to the apical O from Cu, $(h_{\mathrm{O}}^{+}+h_{\mathrm{O}}^{-})/2$, in the LCO (LSCO) side (see Fig. 5). Because the onsite level of $d_{z^{2}}$ orbital is more sensitive to $h_{\mathrm{O}}$ than $d_{x^{2}-y^{2}}$ orbital, a larger $h_{\mathrm{O}}$ value leads to a larger $\Delta E$. Therefore, $\Delta E$ decreases (increases) on the LCO (LSCO) side after the structure relaxation as shown in Fig. 8(a).

The nearest neighbor transfer $t_{x^{2}-y^{2}}$ relatively decreases from the unrelaxed lattice value [Fig. 8(d)] despite a slight increase in the Wannier spread [Fig. 8(i)]. This is presumably because the inplane Cu and O do not align in the flat plane any more but form a zigzag alignment after the structural optimization.

The screened interaction of the antibonding $d_{x^{2}-y^{2}}$-like and $d_{z^{2}}$-like electrons decreases (increases) on the LCO (LSCO) side in comparison to the unrelaxed case. The structure relaxation increases (decreases) the hole concentration on the LCO (LSCO) side. As we see in Sec. III.4, larger hole concentration enhances the screening. Therefore increase (decrease) of the hole concentration may explain the reduction (enhancement) of $U$ values on the LCO (LSCO) side.

In total, the octahedron distortion leads to more gradual layer dependence of all the quantities than the unrelaxed case following the above principle.Such weakened and more gradual layer dependence than that before lattice relaxation seems to follow a general principle: When the spatial gradient becomes strong, the system reacts to weaken it by screening such strong spatial dependence in analogy to Le Chatelier’s priciple in the time dependence. Such screening effects may also work to weaken the effect of impurities, randomness and interface roughness in real interface.

More gradual layer dependence of the onsite level than that before the lattice relaxation may play a role to stabilize the superconductivity. In bulk compounds, there are strong tendency towards the in-plane charge inhomogeneity which suppresses superconductivity Zheng et al. (2017); Ido et al. (2018); Darmawan et al. (2018). However, the interface systems might be able to avoid the in-plane inhomogeneity by making out-of-plane inhomogeneity (interlayer phase separation). In Ref. Misawa et al. (2016), it has been argued that the gradual layer dependence of the onsite level is favorable to realize interlayer phase separation and stabilize superconductivity. It is indeed interesting to analyze the derived interface Hamiltonian to study the stability of the superconductivity.

Finally, the sensitivity of onsite-level, hopping and Coulomb interaction parameters to the lattice distortion suggests nonnegligible coupling between electron and lattice degrees of freedom. Therefore, it is also interesting to investigate the role of electron-phonon couplings in the interface, which might play a role to enhance superconductivity as suggested in other interface systems such as FeSe on SrTiO3 substrate Huang and Hoffman (2017).

VI Conclusion

We have derived the ab initio $E_{g}$ Hamiltonian of La2CuO4/La1.55Sr0.45CuO4 superlattices by performing large-scale cRPA calculations based on DFT. We have shown that the level offset between the $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals, $\Delta E$, and the Coulomb parameters near the interface become smaller than those of the LCO bulk values. This occurs even without performing structure relaxation and can be attributed to the change of the electrostatic potential induced by substituting La with Sr, which increase the energy level of the oxygen $2p$ orbitals more significantly than that of Cu $3d$ orbitals. After the structural relaxation, the layer dependence of the hopping and Coulomb parameters becomes more gradual than the unrelaxed case, which results from the contraction and stretch of CuO6 octahedron in the LCO and LSCO side, respectively. The effect of the structural relaxation is particularly noticeable in $\Delta E$ and the occupation number. Since these parameters as well as the scaled correlation strength $|U/t|$ influence the stability of the superconducting state, the modulation of these parameters at the interface reported in this work should be considered for developing robust understandings of the unique superconducting properties observed in interfaces.

Acknowledgements.

The authors thank Kazuma Nakamura, Motoaki Hirayama, Takahiro Ohgoe, Takahiro Misawa, Youhei Yamaji, Kota Ido, and Kosuke Miyatani for fruitful discussion. This work was partially supported by the MEXT HPCI Strategic Programs, and the Creation of New Functional Devices and High-Performance Materials to Support Next Generation Industries (CDMSI) and a Grant-in-Aid for Scientific Research (No. 16H06345) from Ministry of Education, Culture, Sports, Science and Technology, Japan. The computation in this work has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, The University of Tokyo. Y. N. was supported by Grant-in-Aids for Scientific Research (JSPS KAKENHI) (No. 17K14336 and 18H01158).

Appendix A Comparison of the present cRPA calculation for bulk La2CuO4 with previous works

In the original disentanglement method of Şaşıoğlu et al. Şaşıoğlu et al. (2011), the $d$-$d$ polarization process is excluded via the weight of the $E_{g}$ subspace $w_{\alpha\bm{k}}$. Therefore, the screening process involving the $d_{x^{2}-y^{2}}$-like KS orbitals crossing the Fermi energy can be avoided perfectly if the $w_{\alpha\bm{k}}$ values become exactly one for these KS orbitals. For La2-xSrxCuO4, however, we observed that the weight slightly deviates from one even near the Fermi energy when the KS energy $\epsilon_{\alpha\bm{k}}$ is outside the frozen window. For example, the weight of the $d_{x^{2}-y^{2}}$-like KS orbitals at $\bm{k}=(\frac{\pi}{4},0,0)$ was around 0.97, and the remaining weight of 0.03 origites from the orbitals outside the $E_{g}$ subspace, which contribute to $P_{r}$. Since we want to avoid such a small screening process from the $d_{x^{2}-y^{2}}$-like KS orbitals completely as the disentanglement method of Miyake et al. Miyake et al. (2009) does, we updated the $w_{\alpha\bm{k}}$ values in such a way that those close to one (zero) becomes exactly one (zero) while keeping the total weight $\sum_{\alpha}w_{\alpha\bm{k}}$ unaltered. Such a treatment did not change the results significantly but increased the partially-screened Coulomb parameters by $\sim$5% compared with the original treatment, leading to better agreement with the previous result of $U_{x^{2}-y^{2}}=4.2$ eV Miyatani (2014) which was obtained based on the full potential linearized muffin tin orbital (FP-LMTO) method and the disentanglement method of Miyake and coauthors. The present result of $U_{x^{2}-y^{2}}\sim 3.9$ eV (3.67 eV with 500 bands) is still smaller than 4.2 eV, which can most likely be attributed to the difference in the Wannierization parameters and/or the pseudopotentail adopted in this study. Indeed, the bare Coulomb parameters in this study is also smaller than the FP-LMTO based study by $\sim$15%.

In the previous cRPA studies of LCO, the $U_{x^{2}-y^{2}}$ values of 3.65 eV Werner et al. (2015) and 3.15 eV Jang et al. (2016) have been reported, which agree reasonably well with our result especially given that the $U$ value is rather sensitive to the detail of the Wannierization procedure and the resulting spread of the Wannier orbital. The result of Ref. Werner et al. (2015) was obtained for the one-band Hamiltonian, where the Wannier function was constructed without the frozen window near the fermi energy. Therefore, the shape (spread) of their $d_{x^{2}-y^{2}}$-like Wannier orbital is similar to ours (Fig. 3(b)), resulting in the similar $U_{x^{2}-y^{2}}$ values. If the frozen window is used when constructing the one-band Hamiltonian, the resulting Wannier orbital should be more extended and the $U_{x^{2}-y^{2}}$ value should become smaller because of the strong hybridization between the $E_{g}$ orbitals in LCO. Compared to the results of Ref. Werner et al. (2015) and ours, the $U_{x^{2}-y^{2}}$ value of Ref. Jang et al. (2016) seems somewhat smaller, whose origin is unclear due to the missing details of the Wannierization parameters in Ref. Jang et al. (2016).

Recently, Hirayama et al. Hirayama et al. (2018) has reported ab initio effective Hamiltonians for bulk cuprates, including LCO, obtained within the FP-LMTO and the cGW method supplemented by the self-interaction correction (SIC) of the Hartree term Hirayama et al. (2013, 2017). The onsite Coulomb parameters of the $E_{g}$ Hamiltonian is reported to be $U_{x^{2}-y^{2}}=$ 5.3–5.5 eV, which is $\sim$25–30% larger than the previous cRPA result of 4.2 eV Miyatani (2014). This enhancement in $U$ can be attributed to the refined treatment of the Coulomb interaction by the GW approximation, which makes the band width of the $d_{x^{2}-y^{2}}$-like orbital smaller and thereby increases the bare $V$ value by $\sim$5%. More importantly, the GW calculation increases the level offset between the $d_{x^{2}-y^{2}}$-like orbital and other orbitals in the $r$ subspace, such as the bonding state formed by the Cu $d_{x^{2}-y^{2}}$ and in-plane O $p_{\sigma}$ orbitals, leading to the considerably weaker screening compared with cRPA. While the cGW-SIC scheme is theoretically more refined than cRPA, it is computationally more demanding than cRPA and its application to the LCO/LSCO interface was impractical. Therefore, cRPA is used both for the bulk and heterostructure in this work, which still gives reasonably accurate results and does not change the conclusions of this paper which mainly focuses on the effects of the structural optimization on the $E_{g}$ Hamiltonians of the LCO/LSCO interface.

Appendix B Averaging method for Coulomb parameters

In the present cRPA calculation, we employ the $8\times 8\times 2$ $\bm{k}$ points to use the tetrahedron method for an accurate numerical integration of Eq. (II). However, we found that the original Eq. (12) gave unphysical negative $V$ values around $\bm{R}=(0,0,c)$ as shown in Fig. 11. The negative contribution around $\bm{R}=(0,0,c)$ comes mostly from the component of $\rho_{mn\bm{R}\bm{q}}(\bm{G})$ around $\bm{q}=(0,0,\frac{\pi}{c})$, which is far larger than the largest positive contributions from $\bm{q}=(\frac{\pi}{4a},0,0)$ due to the prefactor $1/|\bm{q}|$ as well as the large anisotropy of the lattice shape, i.e., $a\ll c$. If one can increase the $\bm{k}$-mesh density up to $N\times N\times 2$ so that $\frac{\pi}{c}\simeq\frac{2\pi}{Na}$ is satisfied, the negative $v$ problem could be solved. However, such a calculation is almost infeasible because $N$ must be as large as 28 even for the smallest $(4,4)$ SL to satisfy the condition.

To mitigate this issue, we simply modify the original $\rho_{mn\bm{R}\bm{q}}(\bm{G})$ given in Eq. (14) as

[TABLE]

where $|\bm{q}|_{\mathrm{rep.}}$ is a representative value of the norm around $\bm{q}$ defined as

[TABLE]

Here, $S_{\bm{q}}$ is the Wigner–Seitz cell of the lattice point $\bm{q}$, and $\Omega_{S_{\bm{q}}}$ is its volume. This treatment assumes that a variation of $\braket{\psi_{m\bm{k}+\bm{q}}^{(w)}}{e^{\mathrm{i}(\bm{q}+\bm{G})\cdot\bm{r}}}{\psi_{n\bm{k}}^{(w)}}$ inside $S_{\bm{q}+\bm{G}}$ is small, which is satisfied when the size of $S_{\bm{q}}$ is reasonably small. If we use $\tilde{\rho}_{mn\bm{R}\bm{q}}(\bm{G})$, we can obtain a correct $1/r$ dependence without negative values as shown in Fig. 11.

S1. Layer dependence of off-site Coulomb and exchange parameters

Figure S1 shows the layer dependence of the off-site Coulomb and exchange parameters obtained with the same outer window setting as the main text, i.e., 8 valence bands per CuO2 layer. In comparison with the on-site Coulomb parameters (Fig. 9 of the main text), the layer dependence of the off-site Coulomb parameters is weak as shown in Fig. S1 (a) and (c). The maximum difference of the $V_{n}$ values between the layers is $\sim$ 0.05 eV, which is smaller than 6% of the average value. By contrast, the layer dependence of the exchange parameters is rather significant, and the maximum difference reaches $\sim$ 0.075 eV, which amounts to 15–19% of the absolute $J$ values. The layer dependence of the exchange parameters are similar to that of the on-site Coulomb parameters.

S2. Parameters with wider outer window

In the main text of the paper, we show the hopping and Coulomb parameters based on the rather extended $d_{z^{2}}$ orbital constructed from a narrower outer window that includes 8 valence bands per CuO2 layer. Since the parameters of the low-energy effective Hamiltonian depend on the shape of the Wannier orbital, whose spread is sensitive to the outer window in La2-xSrxCuO4, we also show the computational results obtained when we included 14 valence bands per CuO2 layer in the outer window.

Table S1 shows the calculated hopping and Coulomb parameters of bulk La2-xSrxCuO4, and Figs. S2 and S3 show the layer-dependent hopping and Coulomb parameters of the LCO/LSCO heterostructures, respectively. It can be inferred by comparing these data with Table 2, Figs. 8 and 9 of the main text that parameters involving the $d_{z^{2}}$ orbital changes dramatically with the change of the outer window. For example, the on-site Coulomb parameter of the bulk N2 system increases from 13.83 eV to 22.51 eV for the bare interaction $v$ and from 3.91 eV to 6.85 eV for $U(0)$. These enhancement can be attributed to the smaller Wannier spread of the $d_{z^{2}}$ orbital constructed from the 14 valence bands than that constructed from the 8 valence bands. The shrinkage of the Wannier orbital occurs due to the inclusion of the bonding state formed by the copper $d_{z^{2}}$ and the apical oxygen $p_{z}$ orbitals in the outer window. By contrast, the parameters that only involve the $d_{x^{2}-y^{2}}$ orbital such as $v_{x^{2}-y^{2}}$ and $U_{x^{2}-y^{2}}(0)$ do not change quantitatively because we employed the frozen window in the both outer-window settings.