Examination of permittivity for depolarization field of ferroelectric by ab initio calculation, suggesting hidden mechanisms

Electrostatics of depolarization field Ed in relation to the polarization is studied. In particular, the value of permittivity for Ed (εd) in prototypical situations of ferroelectrics, including Mehta formula, is examined by ab initio calculations. By using spontaneous polarization PS corresponding to accurate experiment ones, we show εd = 1, which suggests that the results of εd ≫ 1 indicate hidden mechanisms; εd = 1 suggests that the effect of Ed is significant to induce intriguing important phenomena overlooked by εd ≫ 1. A bridge between εd = 1 and εd ≫ 1, i.e. the consistency of εd = 1 with conventional results is presented. The exact electrostatic equality of head-to-head–tail-to-tail domains to free-standing ferroelectrics is deduced. Hence, most stoichiometric clean freestanding monodomain ferroelectrics and head-to-head–tail-to-tail domains are shown unstable regardless of size, unless partially metallic. This verifies the previous results in a transparent manner. This conclusion is shown consistent with a recent hyperferroelectric LiBeSb and “freestanding” monolayer ferroelectrics, of which origin is suggested to be adsorbates. In addition, this restriction is suggested to break in externally strained ultrathin ferroelectrics. The macroscopic formulas of Ed are found valid down to a several unit-cells, when electronic and atomic-scale surface effects are unimportant and accurate PS is used.

thin films by using E d = − P S /ε f ε 0 with P S from remnant polarization measurements (7 μC/cm 2 ) and ε d = ε f = 10. In the electrostatic study of Tian et al. 7 , E d in BiFeO 3 thin films was estimated by using ε f = 60 (in Eq. (1) of Ref. 7 ), where this P is initially a total polarization of a single domain state. Kim et al. 8 estimated E d in BaTiO 3 (BTO) ultrathin films through the formula by Mehta 4 with the experimental remnant P obtained from pulse train methods and ε d = ε f = 80 from ε f − E ext curves. In the Ginzburg-Landau-Devonshire (GLD) theory of BTO ultrathin films by Jo et al. 9 , E d was given through the Mehta formula 4 with the same P as P in the GLD equation and ε d = ε f = 80 of Ref. 8 , where P in GLD theory is the total polarization. Schroeder et al. 10 estimated E d in HfO 2 and PZT ultrathin films through the Mehta formula 4 with experimental P S and ε d = ε f = 20-300. Similar analyses with ε d = ε f ≫ 1 are frequently employed 16,17 .
Contrastingly, a primitive considerations show ε d = 1 for P S , i.e. a total P S (≡ P(E f (E ext = 0))) 13 , where E f is the total macroscopic electric filed in FE. For freestanding FEs, for example, E f (E ext ) = E ext − P(E f (E ext ))/ε 0 or E f (E ext ) = E ext /ε f − P(E f (E ext = 0))/ε 0 ("Results"). This implies that ε d = ε f 4-10,16,17 may be double counting, while we note that the electrical measurements of P S 4-10,16,17 are indirect measurements based on induced charge per area in electrodes Q. If ε d = 1 is correct, the successes of the analyses using ε d = ε f are attributed to inappropriate parameters or unidentified screening mechanisms.
We think that the existence of this controversy on ε d is due to explanations based on macroscopic quantities. Because macroscopic explanations are abstract, they are unsuitable to bridge the gap between two conflicting views of ε d . On the other hand, ab initio estimation of ε d is considered as the clearest method for this problem but is not reported to our knowledge. Hence, we clarify ε d in the formulas of E d , by ab initio simulations in which ab initio P S is exactly P(E f (E ext = 0)), which is considered as P S obtained by accurate experiments of ideal samples. Here, the standard theoretical assumptions: pure, stoichiometric, clean FEs are used. E d is related to fundamental issues such as stability of monodomains, critical thicknesses of FE, and the emergence of ferroelectricity in superlattices. Some of these subjects require the consideration of other effects such as strain-induced FE and electronic effect at electrodes 11 . To avoid the complexity, we concentrate on free-standing insulating FE and its electrostatic identicals, i.e. head-to-head-tail-to-tail (HH-TT) domains. Thus, we estimate the value of ε d in the formula of E d in a clear manner.
As expected from the electronic interaction at the electrode 11 , it may be argued that the formula of E d based on electrostatics is not possible for nm-FEs. We resolve this by focusing on the formulas of E d and using ab initio P S in the formulas. Therefore, P S in these formulas contains the effects of the interactions in slabs or superlattices, whereas the absence of these effects in conventional studies has limited the applicability. The use of nonlinear ε f is often better than linear ε f but can be approximated by an average linear ε f 12 . Therefore, the conclusions below are applicable also to the nonlinear ε f (≫ 1).

Model
For simplicity, we discuss 1-dimensional (1D) cases with E ext = 0, where FEs with thickness l f are in periodic slabs (Fig. 1). Here, 1D refers not to the shape of object ( Fig. 2) but to the case where all the properties change only along one coordinate; Fig. 1a,b show FE in vacuum with thickness l V and FE/paraelectric (I adj ) superlattice, respectively, while the latter mimics an inhomogeneous FE. I adj stands for both vacuum and insulator, which is dielectric or FE having different P S. The polarization, thickness, and permittivity of I adj are P I , l I , and ε I , respectively, and the thickness of slab is l SC = l f + l I (l V ). The angles of the polarization P of FE and I adj to the slab direction are θ and θ I , respectively.
The macroscopic and atomic electrostatic potential (ϕ) of these models are represented by Fig. 1c,d, respectively. ab initio E d (E d ab initio ) was obtained from the envelope of the peak tops of atomic electrostatic potential, of which example is Fig. 1d. All these FE/vacuum and FE/paraelectric exhibited the density of states (DOS) of insulators (Fig. 1e,f). Additionally, BaTiO 3 (BTO) capacitors are examined, where metal layers are standard electrode materials for FEs: SrRuO 3 or Pt and ~ 20 Å (Fig. 1g).
Accurate estimation of P S is indispensable for estimating E d correctly and achieved by direct Berry phase calculations. To enable these calculations, we designed special FE slabs and procedures described in "Methods". This is because stable 1D-FEs in vacuum are metallic 14,15 and, hence, direct Berry phase calculations are not possible; Even a two unit-cell thick (~ 8 Å) BTO in vacuum is metallic, when FE is enforced 18 . To achieve insulativity, we used tetragonal (P4mm) SrTiO 3 (STO) of which a-axis lattice constant increased by 0.5% and decreased by 0.01% from that of the theoretical cubic phase. For these a-axis lattice constants, bulk STO was FE 19 . We refer to them as STO1.005 and STO.9999, respectively, of which bulk P S 's were 3.56 μC/cm 2 and 6.15 μC/ cm 2 , respectively, by VASP [19][20][21][22][23][24][25][26][27][28] .
Macroscopic equations of E d are obtained in a following manner. The normal component of P of FE (P ⊥ ) under E d is P ⊥ = P S cos θ + (ε d − 1) ε 0 E d in standard approaches [4][5][6][7][8][9][10]16,17 . The equation of continuity of electric flux is P S cos θ + ε d ε 0 E d = P I cos θ I + ε 0 E I , where E I is the macroscopic electric field in I adj . The validity of this continuity in the presence of peaks at the surfaces (Fig. 1d) is explained in "Methods". The continuity of potential in a periodic boundary condition yields E I l I = − E d l f (Fig. 1c). Therefore, we have for θ = 0 For θ ≠ 0 and θ I ≠ 0, Eq. (1) is E d = − (P S cos θ − P I cos θ I )/ε 0 (ε d + l f /l I ). In the present study, P S and P I in Eq. (1) are given by ab initio calculations that simultaneously yield E d consistent with P S . Therefore, the only unknown quantity is ε d .
When P I = (ε I − 1)ε 0 E I , Eq. (2) for θ = 0 is Equation (2) yields also E d in FE capacitors, because equations of continuity of electric flux similar to the above hold; A short-circuited FE capacitor is modelled as a perfect-metal/insulator(l I /2)/FE/insulator(l I /2)/ perfect-metal 29 , where the perfect metal refers to a metal with zero screening length and the thickness of each screening layer λ is l I /2. Assuming P I = ε I ε 0 E I in screening layer, Eq. (2) is applicable and yields the Mehta formula 4 with ε d = ε f and l I = 2λ. For FE capacitors with θ ≠ 0, the formula beneath Eq. (2) is applicable. Because we neglected the electronic interactions at the metal/FE interface of 1-2 unit-cell, e.g. quantum mechanical smearing 30 , the formula for capacitors may be inaccurate for l f < several unit-cells.
The nominal FE thickness l f is the distance between the center position of a top ion and that of a bottom ion, but twice of the atom radius ~ 0.5 Å × 2 should be added. In case of FE/vacuum, this correction was examined by considering the smear-out of electrons into vacuum 29 ; When λ smear (~ 0.8 Å) is the distance between an outermost electron density and a center of ion position ("Methods"), FE thickness appropriate for electrostatics is l f eff = l f + 2λ smear . As seen below, Eqs. (1) and (2) can be valid down to a few nm l f in case of FE/vacuum. Additionally, surface buckling layer is electrostatically a dipole layer, and its net charge is zero. Therefore, even in the presence of buckling layer, these formulas for 1D are also valid, by regarding buckling layer as dead layer ("Methods"); The effective thickness is l f eff − 2l buckle , where l buckle is the thickness of a buckling layer ~ 1-2 unit-cells.

Results
Estimation of ε d . ε d was examined through the comparison of ab initio E d with E d of Eq. (1) or (2) that uses different values of ε d . For FE/vacuum (P I = 0), P S in Eq. (1) was the rigorously calculated P S of the slab by Berry phase. For FE/insulators and FE capacitors, P S and P I in Eq. (2) were calculated ab initio. Therefore, all the parameters in Eqs. (1) and (2) except for ε d are accurately known. Figure 3 compares E d 's by Eq. (1) with E d ab initio 's, where Eq. (1) uses ε d = 1, 4. Here, ε d = 4 is the lower bound of electronic permittivity of STO 29 . E d ab initio is accurate for a long l f , and data for l V ≫ l f reflects the effect of ε d explicitly because of E d ~ − P S /ε 0 ε d . Large symbols in Fig. 3 show the data points satisfying both conditions and, hence, are important.
In Fig. 3, Eq. (1) with ε d = 1 agrees with E d ab initio within 10% always for λ smear = 0.8 Å and mostly for λ smear = c/2. The difference between E d 's for λ smear = 0.8 Å and = c/2 provides typical error bar and is approximately 10%. Equation (1) with ε d = 4 deviates from E d ab initio 's by more than 140%, and the deviations increase monotonically with ε d . Additionally, Eq. (1) yields the potential difference ϕ max =|E d |l f = P S /ε 0 (ε d /l f + 1/l V ), which, with ε d = 1, quantitatively agrees with bandgap Eg that decreases with l f and l V in Fig. 1e.
For BTO/STO, Fig. 4  . In particular, the agreements are within ± 6% in the results by PBE + U (a method of ab initio calculation ("Methods")). Equation (2) with ε d = 20 deviated by more than 1000% from E d ab initio . The deviation increased monotonically with ε d , whereas ε f > 20 is usual for inorganic FEs 4,5,7-10,16,17 . For capacitors, Fig. 5a Bridge between ε d = 1 and ε f 13 . We showed ε d = 1 and will use it below. By noting that electrical measurements of P S are based on the change of charge per area in electrodes ΔQ induced by E ext , a bridge between ε d = 1 and ε d = ε f will be shown. Because most studies [4][5][6][7][8][9][10]16,17 are for FE capacitors, the followings are for 1D FE capacitors with θ = 0, which can represent FE/vacuum for ε I = 1 and freestanding FE for l I /ε I ≫ l f . Figure 2. Definition of (a) three, (b) two, (c) one-dimensions for FE in this article. Dimensionality is not referred to the shape of an object. θ is the angle between the normal to the surface in (c) and the direction of P S . (d) Typical measurement of P S in a capacitor. This is smaller than bulk P S (E d = 0), because of nonzero screening length in electrodes.   2) is estimated with λ smear = 0.8 Å and λ smear = c/2, respectively. The shape of symbol indicates a slab structure: For STO1.005, orange circles, red triangles, light blue diamonds, and blue inverted triangles correspond to 7-unit-cell-STO with l V = 30 Å, 7-unit-cell-STO with l V = 200 Å, 10-unit-cell-STO with l V = 30 Å, and 10-unit-cell-STO with l V = 100 Å, respectively. For STO0.9999, red squares, green pentagons, and blue 90°-rotated triangles correspond to 5-unit-cell-STO with l V = 300 Å, 6-unit-cell-STO with l V = 30 Å, and 7-unitcell-STO with l V = 30 Å, respectively. www.nature.com/scientificreports/ Using the total electric field in FE E f and the total polarization of FE P, 13 . E I , the field in the screening layer of the electrode, is , suggesting that the measured permittivity is ε f .
We show an example: Additionally, Eq. (2) with ε d = 1 shows D = P S − P S /(1 + ε I l f /l I ), while D = Q. Because l I /ε I ~ 0.1 Å, the difference between the real P S and the measured P S is detectable only for l f < 10 Å. As for potential difference, Eq. (2) is well approximated by E d = − P S l I /(ε I ε 0 l f ) for l f > 10 Å, because l I /ε I is short ~ 0.1 Å. Therefore, the potential difference across the capacitor is independent of the FE thickness l f , when the quality of FE is independent of l f and FE is ideally stoichiometric.
Permittivity for non-polarization field (built-in field). Because the polarization P in standard GLD theories 9,31-34 are formulated with a single total polarization, ε d = 1 should be used in standard GLD theories.
The preceding results have shown that the permittivity that expresses the change of P in response to E ext is ε f (≫ 1) even for FE under E d . By the same logic, the change of P by built-in internal field E bi is also expressed by ε f (≫ 1), where E bi is not due to P or a dipole that is not expressed by P. E bi exists in FEs by various mechanisms such as the diffusion potentials at pn and Schottky junctions and chemical orders, e.g. LaAlO 3 in the polar catastrophe model.
For example, P S = 0 and E d = 0 in a bulk cubic BTO. However, E bi ≠ 0, when the surfaces of a cubic BTO slab are asymmetrically terminated to form a dipole, e.g., BaO/TiO 2 /BaO/…/TiO 2 /BaO/TiO 2 . Hence, to achieve E bi = 0, the present study used chemically symmetric slabs (Fig. 1a,b,g), e.g. BaO/TiO 2 /BaO/…/TiO 2 /BaO. Insulativity condition. For 1D-FE to remain insulating without artifactual screening, el f |E d | < Eg (e: elementary charge), for which Eq. (2) and ε d = 1 yield 1/l f > eP S cos θ/ε 0 Eg − ε I /l I . Therefore, the condition of insulativity is one of the followings where Eg* and P S *(T C ) are bulk Eg normalized by 2 eV and P S of bulk FE at T C normalized by 10 μC/cm 2 , respectively, and the unit of l f and l V is Å. P S of bulk FE at T C approximates the critical P S of FE that is about to become paraelectric by E d 12 .
Equations (3) and (4)   . This suggests that freestanding FEs with normal bulk properties are FEs with metallic layer or insulating paraelectrics 12,15 as explained by the following GLD analysis. This conclusion is valid also for HH-TT domains with θ = 0. Standard GLD theories are based on a single polarization vector P as the order parameter. We approximate the polarization possibly missed in such GLD theory 29 by an extra permittivity ε NG − 1 35,36 , while ε NG is speculatively close to electronic permittivity 29 . The GLD energy F of an insulating FE is F = (T − T 0 )P 2 /2Cε 0 + βP 4 /4 + γP 6 /4 − PE d /2, where T 0 , C, β, γ, and θ are Curie-Weiss temperature, Curie constant, and GLD coefficients, respectively. The effect of strain can be incorporated in T 0 and β [31][32][33][34]37 . Curie temperature T C is T 0 + ΔT, where ΔT = 3β 2 /16γCε 0 . For 2nd order transition, γ = 0, β > 0, and T C = T 0 . The effect of E d = − P S /ε 0 ε NG is the change of T 0 to T 0 − C/ε NG in F, where Eq. (2) with ε d = ε NG and l I = ∞ is used.
These ε f 's of FEs undergoing 2nd and 1st transitions appear too small for experimentally observed bulk metal oxide FEs. Therefore, freestanding insulating FEs satisfying C < ε NG T C are unlikely to exist, unless ε NG is far larger than 5; That is, for freestanding FE materials, there exists virtually one choice between a partial loss of insulativity and a loss of FE.

Design of freestanding insulating FE.
For freestanding insulating FE (θ = 0), C < ε NG T C was shown, while Δϕ ~ P S (E d )l f /ε 0 . Therefore, FE materials having a very large ε NG (≫ 5) can retain FE and remain insulating, when ultrathin. Such FE materials are unlikely to exist. Alternatively, we may consider electrically freestanding FE or FE with clean surface that is not mechanically freestanding. In this case, T 0 (~ T C ) of common FEs increases to T 0 eff by inplane strain, while ab initio calculations shows that T 0 eff is much larger than those of standard GLD theories 37 . Therefore, heavily strained FE materials may retain FE and remain insulating, when ultrathin (Formula estimating an effective T 0 from ab initio P S is Ref. [88] of Ref. 37 ). The above indicates that ε f of such FE is extremely low for E d = 0 but can be large for E d ≠ 0, because the coefficient of the first term GLD energy F is (T − T 0 eff + C/ε NG )/2Cε 0 (freestanding).

LiBeSb. LiBeSb with P6 3 mc symmetry 38 is reported as a hyper-FE that retains both FE and insulativity in FE/
paraelectric superlattices, which may contradict the above conclusion on insulativity and FE stability. Therefore, we ab initio calculated one-unit-cell LiBeSb (l f = 6.08 Å) in vacuum (Fig. 6a). Figure 6b shows a metallic DOS of LiBeSb for l V = 31.7 Å, while metallicity increases with l V . This is consistent with the above conclusion and the previous reports 14,15,[39][40][41] . Equations (3) and (4) explain the insulativity of LiBeSb 38 as the effect of adjacent dielectric. Actually, similar to LiBeSb, BTO/STO superlattices are insulating as in Fig. 1f, while BTO/vacuum is partially metallic 18 .
Freestanding and free-surface FE: hidden mechanism. Mechanically freestanding FE is customarily referred to as freestanding; Ji et al. 42 reported exceptionally intriguing results of the freestanding insulating BiFeO 3 (BFO) that retains FE down to monolayer. This appears to contradict both the reports of metallicity at HH-TT domains of BFO and the present results, esp. the single choice between insulating paraelectric and partially metallic FE. If ε d = ε f = 100, the potential difference Δϕ across freestanding insulating BFO of 1-4 unit-cell thickness with a moderate P S ~ 20 μC/cm 2 is 0.09-0.36 V by Eq. (2) with l I = ∞, which allows this BFO to be insulating in agreement with Ji et al. 42 . For ε d = 1, Δϕ increases by 100 times, by which BFO's have to be partially metallic. www.nature.com/scientificreports/ Hence, we shall look at the measurements of Ji et al. 42 . For freestanding FE, it was shown that the surface or boundary of 1 ~ 2 unit-cell thickness was metal and the rest was insulating FE 14,15,39,40 . So, the metallicity is detectable only by inplane conductance, which is absent in Ji et al. 42 . Second, because the crystallographic properties of FE with these metal layers was shown to be those of FE 14,15,39,40 , the crystallographic measurements of Ji et al. 42 do not exclude metallic layers. Third, because piezoelectric measurements use bottom and top electrode (or tip) and may move ions 43 , those by Ji et al. 42 are not that of freestanding FE. Consequently, all the measurements of Ji et al. 42 do not contradict the conclusion of the present paper.
More importantly, the interdisciplinarity of nano FE hides true mechanisms. In the present case, "freestanding" is defined by surface science and electrostatics. For example, Fong et al. found monodomain FE of 3 unitcell thickness as opposed to E d -limited domain and size effect, which was later attributed to adsorbates 44 . This agrees with recent ab initio study 45 . Further, photoemission spectroscopy in UHV showed that SrTiO 3 surface was covered by adsorbates even in ultrahigh vacuum (UHV) 46 . Actually, the free surface with P S ⊥ surface is insulator-like in air and metallic in UHV when cleaned 14 . Because the insulating freestanding FE 42 was exposed to air and water, we suggest adsorbates as its hidden mechanism.

Conclusion
We studied the electrostatics of E d , especially, the value of permittivity ε d in the formula of E d by ab initio simulations, where ab initio P S corresponded accurately to experimental P S . For this, the standard theoretical assumptions: pure, insulating, stoichiometric, and clean FEs were used. To validate the analyses of E d based on electrostatics, we concentrated on the formulas of E d for accurate ab initio total P(E f (E ext = 0)) that contained various atomic effects and corresponded to experimental P S . Further, we focused on the simplest cases of E d : freestanding 1D-FE, HH-TT domains, and superlattices that mimicked inhomogeneous FE and FE/dielectric. The present ab initio simulations showed ε d = 1 ± 0.06-1 ± 0.2. That is, ε d = 1 should be applied to experimental and standard-GLD P S 's. A contradiction between ε d = 1 and ε d = ε f was resolved by a bridge; Even under E d , the permittivity for E ext and built-in field E bi was ε f . Therefore, if a study requires ε d ≫ 1 4-10,16,17 , the value of P S is incorrect, the values of the parameters are inappropriate, or, most likely, hidden screening mechanisms exist 14,15,[43][44][45][46][47][48] .
For freestanding insulating FEs (l I = ∞), Eq. (2) yields E d = − P S /ε 0 (or E d = − P S cos θ/ε 0 ), while, for HH-TT insulating domains, Eq. (1) with P I = − P S and l I = l f yields the same E d . Therefore, when the effects at surface of 1-2 unit-cell is unimportant, freestanding FEs are electrostatically exactly identical with HH-TT domains.
Consequently, both the electrostatic energy of E d and the FE free energy of insulating freestanding and HH-TT FEs scale linearly with l f . This implies that the stability of 1D-freestanding and HH-TT insulating FEs is independent of size 12,15 , when the energy increase by surface effect and domain walls energy is ignored. A strain effect to overcome this restriction was suggested.
The electrostatic formulas of E d (Eqs. 1 and 2) were valid down to a several unit-cell scale (Figs. 3, 4, 5), when atomic-scale surface effects, e.g. interactions with electrodes 11,28 were unimportant. Even with buckling at FE surfaces, these formulas can be valid by regarding buckling layers as dead layer. presented, because they have Eg wider than Eg of the TiO 2 -terminated STO slabs (Fig. 1e). To enforce FE, the ion-positions in the STO/vacuum slabs were not optimized, because, otherwise, FE disappears (Fig. 7). Therefore, the STO unit-cells in the slabs retained the ion positions of bulk STO1.005 or STO.9999. These calculations of STO/vacuum were only for the examination of E d and ε d and do not correspond to standard experiments. P S of STO1.005 in the slab was typically 1 μC/cm 2 . The models of FE/paraelectric are BTO/STO superlattices. All the calculated forces were < 1 meV/Å after geometry relaxation, and the a-axis lattice constant of STO was expanded by 1.1-1.3%. The bulk STOs that had these a-axis lattice constants were paraelectric 19 . The a-lattice constant of BTO/SrRuO 3 and BTO/Pt capacitor was fixed at the theoretical a-axis lattice constant of cubic STO and bulk tetragonal BTO, respectively, and all other ion positions were relaxed (Fig. 1g). The atomic models of BTO/SrRuO 3 are similar to BTO/STO (Fig. 1b). The use of the theoretical a-lattice constant of cubic STO corresponds to the thin films on STO substrates. The surfaces of the BTO and SrRuO 3 were TiO 2 and SrO, respectively, and Pt atoms at the interface aligned with O atoms of TiO 2 plane.
The present study is about the formulas of E d for given structure parameters. Here, the change of P S by the interactions in the slabs is included consistently in these formulas by the use of ab initio P S in these formulas.
The ab initio calculations with VASP 20-22 used the projector augmented wave method 23 with a Monkhorst-Pack 24 mesh of 8 × 8 × 2 for slabs and an energy cutoff of 650 eV. PBEsol functional 25 was used, unless otherwise mentioned. Ab initio P S was calculated through Berry phase 26 . The results of BTO/STO were reexamined with PBE functional 27 with Hubbard U (PBE + U) 28 , which was used also for BTO/Pt. In the slab calculations, graphic processing units acceleration 49 (1) and (2) using ab initio P S with E d for the ion positions same as those of this P S , the accuracy of P S for given ion positions matters.
Berry phase calculation of P S for given ion positions is accurate but only possible for insulators. For example, the present Berry phase calculations yields P S of bulk BaTiO 3 that agree with experimental P S within 4%, when experimental ion positions and lattice constants of at 303 K are used 37 .
Therefore, to obtain accurate P S , the dipole moment of a whole slab was calculated with Berry phase; We treated these slabs as unit-cells to apply Berry phase calculations directly, unlike conventional approaches. P S was obtained by dividing the dipole moment by the volume of FE part of the slab. These P S 's were referred to "rigorously calculated P S 's of the slab" and obtained for all the FE/vacuum and BTO/STO slabs. Here, STO1.005 and STO.9999 slabs are insulating, allowing accurate Berry phase calculations.
Additionally, P S of the unit-cell that possessed exactly the same ion positions as those in the slab was calculated with Berry phase and, then corrected with atomic polarization by E d by the procedures in Ref. 29 . These P S 's agree perfectly with "rigorously calculated P S 's of the slab", which further confirmed the accuracy of the present P S 's of FE/vacuum and BTO/STO. These corrected P S 's 29 were used for capacitors. Therefore, in the present study, P S 's are accurate total P S 's and self-consistent with E d . Hence, P S 's in Figs. 3, 4 and 5 are accurate. (1) and (2) are applicable to the regions much larger than unit-cell. Here, the peaks of 1.5 Å width at the surface in Fig. 1d may be suspected to invalidate Eqs. (1) and (2). These peaks are due to effective surface dipoles caused by electron tunneling smearout; Surface electrons smear out in vacuum, making positive charge density inside the surface and negative charge density in vacuum.

Validity of continuity of electric flux and surface dipoles. Equations
The heights and shape of the two peaks from the baseline (yellow line in Fig. 1d) are the same (1.53 V). This means σ + R = σ + L and σ − R = σ − L as expected from their origin, where σ + R , σ − R , σ + L , and σ − L are positive and negative charge densities that yield the right and left peak, respectively. Because of the charge neutrality of FE, σ + R + σ + L + σ − R + σ − L = 0, i.e. σ − R + σ + R = 0. Therefore, the continuity of the electric fluxes D FE in FE and D I in I adj (Fig. 1b) is D FE − D I = σ − R + σ + R = 0, i.e. the continuity of electric flux D FE = D I , where D FE = P S cos θ + ε d ε 0 E d and D I = P I cos θ I + ε 0 E I . A clearer example is shown in Fig. 7, which evidently shows the continuity of electric flux and, hence, validates the use of Eqs. (1) and (2).
Because surface buckling in vacuum is electrostatically dipole due to ion displacements, the arguments exactly the same as the above hold. Therefore, the electric flux of the inside D FE and the outside D V of the buckling layer is continuous (D FE = D V ).
Effective l f (l f eff ). For FE/vacuum, the effective l f (l f eff ) was estimated from the planer averaged electron density ρ profiles 29 . Below, z = 0 corresponds to the position of bottom ion. Because ρ at z = − 0.8 Å was same as the minimum ρ of inner part in all the ρ-z curves, the region of z = 0 ~ − 0.8 Å was considered as a part of FE (λ smear = 0.8 Å), and l f eff was l f eff = l f + 2λ smear . In addition, λ smear = c/2 (half unit-cell) was also tested, and l V = l SC − l f eff . For BTO/STO, l f was defined as the distance between the top and bottom Ti ions of BTO (Fig. 1b), and l I = l SC − l f . For FE capacitors, l f eff was l T-B − uc BTO (outermost ions are Ti, uc BTO : length of a BTO unit-cell), for which the quantum mechanical smearing 29 may be responsible. The estimations with l T-B − 1.5uc BTO were also tested. The effective thicknesses of the screening layer, i.e. the effective passive layer l I /2ε I of BTO/SrRuO 3 and BTO/Pt were estimated as 0.1 Å and approximately 0.05 Å, respectively 29

Data availability
The data required to reproduce these findings can be provided upon reasonable requests to the corresponding author.