Ferroelectricity at ferroelectric domain walls
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J un Ferroelectricity at ferroelectric domain walls
Jacek C. Wojde l and Jorge ´I˜niguez
Institut de Ci`encia de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain
We present a first-principles study of model domain walls (DWs) in prototypic ferroelectricPbTiO . At high temperature the DW structure is somewhat trivial, with atoms occupying high-symmetry positions. However, upon cooling the DW undergoes a symmetry-breaking transitioncharacterized by a giant dielectric anomaly and the onset of a large and switchable polarization.Our results thus corroborate previous arguments for the occurrence of ferroic orders at structuralDWs, providing a detailed atomistic picture of a temperature-driven DW-confined transformation.Beyond its relevance to the field of ferroelectrics, our results highlight the interest of these DWs inthe broader areas of low-dimensional physics and phase transitions in strongly-fluctuating systems. PACS numbers: 77.80.B-, 63.70.+h, 71.15.Mb
The structural domain walls (DWs) occurring in ferro-electric (FE) and ferroelastic (FS) materials have becomea focus of attention. Recent studies show that the DWscan present a variety of properties, from conductive [1–4] and optical [5, 6] to magnetic [7–9], that differ fromthose of the neighboring domains, which suggests thatthey could be the active element in nano-technologicalapplications [10, 11]. Elucidating the DW behavior posesmajor experimental challenges, and the origin of most ofthe newly discovered effects remains unclear. In fact,we still lack a detailed structural and dynamical pic-ture of the DWs, and in many cases we can only specu-late about the structure–property relationships at workwithin them. Hence, there is a pressing need for predic-tive theoretical studies tackling the DWs at an atomisticlevel and at the relevant conditions of temperature, etc.The DW structure, and even the possible occurrenceof DW-confined ferroic orders, have been discussed theo-retically for decades, usually in the framework of contin-uum Ginzburg-Landau or phenomenological model theo-ries [12–22]. Materials with competing structural insta-bilities have been a focus of attention, a good examplebeing perovskite SrTiO (STO). STO undergoes a FStransition driven by an anti-ferrodistortive (AFD) modethat involves concerted rotations of the O octahedra inthe perovskite structure. This mode competes with a FEinstability that is suppressed by the onset of the AFDdistortion [23]. Yet, there are both theoretical and ex-perimental indications that a polar order occurs at lowtemperatures within STO’s FS DWs [16, 24–26], i.e., inthe region where the otherwise dominant AFD distor-tions vanish. In this context, it is worth noting recentfirst-principles studies predicting that PbTiO (PTO)[27] and related compounds [28] present a FE-AFD com-petition that is even stronger than the one occurring inSTO. These are the ideal conditions to obtain interestingeffects at structural DWs, and motivated this work. Low-temperature study .– We employed the tools ofRef. 27, which permit large-scale simulations with first-principles predictive power, to investigate an ideal ver-sion of the simplest DWs occurring in PTO, namely, 180 ◦ boundaries separating regions of opposed polariza-tion and being perfectly planar. We used the model po-tential for PTO labeled “ L I ” [27], which we briefly de-scribe in [29]. As shown in Fig. 6(a), we set the polariza-tion of the first domain P I parallel to the [100] directionof the perovskite lattice, and took P II k [¯100] for the sec-ond one; the DW in between was assumed to reside ina (001) plane. Our supercell, which contains 12 × × x -component of the polarization ( P x ) as we move along z . We observe two domains within which PTO adoptsthe structure of its homogeneous ground state, with anassociated polarization of about 0.99 C/m and a cell as-pect ratio of about 1.07. The domains are separated by aDW centered at a PbO plane and presenting a thicknessof about one unit cell.Our DWs do not display any rotations of the O oc-tahedra. This result lends itself to a simple explana-tion: Because the DWs are ultra-thin, hypothetical DW-localized AFD modes would overlap with the neighboringFE distortions, and thus be penalized by the FE-AFDcompetition. As a result, the absence of localized AFDmodes seems rather natural.Nevertheless, the structure of the DWs is far from be-ing trivial. As shown in Fig. 6(c), a non-zero P y polariza-tion appears at the DW plane and rapidly vanishes as wemove into the domains. This DW polarization is switch-able, as evidenced by the hysteresis loop in Fig. 2(a).Further calculations show that the polarizations of neigh-boring DWs couple, and tend to align in a parallel con-figuration when the walls are sufficiently close. How-ever, the DW–DW interaction quickly decreases with theseparation distance; for example, for our 12 × ×
20 su-percell, the energy split between the parallel and anti-parallel states is about 0.01 meV per
DW cell, which isnegligible. Hence, the anti-parallel configuration shownin Fig. 6 is a stable state (i.e., a local minimum of the xyz (a) DW {{ (d) xy z Pb Ti O domain I domain II d o m a i n I d o m a i n II DW (b)(c) PbO planes X TiO planes domain I domain II DW -0.8-0.400.40.8 P x ( C / m ) PBEsol (I)PBEsol (II)
Cell number along z direction -0.8-0.400.40.8 P y ( C / m ) FIG. 1. Panel (a): Sketch of the supercell used in our simulations. The indicated Cartesian axes coincide with the principaldirections of the perovskite lattice. Panels (b) and (c): Polarization profile corresponding to the stable structure of our multi-domain configuration. (Calculation of local polarizations described in Ref. [29].) The “PBEsol (I)” lines show the results ofan unconstrained first-principles structural relaxation [29]; the agreement with the model predictions is essentially perfect forthe P x profile; for the DW polarization we obtain a slightly smaller value. The “PBEsol (II)” lines show the results of afirst-principles relaxation in which the supercell lattice parameters were fixed to match those predicted by our model potential;the agreement for the P y profile is essentially perfect; for the polarization within the domains we get slightly larger values.Panel (d): Views of the atomic structure of our multi-domain configuration. energy), but not more significant than the also stable,quasi-degenerate parallel configuration.Our results show how strain and the reduced dimen-sionality determine the energetics of the DW-confined FEdistortion. The onset of the multi-domain structure for P x implies a strong deformation of the perovskite lat-tice: it becomes tetragonal and acquires an aspect ratioof 1.07, the long lattice vector coinciding with the po-lar axis x . In the case of our simulated system, the xy plane is homogeneously strained throughout the supercell(stretched along x , shrunk along y ). Hence, even if wehave P x = 0 at the DWs of Fig. 6, the strain disfavors theoccurrence of a DW polarization along the y direction,which is subject to a compression. The ensuing effect canbe appreciated in the potential wells of Fig. 2(b): Case Icorresponds to the full development of the FE distor-tion of PTO, as it happens within the domains. Case IIcorresponds to the development of a three-dimensional y -polarized state when we constrain the cell to be strainedas in the x -polarized FE state; the equilibrium polariza-tion and associated energy gain get clearly reduced, andwe obtain a value of P y (about 0.75 C/m ) that is not farfrom our result at the DW center (about 0.65 C/m ). Ad-ditionally, Fig. 2(b) shows a case III corresponding to thecondensation of P y at our DWs; the energy well becomesshallower than in case II, indicating a further weakeningof the polar instability caused by the spatial confinement(i.e., by the truncation of interactions favoring the three-dimensional homogeneous polar state) and the competi- tion with the P x distortion of the neighboring domains.Nevertheless, the obtained well depth (86 meV/cell) issizable, which suggests that the predicted DW instabil-ity should occur at relatively high temperatures.We checked the correctness of our model-potential pre-dictions by running direct first-principles calculations ofour multi-domain structure, using a 1 × ×
20 cell. Asshown in Figs. 6(b) and 6(c), the agreement betweenour model-potential and first-principles results is verygood, and the FE character of PTO’s DWs is confirmed[30]. We should note that there are several first-principlesstudies of the 180 ◦ DWs of PTO in the literature [18, 31–33], and the consensus is that no DW-confined polariza-tion occurs. We cannot be sure about the reasons whythese previous works did not find polarized DWs; somepossibilities are discussed in [29].We checked whether this confined polarization occursin other PTO DWs. We found that 180 ◦ DWs lyingin other planes – e.g., a (011) boundary separating do-mains with P I k [100] and P II k [¯100] – present polardistortions analogous to the one just described. In con-trast, we found that 90 ◦ DWs do not present any FEinstability, a result probably related with the fact thatthese boundaries are considerably more distorted thantheir 180 ◦ counterparts, or to the elastic (epitaxial ten-sile) constraints we had to impose in order to stabilizethem. PTO’s 90 ◦ DWs will be discussed elsewhere.Finally, let us note that the low-temperature configu-ration of our PTO DWs can be described as being Bloch- -4 -2 0 2 4
Electric field along y (10 V/cm) -0.4-0.200.20.4 P y ( C / m ) TotalBulk-like ∆ P y Difference (DW) (a)(b)
DW polarization Fully developedConstrained
Polarization (C/m ) -200-1000100 E n e r gy p e r ce ll ( m e V ) III II I
FIG. 2. Panel (a): Polarization loop computed for the x -polarized multi-domain state of Fig. 6 and for an electricfield along the y direction. For these simulations we useda 1 × ×
20 supercell, thus assuming an homogeneous switch,and performed the calculations in the limit of 0 K (details in[29]). The total response (black line and circles) is split intwo parts: (1) The response of an x -polarized mono-domainstate (dashed blue line); note that this response does not sat-urate, as the P x polarization will eventually rotate to alignwith the field. (2) The difference between the total and themono-domain results (solid red line), which captures the DWresponse. In order to observe the switch at the DWs, weimposed in the simulations the strain of the multi-domainground state, and thus prevented the domain polarizationsfrom fully aligning with the applied field for relatively smallfield values. Such a constraint is similar to the epitaxial onecharacteristic of thin films. Panel (b): Energy wells corre-sponding to FE instabilities in various situations described inthe text. Note that the P = 0 state corresponds to a differentatomistic configuration in each case. In the case of the DWpolarization (red solid line) the energy is given per cell withinthe DW plane. like. First-principles theory has predicted the occurenceof Bloch-like DW configurations in materials like LiNbO [18] and BaTiO in its rhombohedral phase [21]. Behavior with temperature .– We studied PTO DWsas a function of increasing temperature by running MC T ( K ) f r equen cy FIG. 3. Probability distribution ρ ( P DW y ) for the y -componentof the polarization at the center of the DW and as a func-tion of temperature. The three temperature regions men-tioned in the text are marked with different colors: whitefor T ≥
350 K (where we clearly have h P DW y i = 0), red for320 K < T <
350 K (critical region), and blue with stripesfor T ≤
320 K (where we clearly have h P DW y i 6 = 0). Theleft side of the histograms for T = 330 K and T = 325 K isdrawn using dashed black lines; this is to emphasize that, atthese temperatures, we attribute the P DW y fluctuations frompositive to negative values to finite size effects (see [29]). simulations as described in [29]. Figure 3 shows the ob-tained probability density, ρ ( P DW y ), for the y -componentof the polarization at the DW plane. We observe threedistinct regions: (1) For T ≤
320 K the DW presents astable and large polarization, and the equilibrium stateresembles the one discussed above. (2) For T ≥
350 Kthe DW gets disordered and we have a null thermal av-erage h P DW y i = 0. At those temperatures the systempresents a mirror symmetry plane perpendicular to the y axis, and we could say that the DWs are in a paraelectric state [34]. (3) Finally, for 320 K < T <
350 K we havea narrow region in which, during the course of the MCsimulation, P DW y occasionally switches between equiva-lent polar states. Clearly, the finite size of our simulationsupercell is partly responsible for such fluctuations, whichshould be considered a spurious finite size effect below acertain transition temperature T DWC . At any rate, thepresence of a phase transition is obvious from these re-sults, and the analysis of the MC data suggests that wehave T DWC ≈
335 K. (See [29] for details.)The phase transitions in our simulated system areclearly appreciated in Fig. 4, which shows the evolutionof the relevant order parameters and dielectric response.As the system cools down from high T , the diagonal com-ponents of the dielectric tensor increase sharply, reveal-ing a FE transition at T C = 510 K. [The quantitativedisagreement between our computed T C and the exper- P o l a r i za ti on ( C / m )
200 300 400 500 600 700
Temperature (K) ε / ε P x I P y DW ε xx = ε yy = ε zz ε xx < ε yy ~ ε zz ε yy >> ε yy DW dom (a)(b) T C T CDW
FIG. 4. Panel (a): Temperature dependent polarization atthe center of the first domain ( P I x , blue circles) and at thecenter of the DWs ( P DW y , red triangles). Panel (b): Diagonalcomponents of the dielectric permittivity tensor. The discon-tinuity for ǫ yy (triangles) at T = 390 K is related with thefinite-size effects described in [29]. imental result for PTO (760 K) is discussed in Ref. 27.]Note that for T >
510 K we have ǫ xx = ǫ yy = ǫ zz , reflect-ing the cubic symmetry of the paraelectric phase. Then,for T & T DWC the ǫ yy component [triangles in Fig. 4(b)]behaves in an anomalous way. As regards the dielectricresponse along the y direction, our multi-domain con-figuration can be viewed as a set of parallel capacitors,so that the total dielectric constant can be split into do-main and DW contributions, ǫ yy = f dom ǫ dom yy + f DW ǫ DW yy ,where f dom and f DW are the respective volume fractions.Noting that f dom ≈ . f DW ≈ . ǫ dom yy is featureless for T < T C (which we checkedby running MC simulations of the mono-domain case), itis clear that the increase in ǫ yy comes from a very largeDW response ǫ DW yy . Hence, the DW-confined transitionis driven by a FE soft mode and characterized by verystrong fluctuations of the order parameter, as consistentwith the reduced dimensionality. Final remarks .– As mentioned, the possibility of hav-ing polar orders and dielectric anomalies associated withDWs has been discussed in the theoretical literature,typically employing phenomenological continuum mod-els [16, 22]. Our first-principles work corroborates thatsuch effects can indeed occur, providing for the first timea realistic atomistic picture of a T -driven DW-confinedferroic transformation.One might describe the found transition as a changein the DW character, from Bloch to Ising, upon heating; in fact, some authors have discussed similar effects inthese terms [20, 22]. However, we think that our result isbetter described as a proper FE phase transition confinedto the DW, to emphasize that it results in a switchableDW polarization. For the same reasons, we would ratherdenote the low-temperature state as a ferroelectric DW,and not simply as a DW with Bloch-like character. Notethat, as in the case of LiNbO [18], Bloch-like DWs maynot display a net polarization.The present results are the first step in the investi-gation of such an interesting phenomenon. Aspects forfuture work include: the possible critical behavior ofthe transition, the pre-transitional dynamics, the internalstructure of the DWs (can we have multi-domain stateswithin them?), and the role of DW–DW interactions (canit affect the dimensionality and features of the transi-tion?). We thus believe our findings will open excitingresearch avenues in the fields of phase transitions andlow-dimensional physics.It may seem surprising that the predicted effect hasnot been reported experimentally. However, note thatobserving such a FE order may require some uncom-mon measurements. Ideally, one would like to work withsamples presenting a pattern of highly-ordered stripe do-mains separated by 180 ◦ DWs, as displayed by suit-ably grown PTO films [35] and PTO/STO superlattices[36, 37]. High-resolution X-ray measurements at low tem-peratures may reveal the DW order. Additionally, bymeasuring the dielectric response along the in-plane di-rection of the DWs, one should be able to observe a clearfeature around T DWC . We hope our results will motivatefurther experimental work to characterize FE DWs andthe transitions that may occur within them.Work supported by MINECO-Spain (Grants No.MAT2010-18113 and No. CSD2007-00041) and CSIC[JAE-doc program (JCW)]. Some figures were preparedusing VESTA [38]. We thank Ekhard Salje, David Van-derbilt, and Pavlo Zubko for their useful comments.
SUPPLEMENTAL INFORMATION FOR“FERROELECTRIC TRANSITIONS ATFERROELECTRIC DOMAIN WALLS”
In the following we describe the simulation methodsemployed as well as various technicalities pertaining tothe calculation of some quantities. We also comment onthe determination and quantitative accuracy of the tem-perature at which the domain wall-confined ferroelectrictransition takes place, and on previous first-principles lit-erature on the DWs in PbTiO . METHODS
For most simulations we used the first-principles modelfor PTO described in Ref. 27, where it is labeled “ L I ”.This model can be viewed as a Taylor series of the en-ergy, around the ideal cubic perovskite structure, as afunction of all possible atomic distortions and strains.The series was truncated at 4th order and only pairwiseinteraction terms were included. The potential param-eters were computed by using the local density approx-imation (LDA) to density functional theory (DFT). Tocompensate for LDA’s well-known overbinding problem,we simulate the model under the action of a tensile hy-drostatic pressure of 14.9 GPa.We ran some direct first-principles simulations to verifythe model predictions at low temperatures. For that task,we used the PBEsol[39] approach to DFT as implementedin the VASP package.[40] Note that the PBEsol func-tional has been shown to render accurate results for theequilibrium structures of several FE perovskite oxides[41]including PTO (which we checked ourselves). We usedthe projector augmented wave method to represent theionic cores,[42, 43] solving for the following electrons: Ti’s3 p , 3 d , and 4 s ; Pb’s 5 d , 6 s , and 6 p ; and O’s 2 s and2 p . Wave functions were represented in a plane-wavebasis truncated at 500 eV, and a 6 × × k -point gridwas used for integrations within the Brillouin zone corre-sponding to a 1 × ×
20 supercell. To better compare ourmodel-potential and first-principles results, we conductedPBEsol structural relaxations in both an unconstrainedway [“PBEsol (I)” results in Fig. 1 of our manuscript]and by imposing lattice vectors obtained from model-potential calculations [“PBEsol (II)”].For the sake of comparison, we also ran some simu-lations at the LDA level, all other technicalities beingidentical to those of our PBEsol calculations.We computed the electric dipoles caused by the dis-placement of individual atoms within a linear approxima-tion, using the corresponding dynamical charge tensors ascomputed for the cubic phase. To compute Ti-centeredlocal polarizations associated to specific cells, as thoseshown in Fig. 1, we proceeded in the following way: Wesummed up the dipoles associated to the central titaniumand to its nearest-neighboring atoms (i.e., 6 oxygens and8 leads), weighting the contribution of the neighbors ap-propriately (e.g., since an oxygen atom is shared by twoTi atoms, it contributes one half of its associated dipoleto the polarization centered at each of its two neighbor-ing titaniums). We proceeded analogously to computePb-centered polarizations. This is a rather standardscheme[32] related with the local-mode construction inclassic first-principles models of ferroelectrics.[44]Finally, we used a standard Monte Carlo (MC) schemeto account for the effect of temperature in our system,working with a periodically repeated box of 12 × ×
20 PTO unit cells. At all temperatures we initialized ourMC simulations from the multi-domain ground stateconfiguration described in Fig. 1. We typically ran25,000 MC sweeps for thermalization, followed by atleast 25,000 additional sweeps for computing thermal av-erages. In the region of the DW-confined phase tran-sition, marked in red in Fig. 3 of our manuscript, wedid production runs of up to 150,000 MC sweeps, asneeded to obtain acceptable statistics for ρ ( P DW y ). Notethat, because of the quasi-two-dimensional and strongly-fluctuating nature of the polar order at the walls, thiswas an especially challenging task; in particular, as dis-cussed below, our ability to get quantitatively accurateresults in the critical region is limited by the size of theconsidered supercell, which is the largest we can affordwith our current simulation tools. We also ran a sim-ulation of the FE phase transition of PTO, always in amono-domain configuration, using a 12 × ×
12 simula-tion box periodically repeated. In this case, we alwaysstarted the simulations from a perfectly cubic structure;then we ran 10,000 MC sweeps for thermalization, fol-lowed by 50,000 sweeps for computing averages, whichwas enough to obtain well converged results at all tem-peratures. We checked that the temperature dependentmono-domain polarization thus computed is essentiallyidentical to the local polarization obtained at the middleof the domains in our multi-domain simulations.Whenever we needed to run a structural relaxationwith our model potential, we performed MC simulatedannealings in which we reduced the temperature downto essentially 0 K, allowing the simulation to evolve untilthe system would get frozen at an energy minimum. Insome cases we started the annealings at very low tem-peratures, which allowed us to investigate local (meta-stable) minima of the energy, as the lack of thermal ac-tivation prevents the simulated system from evolving toits most stable configuration. This was important, forexample, to compute the hysteresis loops of Fig. 2(a) ofour manuscript, as we needed to follow the evolution ofthe polar states, as a function of the magnitude of anopposing electric field, up to their meta-stability limit.Beyond that limit, the local energy minimum is lost andthe relaxation (i.e., the simulated annealing) naturallyevolves towards the polar state favored by the appliedfield. Using the same strategy, we were also able to quan-tify the energy gap between the states with parallel andanti-parallel DW polarizations, and thus determine thatthe DW–DW interaction is negligible for the separatingdistance (i.e., 9 unit cells or about 35 ˚A) correspondingto our simulation box.
200 250 300 350 400
Temperature (K) P o l a r i za ti on ( C / m ) FIG. 5. Computed temperature dependence of h| P DW y |i , theabsolute value of the y component of the DW polarization.From the inflection point of this curve, obtained numerically,we can estimate T DWC ≈
332 K.
FINITE SIZE EFFECTS IN THE MONTE CARLOSIMULATIONS
As already mentioned in the manuscript, at tempera-tures around the DW-confined transition our MC resultsare affected by the fact that our simulation supercell isfinite in size. More specifically, we simulate DWs thatare composed of 12 ×
12 cells periodically repeated in theDW plane. While we have checked that a lateral size of12 cells is sufficient to describe the bulk FE transition ofPTO with good accuracy, our results suggest that we maynot be converged size-wise as regards the DW-confinedtransformation (see Fig. 3 of our manuscript). When-ever one deals with low-dimensional orders, large ampli-tude fluctuations are more likely to occur, which com-plicates convergence with size. Unfortunately, at presentour software does not allow us to further increase thesize of the simulation box, as the calculations becometoo heavy computationally. Nevertheless, an analysis ofthe results obtained for 12 ×
12 DWs allows us to bracket T DWC quite precisely.First of all, let us stress that the uncertainty causedby these size effects is of little practical importance. Theaffected temperature range is relatively small, approxi-mately from 320 K to 350 K. Hence, from the point ofview of comparing our value of T DWC with an eventualexperimental determination, the error coming from sizeeffects can be considered negligible.There are simple and usually effective strategies to es-timate the temperature dependence of the order param-eter, from our MC data, in ways that are robust againstsize effects. Thus, for example, if we are interested in thethermal average of the DW polarization h P DW y i , it is agood idea to compute the average of its absolute value, h| P DW y |i , instead. In this way, spurious fluctuations (i.e., jumps between equivalent energy wells) of the DW po- larization, permitted by the finite size of the simulationbox, will not result in a vanishing order parameter. Anobvious drawback of this approach is that h| P DW y |i willbe different from zero even in the disordered phase, whichis not correct. Yet, monitoring the temperature depen-dence of this quantity is usually enough to locate a phasetransition, and even permits a quantitative estimation ofthe transition temperature.Figure S1 shows the results for h| P DW y |i from our MCsimulations. As expected, we obtain a smooth temper-ature dependence that clearly points at a transition oc-curring in the range between 320 K and 350 K. Further,we can estimate the transition temperature by numeri-cally calculating the inflection point of this curve, whichhappens at about 332 K.We were not fully satisfied by this simple estimate, andtried to analyze our MC data in more detail, proceedingas follows. Figure S2(a) shows the computed ρ ( P DW y )evaluated at P DW y = 0 and as a function of temperature.There we can appreciate that at T ≤
320 K the DWs donot visit the P DW y = 0 state, as we have ρ (0) ≈
0. Hence,it is clear that our simulations predict T DWC >
320 K.Above that temperature, there is a range in which weget very small but non-zero values of ρ (0), as well as abimodal shape for ρ ( P DW y ). Do we really have disorderedDWs at these temperatures?To gain insight into the behavior of the system inthis temperature range, we considered several ways tocompute the dielectric constant ǫ yy . Let us recall that ǫ yy = ǫ ( χ yy + 1), and that we have ǫ χ yy = ∂P y ∂E y = βV ( h P y i − h P y i ) , (1)where E y is the applied electric field, V is the supercellvolume, β = 1 /k B T , k B is the Boltzmann constant, and h ... i denotes a thermal average as computed from our MCsimulations. Now, instead of using all the sweeps of ourproduction MC runs to compute χ yy , it was instructiveto consider smaller sets of data that we label by i , com-pute the corresponding χ ( i ) yy for each set, and performa statistical analysis of the results. Figure S2(b) sum-marizes the outcome of such an exercise. We show theresults obtained by considering sets of data composed of16,000 MC sweeps each. (We checked that 16,000 sweepsare typically sufficient to get well-converged thermal av-erages of our quantities of interests.) If we approximate χ yy by ¯ χ yy = 1 M M X i =1 χ ( i ) yy , (2)where M is the number of 16,000-sweep sets that we have,we obtain the results shown by a solid black line andtriangles in Fig. S2(b). Interestingly, ¯ χ yy presents two
275 300 325 350 375 400
Temperature (K) ρ ( P y D W = )
200 250 300 350 400
Temperature (K) ε yy / ε averageMinimum -0.6-0.4-0.200.20.40.6 P y ( C / m ) Monte Carlo sweep -0.6-0.4-0.200.20.40.6 P y ( C / m )
335 K320 K (a) (b)
From data sets of16,000 MC sweeps (c)(d)
335 K
335 K320 K
FIG. 6. Panel (a): Probability ρ (0) of having P DW y = 0 as a function of temperature. Panel (b): ǫ yy component of the dielectricpermittivity tensor, as a function of temperature, and estimated in the two different ways described in the text. Panels (c) and(d): Evolution of P y during 50,000 sweeps of the MC runs at 335 K and 320 K, respectively. The shown data is representativeof the behavior of our simulated system during the longer runs. peaks located at 320 K and 335 K, respectively. If, alter-natively, we approximate χ yy by χ min yy = min i χ ( i ) yy , (3)we get the results shown by a dashed black line and circlesin Fig. S2(b). In this case, we get a single broad peakwith a maximum at T ≈
340 K.To understand the origin of these results, it is usefulto look at the evolution of P y during the MC runs, asshown in Figs. S2(c) and S2(d) for 335 K and 320 K,respectively. At 335 K we observe frequent fluctuationsin the value of P y , suggesting that the DWs are not or-dered. (Note that we have two DWs in the simulationcell; at 335 K we get | P y | ≈ . when the twoDW polarizations are parallel, and | P y | ≈ χ yy presents at that temperature. At the sametime, it is possible to find relatively long segments of theMC run in which P y fluctuates around a constant value.Hence, the computed χ min yy is much smaller than ¯ χ yy .Figure S2(d) shows that the situation at 320 K is verydifferent. In this case, our DWs clearly present a sta-ble polarization. Nevertheless, from time to time we doobserve a jump in P y , indicating a switch in the polar-ization of one of the DWs. These infrequent jumps leadto a very large value of ¯ χ yy , which peaks at that temper-ature. In contrast, we find that χ min yy , which is obviouslydetermined by a segment of the MC trajectory that isfree of jumps in P y , takes a relatively small value. Fromthis analysis, we can conclude that the peak observed at320 K in ¯ χ yy is an artifact associated to the finite sizeof our simulation supercell, and that our results suggestthat T DWC is about 335 K or 340 K. To generate the fig-ures in our manuscript, we assumed T DWC = 335 K andprocessed our MC data accordingly.Finally, let us note another feature of our results that is related with the finite size of our simulation box. InFig. 4(b) of our manuscript, the results for ǫ yy with aclear DW contribution extend only up to T ≈
390 K,where a discontinuity can be appreciated. The reasonis that above that temperature, because of the finitesize our supercell, the initial multi-domain configurationevolves during the MC run into a mono-domain state,which is more stable for the elastic (unconstrained) andelectrical (short-circuit) boundary conditions consideredin our simulations. Nevertheless, we do not observe anysort of structural order at the DWs for T & T DWC , andhave no reason to expect a different behavior at highertemperatures. Hence, our current inability to investigatethe DWs up to T C is not troublesome. FIRST-PRINCIPLES AND MODEL-POTENTIALACCURACY
The largest quantitative uncertainty affecting our valuefor T DWC is surely related with the inaccuracies inher-ent to our first-principles model potentials. As shown inRef. 27, our potential for PTO reproduces the energet-ics of the FE instabilities with very good accuracy, andit also accounts well for the lowest-energy perturbationsof the relevant stable structures. Hence, as far as thehomogeneous FE phases of the material are concerned,our potential can be considered an accurate approxima-tion to density functional theory (DFT) results obtainedusing the local density approximation (LDA). However,the computed T C = 510 K turns out to be significantlysmaller than the experimental value of 760 K, which sug-gests there is a problem with the accuracy of the LDAand employed first-principles methods. As discussed inRef. 27, tracking down the origin of this error is not triv-ial; the LDA could be giving us (1) FE instabilities thatare too weak, (2) AFD instabilities that are too strong,(3) a FE-AFD competition that is too strong, or (4) acombination of all these problems. All such errors willaffect our calculation of T DWC in essentially the same wayas they affect the calculation of T C . Hence, we have rea-sons to believe that our model may underestimate signif-icantly the temperature of the DW-confined transition.Additionally, one may wonder whether our model re-produces the DFT energetics associated with the DWs aswell as it captures those related with the bulk FE state.Our model renders a DW energy E DW = 107 mJ/m when we allow for a full structural relaxation of theDW, and E DW = 190 mJ/m when we constrain theDW to remain in the unpolarized high-symmetry state.Using the PBEsol approach to DFT [which we denote“PBEsol (I)” in our manuscript; see Methods section],we get 149 mJ/m and 152 mJ/m , respectively, for thesetwo quantities. When we run the PBEsol relaxation im-posing the cell strain obtained from our model calcula-tions [“PBEsol (II)”], we get 225 mJ/m and 250 mJ/m ,respectively. Note that Meyer and Vanderbilt reported avalue of 132 mJ/m for this quantity,[32] which they com-puted at the LDA level and considering an unpolarizedDW. Thus, we can conclude that our model gives accept-able results for the DW energy, even though this pieceof information was not used when calculating the poten-tial parameters. (For reference, we ran LDA calculationsconsidering an unpolarized DW, as done in Ref. 32, andgot a DW energy of 97 mJ/m .)For us it is more relevant to consider the energy dif-ference between the paraelectric and ferroelectric statesof the DW, i.e., the depth of the energy well corre-sponding to case III in Fig. 2(b) of our manuscript.For this quantity we get 86 meV/cell with our model,4 meV/cell at the “PBEsol (I)” level, and 26 meV/cellat the “PBEsol (II)” level. It is interesting to note that,despite the good agreement for the DW polarization andstructure between our model and the PBEsol calcula-tions [see Fig. 1(c)], the differences in energy are large.This comparison suggests that, with respect to the re-sult one would obtain from a PBEsol simulation, ourmodel may overestimate significantly the value of T DWC .[For reference, we also computed the depth of the en-ergy well corresponding to the DW-confined FE insta-bility at the LDA level, and obtained 6 meV/cell whena full structural relaxation is allowed, and 45 meV/cellwhen the cell obtained from our model-potential simula-tions is imposed. Hence, LDA renders slightly strongerDW-confined instabilities than PBEsol.]How worrisome are the differences between our modeland PBEsol energetics as regards the DW order? Howwell would a PBEsol simulation reproduce the T C of theFE transition occurring in the bulk of PTO? Is PBEsolthe best DFT flavor available for this task? It is not pos-sible to answer those questions at present. As an indica-tion, let us mention that our PBEsol calculations render avalue of 67 meV/cell for the energy difference between thebulk-paraelectric and bulk-ferroelectric phases of PTO. This value corresponds to the depth of the energy wellof case I in Fig. 2(b), for which we get 189 meV/cell us-ing our model. Hence, since our model underestimatesPTO’s experimental T C quite significantly, we can tenta-tively conclude that a PBEsol simulation would lead tovery small Curie temperature for this material. Hence,despite its accuracy at predicting equilibrium structures,it is not clear at all whether PBEsol should be our DFTmethod of choice for calculating the energetics of PTO’sFE instabilities.We hope this discussion illustrates the difficulties in-volved in assessing the accuracy of DFT, and DFT-basedmethods, as regards the calculation of transition temper-atures in materials with non-trivial structural and lattice-dynamical effects like PTO. As mentioned, we have rea-sons to think that our result for T DWC might be smallerthan the experimental value, and we also have reasons tothink it might be larger. At any rate, the existence of aDW-confined polar order and temperature-driven transi-tion remains a solid prediction.
RELATION TO PREVIOUS LITERATURE
As mentioned in our paper, there are several first-principles studies of the 180 ◦ DWs of PTO available inthe literature [18, 31–33]. Interestingly, except for theearliest one,[31] all of them suggest that no polar orderoccurs within the DW plane. We cannot be sure aboutthe reasons why these previous works did not find polar-ized DWs, but some possibilities are worth mentioning.First, it is usual in studies of this sort to assume a sim-plification of the DW structure, typically imposing thatthe symmetry operations that are common to the twoneighboring domains be preserved in the structural re-laxations. Such an approximation, which is adopted toalleviate the computational cost of the simulations, maywell render a null DW polarization. Second, we obtainedour polar DW configurations from MC annealings basedon our PTO potential, which allowed us to explore theconfiguration space very efficiently, and used such struc-tures as the starting point of the first-principles relax-ations. In our experience, one is not guaranteed to reachsuch a distorted structure in a typical first-principles op-timization, which would start from the non-polar DWstate, even if all symmetries are broken. (Note that theinformation on the energetics of the DW instability, givenin Section III above, may be relevant to this issue.)Let us also note that, recently, Wei et al . [45] havefound polar distortions at the structural (antiphase) DWsof non-polar compound PbZrO . Wei et al . present first-principles results showing the bi-stability of the polarDW configuration, and thus argue that such a polariza-tion is switchable. However, the supplemental materialof Ref. 45 suggests that the switch of the DW polariza-tion also involves the reversal of structural distortions inside the domains. This is reminiscent of the so-called hybrid improper ferroelectrics – where the polarizationswitch must be accompanied by the reversal of one addi-tional order parameter [46–49] – and of the polar orderpredicted at the FS walls of CaTiO [17]. In contrast,the polarization of our PTO DWs can switch withoutaffecting the surrounding domains. [1] J. Seidel, L. W. Martin, Q. He, Q. Zhan, Y. H. Chu,A. Rother, M. E. Hawkridge, P. Maksymovych, P. Yu,M. Gajek, N. Balke, S. V. Kalinin, S. Gemming, F. Wang,G. Catalan, J. F. Scott, N. A. Spaldin, J. Orenstein, andR. Ramesh, Nature Materials , 229 (2009).[2] J. Guyonnet, I. Gaponenko, S. Gariglio, and P. Paruch,Advanced Materials , 5377 (2011).[3] S. Farokhipoor and B. Noheda, Physical Review Letters , 127601 (2011).[4] A. Aird and E. K. H. Salje, Journal of Physics: Con-densed Matter , L377 (1998).[5] S. Y. Yang, J. Seidel, S. J. Byrnes, P. Shafer, C. H. Yang,M. D. Rossell, P. Yu, Y. H. Chu, J. F. Scott, J. W. Ager,L. W. Martin, and R. Ramesh, Nature Nanotechnology , 143 (2010).[6] M. Alexe and D. Hesse,Nature Communications (2011), 10.1038/ncomms1261.[7] Q. He, C. H. Yeh, J. C. Yang, G. Singh-Bhalla, C. W.Liang, P. W. Chiu, G. Catalan, L. W. Martin, Y. H. Chu,J. F. Scott, and R. Ramesh, Physical Review Letters , 067203 (2012).[8] M. Daraktchiev, G. Catalan, and J. F. Scott,Physical Review B , 224118 (2010).[9] Z. V. Gareeva and A. K. Zvezdin, Physics of the SolidState , 1714 (2010).[10] E. K. Salje, ChemPhysChem , 940 (2010).[11] G. Catalan, J. Seidel, R. Ramesh, and J. F. Scott, Re-views of Modern Physics , 119 (2012).[12] B. Strukov and A. Levanyuk, Ferroelectric Phenomenain Crystals: Physical Foundations (Springer, 1998).[13] A. Tagantsev, L. E. Cross, and J. Fousek,
Domains inFerroic Crystals and Thin Films (Springer, 2010).[14] J. J. Lajzerowicz and J. J. Niez, Journal de PhysiqueLetters , L165 (1979).[15] B. Houchmandzadeh, J. Lajzerowicz, and E. Salje, Jour-nal of Physics: Condensed Matter , 5163 (1991).[16] A. K. Tagantsev, E. Courtens, and L. Arzel, PhysicalReview B , 224107 (2001).[17] L. Goncalves-Ferreira, S. A. T. Redfern, E. Artacho, andE. K. H. Salje, Physical Review Letters , 097602(2008).[18] D. Lee, R. K. Behera, P. Wu, H. Xu, Y. L. Li, S. B.Sinnott, S. R. Phillpot, L. Q. Chen, and V. Gopalan,Phys. Rev. B , 060102 (2009).[19] S. Conti, S. M¨uller, A. Poliakovsky, and E. K. H. Salje,Journal of Physics: Condensed Matter , 142203 (2011).[20] V. Stepkova, P. Marton, and J. Hlinka, Journal ofPhysics: Condensed Matter , 212201 (2012).[21] M. Taherinejad, D. Vanderbilt, P. Marton, V. Stepkova,and J. Hlinka, Physical Review B , 155138 (2012).[22] P. Marton, V. Stepkova, and J. Hlinka, Phase Transi- tions , 103 (2013).[23] W. Zhong and D. Vanderbilt, Physical Review Letters , 2587 (1995).[24] J. F. Scott, E. K. H. Salje, and M. A. Carpenter, PhysicalReview Letters , 187601 (2012).[25] E. A. Eliseev, A. N. Morozovska, Y. Gu, A. Y. Borisevich,L.-Q. Chen, V. Gopalan, and S. V. Kalinin, PhysicalReview B , 085416 (2012).[26] E. K. H. Salje, O. Aktas, M. A. Car-penter, V. V. Laguta, and J. F. Scott,Physical Review Letters , 247603 (2013).[27] J. C. Wojde l, P. Hermet, M. P. Ljungberg, P. Ghosez,and J. ´I˜niguez, Journal of Physics: Condensed Matter , 305401 (2013).[28] I. A. Kornev, L. Bellaiche, P.-E. Janolin, B. Dkhil, andE. Suard, Physical Review Letters , 157601 (2006).[29] Please see the online Supplemental Materials for a briefdescription of various technicalities of our calculationsand some complementary results and analysis.[30] For the FE phase of PTO, our model potential predictslong and short lattice constants of 4.21 ˚A and 3.94 ˚A,respectively. In contrast, our PBEsol calculations ren-der 4.03 ˚A and 3.87 ˚A for these quantities. These resultsare consistent with those in Figs. 6(b) and 6(c), as thesmaller lattice parameters lead to relatively small polardistortions, both in the domains and within the DWs, inthe PBEsol (I) case.[31] S. P¨oykk¨o and D. J. Chadi, Applied Physics Letters ,2830 (1999).[32] B. Meyer and D. Vanderbilt, Physical Review B ,104111 (2002).[33] L. He and D. Vanderbilt,Physical Review B , 134103 (2003).[34] In this high- T regime, the DW presents all the symme-tries that are common to its neighboring domains.[35] S. K. Streiffer, J. A. Eastman, D. D. Fong,C. Thompson, A. Munkholm, M. V. Ramana Murty,O. Auciello, G. R. Bai, and G. B. Stephenson,Physical Review Letters , 067601 (2002).[36] P. Zubko, N. Stucki, C. Lichtensteiger, and J. M.Triscone, Physical Review Letters , 187601 (2010).[37] P. Zubko, N. Jecklin, A. Torres-Pardo, P. Aguado-Puente, A. Gloter, C. Lichtensteiger, J. Junquera,O. St´ephan, and J. M. Triscone, Nano Letters , 2846(2012).[38] K. Momma and F. Izumi, Journal of Applied Crystallog-raphy , 653 (2008).[39] J. Perdew, A. Ruzsinszky, G. Csonka, O. Vydrov,G. Scuseria, L. Constantin, X. Zhou, and K. Burke,Physical Review Letters , 136406 (2008).[40] G. Kresse and J. Furthm¨uller, Physical Review B ,11169 (1996).[41] R. Wahl, D. Vogtenhuber, and G. Kresse, Physical Re-view B , 104116 (2008).[42] P. E. Bl¨ochl, Physical Review B , 17953 (1994).[43] G. Kresse and D. Joubert, Physical Review B , 1758(1999).[44] W. Zhong, D. Vanderbilt, and K. M. Rabe, PhysicalReview B , 6301 (1995).[45] X.-K. Wei, A. K. Tagantsev, A. Kvasov,K. Roleder, C.-L. Jia, and N. Setter,Nature Communications (2014), 10.1038/ncomms4031.[46] E. Bousquet, M. Dawber, N. Stucki, C. Lichtensteiger, P. Hermet, S. Gariglio, J.-M. Triscone, and P. Ghosez,Nature , 732 (2008).[47] N. A. Benedek and C. J. Fennie, Physical Review Letters , 107204 (2011). [48] Z. Zanolli, J. C. Wojde l, J. ´I˜niguez, and P. Ghosez, Phys-ical Review B , 060102 (2013).[49] Y. Yang, J. ´I˜niguez, A.-J. Mao, and L. Bellaiche, Phys-ical Review Letters112