Flexoelectricity from density-functional perturbation theory

Flexoelectricity from density-functional perturbation theory

Massimiliano Stengel ICREA - Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona, Spain Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain
July 12, 2019

We derive the complete flexoelectric tensor, including electronic and lattice-mediated effects, of an arbitrary insulator in terms of the microscopic linear response of the crystal to atomic displacements. The basic ingredient, which can be readily calculated from first principles in the framework of density-functional perturbation theory, is the quantum-mechanical probability current response to a long-wavelength acoustic phonon. Its second-order Taylor expansion in the wavevector around the () point in the Brillouin zone naturally yields the flexoelectric tensor. At order one in we recover Martin’s theory of piezoelectricity [R. M. Martin, Phys. Rev. B 5, 1607 (1972)], thus providing an alternative derivation thereof. To put our derivations on firm theoretical grounds, we perform a thorough analysis of the nonanalytic behavior of the dynamical matrix and other response functions in a vicinity of . Based on this analysis, we find that there is an ambiguity in the specification of the “zero macroscopic field” condition in the flexoelectric case; such arbitrariness can be related to an analytic band-structure term, in close analogy to the theory of deformation potentials. As a byproduct, we derive a rigorous generalization of the Cochran-Cowley formula [W. Cochran and R. A. Cowley, J. Phys. Chem. Solids , 447 (1962)] to higher orders in . This can be of great utility in building reliable atomistic models of electromechanical phenomena, as well as for improving the accuracy of the calculation of phonon dispersion curves. Finally, we discuss the physical interpretation of the various contributions to the flexoelectric response, either in the static or dynamic regime, and we relate our findings to earlier theoretical works on the subject.

71.15.-m, 77.65.-j, 63.20.dk

I Introduction

Flexoelectricity, the electric polarization linearly induced by an inhomogeneous deformation, Kogan (1964) has become a popular topic in material science during the past few years Cross (2006); Zubko et al. (2007); Lee et al. (2011); Catalan et al. (2011); Lu et al. (2012); Nguyen et al. (2013); Zubko et al. (2013). The interest is motivated by the universality of the flexoelectric effect, which, unlike piezoelectricity, is present in insulators of any symmetry and composition (including simple solids such as crystalline Si and NaCl). 111We exclude from this discussion phenomena that occur in biological systems and liquid crystals, of slightly different nature. While flexoelectricity is not a new discovery (the effect was predicted by Kogan in 1964 Kogan (1964); the bending of a parallel-plate capacitor induced by an applied voltage was experimentally demonstrated in 1968 by Bursian et al. Bursian and Zaikovskii (1968)), it has traditionally been regarded as a very weak effect, hardly detectable in macroscopic samples. Only in the past few years research has taken off on this front, thanks to the breathtaking progress in the design and control of nanoscale structures. From an application point of view, reducing the size of the active elements is crucial to obtaining a sufficiently large response: The uniform strain gradient that can be sustained by a sample before material failure is inversely proportional to its lateral dimensions. Ironically, in the context of perovskite thin films, strain gradients (e.g. occurring during epitaxial growth) have long been regarded as harmful to the operation of ferroelectric memories, Catalan et al. (2005); Dawber et al. (2005) and only later explored as a potentially useful functional property. Several recent experimental breakthroughs Zubko et al. (2007); Catalan et al. (2011); Lu et al. (2012) have convincingly demonstrated that the effect can indeed be giant Lee et al. (2011) in thin films, large enough to rotate Catalan et al. (2011) and/or switch Lu et al. (2012) ferroelectric domains, or to replace conventional piezoelectric materials Fousek et al. (1999); Zhu et al. (2006) in sensors and transducers.

At the level of the theory, advances have been comparatively slow. For a long time, the main reference in the field was the seminal work by Tagantsev Tagantsev (1986), which focused on lattice-mediated responses only, and from a phenomenological perspective. Maranganti and Sharma Maranganti and Sharma (2009) have later applied the method of Ref. Tagantsev, 1986 to the calculation of the flexoelectric coefficients in selected materials. Unfortunately, a considerable spread emerged between the predictions of different microscopic models, hence the need for a more fundamental treatment. It has taken many years before a full first-principles calculation of the flexoelectric coefficients was attempted Hong et al. (2010). More recently Resta Resta (2010), and Hong and Vanderbilt Hong and Vanderbilt (2011) have established the basis for a general formulation of the problem in the context of electronic-structure density functional theory, but a unified approach, encompassing both electronic and lattice-mediated effects, has not emerged yet. Note that most theoretical treatments to date have defined the flexoelectric tensors starting from the real-space moments of localized response functions (either atomic forces induced on neighboring atoms Tagantsev (1986) or multipolar expansions of the charge response to atomic displacements Resta (2010)). This is a drawback in the context of electronic-structure calculations, where working with periodic functions would be preferable, as it would eliminate the need for expensive supercell geometries.

Given the incomplete state of the theory, there are pressing questions coming from the experiments that are still unresolved to date. First, whether the flexoelectric tensor is a well-defined bulk property has been a matter of debate for several years; Resta (2010); Tagantsev (1986) consequently, it is currently unclear if it is at all possible to separate the surface and bulk contributions in a typical experiment. Next, it has been pointed out Zubko et al. (2007) that static measurements alone leave the flexoelectric tensor undetermined – in order to solve for all the independent components one needs to combine static with dynamic data. Is it, however, physically justified to “mix” the two? What do we get as a result, a static or a dynamic quantity? Finally, of particular importance in the area of perovskite oxides (which are by far the best studied and most promising materials for flexoelectric applications) is the interplay of inhomogeneous deformations with the main order parameters (either ferroelectric polarization or antiferrodistortive tilts of the oxygen octahedral network). This has been the subject of several studies in the context of phenomenological Morozovska et al. (2012a, b) and effective Hamiltonian Ponomareva et al. (2012) approaches, but a systematic way to calculate the coupling coefficients (to be used as an input to the higher-level simulations) is still missing.

In a broader context, it worth noting that the interest of flexoelectricity is by no means limited to perovskites: For example, curvature-induced effects are of outstanding relevance in the physics of two-dimensional nanostructures Ortix et al. (2011) such as -bonded crystals (e.g. graphene Kalinin and Meunier (2008) or boron nitride Naumov et al. (2009)). Also, the electrostatic potential induced by deformation fields is a known concern in the performance of optoelectronic quantum-well devices, Lew Yan Voon and Willatzen (2011) especially in the promising area of foldable inorganic light-emitting diodes. Park et al. (2010) The theory of absolute deformation potentials, Van de Walle (1989) intimately related to flexoelectricity, Resta et al. (1990); Resta (2010) is an invaluable tool in the band-gap engineering of these (and other) semiconductor-based systems. Rationalizing these diverse and technologically important phenomena into a unified theory would be, of course, highly desirable from a modeling perspective.

Here we show how to consistently address the above issues by using density functional perturbation theory Baroni et al. (2001) (DFPT) as a methodological framework. By taking the long-wavelength limit of acoustic phonons, we derive the electromechanical tensors (both piezoelectric and flexoelectric) in terms of standard lattice-periodic response functions, which can be readily calculated by means of publicly available first-principles codes. We demonstrate the consistency of our formalism by rederiving already established results, such as Martin’s theory Martin (1972) of piezoelectricity and existing theories Tagantsev (1986); Resta (2010); Hong and Vanderbilt (2011) of flexoelectricity. To substantiate our arguments, we carefully study the nonanalyticities, due to the long-range character of the electrostatic interactions, that plague the electronic response functions in the long-wave regime. In particular, we devise a rigorous strategy to dealing with this issue by suppressing the macroscopic () component of the self-consistent electrostatic potential in the linear response calculations. We find, however, that such a procedure is not unique – there is an inherent ambiguity in the specification of the “zero macroscopic field” condition in the flexoelectric case, which can be traced back to the choice of an arbitrary reference energy in the periodic crystal. We rationalize such ambiguity by establishing a formal link between the present theory of flexoelectricity and the preexisting theory of absolute deformation potentials. Resta et al. (1990) In addition to providing a solid formal basis to our derivations, our treatment of macroscopic electrostatics also yields a rigorous generalization of the Cochran-Cowley formula to higher orders in , which can be of great utility in future lattice-dynamical studies. Finally, based on our findings, (i) we derive an exact sum rule, relating the flexoelectric coefficients to the macroscopic elastic tensor; (ii) we use such a sum rule to demonstrate that the same definition of the flexoelectric tensor is equally well suited to describing static or dynamic phenomena; (iii) we discuss the physical interpretation of the various physical contributions to the flexoelectric tensor, relating them to earlier first-principles Resta et al. (1990) and phenomenological Morozovska et al. (2012a, b) studies.

This work is structured as follows. In Section II we introduce some useful basic concepts of continuum mechanics, and the general strategy that we use to attack the flexoelectric problem. In Section III we introduce the formalism of density-functional perturbation theory, and the basic ingredients that will be used in the remainder of this work. In Section IV we proceed to performing the long-wave analysis of an acoustic phonon, deriving the piezoelectric and flexoelectric response tensors in terms of the basic ingredients defined above. In Section V we discuss several important properties of the electronic response functions (polarization and charge density), and use them to draw a formal connection to Martin’s theory of piezoelectricity, Martin (1972) and to earlier theories Hong and Vanderbilt (2011); Resta (2010); Tagantsev (1986) of flexoelectricity. In Section VI we study the nonanalytic properties of the aforementioned response functions, obtaining (among other results) a higher-order generalization of the lattice-dynamical theory of Pick, Cohen and Martin. Pick et al. (1970) Finally, Section VII is devoted to discussing the physical implications of the derived formulas, while in Section VIII we briefly summarize our main results and conclusions.

Figure 1: Individual components of the strain gradient tensor in a square two-dimensional lattice. The “longitudinal” (), “shear” () and “transverse” (or “bending”, ) components are linearly independent; a fourth pattern of less obvious physical interpretation () is a combination of and . Direction 1 corresponds to the horizontal axis, 2 to the vertical; the polarization vector is oriented along 2.

Ii Preliminaries

ii.1 Strain and strain gradients

In continuum mechanics, a deformation can be expressed as a three-dimensional (3D) vector function, , describing the displacement of a material point from its reference position at to its current location ,

The deformation gradient is defined as the gradient of taken in the reference configuration,


is often indicated in the literature as “unsymmetrized strain tensor”, as it generally contains a proper strain plus a rotation. By symmetrizing its indices one can remove the rotational component, thus obtaining the symmetrized strain tensor,

is a convenient measure of local strain, as it only depends on relative displacements of two adjacent material points, and not on their absolute translation or rotation with respect to some reference configuration.

In this work we shall be primarily concerned with the effects of a spatially inhomogeneous strain. The third-rank strain gradient tensor can be defined in two different ways, both important for the derivations that follow. The first (type-I) form consists in the gradient of the unsymmetrized strain,


Note that , manifestly invariant upon exchange, corresponds to the tensor of Ref. Hong and Vanderbilt, 2011, and to the symbol of Ref. Tagantsev, 1986. Alternatively, the strain gradient tensor can be defined (type-II) as the gradient of the symmetric strain, ,

invariant upon exchange. It is straightforward to verify that the two tensors contain exactly the same number of independent entries, and that a one-to-one relationship can be established to express the former as a function of the latter and viceversa. For example,


In Fig. 1 we illustrate the three independent components of the and tensors on a square two-dimensional (2D) lattice, evidencing analogies and differences. It is clear from the figure that the longitudinal and shear components are elementary objects in both type-I and type-II forms. The main difference between the two representations concerns the third independent component, which assumes the form of a flat displacement pattern in the type-I form, and has the more intuitive interpretation of a pure bending (one can show that , where is the curvature radius) in the type-II form; the latter will be indicated as transverse strain gradient henceforth. In fact, these three components of the strain gradient tensor are, by symmetry, the only types of independent perturbations in a cubic material, and are therefore very important in the context of flexoelectricity.

ii.2 Long-wavelength acoustic phonons

A macroscopic strain gradient breaks the translational symmetry of the crystal lattice. For this reason, the response to such a perturbation cannot be straightforwardly represented in periodic boundary conditions. This makes the theoretical study of flexoelectricity more challenging than other forms of electromechanical couplings, e.g. piezoelectricity. To circumvent this difficulty, we shall base our analysis on the study of long-wavelength acoustic phonons. These perturbations, while generally incommensurate with the crystal lattice, can be conveniently described Baroni et al. (2001), in terms of functions that are lattice-periodic, and therefore are formally and computationally very advantageous.

Figure 2: Displacement fields in longitudinal (top) and transversal (bottom) sound waves.

The direct relationship between an acoustic phonon and a mechanical deformation is clear from Fig. 2: in the longitudinal and transversal waves one can visually identify regions of negative and positive strain gradients, respectively of the longitudinal and shear type. Mathematically, this observation can be formalized by writing (at the lowest order in the wavevector ) an acoustic phonon as a homogeneous displacement of every material point of the type

where is the displacement amplitude and the frequency. Consider now the microscopic polarization currents (these are due to the displacements of the charged particles, electron and nuclei, from their equilibrium positions) induced by the phonon at the linear-response level, 222 We use a tilde symbol in the following equations to indicate the total relaxed-ion polarization. This includes the effects of the internal strains that are dynamically produced by the deformation field, . In the remainder of this work we shall use , without tilde, to indicate the elementary polarization response function, i.e. in absence of internal strains.

Assuming for the moment that is an analytic function of at the point, in a neighborhood of (i.e. in the long-wavelength regime) we can replace it with its second-order Taylor expansion,


Now, by applying Eq. (1) and Eq. (2), we can compute the local deformation gradient and strain gradient that are associated with the acoustic phonon,


A comparison of Eqs. (4), (5) and (6) suggests that the polarization response to a uniform strain (piezoelectricity) and to a strain-gradient (flexoelectricity) are, respectively, related to the first- and second-order Taylor expansion in powers of of the polarization field produced by a sound wave. The latter can be computed by working with lattice-periodic functions only, implying that all the theoretical and computational weaponry developed so far within periodic boundary conditions (e.g. Bloch theorem, plane wave basis set, pseudopotentials, etc.) can be proficiently applied to the flexoelectric problem. This is precisely the approach that we shall take in the remainder of this work.

While this appears conceptually simple, there are a number of important issues that need to be addressed before one can establish the formal link between the electrical properties of sound waves and the macroscopic electromechanical tensors. First of all, it is not clear a priori whether the above strategy is even applicable. Long-wavelength phonons are generally accompanied by macroscopic electric fields. This is a substantial complication from the operational point of view: the longitudinal character of the electrostatic screening causes a nonanalytic behavior of most response functions (e.g. the atomic eigendisplacements and the electronic polarization, see Section VI for a detailed discussion) in a vicinity of , thwarting their expansion in powers of . Whether (and how) these nonanalyticities can be tamed will need to be assessed prior to starting the actual derivations. Second, phonon eigenmodes also contain, in addition to genuine macroscopic deformations, translations and rotations of a given crystal cell with respect to its reference configuration at rest. It will be therefore necessary to show that such rototranslations do not contribute to the macroscopic electrical response. Third, a phonon is inherently a dynamic perturbation, and whether the effects derived for a sound wave are equally applicable to a static deformation will need to be carefully demonstrated. In the following Sections we shall first introduce the basic ingredients that we need in order to derive the total polarization response, Eq. (4); next, we shall proceed to the formal derivation of the electromechanical tensors, and to their validation in relation to the aforementioned sources of concern.

Iii Density-functional perturbation theory

This Section will provide a brief introduction to the DFPT formalism. This is mainly aimed at specifying the general context of our derivations, as well as at pointing out the key modifications to the standard approach Baroni et al. (2001) that are necessary in the context of this work. In particular, we shall put the emphasis on the following three technical points: the treatment of the macroscopic fields; the definition of the microscopic polarization response; the practical calculation of the relevant response functions by means of publicly available codes.

iii.1 Linear response to monochromatic perturbations

Our starting point is an insulating crystal, whose equilibrium configuration is described by the three primitive translation vectors, , and by a basis of atoms located at positions () within the primitive unit cell. Within density-functional theory, the electronic ground state can be written in terms of the self-consistent (SCF) Kohn-Sham equation,

where is the SCF Hamiltonian at the point in the Brillouin zone, and and are respectively the ground-state Bloch orbitals and eigenvalues. In full generality, the Hamiltonian

contains a single-particle kinetic energy operator, , the external potential of the nuclei, and the Hartree and exchange and correlation potential, the latter depending self-consistently on the electronic charge density ,


( is the occupation of the orbital, equal to 2 if spin pairing is assumed.) The total charge density, , includes the contribution of the nuclear point charges,


where is the bare pseudopotential charge (or the atomic number in the case of an all-electron description), and is a Dirac delta function. [Note that is the equilibrium (unperturbed) atomic position in the crystal, and is a cell index.]

Consider now a monochromatic perturbation, where the atoms in the sublattice undergo a small displacement along of the type


The linear response of the crystal to such a perturbation can be readily computed in the framework of density-functional perturbation theory Baroni et al. (2001), by solving the following Sternheimer equation,


Here are the desired first-order wavefunctions, and are the projection operators on the valence and conduction subspaces, and is the sum of the external perturbing potential (due to the atomic displacements) and the linear variation in the Hxc potential due to the rearrangement of the electron cloud. The arbitrary parameter guarantees orthogonality between and and is otherwise irrelevant. Note that Eq. (10) involves lattice-periodic functions only, and thus provides a convenient route to accessing the relevant response functions at an arbitrary wavevector .

In the context of the present work, we need to focus on three basic response functions, all of which are linear in the perturbation amplitude . (To avoid overburdening the notation, from now on we shall omit the “” prefix whenever the linearity of a given response function with respect to is obvious from the context.) The first quantity is the variation of the total charge density,

Similar to , the cell-periodic function can be also decomposed into an electronic and a (trivial) ionic contribution,

where is the solution of Eq. (10).

The second quantity is the microscopic polarization response, defined as the current density, , that is linearly induced when the perturbation is adiabatically switched on via a time-dependent parameter ,

The variation of can also be written as a cell-periodic part multiplied by a phase,

and decomposed into an electronic and ionic part, . The ionic contribution has again a simple expression,


independent of . It is easy to verify that

The electronic contribution, , is a new quantity that is not part of currently available DFPT implementations. Further details on how it can be calculated in practice are provided in Section III.4.

The third and last basic response function that we shall consider in this work is the force induced on the atom along , , whose cell-periodic part is the -space force constant matrix, ,


This is, of course, a central quantity in DFPT, and can be readily computed following the prescriptions of Refs. Baroni et al., 2001, Gonze, 1997 and Gonze and Lee, 1997. (With respect to the procedure described in these works, note that there is an important subtlety related to the phase of the perturbing potentials and response functions, which we shall discuss in Section III.3.1.)

iii.2 Taylor expansion in a vicinity of

As we shall see in Section IV, in order to obtain the long-wave limit of the polarization response to a phonon, Eq. (4), one needs to evaluate a number of intermediate quantities. These are the lowest terms of the Taylor expansion (in space) of the fundamental response functions , and introduced earlier in this Section. Unfortunately, these functions are plagued by a nonanalytic behavior at , which implies that their direct Taylor expansion is not feasible. The nonanaliticity is related to the macroscopic electric fields that occur in response to the perturbation. To clarify this point, it is useful to write the induced electric field as

After expanding the cell-periodic part, into its reciprocal-space coefficients (indexed by the reciprocal-lattice vectors ), it becomes apparent that the term (indicated by a wide bar symbol),

which is purely longitudinal, is problematic for . In fact, one can show (a rigorous derivation is provided in Section VI.2) that, at order zero in , is a direction-dependent constant,


where is the macroscopic dielectric tensor, is the Born dynamical charge tensor, and . Such a nonanalytic behavior of propagates to the charge, polarization and lattice responses, thwarting their Taylor expansion at .

We shall circumvent this difficulty by removing the macroscopic electrostatic component (corresponding to the vector of the reciprocal lattice) from the self-consistent electrostatic potential. This prescription has the effect of screening the longitudinal fields associated with the long-wavelength phonon. Therefore, it corresponds to adopting short-circuit electrical boundary conditions in the calculation of the response functions, which is indeed the standard convention in the definition of the electromechanical coupling coefficients. This way, we have solved two problems at once: i) all the response functions become analytic at and their polynomial expansion is, in principle, well defined at any order in ; ii) we have specified once and for all that the response functions are calculated with a macroscopic electric field kept constant and equal to zero, i.e. in short circuit. A formal demonstration of these claims, based on the dielectric matrix approach, Pick et al. (1970) is provided in Section VI.

By using the aforementioned precautions, it is now formally possible to perform the Tayor expansion of the charge density response (we shall assume from now on that repeated indices are implicitly summed over),


the microscopic polarization,


and the force-constant matrix,


Note the choice of the prefactors, which is motivated by the relationship to the localized real-space representation (see Section V). In practice, at an arbitrary order and for a given response function , we define


This prescription also guarantees that the functions are always real.

iii.3 Practical considerations

iii.3.1 Phase factors

Our definition of the elementary monochromatic perturbations, Eq. (9), differs from that used by Gonze Gonze (1997) and Gonze and Lee Gonze and Lee (1997) (GL),

by a sublattice-dependent (but cell-independent) phase factor,

Such a modification is irrelevant in the calculation of phonon dispersion curves, but is crucial in the context of the long-wave expansion performed here. In fact, it guarantees that the acoustic phonon eigenmodes do not depend on the (arbitrary) assignment of each basis atom to a given cell in the crystal, and therefore we regard it as a very natural choice on general physical grounds.

From the point of view of practical calculations, it should be kept in mind that all the response functions discussed in this work generally differ from the quantities that are computed within the publicly available implementations of DFPT. Given that the modification consists in a trivial phase, however, it is easy to write the correspondence between the response functions defined by GL and those considered here. For example, concerning the charge density response, we have (by using the linearity of the response functions in the perturbation)

In the case of the force-constant matrix, there is an additional phase factor coming from the factorization Eq. (12), which leads to the following correspondence,

Of course, the real-space force constants (i.e. the second derivative of the total energy with respect to the displacements of individual atoms) must be consistent with the definition given by GL,

Therefore, our modification essentially concerns the definition of the Fourier transform that is used to move between direct and reciprocal space. Here we have [compare with Eq. (10) of Ref. Gonze and Lee, 1997]

iii.3.2 Differentiation in -space

To calculate the Taylor expansion of the response functions one can follow two different routes. Ideally, it would be desirable to take the analytical gradients of the Sternheimer equation, Eq. (10), in -space and solve directly for the perturbed wavefunctions at a given order in ,

(The dependence on the sublattice index and the displacement direction has been kept implicit to avoid overburdening the notation.) Then the response functions could be simply calculated from the orbitals at the desired order in by exploiting the linearity of the respective Taylor expansions. For example, the charge density at linear order in would read


We have not implemented the analytic long-wave expansion of the Sternheimer equation here. (The explicit derivation is under way and will be the subject of a future communication.) Instead we propose, for the time being, to extract the needed Taylor-expanded response functions by using a finite-difference approach in -space. 333 Note that the quantities that undergo the differentiation in are first-order response functions in the perturbation amplitude, , to be calculated within DFPT. This has the advantage of allowing the calculation of the flexoelectric tensor in arbitrary solids by means of the existing implementations of DFPT. In practice, it suffices to discretize Eq. (17) (replace with the desired response function), by using an appropriate grid of points surrounding . This procedure is good for performing the long-wave analysis of the charge density and the force-constant matrix, as both quantities are fully implemented in publicly available codes. The calculation of the polarization response deserves a separate comment, as currently available implementations of DFPT provide access only to the macroscopic (cell-averaged) part, and not to the full microscopic current density. In the following Section we shall outline a viable procedure to access the latter quantity.

iii.4 Microscopic polarization response

To derive the electronic contribution to the microscopic polarization response, we shall work in reciprocal space, and write in terms of its Fourier coefficients (we shall omit the superscript “el” in the remainder of this section, as the absence of the ionic contribution is obvious from the context),


We seek an (unknown) operator such that


where and . (Strictly speaking, only the component of the polarization response, corresponding to , is sufficient for the scopes of the present work; we keep the -dependence for the sake of generality.) It is convenient to simplify the notation, and write


where the index runs over valence (occupied) wavefunctions. As the first-order wavefunctions belong, by construction, to the conduction manifold, one can insert a projector ,

where runs over the unoccupied orbitals. Since both and are eigenstates of the unperturbed Hamiltonian, , one can readily write


Now, assuming that does not depend explicitly on time we have, from Ehrenfest theorem,


where stands for the expectation value of the operator . Since in a nonmagnetic insulator , it follows that the commutator in Eq. (22) must correspond to the current density operator,


Hence, we have




are the first-order orbitals induced by as a perturbing operator. These can be conveniently obtained by solving the nonselfconsistent Sternheimer equation,


This result allows us to calculate the cross-gap matrix elements of the unknown microscopic polarization operator, , by means of the more familiar current density operator. The probability current is a fundamental quantum-mechanical observable, and implementing it in an existing DFPT code should not present major conceptual obstacles; such a task will be the topic of a future communication. In this context, it is worth mentioning the work of Umari, Dal Corso and Resta Umari et al. (2001), where the microscopic polarization response to a uniform electric field perturbation was derived and computed; the authors used an approach that is closely related to the one presented here.

Iv Long-wave analysis

In the following, we shall use the long-wave method Born and Huang (1954) to derive the electromechanical response (either piezoelectric or flexoelectric) of the crystal in terms of the elementary ingredients introduced above. We shall first focus on the atomic displacements induced by a “short-circuited” (in the sense specified in the previous Section) acoustic phonon, and later compute the polarization field associated with the deformation.

iv.1 Internal strains

Consider the real-space atomic equation of motion,

where are the displacements, is the real-space force-constant matrix and is the mass of the specie . ( indexes the lattice cell where the atom is located; the atom is located at the cell; and refer to Cartesian directions.) We seek solutions of the type

These are given by the eigenvalue problem


We solve Eq. (28) perturbatively Tagantsev (1986); Born and Huang (1954) in a vicinity of by writing the wavevector as Born and Huang (1954) , where is a dimensionless perturbation parameter. For an acoustic branch, and can be expanded as follows, Born and Huang (1954)


(The expansion of starts with the linear term, as for acoustic waves the frequency approaches zero as .) We shall now proceed to calculating the induced displacements by plugging Eq. (29), Eq. (30) and the Taylor expansion of the force-constant matrix, Eq. (16), into Eq. (28), and by grouping the different terms according to their perturbative order.

At order zero in we have


i.e. the phonon eigenvector must be independent of . In fact, the matrix is the zone-center dynamical matrix of the crystal, and has three zero-frequency eigenmodes corresponding to the rigid translations of the whole lattice along each Cartesian direction. This means that Eq. (31) is identically verified by any real-space vector .

At first order in we have


Solvability requires that

be identically satisfied for any , which is indeed the case Born and Huang (1954). The explicit solution can be written as


Here we have introduced as the pseudoinverse Wu et al. (2005) of the singular matrix (the zero eigenvalues of , corresponding to rigid translations, are mapped into zero eigenvalues of , while the nonsingular remainder of the matrix is inverted), and


is the piezoelectric force-response tensor (following the notation of Ref. Wu et al., 2005). describes the force induced on the sublattice along when the crystal undergoes a homogeneous strain deformation, 444 Strictly speaking, is defined as the response to an unsymmetrized strain, . However, due to the invariance of with respect to exchange, we can readily identify it as the force response to a symmetrized strain, . , and is symmetric with respect to exchange. Born and Huang (1954) The internal-strain tensor Martin (1972), , describes the atomic relaxations induced by , and inherits the invariance from . Note that is specified only modulo an arbitrary -independent (but possibly -dependent) constant, which physically corresponds to a rigid shift of the whole lattice.

At second order in we obtain


where are atomic masses, and we have introduced the type-I flexoelectric force-response tensor as follows,


The square brackets and round brackets are defined (in loose analogy with the notation of Ref. Born and Huang, 1954) as


describes the force induced along on a given atomic sublattice by a “frozen-ion” (in the sense specified in Ref. Hong and Vanderbilt, 2011) strain gradient (in type-I form). describes the additional force produced by the atomic relaxations that are first-order in and is, therefore, only relevant to crystals that have one or more free Wyckoff parameters. Note that the round bracket is a type-II object (i.e. it relates the force along to the type-II strain gradient component ), hence the symmetrization in Eq. 37 (the tensor is a type-I object).

The linear problem, Eq. (36), admits solution if and only if the following condition on is satisfied Born and Huang (1954); Tagantsev (1986),


where is the total mass of the primitive cell, and We shall see in Section IV.3 that the quantity can be considered a “type-I” representation of the macroscopic elastic tensor, 555 Note that in piezoelectrically active materials the long-range electrostatic fields contribute to sound propagation via a nonanalytic term in the elastic tensor. Our neglect of the macroscopic term in the electrostatics implies that such contribution is not accounted for in Eq. (40). One will then easily recognize Eq. (40) as the sound wave equation. Born and Huang (1954) Its solutions depend only on , on the mass density and on the propagation direction , and are therefore characterized by a linear dispersion relation along any given . By combining Eq. (40) with Eq. (36), we readily obtain


where is the type-I flexoelectric internal-strain tensor, and we have also introduced the mass-compensated force-response tensor,


Note that the sum over of the tensor identically vanishes by construction; this is a necessary condition for the linear problem Eq. (42) to be solvable, proving that our derivations are internally consistent.

In summary, the lattice response to a (short-circuited) long-wavelength acoustic phonon can be written as


where and are the desired internal-strain tensors.

Before closing this part, it is useful to briefly comment on the relationship between our derivation and Tagantsev’s. Our approach accurately follows the formalism of Ref. Tagantsev, 1986, except for the procedure to extract the relevant physical quantities from the force constants of the crystal. 666The force-response tensor of this work is closely related to the of Ref. Tagantsev, 1986, with the exception that the latter is not summed over one of the sublattice indices. Note that our Eq. (37) differs from the unnumbered equation after Eq. (A7) of Ref. Tagantsev, 1986: in the former, the contribution that depends on (round brackets) is correctly symmetrized over the indices, while in the latter it is not. Tagantsev wrote the moments of the force constant matrix as lattice sums in real space, whose convergence is not guaranteed unless a specific prescription for dealing with macroscopic electrostatics is formulated. A heuristic treatment of the macroscopic fields might be possible in atomistic models, where the charge response to individual atomic displacements is trivially simple. In the present quantum-mechanical context, the problem is complicated by the presence of higher-order multipolar interactions Resta (2010); Hong and Vanderbilt (2011), whose impact on lattice dynamics might be cumbersome to keep track of. Here we solve this issue by working in space, where a rigorous strategy to suppress the problematic electric fields in the long-wavelength limit is easy to implement once and for all, and does not require any special effort.

iv.2 Macroscopic flexoelectric coefficients in type-I form

In order to write the total polarization response to the long-wavelength phonon (and hence to an arbitrary mechanical deformation), we need to combine the eigendisplacements derived in the previous section with the small- expansion of the induced polarization, Eq. (15). As in this work we are only concerned with the macroscopic response, we shall work on the cell-averaged counterparts of the polarization-response functions, which we indicate by an overline symbol,


By construction, the zero-order term is proportional to the Born effective charge tensor of the specie ,


It follows that, due to the acoustic sum rule Pick et al. (1970), the sum over of the tensor identically vanishes. Thus, the contribution of the rigid translation to the macroscopic vanishes as well (as anticipated above), leaving us with the terms that are linear and quadratic in ,


where the relaxed-ion response tensors are given by


In the above expressions we have used the bar symbol to indicate, following the notation of Ref. Wu et al., 2005, the frozen-ion counterparts of the tensors,


The tensors and correspond to the well-known relaxed-ion and frozen-ion piezoelectric coefficients, and are both symmetric with respect to exchange. (This property of the latter tensor, as well as its relationship to Martin’s theory of piezoelectricity Martin (1972), will be rigorously demonstrated in Sections V.3 and V.4, respectively.) Hence, the unsymmetrized stress tensor, , can be replaced with its symmetric counterpart in Eq. (47), leading to an expression that is fully invariant with respect to either translations or rotations of the original reference,


This formula, which is a central result of this work, allows us to identify and with the sought-after relaxed-ion and frozen-ion flexoelectric tensors, respectively, i.e.


The superscript I indicates that and are “type-I” objects, i.e. they describe the response to a type-I strain gradient tensor . One can, of course, write the flexoelectric tensor, , in type-II form – we shall do this explicitly hereafter, as such a conversion is important for later derivations (and, in particular, for tracing the important link to elasticity that we anticipated in Section IV.1).

iv.3 Type-II form and “elastic sum rule”

From the definition of the type-II strain gradient tensor, , it follows that


where the type-II flexoelectric tensor is related to via a cyclic permutation of the last three indices,


Note that is invariant upon exchange of the last two indices, consistent with the analogous symmetry of the type-II strain gradient tensor. By combining Eq. (55) with Eq. (49) we have


The tensor is defined from via the symmetrization Eq. (55). The internal-strain tensor follows from via an analogous operation on the indices , which can be ultimately traced back to a redefinition of the force-response tensor,


In turn, the tensor can be written explicitly in terms of the type-II flexoelectric force-response tensor,


after separating the mass-dependent part,


Under the assumption that the crystal at rest is free of stresses (see Section 28 of Ref. Born and Huang, 1954),