Research ArticleSUPERCONDUCTORS

Self-optimized superconductivity attainable by interlayer phase separation at cuprate interfaces

See allHide authors and affiliations

Science Advances  29 Jul 2016:
Vol. 2, no. 7, e1600664
DOI: 10.1126/sciadv.1600664

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 La2CuO4 and La2−xSrxCuO4. 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.

Keywords
  • 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 Tc in thin films that are higher than or equal to the maximum Tc of the bulk material were observed (1, 2). These findings suggest the superiority of interfaces for designing high-Tc superconductivity.

We do not have a thorough understanding of cuprate high-Tc 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 Tc as a function of carrier (doping) concentration is a common property irrespective of the compounds. In particular, the optimum Tc is realized only at a specific doping concentration δ around 0.15 per Cu.

In contrast, at the interface of La2CuO4 and La2−xSrxCuO4 (schematically illustrated in Fig. 1A), for x in the range of 0.2 to 0.5, stable superconductivity with Tc ~ 40 K was observed, irrespective of the value of x (2). This discovery was even more unexpected as the value of this “pinned” Tc is very close to the maximum value for the bulk superconductivity in La2−xSrxCuO4, 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.

Fig. 1 Experimental setup and present theoretical model of the cuprate interface.

(A) Schematic experimental setup of the cuprate interface (2). (B) Top: Illustration of the interface model for cuprates. The dotted line denotes the interface between the metallic and the insulating layer. The color schematically illustrates the change in the carrier concentration obtained in the present work. Bottom: Two hypothetical bulk or single-layer phases with charge inhomogeneity within a layer. (C) Layer dependence of onsite energy level chosen to model the interface (red line). In the metallic phase, the onsite energy level is assumed to change linearly. This is an approximation that takes into account the effect of interlayer atomic diffusion [blue curve; taken from the study of Logvenov et al. (20)] combined with effects from the Madelung potential and the spatial extension of the Wannier orbital at the interface. (D) Onsite level of a sharp interface modeled by means of an ab initio calculation for x = 0.4.

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 (311).

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, 1315), 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 byEmbedded Image(1)where Embedded Image (ciσν) is the creation (annihilation) operator of an electron at the ith site on the νth layer with spin σ and Embedded Image 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 tz = 0.05t, and the onsite Coulomb interaction is set to U = 8t. 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 Ns = L × L square lattice stacked as a slab with thickness Llayer (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 Lz = 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.

Fig. 2 Layer dependence of doping concentration around the interface.

(A) Layer and level-slope dependence of carrier density (filled circles and blue surface). At the fourth layer, the green curve is taken from the μbulkbulk relation, and the two horizontal gray sheets show the phase separation boundaries determined in (B). Note that μbulk = μ4 ~ ε4 − 2.4 is satisfied, indicating that the grand canonical ensemble is realized for ν = 4. The phase separation region in the bulk is also evaded around the interface in any layer ν. In contrast, the noninteracting case with the same δ4 plotted for Δε = 0.1 (red line) enters the present phase separation region. (B) Relation between the hole density δbulk = δ and the chemical potential μbulk = μ in the uniform bulk system calculated within the canonical ensemble for a single layer representative of the bulk (10). The Maxwell construction (dashed line) determines the phase separation as the gray region between δbulk ~ 0.2 and 0. (C) Hole density at interfaces δ1 shows pinning against bulk hole density δbulk.

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 tz = 0.05t. 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 ε44 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 tz (= 0.05t) is small, as for the μ-δ relation. The main role of tz 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 Embedded Image-wave symmetry, defined asEmbedded Image(2)where Δν,d denotes the Embedded Image-wave superconducting order parameter at the νth layer, fd(r) is the form factor that describes the Embedded Image-wave symmetry, and δi,j denotes the Kronecker’s delta, with r = (rx,ry) 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 Embedded Image) at the interface layer (ν = 1), implying long-ranged order. Here, Pν,d(r) is similar to the value for the uniform bulk system with a hole density similar to that at the interface. Pν,d(r) 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).

Fig. 3 Superconducting correlations and amplitudes.

(A) Spatial dependence of d-wave superconducting correlations at the interface (ν = 1) for Δε = 0.2 and δbulk ~ 0.32 (blue squares) compared with that of the uniform bulk for a hole density similar to that at the interface (~0.20). The red circles are obtained for the bulk (stacked layers) with uniform chemical potential. The saturation at long distances r indicates long-range order. The data sets are for the linear size in the plane direction, L = 14, for which we confirmed convergence to the thermodynamic limit. (B) Bulk hole density dependence of squared superconducting amplitude at the interface (ν = 1) defined by Embedded Image. Embedded Image hardly depends on the bulk hole densities. (C) Layer dependence of Embedded Image. This function is strongly peaked at the interface ν = 1.

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 fromEmbedded Image(3)for sufficiently large L, where M is the number of vectors satisfying Embedded Image.

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 Embedded Image converges to the square of the order parameter (superconducting amplitude) Embedded Image 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 Tc independent of δbulk observed at the interfaces (2). It is natural that the same order parameter at T = 0 yields the same Tc.

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 (2123) and have long been debated (2433) 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 (3638).

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 Tc ~ 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, ED, there are two contributions—ED1 and ED2—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 ED1 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 tz and the presence of the intermediate LaO layers.

The other energy cost, ED2, 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 asEmbedded Image(4)where Embedded Image and Embedded Image are the Gutzwiller (41) and Jastrow (42, 43) factors, respectively (10, 12). These correlation factors are defined asEmbedded ImageEmbedded Imagewhere gν and vijνμ 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 Embedded Image, 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 asEmbedded Image(5)where fijνμ denotes the variational parameters and 𝒩 represents the total number of electrons. Here, we allowed fijνμ to have a 2 × 2 sublattice structure for each layer (2 × 2 × Llayer 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 104), we simultaneously optimized all the variational parameters by using the stochastic reconfiguration method (39, 45).

In the actual calculations, we took an L × L × Llayer lattice (where L = 10 and Llayer = 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 tz.

Method to determine charge profile

We define the chemical potential of each layer after taking into account many-body effects asEmbedded Image(6)Embedded Image(7)Embedded Image(8)where Embedded Image and Embedded Image (Embedded Image) denotes the total energy (electron number) at the νth layer, when the total electron number of the multilayer slab is Embedded Image. Here, we ignored the negligible contribution from the interlayer kinetic energy as we remark later. Embedded Image should be close to Embedded Image to approximate the derivative by the difference in Eq. 6. In the definition of Embedded Image, 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 tz. We separately confirmed that uniformly stacked layers (slab) coupled by tz = 0.05t 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.

Fig. 4 Relation between chemical potential and hole concentration.

(A) Chemical potentials μ4 (determined from Eq. 6) as a function of the hole density δ4 at the fourth layer for several choices of Δε are plotted as curves with symbols. Here, δν is defined as Embedded Image. For a choice of Δε, μ4 curves are drawn by changing the total electron number in the whole slab. We assume that μ4 converges to the bulk chemical potential μbulk (green curve), which was calculated in the study of Misawa and Imada (10). Therefore, the realistic bulk hole density δbulk is determined from the crossing point between the bulk chemical potential μbulk given by the green curve and μ4 for each choice of Δε. Cases with different δbulk are obtained from different Δε. The nonmonotonic behavior of the green curve signals the existence of a phase separation region. A Maxwell construction (horizontal dashed line) allows us to determine the coexistence region as 0 < δ < δPS ~ 0.2 (gray area). (B) Chemical potentials μν (determined from Eq. 6) as a function of the hole density δν for ν = 1 to 4 for several choices of Δε are plotted as symbols. It follows the bulk behavior shown by the green curve, indicating that each layer behaves as a single layer (or uniform bulk) in the μ-δ relation with negligible effects from tz. (C) Chemical potential difference μν − μPS plotted as a function of the onsite-level difference εν − εPS. The straight bold line shows that the chemical potential at each layer shifts in accordance with the shift of the onsite level, again indicating that the effects of tz is negligible and each layer behaves as the grand canonical ensemble with the hole onsite level εν.

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.

References (4651)

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.
View Abstract

Navigate This Article