The role of surfaces in flexoelectricity

Flexoelectricity describes the electromechanical coupling of strain gradient and polarization.^1–3 First experimentally observed in solids in 1968,⁴ it has only become the focus of significant research over the last 20 years.^1–3 Efforts to understand the flexoelectric effect, largely driven by the presence of large strain gradients at the nanoscale^5–8 and the fact that the coupling is allowed by all space groups,^1–3 have made clear that flexoelectricity has significant implications for a diverse set of fields including biology,^9–11 energy harvesting,^12,13 and triboelectricity.¹⁴ However, utilizing flexoelectricity in these disciplines requires advancing the fundamental physics of flexoelectricity.

Significant progress toward a predictive, microscopic model for flexoelectricity has been made over the last decade with the advent of the first principles theory of flexoelectricity.^15–18 This framework has enabled the ab initio calculation of bulk flexoelectric coefficients, the parameters that describe the linear coupling between polarization and strain gradient in an infinite crystal. Unfortunately, the bulk flexoelectric coefficients predicted by this theory and experimental measurements often differ in both magnitude and sign.^2,15,19 While experiments have demonstrated that the flexoelectric response of a sample is very sensitive to defects such as doping,²⁰ dislocations,²⁰ and grain boundaries,²¹ this theoretical and experimental divide persists even in high-quality, single crystalline samples.^15,19 The prevailing interpretation attributed this deviation to surface contributions in the measured flexoelectric response.

Early ab initio calculations of flexoelectric coefficients supported this interpretation by noting that the flexoelectric response of a finite body possessed a surface dependency that did not tend to zero in the thermodynamic bulk limit.¹⁷ This observation had already been predicted by phenomenological theory^22,23 but was initially controversial.²⁴ Subsequent work has rigorously shown that the polarization of a finite crystal induced by inhomogeneous strain has surface-sensitive contributions, which arise from changes in what has been called a “surface dipole” with strain, do not go to zero in the bulk limit, and are comparable in magnitude to the bulk flexoelectric response.^25,26 It has also been shown that these contributions to the total flexoelectric response of a finite body can be expressed as the change in the mean-inner potential (MIP) with strain.^25,26 (A more rigorous analysis of this will be presented later in this paper, where we will also define what is meant here by the “surface dipole,” as this has a different usage in the surface science community.)

It is well-established that the MIP is sensitive to surface details;^27–29 this is not unexpected since from continuity of the potential, its value in the bulk will depend upon the atomic arrangements as the material transitions from bulk to vacuum. The MIP has been extensively studied due to its role in both electron diffraction and holography,^27,30,31 and there have been significant efforts to measure^32–34 and calculate MIPs^28,29,32 in a wide range of materials. Given the surface sensitivity of the MIP, it is reasonable that the strain derivative of the MIP (which contributes to the total flexoelectric response^25,26) should be similarly surface sensitive.

The surface sensitivity of flexoelectricity was first explicitly demonstrated by a set of density functional theory (DFT) calculations by Stengel²⁶ on bulk truncated SrTiO₃ (100) surfaces, in which it was shown that MIP contributions to flexoelectricity are of the same magnitude as bulk flexoelectric coefficients and lead to a sign reversal in the overall flexoelectric response of bent SrTiO₃ beams. While the work in Ref. 26 represents an important step toward understanding flexoelectricity by demonstrating the importance of surfaces, there has been no systematic consideration to date of the impact of surface structure, chemistry, polarity, and adsorbates on flexoelectricity. This is particularly important because most experimentally observed surfaces, including those in nominally simple systems such as SrTiO₃, are far from simple and deviate significantly from bulk truncations. Beyond over-simplified ideas such as avoidance of unbalanced charges or dipoles, in almost all cases for oxides the bonding is local and needs to be chemically correct: see Andersen et al.³⁵ and references therein for a recent analysis. Unfortunately, the bulk truncated surfaces considered thus far in flexoelectricity²⁶ probably never occur in real systems, as has been suspected in the oxide surface science community for many years and recently highlighted by Sokolović et al.³⁶ As appropriate, we will add further information about this later in the paper.

In this work, we explore the relationship between MIPs and flexoelectricity through DFT calculations on a range of experimentally observed, low energy surfaces on archetypal ionic (MgO), covalent (Si), and mixed ionic-covalent (SrTiO₃) crystals. We consider the impact of surface structure, chemistry, polarity, and adsorbates on the flexoelectric response of bent beams, demonstrating that even minor variations in surface structure, chemistry, and adsorbates lead to significant changes in the total flexoelectric response. The total MIP contribution to the flexoelectric response is found to be dominated by bending-induced volume changes and linearly scales with the unstrained MIP, and the surface-specific component is found to be proportional to the ionization potential. Additionally, we develop a means to estimate MIP contributions to the total flexoelectric response using atomic electron scattering factors based upon the Ibers approximation.³⁷ This approximation is found to estimate MIP contributions to flexoelectricity to within ∼20% of the DFT calculated values. We then use atomic electron scattering factors and bulk flexoelectric coefficients calculated in Ref. 15 to predict the total flexoelectric response for a number of cubic systems. Finally, we discuss strategies to experimentally disentangle bulk and surface contributions.

In the absence of an electric field, the linear coupling between polarization $(Pi)$ and strain gradient $(∂ϵkl/∂xj)$ is described by flexoelectric coefficients $(μijkl)$ defined by

These coefficients can be recast as flexocoupling voltages, which give the gradient of the average Coulomb potential arising from strain gradients, according to

where $χ$ is the dielectric susceptibility and $ε0$ is the permittivity of free space. Flexocoupling voltages are a more convenient description of flexoelectricity for our purposes and will be used throughout the remainder of this paper.

The flexocoupling voltages in Eq. (2) are bulk properties that can be computed from first principles.^15–18 Often, a linear combination of these tensor components known as the effective bulk flexocoupling voltage $fbulk$ is physically relevant for experiments.^15,19 The general effect of $fbulk$ in an infinite, centrosymmetric crystal under a constant strain gradient is shown in Fig. 1(a). Since the scope of this paper is to investigate the role of surfaces in flexoelectricity, we assume that the $fbulk$ values are known (or, at least, can be computed) and where appropriate use values from Ref. 15. Note that test calculations with the all-electron augmented plane wave + local orbitals WIEN2k code³⁸ yielded $fbulk$ values similar to those reported in Ref. 15.

In addition to $fbulk$⁠, the flexoelectric response of a finite body has surface contributions stemming from the change in the average Coulomb potential going from inside the crystal to vacuum as a function of strain. To understand the origins of this effect, consider a finite centrosymmetric crystal with a slab geometry as shown in Fig. 1(b). Such a slab has electric fields near the surface from the potential step but no potential difference across it because it is centrosymmetric. Homogeneously straining the slab modifies the magnitude of electric fields near the surface, but the net effect is still zero potential difference across the slab because homogeneous strain does not break inversion symmetry. Note that this potential step has been referred to as a surface dipole;^25,26 however, we will follow the surface science convention and reserve “surface dipole” for polar entities such as chemisorbed hydroxide.

As shown in Fig. 1(c), the application of a constant strain gradient not only induces a constant electric field from the bulk flexoelectric effect but also breaks the symmetry of the electric fields associated with the potential step between the bulk and vacuum. The potential difference is the sum of these two terms. Importantly, the surface contributions do not go to zero in the thermodynamic bulk limit and are comparable in magnitude to the bulk flexoelectric effect.^25,26 This indicates that, contrary to many other electromechanical couplings, the interpretation of measured flexoelectric responses requires consideration of the surface.

The surface contribution to the total flexocoupling voltage depicted in Fig. 1(c) can be expressed as the change in the mean-inner potential (MIP) with strain,^25,26

The MIP describes the difference between the average Coulomb potential in a solid $(Vavg)$ and the vacuum energy $(Evacuum)$ outside the solid.²⁷ It is only well-defined for finite-sized systems,³⁹ sensitive to surface orientation, structure, chemistry, and adsorbates,^28,29 and directly measurable with electron holography (see Refs. 31 and 33 and references therein). Motivation for Eq. (3) follows from the definitions of $fbulk$ and the MIP: the total potential difference will not only depend on the gradient of the average Coulomb potential $(fbulk)$ but also the difference between the average Coulomb potential and the vacuum energy (MIP). Therefore, if one neglects contributions from defects,^20,21,40 the total flexoelectric response of a finite crystal is

The term in Eq. (3) has been called the surface flexocoupling voltage in the literature,²⁶ but we refer to it as $fMIP$ because, as shown below, only a portion of $dV¯dϵ$ is truly surface sensitive: from a surface science viewpoint, only the local dynamic polarizability at the surface should be considered a “real” surface contribution (e.g., differences in Born charges near the surface⁴¹ or surface dipole contributions from chemisorbed species).

Given the importance of the MIP in flexoelectricity, we now carefully define the MIP as well as the true surface contributions before focusing on its implications for the flexoelectric effect. Figure 2 provides a band diagram description of the MIP. The surface sensitivity of the MIP in an insulator is readily observed by decomposing the MIP into the sum of (1) $Δ$⁠, the difference between $Vavg$ and the bulk valence band maximum $(EVBMbulk)$⁠, and (2) I, the difference between $Evacuum$ and $EVBMbulk$ (i.e., the ionization potential⁴²). The former is a bulk property, and the latter is a surface property. For the purposes of this work, it will be useful to further separate I into the sum of the work function (⁠$ϕ$⁠, the difference between $Evacuum$ and the surface VBM $EVBMsurf$⁠) and the surface valence band offset (⁠$δ$⁠, the difference between $EVBMbulk$ and $EVBMsurf$⁠).

Referring to the decomposition in Fig. 2 and the definition of $fMIP$ in Eq. (3), it is apparent that the MIP contribution to the total flexoelectric response has bulk and surface components,

The bulk component describes the change in the difference between the average Coulomb potential and the bulk VBM with strain. This term is a function of crystallographic orientation (as described in Sec. III), but only depends upon the bulk crystal structure for a given material. The surface component describes the change in the ionization potential with strain, which will be shown to be very sensitive to surface structure, chemistry, and adsorbates. Note that while the VBM is used to separate the MIP and $fMIP$ into bulk and surface terms in this work, this breakdown is not unique: for other materials, it may be more sensible to use another energy to distinguish between bulk and surface contributions, e.g., the conduction band minimum in n-type semiconductors.

Section III describes how $fMIP$ is calculated from first principles, and this quantity is analyzed in Secs. IV and V for several low energy surfaces. Such calculations are useful for understanding the role of surface chemistry, structure, and adsorbates in modifying a flexoelectric response but are often inapplicable to experiments where frequently little is known about the sample surface. In these situations, it is possible to analytically estimate the MIP and $fMIP$ from electron scattering factors. (For reference, electron scattering factors give the scattering amplitude of an electron by an isolated atom/ion and are proportional to the Fourier transform of the isolated atomic/ionic Coulomb potentials.) As shown by Ibers,³⁷ the MIP can be approximated as

where h is Planck's constant, m is the electron rest mass, e is the electron charge, $Ω$ is the unit cell volume, and $fiel(0)$ is the electron scattering factor of species i in the unit cell (atomic and ionic electron scattering factor values are tabulated for each element from Dirac–Fock calculations, e.g., Ref. 43). In the limit of small strains, from Eqs. (3) and (6), $fMIP$ can be approximated from electron scattering factors and volumetric strain according to

where $V¯0$ is the MIP at the equilibrium volume $Ω0$ and $ΔΩΩ0$ is the volumetric strain. The quality of these approximations is assessed in Secs. IV and V.

DFT calculations were performed with the all-electron augmented plane wave + local orbitals WIEN2k code.³⁸ Bulk calculations were first performed to find optimized lattice constants of SrTiO₃ (Pm-3m), MgO (Fm-3m), and Si (Fd-3m). Then conventional slab models using the bulk optimized lattice constants were constructed to simulate each surface of interest. These consisted of 15–20 atomic layers, a separation of 10–15 Å between surfaces, and utilized the highest available symmetry. In addition to performing typical checks on the converged surface slabs (e.g., bond lengths, bond valence sums, etc.), the density of states from the center of the slab was compared to the density of states calculated from a separate bulk calculation to ensure that the slabs were sufficiently large and contained enough vacuum. Additionally, test calculations on surface slabs that contained a larger number of unit cells in the bulk of the slab and slabs that contained a larger amount of vacuum yielded nearly identical MIP and $fMIP$ values (see below for how these are calculated), further confirming that the surface slabs used here were sufficiently large. Atomic positions in the surface slabs were simultaneously converged with the electron density using a quasi-Newton algorithm.⁴⁴ Unless otherwise stated, the local density approximation (LDA)⁴⁵ was used to approximate the exchange–correlation term. As described in Sec. IV, calculations performed with the PBEsol functional⁴⁶ and a PBEsol + on-site hybrid^47,48 approach (only for SrTiO₃, with an on-site hybrid fraction of 0.25 applied to Ti 3d states) yielded qualitatively similar results.

For SrTiO₃, muffin-tin radii of 2.31, 1.54, and 1.40 bohrs were used for Sr, Ti, and O, respectively, with a plane wave expansion parameter RKMAX of 5.98, energy cutoff of −6 Ry, and k-mesh equivalent to 4 × 4 × 4 per bulk conventional unit cell. For MgO, muffin-tin radii of 1.63, 1.20, and 0.60 bohrs were used for Mg, O, and H, respectively, with a RKMAX of 4.42 (2.21 with H), energy cutoff of −6 Ry, and k-mesh equivalent to 10 × 10 × 10 per bulk conventional unit cell. For Si, a muffin-tin radius of 2.00 bohrs was used with a RKMAX of 7, energy cutoff of −8 Ry, and k-mesh equivalent to 7 × 7 × 7 per bulk conventional unit cell. Numerical tests in which the muffin-tin radii, RKMAX, and k-mesh reported above were varied yielded consistent MIP and $fMIP$ values. All calculations used a Mermin functional at room temperature. For bulk calculations, convergence criteria of 10⁻⁴ e and 10⁻³ Ry were used. This yielded optimized lattice parameters of 7.290, 7.874, and 10.208 bohrs for SrTiO₃, MgO, and Si, respectively. For surface calculations, convergence criteria of 10⁻⁴ e, 10⁻³ Ry, and 10⁻³ mRy/bohr were used with a force tolerance of 0.1 mRy/bohr. Surface calculations were found to be particularly sensitive to force convergence and tolerance.

The MIP for each surface was calculated via a core-level approach^49,50 according to

This approach avoids macroscopic averaging that can have numerical problems; with infinite precision, they are identical. A bulk calculation was used to determine the difference between a deep core eigenvalue, $Ecore$⁠, and the average Coulomb potential, $Vavg$⁠. Then, a surface slab calculation was used to determine the difference between the Coulomb potential in the center of vacuum, $Evacuum$⁠, and the same deep core eigenvalue from the innermost atomic plane of the slab. Checks were made to ensure that the slab and vacuum were large enough to minimize oscillations in the average Coulomb potential in vacuum and recover the bulk structure in the innermost slab layer. Tests were also performed using different deep core eigenvalues for a particular structure to confirm consistency in $V¯$ calculations. The computational parameters and method used here yielded uncertainties in $V¯∼0.2V$⁠. Once $V¯$ was known, the energies of the bulk and surface VBM were used to decompose $V¯$ according to Fig. 2.

The $fMIP$ for each surface was calculated with a core-level approach analogous to the method used to calculate the MIP,

The quantity $(d(Ecore−Evacuum)dϵ)slab$ was determined by calculating the linear variation in the difference between $Evacuum$ and $Ecore$ from the innermost atomic plane of the slab with clamped plate bending strains of the form

In Eq. (10), the z axis is normal to the surface, $ν=c12c11$ for (100) cubic surfaces,⁵¹ and $ν=c11+2(c12−c44)c11+2(c12+2c44)$ for (111) cubic surfaces.⁵¹ Elastic constants for SrTiO₃, MgO, and Si were taken from Refs. 52–54, respectively. The quantity $(d(Vavg−Ecore)dϵ)bulk$ was determined from bulk calculations by calculating and scaling the hydrostatic deformation potential of $Ecore$ (i.e., the linear variation in $Vavg−Ecore$ with uniform strain⁵⁵) by the trace of Eq. (10). In all cases, at least five strains were used to determine slopes, and checks were performed with multiple core eigenvalues to ensure consistency. The quality of the linear fits is demonstrated in Fig. 3 with the Mg1s eigenvalue as $Ecore$⁠. Once $fMIP$ were calculated from Eq. (9), the energies of the surface and bulk VBM were used to decompose $fMIP$ into $fMIPbulk$ and $fMIPsurf$ as defined in Eq. (5). Overall, the parameters and methods outlined above yielded an uncertainty in $fMIP∼0.2V$⁠.

While the clamped plate bending in Eq. (10) was computationally advantageous [many of the surface slabs investigated here have tetragonal symmetry, which is preserved by Eq. (10)], beam bending is most common in experimental flexoelectric characterization.¹⁹ To facilitate comparison with experiment, all reported $fMIP$ values throughout the paper have been converted to pure beam bending values using the relation

As described in Sec. II, the MIP matters in the context of flexoelectricity because its strain derivative encapsulates the surface contributions to the flexoelectric response^25,26 in a manner that is consistent with the standard definition of bulk flexoelectric coefficients.^15,18 Given the intimate relationship between flexoelectricity and the MIP, we first examine the impact of surfaces on the MIP before investigating their impact on flexoelectricity. Though there have been some studies using DFT to investigate MIP variations arising from differences in surfaces (e.g., Ref. 28), such analysis is included here because the MIP has not been studied for many of the surfaces we consider. Table I includes the DFT calculated MIP for a wide range of known surfaces on SrTiO₃,^56–68 MgO,^69–71 and Si,⁷² which are archetypal mixed ionic-covalent, ionic, and covalent crystals, respectively. The DFT relaxed surface structures are depicted in Figs. 4–6; in most cases, they can also be found in the source references. For this work, we focus on surfaces that have been experimentally observed and/or lie near the theoretical convex hull.

First, we examine the MIP of different (100) SrTiO₃ surfaces to ascertain the impact of surface chemistry and structure on the MIP. There are two possible bulk truncations for a cubic (100) perovskite surface yielding (1 × 1) TiO₂ and SrO single-layer (SL) terminations for SrTiO₃. The (100) surface of SrTiO₃ is also known to possess many TiO₂ double-layer (DL) reconstructions (e.g., Ref. 35 and references therein). The (1 × 1), (2 × 2)A, (2 × 2)C, (2 × 1), and c(4 × 2) DL reconstructions studied here and shown in Fig. 4 possess identical surface stoichiometries with structural differences, thus providing a means to isolate the impact of surface structure on the MIP while keeping surface chemistry fixed.

Beginning with surface chemistry, Table I demonstrates that the MIP of (100) SrTiO₃ surfaces increases with the amount of excess TiO₂ on the surface from a value of 15.2 V for bulk truncated SrO (−0.5 TiO₂/1 × 1) to 17.7 V for bulk truncated TiO₂ (0.5 TiO₂/1 × 1) to an average of 18.4 V for TiO₂ DL reconstructions (1.5 TiO₂/1 × 1). This indicates that surface chemistry is an important factor affecting MIPs. The data in Table I also show that the MIP has a comparable sensitivity to surface structure. The MIPs of the stoichiometrically identical DL reconstructions have a 1.7 V range, which reflects the impact of differences in the coordination environment of the surface Ti–O polyhedra [e.g., the (2 × 2)A, (2 × 2)C, and c(4 × 2) DL reconstructions are similarly coordinated⁷⁴ and have nearly identical MIP]. The range in MIP for the DL reconstructions is comparable to the MIP difference between the chemically distinct SrO and TiO₂ SL bulk terminations, indicating that both surface chemistry and structure play a large role in determining the MIP of a surface.

To our knowledge, there has been one experimental MIP measurement for bulk terminations of (100) SrTiO₃⁷³ and no measurements for DL reconstructions, hampering our ability to make meaningful comparisons and highlighting an area for future work. The measurements on bulk terminations of (100) SrTiO₃ surfaces found MIPs of 13.3 and 14.6 V for SrO and TiO₂ terminations, respectively. In contrast, we calculated MIPs of 15.2 and 17.7 V for SrO and TiO₂ terminations, respectively. While our calculations capture the qualitative relationship between the MIP of the bulk terminations indicated by these measurements, the source of the discrepancy in the magnitude of the values is unclear. (Note that as mentioned earlier, it is unlikely that pure bulk-like SrO and TiO₂ terminated surfaces actually exist in reality.³⁶) As shown in Table II, additional calculations with the PBEsol functional and a PBEsol + hybrid approach yielded qualitatively similar MIP to those calculated with the LDA functional. The product of the MIP and the optimized unit cell volume is also found to be constant across all functionals indicating that the deviation from experiment is not simply a difference in lattice parameter. Furthermore, the work functions we have calculated for SrO and TiO₂ SL surfaces are in good agreement with other DFT calculated work functions.⁷⁵ This consistency coupled with the good agreement between our calculations and electron holography MIP measurements for MgO and Si surfaces (discussed below) suggests that (1) the functionals used here are inadequate for SrTiO₃ surfaces, (2) there is a difference between MIP measured with different techniques (electron holography³³ vs reflection high-energy electron diffraction⁷³), or (3) some contamination affected the MIP measurements reported in Ref. 73.

Next, we focus on the MIP of MgO (100) and (111) surfaces. MgO was selected for its well-studied “polar” surfaces: bulk truncations of the (111) surface of a rock salt structure are not valence-neutral, and so these surfaces must stabilize via reconstruction, metallization, or adsorption.⁶⁹ MIP calculations were first performed on bulk truncations of MgO (100) to serve as a baseline for comparison with experiment and different MgO (111) reconstructions. The calculated MIP of 13.4 V for the MgO (100) surface is close to the MIP of 13.01 V measured with electron holography.³³ Using the experimental and theoretical MIPs and lattice parameters, we find $V¯LDA⋅ΩLDA=0.969Vnm3$ and $V¯exp⋅Ωexp=0.972Vnm3$ indicating that, unlike SrTiO₃, the difference between the calculated and experimental MIP for MgO is fully accounted for by considering the difference between the DFT optimized and experimental lattice parameters.

Turning now to unhydroxylated MgO (111) surfaces, we first analyze the MIP of the (2 × 2) octapolar reconstructions.⁷⁶ While these reconstructions have not been unequivocally observed as an isolated phase, they represent canonical, low energy structures.⁶⁹ We find that surface chemistry effects on the MIP of polar surfaces are similar to those on non-polar surfaces: there is a MIP difference of 2.5 V between the Mg-terminated and O-terminated octapolar structures, which is comparable to the difference between the SrO and TiO₂ bulk terminations of SrTiO₃. This indicates that there is nothing “special” about polar surfaces with regard to the MIP other than the fact that they must resolve their polarity through reconstruction or other means to exist. We also note that the MIP of the MgO (100) bulk truncation is roughly the average of the MIP of the octapolar structures, which follows from the difference in electron affinity of Mg and O.

Though the (2 × 2) octapolar surfaces have not been definitively experimentally observed, the closely related (2 × 2)-α structures have been.^69,70 The three variants of the unhydroxylated (2 × 2)-α structure have identical surface chemistries with minor structural differences corresponding to different combinations of two out of three surface sites occupied by oxygen [the generic (2 × 2)-α is shown in Fig. 5 with the three surface sites indicated by an orange hexagon, a green triangle, and a purple triangle]. Given the structural similarities, one would expect the three variants of the (2 × 2)-α structures to have similar MIP. While the (2 × 2)-α-O1 and (2 × 2)-α-O2 structures have nearly identical MIPs with values of 18.5 and 18.8 V, respectively, the (2 × 2)-α-O3 structure has a MIP of 17.1 V. This is likely owing to the relative stability of these three structures: the (2 × 2)-α-O3 is ∼0.5 eV per (1 × 1) unit cell lower in energy than the other two structures.⁶⁹ These results further support the conclusion that structural differences have a comparable impact to chemical differences on the MIP.

Next, the role of adsorbates was investigated by examining the MIP of hydroxylated MgO (111) surface reconstructions.⁶⁹ On average, the introduction of hydrogen reduces the MIP from 16.5 to 14.1 V (averaged over all unhydroxylated and hydroxylated structures) with the (1 × 1) hydroxylated structure exhibiting the lowest MIP (11.7 V) of all the investigated surfaces. The three variants of the hydroxylated (2 × 2)-α structures have smaller MIPs than their unhydroxylated counterparts, but their MIPs are reduced by different amounts, which further emphasizes the importance of surface structure: the MIP of (2 × 2)-α-OH1, (2 × 2)-α-OH2, and (2 × 2)-α-OH3 are decreased by 3.4, 3.9, and 1.4 V relative to (2 × 2)-α-O1, (2 × 2)-α-O2, and (2 × 2)-α-O3, respectively. The similarity between (2 × 2)-α-OH1 and (2 × 2)-α-OH2 reflects their relative energetic stability.⁶⁹ Finally, the Mg and O variants of the hydroxylated Rt3 structure⁶⁹ exhibit similar behavior to the (2 × 2)-octapolar structures, where the O-rich variant has a higher MIP. This result follows from the difference in the electron affinity of Mg and O.

To our knowledge, the MIP has not been measured for any (111) MgO surface, precluding a comparison with our calculations. Before moving onto Si, we will note that as with the SrTiO₃ structures, changing the DFT functional changed the MIP and optimized lattice parameter such that $V¯⋅Ω$ was constant for a given structure, independent of the functional used.

Finally, the (100) surface of Si was studied as an archetypal covalent material. The two surfaces studied in this material system were the (1 × 1) bulk truncation and asymmetric dimerized (2 × 1) reconstruction.⁷² Both are depicted in Fig. 6. These surfaces have identical chemistry and minor structural differences. Our calculations indicate a difference of 0.5 V between their MIPs which is smaller but comparable to the differences in the MIP of the chemically identical but structurally distinct SrTiO₃ and MgO surfaces investigated above. Our calculated MIP values are similar to some measured MIPs for Si;^28,32,33 however, we cannot make a quantitative comparison because we could not find MIP measurements for (100) Si surfaces, and there is a sizable spread in the available experimental values of the MIP for other Si surfaces.^28,32,33

Our calculations demonstrate that there is sizable variation in the MIP of nominally similar structures because the MIP is sensitive to surface structure, chemistry, and adsorbates. While calculations such as the ones presented here are useful for determining a MIP, they require the atomic structure of a surface to be known. This is frequently not the case in experiment. Therefore, we now assess the quality of predicting the MIP using electron scattering factors and experimental lattice parameters^77–79 according to the Ibers approximation given in Eq. (6). This method of determining the MIP is often applicable to experiments as it only requires knowledge of the bulk structure.

Using the appropriate atomic (ionic) electron scattering factors values from Ref. 43 in conjunction with Eq. (6) leads to predicted MIPs of 18.4 V (12.6 V), 15.1 V (22.3 V), and 13.8 V for MgO, SrTiO₃, and Si, respectively. Figure 7 compares the MIP predicted by electron scattering factors to those calculated with DFT for the surfaces analyzed above. We find both electron scattering factors do a comparable job in predicting the MIP, with mean absolute errors (MAEs) of 3.0 V for atomic electron scattering factors and 3.3 V for ionic electron scattering factors across all surfaces. Since MIP values are ∼15 V, this MAE indicates that the Ibers approximation is good to ∼20%. For individual materials, there are slight differences between the MIP predicted with ionic and atomic electron scattering factors. The MIP of MgO is better predicted by ionic electron scattering factors (MAE = 2.7 V) than atomic scattering factors (MAE = 3.5 V), whereas the MIP of SrTiO₃ is better predicted by atomic scattering factors (MAE = 3.0 V) than ionic scattering factors (MAE = 3.3 V). The predicted MIP of Si is comparable to the other two materials with a MAE of 3.0 V using atomic scattering factors (there are no ionic scattering factors for Si).

The results from Sec. IV demonstrate that the MIP is sensitive to the details of the surface including its structure, chemistry, and adsorbates. Using the same SrTiO₃, MgO, and Si surfaces shown in Figs. 4–6, we now explore the sensitivity of $fMIP$⁠, i.e., the MIP contribution to the flexoelectric response of a finite body, to these factors. For each of the investigated surfaces, Table III includes the DFT calculated values of the MIP contributions to the total flexoelectric response $(fMIP)$⁠, and the bulk $(fMIPbulk)$ and surface $(fMIPsurf)$ components of $fMIP$ defined by Eq. (5). Table III also lists the total flexoelectric response $(ftotal=fMIP+fbulk)$ of a sample with each surface with values for $fbulk$ taken from Ref. 15; $ftotal$ is analyzed in Sec. VI. Note that all values in Table III correspond to the beam bending strains described in Sec. III.

Beginning with the (1 × 1) SrO and TiO₂ SL structures on (100) SrTiO₃, we find that the surface-dependent contributions modify $fMIP$ by upward of 30% and $fMIPsurf$ values, which differ in sign and by 3.7 V. These results directly demonstrate that $fMIPsurf$ is sizable and reflect the importance of surface chemistry in determining the MIP contributions to the total flexoelectric response. Our calculations also corroborate the work by Stengel,²⁶ which found $fMIP$ of 6.6 and 9.4 V for SrO and TiO₂ SL structures, respectively.

To isolate the effects of surface structure on $fMIP$ in SrTiO₃ from those of surface chemistry, we analyze the (1 × 1), (2 × 2)A, (2 × 2)C, (2 × 1), and c(4 × 2) DL reconstructions shown in Fig. 4 because they have identical stoichiometry. Our calculations indicate that DL reconstructions have an average $fMIP$ of 9.0 V with a range of 1.6 V. This range is smaller, but comparable to the difference in $fMIP$ between SrO and TiO₂ terminations, which suggests that, like the MIP, $fMIP$ is sensitive to the details of surface structure. However, unlike what was found in the MIP calculations, there does not appear to be a distinct relationship between $fMIP$ and Ti–O polyhedra coordination or excess TiO₂ coverage. Instead, the surface contributions to $fMIP$ for all TiO₂-rich terminations are qualitatively similar with $fMIPsurf≥0V$⁠, which is in stark contrast to $fMIPsurf<0V$ for the SrO termination. This change in the sign and magnitude of $fMIPsurf$ for the TiO₂-rich and SrO terminations is consistent with the relative acidity and basicity of the surfaces, a trend which is also found in the MgO surfaces analyzed below.

Taken together, these results indicate that while surface structure is important in dictating the magnitude of $fMIPsurf$⁠, it is a secondary consideration to surface chemistry that controls the sign. This is most clearly exemplified by studying $fMIP$ in the (1 × 1) and (2 × 1) surfaces of Si because these two surfaces have identical surface chemistry and very minor structural differences compared to the structural differences in the SrTiO₃ DL reconstructions. As shown in Table III, the (1 × 1) and (2 × 1) have identical $fMIP$ within the uncertainty of these calculations.

To further assess the importance of surface chemistry and structure, and investigate the role of polar surfaces and adsorbates, we turn to MgO. The (100) and (111) surfaces of MgO studied here also demonstrate the existence of substantial surface-dependent modifications to $fMIP$⁠. As with the SrTiO₃ surfaces, we find the most important factor in determining $fMIPsurf$ to be surface chemistry, which dictates the sign. For example, there is a 2.5 V difference in $fMIPsurf$ of the Mg-terminated and O-terminated (2 × 2) octapolar structures as well as a change in sign, which is similar to the change in magnitude and sign of $fMIPsurf$ between the (1 × 1) bulk SL terminations of (100) SrTiO₃. In general, we find the Mg-rich MgO surfaces to have $fMIPsurf≤0V$⁠, O-rich MgO surfaces to have $fMIPsurf≥0V$⁠, and mixed termination surfaces (e.g., the (100) bulk-terminated MgO surface) to have $fMIPsurf≈0V$⁠. As with the SrTiO₃ surfaces, $fMIPsurf$ values for MgO surfaces track with the relative acidity and basicity of the surfaces.

Surface structure is also found to be an important consideration in MgO, particularly for surfaces with the same surface chemistry. For example, $fMIPsurf$ of the (2 × 2)-α-O3 structure is approximately double that of the (2 × 2)-α-O1 and (2 × 2)-α-O2 structures, even though these three structures have identical stoichiometries and only minor variations in surface site occupancy. The large differences in $fMIPsurf$ of the hydroxylated (2 × 2)-α structures similarly highlight the importance of surface structure. Note that as with the MIP, the similarities between the $fMIPsurf$ of the (2 × 2)-α-O1/O2 and (2 × 2)-α-OH1/OH2 structures is anticipated based on their similar energetic stability compared to the (2 × 2)-α-O3 and (2 × 2)-α-OH3 structures.⁶⁹ 

Regarding polar surfaces, our calculations indicate that the magnitudes of $fMIP$ and $fMIPsurf$ are similar for both polar and non-polar surfaces. This suggests that while polar surfaces may exhibit a wider range of $fMIP$ and $fMIPsurf$ values owing to their need to reconstruct, metallicize, or adsorb species; the flexoelectric response of polar surfaces is not intrinsically different than non-polar surfaces. On the contrary, adsorbates are found to have a rather dramatic effect on $fMIP$⁠, with hydroxylation tending to suppress the strain derivative of the MIP for MgO. On average, hydroxylation reduces $fMIP$ from 10.7 to 9.4 V (averaged over all unhydroxylated and hydroxylated structures). This behavior is most clearly demonstrated by comparing the $fMIP$ of the hydroxylated and unhydroxylated (2 × 2)-α structures, where the magnitude of $fMIP$ is reduced by 1.4, 0.9, and 5.0 V upon hydroxylation of (2 × 2)-α-O1, (2 × 2)-α-O2, and (2 × 2)-α-O3, respectively.

In general, though there are sizable variations in $fMIP$ across materials systems, even for nominally similar surfaces, we find that a good predictor (r = 0.87) for $fMIP$ is $V¯0ΔΩΩ0$⁠, where $V¯0$ is the DFT value for the MIP of the equilibrium structure and $ΔΩΩ0$ is the volumetric strain arising from bending. This dependency follows from the expression for $fMIP$ given in Eq. (7). Since the slope of $fMIP$ and $V¯0ΔΩΩ0$ across all investigated surfaces and material systems shown in Fig. 8(a) is very close to 1 (0.96), the volume change associated with bending appears to be the dominant contribution to strain-induced MIP changes.

Similarly, we find that variations in the ionization potential (and to a lesser extent, the work function) serve as a good predictor for $fMIPsurf$⁠. As shown in Fig. 8(b), there is a strong correlation (r = 0.79) between ionization potential and $fMIPsurf$⁠. This result indicates that materials with either high or low ionization potentials have strong surface-dependent contributions to flexoelectricity. Specifically, materials with high ionization potentials tend to have large, positive $fMIPsurf$ values, and those with low ionization potentials have large, negative $fMIPsurf$ values, with $I≈5.7V$ corresponding to the crossover from negative to positive $fMIPsurf$ values. Since $fMIPbulk≥0$ for beam bending, the fact that $fMIPsurf$ can be positive or negative means that $fMIPsurf$ (and in turn the ionization potential) plays a deciding role in determining the sign of the overall flexoelectric response. This suggests the prediction that the flexoelectric response of a (100) SrTiO₃ beam with an SrO termination will have a different sign than the flexoelectric response of a (100) SrTiO₃ beam with a TiO₂ termination²⁶ should be replicable in other systems.

As was noted at the end of Secs. II and IV, it is often experimentally the case that the detailed nature of a sample's surface is unknown making it impossible to disentangle bulk and surface contributions to a measured flexoelectric response using DFT calculations. In these cases, it is possible to approximate the MIP using electron scattering factors according to the Ibers approximation (which only relies upon knowledge of the bulk structure) and then estimate $fMIP$ using $fMIP≈V¯IbersΔΩΩ0$ (where $ΔΩΩ0$ is calculated from experimentally measured bending strains and elastic constants). In Fig. 9, we assess the quality of this approximation by comparing $fMIP$ calculated using DFT for the SrTiO₃, MgO, and Si surfaces shown in Figs. 4–6 with $fMIP$ calculated using atomic and ionic electron scattering factors. The MAE across all surfaces is 1.7 and 3.3 V using atomic and ionic electron scattering factors, respectively. The atomic electron scattering factors also outperform the ionic scattering factors for each individual material system with MAE of 2.0 (3.4) V for MgO, 1.3 (2.9) V for SrTiO₃, and 1.0 V for Si using atomic (ionic) electron scattering factors. This indicates that approximating $fMIP$ with atomic electron scattering factors is more accurate than with ionic electron scattering factors, even though both scattering factors performed similarly in approximating the MIP. Since the value of $fMIP$ across all surfaces was ∼10 V, utilizing the Ibers approximation with atomic electron scattering factors to estimate MIP contributions to flexoelectricity is good to ∼20%.

It is necessary to also include the effects of the bulk flexocoupling voltage, $fbulk$⁠, to obtain the total flexoelectric response one would experimentally measure for a bent beam with the surfaces investigated here. Combining the bulk flexocoupling voltage tensor components calculated from first principles by Hong and Vanderbilt¹⁵ with the expressions for the effective flexocoupling voltage of cubic beams with (100) and (111) surfaces described in Ref. 51 yields $fbulk$ values of −7.7 V for (100) SrTiO₃, −2.4 V for (100) MgO, −4.7 V for (111) MgO, and 0.8 V for (100) Si. The sum of $fbulk$ and $fMIP$ for each of the surfaces studied in this work is provided in Table III under the column $ftotal$⁠. Note that test calculations of bulk flexocoupling voltage tensor components with WIEN2k yielded similar values to those reported in Ref. 15. It is clear from these calculations that there is a large variation in the magnitude of the total flexoelectric response owing to surface-sensitive MIP contributions. For example, the total flexocoupling voltage of the MgO (111) surfaces investigated here has a 7.9 V range.

Beyond predicting values for samples with specific surfaces, known bulk flexocoupling voltage tensor components can be used in conjunction with the Ibers approximation to $fMIP$ to predict total flexoelectric responses. As an example, Table IV includes the predicted total flexocoupling voltages for all the cubic materials whose bulk flexocoupling voltage tensor components were calculated in Ref. 15. Beam geometries with {100}-type faces are assumed in these calculations. We find good agreement with the experimental SrTiO₃ flexocoupling voltages of 2.2 V¹⁹ and 4.5 V,⁵ but poor agreement with the experimental BaTiO₃ flexocoupling voltage of 22 V.⁵¹ Note that the disagreement with BaTiO₃ likely stems from the presence of precursor ferroelectric domains in experiments^51,80,81 and does not necessarily invalidate our application of the Ibers approximation in this context. Unfortunately, a more complete comparison to experimental measurements is not possible because the overlap between calculated and measured flexoelectric properties is rather small.

Among the many subtleties associated with the flexoelectric effect, perhaps the most counter-intuitive is the seemingly disproportionate impact of surfaces. As we have shown in this paper, the surface sensitivity of flexoelectricity follows naturally from recasting the problem in terms of the MIP and its strain derivative. For decades, it has been known that the MIP has a bulk component common to all surfaces of a given material and a sizable surface-specific component with a magnitude and sign determined by the detailed nature of the surface. Through our calculations on a range of low energy surfaces from archetypal ionic, covalent, and mixed ionic-covalent systems, we demonstrate that the total flexoelectric response of a finite sample exhibits a similar surface sensitivity owing to the impact of surface chemistry, structure, and adsorbates on the strain derivative of the MIP. Additionally, large variations in the signs and magnitudes of $fMIPsurf$ for both polar and non-polar surfaces lead to sizable differences in the total flexoelectric response of nominally similar surfaces. These results could explain the variation in the sign and magnitude of experimentally measured flexoelectric coefficients, though there may be additional interfacial effects from electrodes and/or adsorbed hydrocarbons in flexoelectric experiments. Moreover, since MIPs are affected by doping and defect states (e.g., Ref. 82), the framework introduced in this paper should be amenable to treating these effects in flexoelectricity.

Short of independently measuring $fMIP$⁠, our work also provides a series of approximations one may make to interpret flexoelectric measurements and separate $fbulk$ and $fMIP$ contributions in measured flexoelectric responses. The first level of approximation is to use the Ibers approximation with atomic electron scattering factors to estimate $fMIP$ according to Eqs. (6) and (7). While this is a simple approximation, it performs relatively well (e.g., Fig. 9) and is often the most appropriate since frequently too little is known about the surface of a sample to calculate $fMIP$ with DFT. The next level of approximation is to measure the MIP for a sample and then use the volumetric strain approximation given in Eq. (7) to estimate $fMIP$⁠. Note that some care should be taken if the MIP is measured with electron microscopy as phase contrast techniques are sensitive to electric field instead of polarization. Combining this value for $fMIP$ with the total measured flexoelectric response should allow one to isolate bulk and surface flexoelectric contributions and make direct comparisons to DFT calculations of $fbulk$⁠.

This work demonstrates that surfaces matter in flexoelectricity and emphasizes the applicability of surface science beyond traditional research domains such as catalysis. Our DFT calculations establish that surface chemistry, structure, and adsorbates all play a significant role in determining the MIP and its strain derivative and provide predictions of flexoelectric responses for a wide range of materials systems and surfaces. Furthermore, we establish approaches to experimentally disentangle surface and bulk contributions to flexoelectricity. Together, this suggests the enticing possibility of tailoring flexoelectric responses through procedures such as surface treatments but makes clear the need for flexoelectric experiments on single crystals with well-defined surfaces.

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-FG02-01ER45945, and the National Science Foundation under Award No. DMR-1507101. C.A.M. performed the DFT calculations and analysis under the supervision of L.D.M. Both authors contributed to the writing of the paper.

The data that support the findings of this study are available from the corresponding author upon reasonable request.