Meta-screening and permanence of polar distortion in metallized ferroelectrics
Ferroelectric materials are characterized by a spontaneous polar distortion. The behavior of such distortions in the presence of free charge is the key to the physics of metallized ferroelectrics in particular, and of structurally-polar metals more generally. Using first-principles simulations, here we show that a polar distortion resists metallization and the attendant suppression of long-range dipolar interactions in the vast majority of a sample of 11 representative ferroelectrics. We identify a meta-screening effect, occurring in the doped compounds as a consequence of the charge rearrangements associated to electrostatic screening, as the main factor determining the survival of a non-centrosymmetric phase. Our findings advance greatly our understanding of the essentials of structurally-polar metals, and offer guidelines on the behavior of ferroelectrics upon field-effect charge injection or proximity to conductive device elements.
In many materials, spontaneous structural distortions occur that break the inversion symmetry of a parent centrosymmetric (CS) structure. These are usually named polar distortions (PDs), since they enable the existence of non-zero polar-vector observables, such as spontaneous electric polarization. Ferroelectrics (FEs) display just such a PD and consequently possess a spontaneous polarization. By definition Lines and Glass (1977), in a FE polarization must be switchable by an external field (non-switchable polarized materials do exist, named pyroelectrics Ambacher et al. (2002)). Because of this requirement, ferroelectrics should be insulators or semiconductors, as opposed to metals, so that they can be acted upon with an external bias. However, it is not a priori obvious that the insulating character itself is necessary for a PD to occur: could it not Anderson and Blount (1965) happen in a metal?
Our general understanding of basic ferroelectric phenomena – largely based on empirical Lines and Glass (1977); Strukov and Levanyuk (1998) and early first-principles Cohen (1992); Posternak et al. (1994); Zhong et al. (1994a); Rabe et al. (2007) studies of perovskite oxides such as BaTiO, PbTiO, or KNbO, centers on the role of electrostatic dipole-dipole couplings as the driving force of the long-range polar order. As a result, free carriers and the attendant electrostatic screening are usually regarded as incompatible with the existence of PDs. Hence, at least among perovskite oxides Benedek and Birol (2016), non-centrosymmetric metals (NCSMs) are usually deemed exotic. This viewpoint has been supported by theoretical work on BaTiO Iwazaki et al. (2012); Wang et al. (2012), whose results seem to be taken as a general rule.
NCSMs are currently a hot topic for obvious reasons of fundamental understanding, but also because of the possible occurrence of quantum phenomena in the context of superconductivity Yildirim (2013); Bauer and Sigrist (2012), and of course their technological relevance to devices involving conductive and FE elements. Indeed, considerable efforts Shi et al. (2013); Filippetti et al. (2016); Puggioni and Rondinelli (2014); Kim et al. (2016) are currently focused on the experimental discovery and first-principles prediction of NCSM compounds, and are yielding experimental Shi et al. (2013), and very recently theoretical Filippetti et al. (2016); Benedek and Birol (2016); He et al. (2016); He and Jin (2016); Shimada et al. (2017), results that question the common wisdom that metallization is incompatible with the occurrence of a PD. For example, first-principles studies have recently suggested that the PD of materials like PbTiO and BiFeO is not strongly affected by the presence of free carriers He et al. (2016); Shimada et al. (2017); He and Jin (2016). Further, some of us took advantage of the chemical origin of ferroelectricity in Bi-based compounds to predict a switchable polar order in BiTiO, a layered perovskite that is metallic Filippetti et al. (2016). A careful examination and rationalization of the compatibility between PDs and free carriers is thus certainly warranted, both to buttress our fundamental understanding and to suggest practical routes to obtain NCSMs, for example by the metallization of a known ferroelectric compound (e.g., by suitable chemical doping or field-effect charge injection).
Here we analyze the effect of doping on PDs by studying from first principles a collection of diverse and representative FE materials. We find that the PD coexists with metallicity in most of the considered compounds. We discuss the atomistic interactions responsible for the observed behaviors, revealing a largely universal meta-screening effect that favors polar distortions upon doping. As a by-product of our work, we come up with obvious prescriptions to obtain FE materials that should yield non-centrosymmetric metals upon doping. Other implications of our results, e.g., as regards hyperferroelectric effects, are also briefly discussed.
Ii Results and discussion
We consider a total of 11 ferroelectric compounds that represent different families owing their FE order to different physical and chemical mechanisms. More specifically, we have LiNbO (LNO), several perovskites (BaTiO or BTO, KNbO or KNO, PbTiO or PTO, BiFeO or BFO, BaMnO or BMO, and BiAlO or BAO), and layered perovskites (LaTiO or LTO227, SrNbO or SNO227, and CaTiO or CTO327), and a (001)-oriented superlattice formed by LaFeO and YFeO perovskite layers that are one unit cell thick (LFO/YFO). Beyond these, we also consider other paraelectric perovskite compounds (LaAlO or LAO), and even metals (Cr and V) and Zintl semiconductors (KSnSb or KSS), to run additional calculations that aid our discussion. Most of our calculations take the ground state structure of these materials, which in all cases is known from the literature, as a starting point to study their behavior upon doping. In a few cases we consider (or identify) additional phases that are stabilized upon doping, and which we introduce in due course. Further details on our calculations are in Appendix A.
ii.1 Polar distortions under doping
We begin by discussing the behavior of PDs in our sample of FE compounds as a function of doping. We adopt the convention that a positive carrier density corresponds to extra electrons (i.e., doping), while negative values indicate hole () doping. We relax all structures as a function of carrier concentration, and monitor the evolution of the PD normalized to its value in the undoped case (see Appendix A for details).
In Fig. 1(a), we present the results obtained under the constraint that the unit-cell volume be fixed and equal to the value obtained in the undoped case. In Fig. 1(b), we show instead the corresponding data when a full volume relaxation is permitted. Figures 1(a) and 1(b) display the same qualitative behavior; the distinction is relevant for reasons to be discussed below.
Figure 1 yields one clear main message: the PDs survives metallization in the vast majority of the considered FE compounds. The PD is unaffected or reinforced in materials in which ferroelectricity is mainly driven by chemical or steric effects (as in PbTiO, BiFeO, BiAlO, and LiNbO), caused by a particular lattice topology or geometry (as in LaTiO and SrTiO López-Pérez and Íñiguez (2011)), or an improper effect triggered by a different primary order parameter (as in CaTiO Oh et al. (2015); Benedek and Fennie (2011) and LaFeO/YFeO Zanolli et al. (2013); Rondinelli and Fennie (2012) superlattices). In fact, in our doping range, the PD disappears only for BTO, BMO, and KNO under doping, and even then, it does take quite some free charge (well above 10 cm) to kill it.
In our description (see also Appendix A) of doping, charge localization, e.g., into narrow gap states, is excluded since we work with perfect crystals, the periodic unit being that of the undoped compound. Hence, the doping charges occupy itinerant Bloch states at the conduction band bottom (electrons) or valence band top (holes), as illustrated by the density of states of BaTiO in Fig. 2, which is representative of all materials.
ii.2 Screening and interactions under doping
To better understand how doping affects the PD, we inspect the effect of the carriers on the relevant interatomic interactions. We specifically analyze the behavior of BTO, BMO, PTO, and BFO, four perovskites that share some similarities, but also present key differences. For example, in both BTO and BMO the PD is mainly driven by the off centering of the B cations, and is known to rely strongly on dipole-dipole interactions Zhong et al. (1994a, 1995); Bhattacharjee et al. (2009). However, Ti has a 3 electronic configuration, while Mn presents a 3 state; hence, the doping electrons and holes occupy different types of orbitals in these two compounds. On the other hand, BFO is a material in which the (very large) PD is driven by the A cation and has a widely accepted chemical origin (Bi’s lone pair) Ravindran et al. (2006); Diéguez et al. (2011). Finally, PTO is a material that shares features of BTO (Ti in a 3 state, with large dipole-dipole interactions) and BFO (Pb’s lone pair).
ii.2.1 BaTiO: raw results
We first focus on BTO, the material where the PD is the least robust of all. To visualize the interactions responsible for the FE instability of BTO, we run the following simulations. We consider the long supercell sketched in Fig. 3(a), which comprises 1120 elemental 5-atom units, with the atoms in their high-symmetry (cubic phase) positions. Then, we displace by 0.05 Å along the Ti atom in the first cell, noting that, because we work with a periodically-repeated supercell, this amounts to creating an array of planes of -polarized dipoles, separated by 19 unit cells (about 76 Å) from each other. Then we compute the forces, considering the undoped case as well as representative doping values. The results are summarized in Figs. 4, 5, and 6.
ii.2.2 Undoped BaTiO
In the undoped case, we find that the force acting on the displaced Ti atom is large and negative. This is a restoring force resulting from two types of interactions: one, short-range repulsive coupling between the Ti and its neighboring oxygens; two, long-range interactions between dipoles within the Å plane, as well as with their periodic images. As sketched in Fig. 3, the lateral interactions between the dipoles in a given plane favor an antipolar order, i.e., they add to the restoring force acting on our displaced Ti. In particular, by performing the corresponding Ewald sum, we estimate this dipole-dipole contribution to be about 0.35 eV/Å in the present case, which is about 25 % of the total force of 1.37 eV/Å obtained in our calculation. (The dominant interactions are those between dipoles in the same plane; the coupling with periodic-image dipole planes is very small.)
If we now move to the two apical oxygens [labeled O(3) in Fig. 3] that lie closest to the displaced Ti, we find relatively large and positive forces acting on them. If we try to understand such forces as the result of short- and long-range interactions, it becomes apparent that they must be dominated by the former kind. Note that the positive dipoles created by the plane of displaced Ti atoms yield a net positive electric field on these O(3) oxygens, which should result in negative dipole-dipole forces. [The relevant dynamical charges are 7.73 for Ti and 6.15 for O(3).] Hence, the computed positive forces must thus be the result of a stronger and repulsive short-range interaction between the Ti and O(3) atoms; this interaction can be seen as tending to preserve an optimal Ti–O(3) distance. Note also that the force computed for the O(3) on the left of the displaced Ti is different from that of the O(3) on the right; this is quite natural, as these two O(3) atoms are not related by symmetry in the distorted configuration; in fact, this difference reflects anharmonic interactions that have an effect even though the considered displacement of the Ti atom (0.05 Å) is relatively small.
As regards the equatorial oxygens [O(1) and O(2)] that are nearest neighbors from the displaced Ti, the obtained positive forces are not a surprise, as both short-range [which will tend to preserve the optimum Ti–O(1) distance in the cubic phase] and long-range [the dipole field in the dipole plane is negative] interactions give a positive contribution. [In this case, the relevant dynamical charge for O(1) and O(2) is about 2.15 .] As regards the Ba atoms, we obtain relative small forces that we do not discuss here.
Interestingly, none of the forces just mentioned, which act on atoms close to the dipole plane, tend to stabilize the polar distortion. Indeed, they are all restoring forces, and it seems safe to interpret them as dominated by short-range (repulsive) couplings favoring the high-symmetry cubic structure. (Short-range interactions are indeed often mentioned in the literature as detrimental to ferroelectricity in BTO Cohen (1992).) However, the situation changes drastically for atoms far from the dipole plane. For those, we obtain finite forces saturating to a non-zero value at around 8 Å from the displaced Ti: in that region, we observe positive forces of about 0.06 and 0.02 eV/Å acting on the Ti and Ba atoms, respectively; and negative forces of about 0.02 eV/Å and 0.05 eV/Å, respectively, acting on the O(1,2) and O(3) anions. Such forces are the result of the quasi-homogeneous field that the dipole planes create in the intermediate region of the supercell; as shown in Supplemental Material (Note 1 and Fig. 1) SupMat (), they can be easily recovered from the potential (Fig. 5) and dynamical charges obtained from our simulations. (By performing the corresponding Ewald sums Allen and Tildesley (1989) for our periodic planes of spaced dipoles, we checked explicitly that, for the situation here considered, a nearly constant field must indeed appear in the intermediate regions. As the separation between dipole planes increases, the field develops small spatial inhomogeneities and eventually decays to zero away from the dipole planes.) These dipole-dipole forces push the cations and anions to move against each other, and thus tend to stabilize a PD. Hence, this is a manifestation of the dipole-dipole interactions responsible for the PD of ferroelectrics like BTO. From a related perspective, since there is no free charge, the equilibrium state of the material should satisfy the Maxwell relation for the electric displacement field . Thus the computed forces in the intermediate regions capture the response of the compound aiming at an homogeneous state of constant when a dipole plane is created.
ii.2.3 Doped BaTiO: electrostatic screening
Let us now discuss the results obtained under doping. One obvious difference with the undoped case is that the forces vanish in the regions away from the dipole plane. Correspondingly, as shown in Fig. 5, the computed potential is flat in those areas. Hence, as expected, the presence of dopants, positive or negative, renders a metallic system and permits the screening of the dipole-dipole interactions. Naturally, this effect goes against the onset of a PD.
We can appreciate how the screening comes about by comparing the DFT results for the non-polar (cubic) and polar (Ti-displaced) structures, as shown in Fig. 6. For example, our results for doping show that an excess of electrons appears in a region within 8 Å to the right of the dipole plane, while an excess of holes occurs in a region of about 12 Å on the left side.
The fact that these two regions are not symmetric makes physical sense: In the cubic structure, the -dopants occupy the Ti- levels, and distribute homogeneously throughout the supercell. Upon displacement of the Ti atom at Å, we essentially have a transfer of mobile electrons from the Ti’s on the left of the dipole plane to the Ti’s on the right side of it. Since the doping level is low, the amount of mobile electrons available in the left-side Ti’s is small, and a relatively large number of atoms are required to provide sufficient charges; in contrast, there are plenty of empty 3 orbitals in the Ti’s on the right, and the excess electrons can be accommodated in a relatively small number of atoms. In the case of -doping [Fig. 6(b)] we observe the same kind of electron depletion (on the left) and accumulation (on the right), and a similarly efficient electrostatic screening (Fig. 4); yet, the details are different, reflecting the different orbitals involved in the charge redistribution. Indeed, in this case the left-side electron donors are O-2 orbitals, and it is also O-2 orbitals that mainly receive electrons on the right.
In accordance with these findings, we observe that electrostatic screening reduces the restoring force on the displaced Ti, as a result of the reduced lateral dipole-dipole interactions within the dipole plane. As Fig. 4 shows, the decrease of the on-site repulsive force is of the order of our ideal estimate of it (i.e., about 0.35 eV/Å). Therefore, in this specific regard, screening favors the occurrence of the polar distortion.
ii.2.4 Doped BaTiO: Short-range effects, meta-screening
Understandably, most discussions of free-carrier effects in the ferroelectrics literature focus on the suppression of the long-range electrostatic interactions. However, our results reveal another important – even dominant – effect in the doped materials, one that is largely independent of the doping type. It is a short-range, screening-related effect that we term meta-screening, which enhances the tendency of the material to display polar distortions.
Compared to the undoped ones, the doped systems exhibit (Figs. 4 and 6) significantly modified forces on atoms close to the dipole planes. These changes happen concurrently with the accumulation of screening electrons and holes (e.g., in the regions marked in Figs. 4 and 6), and follow their variation in width as a function of doping. For atoms in those regions, the forces in the undoped case had an obvious electrostatic character. But, surprisingly, such forces become significantly stronger upon doping, e.g., increasing by a factor of 2, from 0.15 eV/Å to about 0.35 eV/Å for = 0.01 /f.u. on the Ti’s marked with horizontal arrows in Fig. 4(b). Since the dipole-dipole interactions essentially vanish in the doped case, these stronger forces have a different origin, and fall within the general category of short-range interactions. This effect is associated to the electrostatic screening, since it occurs in response to the spatial modulation of the accumulated screening charge (almost irrespective of its sign) around the dipole plane; yet, it clearly transcends the screening of long-range dipolar couplings. We thus term it meta-screening, i.e., occurring along with, but beyond, normal screening.
While a complete discussion of this meta-screening will require further work, its central features lend themselves to simple interpretations. For example, upon doping, the forces acting on the apical O(3) closest to the displaced Ti [marked with arrows in Fig. 4(d)] are positive and significantly smaller than in the undoped case. Hence, it appears that we see in action the repulsive interactions invoked above to rationalize these forces in absence of doping. However, in the doped cases, the accumulation of electrons in the Ti at z6 Å may itself repel the O(3) anion at z4 Å and result in a smaller positive force than in the undoped case; similarly, the accumulation of holes in the Ti at z Å may attract the negatively charged O(3) at and result in relatively small positive force acting on that oxygen. Such considerations apply as well to the forces obtained for the Ti atoms in the immediate vicinity of the dipole plane [marked with horizontal green arrows in Fig. 4(b)]. The one on the right is strongly populated with screening electrons; the obtained positive force would tend to separate it from the displaced Ti, thus expanding the lattice as required to accommodate such an electron excess. The one on the left is in an electron-depleted region, and the obtained positive force would tend to shrink the lattice on that side. Interestingly, this interpretation is consistent with the doping-driven pressure-like effects reported below.
Now, it is important to note that the largest effects observed, especially those pertaining to the Ti atoms closest to the dipole plane, tend to favor the onset of a PD parallel to the imposed dipoles. Indeed, in the accumulation and depletion regions, the computed forces are positive on the cations and negative on the oxygens, and will yield a PD that is qualitatively similar to the FE mode of undoped BTO. It is tempting to interpret the forces obtained under doping as a consequence of imperfect screening, and a signature of how the material tries to reduce the inhomogeneity in the displacement field via a PD. However, as emphasized above, such an electrostatic effect should be strongest in the undoped compound, while we find the largest PD-favoring short-range forces in the doped cases.
Hence, we conclude that the dominant mechanism causing the strongest changes in the short-range forces under doping is a local lattice response accommodating the screening electrons and holes. Incidentally, the similarity between the meta-screening-induced relaxation and BTO’s soft FE mode – both of which are essentially characterized by the relative displacement of Ti-O(3) pairs – is not surprising: upon a local perturbation (i.e., our imposed dipole planes), the lattice response will typically be dominated by the lowest-energy distortions that become activated by the perturbation; in our case, such distortions are the soft polar modes, which continue to be rather low in energy in BTO even upon doping [this is obvious from Fig. 7(a), discussed below].
In summary, we have evidence for a previously unnoticed, short-range meta-screening effect, which is a by-product of the electronic screening and favors polar distortions for both and doping. As shown below, meta-screening occurs in all the considered perovskite oxides, hence it is likely to be a general phenomenon.
ii.2.5 Soft modes under doping
To address the (in)stability of cubic BTO against polar distortions and its dependence on doping, we compute the force-constant matrix at the point (Brillouin zone center) via standard finite-displacement methods in our 1120 supercell. We focus on the -polarized instability, and displace the atoms by 0.01 Å from their ideal cubic positions. The -point force-constant matrix is trivially derived from the computed forces by a supercell average. While the same -point matrix can be easily obtained in the five-atom BTO unit cell, using the long supercell we can monitor the various interactions in real space, and modify them by hand to test their individual effects. Note also that this force-constant matrix yields the zone-center dynamical matrix just by introducing suitable mass factors. Any soft-mode instability of the cubic structure results in both matrices having (at least) one negative eigenvalue, corresponding to a negative force constant (energy curvature) in the former case, and to an imaginary frequency in the latter.
Figure 7(a) shows our basic result, i.e., the evolution of the force constant (or stiffness) of the soft polar mode, , as a function of doping. As expected, we find that electron doping eliminates the polar instability at /f.u., which roughly agrees with the results in Fig. 1. (Slight quantitative differences are due to volume effects, because in Fig. 7 we work with the optimized undoped cubic cell, while in Fig. 1 we optimize the cell of the polar structure.) In contrast, the polar instability survives when the doping is with holes. Let us stress that our supercell calculations only involve displacements of atoms in the unit cell at the origin, so the settings are identical (except for the use of smaller displacements, to make sure we are in the harmonic regime) to those used in the dipole-plane simulations described above. Hence, all the electronic effects discussed earlier in this paper are obviously active in the simulations, and contribute to the obtained evolution of .
We have seen above that the long-range dipole-dipole interactions – well, established to be the driving force for ferroelectricity in undoped BTO, are all but gone as soon as some dopants are introduced in the material. It is thus surprising that doped BTO retains a polar soft mode in some doping ranges. Incomplete electrostatic screening might be a tempting explanation for the case of small doping, but it most certainly does not apply to the results for large doping. Instead, it seems more reasonable to turn our attention to the meta-screening effects revealed above as a possible origin for the observed behavior. Let us focus on the most obvious one, i.e., the strong coupling between first-nearest-neighboring Ti atoms that renders the very large forces marked with green horizontal arrows in Fig. 4(b). To test whether such an interaction may explain the polar instability in doped BTO, we run the following computational experiment.
The -point force-constant matrix and the soft polar mode obtained from its diagonalization satisfy
where and run over the atoms in the unit cell and spatial directions, and is the soft-mode force constant, depicted in Fig. 7(a). Naturally, all these quantities depend implicitly on . We now test how the stiffness constant of the soft mode changes if we modify some key interactions. To do this, we construct a new force-constant matrix that is identical to except that we impose the coupling between first-nearest-neighboring Ti atoms be always that of the undoped case, independently of the doping level. We thus remove the most prominent meta-screening effect revealed above. The modified stiffness
is shown as function of doping in Fig. 7(a) (dashed red lines). It is obvious that once the meta-screening effect is removed, BTO instantly loses its polar instability upon doping, irrespective of the sign of the extra charges. Hence, the meta-screening effect is the driving force for the polar instability of doped BTO.
Note that in the past, e.g., in the important work of Wang et al. Wang et al. (2012), short-range forces have generally been assumed to be independent of doping. Based on this (incorrect) assumption, it is most natural to attribute the persistence of the PD in metallized BTO to the action of screened, but strong enough, Coulomb interactions. Our present results clearly show that this is not the case.
There is a clear - asymmetry in Fig. 7(a), evidenced e.g. by the slope discontinuity of around . This is a direct consequence of the existence of a band gap in the material, and of the different character of the states occupied by the doping electrons (Ti’s 3) and holes (O’s 2). Further, while the meta-screening effect is sufficient to preserve the polar instability in doped BTO in this range, it is overcome by some other interaction in the doped compound, where the PD eventually disappears ( for /f.u.). The largest and most relevant differences between and doping do not pertain to electrostatic screening, which is very efficient in both cases and causes similar meta-screening effects. Instead, the greatest differences pertain to the shortest-range interactions; most importantly, the results in Fig. 4 show that the restoring forces are systematically weaker for -doping.
This result can be understood by recalling the usual picture of the Ti–O electronic hybridizations in BTO, which emphasizes the key role of second-order Jahn-Teller effects to permit the FE distortion of this material. In essence, the energy of the compound can be reduced by the hybridization of (empty) Ti-3 and (occupied) O-2 states, which is prompted by the onset of the PD and associated reduction of the Ti–O(3) distance. Additional electrons would tend to occupy the empty orbitals above the band gap, and thus increase the energy significantly; in contrast, additional holes would occupy filled valence states, and result in a relatively moderate energy increase. Hence, it naturally follows that short-range restoring (repulsive) forces will be stronger for the -doping case, which is consistent with the observed suppression of the PD only upon electron doping.
To test the effect of these different forces, we run another computational experiment along the lines of the one just described. We construct modified force-constant matrices in the following way: For a certain -doping (-doping) given by , we substitute the self-interaction of the Ti atom [responsible for the largest restoring force, marked with a gray arrow in Fig. 4(b)] by the value obtained for the corresponding doping ( doping). We thus obtain a second modified stiffness ; the results are in Fig. 7(a), green dotted lines. We observe a notable degradation of the polar instability under -doping, and a sizable strengthening upon -doping. (The irregular behavior of near reflects the qualitatively different effects of - and -doping on the short-range interactions, due to the band gap. Similarly, the occurrence of a minimum of for is a by-product of the artificial way in which we construct , and not worth discussing.) These results thus indicate that the main difference between electron and hole doping lies in their effect on the short-range repulsive couplings.
ii.2.6 Other materials
Having discussed in detail BTO’s case, our findings for BMO, PTO and BFO are easy to present. Figure 8 summarizes the results from our supercell simulations with imposed dipole planes, which we create by displacing Ti and Pb atoms in the case of PTO [Figs. 8(b) and (c), respectively], Bi atoms in the case of BFO [Fig. 8(d)], and Mn atoms in the case of BMO [Fig. 8(e)]. We also include in Fig. 8(a) the results for BTO, for an easier comparison. Remarkably, the computed forces exhibit the same essential features discussed above for BTO.
Most importantly, we emphasize that meta-screening, i.e. the enhancement of short-range interactions upon doping, occurs in all the considered materials, and is thus very likely to be a general phenomenon. Moreover, in all cases, meta-screening favors again polar distortions. (Figure 8 shows positive forces on the key cations; the forces on the oxygens, not shown here, are negative.)
To drive this point home, we show in Fig. 7(b) three versions of the stiffness constant of the soft mode of PTO as a function of doping. Similarly to BTO, we present the stiffness obtained from the -point force-constant matrix, along with two other quantities: one is (Pb), obtained from Eq. (2) for the same matrix, except for the strongest meta-screening forces acting on Pb ions [marked with arrows in Fig. 8(c)] being replaced by the corresponding values in the undoped case. If we also similarly modify the forces acting on the Ti ions [marked with arrows in Fig. 8(b)] we obtain by the same procedure a third stiffness variant, (Pb&Ti). Essentially, when the system is purged of the meta-screening couplings, the soft modes are much less soft, i.e., their force constants are much less negative, in accordance with our previous conclusion that meta-screening is the main driver of the permanence of PDs in doped FEs.
We should note that, from the evidence at hand, we cannot tell whether the meta-screening mechanism is a necessary condition for the PD to occur in a compound like PTO. To elucidate that question, we would need an accurate quantification of the meta-screening contribution to the forces, so that such effects can be clearly disentangled from other (steric/chemical) factors. This poses an interesting and non-trivial challenge to electronic-structure theory, and remains for future work.
The results in Fig. 8 offer other interesting insights. For example, it is apparent that the restoring forces are relatively small for the Pb (in PTO) and Bi (in BFO) cations, and relatively large for Ti (in both BTO and PTO) and Mn (in BMO). We think this difference can be partly attributed to the stereochemical activity of Pb and Bi’s lone pairs, which tends to compensate the electronic repulsion between ionic cores.
It is also interesting to note that the restoring force acting on the displaced Mn (3) cation in BMO is significantly smaller than that on displaced Ti (3) cation in both BTO and PTO. This may seem at odds with the usual view that empty 3 orbitals are indispensable for B-site driven ferroelectricity to occur. Yet, one should note that, as regards the possibility that a Mn cation in an O environment drives ferroelectricity, the most relevant 3 orbitals are those with symmetry, which are directed towards the oxygen anions and are empty in this case. Hence, ferroelectricity in BMO should not be penalized by strong repulsive forces associated to the Mn-3 configuration Filippetti and Hill (2002); Bhattacharjee et al. (2009). Having said this, to explain why the restoring forces acting on BMO’s Mn cation are significantly smaller than those obtained for BTO’s Ti, we probably should resort to simple steric arguments. Indeed, the ionic radii of Ti and Mn in an octahedral O environment are 0.605 Å and 0.53 Å, respectively Shannon (1976); then, noting that BTO and BMO share the same A-site cation, size considerations suggest that it will be easier for the smaller Mn to move off-center, which is clearly consistent with the relatively weak restoring force obtained in our calculations.
Finally, let us remark the striking similarity between our results for the Ti forces in BTO [Fig. 8(a)] and the corresponding ones in PTO [Fig. 8(b)]; this suggests that interactions between same atom pairs are relatively unaffected by the different chemical environment in different perovskite oxides, an observation that is in line with previous first-principles studies Ghosez et al. (1999). Additionally, note that the results for the Pb forces in PTO [Fig. 8(c)] and the Bi forces in BFO [Fig. 8(d)] are quite similar as well. While we do not want to overinterpret these observations, they are clearly suggestive of the hybrid nature of ferroelectricity in PTO, as the polar soft mode of this material is obviously participated by both the A and B cationic sublattices; in contrast, BFO and BTO are textbook examples of compounds in which ferroelectricity is driven by only one cation sublattice, respectively A and B.
ii.3 Additional remarks
ii.3.1 Volume changes and transitions under doping
As shown in Fig. 9(a), our simulations yield a universal behavior regarding the volume of the doped materials: additional electrons cause an expansion, while additional holes cause a contraction. Such an effect had already been observed in the past, in independent investigations of BaTiO Iwazaki et al. (2012), BiFeO He et al. (2016) and PbTiO He and Jin (2016). Our present work confirms this behavior and shows that it pertains to all the diverse ferroelectrics here considered.
One may wonder whether this volume effect has any influence on the survival, or disappearance, of the PD upon doping. To check this, in Fig. 1 we compared the results obtained for constant volume [Fig. 1(a)] and relaxed volume [Fig. 1(b)], noting that in the considered doping range the volume changes can be up to %. Our results show that FEs conserve their PD irrespective of whether we allow the volume to relax or not (with the partial exception of -doped BTO, KNO and BMO). This suggests that the effects discussed above, responsible for the disappearance (screening) or survival (meta-screening) of the PD, are not much affected by even fairly substantial volume changes.
Naturally, we do find some differences when volume relaxation is allowed. For example, it is apparent that the contraction associated to doping is detrimental to the PD of BTO and BMO. This result lends itself to a simple interpretation, as it is well-known that a compression tends to weaken ferroelectricity in conventional perovskite oxides like BTO Ishidate et al. (1997); Íñiguez and Vanderbilt (2002).
As emphasized by other authors He et al. (2016), the doping-driven volume changes operate in essentially the same way as a hydrostatic pressure would, and can potentially induce structural phase transitions beyond those (polar to non-polar) discussed above. As an example, in Fig. 9 we show the behavior of PTO under doping and under a negative pressure [Figs. 9(b) and (c), respectively]. In both cases, the volume increase causes a transition into a so-called super-tetragonal phase with giant aspect ratio Tinte et al. (2003); Wang et al. (2015). The analogy between doping and pressure is further ratified by our studies of BiFeO and LaAlO [see Supplemental material (Note 2 and Figs. 2–4) SupMat ()], and suggests that non-trivial structural effects may occur, to some extent at least, whenever dopants stay spatially delocalized.
We can try to rationalize the volume changes in terms of the bonding/anti-bonding character of the electronic states affected by the doping. As described in the Supplemental Material (Note 3 and Figs. 5–12) SupMat (), some of our results are straightforwardly interpreted (e.g., dopants occupy anti-bonding states in our insulating oxides, which suggests a lattice expansion consistent with our calculations), and others can be explained by invoking plausible second-order orbital mixing effects. Yet, we also find examples (in particular, for the non-oxidic materials V, Cr and KSS) where such bonding arguments clearly fail, which questions their general validity. We are thus inclined to believe that the obtained volume effects may be the consequence of a rather crude steric mechanism of sorts (grossly speaking: electrons do occupy space), which prevails over the bonding characteristics of the (de)populated states.
We also note that our way of simulating doping is not expected to reproduce polarons. Since previous work suggests that in some cases volume changes are suppressed when chemical dopants Iwazaki et al. (2012) or self-trapped electrons and holes He et al. (2016) are considered explicitly, the doping-driven volume changes just reported should be considered realistic insofar as the free charges remain extended. Since localization is frequent in oxides, our volume changes may be considered an upper limit when compared with experiment, but should apply fairly closely when the injected charge is delocalized, as at metal/ferroelectric interfaces (where some charge spillage always occurs) and in the case of field-effect injection or electrostatic doping.
Hyperferroelectric compounds Garrity et al. (2014) are soft-mode ferroelectrics whose paraelectric phase displays an unstable longitudinal-optical (LO) polar phonon band. To obtain such an exotic property, which suggests, e.g., that an hyperferroelectric can form (meta)stable FE domain walls that would be formally charged, it is mandatory to have unstable transversal-optical (TO) polar phonons and a relatively small LO-TO splitting. The latter is typical of materials with large high-frequency dielectric permittivity , i.e., materials with a very efficient electrostatic screening. Hence, whenever we have a hyperferroelectric that displays regular (TO) FE instabilities in spite of weak dipole-dipole interactions, that is a good candidate to remain polar when such couplings are totally screened ( diverges upon doping). Conversely, materials that remain polar upon metallization may in principle be good candidates for hyperferroelectricity.
To investigate this connection, we looked for hyperferroelectricity in a subset of our considered FE materials, by running straightforward phonon and perturbative calculations that allow us to compute the LO-TO splitting [see details in Supplemental Material (Note 4 and Table I) SupMat ()]. To our surprise, we find that only four compounds (LNO, LTO227, SNO227 and CTO327) are hyperferroelectric, while most of the materials displaying a strong and robust PD upon doping are not. Indeed, in materials like PTO and BFO, while the zone-center (TO) polar instability of the cubic phase is very strong, the LO-TO splitting is even stronger, yielding a stable LO band. Note that the very large LO-TO splitting that is typical of FE perovskite oxides can be traced back to the anomalously large polarity of the soft modes (which in turn reflects unusually large dynamical charges Zhong et al. (1994b)) and their relatively small .
In conclusion, our first-principles study of diverse ferroelectrics shows that their characteristic polar distortion is generally stable upon charge doping. Remarkably, our results reveal a previously unnoticed meta-screening effect that is essential to the permanence of the non-centrosymmetric phase. This seemingly-universal meta-screening mechanism is triggered by the rearrangement of mobile electrons and holes associated to the screening of dipolar interactions, is essentially independent of the sign of the doping charges, and results in short-range couplings favoring a polar distortion. Our results thus provide unprecedented insight into the behavior of metallized ferroelectrics, potential implications ranging from the discovery of new polar metals to the design of metal/ferroelectric interfaces or charge-injection effects in these compounds.
Acknowledgements.Work supported by the Luxembourg National Research Fund through Grants P12/4853155 COFERMAT (H.J.Z. and J.I.), INTER/MOBILITY/15/9890527 GREENOX (L.B., H.J.Z. and J.I.), and AFR Grant No. 9934186 (C.E.S.). Additionally, V.F. was supported by Progetto biennale di ateneo UniCA/FdS/RAS 2016; E.C. by the Spanish MINECO through the Severo Ochoa Centers of Excellence Program under Grant SEV-2015-0496, as well as through Grant FIS2015-64886-C5-4-P, and by Generalitat de Catalunya (2017SGR1506); and L.B. by the Air Force Office of Scientific Research under Grant No. FA9550-16-1-0065. Computational resources have been provided by the PRACE-3IP DECI-13 grant 13DECI0270 INTERPHON (Salomon cluster at the Czech National Supercomputing Center), the CRS4 Computing Center (Piscina Manna, Pula, Italy), and the Arkansas High Performance Computing Center. Also, we are grateful to M. Stengel and P. Zubko for fruitful discussions.
Appendix A Methods
We use density functional theory (DFT) within the generalized gradient approximation (PBEsol functional Perdew et al. (2008)) as implemented in the software package VASP Kresse and Furthmüller (1996); Kresse and Joubert (1999). For all considered compounds, the electronic wave functions are represented in a basis of plane waves truncated at 500 eV. Reciprocal space integrals are computed using -point grids that are equivalent to (or denser than) a sampling of the Brillouin zone of an elemental five-atom perovskite cell. The interaction between ionic cores and electrons is treated within the so-called plane augmented wave (PAW) approach Blöchl (1994), solving explicitly for the following electrons: O’s 2 and 2; Li’s 2; K’s 3, 3, and 4; Ba’s 5, 5 and 6; Pb’s 6 and 6; Ca’s 3 and 4; Sr’s 4, 4, and 5; Bi’s 6 and 6; La’s 5, 5, 5, and 6; Y’s 4, 4, 4 and 5; Al’s 3 and 3; Ti’s 3 and 4; Mn’s 3 and 4; Fe’s 3 and 4; Nb’s 4, 4, 4, and 5; Sn’s 5 and 5; Sb’s 5 and 5; Cr’s 3 and 4; and V’s 3 and 4. For Fe’s 3 electrons we use the “Hubbard correction” introduced by Dudarev et al. Dudarev et al. (1998) with eV; for Mn’s 3 electrons we use the correction introduced by Liechtenstein et al. Liechtenstein et al. (1995) with eV and eV. (In the case of BiFeO, we explicitly verified that our results for the persistence of the PD upon doping remain essentially the same for values between 3 and 5 eV.) Structural relaxations are run until residual forces and stresses fall below 0.005 eV/Å and 0.05 GPa, respectively. These calculations conditions were checked to render sufficiently converged results.
We simulate the effect of doping by varying the number of electrons in the cell, and adding a neutralizing homogeneous charge background. This approach, the standard one employed in most of the previous works on this problem Iwazaki et al. (2012); Wang et al. (2012); Benedek and Birol (2016); He et al. (2016); He and Jin (2016), does not describe the doping species explicitly, which greatly simplifies the calculation. Further, we use the smallest cells describing the equilibrium structures of the undoped material, namely a 5-atom cell for perovskites like BaTiO and PbTiO, a 10-atom cell for a material like BiFeO, etc. Such settings impose restrictions on the possible arrangements of added electrons or holes, such as for example polaron states (we note in passing that standard semi-local density functional methods are a priori not expected to yield stable states of that type). We thus expect that our simulations will tend to exaggerate the tendency towards metallization and the effectiveness of doping in producing screening, as well as in modifying the structure. Nevertheless, as evidenced by the results here reported, these idealized conditions are relevant to better understand the intrinsic response of FE materials to carrier doping. On the other hand, our results are directly relevant to situations that are typical of ferroelectric nanostructures, e.g., whenever the ferroelectric material is partly metallized near the interface with an electrode, or extra carriers are injected by electrostatic doping, etc.
For the ferrites (BiFeO and LaFeO/YFeO) and manganite (BaMnO), we use the well-known lowest-energy spin arrangement (anti-ferromagnetic with anti-parallel nearest-neighboring spins) and the standard scalar-magnetism (collinear) approximation. Note that, according to previous studies Wojdeł and Íñiguez (2009, 2010), non-collinear magnetism and spin-orbit interactions are expected to have a negligible impact on the FE instabilities of these compounds; hence, we do not consider them here.
We use standard analysis tools to study the doping-induced effects. In particular, we use the FINDSYM Stokes and Hatch (2005) and AMPLIMODES Orobengoa et al. (2009); Perez-Mato et al. (2010) codes to determine the space group of our doped structures and to calculate the mode-resolved distortion amplitudes, respectively. When computing the distortion amplitudes with AMPLIMODES, the undoped high-symmetry phase ( for simple perovskites, for layered perovskite CaTiO, for layered perovskites LaTiO and SrNbO, and for superlattice LaFeO/YFeO) is taken as the reference structure. Note that when AMPLIMODES compares a reference CS structure with a polar one (doped or undoped), it will in general yield a collection of amplitudes corresponding to modes of different symmetries; from those, we retain the result corresponding to the polar mode (which e.g. corresponds to the irreducible representation in the case of simple perovskites) to quantify the CS-breaking distortion.
Finally, we also use the ASE tools Larsen et al. (2017); Hunter (2007) and VESTA Momma and Izumi (2011) for analysis and visualization of our results, as well as the LOBSTER code Dronskowski and Bloechl (1993); Deringer et al. (2011); Maintz et al. (2013, 2016a, 2016b) to characterize the bonds and electronic structure via a standard COHP (crystal orbital hamilton population) analysis.
- Lines and Glass (1977) M. E. Lines and A. M. Glass, Principles and Applications of Ferroelectrics and Related Materials, Oxford Classic Texts in the Physical Sciences (Clarendon Press, Oxford, 1977).
- Ambacher et al. (2002) O. Ambacher, J. Majewski, C. Miskys, A. Link, M. Hermann, M. Eickhoff, M. Stutzmann, F. Bernardini, V. Fiorentini, V. Tilak, B. Schaff, and L. F. Eastman, J. Phys.: Condens. Matter 14, 3399 (2002).
- Anderson and Blount (1965) P. W. Anderson and E. I. Blount, Physical Review Letters 14, 217 (1965).
- Strukov and Levanyuk (1998) B. Strukov and A. Levanyuk, Ferroelectric Phenomena in Crystals: Physical Foundations (Springer, 1998).
- Cohen (1992) R. E. Cohen, Nature 358, 136 (1992).
- Posternak et al. (1994) M. Posternak, R. Resta, and A. Baldereschi, Physical Review B 50, 8911 (1994).
- Zhong et al. (1994a) W. Zhong, D. Vanderbilt, and K. M. Rabe, Physical Review Letters 73, 1861 (1994a).
- Rabe et al. (2007) K. M. Rabe, C. H. Ahn, and J. Triscone, eds., Physics of Ferroelectrics: A Modern Perspective (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007).
- Benedek and Birol (2016) N. A. Benedek and T. Birol, J. Mater. Chem. C 4, 4000 (2016).
- Iwazaki et al. (2012) Y. Iwazaki, T. Suzuki, Y. Mizuno, and S. Tsuneyuki, Physical Review B 86, 214103 (2012).
- Wang et al. (2012) Y. Wang, X. Liu, J. D. Burton, S. S. Jaswal, and E. Y. Tsymbal, Physical Review Letters 109, 247601 (2012).
- Yildirim (2013) T. Yildirim, Physical Review B 87, 020506 (2013).
- Bauer and Sigrist (2012) E. Bauer and M. Sigrist, eds., Non-Centrosymmetric Superconductors: Introduction and Overview (Springer, 2012).
- Shi et al. (2013) Y. Shi, Y. Guo, X. Wang, A. J. Princep, D. Khalyavin, P. Manuel, Y. Michiue, A. Sato, K. Tsuda, S. Yu, M. Arai, Y. Shirako, M. Akaogi, N. Wang, K. Yamaura, and A. T. Boothroyd, Nature Materials 12, 1024 (2013).
- Filippetti et al. (2016) A. Filippetti, V. Fiorentini, F. Ricci, P. Delugas, and J. Íñiguez, Nature Communications 7, 11211 (2016).
- Puggioni and Rondinelli (2014) D. Puggioni and J. M. Rondinelli, Nature Communications 5, 3432 (2014).
- Kim et al. (2016) T. H. Kim, D. Puggioni, Y. Yuan, L. Xie, H. Zhou, N. Campbell, P. J. Ryan, Y. Choi, J.-W. Kim, J. R. Patzner, S. Ryu, J. P. Podkaminer, J. Irwin, Y. Ma, C. J. Fennie, M. S. Rzchowski, X. Q. Pan, V. Gopalan, J. M. Rondinelli, and C. B. Eom, Nature 533, 68 (2016).
- He et al. (2016) X. He, K.-J. Jin, H.-z. Guo, and C. Ge, Physical Review B 93, 174110 (2016).
- He and Jin (2016) X. He and K.-J. Jin, Physical Review B 94, 224107 (2016).
- Shimada et al. (2017) T. Shimada, T. Xu, Y. Araki, J. Wang, and T. Kitamura, Advanced Electronic Materials 3, 1700134 (2017).
- López-Pérez and Íñiguez (2011) J. López-Pérez and J. Íñiguez, Physical Review B 84, 075121 (2011).
- Oh et al. (2015) Y. S. Oh, X. Luo, F.-T. Huang, Y. Wang, and S.-W. Cheong, Nature Materials 14, 407 (2015).
- Benedek and Fennie (2011) N. A. Benedek and C. J. Fennie, Physical Review Letters 106, 107204 (2011).
- Zanolli et al. (2013) Z. Zanolli, J. C. Wojdeł, J. Íñiguez, and P. Ghosez, Physical Review B 88, 060102 (2013).
- Rondinelli and Fennie (2012) J. M. Rondinelli and C. J. Fennie, Advanced Materials 24, 1961 (2012).
- Zhong et al. (1995) W. Zhong, D. Vanderbilt, and K. M. Rabe, Physical Review B 52, 6301 (1995).
- Bhattacharjee et al. (2009) S. Bhattacharjee, E. Bousquet, and P. Ghosez, Physical Review Letters 102, 117602 (2009).
- Ravindran et al. (2006) P. Ravindran, R. Vidya, A. Kjekshus, H. Fjellvåg, and O. Eriksson, Physical Review B 74, 224412 (2006).
- Diéguez et al. (2011) O. Diéguez, O. E. González-Vázquez, J. C. Wojdeł, and J. Íñiguez, Physical Review B 83, 094105 (2011).
- (30) See Supplemental Material for additional details on various aspects of this work.
- Allen and Tildesley (1989) P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications (Clarendon Press, 1989).
- Filippetti and Hill (2002) A. Filippetti and N. A. Hill, Physical Review B 65, 195120 (2002).
- Shannon (1976) R. D. Shannon, Acta Crystallographica Section A 32, 751 (1976).
- Ghosez et al. (1999) P. Ghosez, E. Cockayne, U. V. Waghmare, and K. M. Rabe, Physical Review B 60, 836 (1999).
- Ishidate et al. (1997) T. Ishidate, S. Abe, H. Takahashi, and N. Môri, Physical Review Letters 78, 2397 (1997).
- Íñiguez and Vanderbilt (2002) J. Íñiguez and D. Vanderbilt, Physical Review Letters 89, 115503 (2002).
- Tinte et al. (2003) S. Tinte, K. M. Rabe, and D. Vanderbilt, Physical Review B 68, 144105 (2003).
- Wang et al. (2015) J. Wang, B. W. Eerd, T. Sluka, C. Sandu, M. Cantoni, X.-K. Wei, A. Kvasov, L. J. McGilly, P. Gemeiner, B. Dkhil, A. Tagantsev, J. Trodahl, and N. Setter, Nature Materials 14, 985 (2015).
- Garrity et al. (2014) K. F. Garrity, K. M. Rabe, and D. Vanderbilt, Physical Review Letters 112, 127601 (2014).
- Zhong et al. (1994b) W. Zhong, R. D. King-Smith, and D. Vanderbilt, Physical Review Letters 72, 3618 (1994b).
- Perdew et al. (2008) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Physical Review Letters 100, 136406 (2008).
- Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
- Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
- Blöchl (1994) P. Blöchl, Physical Review B 50, 17953 (1994).
- Dudarev et al. (1998) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Physical Review B 57, 1505 (1998).
- Liechtenstein et al. (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Physical Review B 52, R5467 (1995).
- Wojdeł and Íñiguez (2009) J. C. Wojdeł and J. Íñiguez, Physical Review Letters 103, 267205 (2009).
- Wojdeł and Íñiguez (2010) J. C. Wojdeł and J. Íñiguez, Physical Review Letters 105, 037208 (2010).
- Stokes and Hatch (2005) H. T. Stokes and D. M. Hatch, Journal of Applied Crystallography 38, 237 (2005).
- Orobengoa et al. (2009) D. Orobengoa, C. Capillas, M. I. Aroyo, and J. M. Perez-Mato, Journal of Applied Crystallography 42, 820 (2009).
- Perez-Mato et al. (2010) J. M. Perez-Mato, D. Orobengoa, and M. I. Aroyo, Acta Crystallographica A66, 558 (2010).
- Larsen et al. (2017) A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. DuÅak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. SchiÃ¸tz, O. SchÃ¼tt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, Journal of Physics: Condensed Matter 29, 273002 (2017).
- Hunter (2007) J. D. Hunter, Computing In Science & Engineering 9, 90 (2007).
- Momma and Izumi (2011) K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).
- Dronskowski and Bloechl (1993) R. Dronskowski and P. E. Bloechl, The Journal of Physical Chemistry 97, 8617 (1993).
- Deringer et al. (2011) V. L. Deringer, A. L. TchougrÃ©ef, and R. Dronskowski, The Journal of Physical Chemistry A 115, 5461 (2011).
- Maintz et al. (2013) S. Maintz, V. L. Deringer, A. L. TchougrÃ©eff, and R. Dronskowski, Journal of Computational Chemistry 34, 2557 (2013).
- Maintz et al. (2016a) S. Maintz, V. L. Deringer, A. L. TchougrÃ©eff, and R. Dronskowski, Journal of Computational Chemistry 37, 1030 (2016a).
- Maintz et al. (2016b) S. Maintz, M. Esser, and R. Dronskowski, Acta Phys. Pol. B 47, 1165 (2016b).