From flexoelectricity to absolute deformation potentials: The case of SrTiO
Abstract
Based on recent developments in the firstprinciples theory of flexoelectricity, we generalize the concept of absolute deformation potential to arbitrary nonpiezoelectric insulators and deformation fields. To demonstrate our formalism, we calculate the response of the band edges of SrTiO to both dynamic (sound waves) and static (bending) mechanical loads, respectively at the bulk level and in a slab geometry. Our results have important implications for the understanding of straingradientrelated phenomena in crystalline insulators, formally unifying the description of bandstructure and electrostatic effects.
pacs:
71.15.m, 77.65.j, 63.20.dkI Introduction
The socalled deformation potentials, introduced in 1950 by Bardeen and Shockley, Bardeen and Shockley (1950) describe the shifts in the singleparticle electronic energy levels that are induced by a strain field. Originally aimed at estimating the intrinsic contributions to carrier mobility in semiconductors, this concept is nowadays important for a wide range of phenomena of fundamental and technological interest. For example, deformation potentials were invoked to explain band offset formation at strained semiconductor heterojunctions Peressi et al. (1998), carrier confinement in optoelectronic devices, Janotti and Van de Walle (2007) and the impact of dopants and impurities on the lattice parameter of insulators Bruneval et al. (2015); Janotti et al. (2012). In all the above contexts, one needs to know how the energy of a given electronic state is modulated by an inhomogeneous strain field, which can be produced by an acoustic phonon, by a compositional gradient or by an applied external load. This physical property goes now by the name of absolute deformation potential (ADP), where the qualification “absolute” was introduced Van de Walle and Martin (1989) to emphasize the focus on the single energy levels, rather than their relative differences.
For theoretical analysis and computational purposes, ADPs are typically split into two separate contributions, a macroscopic electrostatic (ME) term and a bandstructure (BS) term. The latter consists in the strain derivative of the Bloch eigenvalues, where the (arbitrary) energy reference is fixed to the average electrostatic potential (or to the Fermi level in a metal). This is manifestly a bulk property, as it can be readily obtained by performing periodic calculations that involve the primitive unit cell of the crystal only. The former, which is by far the most challenging to calculate, consists in the strain derivative of the electrostatic reference itself. (This term is due to longrange electrostatic interactions, and therefore it is only present in insulators.) As such reference is generally illdefined in an infinite crystal, Kleinman (1981) one needs to explicitly consider an inhomogeneous strain field, where the cell parameters gradually evolve between two different configurations A and B. The integral of the electric field along the AB path yields then the soughtafter voltage lineup. The nonanaliticity of electrostatic interactions implies that the lineup depends not only on the end points, A and B, but also on the specific path connecting them. Yet, once the direction of such a path is known (assuming it is rectilinear), this term can also be expressed as a bulk property, Resta et al. (1990); Resta (1991) and explicit firstprinciples calculations have been reported for a number of materials. Van de Walle and Martin (1989); Resta et al. (1990); Janotti and Van de Walle (2007)
In all the aforementioned works, the focus has been mainly restricted to a particular class of deformations (longitudinal acoustic phonons) and dielectrics (nonpolar semiconductors like Si or Ge). These essentially reflect the established range of validity of the current ADP theory. Resta et al. (1990) There are increasingly good reasons, however, to seek a more general approach to the problem in order to overcome the above limitations. Regarding the materials issue, there has been recently increasing interest in semiconducting oxides: ZnO, SrTiO, etc. with highly confined 2D electron gases. Acoustic phonon scattering has been identified as one of the main factors limiting mobility, e.g., at the LaAlO/SrTiO interface, either in the superconducting state Klimin et al. (2014), or in thermoelectric applications Pallecchi et al. (2015); in other cases, ADPs were invoked to explain the carrier confinement mechanism itself. Janotti and Van de Walle (2007) Polar materials (in this work we shall indicate as “polar” all crystals where the Born charge tensor does not vanish) presently fall outside the scopes of ADP theory, Resta et al. (1990) and this poses a clear problem for the interpretation of the above phenomena. Regarding the “geometric” issue, there has been recently a strong interest in the electronic properties of bent nanostructures, and in particular in determining how curvature affects the effective potential felt by quantum particles. Ortix et al. (2011) Accurate knowledge of the bendingrelated ADPs is necessary to obtain a quantitative solution to this problem, Ortix et al. (2011) but bending is a type of deformation that cannot be described as a longitudinal acoustic phonon; this seriously limits the applicability of firstprinciples methods in this context.
Here we show that the theory of flexoelectricity, which has undergone an impressive development in the past few years, is now mature enough to provide a rigorous answer to both questions. Flexoelectricity (FxE) describes the electrical polarization that is linearly induced by a strain gradient, and is therefore ideally suited to tackle the macroscopic electrostatic contribution to the ADPs in the most general case. (Indeed, the longitudinal components of the macroscopic FxE tensor reduce Resta (2010) to the ME term of Ref. Resta et al., 1990 in a nonpolar crystal.) In particular, we argue that (i) the atomic relaxations (internal strains) that occur in a straingradient field are bulk properties of the material (just like the purely electronic effects that were considered in earlier works Resta (1991); Resta et al. (1990)), thereby extending the ADP theory to polar insulators; and (ii) the transverse components of the FxE tensor naturally yield the desired information on the bendinginduced electrostatic effects, which allows us to formally extend the ADP theory to arbitrary deformation fields. As a practical demonstration, we provide a full calculation of the ADPs in cubic SrTiO, considering both the dynamic regime of sound waves and the static one of a bent slab.
Apart from the immediate relevance to the physics of SrTiO, these results highlight a crucial aspect of flexoelectricity that is absent in other types of electromechanical couplings: As each energy level experiences a different tilt, in a macroscopic strain gradient the very notion of macroscopic electric field is inherently ambiguous, and depends on the physical nature (electron, hole, classical “test” charge, etc.) of the charged particle that is used to define it. An elegant solution to this conceptual puzzle consists in identifying the flexovoltage coefficients with the ADPs, and thereby rationalizing the aforementioned ambiguity via the bandenergy dependence of the deformation potential. This points to a formal unification of ADP and flexoelectricity theories, with important implications for both the interpretation of the experiments and the theoretical modeling of straingradient related phenomena.
Ii Theory
Consider a macroscopic deformation where a given component of the symmetric stress tensor, , undergoes a linear increase along the direction . The absolute deformation potential of the th band at the point in the Brillouin zone can be written as
(1) 
Here is the energy eigenvalue referred to the mean electron potential energy,
(2) 
is the typeII flexoelectric tensor (referring to the derivative of the polarization with respect to the strain gradient component , where the latter is the gradient of along the Cartesian direction ), is the relative permittivity along , and is the vacuum permittivity. The first term on the rhs of Eq. (1) is the bandstructure term, and is independent of the direction ; following earlier works, Van de Walle and Martin (1989); Resta et al. (1990) we shall focus on the valence () and conduction () band edges, and indicate the corresponding BS terms as . The second term in Eq. (1) is minus the electric field that is generated by the strain gradient when opencircuit boundary conditions are imposed along ; this corresponds precisely to the macroscopic electrostatic term in the ADP. (Note the factor, reflecting the negative electronic charge.)
To see that the present theory recovers earlier treatments Resta et al. (1990) as special cases, consider the macroscopic displacement pattern associated with an acoustic phonon of wavevector ,
(3) 
where is a cell index, is a sublattice index, and is a realspace vector indicating the displacement amplitude and direction. For a longitudinal phonon, the symmetric strain is
(4) 
where and . We readily obtain, for the ME part,
(5) 
It is then straightforward to show Stengel (2013a) that, in the case of a nonpolar insulator, reduces to the expression of Ref. Resta et al., 1990, written in terms of the dynamical quadrupole and octupole tensors. In a material with cubic symmetry this can be written, in turn, as Resta et al. (1990)
(6) 
where the anisotropy function varies between 0 and 1 depending on the direction. Then, the final result is
(7) 
where .
Having shown the consistency of Eq. (1) with earlier works, we can now move on to clarifying in what respects such expression generalizes the scopes of the ADP theory to a broader range of phenomena. In a nutshell, by using Eq. (1) to define the ADPs, one can in principle tackle all materials (including polar crystals) and deformation fields (including bending) for which the flexoelectric tensor is well defined. To appreciate the practical implications of this statement, it is useful to discuss in some detail the example of a polar (but nonpiezoelectric Stengel (2013a)) insulator. In Ref. Resta et al., 1990, the authors argued that ADPs cannot be written as bulk material properties in this case, as the displacement of a plane of atoms induces a net dipole density, and hence an arbitrary shift in the electrostatic lineup. The theory of flexoelectricity, however, shows that such displacements are not arbitrary, but can be themselves written as bulk material properties, via a higherorder internalstrain tensor (indicated as or in Ref. Stengel, 2013a). This allows one to write, in full generality, the latticemediated (LM) contribution to the ADPs, which for a LA phonon in a cubic insulator has the same form as Eq. (6). (The frozenion and relaxedion ADPs retain the same functional form, only with different values of the coefficients and .)
It is worth making, in this context, an additional comment on the physical nature of the LM contribution. Such contribution depends, as we said, on the flexoelectric internalstrain tensor components, and this is an inherently dynamic quantity, i.e. it depends explicitly on the atomic masses; Stengel (2013a) this characteristic obviously propagates to the relaxedion ADPs. This is not a concern when studying, e.g., the impact of acoustic vibrations (a dynamic effect) on carrier mobility, but questions, at first sight, the applicability of the present theory to genuinely static cases, such as a bent slab under an external mechanical load. Note, however, that the condition of static equilibrium imposes a precise linear relationship between the straingradient components that are present in a given region of the sample interior. (For example, in a bent plate the “primary” transverse deformation is always accompanied by a “secondary” longitudinal strain gradient that is oriented along the surface normal. Stengel (2013b, 2014)) Once these elastic relaxation effects are appropriately incorporated, the massdependent terms that are present in the internalstrain tensor exactly cancel out. Stengel (2013a) This clears up the above worries regarding latticemediated effects: the ADPs described by Eq. (1) are indeed static quantities in the slabbending case, while they are dynamic in the context of LA phonons, consistent with the physical nature of the phenomenon under study.
Iii Results
iii.1 Longitudinal acoustic phonons
Frozenion  16.15  0.33  1.40  1.29 

Total  18.68  2.20  3.93  0.98 
In the following, we shall illustrate the above arguments by explicitly calculating the relaxedion ADPs in SrTiO. Our calculations are performed within the localdensity approximation Perdew and Wang (1992) to densityfunctional theory. We use the same pseudopotentials and computational parameters as in Ref. Stengel, 2014, from which we also take the numerical values of the relevant bulk flexoelectric tensor components. These data are then processed to yield the macroscopic electrostatic contributions to the ADPs. The bandstructure terms, are calculated by finite differences, considering an isotropic volume expansion (or compression) of the cubic SrTiO cell. This corresponds to neglecting the splittings induced by anisotropic deformations, consistent with earlier works. Van de Walle and Martin (1989); Resta et al. (1990) The reported values correspond to the valenceband top at the point of the Brillouin zone, and to the conductionband bottom at .
We shall start by discussing our results for the “acoustic” ADPs, which are described by Eq. (6). Our calculated numerical values for , , and , both at the frozenion and relaxedion levels, are reported in Table 1. The case of a [100]oriented LA phonon is also schematically illustrated in Fig. 1.
First, note the large difference in absolute value between and , which points to a substantial cancellation between the (positive) macroscopic and (negative) bandstructure terms. A closer look at the numerical data indicates that such cancellation occurs at the purely electronic (frozenion) level. To rationalize this outcome (which is a common feature of ADP calculations Resta et al. (1990)), recall that the frozenion value of the longitudinal flexoelectric coefficients (and hence of ) are related Hong and Vanderbilt (2013); Stengel (2013a) to the octupolar moments of the electronic charge that each atom “drags” along during the deformation. Such moments are typically dominated by a spherical term, which originates from the inner orbitals that rigidly follow each nucleus along its distortion path. These orbitals, in turn, are chemically inert, meaning that their displacement contribute to but not to the band energies. (In a ideally ionic solid, e.g. like the noninteracting rigidion model of Ref. Stengel, 2013b, the cancellation would be complete, and yield vanishing values for .)
Next, it is interesting to note the relatively large and positive contribution from latticemediated effects. (Once relaxation effects are accounted for, both bands undergo an upward shift upon uniaxial tension, see Fig. 1.) This observation can be readily explained by recalling the dynamical nature of the latticemediated flexoelectric effect in sound waves. Indeed, a strain gradient dynamically induces (in addition to the massindependent contributions from the interatomic force constants) a force on each ionic sublattice that is directly proportional to its nuclear mass Stengel (2013a). Since the cations (Sr, Ti) are much heavier than the anions (O), this produces a systematic bias on the sign of the flexoelectric coefficients, consistent with the results of Table 1. To verify such a hypothesis, we have performed a computational experiment, where we have recalculated the latticemediated contribution to , , while setting all nuclear masses to a uniform value. As a result, changed from 2.53 eV to 1.88 eV, confirming our arguments. (A closely related analysis was also performed in Ref. Hong and Vanderbilt, 2013, with qualitatively similar conclusions, see Table IV therein.)
Finally, note the relatively small anisotropy, which is only marginally affected by latticemediated effects, and the comparably larger value of respect to . The latter observation indicates that the gap of SrTiO shrinks upon volume expansion, consistent with the dominant ionic character of the bonding. (Our calculated gap deformation potential is eV.)
Our frozenion results can be directly compared with those that were recently reported by Janotti et al. Janotti et al. (2012) There appears to be a significant discrepancy ( eV, eV, eV in Ref. Janotti et al., 2012), whose precise origin is unclear at present. One possible source of disagreement may stem from the different exchange and correlation functional that was used therein (HSE). Heyd et al. (2003) This, however, would imply a worrisome dependence of the ADPs on the specific details of the computational model; further investigations will be necessary in order to clarify this important point.
iii.2 Static bending of a (100)oriented plate
Frozenion  10.37  1.31  0.08 

Total  10.81  0.87  0.36 
Our flexoelectricitybased theory of ADPs can readily answer another important physical question Ortix et al. (2011) that was so far unaddressed at the firstprinciples level. This concerns the modifications to the electronic structure of a slab that are induced by a bending deformation. The validity of the ADP concept in this specific case rests on the demonstration, provided in Ref. Stengel, 2013b, that the opencircuit internal field of a slab subjected to bending is a bulk material property. As we mentioned earlier, it is most appropriate here to assume mechanical equilibrium (under the action of an external load), and therefore consider static ADPs. These can be expressed as effective linear combination of the transverse and longitudinal flexovoltage coefficients, and decomposed according to Eq. (1) into a macroscopic and a bandstructure term,
(8) 
Here corresponds to minus the effective flexovoltage coefficient, Stengel (2014) , scaled by the electron charge. The bandstructure term, on the other hand, consists in the usual isotropic volume shift with a prefactor, Stengel (2014) , that accounts for Poisson’s ratio effects,
(9) 
(, where is the bulk elastic tensor in Voigt notation.)
As a testcase, we considered a (100)oriented slab of SrTiO in the platebending regime, analogous to Ref. Stengel, 2014. (A more extensive analysis, involving slabs of different orientations and a comparison of the beambending versus platebending limits, is reported in Section III.3.) In Table 2 we show our results for and . Again, we calculated both quantities at the frozenion level first (note that here “frozenion” implies neglect of the internal strains, but inclusion of the aforementioned Poisson’s ratio effects), and later incorporating full atomic relaxation. As in the phonon case discussed earlier, there is a large cancellation between the macroscopic and the bandstructure terms at the frozenion level. Due to this cancellation, the relative impact of the latticemediated contribution, small on the scale of the large macroscopic term, becomes significant in the context of the ADPs (the latter are about an order of magnitude smaller than ); for example, the valenceband parameter changes sign after relaxation. Compared to the phonon ADPs of Table 1, here the relaxedion values are significantly smaller in magnitude; also, the conductionband and valenceband ADP have opposite signs. This corroborates our interpretation that the large and positive values of the phonon ADPs are mainly due to dynamical effects, which are absent in the present slabbending case. The above results are schematically illustrated in Fig. 1, where we also have incorporated the results of Ref. Stengel, 2014 regarding the surface contributions to the opencircuit flexovoltage.
iii.3 Back to flexoelectricity
The energy level diagram of Fig. 1 may appear, at first sight, confusing: in the slab interior we have represented bandstructure effects (which are the realm of ADP theory), while in the vacuum regions we have plotted the electromechanical response coefficients of Ref. Stengel, 2014 (which were derived in the context of flexoelectricity). Such a juxtaposition was made on purpose – as we shall explain shortly, there is an even more intimate link between flexoelectricity and ADPs than the earlier sections of this work suggest.
To see this, consider the total opencircuit voltage induced by bending, , which in Fig. 2 corresponds to (minus) the net shift between the (constant) vacuum potentials located at either side of the bent slab. (This is the physical quantity that flexoelectric experiments typically focus on.) In Ref. Stengel, 2014, such a quantity was obtained by summing up the relevant bulk flexovoltage coefficient (corresponding to i.e. the tilt of the macroscopic electrostatic potential in the slab interior) and the surface dipolar contribution, ,
(10) 
Fig. 2 shows that the same quantity can be, equivalently, obtained by summing up the bulk ADP (either the valenceband or the conductionband ), and the appropriate surface deformation potential Zhou et al. (2013), ,
(11) 
(The latter is the strain derivative of the offset between the conduction or valence band edges and the vacuum level; such offsets are commonly known in the literature by the name of electron affinity and ionization potential, respectively.) In our slab models of SrTiO, we find eV, and eV, where the superscript refers to the two types of surface terminations considered in Ref. Stengel, 2014. (The corresponding values for the conduction band surface deformation potentials can be trivially determined by using the data of Table 2.)
The fact that, by summing up either the bulk and surface flexovoltage coefficients or the corresponding valenceband ADPs, we have obtained the same number (the total flexovoltage response of the slab) is not fortuitous. Flexovoltage coefficients really are ADPs in disguise (they can be thought as “electrostatic potential ADPs”) and, similarly, ADPs can be perfectly well regarded as flexovoltage coefficients. The former statement is an obvious consequence of Eq. (1): by setting the bandstructure term to zero, and by recalling the relationship between flexovoltage and flexoelectric coefficients Stengel (2014), one trivially obtains . Accepting the latter statement, on the other hand, may seem awkward at first sight, as it implies that the bulk flexoelectric tensor is not uniquely defined. (ADPs depend on the specific band feature that is being considered.) This is, however, perfectly consistent with the fundamental theory of flexoelectricity, where such a nonuniqueness emerges as a natural consequence of the formalism. Stengel (2013a)
To appreciate the physical origin of this arbitrariness, one does not necessarily need to go through pages of complicated algebra. In fact, it suffices to take a closer look at Fig. 2 and Table 2 to get a reasonably clear idea of what’s going on: The valence band, conduction band and mean electrostatic potential all undergo a different tilt. This means that a conduction electron, a valence hole and a classical “test” charge will not feel the same electrical force. Otherwise stated, in a uniform strain gradient the very definition of “macroscopic electric field”, , depends on the physical nature of the charged particle that is used to probe it. As the field in the interior of a bent slab is, in the above sense, ambiguous, the notion of “shortcircuit boundary conditions” (a necessary ingredient for defining the bulk flexoelectric tensor) is equally ambiguous therein – it depends on the (arbitrary) energy reference that we choose when imposing the (i.e. flatband) condition. ^{1}^{1}1Of course, one can always think of depositing metal electrodes onto the slab surfaces – this way, the shortcircuit condition is no longer arbitrary. Still, the precise way that a “zero external potential” translates into a “vanishing internal field” depends on the properties (i.e. on the interface deformation potential) of the filmelectrode interface. In some sense, the interface/surface “selects” a given energy reference out of the infinite possibilities that exist at the bulk level, thereby lifting the aforementioned arbitrariness.
Now, while it is true that selecting one or the other reference potential to define the internal electric field of the slab does not affect the overall result (the total opencircuit voltage), this choice does modify the way the effect is split into surface and bulk contributions. In practice, such a freedom can be exploited to achieve a more meaningful physical description of the two individual pieces. In the present case of the SrTiO slab, for example, it appears tempting to identify the relevant flexovoltage coefficients with the valenceband (VB) ADPs,
(12) 
rather than with the gradient of the bare electrostatic potential as it was done earlier Stengel (2014). This way, we obtain a partition between bulk and surface effects where both contributions are small in magnitude (i.e. of the same order as the total opencircuit flexovoltage of the slab), thereby facilitating the identification of the physical mechanisms that play a dominant role in either context. (Recall that the “electrostatic” flexovoltages are typically characterized by a large cancellation between bulk and surface contributions, Stengel (2014) which complicates the physical interpretation of the two individual terms.)
iii.4 Bending anisotropy of SrTiO slabs
Conduction band  

(100)  (110)  (111)  
Plate  0.87  1.97  1.78 
Beam  0.68  1.73  1.33 


Valence band  
(100)  (110)  (111)  
Plate  0.35  0.83  0.63 
Beam  0.27  0.85  0.47 
To illustrate the above ideas, it is an insightful excercise to compute the bulk ADPs under static slab bending by considering different surface orientations and mechanical boundary conditions. (The two extremes in the latter context are represented by the platebending and beambending limits.) This is especially interesting in light of the experimental results of Ref. Zubko et al., 2007, which concern the flexoelectric response (under bending) of SrTiO slabs with different orientations. Zubko et al. Zubko et al. (2007) argued, based on their analysis, that the measured coefficients can be interpreted reasonably well by assuming a purely bulk flexoelectric response, i.e. by neglecting possible surface effects. It would then be desirable to compare the reported values with the existing theoretical estimates of the bulk flexoelectric coefficients. Interestingly, reliable ab initio calculations of these quantities have been recently performed Stengel (2014); Hong and Vanderbilt (2013), but the reported values show a marked discrepancy with the experimental data: The firstprinciples results for the bulk flexovoltage coefficients are of the order of V, and are systematically negative, while the measurements cluster Zubko et al. (2013) around V, with a positive or negative sign depending on sample orientation. To account for such a discrepancy, the contribution of surfaces is typically invoked. Stengel (2014) The concepts developed in this work, however, suggest that an alternative interpretation is possible: The disagreement may be largely due to an unfortunate choice of the energy reference when calculating the bulk contribution, and only to a lesser extent to the aforementioned surface effects.
Following up on this speculation, we report in Table 3 a complete overview of the (bulk) conductionband and valenceband ADPs in a SrTiO slab subjected to static bending. (The three surface orientations and the choices of the bending axes are consistent with the experimental setup of Ref. Zubko et al., 2007; for completeness we report the results for both the platebending and beambending regimes.) Clearly, the calculated ADPs are much closer to the experimentally measured flexoelectric data of Zubko et al. than the bulk coefficients that were quoted in Refs. Stengel (2014); Hong and Vanderbilt (2013), corroborating our point. (The values of Table 3 are of the order of 1 V; both negative and positive values are present.) Note that the difference between the platebending and beambending results is minor, which suggests that assuming one or the other regime might be of relatively little importance for a qualitatively correct interpretation of the experiments. Note also the relatively small orientation dependence, which amounts to about 1 V; this supports the idea, recently proposed in Ref. Narvaez et al., 2015, that an unusually large anisotropy (as measured by Narvaez et al. in BaTiO samples) is a clear signature of other effects (i.e. not bulklike in origin) being at play.
Iv Conclusions and outlook
We have proposed a unified perspective on the theory of absolute deformation potentials and flexoelectricity. To illustrate our ideas, we have used cubic SrTiO (either in bulk or slab form) as a testcase. The concepts developed here have important implications, both for the interpretation of the experimental measurements and for the macroscopic modeling of electromechanical phenomena. In the latter context, our work highlights (and corroborates with quantitative examples) both the arbitrariness of the bulk flexoelectric tensor, and the need for a consistent treatment of bulk and surface effects in order to achieve a sound physical picture. On a positive note, our work opens the way to the firstprinciples (and firstprinciplesbased) study of an essentially unlimited range of phenomena related to inhomogeneous deformation fields, complementing longrange electrostatics with bandstructure effects.
Acknowledgments
This work was supported by MINECOSpain (Grant No. FIS201348668C22P), and Generalitat de Catalunya (2014 SGR301). We thankfully acknowledge the computer resources, technical expertise and assistance provided by the Supercomputing Center of Galicia (CESGA).
References
 Bardeen and Shockley (1950) J. Bardeen and W. Shockley, Phys. Rev. 80, 72 (1950).
 Peressi et al. (1998) M. Peressi, N. Binggeli, and A. Baldereschi, J. Phys. D: Appl. Phys. 31, 1273 (1998).
 Janotti and Van de Walle (2007) A. Janotti and C. G. Van de Walle, Phys. Rev. B 75, 121201(R) (2007).
 Bruneval et al. (2015) F. Bruneval, C. Varvenne, J.P. Crocombette, and E. Clouet, Phys. Rev. B 91, 024107 (2015).
 Janotti et al. (2012) A. Janotti, B. Jalan, S. Stemmer, and C. G. Van de Walle, Applied Physics Letters 100, 262104 (2012).
 Van de Walle and Martin (1989) C. G. Van de Walle and R. M. Martin, Phys. Rev. Lett. 62, 2028 (1989).
 Kleinman (1981) L. Kleinman, Phys. Rev. B 24, 7412 (1981).
 Resta et al. (1990) R. Resta, L. Colombo, and S. Baroni, Phys. Rev. B 41, 12358 (1990).
 Resta (1991) R. Resta, Phys. Rev. B 44, 11035 (1991).
 Klimin et al. (2014) S. N. Klimin, J. Tempere, J. T. Devreese, and D. van der Marel, Phys. Rev. B 89, 184514 (2014).
 Pallecchi et al. (2015) I. Pallecchi, F. Telesio, D. Li, A. Fête, S. Gariglio, J.M. Triscone, A. Filippetti, P. Delugas, V. Fiorentini, and D. Marré, Nature Communications 6, 6678 (2015).
 Ortix et al. (2011) C. Ortix, S. Kiravittaya, O. G. Schmidt, and J. van den Brink, Phys. Rev. B 84, 045438 (2011).
 Resta (2010) R. Resta, Phys. Rev. Lett. 105, 127601 (2010).
 Stengel (2013a) M. Stengel, Phys. Rev. B 88, 174106 (2013a).
 Stengel (2014) M. Stengel, Phys. Rev. B 90, 201112(R) (2014).
 Stengel (2013b) M. Stengel, Nature Communications 4, 2693 (2013b).
 Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
 Hong and Vanderbilt (2013) J. Hong and D. Vanderbilt, Phys. Rev. B 88, 174107 (2013).
 Heyd et al. (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118, 8207 (2003).
 Zhou et al. (2013) M. Zhou, Z. Liu, Z. Wang, Z. Bai, Y. Feng, M. G. Lagally, and F. Liu, Phys. Rev. Lett. 111, 246801 (2013).
 Zubko et al. (2007) P. Zubko, G. Catalan, A. Buckley, P. R. L. Welche, and J. F. Scott, Phys. Rev. Lett. 99, 167601 (2007).
 Zubko et al. (2013) P. Zubko, G. Catalan, and A. K. Tagantsev, Annu. Rev. Mater. Res. 43, 387 (2013).
 Narvaez et al. (2015) J. Narvaez, S. Saremi, J. Hong, M. Stengel, and G. Catalan, Phys. Rev. Lett. 115, 037601 (2015).