## Abstract

Stabilizing superconductivity at high temperatures and elucidating its mechanism have long been major challenges of materials research in condensed matter physics. Meanwhile, recent progress in nanostructuring offers unprecedented possibilities for designing novel functionalities. Above all, thin films of cuprate and iron-based high-temperature superconductors exhibit remarkably better superconducting characteristics (for example, higher critical temperatures) than in the bulk, but the underlying mechanism is still not understood. Solving microscopic models suitable for cuprates, we demonstrate that, at an interface between a Mott insulator and an overdoped nonsuperconducting metal, the superconducting amplitude is always pinned at the optimum achieved in the bulk, independently of the carrier concentration in the metal. This is in contrast to the dome-like dependence in bulk superconductors but consistent with the astonishing independence of the critical temperature from the carrier density *x* observed at the interfaces of La_{2}CuO_{4} and La_{2−x}Sr_{x}CuO_{4}. Furthermore, we identify a self-organization mechanism as responsible for the pinning at the optimum amplitude: An emergent electronic structure induced by interlayer phase separation eludes bulk phase separation and inhomogeneities that would kill superconductivity in the bulk. Thus, interfaces provide an ideal tool to enhance and stabilize superconductivity. This interfacial example opens up further ways of shaping superconductivity by suppressing competing instabilities, with direct perspectives for designing devices.

- Thin films
- copper oxide superconductors
- interface
- variational Monte Carlo method
- first-principles calculation
- phase separation
- inhomogeneity
- self-organization
*d*-wave symmetry- doped Mott insulator

## INTRODUCTION

Thin films and interfaces offer unique platforms for designing materials functions, beyond what is possible in the bulk. Above all, superconductivity at interfaces was observed even in cases where the bulk compounds sandwiching the interface are both nonsuperconducting. Furthermore, critical temperatures *T*_{c} in thin films that are higher than or equal to the maximum *T*_{c} of the bulk material were observed (*1*, *2*). These findings suggest the superiority of interfaces for designing high-*T*_{c} superconductivity.

We do not have a thorough understanding of cuprate high-*T*_{c} superconductors in the bulk, and our knowledge stemming from experimental measurements has been constantly updated since the discovery of superconductivity in the cuprates. However, the dome structure of *T*_{c} as a function of carrier (doping) concentration is a common property irrespective of the compounds. In particular, the optimum *T*_{c} is realized only at a specific doping concentration δ around 0.15 per Cu.

In contrast, at the interface of La_{2}CuO_{4} and La_{2−x}Sr_{x}CuO_{4} (schematically illustrated in Fig. 1A), for *x* in the range of 0.2 to 0.5, stable superconductivity with *T*_{c} ~ 40 K was observed, irrespective of the value of *x* (*2*). This discovery was even more unexpected as the value of this “pinned” *T*_{c} is very close to the maximum value for the bulk superconductivity in La_{2−x}Sr_{x}CuO_{4}, which was realized in the bulk only for the optimum carrier doping concentration of δ = *x* ~ 0.15. If the mechanism of this superiority and stability at the interface is understood, we may gain insight not only into the unsolved mechanism of superconductivity but also into how to reach higher critical temperatures in elaborately designed devices.

For the bulk superconductivity of the cuprates, noteworthy theoretical progress was recently made: Numerical calculations using various tools have been able to reproduce the basic experimental characteristics, particularly the *d*-wave symmetry of the gap and the dome structure (*3*–*11*).

Using cutting-edge variational Monte Carlo simulations (*10*, *12*) for a stacked layer model shown in Fig. 1B (top panel), we show here that superconductivity emerges dominantly at a single layer of the interface between a Mott insulator and an overdoped metal and that its amplitude is independent of the carrier density in the metallic side. The amplitude is indeed pinned at the maximum of the dome structure in the bulk in perfect agreement with the experiment.

Our numerical result shows that this pinning originates from the underlying electronic phase separation in the bulk (*10*, *13*–*15*), which by itself would destroy the superconductivity in the bulk but is now replaced by an interlayer phase separation around the interface instead of the phase separation within a layer (schematically shown in the bottom panels in Fig. 1B). In general, strong coupling superconductivity with a high critical temperature would require a strong effective attraction between electrons. However, this strong attraction works like a “double-edged sword.” Namely, it also makes the system prone to charge inhomogeneity that destroys superconductivity. The interface cleverly eludes this trade-off.

## RESULTS

### Theoretical model

We analyze the multilayer single-band Hubbard model, which is suitable for studying interfaces of the cuprates, defined by(1)where (*c*_{iσν}) is the creation (annihilation) operator of an electron at the *i*th site on the νth layer with spin σ and is the corresponding number operator. For simplicity, we only consider the nearest-neighbor pair 〈*i*, *j*〉 for the intralayer transfer *t*. For the interlayer transfer, we take *t*_{z} = 0.05*t*, and the onsite Coulomb interaction is set to *U* = 8*t*. These are realistic values in terms of first-principles and numerical estimates (*10*, *16*, *17*) compared to the experimental optical gap and transport measurements (*18*, *19*). Hereafter, we set the energy unit *t* = 1 (~0.5 eV in the cuprates). The layer-dependent onsite hole level is represented by ε_{ν}. We confirmed that details of the parameter values do not alter our results.

We perform high-accuracy many-variable variational Monte Carlo (mVMC) calculations at temperature *T* = 0 for an *N*_{s} = *L* × *L* square lattice stacked as a slab with thickness *L*_{layer} (see Materials and Methods for details of the model and method). The mVMC method is capable of describing quantum and spatial fluctuations (*10*, *12*), allowing for an accurate estimate of the superconducting stability among the competing orders.

The experimental interface illustrated in Fig. 1A has a transient region caused by the interlayer diffusion and exchange between La and Sr atoms (blue line in Fig. 1C) (*2*, *20*). To realistically mimic the interlayer diffusion effect that gradually changes the onsite energy level within a few layers, we construct a slab around the interface, with the layer-dependent onsite level as ε_{ν+1} = ε_{ν} − Δε (3 ≥ ν ≥ 1) with a constant Δε (see the red line in Fig. 1C). The 0th layer is assumed to be insulating, and the other layers (ν ≥ 1) become metallic. For the 0th layer, we use ε_{0} = ε_{1} + 1, ensuring the insulating nature. On the other hand, density functional theory calculations for a sharp interface predict a more abrupt change in the onsite energy level (see Fig. 1D and the Supplementary Materials for the first-principles estimate).

Note that the properties around the interface embedded in a sufficiently thick slab with ε_{ν} = ε_{4} for ν ≥ 4 can be well simulated by a slab with a total thickness of *L*_{z} = 5. This is because the transient region near the interface is confined to the region ν < 4 if the Sr concentration is similar to that illustrated in Fig. 1C. The density at ν ≥ 4 converges to a constant corresponding to the bulk value in the overdoped metallic side. In practical calculations, the bulk hole densities at ν = 4 are controlled by changing Δε and the total electron number in the slab.

### Pinning of electron density at the interface

In Fig. 2A, we plot the layer dependence of the hole density δ_{ν} defined by δ_{ν} = 1 − *N*_{ν}/(*L* × *L*), where *N*_{ν} is the average electron number in the νth layer (see Materials and Methods for the method used to determine the charge profile). The bulk hole density, δ_{bulk} = δ_{4}, monotonically increases with increasing Δε. Experimentally, this corresponds to *x* in the metallic side of the interface. Even if δ_{bulk} changes substantially, at the interface, δ_{1} is pinned.

To understand this counterintuitive pinning, we show the calculation for the μ-δ relation of a single layer in Fig. 2B, where μ is the chemical potential (*10*). We essentially find the same μ-δ relation for the uniform bulk (μ = μ_{bulk}, δ = δ_{bulk}) consisting of stacked layers with the same single-particle level coupled by the small interlayer transfer *t*_{z} = 0.05*t*. Nonmonotonic δ dependence of μ leads to a thermodynamic instability with the phase separation for 0 < δ = δ_{bulk} ≲ δ_{PS} ~ 0.20 (see also “Method to determine charge profile” in Materials and Methods). The ε_{4}-δ_{4} relation traces the μ-δ relation by mapping ε_{4} ↔ 4Δεμ and δ_{4} ↔ δ (green curve in Fig. 2A). Remarkably, the μ_{ν}-δ_{ν} and ε_{ν}-δ_{ν} relations at all the νth layers trace the same relation.

This indicates that each layer is well represented by the single-layer model, and the effect of *t*_{z} (= 0.05*t*) is small, as for the μ-δ relation. The main role of *t*_{z} is to distribute the holes across the layers, where neighboring layers simply work as carrier reservoirs. A doping concentration that would lie in the region 0 < δ < δ_{PS} is prohibited at any layer. Consequently, the first layer that would lie in the phase separation region is in reality pinned at the border δ_{PS} as in Fig. 2C. The consequences of the pinning of δ_{1} at δ_{PS} are further discussed later.

### Pinned superconducting order at the interface

To investigate the superconducting properties, we calculate the layer-dependent equal-time superconducting correlations of -wave symmetry, defined as(2)where Δ_{ν,d} denotes the -wave superconducting order parameter at the νth layer, *f*_{d}(** r**) is the form factor that describes the -wave symmetry, and δ

_{i,j}denotes the Kronecker’s delta, with

**= (**

*r**r*

_{x},

*r*

_{y}) being the two-dimensional lattice coordinate scaled by the lattice constants of the square lattice.

In Fig. 3A, we plot *P*_{ν,d}(** r**) for ν = 1 at Δε = 0.2 (blue squares). The superconducting correlation becomes a nonzero constant at the long-ranged part (essentially for ) at the interface layer (ν = 1), implying long-ranged order. Here,

*P*

_{ν,d}(

**) is similar to the value for the uniform bulk system with a hole density similar to that at the interface.**

*r**P*

_{ν,d}(

**) for the bulk system (red circles) is calculated for uniformly stacked layers that have the same single-particle level for all the layers (see Materials and Methods).**

*r*In Fig. 3B, we plot the δ_{bulk} (metallic bulk density) dependence of the superconducting correlations in the long-range limit, which is, in practice, calculated from(3)for sufficiently large *L*, where *M* is the number of vectors satisfying .

In previous work (*10*), this quantity was shown to allow for a practical estimate of the long-range order, and Fig. 3A also supports this criterion. Note that converges to the square of the order parameter (superconducting amplitude) in the thermodynamic limit. We find in Fig. 3B that the squared superconducting order parameter is pinned irrespective of the bulk hole densities, in accordance with the pinned δ_{PS} ≃ 0.20 at the interface. Note that the pinned superconducting amplitude is equal to the maximum value achievable in the stable uniform bulk as we discuss below. This pinning at the maximum is a central result of the present report. The pinned superconducting order parameter is consistent with the anomalous pinning of *T*_{c} independent of δ_{bulk} observed at the interfaces (*2*). It is natural that the same order parameter at *T* = 0 yields the same *T*_{c}.

A question arises on how robust the results are when the atomic interlayer diffusion is absent, where the density functional calculation indicates that the onsite level varies relatively suddenly at the interface (Fig. 1C). We show in the Supplementary Materials that the pinning still exists.

## DISCUSSION

### Relation between intralayer and interlayer phase separations

In the bulk system, the phase separation and enhancement of charge susceptibility near the Mott insulator were first theoretically pointed out (*21*–*23*) and have long been debated (*24*–*33*) in experiments and theories. Even in the simple Hubbard model on the square lattice, the exact solution is not available in terms of the existence of the phase separation and superconductivity [see the study of Misawa and Imada (*10*) for detailed comparisons of the occasionally controversial theoretical results and their accuracies]. Above all, many accurate numerical results have indicated the existence of an extended region of phase separation region for large *U*/*t*. Furthermore, it was shown that the superconducting correlation exhibits its maximum inside the phase separation region, if we allow for metastable states (*10*). However, as a thermodynamically stable state, the maximum emerges at the phase separation border δ = δ_{PS}. The region 0 < δ < δ_{PS} is subtle: Here, the long-ranged Coulomb repulsion ignored in the Hubbard model would lead to a diverging electrostatic energy, if the macroscopic phase separation occurs, which is prohibited in reality. Consequently, the true ground state is replaced by mesoscale inhomogeneous states or long-period charge order, thus compromising the long-ranged Coulomb force as was previously observed (*34*). Even for the Hubbard model in the absence of long-range interactions, stripe-type charge ordering is nearly degenerate with the phase-separated state (*10*, *28*, *29*, *35*). These inhomogeneities introduce pair breaking and suppress superconductivity, allowing for the maximum superconducting order only at the pinpoint of δ ~ δ_{PS}. Even when the stripe (charge density wave) is perfectly ordered and clean, this suppression occurs where the superconductivity at the optimum carrier concentration is connected by the Josephson tunneling through the nonoptimized density region, which has a similarity to continuous superconductor-insulator transition caused by randomness (*36*–*38*).

In contrast, around the interface, our result indicates that the charge inhomogeneity is circumvented by transferring holes between the neighboring metallic layers and the interface to avoid the energy cost due to the intralayer charge inhomogeneity. This transfer violates the charge neutrality of each layer, but the electrostatic energy remains small, because it corresponds to the formation of a capacitor, where the electric field is confined only within the capacitor. This is a remarkable way of eluding the harmful electronic inhomogeneity that is unavoidable in the bulk. At the interface, the inhomogeneity is dissolved into an imbalanced density between the neighboring layers, pinned at both ends of the phase separation region, δ = δ_{PS} and 0. The stable hole density at the interface at δ_{PS} ensures the maximum superconducting amplitude ever realized in the bulk, consistently with the pinned *T*_{c} ~ 40 K (*2*).

We further discuss our intuitive interpretation of why the interlayer phase separation is more stable than a state with intralayer charge inhomogeneities. Although the divergence of the electrostatic energy is avoided even by the introduction of mesoscale inhomogeneities within a layer, the formation of stripes or puddles costs a boundary energy proportional to the length of the domain wall within the layer. From the energy cost due to the domain-wall formation, *E*_{D}, there are two contributions—*E*_{D1} and *E*_{D2}—that will be crucially different between the intralayer and interlayer domain-wall formations.

Because it is in the phase separation region, the energy as a function of doping concentration has a double-well structure, whose two minima are realized at the two phase-separated densities. Forming a domain wall within a layer costs energy *E*_{D1} because in the transient region at the domain wall, the charge density crosses through the maximum in the center of this double-well structure. On the other hand, if the domain wall is located between two layers, this energy cost can be largely avoided, because the charge density can jump from low to high values, owing to the small *t*_{z} and the presence of the intermediate LaO layers.

The other energy cost, *E*_{D2}, is that arising from the spatial dependence of the charge density. Because the Coulomb contribution that arises from the long-ranged part to this spatial dependence is material- and model-dependent, we do not discuss it in detail. A crucial difference between the interlayer and intralayer phase separations for the spatial-dependent part of the energies arises from the kinetic energy: The kinetic energy is clearly lost in the presence of the domain wall because of the carrier confinement in the carrier-rich region. The kinetic energy is dominated by the intralayer hopping contribution; therefore, the domain wall within the layer costs more kinetic energy than the interlayer domain wall.

Our finding offers possible ways for enhancing and stabilizing the superconducting amplitude by making use of the translational symmetry breaking in the interlayer direction. Controlling the carrier density to reach the end point of the phase separation in the bulk is the best way to optimize superconductivity in the uniform bulk. However, it requires careful tuning. At the interface, the situation is much more robust, because the optimal value is automatically reached in a self-organized manner. One can therefore expect easier routes for materials preparation than the careful tuning needed in the bulk.

Furthermore, elucidating the pinning mechanism provides guidelines for the design of materials and devices with enhanced superconductivity: A likely strategy is to attempt interface engineering; an example is to keep the carrier density even in the metastable region inside the phase separation region. Here, the superconducting amplitude would be even larger than that at the end point of the phase separation region.

A related future issue is the mechanism in multilayer superconductors (*39*, *40*). The charge inhomogeneity and the resultant suppression of the superconductivity can be avoided by the external breaking of the translational symmetry as in the case of the interface and the multilayer systems. Because both the iron-based and cuprate superconductors are on the verge of phase separation (*10*, *14*), this strategy may universally apply to materials where high-temperature superconductivity is driven by electron correlation effects.

## MATERIALS AND METHODS

### Numerical methods

To analyze the multilayer Hubbard model in the ground state, we used the mVMC method. Here, we briefly summarize this method. Details can be found in the study of Tahara and Imada (*12*). The variational wave function for the ground state is defined as(4)where and are the Gutzwiller (*41*) and Jastrow (*42*, *43*) factors, respectively (*10*, *12*). These correlation factors are defined aswhere *g*_{ν} and *v*_{ijνμ} are variational parameters. These factors express many-body correlations beyond the mean-field starting point. To restore the symmetry of the Hamiltonian, we used the quantum number projection method (*44*). Here, we used the total spin quantum number projection operator , which restores *SU*(2) spin symmetry with the total spin *S*, where *S* = 0. The one-body part |φ_{pair}〉 is the generalized pairing wave function defined as(5)where *f*_{ijνμ} denotes the variational parameters and 𝒩 represents the total number of electrons. Here, we allowed *f*_{ijνμ} to have a 2 × 2 sublattice structure for each layer (2 × 2 × *L*_{layer} sites exist in the unit cell). We note that the variational wave function |ψ〉 defined in Eq. 4 can flexibly describe different phases, such as the antiferromagnetic, superconducting, and correlated paramagnetic phases. This flexibility is necessary to analyze the multilayer model where the competitions and/or the coexistence of different phases appears. Although the number of variational parameters becomes large to allow flexibility (in this calculation, the number of variational parameters is more than 10^{4}), we simultaneously optimized all the variational parameters by using the stochastic reconfiguration method (*39*, *45*).

In the actual calculations, we took an *L* × *L* × *L*_{layer} lattice (where *L* = 10 and *L*_{layer} = 5) with antiperiodic-periodic boundary conditions in each layer and with open boundary conditions at the two end layers in the direction perpendicular to the layers. The system size was sufficiently large even when one wishes to examine the long-range order of the superconductivity: We confirmed the saturation of the superconducting correlation at long distances when the superconductivity emerged. The obtained superconducting correlations at each layer were close to those obtained for the uniform bulk simulation for the same hole density with that layer. The superconducting correlation of the uniform bulk does not appreciably depend on the thickness of the uniformly stacked layers if the thickness exceeds three layers, although it is slightly smaller than the single-layer result. The small difference originates from the small interlayer hopping *t*_{z}.

### Method to determine charge profile

We define the chemical potential of each layer after taking into account many-body effects as(6)(7)(8)where and () denotes the total energy (electron number) at the νth layer, when the total electron number of the multilayer slab is . Here, we ignored the negligible contribution from the interlayer kinetic energy as we remark later. should be close to to approximate the derivative by the difference in Eq. 6. In the definition of , the site indices *i* and *j* run over the sites contained within the νth layer.

For several choices of Δε, we show δ_{ν} (hole density at the νth layer) dependence of μ_{ν} in Fig. 4A for ν = 4, which was obtained by changing the total electron number in the canonical ensemble of the slab. Here, the hole density and the chemical potential in the bulk layer at ν = 4, δ_{4} and μ_{4}, respectively, have to satisfy the relation between the bulk hole density (δ_{bulk}) and the bulk chemical potential (μ_{bulk}) calculated independently in the uniform bulk system. For the latter, we used the result of the single layer (*10*) because of the periodicity of the bulk and the negligible contribution of *t*_{z}. We separately confirmed that uniformly stacked layers (slab) coupled by *t*_{z} = 0.05*t* did not provide a difference in δ dependence of μ regardless of the layer thickness of the slab. The μ-δ relation is shown as the green solid curve without symbols in Fig. 4A. This posed a constraint wherein the total electron number in the canonical ensemble of the slab was uniquely determined when we fixed Δε. Namely, the point where the doping dependence of μ_{4} crosses with the chemical potential of the bulk (μ_{bulk}) represents the true bulk hole density for a given Δε. For instance, for Δε = 0.2, μ_{4} crosses with the μ_{bulk} around δ_{4} ~ 0.32. We then used δ_{4} ~ 0.32 as the bulk hole density δ_{bulk} for Δε = 0.2. The results shown in the main text were obtained from the calculations that satisfy this constraint. Figure 4 (B and C) shows that the relation between the chemical potential μ and the hole density δ in each layer follows the relation for the uniform bulk, confirming that the interlayer transfer does not change this relation.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/7/e1600664/DC1

Text (sections A and B)

fig. S1. First-principles estimate of electronic structure around a sharp interface.

fig. S2. Chemical potential in bulk metal and saturated superconducting correlation at the interface as functions of bulk hole concentration in the case of a sharp interface.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Funding:**We thank the computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the High Performance Computing Infrastructure System Research projects (hp140215, hp150211, hp150173, and hp160201) supported by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan. This work was also supported by Grant-in-Aid for Scientific Research (16H06345, 16K17746) from the MEXT of Japan. We also thank numerical resources from the Supercomputer Center of the Institute for Solid State Physics at the University of Tokyo. This work was further supported by the European Research Council under its Consolidator Grant scheme (project number 617196) and by the Institut du développement et des ressources en informatique scientifique/Grand Équipement National de Calcul Intensif Orsay under project t2016091393.

**Author contributions:**M.I. envisioned the initial conception and supervised the whole project. T.M. performed the detailed variational Monte Carlo calculations. The results were analyzed and the paper was first written by T.M. and M.I. Y.N. performed first-principles calculation under the supervision of S.B. The article was revised and reflects the contributions of all authors.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2016, The Authors