Electrostatic engineering of strained ferroelectric perovskites from first-principles
Abstract
Design of novel artificial materials based on ferroelectric perovskites relies on the basic principles of electrostatic coupling and in-plane lattice matching. These rules state that the out-of-plane component of the electric displacement field and the in-plane components of the strain are preserved across a layered superlattice, provided that certain growth conditions are respected. Intense research is currently directed at optimizing materials functionalities based on these guidelines, often with remarkable success. Such principles, however, are of limited practical use unless one disposes of reliable data on how a given material behaves under arbitrary electrical and mechanical boundary conditions. Here we demonstrate, by focusing on the prototypical ferroelectrics PbTiO and BiFeO as testcases, how such information can be calculated from first principles in a systematic and efficient way. In particular, we construct a series of two-dimensional maps that describe the behavior of either compound (e.g. concerning the ferroelectric polarization and antiferrodistortive instabilities) at any conceivable choice of the in-plane lattice parameter, , and out-of-plane electric displacement, . In addition to being of immediate practical applicability to superlattice design, our results bring new insight into the complex interplay of competing degrees of freedom in perovskite materials, and reveal some notable instances where the behavior of these materials depart from what naively is expected.
pacs:
77.55.Nv, 61.50.Ks, 68.65.Cd, 77.80.bnI Introduction
Perovskite-structure oxides with ABO formula display a wide range of functional properties that are sought after for applications in energy and nanoelectronics. Such functionalities are, to a large extent, governed by a number of symmetry-lowering lattice instabilities, which naturally occur in the crystals depending on the chemical nature of the A and B cations. Ferroelectric (FE) and antiferrodistortive (AFD) modes are by far the most common distortions in ABO compounds. The former are typically associated with “soft” (i.e. unstable) zone-center phonons of the reference cubic phase (see Fig. 1a) cohen90 ; cohen92 . Whenever they dominate, the crystal possesses a spontaneous (and switchable) electric polarization, which can be exploited for information storage and (via the piezoelectric effect) for electromechanical transduction. The latter, on the other hand, are related to zone-boundary instabilities, whereby the oxygen octahedral cages surrounding the B cation cooperatively rotate (either in antiphase or in phase) about a particular axis (see Fig. 1a). While AFD modes are non-polar, they are extraordinarily important for the physics of magnetic perovskites, where the couplings between neighboring transition-metal sites are sensitive to the BOB bond angle. Cases where FE polarization and AFD distortions are present and mutually coupled in the same phase are especially intriguing, both from a practical and a fundamental perspective. Indeed, recent research has unvealed a number of notable examples where octahedral rotations can strongly enhance, or even induce, a ferroelectric distortion in materials that would be otherwise nonpolar. junquera12 ; bousquet08 Most importantly, such unconventional forms of ferroelectricity, often referred to as “improper” levanyuk74 or “hybrid improper” benedek11 ; benedek12 , appear as a promising route for achieving electrical control of magnetism in multiferroics and magnetoelectrics, a tantalizing goal that has remained so far elusive in spite of the numerous recent breakthroughs heron14 ; fusil14 .
Practical ways to manipulate and design the aforementioned functionalities are currently being investigated along two main avenues. These consist in either synthesizing new perovskite materials (e.g. by playing with the A- and B-site chemistry), or in altering the behavior of existing ones via external perturbations. In the context of the latter strategy, a large part of the recent efforts have been directed at appropriately modifying the strength of the FE and AFD instabilities (and/or their mutual interplay) via the imposed mechanical hatt10 ; rondinelli11 ; rondinelli11b ; may10 and electrostatic boundary conditions wu12 ; stengel12b ; hong13 (see Fig. 1b for a summary of the most common trends).
Strain engineering, for instance, consists in growing the material in a thin-film form, where bonding to the substrate favors coherency in the in-plane registry. Thanks to the availability of numerous substrates with a variety of lattice parameters, such a practice has become very popular in the past few years dawber05b ; choi04 . Polar (FE) and structural (AFD) degrees of freedom are generally sensitive to the applied strain, and this often allows for the stabilization of phases with enhanced functional properties which would be otherwise inaccessible at the bulk level bea09 ; warusawithana09 . (The regions of parameter space that lie close to phase boundaries are particularly responsive to external perturbations because of their structural “softness” jacek10 .)
Control of electrostatics in perovskites, on the other hand, is comparatively more subtle, as it can be easily thwarted by a redistribution of the mobile charges (e.g. electron/hole carriers, ionic defects, etc.) or the formation of ferroelectric domains. Also, leakage currents often limit the magnitude of the electric field that can be realistically applied via an external voltage source junquera11 . Yet, the so-called electrostatic coupling junquera11 ; ghosez08 ; zubko12 ; dawber12 ; dawber05 ; dawber07 is one of the most valuable design principles in ferroelectric superlattices based on perovskites. It states that, at the boundary between two insulating layers, the normal component of the electric displacement field must be preserved; this effectively constrains the out-of-plane polarization state of the stack to a common value that, depending on the materials choice and conditions, might be markedly different than the spontaneous polarization of the parent compounds. In the vast majority of cases, electrostatics and strain work together, and need to be simultaneously taken into account when interpreting the experiments with quantitatively predictive models.
First-principles simulations within the framework of density-functional theory (DFT) have been of invaluable help in rationalizing and guiding the experimental efforts towards synthesis of new materials with enhanced properties. Until few years ago DFT methods were restricted to zero-electric-field simulations, which complicated the description of electrostatics in multicomponent oxide systems. Such effects have traditionally been described by means of phenomenological models, where the DFT input was limited to the computation of the relevant coupling coefficients. Recent and timely advances in first-principles computational methods have now enabled full control over the macroscopic electrical variable of interest (i.e., electric field, polarization and electric displacement) in the simulation of ferroelectric, magnetoelectric, and piezoelectric materials souza02 ; dieguez06 ; stengel09c . These methods now allow to systematically and accurately explore both strain and electrostatic effects at the level of the individual bulk material. In particular, by calculating the ground-state of a given compound as a function of the in-plane strain and out-of-plane electric displacement, one can trace two-dimensional maps that comprehensively describe all configurations that are, in principle, accessible in a superlattice geometry, i.e. at arbitrary electrical and mechanical boundary conditions.
Among the variety of materials that have been considered as building blocks for ferroelectric superlattices, perovskite-structure PbTiO (PTO) and BiFeO (BFO) occupy a special place. Indeed, the former is the end member of the PbZrTiO family, the technologically most important piezoelectric material, while the latter is the most intensely studied of all naturally occurring multiferroics. Interestingly, even if both materials present active (and strong) FE and AFD instabilities in their cubic reference phase, their ground state is remarkably different. On one hand, PTO condenses in a purely ferroelectric phase, where the AFD instability has been suppressed by the large polar distortion. (FE and AFD modes generally, although not always benedek11 ; benedek12 , compete.) Note that the AFD modes, even if they are absent in the low-temperature bulk phase, have been shown to play an important role in the ferroelectric paraelectric phase transition wojdel13 , and are predicted to become active in thin films at medium to large values of the epitaxial strain yang12 . Bulk BFO, on the other hand, adopts a rhombohedral () structure where a large polarization and antiphase AFD distortions [both oriented along the pseudocubic (111) direction] coexist. The balance between polar and non-polar instabilities is extraordinarily subtle in this latter material, leading to a rich and complex energy landscape with many local minima dieguez11 . In fact, numerous phase transitions driven by pressure guennou11 , temperature prosandeev12 ; rahmedov12 ; cazorla13 , and expitaxial strain bea09 ; zeches09 have been recently predicted and observed in BFO. Given the extraordinarily rich behavior of both materials, it is hardly surprising that their combination results in an even broader variety of original phenomena. Indeed, recent ab initio simulations have shown that epitaxial short-period superlattices made of PbTiO and BiFeO layers may exhibit strain-driven isostructural phase transitions, exotic metastable phases, and unusually large oxygen octahedra tilts yang12 ; stengel12 ; cazorla14 . In view of the great fundamental and technological interest of bulk PTO, BFO, and layered heterostructures comprising both species, a detailed knowledge of how mutually coupled FE and AFD degrees of freedom in these materials respond to external electric fields and/or mechanical constraints, is highly desirable.
In this article, we present a systematic first-principles study of the energy, structural and dielectric properties of ‘‘bulk’’ ^{1}^{1}1We perform bulk calculations, but bearing in mind a thin-film and/or superlattice geometry, where such conditions of strain and electric field can be experimentally achieved. PTO and BFO as a function of both in-plane epitaxial strain () and out-of-plane electric displacement field (). Our computational results are presented in the form of two-dimensional maps, which describe the structural and electronic ground state of the material at an arbitrary choice of the electrical and mechanical boundary conditions. These maps provide a comprehensive overview of all the phases of BFO and PTO that can be realized via strain engineering and electrostatic coupling, i.e. by growing the material in a thin-film or superlattice form. Most importantly, these maps illustrate how the main functional properties (e.g. dielectric, piezoelectric) of either material evolve as a function of , and thus constitute an invaluable guide in the experimental search for new artificial materials. For instance, we identify the specific regions of the phase diagram where PTO or BTO display an unusual dielectric and structural softness, leading to very large dielectric and piezoelectric responses. For PTO, this region is located at small positive strain, where the polarization abruptly rotates from out-of-plane to in-plane. (Interestingly, in a very small range of in-plane strain values, we find a stable monoclinic structure, where the polarization has both in-plane and out-of-plane components and is accompanied by an out-of-phase AFD tilt mode.) In the case of BFO, such a structural softness occurs in the orthorhombic phase that becomes stable at high tensile strain, and whose existence was recently predicted by means of first-principles calculations yang12 ; stengel12 .
In addition to their practical interest, our results also provide a wealth of valuable insight into the physics of either material, especially regarding the mutual interplay of strain, polarization and AFD tilts. In this context, parallel to rationalizing the results of earlier works on PTO- and BFO-based systems junquera12 ; yang12 ; stengel12 we could identify a notable case where the main order parameters do not follow the typical behavior that has been traditionally assumed in phenomenological theories. In particular, both PTO and BFO are characterized, at high tensile strain, by an in-phase out-of-plane AFD mode, whose amplitude grows, rather than decreasing, when an out-of-plane electric field (of arbitrary strength) is applied. This is in stark constrast with the expected trends (AFD modes usually tend to compete with the polar degrees of freedom) represented in Fig. 1b.
This work is organized as follows: In Sec. II we provide the details of the computational strategy and ab initio methods that we used in the construction of the diagrams. In Sec. III, we present our results and discuss them extensively. Finally, we conclude by summarizing our main findings in Sec. IV.
Ii Computational Details
ii.1 Construction of the maps
Our calculations are performed within the local spin density approximation to density-functional theory, as implemented in the “in-house” LAUTREC code. We apply a Hubbard eV to the Fe ions yang12 ; kornev07 . In all the simulations, we use the 20-atom simulation cell depicted in Fig. 1c, which allows for reproducing the ferroelectric and antiferrodistortive distortions of interest (i.e., in-phase AFD and out-of-phase AFD, see Fig. 1a).
It is worth noting that working at fixed electric displacement is in many respects analogous to working at fixed electric field . At fixed , however, a much broader range of configurations can be explored in FE materials, including regions of parameter space where the system is unstable with respect to a polar distortion stengel09c ; stengel09 . Controlling is akin to controlling , and therefore provides an intuitive and direct conceptual link to phenomenological descriptions (i.e., Landau theory) where is treated as the main order parameter.
Atomic and cell relaxations are performed by constraining the value of the out-of-plane component of stengel09c and the two in-plane lattice vectors and (see Fig. 1). In the remainder of this article, will stand for the out-of-plane component of the electric displacement vector. Small tolerances of eV/Å and Kbar are imposed on the atomic forces and stresses. We correct for the bias deriving from the Pulay stress in the out-of-plane direction, , by introducing a compensating term (where is the cell volume) in the expression of the internal energy. The calculations are repeated in order to span the physically relevant space defined by the set of parameters .
In order to achieve an accurate description of the energy and structural properties of PTO and BFO in the Å and intervals, we carry out explicit calculations on a fine regular grid of points, spaced by Å and (where is the electron charge and the surface of the five-atom unit cell). Spline interpolations are subsequently applied to the calculated data points, in order to obtain smooth two-dimensional maps (see Figs. 3-7).
ii.2 Electrostatic and distorsion analysis
The components of the total polarization in each relaxed configuration, i.e., in-plane and out-of-plane , were computed with the formula
(1) |
where is the volume of the cell, runs over all the atoms there, represent Cartesian directions, is the displacement vector of the -th atom as referred to the cubic perovskite phase, and the Born effective charge tensor calculated in the paraelectric reference state. We note that since both PTO and BFO materials are good ferroelectrics, one can regard the out-of-plane component of the total polarization, , to be approximately equal to .
The magnitude of the antiferrodistortive AFD, AFD (both with the rotation axis oriented along ) and AFD (with axis) rotations were evaluated by considering the projection of the vectors on the corresponding zone-boundary eigenmode, see Fig 1a. (We use a “” or a “” symbol to indicate in-phase and out-of-phase tilt modes, respectively.)
ii.3 Computation of functional properties
The availability of the full maps for a given material allows us to extract important functional properties (e.g. dielectric and piezoelectric coefficients) in a straightforward way as a by-product of our calculations stengel09c . In particular, the inverse dielectric constant of the system, (by without subscripts we shall always imply the 33 component of the dielectric tensor), adopts the simple form
(2) |
where is the internal energy, the unit cell volume, the unit cell surface, and the unit cell out-of-plane parameter. Also, the longitudinal piezoelectric coefficient is readily obtained as
(3) |
where direction [stengel09c, ]. Likewise, the shear piezoelectric coefficient is expressed within our formalism as is the potential drop across the unit cell in the
(4) |
where corresponds to the projection of the lattice vector on the basal plane defined by and (see Fig. 1). (The Voigt notation of the shear strain, therefore, refers to the rotated axes of the 20-atom cell, and not to the pseudocubic axes.)
The response coefficients defined in Eq. (2), Eq. (3) and Eq. (4) are directly related to those discussed by Wu, Vanderbilt and Hamann (WVM) wu05 . In particular, the dielectric constant of Eq. (2) corresponds to the free-stress permittivity of WVM, . Regarding the tensor, our notation is consistent with WVM. It is useful for later use (again, following WVM) to introduce the fixed-, free-stress piezoelectric coefficients
(5) |
which trivially relate to the corresponding component of the -tensor as .
Note that the above definitions can be readily used to predict the dielectric and piezoelectric properties of a multicomponent system (e.g. a layered superlattice) based on the bulk properties of the constituents. For example, consider a superlattice (SL) with a total thickness of material A and of material B. In this case, one has
(6) |
As above, one can then easily recover the overall -coefficients as .
Iii Results and Discussion
iii.1 Bulk PbTiO
In Fig. 2, we show a schematic phase diagram of PTO, calculated as a function of , while in Fig. 3 we report a series of detailed 2D maps that illustrate the evolution of most relevant structural and electrical degrees of freedom within the same range. For each region of the phase diagram (Fig. 2) we indicate the space group of the stable crystal phase therein and we mark the resulting two-phase boundaries as a dashed line. (The same phase boundaries are also shown in the maps of Fig. 3.) Moving from smaller to larger values of , the dashed lines correspond, respectively, to the suppression of AFD tilts, the onset of and the onset of AFD rotations (see Fig. 3). We note that, according to our calculations, all the phase transitions occurring in PTO are of second-order type, i.e., the energy of the crystal varies continuously while specific distortions – either FE, AFD, or both simultaneously – emerge or are suppressed. The solid red line in Figs. 2 and 3 indicates, for each value of , the value of that minimizes the energy; the points of this curve therefore correspond to the electrostatic equilibrium configurations as a function of the in-plane strain. Note that exactly coincides with along such path, as at the stationary points of the curve the macroscopic electric field vanishes by definition. The curve mostly follows the known trends for PTO thin films: yang12 ; stengel12 (i) A zero or negative in-plane strain favors a tetragonal state with a strong out-of-plane component of the polarization, . PTO remains in the tetragonal phase down to Å, where the resulting ratio and out-of-plane polarization become very large ( and C/m, respectively). (ii) is suppressed upon application of a tensile strain, whereby the polarization aligns with the in-plane direction (the orientation of the polarization abruptly changes from to at Å) and the system transitions to the orthorhombic phase. This is analogous to the phase occurring in BaTiO under similar conditions dieguez04 . Finally, (iii) at a large enough value of the tensile strain a second phase transition occurs, where an AFD mode appears; as a consequence, the symmetry of the system reduces to as recently reported in Refs. yang12 ; stengel12 .
While the behavior described above is fully consistent with the results of earlier studies we found, however, some interesting surprises as well. First, we find a stable monoclinic phase in the narrow Å interval. This phase is characterized by and (see Fig. 4) and, to our knowledge, has never been predicted or observed previously. (The energy difference between the monoclinic and orthorhombic phases is very small, i.e., mev/f.u., which explains why the monoclinic structure might have been missed in earlier studies.) Second, our calculated phase diagram of Fig. 2 clearly shows that AFD distortions do not occur in PTO at any value of or . Earlier literature studies yang12 ; stengel12 ; blok11 predicted such distortions to occur in the tensile-strain region, where the polarization vector has an in-plane orientation. (These structures were, therefore, identified as orthorhombic , rather than .) This discrepancy can be rationalized by observing that the AFD mode amplitude was found yang12 ; stengel12 ; blok11 to be very small (i.e., ), which implies that that the and phases are essentially degenerate in energy. Whether they show up or not in a calculation might therefore depend on subtle details of the numerical implementation.
Now, we shall discuss more in detail the results that we obtained by exploiting the full capabilities of our computational approach, i.e., by imposing arbitrary , and including those configurations that are characterized by a nonzero internal field along the direction.
First of all, it is interesting to note the unusually flat potential-energy landscape in a vicinity of the zero-strain equilibrium state (see Fig. 3); this indicates that the system is extremely sensitive to external perturbations in this region of the phase diagram. Indeed, recall that in the orthorhombic tetragonal transition the polarization abruptly rotates from in-plane to out-of-plane; thus, large dielectric and piezoelectric responses are expected to occur near Å(we will comment more extensively on this point in Sec. III.3).
The structurally most complex region of the phase diagram occurs in the moderately compressive strain regime ( Å) and values that are smaller than the equilibrium . Here all the most important structural modes of PTO (, and AFD distortions) coexist, and the symmetry of the crystal is monoclinic (see Figs. 2 and 3). (It reduces to orthorhombic at .) The occurrence of the and AFD distortions, which are active instabilities of the reference structure but are absent in the zero-field equilibrium phase within this strain interval, may appear surprising. This, however, can be easily rationalized in terms of the mutual competition between different instabilities in PTO. At zero electric field (i.e., red curves in Figs. 2 and 3), is large enough to prevail over the other degrees of freedom, which are completely suppressed. By reducing (as can occur, in practical situations, as a consequence of depolarizing field effects), the aforementioned cross-couplings are weakened, allowing the competing and AFD distortions to reenter the ground-state structure. This finding explains the results of a recent study junquera12 , where a superlattice geometry was found to induce a monoclinic structure analogous to that described here in the PTO layers. Our results indicate that the electrical boundary conditions are, in fact, likely to be the main cause for such an outcome. Higher in-plane compression ( Å) suppress the in-plane polarization instabilities completely, and at the same time increase the strength of AFD rotations; the symmetry of the crystal is tetragonal (tetragonal at ) Our data suggest that the AFD rotations will consistently show up at arbitrary (moderate to strong) compressive strain values, provided that a depolarizing effect is present.
In the moderate tensile strain regime ( Å), all oxygen octahedral rotations are suppressed and the symmetry of the system is monoclinic . (This corresponds to a polarization oriented along , with , where the out-of-plane component can be induced by applying an out-of-plane electric field to the equilibrium phase.) Finally, in the Å interval, the monoclinic angle increasingly deviates from as is increased, and the AFD rotations become more prominent. Note that, as increases, the in-plane component of the polarization becomes smaller (although full suppression of is not observed). These features are consistent with a monoclinic phase that is closely related to the equilibrium orthorhombic found at . We note that the trend observed in this region is quite surprising: as , or equivalently , increases the AFD rotations become larger. This contrasts with the typical trends shown in Fig. 1. (We already pointed out that such a trend was indeed followed by the AFD / modes in the compressive regime.) Thus, a cooperative rather than a competitive coupling appears to be operating between the and AFD degrees of freedom. This observation points to new interesting opportunities for tuning of the AFD oxygen rotations via application of an external electric field.
iii.2 Bulk BiFeO
For the case of BFO, we shall present our results by following a similar strategy to what we did in the PTO case: The overall phase diagram is shown in Fig. 5, while Figs. 6 and 7 illustrate the detailed behavior of the individual degrees of freedom. (We split the relevant parameter space into two parts here, to better describe the large compressive strain region, where the spontaneous polarization adopts very large values.) We use the same convention as in Figs. 2 and 3 regarding the electrostatic equilibrium path and the two-phase boundaries, with the only difference that, contrary to the PTO case, the latter are almost always shown as solid lines here. This reflects an important outcome of our calculations: most phase boundaries in BFO turn out to be of the first-order type. (The only exception is the transition to the state at very small values of , marked in the upper left corner of Fig. 5.)
As before, the path described by the states (minimum energy configurations at a given value of ) is consistent with the results of earlier works yang12 ; stengel12 : (i) For a large range of in-plane strain values around the equilibrium parameter, BFO adopts a monoclinic (we shall indicate this phase as -I) state that is closely related to the bulk rhombohedral ground state. Here, both the polarization and the AFD vector are large and roughly oriented along the (111) pseudocubic direction. Their evolution with decreasing follows the usual trend: and AFD display a progressive enhancement, while the in-plane components of both degrees of freedom simultaneously decrease in amplitude. (ii) At a sufficiently large compressive strain ( Å), the AFD distortion amplitude peaks and then rapidly falls to zero, while both the ratio and undergo a drastic increase. (At Å, for instance, the value of is C/m and .) We identify this abrupt change with the isosymmetric hatt10b transition to the well-known supertetragonal phase of BFO (also known in the literature as “T-like”) that was experimentally observed in epitaxial films under a large compressive strain bea09 ; zeches09 . (iii) At the other extreme end of the phase diagram (i.e. in the strongly tensile regime, Å) BFO stabilizes in an orthorhombic phase, in all respects analogous to the structure that occurs in PTO at similarly large strain values (see previous section).
When we move away from the electrostatic equilibrium curve, BFO reveals a number of remarkable (and potentially useful) features, as we shall illustrate in the following. In the zero or mildly compressive strain regime ( Å) we observe at small values a region where BFO adopts an orthorhombic structure. This state exhibits very large AFD and AFD O rotations, and a nonzero . The latter, however, is not “ferroelectric” in the usual sense: the system here behaves like a nonlinear dielectric, with a metastable local minimum at , where the symmetry becomes orthorhombic . We can thus expect that under a strong depolarizing field (and conditions that prevent domain formation) a BFO film might adopt a ground state, which was so far only observed in bulk at high pressure haumont09 . As is decreased, the -interval over which this orthorhombic phase occurs becomes broader, while the amplitude of the AFD and AFD oxygen rotations remain practically unchanged. At even larger compressive strains ( Å) we find a new region, this time at large values of , where the AFD oxygen octahedral rotations and in-plane polarization completely disappear and the crystal becomes tetragonal . According to our calculations, the electric field-induced monoclinic -II tetragonal transformation is of second-order type. Note that this same phase could be, in principle, obtained even without the application of an electric field, by further compressing the crystal in-plane [stengel12, ]. However, the strain that one would have to apply is enormous ( Å) – our calculations show that, by combining the effects of in-plane compression and out-of-plane electric field, one might have better chances of experimentally stabilizing the high-energy phase, which is probably out of reach if one relies solely on mechanical means. Finally, moving now to the large tensile strain region ( Å), BFO adopts the same monoclinic phase that we have found in PTO under similar conditions (see Sec. III.1). Also in this case it is interesting to note that the AFD O rotations become larger in amplitude as , and consequently , is increased. Remarkably, the potential energy surface around appears to be rather flat, which indicates a very large dielectric polarizability; this has important implications for practical applications, as we shall discuss shortly in Sec. III.3.
iii.3 Functional properties
The results described above have direct implications over the functional properties of the two materials, as we shall see in the following. We shall focus more specifically on the dielectric and piezoelectric properties, whose dependence on the external parameters will be again illustrated in the form of 2D maps. We shall discuss the two materials separately, starting with PbTiO.
iii.3.1 Bulk PbTiO
In Fig. 8a, we show the calculated inverse dielectric constant, , of bulk PTO in the region of parameter space spanned by Å and . First of all, note the blank area in the lower left corner of the graph; it corresponds to points where adopts negative values, pointing to an instability of the system towards a ferroelectric distortion. As these points are experimentally inaccessible, we have decided to exclude them from our analysis. The most interesting region in the map is that where the crystal, under zero-field conditions, transitions (for increasing ) from the orthorhombic to the tetragonal transition via the intermediate monoclinic phase (see Sec. III.1). Here In the Å interval, for instance, adopts values which are higher than (in units of , that is, the dielectric constant of vacuum) in the corresponding zero-field states (indicated with a magenta line). The phase transition corresponds to an abrupt rotation of the polarization vector from in-plane to out-of-plane, and thus the large dielectric response of the system in that region stems from the flatness of the potential landscape that governs . Indeed, in the tensile-strain regime the system behaves as a normal dielectric along , i.e. the internal energy of the system depends quadratically on . Conversely, in the compressive regime, adopts a symmetric double-well form (see Fig. 4). In crossing the boundary, flattens significantly in the region near , where its second derivative, and consequently , becomes very small.
As they are directly proportional to , the and piezoelectric coefficients of bulk PTO are expected to be giant (see Eqs. 3 and 4 in Sec. II.3) in the transition region described above. Of course, a necessary condition is that the corresponding and coefficients do not vanish. Regarding the longitudinal coefficient, Fig. 8b shows a relatively simple monotonous behavior, with a roughly linear increase from zero (at the state, where it vanishes by symmetry) to (depending on the strain state) at the highest values of . (Interestingly, the rate of increase of with , which can be regarded as an effective electrostrictive coefficient, seems to be rather sensitive to the in-plane strain, being largest in the extreme compressive regime, and milder under tensile conditions.) By combining the information of and , one can draw the following conclusions: (i) vanishes at tensile strain by symmetry, unless a sufficiently large electric field is applied to the crystal. (ii) is nonzero but relatively small in the compressive strain regime, where the zero-field curve lies far from the region where diverges. (One can, in principle, try to approach that region by applying an -field antiparallel to , although the crystal will likely switch to the opposite polarization state well before reaching it.) (iii) In the intermediate region where PTO is either monoclinic or tetragonal (but close to the phase boundary), becomes giant as speculated above, due to the structural and dielectric softness of the system therein.
As for the shear piezoelectric coefficient, (see Fig. 8c), we find a very rich behaviour that is remarkably sensitive to the mechanical and electrical boundary conditions. First, note that whenever vanishes (i.e. at compressive strain), vanishes because of symmetry (the crystal axes form orthogonal angles with no shear). The physically interesting region occurs from zero to moderately tensile strain, where the in-plane components of the polarization start being active. Surprisingly, the coefficient shows a rather peculiar behavior: it is positive in most of the region of parameter space where is nonzero (this includes all equilibrium zero-field states with symmetry, occurring at Å), except for a narrow region close to the phase boundary where it becomes negative.
In order to understand the physical origin of this result, it is useful to recall the form of the coupling terms in the free energy that involve polarization and shear strain. At the lowest order, in PTO such a coupling is governed by pertsev98
(7) |
where is the corresponding component of the stress tensor (note that the shear occurs between the out-of-plane [001] and in-plane [110] axes), and is the coupling coefficient. The piezoelectric coefficient can be roughly estimated by deriving the above coupling term twice with respect to (which in the context of our calculations essentially corresponds to ) and ,
(8) |
It is useful, for analysis purposes, to break down the above equation into two intermediate steps, by introducing the shear strain, ,
(9) |
One has, trivially,
(10) |
This behavior of , which corresponds to the deviation of the monoclinic angle from 90, is qualitatively well reproduced in our calculations (see Fig.3). From the evolution of as a function of (for a fixed ), one can immediately appreciate the origin of the sign change in : close to , grows linearly with (i.e. its derivative with respect to is positive), reaches a peak and then rapidly drops back (negative region) to zero when approaching the phase boundary.
In a close vicinity of such boundary, has a finite value, and therefore can be approximated by a constant; one should then have
(11) |
Now, since the transition is of second-order type, for a given one expects a critical behavior of with respect to of the form
(12) |
where describes the phase boundary. By combining Eq. (11) and Eq. (12), and by further taking the first derivative in , one obtains that the piezoelectric coefficient diverges at the boundary as
(13) |
Unfortunately, the divergence does not show up in our graphs because of the relatively coarse grid of points that we have used to sample our parameter space. Nevertheless, the above conclusion is solid, as it stems from basic symmetry arguments and from the second-order nature of the phase transition. (The latter is unambiguously demonstrated by our numerical results.) This points to a promising route to achieving a large piezoelectric response by tuning the -coefficient , rather than the dielectric constant. We expect this mechanism to be especially useful in superlattices, because of the additive property [see Eq. (6)] of the -coefficients. [Conversely, the dielectric constant obeys the typical series-capacitor model, Eq. (6), making it harder to enhance this latter functionality in a layered system.]
iii.3.2 Bulk BiFeO
The calculated inverse dielectric constant of BFO, expressed as a function of the in-plane lattice parameter and out-of-plane electric displacement, is shown in Fig. 9a. As in the PTO case, we have excluded from the graph the region where the system is unstable against a polar distortion, and hence has a negative permittivity (white area). (In the region where BFO adopts the phase is in fact positive however, for the sake of simplicity and because that region is experimentally difficult to access, we ignore it in the following analysis). The data points that lie closest to the white area are characterized by a vanishing value of , and hence by a diverging dielectric constant, which could be interesting for applications. At large tensile strain, presents some discontinuities, which occur more specifically along the -I phase boundary. (Recall that the involved phase transformation is of first-order type, and is characterized by a drastic structural change; this explains the abrupt jump in the response properties of the crystal.) Interestingly, at Åand for moderate values of ( ), where the system is in a monoclinic phase (orthorhombic at ), adopts very large values. In fact, the high-tensile strain phase appears to be on the verge of having a polar instability along . This fact reflects itself in an unusual flatness of the energy curves around , which in turn results into a very large dielectric (and piezoelectric, as we shall see shortly) response of the system (see Eq. 2).
Regarding the longitudinal piezoelectric coefficient (see Fig. 9b), , BFO behaves roughly like PTO: is zero at , grows linearly for increasing and is largest in the upper left corner of the plot, i.e. at large compressive strain and large out-of-plane polarization. (Of course, the first-order phase transition lines and the corresponding discontinuities in are absent in PTO; however, these do not produce a profound qualitative change in the overall picture.) Note that, at difference with PTO, the system does not cross a state of marked structural softness at small applied strains (the zero-field line stays far away from the high- region unless a very large tensile strain is applied, while in PTO it crosses through it at Å, see Fig. 8a). Thus, as in the large tensile strain region (and hence ) vanishes by symmetry at zero electric field, we conclude that the transition from -I to -II (which occurs at large compressive strain) constitutes the most promising configuration for obtaining an enhanced value of in BFO.
The shear piezoelectric coefficient, on the other hand, displays a complex behavior that differs qualitatively from the PTO case and that does not lend itself to a comparatively simple physical interpretation (see Fig. 9c). As a matter of fact in BFO we have, in addition to the polarization-strain coupling that we mentioned in the previous Section, a further term in the free energy that couples the strain to the AFD modes,
(14) |
where is a second coupling coefficient and are the AFD distortion amplitudes. It is difficult, in the context of our calculations, to disentangle the effects of the two different coupling terms; for this reason, in the following we shall limit ourselves to commenting on the general trends that emerge from our data.
To start with, we note that the system in the proximity of zero epitaxial strain is characterized by small values (the zero-field curve stay close to the isoline for a relatively wide range of in-plane lattice parameters. Again, only at rather extreme values of the epitaxial strain does the system become potentially interesting for piezoelectric applications. For example, at high compressive strains the shear angle shows a steep increase as a function of (see Fig. 7), which results in a relatively large coefficient. (Recall that in similar conditions the out-of-plane parameter is also characterized by a strong dependence on – this yields a large value of the longitudinal coefficient, as we pointed out earlier.) The most interesting region of parameter space in the present context, however, is by far the tensile strain one, where BFO adopts the phase. Here a relatively large value of is combined with a giant value of the out-of-plane dielectric constant, ; this results in a coefficient in excess of pm/V at zero applied field.
Iv Summary
We have studied the physical behaviour of archetypal ferroelectrics PbTiO and BiFeO under arbitrary electrical and mechanical boundary conditions, by means of state-of-the-art ab initio methods. Our calculations provide a detailed description of the complex interplay between polar and structural degrees of freedom in these materials over a comprehensive range of in-plane strain and out-of-plane electric displacement field values. Of particular relevance is the assessment and rationalisation of the related dielectric and piezoelectric response functions, which can be naturally accessed as a by-product of our calculations.
In bulk PTO, we have predicted the existence of a new monoclinic phase that is stable in a narrow strain interval, and that bears important similarities to a previously reported structure in PTO-based superlattices. We have identified the regions in space wherein PTO exhibits a giant dielectric and piezoelectric response. In particular, we have unravelled a coupling mechanism by which its shear piezoelectric coefficient diverges in reaching the monoclinic –tetragonal phase boundary. In BFO, we have shown that a complex interplay between polar and antipolar distortions reigns in the region of high tensile strains wherein an orthorhombic phase becomes stable. In that same region, we find that both the dielectric and shear piezoelectric responses of the system are extremely large.
While the present analysis has been based on bulk calculations, we emphasise that our findings are directly relevant to the the rational design of multicomponent ferroelectric superlattices, as our methodology allows one to accurately predict their functional properties based on the maps of the bulk building blocks. In this sense, our study opens a promising new avenue in the ongoing search for new functional oxide systems with enhanced functionalities.
Acknowledgements.
This research was supported by the Australian Research Council’s Future Fellowship funding scheme (project numbers RG134363 and RG151175), MINECO-Spain (Grants No. MAT2010-18113, No. CSD2007-00041, and No. FIS2013-48668-C2-2-P), and Generalitat de Catalunya (2014 SGR301). We thankfully acknowledge the computer resources, technical expertise and assistance provided by RES, CESGA, and the National Computational Infrastructure supported by the Australian Government.References
- (1) R. E. Cohen and H. Krakahuer, Phys. Rev. B 42, 6416 (1990).
- (2) R. E. Cohen, Nature 358, 136 (1992).
- (3) P. Aguado-Puente, P. García-Fernández, and J. Junquera, Phys. Rev. Lett. 107, 217601 (2011).
- (4) E. Bousquet, M. Dawber, N. Stucki, C. Lichtensteiger, P. Hermet, S. Gariglio, J.-M. Triscone, and P. Ghosez, Nature (London) 452, 732 (2008).
- (5) A. P. Levanyuk and D. G. Sannikov, Sov. Phys. Uspekhi, 17, 199 (1974).
- (6) N. A. Benedek and C. J. Fennie, Phys. Rev. Lett. 106, 107204 (2011)
- (7) N. A. Benedek, A. T. Mulder, and C. J. Fennie, J. Sol. Stat. Chem. 195, 11 (2012).
- (8) J. T. Heron et al., Nature 516, 370 (2014).
- (9) S. Fusil, V. García, A. Barthélémy, and M. Bibes, Annu. Rev. Mater. Res. 44, 91 (2014).
- (10) A. J. Hatt and N. A. Spaldin, Phys. Rev. B 82, 195402 (2010).
- (11) J. M. Rondinelli and N. A. Spaldin, Adv. Mater. 23, 3363 (2011).
- (12) J. M. Rondinelli and S. Coh, Phys. Rev. Lett. 106, 235502(2011).
- (13) S. J. May, J. -W. Kim, J. M. Rondinelli, E. Karapetrova, N. A. Spaldin, A. Bhattacharya, and P. J. Ryan, Phys. Rev. B 82, 014110 (2010).
- (14) C. W. Swartz and X. Wu, Phys. Rev. B 85, 054102 (2012).
- (15) M. Stengel, C. J. Fennie, Ph. Ghosez, Phys. Rev. B 86, 094112 (2012).
- (16) J. Hong and D. Vanderbilt, Phys. Rev. B 87, 064104 (2013).
- (17) M. Dawber, K. M. Rabe, and J. F. Scott, Rev. Mod. Phys. 77, 1083 (2005).
- (18) K. J. Choi et al., Science 306, 5698 (2004).
- (19) H. Ba, B. Dup, S. Fusil, R. Mattana, E. Jacquet, B. Warot-Fonrose, F. Wilhelm, A. Rogalev, S. Petit, V. Cros, A. Anane, F. Petroff, K. Bouzehouane, G. Geneste, B. Dkhil, S. Lisenkov, I. Ponomareva, L. Bellaiche, M. Bibes, and A. Barthlmy, Phys. Rev. Lett. 102, 217603 (2009).
- (20) M. P. Warusawithana et al., Science 324, 367 (2009).
- (21) J. C. Wojdełand J. iguez, Phys. Rev. Lett. 105, 037208 (2010).
- (22) C. Lichtensteiger et al., in Oxides Ultrathin Films: Science and Technology, edited by G. Pacchioni and S. Valeri, Ch. 12, 265 (Wiley-VCH, Germany, 2011).
- (23) P. Ghosez and J. Junquera, J. Comp. Theor. Nanosci. 5, 2071 (2008).
- (24) P. Zubko et al., Nano Letters 12, 2846 (2012).
- (25) J. Sinsheimer et al., Phys. Rev. Lett. 109, 167601 (2012).
- (26) M. Dawber et al., Phys. Rev. Lett. 95, 177601 (2005).
- (27) M. Dawber et al., Adv. Mater. 19, 4153 (2007).
- (28) I. Souza, J. iguez, and D. Vanderbilt, Phys. Rev. Lett. 89, 117602 (2002).
- (29) O. Diéguez and D. Vanderbilt, Phys. Rev. Lett. 96, 056401 (2006).
- (30) M. Stengel, N. A. Spaldin, and D. Vanderbilt, Nature Physics 5, 304 (2009).
- (31) J. C. Wojdeł, P. Hermet, M. P. Ljungberg, P. Ghosez, and J. iguez, J. Phys.:Condens. Matter 25, 305401 (2013). Phys. Rev. Lett. 105, 037208 (2010).
- (32) Y. Yang, W. Ren, M. Stengel, X. H. Yan, and L. Bellaiche, Phys. Rev. Lett. 109, 057602 (2012).
- (33) O. Diéguez, O. E. González-Vázques, J. C. Wojdeł, and J. iguez, Phys. Rev. B 83, 094105 (2011).
- (34) M. Guennou, P. Bouvier, G. S. Chen, B. Dkhil, R. Haumont, G. Garbarino, and J. Kreisel, Phys. Rev. B 84, 174107 (2011).
- (35) S. Prosandeev, D. Wang, W. Ren, J. iguez, and L. Bellaiche, Adv. Funct. Mater. 23, 234 (2013).
- (36) D. Rahmedov, D. Wang, J. iguez, and L. Bellaiche, Phys. Rev. Lett. 109, 037207 (2012).
- (37) C. Cazorla and J. iguez, Phys. Rev. B 88, 214430 (2013).
- (38) R. J. Zeches et al., Science 326, 977 (2009).
- (39) Y. Yang, M. Stengel, W. Ren, X. H. Yan, and L. Bellaiche, Phys. Rev. B 86, 144114 (2012).
- (40) C. Cazorla and M. Stengel, Phys. Rev. B 90, 020101(R) (2014).
- (41) I. A. Kornev, S. Lisenkov, R. Haumont, B. Dkhil, and L. Bellaiche, Phys. Rev. Lett. 99, 227602 (2007).
- (42) M. Stengel, D. Vanderbilt, and N. A. Spaldin, Phys. Rev. B 80, 224110 (2009).
- (43) X. Wu, D. Vanderbilt, and D. R. Hamann, Phys. Rev. B 72, 035105 (2005).
- (44) O. Diéguez, S. Tinte, A. Antons, C. Bungaro, J. B. Neaton, K. M. Rabe, and D. Vanderbilt, Phys. Rev. B 69, 212101 (2004).
- (45) J. L. Blok, D. H. A. Blank, G. Rijnders, K. Rabe and D. Vanderbilt, Phys. Rev. B 84, 205413 (2011).
- (46) A. J. Hatt, N. A. Spaldin, and C. Ederer, Phys. Rev. B 81, 054109 (2010).
- (47) R. Haumont, P. Bouvier, A. Pashkin, K. Rabia, S. Frank, B. Dkhil, W. A. Crichton, C. A. Kuntscher, and J. Kreisel, Phys. Rev. B 79, 184110 (2009).
- (48) N. A. Pertsev, A. G. Zembilgotov, and A. K. Tagantsev, Phys. Rev. Lett. 80, 1988 (1998).