Switching a polar metal via strain gradients
SSwitching a polar metal via strain gradients
Asier Zabalo and Massimiliano Stengel
1, 2 Institut de Ci`encia de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain ICREA-Instituci´o Catalana de Recerca i Estudis Avanats, 08010 Barcelona, Spain (Dated: July 24, 2020)Although rare, spontaneous breakdown of inversion symmetry sometimes occurs in a materialwhich is metallic: these are commonly known as polar metals or ferroelectric metals. Their polar-ization , however, cannot be switched via an electric field, which limits the experimental control overband topology. Here we shall investigate, via first-principles theory, flexoelectricity as a possibleway around this obstacle with the well known polar metal LiOsO . The flexocoupling coefficientsare computed for this metal with high accuracy with a completely new approach based on real-space sums of the inter-atomic force constants. A Landau-Ginzburg-Devonshire-type first-principlesHamiltonian is built and a critical bending radius to switch the material is estimated, whose orderof magnitude is comparable to that of BaTiO . The so-called polar or ferroelectric metals [1], first pro-posed by Anderson more than half a century ago in thecontext of martensitic transformations [2], have been at-tracting increasing attention recently. The prototypi-cal (and historically the first experimentally known) ma-terial realization is lithium osmate, which undergoes aferroelectric-like transition at 140 K from the centrosym-metric R ¯3 c to the non-centrosymmetric R c space group[3]. Since its discovery, the list of known polar metalshas been steadily growing [4]. Their interest lies on theunusual physics that may emerge from the coexistenceof metallicity and polarity, two properties that were ini-tially regarded as contraindicated. For instance, theyprovide excellent opportunities to study exotic quantumphenomena, like non-centrosymmetric superconductivity[5, 6] or spin-polarized currents [7]. In spite of consid-erable progress, however, a long-standing issue still re-mains, and concerns the ability to control polarity via anappropriate external field. Indeed, due to the presenceof free carriers in the bulk the most obvious means ofswitching polarity in ferroelectrics, i.e. an external elec-tric field, is ruled out. Such a control would help shedsome light on their fundamental physics, and possibly de-vise some applications, e.g., in nanoscale electronic andthermoelectric devices [8]. Our goal is to demonstratethat flexoelectricity can solve this issue.Flexoelectricity describes the coupling between a straingradient and the macroscopic polarization and, unlike itshomogeneous counterpart (piezoelectricity), it does notrequire any particular space group to be present [9–11].While flexoelectricity is hardly a new discovery [12], itspractical relevance was demonstrated only recently, thusreviving this field from both the experimental [13–15] andtheoretical [16] points of view. Of course, the electricalpolarization can only be defined in insulating crystals, sothe macroscopic flexoelectric coefficient of a polar metalvanishes identically. Hovever, as we shall demonstrateshortly, the coupling between polar lattice modes anda strain gradient does exist even in metals. Since elas-tic fields, unlike electric fields, are not screened by free carriers, this constitutes, in principle, a viable means ofcontrolling polarity. Still, whether the relevant couplingsare strong enough for such a mechanism to be experimen-tally accessible, is currently unknown. First-principlescalculations could be very helpful in this context, andindeed electronic-structure methods to study flexoelec-tricity have seen an impressive progress in recent years[16, 17]. However, the calculation of the flexocouplingcoefficients remains a subtle task even with insulatorsand treating metallic crystals falls outside the present ca-pabilities of the density functional theory (DFT) basedcodes.Here we overcome such limitations by developing anaccurate and general method to calculate flexocouplingcoefficients in metals, which is based on real space sumsof the inter-atomic force constants (IFC-s). We demon-strate our computational strategy by calculating the flex-ocoupling coefficients for LiOsO as a test case and wecompare them with the ones of BaTiO , probably oneof the most studied ferroelectric materials. Finally, weuse the aforementioned values, in combination with afirst-principles based effective Hamiltonian that we haveconstructed by expanding the energy around the cen-trosymmetric cubic phase, to estimate the critical bend-ing radius of LiOsO . We find the values in line withthose calculated for BaTiO , a material where flexoelec-tric switching of the polar domains has been experimen-tally demonstrated already. Based on these results, me-chanical switching of LiOsO mediated by flexoelectricityappears well within experimental reach.To start with, we shall consider a setup as illustratedin Fig. 1, i.e. of a LiOsO (or n -doped BaTiO ) samplethat is cut along some crystallographic direction ˆ q andmechanically bent via some external load. Within the in-terior of the film, the polar order parameter is assumed tobe homogeneous, and its amplitude is described by somethree-dimensional vector u with the physical dimensionof length. In the following we shall quantify the coercivebending radius , i.e. the radius of curvature that needs tobe applied in order to switch the polar order parameter a r X i v : . [ c ond - m a t . m t r l - s c i ] J u l FIG. 1. A bending type strain gradient is applied to a macro-scopic crystal along the direction ˆ q . The external strain gra-dient couples to the polar modes resulting in a displacementof the atoms and, as a consequence, the structure evolves toanother symmetrically equivalent ferroelectric state. between two neighboring local minima, which are degen-erate at mechanical equilibrium. We shall calculate thecritical radius via the following formula, R crit = f eff F coerc , (1)where f eff is the effective flexocoupling coefficient asso-ciated with the flexural deformation, and F coerc is theminimal generalized force that is required for the mode u to cross the energy barrier between two minima. Thus,the problem can be divided into two separate tasks: (i)determining the coupling between a flexural deformationand the polar mode, described by f eff , as a function of thecrystallographic orientation, and (ii) identifying the mostlikely switching paths and the corresponding energetics.From now on, we shall assume a Landau-like expansionof the energy around the high-symmetry cubic structureas a function of the relevant parameters, following the es-tablished common practice in theoretical studies of per-ovskite ferroelectrics. In this context, task (ii) entails noconceptual difficulties, as it consists in mapping the po-tential energy surface of the crystal as a function of therelevant lattice degrees of freedom – such a procedure hasbeen successfully carried out for a wide range of mate-rials already. The main technical obstacle resides in (i),since no established methods exist for the calculation of f eff in metals. Given the novelty, we shall focus on thispoint in the following.In full generality, the flexocoupling tensor in a “soft- mode” material can be defined as follows [18], f I αβ,γλ = − M h P α | ∂ D ( q ) ∂q γ ∂q λ (cid:12)(cid:12)(cid:12)(cid:12) q =0 | A β i (2)where D ( q ) is the dynamical matrix of the crystal at acertain wavevector q , and the bra and kets represent theTO and acoustic eigenvectors at the Γ ( q = 0) point ofthe Brillouin zone. This tensor describes the force on u α that is produced by a macroscopic strain gradient η β,γλ ,the latter expressed in “type-I” form (hence the super-script “I”), i.e., as the second gradient of the displace-ment field. For our present scopes, it is more convenientto work in type-II form, which can be recovered via thefollowing transformation, f II αλ,βγ = f I αβ,γλ + f Iαγ,λβ − f I αλ,βγ . (3)The main technical challenge from a computational pointof view consists in taking the second gradient with re-spect to q of the dynamical matrix, D ( q ). In insulators,this task is already delicate at the formal level, since D ( q )has a nonanalytic behavior in vicinity of Γ; this requiresa careful treatment of the macroscopic electric fields be-fore performing the perturbative long-wave expansion. Inmetals at finite temperature things appear simpler con-ceptually, since the adiabatic dynamical matrix D ( q ) isan analytic function of the wavevector q over the wholeBrillouin zone. This means that the second q -gradientappearing in Eq. (2) is always well defined without tak-ing any further precaution. However, this methodologyhas not been generalized to metals yet. To circumventthis obstacle, we define the long-wave expansion of thedynamical matrix as the real-space moments of the inter-atomic force constants, following the method describedin [11, 18]. In particular, the IFC’s are first defined asthe second derivative of the total energy with respect toatomic displacements,Φ lκα,κ β = ∂ E∂u κα ∂u lκ β . (4)Then, we write ∂ D ( q ) κα,κ β ∂q γ ∂q λ (cid:12)(cid:12)(cid:12)(cid:12) q =0 = X l Φ lκα,κ β √ m κ m κ ( d lκκ ) γ ( d lκκ ) λ , (5)where m κ is the mass of atom κ , d lκκ = R l + τ κ − τ κ , R l is the Bravais lattice vector indicating the location of the l -th cell and τ κ is the position of atom κ within the unitcell l = 0. The short-range nature of the interatomicforces guarantees that the lattice sums of Eq. (5) willeventually converge to the correct physical value whena dense enough q -point mesh is used to calculate thereal-space force constants of Eq. (4). Interestingly, thecoupling between two acoustic modes is directly related TABLE I. Independent components of the elastic (in GPa)and flexocoupling (in eV) tensor with a 16 × × q points for the IFC-s. The n-type flexocoupling coefficientsof BaTiO are shown here. C C C f II11 f II12 f II LiOsO Long-wave 364.7 129.5 44.3 − Long-wave 346.1 121.7 134.5 − − (in a crystal that is free of stresses) to the elastic tensorcomponents via [19] C αγ,βλ + C αλ,βγ − M h A α | ∂ D ( q ) ∂q γ ∂q λ (cid:12)(cid:12)(cid:12)(cid:12) q =0 | A β i . (6)This is a useful consistency check: one can then comparethe results with a more conventional calculation of theelastic tensor [20] to gauge the reliability of the flexocou-pling coefficients as determined via Eq. (2). Note that theelastic tensor components are themselves a crucial ingre-dient for calculating the effective flexocoupling of Eq. (1)starting from the flexocoupling tensor f II ; therefore, it isimportant to ensure that the two physical quantities arecalculated with consistent accuracy.Our first principles calculations are performed withthe open-source abinit [21, 22] package. (Details of thecomputational parameters are provided in the Supple-mentary Material.) Numerical results for both BaTiO and LiOsO are shown in Table I. Clearly, the largestflexocouplings are f in LiOsO and f for BaTiO .(The latter material behaves very similarly to SrTiO [18], which is natural to expect given the affinities in theelectronic and atomic structure.) Their absolute valuesare similar overall, which provides a first indication thatthe flexocoupling is comparably strong in these two mate-rials. Note that the discrepancy in the elastic constantscalculated via the two different methods is less than a1 % for the three independent components of LiOsO ,which confirms the excellent quality of the calculations.We also show in Fig. S3 (and Table S4) the convergenceof the numerical results for both the elastic and flexocou-pling constants as a function of the q -point mesh, furthercorroborating this point.To make further progress, we use the value of Table I tocompute the effective flexocoupling coefficients for threerepresentative orientations of the sample ([100], [110] and[111]), either in the beam-bending or the plate-bendinglimit. (We focus on the beam-bending limit following thedefinitions of Ref. [23]; explicit formulas are reported inthe Supplementary Material.) The results, shown in Ta-ble II, indicate that [100] is by far the bending directionthat produces the largest flexocoupling in LiOsO . Thesituation in BaTiO seems to be more balanced overall, TABLE II. Effective flexocoupling coefficients (absolute val-ues) for 100, 110 and 111 oriented samples (in eV units) forthe beam-bending limit. f f f LiOsO with a slight preference for [110] and [111] directions over[100]. Note, however, that for each surface orientationˆ q , the effective flexocoupling describes the flexo-inducedforce acting on the polar mode along ˆ q . Depending onthe switching path, such force might not be parallel tothe direction along which the polar mode evolves dur-ing switching, ˆ s ; in such cases the effective flexocouplingneeds to be scaled by the projection ˆ q · ˆ s . Since therelevant paths in BaTiO (see next paragraph) involve[100]-oriented switching, such geometrical factor reducesthe [110] and [111] coefficients by √ √ f eff to a similar magnitude.Having calculated the values of f eff , we now need theinformation about the switching path to obtain R crit ac-cording to Eq. (1). To this end, we construct a Landau-Ginzburg-Devonshire-type first-principles Hamiltonianby expanding the energy around the reference cubicphase, of P m ¯3 m symmetry. The Hamiltonian includesthe most important degrees of freedom of the structure:the strain s ij , the tilts of the oxygen octahedra q i , where q i represents the displacements of the oxygen atoms per-pendicular to the rotation axis and the polar modes u i .(Details on the model can be found in the Supplemen-tary Material.) First, we validate our effective Hamil-tonian H eff by calculating the energetics of the relevantphases (Table S2) and their variation as a function ofexternal pressure (Fig. S1); in both cases we obtain ex-cellent agreement to the first-principles results. Next,we proceed to calculating the most favorable switchingpaths by constraining one component of the polar vec-tor [24, 25] and numerically minimizing (simulated an-nealing) the energy functional with respect to the otherparameters.The resulting double-well potential curves of LiOsO and BaTiO are shown in Fig. 2. Before commentingon the results, it is useful to recall the structural proper-ties of each of the two materials. The structural groundstate of LiOsO has R c symmetry, containing both po-lar distortions and antiphase octahedral tilts ( a − a − a − in Glazer notation) oriented along the [111] pseudocu-bic direction. Since the energy scale associated to theAFD distortions is an order of magnitude larger thanthat associated to u , they are unlikely to be affected bya weak elastic field; in practice, u can only switch be-tween the [111] and [¯1¯1¯1] states. Regarding the actualswitching path, two scenarios are in principle possible. If − − − − − − − − − . − .
04 0 0 .
04 0 . (a) E n e r g y ( m e V / f . u . ) u (Å) u = u = u − − − − − − − − − . − . − .
02 0 0 .
02 0 .
04 0 . (b) E n e r g y ( m e V / f . u . ) u (Å) u = u = 0 FIG. 2. Potential energy landscape for LiOsO (a) andBaTiO (b) from our first principles effective Hamiltonians,obtained by minimizing the energy at fixed u . For LiOsO ,a double-well like curve is obtained when u = u = u is en-forced (dashed line) and a butterfly-like diagram is obtainedwhen all the parameters are allowed to evolve freely (solidline). For BaTiO , the dashed line represents the study un-der the u = u = 0 constraint, and the solid line the casewith no constraints. the non-polar R ¯3 c structure were stable under the con-straint u = 0, the polar modes would be forced to evolvealong the same pseudocubic [111] direction even underthe action of a [100]-oriented external force. However,previous first-principles calculations have shown [26] thatthe R ¯3 c phase has more than one imaginary mode at Γ,which means that R ¯3 c is unlikely to be the saddle point.This suspicion is nicely confirmed by the results of oureffective Hamiltonian: indeed, the “butterfly diagram”of Fig. 2 clearly reflects the presence of a switchablein-plane polarization at u = 0; the resulting coercivefield is F coerc =0.34 eV/˚A. To quantify how much thesystem gains by circumnavigating the energy barrier, weattempted the same computational experiment while im-posing u = u = u along the path; as expected, weobtain a substantially larger critical field of F coerc =0.69 eV/˚A, assuming that the field is still applied along [100].For BaTiO the polarization cannot be constrained bythe tilts, since the latter are absent in this material.At low-temperature BaTiO has R m symmetry, andwe find that the lowest switching barrier occurs whenthe polarization continuously rotates from [111] to [¯111]by passing through an orthorhombic [110] saddle point.(The path is roughly oriented along [100]). We find a crit-ical coercive field of F coerc =0.14 eV/˚A for such a switch-ing path. For comparison to room-temperature experi-ments, where BaTiO adopts a tetragonal structure, wealso calculate the hypothetical barrier that one would ob-tain by constraining P k [100] (i.e. by setting the in-planecomponents of P to zero). We find F coerc =0.29 eV/˚A.This is a substantially larger value than the aforemen-tioned threshold for polarization rotation, in line withliterature results.We are now ready to answer the main physical questionwe asked to ourselves at the beginning: how much do weneed to bend a LiOsO sample to reverse its polar latticedistortion? By means of Eq. (1) we can compute thecritical bending radius for both materials. The obtainedvalues are R crit ∼
118 ˚A for LiOsO and R crit ∼
125 ˚Afor rhombohedral BaTiO and R crit ∼
60 ˚A for tetrago-nal BaTiO . Remarkably, the calculated critical bendingradius of LiOsO is twice as large as that of tetrago-nal BaTiO , essentially matching the calculated value ofrhombohedral BaTiO . Since mechanical switching ofpolar domains in tetragonal barium titanate has alreadybeen experimentally achieved [27] via strain gradients,our results indicate that this is very likely to be feasiblein LiOsO as well.Our results for BaTiO are in good agreement with theones reported in Ref. [28] where a critical bending radiusof 110 ˚A was estimated. This is substantially smallerthan the available experimental estimates (a value of R crit ∼
300 ˚A was observed in BaTiO [27]). This isexpected: theoretical estimations of coercive fields in fer-roelectrics that are based on the homogeneous Landaupotential are typically overestimated by one or two ordersof magnitude [29]. Consideration of more realistic mech-anisms (e.g., domain wall nucleation and motion) woulddrastically complicate our study, and bring us far fromour main scopes. We stress in any case, that our underes-timation of the critical bending radii compared to exper-iments should be ascribed to an overestimation of F crit ,while we regard our calculation of the flexocouplings asaccurate. (The contribution of the oxygen octahedral tiltgradients to the flexocoupling in lithium osmate was ne-glected in this work. While certainly present, we considerit unlikely to qualitatively affect our conclusions; it willbe an interesting topic for follow-up studies.)We expect our results to significantly broaden thescopes of both flexoelectricity and the ongoing search fornew functionalities based on polar metals. In this sense,we believe that our work may open several unexploredresearch directions. First and foremost, we regard anexperimental verification of our predictions as the mostpressing priority. Second, it will be interesting to esti-mate the magnitude of the flexocouplings in a broaderrange of polar metals, and identify candidates where theeffect is especially strong. Finally, from the point of viewof the theory, developing the methodological tools to as-sess the impact of tilt gradients on the calculated coef-ficients is another topic that we regard as promising forfuture studies.We acknowledge the support of Ministerio de Econo-mia, Industria y Competitividad (MINECO-Spain)through Grants No. MAT2016-77100-C2-2-P and No.SEV-2015-0496, and of Generalitat de Catalunya (GrantNo. 2017 SGR1506). This project has received fundingfrom the European Research Council (ERC) under theEuropean Union’s Horizon 2020 research and innovationprogram (Grant Agreement No. 724529). Part of the cal-culations were performed at the Supercomputing Centerof Galicia (CESGA). [1] T. Kim, D. Puggioni, Y. Yuan, L. Xie, H. Zhou,N. Campbell, P. Ryan, Y. Choi, J.-W. Kim, J. Patzner, et al. , Nature , 68 (2016).[2] P. W. Anderson and E. Blount, Physical Review Letters , 217 (1965).[3] Y. Shi, Y. Guo, X. Wang, A. J. Princep, D. Khalyavin,P. Manuel, Y. Michiue, A. Sato, K. Tsuda, S. Yu, et al. ,Nature materials , 1024 (2013).[4] N. A. Benedek and T. Birol, Journal of Materials Chem-istry C , 4000 (2016).[5] E. Bauer and M. Sigrist, Non-centrosymmetric supercon-ductors: introduction and overview , Vol. 847 (SpringerScience & Business Media, 2012).[6] S. Yip, Annu. Rev. Condens. Matter Phys. , 15 (2014).[7] C.-K. Lu and S. Yip, Physical Review B , 104501(2010).[8] C. Ma and K. Jin, SCIENCE CHINA Physics, Mechanics& Astronomy , 97011 (2018). [9] M. Stengel and D. Vanderbilt, (2016).[10] P. Yudin and A. Tagantsev, Nanotechnology , 432001(2013).[11] M. Stengel, Physical Review B , 174106 (2013).[12] S. M. Kogan, Soviet Physics-Solid State , 2069 (1964).[13] F. Vasquez-Sancho, A. Abdollahi, D. Damjanovic, andG. Catalan, Advanced materials , 1705316 (2018).[14] J. Narvaez, F. Vasquez-Sancho, and G. Catalan, Nature , 219 (2016).[15] H. Lu, C.-W. Bark, D. E. De Los Ojos, J. Alcala, C.-B.Eom, G. Catalan, and A. Gruverman, Science , 59(2012).[16] M. Royo and M. Stengel, Physical Review X , 021050(2019).[17] M. Stengel and D. Vanderbilt, in Flexoelectricity inSolids: From Theory to Applications (World Scientific,2017) pp. 31–110.[18] M. Stengel, Physical Review B , 245107 (2016).[19] M. Born and K. Huang, Dynamical theory of crystal lat-tices (Clarendon press, 1954).[20] X. Wu, D. Vanderbilt, and D. Hamann, Physical ReviewB , 035105 (2005).[21] X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken,F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Cara-cas, M. Cˆot´e, et al. , Computer Physics Communications , 2582 (2009).[22] X. Gonze, B. Amadon, G. Antonius, F. Arnardi,L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin,J. Bouchet, E. Bousquet, et al. , Computer Physics Com-munications , 107042 (2020).[23] J. Narvaez, S. Saremi, J. Hong, M. Stengel, and G. Cata-lan, Physical review letters , 037601 (2015).[24] O. Di´eguez and D. Vanderbilt, Physical review letters ,056401 (2006).[25] M. Stengel, N. A. Spaldin, and D. Vanderbilt, NaturePhysics , 304 (2009).[26] H. Sim and B. G. Kim, Physical Review B , 201107(2014).[27] J. Oˇcen´aˇsek, H. Lu, C. Bark, C.-B. Eom, J. Alcal´a,G. Catalan, and A. Gruverman, Physical Review B ,035417 (2015).[28] G. Li, X. Huang, J. Hu, and W. Zhang, Physical ReviewB , 144111 (2017).[29] V. Shelke, D. Mazumdar, G. Srinivasan, A. Kumar,S. Jesse, S. Kalinin, A. Baddorf, and A. Gupta, Ad-vanced Materials , 669 (2011). upplementary notes for ”Switching a polar metal via strain gradients” Asier Zabalo and Massimiliano Stengel
1, 2 Institut de Ci`encia de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain ICREA-Instituci´o Catalana de Recerca i Estudis Avanats, 08010 Barcelona, Spain
I. COMPUTATIONAL DETAILS
Calculations on LiOsO are carried out with pseudopotentials under the generalized gradient approximation (GGA)with Perdew-Burke-Ernzerhof (PBE) exchange-correlation functionals [1]. A dense 16 × × k -point mesh is usedto sample the Brillouin zone with a plane-wave cutoff of 60 Ha. A Gaussian smearing of 0.01 Ha is applied andeffects of spin-orbit coupling are neglected [2]. The dynamical properties, i.e., phonon frequencies and eigenmodes,are calculated using DFPT as implemented in abinit . The independent components of both the elastic and theflexocoupling tensor are computed with Eq. (2) and Eq. (6) of the main text and convergence is tested with respectto the two-dimensional grid of q points. For cubic systems a two dimensional grid is enough to calculate the threeindependent components of the mentioned tensors, commonly denoted as longitudinal ( C , f ), transverse ( C , f )and shear ( C , f ). For a 16 × × q points convergence was clearly reached for the elastic constants, andas a consequence, for the flexocoupling tensor. (Remember that both calculations are carried out using exactly thesame method, so a similar convergence rate and errors are expected.) Calculations for BaTiO are performed by usingLDA and a 12 × × k -point mesh, and the flexocoupling coefficients are extracted via the method described inRef. [3]. For a meaningful comparison with metallic LiOsO we use the conduction band bottom as energy reference,which yields the n -type flexocoupling coefficients as defined in Ref. [3]. (Note that the f component is independentof such a choice.) II. FIRST-PRINCIPLES EFFECTIVE HAMILTONIAN OF LIOSO Expression and parameters
The explicit expression for the first-principles Hamiltonian of LiOsO , which givesthe energy of a 10-atom cell, is the following: H = Ω2 C ( s + s + s ) + Ω C ( s s + s s + s s ) + Ω2 C ( s + s + s )+ β ( q + q + q ) + β ( q + q + q ) + β ( q q + q q + q q )+ ζ ( u + u + u ) + ζ ( u + u + u ) + ζ ( u u + u u + u u )+ ζ ( u + u + u ) + ζ [ u ( u + u ) + u ( u + u ) + u ( u + u )]+ ζ ( u + u + u )+ λ ( s q + s q + s q ) + λ [ s ( q + q ) + s ( q + q ) + s ( q + q )]+ λ ( s q q + s q q + s q q ) + λ ( s q + s q + s q )+ λ ( s q q + s q q + s q q )+ ρ ( s u + s u + s u ) + ρ [ s ( u + u ) + s ( u + u ) + s ( u + u )]+ ρ ( s u u + s u u + s u u ) + ρ ( s u + s u + s u )+ ρ ( s u u + s u u + s u u )+ ρ [ s ( u + u ) + s ( u + u ) + s ( u + u )]+ γ ( u q + u q + u q ) + γ ( u q + u q + u q ) . (1)where the strain is given in the matrix Voigt notation. The calculated model parameters are given in Table S1. a r X i v : . [ c ond - m a t . m t r l - s c i ] J u l TABLE S1. Calculated model parameters for the first-principles Hamiltonian of LiOsO .Parameter Value Units Parameter Value UnitsΩ 112.94 ˚A λ C λ C λ -12.3815 eV/˚A C λ -32.8244 eV/˚A β -4.31776 eV/˚A λ β ρ β ρ -136.627 eV/˚A ζ -35.4382 eV/˚A ρ -110.971 eV/˚A ζ ρ -308.951 eV/˚A ζ ρ -14380.5 eV/˚A ζ -6978.34 eV/˚A ρ -11808.3 eV/˚A ζ γ ζ γ -990.2 eV/˚A Validation of the model
As a first test of the model we compute the energy of the most relevant phases that can beobtained by combining antiferrodistortive and polar instabilities, compared with the reference
P m ¯3 m structure. Thisis, of course, the structural ground state of R c symmetry, containing both polar distortions and antiphase octahedraltilts ( a − a − a − in Glazer notation) oriented along the [111] pseudocubic direction. Other important phases are thecentrosymmetric R ¯3 c and the polar R m , where either the polar or the antiferrodistortive modes are suppressed. Aswe can see from the calculated values (Table S2), the model accurately reproduces the first-principles results obtainedby fully relaxing the crystal cell within the given symmetry. In particular, this first-principles Hamiltonian almostperfectly describes the energy difference between the lowest-energy R ¯3 c and the R c phases, which are the mostrelevant for ferroelectric switching: -22 meV/f.u. with DFT and -23 meV/f.u. with our Hamiltonian, were f.u. standsfor perovskite formula unit. As a further test, we study the phase diagram of LiOsO as a function of cell volume.The dots shown in Fig. S1 were obtained via explicit DFT calculations by fully relaxing the ionic positions for eachfixed cell volume. These results are in excellent agreement with experimental values [4] and previous first-principlesstudies [5]. Solid lines represent the energies obtained with our model, again showing excellent agreement with thefirst-principles values. TABLE S2. Volume and energy comparison between first-principles calculations and our model. Energies are in meV /.f.u.units. DFT Effective Hamiltonian∆ E ∆ V (%) ∆ E ∆ V (%) R m -149 1.1 -157 1.3 R ¯3 c -908 -11.4 -906 -11.5 R c -930 -12.0 -929 -11.9 − − . − . − . − . . E n e r g y ( e V / f . u . ) Volume (Å / f.u.) P m ¯3 mR mR ¯3 cR c FIG. S1. Energy per formula unit as a function of the volume of the perovskite formula unit of
P m ¯3 m , R m , R ¯3 c , and R c phases. Solid lines are obtained from our effective Hamiltonian. III. FIRST-PRINCIPLES EFFECTIVE HAMILTONIAN OF BaTiO Expression and parameters
The explicit expression for the first-principles Hamiltonian of BaTiO , which givesthe energy of a 5-atom cell, is the following: H = Ω2 C ( s + s + s ) + Ω C ( s s + s s + s s ) + Ω2 C ( s + s + s )+ ζ ( u + u + u ) + ζ ( u + u + u ) + ζ ( u u + u u + u u )+ ρ ( s u + s u + s u ) + ρ [ s ( u + u ) + s ( u + u ) + s ( u + u )]+ ρ ( s u u + s u u + s u u ) + ρ ( s u + s u + s u )+ ρ ( s u u + s u u + s u u )+ ρ [ s ( u + u ) + s ( u + u ) + s ( u + u )] (2)The calculated model parameters are given in Table S3. TABLE S3. Calculated model parameters for the first-principles Hamiltonian of BaTiO .Parameter Value Units Parameter Value UnitsΩ 60.745 ˚A ρ -911.228 eV/˚A C ρ -97.0643 eV/˚A C ρ -99.2179 eV/˚A C ρ -29262.0 eV/˚A ζ -11.5402 eV/˚A ρ ζ ρ ζ — — — Validation of the model
As with lithium osmate, we first compute the energy of the relevant phases. Since thereare no tilts, this reduces to calculating the energy of the ground state R m phase. Fig. S2 shows the phase diagramof the cubic and R m phases as a function of cell volume, showing an excellent agreement with the first principlescalculations. − . − . . . . . E n e r g y ( e V / f . u . ) Volume (Å / f.u.) P m ¯3 mR m FIG. S2. Energy per formula unit as a function of the volume of the perovskite formula unit of
P m ¯3 m and R m phases ofBaTiO . Solid lines are obtained from our effective Hamiltonian. IV. FLEXOCOUPLING COEFFICIENTS
Convergence of flexocoupling coefficients
The convergence of the flexocoupling coefficients (in the type-IIrepresentation) of LiOsO with respect to the q point mesh, for a fixed k point mesh, is analyzed. Numerical valuesare shown in Table S4 and graphical results are given in Fig. S3. E l a s t i cc o n s t a n t s ( G P a ) q points C C C − F l e x o c o up li n g c o e ffi c i e n t s ( e V ) q points f f f FIG. S3. Convergence of the elastic tensor and the flexocoupling tensor with respect to the two-dimensional grid of q points( q × q × q points 2 × × × × × × × × C C C f II11 -23.84 -4.53 -15.62 -13.82 f II12 f II44 -3.43 5.43 3.24 3.27
Effective flexocoupling coefficients
In order to have a self contained manuscript, here we give the expressions ofthe effective flexocoupling coefficients for the [100], [110] and [111] orientations for the beam-bending limit followingthe definitions of Ref. [6]. f = − C C + C f II11 + C C + C f II12 ,f = Af II11 + Bf II12 − − A ) f II44 ,f = C C + 2 C + C ( f II11 + 2 f II12 ) − C C C + 2 C + C f II44 , (3)where A = 2 C C ( C − C )( C + 2 C ) + 2 C C ,B = 2( C − C ) C ( C − C )( C + 2 C ) + 2 C C . (4) [1] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters , 3865 (1996).[2] N. A. Benedek and T. Birol, Journal of Materials Chemistry C , 4000 (2016).[3] M. Stengel, Physical Review B , 245107 (2016).[4] Y. Shi, Y. Guo, X. Wang, A. J. Princep, D. Khalyavin, P. Manuel, Y. Michiue, A. Sato, K. Tsuda, S. Yu, et al. , Naturematerials , 1024 (2013).[5] H. Sim and B. G. Kim, Physical Review B , 201107 (2014).[6] J. Narvaez, S. Saremi, J. Hong, M. Stengel, and G. Catalan, Physical review letters115