PHYSICAL REVIEW X 12, 031002 (2022) Translational Covariance of Flexoelectricity at Ferroelectric Domain Walls 1,2,* 3,4,† Oswaldo Diéguez and Massimiliano Stengel 1 Department of Materials Science and Engineering, Faculty of Engineering, Tel Aviv University, Tel Aviv 6997801, Israel 2 Theory and Simulation Group, Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC, BIST, Campus UAB, Bellaterra, Barcelona 08193, Spain 3 Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain 4 ICREA - Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona, Spain (Received 1 February 2022; revised 23 May 2022; accepted 3 June 2022; published 5 July 2022) Macroscopic descriptions of ferroelectrics have an obvious appeal in terms of efficiency and physical intuition. Their predictive power, however, has often been thwarted by the lack of a systematic procedure to extract the relevant materials parameters from the microscopics. Here we address this limitation by establishing an unambiguous two-way mapping between spatially inhomogeneous fields and discrete lattice modes. This yields a natural treatment of gradient couplings in the macroscopic regime via a long- wavelength expansion of the crystal Hamiltonian. Our analysis reveals an inherent arbitrariness in both the flexoelectric and polarization gradient coefficients, which we ascribe to a translational freedom in the definition of the polar distortion pattern. Remarkably, such arbitrariness cancels out in all physically measurable properties (relaxed atomic structure and energetics) derived from the model, pointing to a generalized translational covariance in the continuum description of inhomogeneous ferroelectric structures. We demonstrate our claims with extensive numerical tests on 180° domain walls in common ferroelectric perovskites, finding excellent agreement between the continuum model and direct first- principles calculations. DOI: 10.1103/PhysRevX.12.031002 Subject Areas: Computational Physics Condensed Matter Physics Materials Science I. INTRODUCTION sufficiently strong [7–9]. On the other hand, flexoelectricity endows the spatial gradients of the main order parameters Spatially inhomogeneous structures in ferroelectrics with potentially useful functionalities, e.g., a spontaneous such as domain walls, vortices, etc., have been the subject polarization at ferroelastic twin boundaries [10,11], and a of intense research in the past few years [1,2], because of spontaneous strain at ferroelectric walls via the converse their emerging physical properties and nontrivial topology. effect [12,13]. Considerable efforts are currently directed at identifying In light of these findings, improving our understanding the physical mechanisms that govern the stability of the of the interplay between flexoelectricity and ferroelectricity observed patterns and their response to external probes. In appears as essential for future progress. Given the advances addition to the well-known factors related to the electrical that the first-principles theory of flexoelectricity has made and mechanical boundary conditions, flexoelectricity since the pioneering works of Resta [14] and Hong et al. (describing the coupling between polarization and strain [15], such a goal appears now well within reach. As of gradients [3–6]) has been receiving increasing attention in early 2020, a complete calculation of the bulk flexoelectric this context. On one hand, the flexoelectric coupling tensor can be carried out [16,17] with the latest release of the contributes substantially to the gradient energy, to the point publicly distributed ABINIT [18,19] package, providing in that a transition to a modulated phase may occur if principle a solid theoretical reference for the interpretation of the experimental data [13,20]. Unfortunately, these studies * dieguez@post.tau.ac.il have also revealed that the bulk flexoelectric coefficients are † mstengel@icmab.es ill defined as stand-alone material properties. More specifi- cally, their definition is plagued by two distinct ambiguities, Published by the American Physical Society under the terms of which are respectively related to the treatment of elastic the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to and electrostatic fields in the long-wavelength expansion the author(s) and the published article’s title, journal citation, [17,21,22]. This fundamental limitation prevents a straight- and DOI. forward comparison between theoretical and experimental 2160-3308=22=12(3)=031002(19) 031002-1 Published by the American Physical Society OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022) results, as further considerations are needed to make sure that yields unique answers for the domain-wall structure once the calculated values relate to what is being measured in a the local field amplitudes are converted back into atomic physically meaningful way. distortions, enabling a straightforward validation of the If they only concerned the specifics of how the flexo- method against direct density-functional theory (DFT) electric effect is defined and measured, these issues would calculations. We illustrate this point by calculating 180° be of limited importance. On the contrary, the arbitrariness ferroelectric walls in six different perovskite materials, of the coupling coefficients is problematic in a much finding answers that are within 10%–20% of the “exact” broader context, as it questions the validity of the widely result. Given the extreme (one-cell thick) abruptness of the popular Landau-Ginzburg-Devonshire (LGD) theories of structures, we regard this as a severe test for a continuum ferroelectrics, a cornerstone of the theoretical understand- approach, and such an accuracy is exceptionally good. By ing of inhomogeneous polar structures for several decades. calculating domain walls under hydrostatic pressure in At a domain wall, standard LGD models predict [12,13] a BaTiO3 we also demonstrate the exactness of our theory in dependence of both the energy and structure on the the limit of smooth domain walls. flexocoupling coefficient via the converse effect, which From our results, a new paradigm emerges in the associates a uniform strain with a gradient of the polar order construction of continuum models of ferroics: the invari- parameter. The obvious question is then, how can we trust ance of the Landau-Ginzburg-Devonshire free energy with such physical predictions once we know that one of the respect to a number of generalized gauge transformations materials properties on which they depend is ill (or at least of the parameters and fields. This implies abandoning the nonuniquely) defined? widespread belief that such parameters (e.g., the flexo- To formulate the problem on firm theoretical grounds, electric coefficient) and fields (e.g., the local strain) be the first challenge consists in establishing a rigorous well-defined physical properties of the crystal. It also two-way mapping between microscopic degrees of free- emphasizes the need for an intimate connection between dom and macroscopic order parameters. In the case of microscopics and macroscopics in order to achieve a spatially homogeneous crystal phases, such a task poses qualitatively sound picture. limited conceptual issues: Building effective low-energy This work is organized as follows. In Sec. II we address Hamiltonians in terms of the physically relevant lattice the theoretical issues that arise in the design of continuum distortions (in perovskite crystals these typically include models, focusing on the aforementioned ambiguities in polarization, strain, and antiferrodistortive oxygen tilts) is the definition of local strains, polar distortion patterns, now common practice [23–26] within the ab initio com- and flexoelectric coefficients. In Sec. III we present our munity. Whenever these degrees of freedom are no longer numerical tests on 180° domain walls in perovskite oxides, constant over space, however, many subtleties arise, and the alongside a detailed validation against the results of direct partition of the energy into different macroscopic contri- butions generally becomes nonunique [27]. The question, first-principles calculations. In Sec. IV we discuss the then, is, are there specific criteria for ensuring that the result implications of our findings in the context of the relevant is physically meaningful? And, once we have solved the literature. We summarize our work and present our con- continuum equations, how can we verify that our solution clusions in Sec. V. is consistent with the “training model,” i.e., our first- principles engine? II. THEORY Here we show, by deriving the continuum equations and parameters via a rigorous long-wavelength approximation A. Statement of the problem of the first-principles lattice Hamiltonian, that the above Consider a 180° domain wall in a ferroelectric crystal as difficulties can be traced back to a translational freedom in schematically illustrated in Fig. 1. The outer extremes of the definition of the polar distortion pattern. As a conse- both panels correspond to the oppositely oriented ferro- quence of such freedom, not only the flexocoupling f but electric domains, where the polarization (P) saturates to its also the polarization gradient coefficient G (entering the bulk value; in the central domain-wall region P transitions continuum functional via the squared gradient of the from negative to positive values when moving from left to polarization field) is affected by an inherent arbitrariness right. The atomic structure far from the wall is well in its definition. Crucially, we find that the respective described in terms of a periodically repeated crystalline ambiguities in f and G cancel out exactly in any physical cell. Such a structure can be readily obtained from a bulk prediction of the continuum model, implying that a con- calculation: one typically starts from the reference con- sistent treatment of both terms is essential for the overall figuration, breaks the centrosymmetry by hand (e.g., by theory to work. displacing atom A upward or atom O downward by a small Of particular note, the aforementioned arbitrariness amount), and lets the atoms relax to their polar ground state. directly affects the definition of the elastic displacement Note that one can perform the relaxation while fixing either field (and hence the strain), which we find to be nonunique. A or O to their original locations; while the resulting Our long-wave approach to continuum theory, however, distortions dA;O differ [compare Figs. 1(a) and 1(b)], the 031002-2 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) (a) (quantified here by dκ, with κ ¼ A; O) is ill defined, the u dA = 0 dO amount of elastic displacement Δu that we must incorpo- rate to obtain the correct registry between the left and right half-lattices is also ill defined. This is obvious by looking at Figs. 1(a) and 1(b), where we compare two different choices for dκ and Δu. Clearly, the total atomic distortions within each domain, given by ðdκ þ Δu=2Þ, are the same in Figs. 1(a) and 1(b). And yet, what we understand as the “macroscopic elastic offset” between the domains, Δu, markedly differs. This situation is paradoxical in light of (b) the widespread assumption that the local elastic displace- u ment field (and hence the local strain) be a physically dO = 0 unambiguous degree of freedom of the crystal. The sketch of Fig. 1 seems to disprove such an interpretation. As we dA have anticipated in the Introduction, abandoning such a long-established paradigm brings about a number of con- ceptual troubles; we address them in the remainder of this work. FIG. 1. Schematic model of a [100]-oriented ferroelectric B. Macroscopic theory domain wall in a perovskite crystal. For simplicity, only atoms To frame our discussion, in this section we recap the in AO planes perpendicular to the wall are shown as red circles established [12] macroscopic theory for an Ising-like 180° (A) and black squares (O), respectively. The elastic displacement ferroelectric domain wall as illustrated in Fig. 1. For of the crystal cells (dashed lines) is also shown. The ferroelectric simplicity, we shall exclusively focus on the transverse distortion of the lattice deep within the domain is indicated by red and black arrows, respectively, for the displacements (dκ ) of the A (y) component of the polarization: Longitudinal (x) com- and O sublattice with respect to their high-symmetry positions. ponents are typically small in perovskite ferroelectrics, and The macroscopic shift of the crystal cell Δu is also shown. would require an explicit treatment of the electrostatic Displacements are exaggerated for illustration purposes. Panels energy, which is nontrivial in the flexoelectric case [27]. We (a) and (b) illustrate how two different choices for dκ and Δu can also restrict our attention to a single strain component, the be made for the same atomic structure. xy shear, which is the most relevant one to our present scope. (This implies neglecting the tetragonal distortion of the crystal cells deep within the domains.) Within these two structures are related by a rigid translation of the whole assumptions, the simplest free-energy functional to lattice, and are therefore degenerate. describe the problem is Once the domain walls form, the two oppositely oriented domains no longer enjoy translational invariance sepa-   C 2 A 2 B 4 ∂P G ∂P 2 rately: the wall lifts the indeterminacy modulo a rigid shift FðP; εÞ ¼ ε þ P þ P þ fε þ : ð1Þ of the cell, and uniquely sets the registry between the two 2 2 4 ∂x 2 ∂x oppositely polarized half-lattices. The shift that we must add to the calculated bulk atomic structure to correctly align Here P is the parallel (y) component of the polarization, the two semi-infinite regions can conveniently be ration- while the normal to the wall is indicated as x; alized in terms of the elastic displacement field uðxÞ, which undergoes a jump Δu when moving across the wall along ∂u ε¼ ð2Þ the normal direction (x). Thus, the spontaneous alignment ∂x between the two domains can be understood physically as an electromechanical effect, where the elastic degrees of is the shear (xy) component of the strain, defined as the x freedom emerge as a secondary consequence of the primary derivative of the parallel (y) component of the displacement polar distortion of the lattice. And indeed, recent works field u, C is the elastic constant (C44 component), A and B [12,13] have clarified that the net elastic displacement is are the standard parameters of the homogeneous Landau due to converse flexoelectricity, and Δu can be related potential, f is the flexocoupling coefficient, and G is the linearly to the flexocoupling coefficient(s) of the crystal in polarization gradient coefficient. Note that the flexoelectric its cubic reference phase. coupling is written here in terms of the converse effect Earlier works, however, have overlooked the central (uniform strain in response to a P gradient); it is related to conceptual issue with the above interpretation. Since the the direct effect and to the standard Lifshitz-invariant ferroelectric distortion pattern within the bulk domains formula via simple integrations by parts: 031002-3 OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022)   ∂P ∂ε f ∂P ∂ε The domain-wall energy per unit area can be obtained fε ≃ −f P ≃ ε − P : ð3Þ ∂x ∂x 2 ∂x ∂x by integrating the free-energy density along the normal to the wall, (The difference between the three expressions consists in Z surface terms, irrelevant for the present bulk theory.) W¼ dx½F(PðxÞ) − F0 ; ð11Þ By imposing the stationary condition, ∂F where F0 is the energy density of the monodomain ground ¼ 0; ð4Þ ∂ε state. One arrives at Z the strain can be integrated out, which immediately leads to A2 þ∞ A2 4 8 W¼ ξ dx sech4 ðxÞ ¼ ξ ¼ jF jξ; ð12Þ the following result [12]: 2B −∞ 2B 3 3 0 f ∂P f where F0 ¼ −A2 =4B is the bulk energy density. [Since F0 ε¼− → u ¼ − P: ð5Þ C ∂x C is an energy per unit volume, Eq. (12) correctly describes W in units of energy per length squared]. (We have eliminated a trivial integration constant in u by An interesting consequence of the above derivations is imposing that u vanishes whenever P ¼ 0.) Therefore, the that the main physical properties of the wall (the thickness ξ displacement field at a ferroelectric domain wall adopts and the energy W) depend on the inhomogeneous coef- the exact same spatial profile as the polarization, except ficients only via the renormalized gradient coupling G̃, for the −f=C scaling factor. From Eq. (5) one then can while other features (the elastic offset Δu) explicitly extract the net elastic offset Δu that we introduced in the depend on the flexocoupling f. This means that we can previous section: replace f and G with arbitrary numbers, provided that P0 f G − f 2 =C retains the original value, and extract the same Δu ¼ −2 : ð6Þ physical answers for ξ and W; Δu, on the other hand, is not C invariant with respect to such a transformation. This After eliminating the strain, we obtain the following property of Eq. (1) may appear at first sight as a math- simpler expression for the free energy: ematical curiosity, but has profound physical implications in relation to the paradox illustrated in Fig. 1; we explore   A 2 B 4 G̃ ∂P 2 f2 them in the following sections. FðPÞ ¼ P þ P þ ; G̃ ¼ G − : ð7Þ 2 4 2 ∂x C C. Direct mapping to the microscopics The condition for stability is that the renormalized polari- To test the validity of Eq. (1) in real systems, we need to zation gradient coefficient be positive, GC > f 2 . (This establish a microscopic interpretation of the order param- criterion is well known: if f is large enough, the system eters entering the continuum functional. In full generality, becomes unstable and an incommensurate transition to a we use the following two-step procedure. First, we express modulated state may occur [7–9].) The equation of state is the individual atomic displacements as continuum func- given by the stationary condition with respect to P: tions of the real-space coordinates r via a linear trans- formation of the relevant vector fields vνβ ðrÞ: ∂2 P AP þ BP3 − G̃ ¼ 0: ð8Þ X ∂x2 uκα ðrÞ ¼ vνβ ðrÞT νβ;κα : ð13Þ νβ We shall attempt a trial solution of the type   (T νβ;κα is the transformation matrix of the mapping, describ- x PðxÞ ¼ P0 tanh : ð9Þ ing the displacement of the sublattice κ along α for a unit ξ amplitude of vνβ , where β runs over the Cartesian compo- After a few steps of straightforward algebra, we arrive at nents of the vector field vν .) Second, we write the displace- ment of the atom in the lth cell by sampling the atomic ð0Þ A 2G̃ displacement fields at the undistorted lattice sites Rlκ : P20 ¼ − ; ξ2 ¼ : ð10Þ B jAj ð0Þ ulκα ¼ uκα ðRlκ Þ: ð14Þ P0 is determined by the bulk Landau potential, while ξ is a length, and has the obvious physical meaning of domain- The combination of Eqs. (13) and (14) endows the wall thickness. order parameters with the meaning of envelope functions, 031002-4 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) generally smooth on the scale of the interatomic spacings, D. Calculation of the coupling coefficients that modulate a cell-periodic displacement pattern of the The homogeneous coefficients A and B are easy to atoms over the volume of the crystal. extract from a first-principles calculation: they are readily We assume a perovskite-structure lattice henceforth and given by a quartic fit of the energy of the primitive cell as a focus on two specific order parameters, ν ¼ u, P, corre- function of the distortion amplitude along the direction (in sponding to the elastic displacement and ferroelectric ðβÞ polarization. We identify their respective blocks of the configuration space) spanned by pκ . The gradient terms transformation matrix T with the threefold degenerate (especially f and G) are technically more challenging to acoustic and “soft” transverse modes of the undistorted calculate, in that they are defined in terms of spatially cubic structure at the zone center. Within the standard modulated (and hence non-cell-periodic) atomic distortion choice of the coordinate axes, both blocks are diagonal on patterns. Recent developments [16,22,27] in density- the Cartesian indices, functional perturbation theory have overcome these diffi- culties by applying the long-wavelength method to the ðβÞ ðβÞ phonon problem. We show in the following that Eq. (13) T uβ;κα ¼ uκ δαβ ; T Pβ;κα ¼ pκ δαβ : ð15Þ directly connects to the formalism of Refs. [22,27], and hence leads to a physically sound definition of C, f, and G. ðβÞ ðβÞ Since all gradient terms are harmonic, we consider a uκ and pκ are sets of three five-dimensional basis vectors, each forming a T 1u irreducible representation of linear-response regime in the field amplitude with respect to the high-symmetry cubic phase. To capture the spatial the Pm3̄m point group. Regardless of the microscopics, ðβÞ modulation, it is convenient to work in Fourier space and uκ depends neither on the sublattice index κ nor on the express the relevant perturbations of the continuum fields Cartesian index β, as it describes a rigid shift of the cell that as a constant times a complex phase, e.g., PðrÞ ¼ Pq eiq·r . is collinear with uβ [27]: Via Eq. (14), the corresponding lattice distortions can be written as linear combinations of monochromatic displace- ðβÞ uκ ¼ 1: ð16Þ ment patterns of the atoms: ð0Þ Equation (13) reduces then to the following simplified ulκα ¼ uqκα eiq·Rlκ : ð19Þ expression, The second derivatives of the energy with respect to uqκα ðαÞ define [22] the force-constants matrix at the specified wave uκα ðrÞ ¼ uα ðrÞ þ Pα ðrÞpκ ; ð17Þ vector q, ðαÞ where the only remaining task consists in specifying pκ . ∂2 E Φqκα;κ0 β ¼ ; ð20Þ As detailed in Sec. II F, we require that the homogeneous ∂u−q q κα ∂uκ0 β solution of the Landau potential, Eq. (1), reproduce the ðβÞ spontaneous bulk ferroelectric distortion pattern dκ via which provides the formal link to the established density- Eq. (17): functional perturbation theory framework [22]. The macroscopic limit of Eq. (20) is taken via a long- ðβÞ dκ ¼ P0 pκ : ðβÞ ð18Þ wave expansion [17,22] in powers of q: ð0Þ ð1;γÞ qγ qδ ð2;γδÞ ðβÞ Φqκα;κ0 β ¼ Φκα;κ0 β − iqγ Φκα;κ0 β − Φκα;κ0 β þ    : ð21Þ This condition implies that pκ has four independent 2 entries and depends on β by a permutation of the oxygen indices [21]. Equation (18) does not lead to a unique The zeroth order term is the usual zone-center force- ðβÞ constants matrix in short-circuit electrical boundary con- solution for pκ , though: (i) There is a (trivial) freedom in ðβÞ ditions. (It may be used to compute the homogeneous the choice of the unit in which both PðrÞ and pκ are quadratic coefficient A.) The first-order term vanishes in expressed, and (ii) because of the translational invariance the cubic perovskite reference structure. Finally, the sec- ðβÞ that we have mentioned earlier, dκ is only defined modulo ond-order term allows one to extract the sought-after an arbirary shift of the whole lattice. In the following, we information about flexoelectricity, polarization gradient, assume that some choice has been made for (i) and (ii) and and elasticity via a projection onto the elastic and polar proceed to deriving all the coefficients entering Eq. (1) in displacement patterns. In particular, in our specific context terms of microscopic quantities; later on, we discuss the of the [100]-oriented wall with the polarization oriented implications of (ii) in regards to the apparent paradox along [010], Eqs. (16) and (15) lead to the following of Fig. 1. explicit formulas, 031002-5 OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022) 1 X ð2;xxÞ coefficients within microscopic theory, and is directly C¼− Φ 0 ; ð22aÞ 2Ω κκ0 κy;κ y related to the paradox of Fig. 1. After a quick glance at Eq. (22), it is not difficult to see 1 X ð2;xxÞ ðyÞ where this arbitrariness may come from. Of the two f¼− Φ 0 p0 ; ð22bÞ ingredients that enter the definition of C, f, and G, the 2Ω κκ0 κy;κ y κ Φð2;xxÞ matrix is well defined, as it directly emerges from a 1 X ð2;xxÞ ðyÞ ðyÞ long-wave expansion of the force constants. As we have G¼− Φ 0 pκ pκ0 ; ð22cÞ ðyÞ 2Ω κκ0 κy;κ y anticipated earlier, however, the basis vector pκ is only defined modulo a rigid displacement of the whole lattice ðyÞ where Ω is the volume of the undistorted primitive cell. (Fig. 1). In other words, we can always replace pκ with ðyÞ To connect with the existing first-principles theory of any vector that differs from pκ by a κ-independent flexoelectricity, it is useful to recall the definition [17,22] of constant. There is no fundamental symmetry principle that the force-response coefficient, favors one choice over the other—it is entirely a matter of convention. Let us see what happens if we operate such a 0 1 X ð2;xxÞ fκ ¼ − Φ 0 ; ð23Þ transformation, by defining a new basis vector as 2 κ κy;κ y ðyÞ 0 ðyÞ pκ ¼ pκ þ λ: ð25Þ as the force on the sublattice κ 0 produced by a gradient of the shear strain. A comparison between Eqs. (22) and (23) Evidently, the homogeneous coupling coefficients are shows that the flexocoupling coefficient is properly defined unaffected, since the energetics of the uniform phase is here as the geometrical force on the polar mode produced P 0 ðyÞ insensitive to the choice of the origin. This is not the case by a strain gradient, f ¼ κ0 f κ pκ0 =Ω, consistent with for the gradient coefficients, whose transformation rules earlier works [27,28]. Also, the definition of the elastic can be straightforwardly derived by plugging Eq. (25) into constant C is consistent with the classic result of Born and Eq. (22): Huang [29], as revisited recently in a modern electronic- structure context [17,22]. Note that f κ and C comply with f 0 ¼ f þ λC; ð26aÞ the “elastic sum rule” [17,22], 1 X κ0 G0 ¼ G þ 2λf þ λ2 C: ð26bÞ C¼ f ; ð24Þ Ω κ0 The fact that both the flexoelectric and gradient coefficients as can be easily verified via Eq. (22). depend on an arbitrary constant λ is actually easy to The fact that we obtain the elastic tensor, which is rationalize on elementary physical grounds. Introducing usually regarded as a homogeneous coupling parameter, via a shift in the polar distortion is harmless in the homo- a similar procedure as spatial dispersion coefficients such geneous case, but the gradient of P comes with an extra as G and f, might come as a surprise. This is justified by strain field, which contributes to both the elastic and the fact that a uniform strain is a gradient of the elastic flexoelectric terms in the free energy. displacement field, and therefore it formally enters the One might wonder, at this point, whether Eq. (1) can be long-wave expansion of the dynamical matrix on the same trusted at all, given the aforementioned arbitrariness. To footing as a gradient of the polar mode: elasticity, flex- answer this question, it is useful to understand the impact oelectricity, and polarization gradient coupling all occur at of Eqs. (25) and (26) on the domain-wall solution. second order in the wave vector q. Crucially, the renormalized gradient coefficient G̃ remains unchanged, E. Covariance principle As expressed via Eqs. (22), it might appear that the three G̃0 ¼ G̃; ð27Þ dispersion coefficients C, f, and G are well-defined (and measurable) physical properties of the crystal; most authors as the respective λ-dependent contributions to f and G have indeed used such an assumption in the past, either exactly cancel out. This means that all physically measur- implicitly or explicitly. And yet, while C is indeed a well- able properties of the domain wall, i.e., the thickness ξ and defined crystal property, neither f (as pointed out in earlier the energy W, are well defined regardless of the specific ðyÞ works [21]) nor G (as we demonstrate in the following) are; convention that we choose for pκ ; the equilibrium solution on the contrary, they both suffer from an unavoidable for PðxÞ is also unaffected. The only feature that changes arbitrariness. We anticipate that their ambiguity is physical, with λ is the equilibrium solution for the elastic displace- i.e., it is not specific to the method one uses to calculate the ment field, and hence the local strain, 031002-6 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) f0 active lattice modes, i.e., those associated to the fields vν ðrÞ u0 ðxÞ ¼ − PðxÞ ¼ uðxÞ − λPðxÞ: ð28Þ C via the transformation matrix T. Step (i) does not involve any ambiguity as long as the This result looks, at first sight, surprising: by modifying the atomic displacement pattern is mesoscopic in nature (i.e., definition of the polar distortion we have obtained the same the distortion amplitudes vary on a length scale that is much profile for PðxÞ, exactly the same energy, but a different larger than the interatomic spacings). Such an assumption solution for the strain field. The solution to this puzzle implies that, if we express the distortion in reciprocal space, resides in Eq. (14), which is our gateway from the all phonon amplitudes uκα ðqÞ vanish at the zone boundary. continuum solution back to the microscopics. And indeed, This is a sufficient condition for the Fourier continuation of one can quickly verify that the change in the strain field the atomic displacements from discrete to continuum, exactly cancels out with the change in the displacements ðyÞ ð0Þ that are associated with the redefinition of pκ , leaving the ulκα ¼ uκα ðRlκ Þ → uκα ðrÞ; ð30Þ equilibrium solution for the individual atomic displace- ments well defined (that is, λ independent) and unique. This to be uniquely defined [27]. (Whenever the above condition provides, in a nutshell, the solution to the paradox of Fig. 1, breaks down one can still extract continuum fields from the and constitutes one of our main formal results. atomistics by applying standard “macroscopic averaging” That Eq. (1) behaves this way is not a coincidence, [30,31] techniques). but rather the consequence of a more general covariance Regarding step (ii), we shall define the continuum fields principle. Suppose we operate the following transformation by inverting Eq. (13), of the fields and the distortion vectors, X ðαÞ 0 ðαÞ vνβ ðrÞ ¼ uκα ðrÞT̃ κα;νβ ; ð31Þ pκ ¼ pκ þ λ; ð29aÞ κα u0 ðrÞ ¼ uðrÞ − λPðrÞ; ð29bÞ where T̃ κα;νβ is the converse transformation matrix, which satisfies the condition T̃T ¼ I. Similarly to the direct one, where λ is an arbitrary dimensionless scalar. The atomic T̃ is diagonal on the Cartesian indices, T̃ κα;νβ ¼ δαβ T̃ κα;να . distortions associated with the displacement and polariza- The columns referring to ν ¼ u, P, which we indicate as tion fields via Eq. (17) are manifestly invariant with respect ðαÞ ðαÞ to Eq. (29). This means that the original and primed ũκ ¼ T̃ κα;uα and p̃κ ¼ T̃ κα;Pα , are the duals to the direct ðβÞ ðβÞ quantities refer to the same configuration of the system; basis vectors, respectively, uκ and pκ , that we introduced i.e., they are physically equivalent macroscopic represen- in Sec. II C. The basic requirement on the direct and dual tations of the same distorted crystal structure. It is natural vectors is that they form an orthonormal set: then to require a priori from any continuum functional of u X ðαÞ ðβÞ X ðαÞ ðβÞ and P to be covariant with respect to the choice of λ, i.e., ũκ uκ ¼ p̃κ pκ ¼ δαβ ; that it transforms as Eq. (29). The results of this section κ κ demonstrate that Eq. (1) complies with such a requirement, X ðαÞ ðβÞ X ðαÞ ðβÞ provided that both flexoelectricity and polarization gradient ũκ pκ ¼ p̃κ uκ ¼ 0: ð32Þ κ κ coefficients are consistently calculated. This is a necessary condition to ensure consistency, e.g., F. Converse mapping to the continuum fields that a subsequent application of the direct and converse It is ironic, in light of these results, to realize that mapping recovers the initial values of the continuum fields. macroscopic theory is far better suited to predicting The most general choice that satisfies these constraints equilibrium atomic positions rather than “traditional” consists in introducing a set ofPsublattice-dependent macroscopic quantities, such as the strain. To rationalize weights wκ whose sum is unity, κ wκ ¼ 1. Then, we such an outcome, and get convinced that “it cannot be define the duals as otherwise,” it is illuminating to consider the converse ðβÞ mapping between continuum fields and microscopics, ũκ ¼ wκ ; ð33aÞ i.e., the procedure that allows one to extract the values of PðrÞ and uðrÞ given a distorted configuration of the ðβÞ ðβÞ crystal. We shall follow the same two-step procedure as in p̃κ ¼ wκ pκ : ð33bÞ Sec. II C, but taken in reverse order: (i) transform the discrete sublattice distortions into continuum functions of This way, the orthonormality of the elastic displacement all space, (ii) perform a local projection of the individual vectors is enforced by construction, while the remainder of ðβÞ atomic displacements onto the subspace spanned by the Eq. (32) leads to the following two conditions on pκ : 031002-7 OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022) X ðβÞ wκ pκ ¼ 0; ð34aÞ G. Arbitrariness of the strain field κ With the above derivations, we have established the continuum displacement field as a weighted average over X ðβÞ all sublattices: wκ ½pκ 2 ¼ 1: ð34bÞ κ X uðrÞ ¼ wκ uκ ðrÞ: ð37Þ κ Equation (34a) is a “hard” requirement on pκ , and must always be enforced after some choice of weights is made. The issue with this formula, which is otherwise rather This condition lifts the indeterminacy of the polar distortion trivial, is the fact that the weights are completely arbitrary. vector that we illustrated in Fig. 1, and clarifies the role of There may be, of course, some choices that are preferable wκ in subtracting the (weighted) average displacement of over others, for different reasons. Several authors [13,27], the cell from the polar mode. Doing so is consistent with for example, advocate the use of the physical masses of the physical intuition: the polarization, by its nature, is a atoms as weights; this is convenient for dynamical prob- distortion of the lattice that does not move the unit cell lems, where masses indeed play a role, and provides the of the crystal as a whole. Thus, the translational freedom physically intuitive interpretation of the displacement field that we have described in the earlier sections can be as the displacement of the local center of mass. Then, equivalently expressed as a weight freedom in the converse microscopists routinely use the positions of the heaviest mapping to the macroscopics, which provides an even more ions to define the local strain, as they correspond to the direct connection to the theory of Ref. [21]. brightest spots in the images; this implies setting their Equation (34b), on the other hand, is a consequence of weight to unity, and the others to zero. Simply taking the Eq. (33b), which is to some extent arbitrary. Indeed, one average displacement of the cell (with equal weights) is not can always multiply Eq. (34b) by a constant factor; such uncommon, either. freedom has to do with the choice of units that we use to The key point is that there is no fundamentally right (or measure the polar distortion amplitude. (For example, one wrong) choice: since spatial inversion is broken in the polar could require P0 to coincide with the spontaneous polari- structure, within the bulk domains the relation between the zation of the ferroelectric crystal, as customary in macro- cell origin and the atomic positions cannot be fixed by scopic theories.) The present convention, which consists in symmetry (see Fig. 1). Yet, the “covariance” of Eq. (1) with measuring PðxÞ in length units, has the drawback that the respect to the weight arbitrariness guarantees that the normalization condition [Eq. (34b)], and hence the values physics is uniquely described, even if the strain field of all coefficients of Eq. (1), depends on the choice of (and hence the net elastic offset across the wall, Δu) weights. Still, we prefer it here because it bears a direct depends on such choice. Note that this result, which has formal link to the eigenvectors of the dynamical matrix been established here for a static domain-wall structure, (see Appendix A), and for consistency with earlier works holds in full generality: In Appendix B we generalize it to [10,27,28]. the time-dependent regime, and use it to reconcile the Equations (34a) and (34b), together with the prescrip- existing controversies around the so-called “dynamical tions of Sec. II C, yield a well-defined procedure to flexoelectric effect.” construct the eigendisplacement vectors, and hence the An important consequence of the formalism developed model parameters, given a set of weights wκ . Starting from here is that the definitions of strain and polarization are the atomic distortions in the relaxed ferroelectric structure intimately related: they are both bound, respectively, via ðβÞ with the polarization oriented along β, dκ , we first of all Eqs. (37) and (34), to the same weight choice ambiguity. enforce Eq. (34) via And indeed, Eq. (29) shows that uðrÞ is ambiguous only in the presence of a spatially nonuniform polarization, and is ðβÞ ðβÞ X ðβÞ d̄κ ¼ dκ − wκ dκ : ð35Þ uniquely defined otherwise. To see this, recall that the strain κ is defined as the first gradient of the displacement field. If PðrÞ vanishes or is constant over space, Eq. (29) yields Then, we define the amplitude of the spontaneous dis- ∂uα =∂rβ ¼ ∂u0α =∂rβ , independent of λ; i.e., the strain tortion as becomes a well-defined quantity. This is manifestly con- sistent with Eq. (37): if all the uκ ðrÞ are equal modulo a rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X constant (which is true if the polarization is uniform), their ðβÞ P0 ¼ wκ ½d̄κ 2 : ð36Þ spatial gradients coincide; then, any choice of the weights κ yields the same result for the strain. We have achieved, therefore, a complete physical pic- Finally, we enforce Eq. (34b) by defining the dimensionless ture. There are three, at first sight unrelated, ambiguities in ðβÞ ðβÞ eigendisplacement vector as pκ ¼ d̄κ =P0 . the mapping from continuum to atomistics and vice versa, 031002-8 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) and concern (i) the definition of the local strain, (ii) the the equilibrium lattice parameter a0 of the cubic reference definition of the polar distortion, and (iii) the definition of structure by fitting the energy to the Murnaghan equation of the flexocoupling and gradient coefficients in the free state. (The resulting values are listed in Table I.) Using the energy. We have shown that (i)–(iii) share the same formal corresponding a0 × a0 × a0 five-atom cells, we then dis- root, and can be expressed as a freedom in the choice of a place the atoms along one of the main axes and reoptimize set of atomic weights wκ . This choice should be made once their positions, leading to tetragonal configurations with and for all at the beginning, and consistently respected lower energy. The displacements from the cubic positions throughout the calculation of all free-energy coefficients; dκ and the energy differences ΔE and spontaneous polari- then, the physical answers that we extract from Eq. (1) zation P (calculated via the method of Ref. [36]) are should not depend on the specific set of wκ that we use. The reported in Table I. Based on the calculated values of ΔE variational solution of the continuum differential equations and equilibrium polar distortion amplitude P0 , we then does depend, in general, on wκ , but it must be this way: if calculate A and B as we are asking, for example, “what is the local displacement of the cell at the point r,” the answer inevitably depends on 4ΔE 4ΔE A¼ ; B¼− : ð38Þ how we define such a displacement via Eq. (37). Similar ΩP20 ΩP40 considerations hold whenever we use the information on uðrÞ, extracted from experimental or theoretical domain- In order to calculate the second-order term in the long- wall structures via Eq. (37), to estimate the flexocoupling wave expansion of the force-constant matrix, Eq. (21), we coefficient by inverting Eq. (6). use the real-space supercell approach of Ref. [21]. This implies calculating the second-order moments of the III. RESULTS interatomic force constants according to A. Computational parameters ð2;xxÞ X Φκy;κ0 y ¼ Φlκy;κ0 y ðRl þ τ κ0 − τ κ Þ2x : ð39Þ Our calculations are performed in the framework l of DFT as implemented in the “in-house” LAUTREC code [32]. We use the local-density approximation [33], the [Because of inversion symmetry, displacements of the projector augmented wave method (PAW) [34], and a atoms along y do not generate electric fields along x; this plane-wave basis set with a kinetic energy cutoff of implies that the interatomic forces decay exponentially with 50 Ry in all our calculations. The PAW datasets are distance, and the lattice sums in Eq. (39) converge to a generated by treating the following orbitals as valence unique, well-defined value.] To compute Φlκy;κ0 y , we carry electrons: Bað5s2 5p6 6s2 Þ, Cað3s2 3p6 4s2 Þ, Kð3s2 3p6 4s1 Þ, out calculations in which we displace one atom at the time Nað2s2 2p6 3s1 Þ, Nbð4s2 4p6 5s2 4d3 Þ, Oð2s2 2p4 Þ, by 0.005 a.u. along y and extract the resulting forces. In Pbð6s 5d 6p Þ, Srð4s 4p 5s Þ, Tið3s 3p 4s2 3d2 Þ, and 2 10 2 2 6 2 2 6 practice, we use a 8a0 × a0 × a0 supercell, with the same Zrð4s2 4p6 5s2 4d2 Þ. In all calculations, the Brillouin zone is geometry as in Fig. 2, except that we use the centrosym- sampled using Monkhorst-Pack [35] meshes that are metric paraelectric structure as reference. (The resulting equivalent (or better) to an 8 × 8 × 8 k-point grid in the matrix elements are given in Appendix A.) five-atom bulk cell. Finally, to validate the model results against full DFT For each of six ABO3 perovskites (BaTiO3 , CaTiO3 , calculations, we prepare two domains with opposite polari- KNbO3 , NaNbO3 , PbTiO3 , and PbZrO3 ) we first calculate zation in a long (100)-oriented supercell (as in Fig. 2, for TABLE I. Data computed using DFT for six perovskite oxides: simple-cubic lattice parameter a0 (in Å); displacements (in Å) of the atoms from their high-symmetry positions in the distorted tetragonal structure at fixed lattice constants a0 , dκ , for the A cation (dA ), B cation (dB ), apical anion (dO1 ), and equatorial anions (dO2;3 ); energy of this relaxed tetragonal configuration with respect to the simple cubic one (in meV/f.u.); and polarization of the tetragonal phase (in C=m2 ). (*) indicates calculations of BaTiO3 under hydrostatic pressure; see Sec. III D. a0 dA dB dO1 dO2;3 ΔE P BaTiO3 3.935 þ0.0186 þ0.0538 −0.0359 −0.0183 −2.87 0.188 ðÞ BaTiO3 3.904 þ0.0087 þ0.0237 −0.0117 −0.0064 −0.083 0.078 CaTiO3 3.799 þ0.2369 þ0.0584 −0.0432 −0.1261 −36.72 0.508 KNbO3 3.947 þ0.0244 þ0.0553 −0.0229 −0.0287 −3.77 0.205 NaNbO3 3.907 þ0.2097 þ0.0498 −0.0480 −0.1057 −18.75 0.359 PbTiO3 3.879 þ0.1894 þ0.0793 −0.0524 −0.1081 −30.56 0.578 PbZrO3 4.102 þ0.3674 þ0.0716 −0.0241 −0.2075 −187.50 0.647 031002-9 OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022) FIG. 2. Example of unit cell used in some of our full-DFT calculations (here, N ¼ 8, and domains with opposite orientations give rise to an AO-type domain wall). different values of N), and we allow the atoms to relax along y focus our attention on BaTiO3 , PbTiO3 , and KNbO3 , and till the forces on them are negligible. We consider both AO- present the data on other materials for comparison purposes centered and BO2 -centered wall types. Note that we neglect and future reference. octahedra rotations, strain relaxations, or other wall orienta- tions that might result in energetically more favorable B. DFT calculations of domain walls structures. Our main goal here is testing the continuum We start by discussing our direct DFT calculation of approximation on a minimal Landau model of a ferroelectric the domain-wall structures. Figure 3 shows the resulting wall, and discussing the subtleties related to the treatment of atomic displacements from the high-symmetry positions. gradient effects. In this sense, incorporating additional Two features are common to all materials: the domain walls degrees of freedom to achieve a more realistic picture would are atomically thin and the atomic positions at the center of have constituted an unnecessary complication. In some each domain depend only on the material, and not on members of our materials set such a simplified model does the type of wall (AO or BO2 ). The six oxides considered not yield a physically meaningful description of the bulk or here can be roughly classified into two categories: those domain-wall structure (or both). For this reason, we primarily for which the relative displacement between A cations in AO walls BO2 walls 0.05 0 Atomic displacement parallel to wall (Å) -0.05 BaTiO3 BaTiO3 0.05 0 -0.05 KNbO3 KNbO3 0.2 0.1 0 -0.1 -0.2 PbTiO3 PbTiO3 -5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 4 5 Distance from wall (a0) Distance from wall (a0) FIG. 3. Empty symbols: relaxed displacements of the atoms in the direction parallel to the domain wall as a function of the normal coordinate, extracted from full DFT calculations (circles, squares, diamonds, and triangles represent, respectively, A, B, apical O, and equatorial O atoms). Continuous curves: atomic displacement profiles predicted by the Landau model. Two possible types of domain wall (either AO- or BO2 -centered) are shown for three of the six perovskites considered in this work. 031002-10 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) TABLE II. Domain-wall properties for six perovskite oxides: example concerns the prediction of secondary Bloch-like domain-wall energy W as computed using DFT (for the two components in PbTiO3 [40]. The contribution of secondary possible AO and BO2 domain walls) and as computed using the (or coprimary) antiferrodistortive modes (involving rota- Landau model; and domain-wall thickness ξ as computed using tions of the O6 octahedra) to the domain-wall energy and the Landau model. structure has also been studied in selected cases [41,42]. W (mJ=m2 ) ξðÅÞ These works clearly indicate that our Eq. (1) is too simplified to provide a realistic picture in many materials; DFT (AO) DFT (BO2 ) Eq. (12) Eq. (10) a follow-up work is currently under way to generalize our BaTiO3 2.95 4.59 4.43 2.204 model to more complex geometries and boundary con- ðÞ 0.321 0.322 0.324 5.441 BaTiO3 ditions. Although experimental probing of the structure and CaTiO3 63.2 60.8 70.5 2.463 energetics of ferroelectric domain walls is still challenging, KNbO3 4.46 6.61 6.03 2.303 by now the characterization methods are mature enough to NaNbO3 26.6 32.7 34.3 2.557 allow for a meaningful comparison to theoretical results. PbTiO3 54.9 58.1 63.8 2.855 For example, it is now widely accepted that domain walls PbZrO3 200 133 252 2.170 in perovskite oxides can be atomically thin (see, for example, Ref. [43] and other references therein) as we have found here. adjacent domains is significantly larger than the one between B cations (CaTiO3 , NaNbO3 , PbTiO3 , and PbZrO3 ), and those for which this is not the case (BaTiO3 and KNbO3 ). C. Landau model calculations of domain walls Such an outcome reflects the bulk distortion patterns quoted In this section we develop Landau models for each of our in Table I, which indeed shows that BaTiO3 and KNbO3 have six perovskite oxides according to the guidelines specified similar properties, e.g., regarding the small displacements of earlier. We start from the calculated distortion pattern (dκ ) their A cations. A similar classification also applies to the and energy gain (ΔE) of the tetragonal ferroelectric phase, domain-wall energies, listed in Table II: for BaTiO3 and as reported in Table I. After processing the latter values via KNbO3 the domain-wall energies are significantly smaller the procedure described in Secs. II C and II F, we readily than in other oxides (about one order of magnitude smaller obtain the values of A, B, and pκ for each bulk material. The than in PbTiO3 ), with the BO2 -type wall energy approx- gradient-mediated coefficients of the model (C, f, G) are imately 50% higher than the AO-type value. In other ð2;xxÞ then computed from pκ and the calculated Φκy;κ0 y via oxides the energies are larger, and for both types of wall they are within 25% of each other. (The BO2 type becomes Eq. (22). To illustrate our arguments, we use two different favored for CaTiO3 and PbZrO3 .) This picture is consistent choices of weights for defining pκ (and hence the A, B, f, with earlier calculations [37] in PbTiO3 and BaTiO3 , and G parameters of the model), by setting wκ either to even though in our calculations the relaxation of the cell equal values or to the physical masses of the atoms. We parameters is not allowed. shall compare the results and demonstrate their mutual The atomic configurations and energies of some per- consistency in the following. ovskite oxide domain walls have been studied in the past In Table III we summarize our results for the calculated using DFT-based methods. Padilla et al. [38] carried out a model parameters depending on the weight choice. As pioneering study on 180° domain walls in BaTiO3 using an expected, all parameters (with the exception of the elastic effective Hamiltonian built from DFT results; they reported coefficient) considerably differ between the equal-mass and that the walls are atomically thin and centered at the Ba the standard-mass convention. Recall that this difference is atoms, consistent with our results, and that the domain-wall twofold. First, there is a trivial scale factor (proportional to energies are of the order of 10 mJ=m2 . Full DFT studies of some power of P0 ) that is due to the normalization of the 180° domain walls in PbTiO3 , first by Pöykkö and Chadi polar eigendisplacements pκ . The present choice of meas- [39] and later by Meyer and Vanderbilt [37], reached the uring the polar order parameter via the norm (in length same conclusion regarding thickness, and predicted units) of the atomic distortion amplitude differs from the domain-wall energies of 100 to 200 mJ=m2 . In this case, usual convention of macroscopic theories. To connect with the most favorable domain walls are found to be centered the latter we provide in Table V (Appendix C) the same on the Pb atoms, which is again consistent with our coefficients in SI units, where we have set P0 to the findings. To the best of our knowledge, Ref. [37] was spontaneous polarization of the bulk crystal, in C=m2 . In the first to point out the geometrical offset of the atomic this case, the A and B coefficients agree between different rows between the oppositely polarized domains; we discuss choices of the weights, while a discrepancy remains in both this point extensively in the next section. f and G. Such a dependence of f and G on the weight More recent studies have revealed that surprises may be choice relates to the ambiguity in the definition of the in store even in systems that were hitherto believed to be center of mass, which is subtracted by construction from pκ simple and well understood. A particularly illuminating via Eq. (34). As we said, either source of arbitrariness must 031002-11 OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022) TABLE III. Coefficients of the Landau model of Eq. (1) 0.02 computed from DFT calculations. The values are calculated by expressing both u and P fields in Bohr units of length. Then, A is in 10−3 Ha bohr−5 ; B in 10−3 Ha bohr−7 ; f, G, and C are all in 0 atomic units of pressure, i.e., 10−3 Ha bohr−3 . In the top block Elastic displacement field u (Å) unit weights are used. In the bottom block physical atomic masses -0.02 BaTiO3 are used as weights. A B C f G 0.04 BaTiO3 −0.276 74.6 4.61 −0.560 2.47 0 ðÞ BaTiO3 −0.0503 83.3 4.77 −0.590 2.66 CaTiO3 −0.219 3.29 3.47 −0.861 2.59 -0.04 KNbO3 KNbO3 −0.322 77.8 3.27 −1.62 3.86 NaNbO3 −0.135 2.65 2.56 −1.36 2.30 0.2 PbTiO3 −0.234 4.79 3.45 −1.06 3.73 PbZrO3 −0.365 2.26 2.09 −0.386 3.14 0.1 BaTiO3 −0.448 195 4.61 −3.82 7.04 0 ðÞ BaTiO3 −0.0801 212 4.77 −3.86 7.35 -0.1 CaTiO3 −0.219 3.30 3.47 −2.28 3.87 PbTiO3 -0.2 KNbO3 −0.316 74.9 3.27 −4.15 8.26 NaNbO3 −0.215 6.72 2.56 −2.60 5.14 -5 -4 -3 -2 -1 0 1 2 3 4 5 Distance from wall (a0) PbTiO3 −0.302 8.02 3.45 −5.49 13.1 PbZrO3 −0.431 3.14 2.09 −2.75 7.24 FIG. 4. Value of the elastic displacement field as a function of distance to the domain wall, for three perovskite oxides. The lines correspond to Landau model predictions using unit weights have no impact on the physical results that we extract from (broken lines) or atomic weights (continuous lines). The symbols the model—we use this criterion to validate the internal correspond to local center-of-mass calculations computed with consistency of our theory. full DFT, for different cells centered on AO or BO2 layers, and for As a first test of such a claim, we use the calculated those two different weights. values of the coefficients to compute the energy and the width of the domain walls following Eqs. (10) and (12). We weight convention. (For an AO wall, we compute this by find that the results are indeed consistent (to machine considering the positions of the A and O atoms in a AO precision) between the two weight choices. Table II shows layer together with an average of the B and O atoms in that the domain-wall energies extracted from the model are adjacent BO2 layers—we follow an equivalent procedure also consistent with the results of our full DFT calculations; for the BO2 walls.) In all materials considered here, taking the level of agreement is remarkable considering the unit weights (discontinuous lines) leads to flatter profiles simplicity of the Landau model. The wall width, in than taking atomic weights (continuous lines), but with particular, is of the order of the cubic lattice parameter neither choice the centers of mass is fully aligned across in all cases, which provides a rather difficult test for the the wall. Any choice of weights, on the other hand, leads to continuum approximation. (The latter is expected to break the same values of atomic displacements predicted by the down at length scales that are comparable with the lattice Landau model via Eq. (13). The resulting curves are plotted periodicity.) Note that the direct DFT calculations yield a as continuous lines next to the full DFT results for the marked dependence of the domain-wall energy on the wall domain-wall structures in Fig. 3. The agreement between location with respect to the underlying atomic structure; the displacements predicted by the Landau model and those such a dependence is obviously missing within the con- from full DFT calculations is obvious. tinuum description. For a more quantitative analysis, we extract the overall Next, we compute the elastic displacement field by displacement of the center of mass across the wall Δu from combining Eqs. (5) and (9): the relaxed domain-wall structures, and use it to estimate f f x the effective flexocoupling coefficient by inverting Eq. (6): uðxÞ ¼ − PðxÞ ¼ − P0 tanh : ð40Þ C C ξ Δu f eff ¼ −C : ð41Þ The resulting profile, shown in Fig. 4, drastically differs 2P0 depending on the weight choice, but is in excellent agree- ment with the “local center of mass” that we extract from The results for both choices of weights are compared in our explicit domain-wall calculations by using the same Table IV to the values predicted by the macroscopic model. 031002-12 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) TABLE IV. Value of feff ¼ −CΔu=2P0 , in V, as computed with external parameter to bring the material closer to the phase DFT and with our Landau model. In the top block unit masses are transition. [The domain-wall thickness, also known as used. In the bottom block actual atomic masses are used. correlation length, diverges in a vicinity of a ferroelectric DFT (AO) DFT (BO2 ) Landau model phase transition, because of the vanishing A coefficient at the denominator of Eq. (10).] Temperature is the most BaTiO3 −0.203 −0.198 −0.283 obvious choice in an experimental context; this is, however, ðÞ BaTiO3 −0.279 −0.279 −0.288 impractical in the context of direct DFT simulations. A CaTiO3 −0.409 −0.454 −0.681 simpler alternative consists in applying a hydrostatic KNbO3 −0.676 −0.670 −0.793 pressure p to our simulation cells; in many ferroelectrics, NaNbO3 −1.07 −1.09 −1.33 this results in a suppression of the ferroelectric instability PbTiO3 −0.395 −0.400 −0.628 PbZrO3 0.038 −0.124 −0.373 already at moderate values of p, thus mimicking the effect of increasing temperature in real experiments. BaTiO3 −1.43 −1.43 −1.52 We apply this strategy to BaTiO3 by repeating all the BaTiO3  −1.47 −1.48 −1.49 calculations under hydrostatic pressure; the corresponding CaTiO3 −1.53 −1.57 −1.80 results are indicated in the tables with a (*) symbol. We find KNbO3 −1.93 −1.92 −2.05 NaNbO3 −1.75 −1.76 −2.01 that, at a lattice constant of 99.2% of the equilibrium a0 (see PbTiO3 −2.63 −2.64 −2.87 Table I), the A parameter is still negative, but its absolute PbZrO3 −2.04 −2.20 −2.45 value is about an order of magnitude smaller (see Table III) than in standard conditions. Consistently, as the soft mode eigenvalue approaches zero the atomic displacements One can note that the equal-weight convention yields around the wall become smoother, and the value of ξ values of f and f eff that are systematically smaller than increases to 5.441 Å. The profile of atomic displacements those obtained by setting wκ to the physical atomic masses. for these walls is shown in Fig. 5: the results of the Landau Overall, the values of f eff and f nicely agree, differing at model (continuous lines; parameters are reported in most by few tenths of a volt in all cases. Table III) are essentially in perfect agreement with the Interestingly, the continuum approach yields a slight full DFT results in a 20a0 × a0 × a0 unit cell (symbols), overestimation of the converse flexoelectric effect at the both for the AO and for the BO2 types of domain wall. The wall; such a feature appears to be systematic across all the agreement regarding domain-wall energy is also extremely materials set and irrespective of the weight convention being used. A possible explanation might lie in the BaTiO3 (under pressure) atomically sharp nature of the domain-wall structures, 0.03 Atomic disp (Å) AO DOMAIN WALL which clearly challenges the continuum description. The 0.02 B 0.01 A estimated values of f eff , however, appear to be largely O 2 O3 0 insensitive to the atomistic details of the wall (AO- and -0.01 O1 BO2 -type walls yield very similar values for most materi- -0.02 als); therefore, it is unlikely that the aforementioned -0.03 discrepancy originates from the continuum approximation 0.03 Atomic disp (Å) BO2 DOMAIN WALL itself. We suspect that such an effect may depend on higher- 0.02 0.01 order terms (either in the gradient expansion or in the polar 0 distortion amplitude) that we neglect in Eq. (1). In any case, -0.01 the accuracy of the present theory is more than sufficient for -0.02 a quantitative comparison of first-principles calculations -0.03 and experiments, e.g., along the lines of Ref. [13]. 0.01 AO-centered, atomic m Elastic disp (Å) AO-centered, unit m BO2-centered, atomic m BO2-centered, unit m D. “Soft” ferroelectric walls: BaTiO3 under pressure 0 The values of ξ in Table II are smaller than one lattice parameter in all cases, which agrees with the usual -0.01 perception that 180° ferroelectric domain walls are very -5 -4 -3 -2 -1 0 1 2 3 4 5 Distance from domain wall center (0.992a0) thin. In order to evaluate the accuracy of the polarization profile of Eq. (9) it would be desirable to study thicker FIG. 5. Atomic displacements parallel to an AO domain wall ferroelectric domain walls with a smaller distortion ampli- (top panel) and to a BO2 wall (middle panel), and elastic tude; in such a limit we expect the free-energy functional of displacement (bottom panel), as a function of the distance from Eq. (1) to match the results of direct DFT calculations the domain wall. In all cases the system studied is compressed exactly. One way to access this regime consists in using an BaTiO3 , built using supercells with a lattice parameter of 0.992a0 . 031002-13 OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022) good (see Table II): Equation (12) yields 0.324 mJ=m2 , From the point of view of the physics, the main conceptual while from full DFT calculations we obtain 0.321 mJ=m2 advance of our work can be summarized as follows. The (for the AO type) and 0.322 mJ=m2 (for the BO2 type). contribution of flexoelectricity to the domain-wall properties Note that the direct DFT calculations of the domain-wall is ill defined; however, the contribution of the polarization energies and atomic structures are numerically much more gradient energy is ill defined as well, and including both challenging than in the zero pressure case, because of the terms is essential for guaranteeing that their respective extreme softness of the ferroelectric instability. To achieve a arbitrariness cancels out. The obvious question then is, is reasonable level of accuracy, and hence allow for a mean- it really necessary to consider both terms explicitly? Or can ingful comparison to the model results, we had to be we choose the weights in such a way that the contribution of unusually careful with the usual computational parameters: flexoelectricity vanishes identically, leaving only the gra- (i) the domain-wall calculations were performed with a dient terms? A quick glance at Fig. 4 suggests, at first sight, larger 20a0 × a0 × a0 unit cell, to accommodate the thicker that the answer to the latter question be “yes”: for the equal- wall structure, (ii) the five-atom unit cell calculations were weights case, the elastic offset between the domains is done using a 20 × 8 × 8 Monhkhorst-Pack grid, exactly already small—by weighing the cations slightly less than reproducing the folded Brillouin zone of the supercell (the the oxygens, one could certainly make Δu to vanish exactly. ferroelectric distortion amplitude and double-well potential However, this would “renormalize out” flexoelectricity only depth is remarkably sensitive to the k-point mesh), (iii) the for a [100]-oriented wall; the same choice of weights would not yield a vanishing Δu for a different (say, [110]) tolerance on residual forces was set to 0.0001 eV=Å, orientation. This means that, for a truly isotropic material, approaching the inherent precision limits of the numerical flexoelectricity can be reabsorbed into the polarization algorithms. On the other hand, the calculation of the gradient energy, while anisotropic crystals generally require Landau model parameters required a similar computational its explicit treatment. In other words, one can equivalently effort as in the zero-pressure case, highlighting the obvious state that the renormalized gradient coefficients G̃ contain a advantages of our multiscale approach in softer ferroic nonanalytic elastic contribution, which generally prevents materials. their representation in a straightworward tensorial form. This adds up to the more conventional sources of nonanaliticity, IV. DISCUSSION due to the long-range electrostatic interactions, which we have not covered in the present work. Interestingly, electro- The impact of flexoelectricity on the properties of static and elastic interactions mediated by flexoelectricity ferroelectric domain walls was studied in several recent share many similarities, as the former also enjoy a covariance works [12,13,44–46]. Yudin and Tagantsev [12] established principle [27] and become analytic (i.e., short-ranged) in the role of flexoelectricity in determining the elastic offset isotropic media [27]. at the wall, as well as its impact on domain-wall energy and Another conceptual advance of this work consists in thickness via the renormalization of the polarization gra- establishing a first-principles theory of the polarization dient coefficient G. These formal results, which we have gradient coefficient G. The established procedure to infer largely built upon in our present work, were applied by the value of G from the phenomenological point of view Wang et al. [13] to domain-wall structures in PbTiO3 that consists in analyzing the dispersion of the phonon band were obtained via either first-principles calculations or associated with the soft mode [47]. The conceptual diffi- experimental microscopy images. These authors, however, culties of defining G in a microscopic context, and the often assumed that an “authentic” [13] definition of the strain counterintuitive consequences of combining LGD models exists, and overlooked (as did earlier works on this topic) with density-functional theory, have been emphasized very the inherent arbitrariness that we highlight here. The recently [48]. (For example, an estimation based on the definition used by Ref. [13] corresponds to using physical calculated spontaneous polarization and domain-wall masses as weights within our formalism. Within such a energy resulted [48] in a large variation of G as a function convention, however, our result of f 44 ¼ −2.87 V (or of the applied pressure.) In contrast with the results of 44 ¼ −2.64 V as extracted from the direct calculation f eff Ref. [48], the G coefficient as defined in our work changes of the BO2 -centered domain-wall structure) in PbTiO3 little (2%–3% deviation for an isotropic compressive strain disagrees with the value of f 44 ¼ −5.4 V quoted by of 1%) with external pressure in BaTiO3 , consistent with Wang et al. by a factor of 2. We believe that the source the usual assumptions of LGD theory. This, together with of the disagreement lies in Eq. (9b) of Wang et al. [13], the excellent accuracy of the domain-wall energies and where an additional factor of 2 is indeed present compared structures calculated within our continuum approach, to our Eq. (41). Note that our calculated C ¼ 101 GPa is in allows us to reach a more optimistic conclusion (compared excellent agreement with the value of C44 ¼ 100.8 GPa to Ref. [48]) regarding the suitability of LGD equations for quoted therein, ruling out possible issues in the definition the description of realistic domain-wall structures. of the shear elastic constant. Once the result of Ref. [13] is Before closing, it is useful to place the methodology divided by 2, it is in excellent agreement with ours. that we have developed here in the context of the existing 031002-14 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) literature. Our strategy differs in spirit from Landau- predictive power in the study of spatially inhomogeneous Ginzburg-Devonshire theories in that the latter have been structures (e.g., domain walls) in ferroics. Our formalism traditionally approached with a phenomenological mindset. demonstrates the necessity of abandoning some widespread Granted, combining first-principles techniques with LGD beliefs in continuum theory, for example, that the local theories is not new: the former are being increasingly used strain field be a physically well-defined degree of freedom [49] to estimate a subset, or even the entirety of the LGD of the crystal. While the arbitrariness in the definition of the model coefficients. However, such a practice is seldom elastic strain has a profound impact on the continuum supported by a formal link between the macroscopic and model coefficients, we demonstrate that the physical microscopic degrees of freedom, which has thwarted so far a answers derived from the model are robust against the quantitative validation of LGD models against the ab initio specific convention that is being used. results. (At the qualitative level, successful comparisons do On one hand, this results in a fundamental principle of exist; see Ref. [45] for an example that is relevant to the topics invariance that we deem of great practical utility in of this work.) The conceptual novelty of our approach validating the internal consistency of the continuum equa- consists in deriving the macroscopic field equations as a tions. On the other hand, and most importantly, our results rigorous long-wave approximation to the discrete lattice provide a stringent benchmark to determine what are the Hamiltonian. This constitutes a much more intimate synergy physically sound questions that one can ask, and what between the two levels of theory, which allows for quanti- are not. As illustrated by our practical tests, examples of tatively accurate predictions (compared to the “training” first- well-posed questions concern the domain-wall energy, or principles model) of the relevant physical properties. Most the atomic positions; we show that macroscopic theory importantly, our approach leads to a deeper awareness of the (within the validity range of the continuum approximation) internal structure of the theory, the mutual relation between can be an excellent tool to predict both—with comparable the many physical quantities involved, and their potential accuracy to the full first-principles “training model.” dependence on some arbitrary choices that are inevitable in Conversely, questions of the type “what is the contribution the mapping of the problem onto continuum fields. of flexoelectricity to the domain-wall energy?” or “what is For the above reasons, our “first-principles macro- the impact of the dynamic flexoelectric effect on the scopic theory” belongs to the class of methods that are acoustic phonon dispersion?” are physically meaningless, commonly known as second principles. The latter are as the answer can be about “anything,” depending on some obtained via an appropriate coarse graining of the first- (necessarily) arbitrary choices that one makes along principles Hamiltonian, whose physically relevant low- the way. energy degrees of freedom are treated explicitly, while most In the present work we have voluntarily chosen, for the of the original complexity is integrated out. In the case of sake of clarity, a minimal model of macroscopic phenom- ferroelectrics, the reference in this context is the “effective ena where the above ideas have a nontrivial impact. (For Hamiltonian” method, pionereed by Zhong et al. [24]; example, we have neglected most components of the strain many improvements and generalizations thereof have been tensor, as well as competing antiferrodistortive modes.) introduced over the years [50,51]. Performing a fair This choice has inevitably limited the predictive power of comparison of the virtues and shortcomings of either our study in some materials. An obvious future develop- strategy would require substantial additional work, and ment of this work consists in extending the scopes of our will be best carried out in the framework of a separate first-principles continuum approach to more complex publication. Here we limit ourselves to observe that, once structures, e.g., involving a higher dimensionality or a the continuum differential equations are discretized on a broader range of degrees of freedom. For example, it will be regular mesh corresponding to the unit cells of the original interesting to clarify whether the principle of invariance perovskite lattice, our approach essentially reduces to an established here also applies to the bichiral domain-wall effective Hamiltonian—at least for the simple domain-wall structures described in Refs. [44,45], or to the secondary geometry considered here. [A generalization to the full Bloch-like components that were theoretically predicted in three-dimensional case appears feasible, too, by writing all PbTiO3 [40]. Also, incorporating the effect of octahedral the couplings of Eq. (1) in a tensorial form and by explicitly tilts appears especially promising, in light of the results of treating the electrostatic energy.] From this point of view, Ref. [10]. Generalizing the ideas developed here in such the present method effectively bridges the gap between directions will be an exciting avenue for further study. atomistic and continuum approaches, while preserving an exact limit at length scales that are large compared to the ACKNOWLEDGMENTS interatomic spacings. O. D. acknowledges funding from the Israel Science Foundation under Grants No. 1814/14 and No. 2143/14. V. CONCLUSIONS M. S. acknowledges the support of Ministerio de Ciencia e We have established a formal mapping between con- Innovación (MICINN-Spain) through Grant No. PID2019– tinuum fields and atomistics, and demonstrated its 108573 GB-C22 and Severo Ochoa FUNFUTURE center 031002-15 OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022) of excellence (CEX2019-000917-S), of Generalitat de Eq. (A3), we find that both Eqs. (34) and (34b) are Catalunya (Grant No. 2017 SGR1506), and of the automatically satisfied by all polar modes. We can, thus, European Research Council (ERC) under the European identify the symbol hκβjναi with the transformation matrix Union’s Horizon 2020 research and innovation program introduced in Sec. II C: (Grant Agreement No. 724529). hκβjναi ¼ T κβ;να : ðA5Þ APPENDIX A: RELATION TO THE EIGENMODE REPRESENTATION In particular, the elastic displacement and polarization basis vectors are defined by In this appendix, we shall link the formalism developed in the main text to the prescriptions of Ref. [27] for the ðαÞ hκβjuαi ¼ δαβ ; hκβjTO1 αi ¼ δαβ pκ ; ðA6Þ construction of the continuum Hamiltonian. Their proposed strategy consists in identifying the distortion patterns where we have assumed that TO1 corresponds to the associated with the mechanical displacement and polari- ferroelectric “soft mode.” zation fields with, respectively, the acoustic and “soft” There is a slight drawback with such a procedure: polar eigenmodes of the zone-center dynamical matrix D. different conventions for the weights lead to definitions We adopt here a slightly more general definition by ðαÞ introducing the following operator, of pκ that are generally not related via Eq. (25). Indeed, ðαÞ the configuration space spanned by pκ changes depending 1 ð0Þ on the weights, as the three polar optical modes of the D̃κα;κ0 β ¼ pffiffiffiffiffiffiffiffiffiffiffi Φκα;κ0 β ; ðA1Þ wκ wκ0 perovskite structure can mix. (This is a well-known issue in the construction of effective Hamiltonian for ferroelectrics, which reduces to the standard definition of D when the where typically the lowest eigenvector of the force-constant weights are set to the physical atomic masses mκ . (In our ðαÞ matrix is used for pκ ; this corresponds to choosing equal formalism the weights are dimensionless, so the appropriate weights in the context of our formalism.) In practical cases, choice in this case consists in setting them to the fractional ðαÞ pκ might not reproduce the correct distortion pattern (and mass ofPthe sublattice; then, D̃ reduces to M tot D, where energetics) of the bulk ferroelectric ground state, which is Mtot ¼ κ mκ is the total mass of the cell.) Following undesirable in the study of a domain wall. To avoid this Ref. [27], we define the eigendisplacement pattern jναi issue, in this work we have followed the prescriptions of associated to a given normal mode να (ν runs over the four ðαÞ T 1u irreps of the cubic perovskite structure, including Sec. II F and defined pκ starting from the relaxed acoustic and optical TO1 –TO3 modes; α is a Cartesian distortion pattern of the bulk ferroelectric crystal instead. direction) as In the limit of a weak ferroelectric instability, one can show ðαÞ that this definition of pκ exactly matches the eigenvector sffiffiffiffiffiffi representation provided by Eq. (A6) regardless of the 1 ðνÞ hκβjναi ¼ hκβjvα i; ðA2Þ choice of the weights. wκ ðνÞ APPENDIX B: DYNAMICAL EQUATIONS where jvα i are the normalized eigenvectors of D̃. (The OF MOTION present use of the bra-ket notation follows the conventions established in earlier works [27,28]. Both bras and kets are Tagantsev [52] and Kvasov and Tagantsev [53] claimed real vectors in a space of dimension 3N at , where N at is the that there are two well-defined contributions to the bulk total number of atoms in the unit cell; hκβj form a complete flexoelectric tensor, static and dynamic in nature, and that orthonormal basis.) This definition results in a generalized they are, in principle, separately measurable. Later works orthonormality condition for the eigendisplacements, [22] clarified that such a partition is arbitrary, and that (i) the total flexoelectric coefficient is meaningful for hναjWjν0 βi ¼ δνν0 δαβ ; ðA3Þ dynamical problems and (ii) either the total or the “static” flexoelectric tensor yield identical answers at mechanical where the “overlap operator” W is defined as equilibrium. To firm up our arguments, we shall revisit this long-standing debate on the static versus dynamic contri- W κα;κ0 β ¼ wκ δκκ0 δαβ : ðA4Þ bution to the flexoelectric tensor in light of the results presented so far. We shall see that it bears strong con- One can show [27] that, within such prescriptions, the nections to the aforementioned ambiguities in the defini- eigendisplacements associated to the acoustic mode tions of the continuum fields, and that our formalism (ν ¼ u) are independent of the weights, thus recovering resolves once and for all the existing confusion around Eq. (16). Then, by setting ν ¼ u and ν0 ¼ TO1 –TO3 in this topic. The coupled dynamics of the polar and acoustic 031002-16 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) degrees of freedom is, of course, irrelevant to the study of free-energy functional coefficients. If we do so, the matrix static structures, such as the domain walls that we consider element MuP then reduces to their definition of the in this work. Still, it is interesting to discuss this topic here, “dynamic” contribution. If we made a different choice, as it provides an additional proof of the internal consistency the partition between the two would change arbitrarily— of our arguments. and yet both the dynamical (phonon frequencies and dis- To describe the dynamical evolution of the mechanical persions) and static properties (domain-wall energy, equi- and acoustic degrees of freedom, we need to work out the librium atomic positions) predicted by our Lagrangian kinetic energy density in terms of the mode velocities. We would be exactly the same. It is interesting to consider shall write it as the special case where the weights are set to the physical masses of the atoms P divided by the total mass of the 1X cell, wκ ¼ mκ =ð κ0 mκ0 Þ. The orthogonality condition, T ðrÞ ¼ v_ ðrÞM να;ν0 β v_ ν0 β ðrÞ; ðB1Þ Eq. (A3), immediately leads then to MuP ¼ 0; i.e., the 2 νν0 αβ να “dynamical flexoelectric effect” disappears altogether. Given that the magnitude, and even the very existence, where αβ run over all the Cartesian components of the of such an effect depends on some arbitrary convention we vector fields indexed by νν0, and the matrix Mνα;ν0 β , of the have made along the way in order to map our lattice- dimension of a mass density, is the normal mode repre- dynamical problem onto a continuum Lagrangian density, sentation of the “mass operator” M: we are forced to conclude that such an effect is not measurable. Still, we find that incorporating a mass cross 1 term, as suggested in Ref. [12], is necessary to guarantee Mνα;ν0 β ¼hναjMjν0 βi; Ω that the physical predictions of the theory are unaffected by hκαjMjκ0 βi ¼ mκ δκκ0 δαβ : ðB2Þ such ambiguities. (As usual, κκ0 are sublattice indices; mκ are atomic masses.) APPENDIX C: SUPPORTING NUMERICAL DATA The off-diagonal kinetic term, coupling the acoustic and In Table V we provide the complete list of the calculated optical mode velocities, has been identified as a dynamical model coefficients for all materials. This is essentially the contribution to the bulk flexoelectric tensor by Tagantsev same data as in Table III, only expressed in SI units while and co-workers [12,52,53]. setting P0 to the spontaneous ferroelectric polarization of Based on the arguments of the earlier sections, it is clear the bulk crystal. This conversion is useful for two purposes. that the magnitude of such contribution, and even whether First, it shows that the A and B coefficients are consistent it exists at all, depends on the choice we make for the between different weight choices, provided that the respec- weights wκ . The definition given by Tagantsev of the tive distortion vectors pκ are related via Eq. (25); i.e., they “static” flexoelectric tensor corresponds, within our for- only differ by a sublattice-independent constant. [This is malism, to using equal weights in the construction of our obviously the case if the electrical polarization is chosen as TABLE V. Coefficients (in SI units) of the Landau model of Eq. (1) computed from DFT calculations, when P0 is the spontaneous polarization of the bulk crystal (measured in C=m2 ). In the top block unit weights are used. In the bottom block physical atomic masses are used. A (109 m=F) B (1010 Ω2 m3 =kg) C (1011 Pa) f (V) G (10−10 m3 =F) G̃ (10−10 m3 =F) BaTiO3 −0.857 2.43 1.36 −0.283 0.214 0.208 ðÞ BaTiO3 −0.145 2.36 1.40 −0.288 0.221 0.215 CaTiO3 −1.66 0.645 1.02 −0.681 0.550 0.505 KNbO3 −0.934 2.22 0.961 −0.793 0.313 0.248 NaNbO3 −1.56 1.21 0.755 −1.33 0.746 0.510 PbTiO3 −1.00 0.299 1.01 −0.628 0.447 0.409 PbZrO3 −4.15 0.991 0.616 −0.373 1.00 0.979 BaTiO3 −0.857 2.43 1.36 −1.52 0.378 0.208 ðÞ BaTiO3 −0.145 2.36 1.40 −1.49 0.373 0.215 CaTiO3 −1.66 0.645 1.02 −1.80 0.822 0.505 KNbO3 −0.934 2.22 0.961 −2.05 0.684 0.248 NaNbO3 −1.56 1.21 0.755 −2.01 1.05 0.510 PbTiO3 −1.00 0.299 1.01 −2.87 1.22 0.409 PbZrO3 −4.15 0.991 0.616 −2.45 1.95 0.979 031002-17 OSWALDO DIÉGUEZ and MASSIMILIANO STENGEL PHYS. REV. X 12, 031002 (2022) TABLE VI. Elements of the second-moment interatomic force- [1] J. Seidel, Topological Structures in Ferroic Materials constant matrix (in atomic units), as computed using DFT (Springer, Cham, 2016), p. 826. calculations in 40-atom 8a0 × a0 × a0 unit cells, for six per- [2] Shanquan Chen, Shuai Yuan, Zhipeng Hou, Yunlong Tang, ovskite oxides. Jinping Zhang, Tao Wang, Kang Li, Weiwei Zhao, Xingjun Liu, Lang Chen, Lane W. Martin, and Zuhuang Chen, Φð2;xxÞ Recent Progress on Topological Structures in Ferroic Thin BaTiO3 0.0920 −0.5561 −0.2022 −0.0383 −0.1299 Films and Heterostructures, Adv. Mater. 33, 2000857 −0.5561 −0.0191 0.0539 −0.0259 −0.0202 (2021). −0.2022 0.0539 −0.0847 −0.7573 0.0142 [3] P. Zubko, G. Catalan, and A. K. Tagantsev, Flexoelectric −0.0383 −0.0258 −0.7573 −0.4640 0.1104 Effect in Solids, Annu. Rev. Mater. Res. 43, 387 (2013). −0.1299 −0.0202 0.0142 0.1104 −0.2174 [4] G. Catalan, A. Lubk, A. H. G. Vlooswijk, E. Snoeck, C. Magen, A. Janssens, G. Rispens, G. Rijnders, D. H. A. BaTiO3 ðÞ 0.0915 −0.5773 −0.1854 −0.0319 −0.1605 Blank, and B. Noheda, Flexoelectric Rotation of Polariza- −0.5773 −0.0230 0.0580 0.0177 −0.0231 tion in Ferroelectric Thin Films, Nat. Mater. 10, 963 (2011). −0.1854 0.0580 −0.0799 −0.7699 0.0142 [5] H. Lu, C.-W. Bark, D. Esque de los Ojos, J. Alcala, C. B. −0.0319 0.0178 −0.7699 −0.4921 0.1119 Eom, G. Catalan, and A. Gruverman, Mechanical Writing of −0.1605 −0.0231 0.0142 0.1119 −0.2326 Ferroelectric Polarization, Science 336, 59 (2012). CaTiO3 0.0082 −0.3563 −0.1513 −0.0596 0.1967 [6] L. E. Cross, Flexoelectric Effects: Charge Separation in −0.3563 −0.0374 0.1251 0.2363 −0.0124 Insulating Solids Subjected to Elastic Strain Gradients, −0.1513 0.1252 −0.1968 −0.9531 0.0092 J. Mater. Sci. 41, 53 (2006). −0.0596 0.2363 −0.9531 −0.7121 0.1835 [7] J. D. Axe, J. Harada, and G. Shirane, Anomalous Acoustic 0.1967 −0.0124 0.0092 0.1835 −0.0636 Dispersion in Centrosymmetric Crystals with Soft Optic KNbO3 0.0195 −0.3345 −0.1154 −0.0194 −0.0530 Phonons, Phys. Rev. B 1, 1227 (1970). −0.3345 −0.0671 0.1619 0.2763 −0.0216 [8] Henning Pttker and Ekhard K. H. Salje, Flexoelectricity, −0.1154 0.1619 −0.1782 −0.9074 0.0163 Incommensurate Phases and the Lifshitz Point, J. Phys. −0.0194 0.2763 −0.9075 −0.7123 0.1230 Condens. Matter 28, 075902 (2016). −0.0530 −0.0216 0.0163 0.1230 −0.0259 [9] A. K. Tagantsev, K. Vaideeswaran, S. B. Vakhrushev, A. V. Filimonov, R. G. Burkovsky, A. Shaganov, D. Andronikova, NaNbO3 0.0022 −0.2023 −0.0782 −0.0274 0.0916 A. I. Rudskoy, A. Q. R. Baron, H. Uchiyama, D. Chernyshov, −0.2023 −0.0821 0.2043 0.3915 −0.0237 A. Bosak, Z. Ujma, K. Roleder, A. Majchrowski, J. H. Ko, −0.0782 0.2043 −0.2783 −0.9913 0.0250 and N. Setter, The Origin of Antiferroelectricity in PbZrO3 , −0.0274 0.3916 −0.9913 −0.8334 0.1769 Nat. Commun. 4, 2229 (2013). 0.0916 −0.0237 0.0250 0.1769 −0.0063 [10] Andrea Schiaffino and Massimiliano Stengel, Macroscopic PbTiO3 −0.0529 −0.3944 −0.0856 −0.0003 0.1969 Polarization from Antiferrodistortive Cycloids in Ferroe- −0.3944 −0.0254 0.0977 0.1092 −0.0240 lastic SrTiO3 , Phys. Rev. Lett. 119, 137601 (2017). −0.0856 0.0978 −0.1416 −0.7698 −0.0300 [11] Ekhard K. H. Salje, Suzhi Li, Massimiliano Stengel, Peter −0.0003 0.1092 −0.7698 −0.4886 0.0226 Gumbsch, and Xiangdong Ding, Flexoelectricity and the 0.1969 −0.0240 −0.0300 0.0226 −0.2507 Polarity of Complex Ferroelastic Twin Patterns, Phys. Rev. B 94, 024114 (2016). PbZrO3 −0.1344 −0.3815 −0.1047 0.0158 0.3713 [12] P. V. Yudin and A. K. Tagantsev, Fundamentals of Flexoe- −0.3815 0.0260 0.0573 0.1440 −0.0367 lectricity in Solids, Nanotechnology 24, 432001 (2013). −0.1047 0.0573 −0.1189 −0.5866 −0.0134 [13] Y. J. Wang, Y. L. Tang, Y. L. Zhu, Y. P. Feng, and X. L. Ma, 0.0158 0.1440 −0.5866 −0.5552 0.0543 Converse Flexoelectricity around Ferroelectric Domain 0.3713 −0.0367 −0.0134 0.0543 −0.2071 Walls, Acta Mater. 191, 158 (2020). [14] R. Resta, Towards a Bulk Theory of Flexoelectricity, Phys. Rev. Lett. 105, 127601 (2010). a measure of the atomic distortion, but not when the total [15] J. Hong, G. Catalan, J. F. Scott, and E. Artacho, The norm of the distortion is used; in the latter case there is Flexoelectricity of Barium and Strontium Titanates from generally an overall scaling factor that originates from First Principles, J. Phys. Condens. Matter 22, 112201 Eq. (34b).] Second, the coefficients are now expressed in (2010). [16] Miquel Royo and Massimiliano Stengel, First-Principles the same units as in conventional macroscopic theory, Theory of Spatial Dispersion: Dynamical Quadrupoles and allowing for a direct comparison. In the case of BaTiO3 , for Flexoelectricity, Phys. Rev. X 9, 021050 (2019). example, our calculated polarization gradient coefficient G [17] Miquel Royo and Massimiliano Stengel, Lattice-Mediated is in good agreement with the phenomenological value of Bulk Flexoelectricity from First Principles, Phys. Rev. B G ¼ 0.2 × 10−10 m3 =F reported in Ref. [45]. 105, 064101 (2022). The elements of the second-moment interatomic force- [18] X. Gonze et al., ABINIT: First-Principles Approach to constant matrix, calculated via Eq. (39), are reported in Material and Nanosystem Properties, Comput. Phys. Com- Table VI for reference. mun. 180, 2582 (2009). 031002-18 TRANSLATIONAL COVARIANCE OF FLEXOELECTRICITY AT … PHYS. REV. X 12, 031002 (2022) [19] Aldo H. Romero et al., ABINIT: Overview and Focus [38] J. Padilla, W. Zhong, and D. Vanderbilt, First-Principles on Selected Capabilities, J. Chem. Phys. 152, 124102 Investigation of 180° Domain Walls in BaTiO3 , Phys. Rev. (2020). B 53, R5969 (1996). [20] Q. Li, C. T. Nelson, S. L. Hsu, A. R. Damodaran, L. L. Li, [39] S. Pykk and D. J. Chadi, Ab Initio Study of 180° Domain A. K. Yadav, M. McCarter, L. W. Martin, R. Ramesh, and Wall Energy and Structure in PbTiO3 , Appl. Phys. Lett. 75, S. V. Kalinin, Quantification of Flexoelectricity in 2830 (1999). PbTiO3 =SrTiO3 Superlattice Polar Vortices Using Machine [40] Jacek C. Wojdeł and Jorge Íñiguez, Ferroelectric Transi- Learning and Phase-Field Modeling, Nat. Commun. 8, tions at Ferroelectric Domain Walls Found from First 1468 (2017). Principles, Phys. Rev. Lett. 112, 247603 (2014). [21] J. Hong and D. Vanderbilt, First-Principles Theory and [41] Axel Lubk, S. Gemming, and N. Spaldin, First-Principles Calculation of Flexoelectricity, Phys. Rev. B 88, 174107 Study of Ferroelectric Domain Walls in Multiferroic Bis- (2013). muth Ferrite, Phys. Rev. B 80, 104110 (2009). [22] M. Stengel, Flexoelectricity from Density-Functional Per- [42] Oswaldo Diguez, Pablo Aguado-Puente, Javier Junquera, turbation Theory, Phys. Rev. B 88, 174106 (2013). and Jorge Iniguez, Domain Walls in a Perovskite Oxide [23] W. Zhong, David Vanderbilt, and K. M. Rabe, Phase with Two Primary Structural Order Parameters: First- Transitions in BaTiO3 from First Principles, Phys. Rev. Principles Study of BiFeO3 , Phys. Rev. B 87, 024102 Lett. 73, 1861 (1994). (2013). [24] W. Zhong, David Vanderbilt, and K. M. Rabe, First- [43] Donald M. Evans, Vincent Garcia, Dennis Meier, and Principles Theory of Ferroelectric Phase Transitions for Manuel Bibes, Domains and Domain Walls in Multifer- Perovskites: The Case of BaTiO3 , Phys. Rev. B 52, 6301 roics, Phys. Sci. Rev. 5, 20190067 (2020). (1995). [44] P. V. Yudin, A. K. Tagantsev, E. A. Eliseev, A. N. [25] Ph. Ghosez and K. M. Rabe, Microscopic Model of Ferro- Morozovska, and N. Setter, Bichiral Structure of Ferro- electricity in Stress-Free PbTiO3 Ultrathin Films, Appl. electric Domain Walls Driven by Flexoelectricity, Phys. Phys. Lett. 76, 2767 (2000). Rev. B 86, 134102 (2012). [26] Igor A. Kornev, S. Lisenkov, R. Haumont, B. Dkhil, and L. [45] Yijia Gu, Menglei Li, Anna N. Morozovska, Yi Wang, Bellaiche, Finite-Temperature Properties of Multiferroic Eugene A. Eliseev, Venkatraman Gopalan, and Long-Qing Chen, Flexoelectricity and Ferroelectric Domain Wall BiFeO3 , Phys. Rev. Lett. 99, 227602 (2007). Structures: Phase-Field Modeling and DFT Calculations, [27] Massimiliano Stengel, Unified Ab Initio Formulation of Phys. Rev. B 89, 174111 (2014). Flexoelectricity and Strain-Gradient Elasticity, Phys. Rev. [46] Bo Wang, Yijia Gu, Shujun Zhang, and Long-Qing Chen, B 93, 245107 (2016). Flexoelectricity in Solids: Progress, Challenges, and Per- [28] Asier Zabalo and Massimiliano Stengel, Switching a Polar spectives, Prog. Mater. Sci. 106, 100570 (2019). Metal via Strain Gradients, Phys. Rev. Lett. 126, 127601 [47] J. Hlinka and P. Márton, Phenomenological Model of a 90° (2021). Domain Wall in BaTiO3 -Type Ferroelectrics, Phys. Rev. B [29] Max Born and Kun Huang, Dynamical Theory of Crystal 74, 104104 (2006). Lattices (Oxford University Press, Oxford, 1954). [48] Atanu Samanta, Suhas Yadav, Zongquan Gu, Cedric J. G. [30] A. Baldereschi, S. Baroni, and R. Resta, Band Offsets in Meyers, Liyan Wu, Dongfang Chen, Shishir Pandya, Robert Lattice-Matched Heterojunctions: A Model and First- A. York, Lane W. Martin, Jonathan E. Spanier, and Ilya Principles Calculations for GaAs/AlAs, Phys. Rev. Lett. Grinberg, A Predictive Theory for Domain Walls in Oxide 61, 734 (1988). Ferroelectrics Based on Interatomic Interactions and Its [31] J. Junquera, M. H. Cohen, and K. M. Rabe, Nanoscale Implications for Collective Material Properties, Adv. Mater. Smoothing and the Analyis of Interfacial Charge and 34, 2106021 (2022). Dipolar Densities, J. Phys. Condens. Matter 19, 213203 [49] Sergey Artyukhin, Kris T. Delaney, Nicola A. Spaldin, and (2007). Maxim Mostovoy, Landau Theory of Topological Defects in [32] LAUTREC is a massively parallel “in-house” electronic- Multiferroic Hexagonal Manganites, Nat. Mater. 13, 42 structure code currently maintained by M. Stengel. (2014). [33] J. P. Perdew and Y. Wang, Accurate and Simple Analytic [50] Jacek C. Wojdeł, Patrick Hermet, Mathias P. Ljungberg, Representation of the Electron-Gas Correlation Energy, Philippe Ghosez, and Jorge Íñiguez, First-Principles Model Phys. Rev. B 45, 13244 (1992). Potentials for Lattice-Dynamical Studies: General Meth- [34] P. E. Blöchl, Projector Augented-Wave Method, Phys. Rev. odology and Example of Application to Ferroic Perovskite B 50, 17953 (1994). Oxides, J. Phys. Condens. Matter 25, 305401 (2013). [35] H. J. Monkhorst and J. D. Pack, Special Points for Brillouin- [51] Philippe Ghosez and Javier Junquera, Modeling of Ferro- Zone Integrations, Phys. Rev. B 13, 5188 (1976). electric Oxide Perovskites: From First to Second Principles, [36] M. Stengel and N. A. Spaldin, Accurate Polarization within Annu. Rev. Condens. Matter Phys. 13, 325 (2022). a Unified Wannier Function Formalism, Phys. Rev. B 73, [52] A. K. Tagantsev, Piezoelectricity and Flexoelectricity in 075121 (2006). Crystalline Dielectrics, Phys. Rev. B 34, 5883 (1986). [37] B. Meyer and David Vanderbilt, Ab Initio Study of Ferro- [53] Alexander Kvasov and Alexander K. Tagantsev, Dynamic electric Domain Walls in PbTiO3 , Phys. Rev. B 65, 104111 Flexoelectric Effect in Perovskites from First-Principles (2002). Calculations, Phys. Rev. B 92, 054104 (2015). 031002-19