Interplay of epitaxial strain and rotations in PbTiO_{3}/PbZrO_{3} superlattices from first principles

Jeroen L. Blok,^{∗} Dave H.A. Blank, and Guus Rijnders

Faculty of Science and Technology and MESA+ Institute for Nanotechnology, University of Twente, 7500 AE Enschede, The Netherlands

Karin M. Rabe and David Vanderbilt

Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854-8019, USA (Dated: August 16, 2011)

We present first-principles calculations of the structural phase behavior of the [1:1] PbTiO_{3}/PbZrO_{3} superlattice and the PbTiO_{3} and PbZrO_{3} parent compounds as a function of in-plane epitaxial strain. A symmetry analysis is used to identify the phases and clarify how they arise from an interplay between different kinds of structural distortions, including out-of-plane and in-plane polar modes, rotation of oxygen octahedra around out-of-plane or in-plane axes, and an anti-polar mode. Symmetry-allowed intermode couplings are identified and used to elucidate the nature of the observed phase transitions. For the minimum-period [1:1] PbTiO_{3}/PbZrO_{3} superlattice, we identify a sequence of three transitions that occur as the in-plane lattice constant is increased. All four of the phases involve substantial oxygen octahedral rotations, and an antipolar distortion is important in the high-tensile-strain phase. Inclusion of these distortions is found to be crucial for an accurate determination of the phase boundaries.

PACS numbers: 77.80.bn,77.55.hj,77.55.Px,77.84.Cg,68.65.Cd

I. INTRODUCTION

Ferroelectric materials have long been the subject of intense study both because of their fundamental scientific interest and because their properties are promising for device applications. Experimental and theoretical studies of the relationships between their symmetries, structures and properties are essential for the development of improved materials. One of the most widely used and thoroughly studied materials is PbZr_{x}Ti_{1−x}O_{3} (PZT), which in bulk shows an enhanced piezoelectric effect near the rhombohedral-tetragonal morphotropic phase boundary (MPB) around x = 0.5. Noheda et al.^{1} explained this anomalous response in connection with their discovery of a narrow wedge of monoclinic phase, bridging between tetragonal and rhombohedral phases at the MPB.

Modern growth techniques such as pulsed laser deposition and molecular-beam epitaxy now allow for the outof-equilibrium synthesis of atomic-scale perovskite superlattices of arbitrary layer sequences, and control of the in-plane lattice constant via coherent epitaxy on chosen perovskite substrates.^{2} As the behavior of a bulk solid solution and a superlattice having the same overall composition can in principle be quite different, these developments offer an exciting opportunity for the development of improved new materials. However, given the enormous number of possible structures that can be constructed in this way, it is clear that theory has an important role to play in guiding the search.

Not surprisingly, then, there has been a simultaneous explosive growth in the application of first-principles theoretical methods to the study of perovskite superlattice structures,^{3} specifically to the prediction of structure, polarization and polarization-related properties such as the dielectric and piezoelectric response. A previous study by Bungaro and Rabe^{4} considered [1:1] superlattices (alternating single unit cells) of PTO and PZO, finding a sequence of structural phase transitions as a function of epitaxial in-plane lattice constant. However, the possibility of oxygen octahedron rotations was not taken into account. Such rotations, which arise from antiferrodistortive instabilities in many perovskite materials, may in many cases compete with the ferroelectric distortions, as shown by the early work of Zhong and Vanderbilt.^{5} In subsequent first-principles work, Ghosez et al.^{6} identified antiferrodistortive instabilities in a number of ferroelectric perovskites, while Fornari and Singh^{7} showed that octahedral rotations in PZT play an important role in the morphotropic phase boundary. Most recently, combinations of rotations at the interfaces were found in firstprinciples investigations to generate improper ferroelectricity in PT/ST superlattices^{8} and in layered perovskite compounds.^{9}

In the present work, we apply first-principles DFT methods to determine the ground-state phases of the [1:1] PTO/PZO superlattice as a function of epitaxial strain, as well as those of pure PTO and PZO for reference. We go beyond previous studies by including octahedral rotations and tilts (both of which shall henceforth be referred to as “rotations”), which can be accommodated in a \sqrt{2}×\sqrt{2}×2 supercell; the types of distortions identified as playing an important role are shown in Fig. 1. We find that octahedral rotations are ubiquitous throughout the phase diagram of the PTO/PZO superlattice, and have an significant influence on the polarization and the locations of the phase boundaries. In particular, they help stabilize a monoclinic phase over a very broad range of epitaxial strain. We relate these results to the behavior

FIG. 1: Sketch of distortions of type (a) P_{001}, (b) P_{110}, (c) R_{001}, (d) R_{110} and (e) A_{¯110} in the PTO/PZO superlattice.

of the PTO and PZO constituents and analyze the results in terms of the symmetries of the various phases. Our results clarify the role played by the coupling between these degrees of freedom in determining the structure and polarization of the superlattice structure under consideration.

This paper is organized as follows. In Sec. II we describe the methods used to apply DFT to bulk PTO and PZO and to the PTO/PZO superlattice. The results of these calculations are presented is Sec. III, and a discussion is given in Sec. IV. We summarize and conclude in Sec. V.

II. METHODS

A. Computational details

The effect of rotations on the PTO/PZO superlattice was studied using density-functional theory (DFT) calculations within the local density-approximation (LDA).^{10,11} The calculations were performed using the Vienna Ab initio Simulation Package (VASP).^{12,13} The Brillouin zone was sampled with a 4x4x4 Monkhorst-Pack grid.^{14} The projector augmented wave method (PAW)^{15,16} was used with an 800 eV plane-wave cutoff. The polarization for all systems was calculated using the Berry-phase method^{17} as implemented in VASP.

For comparison with the PTO/PZO superlattice, we performed analogous studies of PTO and PZO with the same epitaxial constraints. For all three systems, the structure was first initialized by imposing selected polar and/or rotational distortions, not limited to those shown in Fig. 1, and the system was then relaxed subject to the constraint of fixed in-plane lattice parameters, corresponding to the epitaxial boundary condition. While the symmetry of some of the phases allowed for a tilting of the c axis away from the normal to the a-b plane, tests described below in Sec. II D show that this effect is small, and the calculations reported here were all carried out with the c axis kept parallel to ẑ. The in-plane lattice parameters were varied over the range 3.80 Å < a < 4.10 Å for PTO, 3.90 Å < a < 4.20 Å for PZO and 3.90 Å < a < 4.25 Å for the PTO/PZO superlattices.

In order to consider the effects of rotations of oxygen octahedra, the calculations were carried out on a 20-atom \sqrt{2} × \sqrt{2} × 2 unit cell. This supercell allows distortions corresponding to rotational phonon modes at the R point (π/a, π/a, π/a) and at the M point (π/a, π/a, 0), as well to polar zone-center modes associated with ferroelectric distortions (see Fig. 1). In all cases considered, we found that structures with P_{100} > 0 had higher energies than those with P_{110}> 0. Therefore in the following we limit our consideration to structures with nonzero P_{001} and/or P_{110} as shown in Figs. 1(a-b). Furthermore, when comparing PTO and PZO phases incorporating R-point (outof-phase) rotations to those with M-point (in-phase) rotations, we always found the R-point patterns to be lower in energy. Similarly, with one exception that will be noted in due course, we always found the rotations of the ZrO_{6} and TiO_{6} octahedra to have opposite signs in the PTO/PZO superlattices. In what follows, therefore, out-of-phase rotations of the two layers stacked along ẑ are to be understood unless otherwise specified. Rotations about [001] and [110] axes will be denoted by R_{001} and R_{110} respectively, as in Figs. 1(c-d). The 20-atom supercell also accommodates antipolar distortions of the type shown in Fig. 1(e), denoted as A_{¯110}. Structural energies for the PTO, PZO and PTO/PZO systems are given relative to the energy of the corresponding ideal high-symmetry P 4/mmm structure (i.e., with no rotations or polarization) at a reference lattice constant of 3.99 Å.

B. Small distortions

Typically we started our structural minimizations from an initial geometry of relatively low symmetry, corresponding to a mixture of several kinds of distortion, so as to be conservative in our assumptions about what kind of structural distortions would actually occur. For example, in the PTO calculations to be described below, we started from a system with the P_{001}, P_{110}, R_{001} and R_{110} distortions of Figs. 1(a-d) present. After structural relaxation, some of these components (especially the rotational ones) would become very small but not disappear completely. To determine whether these small distortions were physically meaningful, or just a result of incomplete convergence or related numerical difficulties, the distortions were manually turned off (i.e., the structure was symmetrized) and then the total energy was calculated as the distortions were added back in in small increments. If the resulting energy curve showed a minimum away from zero distortion amplitude, the distortion was taken to be physically meaningful. This analysis proved useful, in particular, in clarifying the phase diagram for PTO in Sec. III A.

C. Location of phase boundaries

As we shall show later, the case of the PTO/PZO superlattice presents a related difficulty. Here, we find two second-order transitions at which one or more order parameters go to zero in an apparent square-root singularity. In order to get a more precise location for the critical epitaxial lattice constant at which this occurs, we have carried out the following procedure. We describe this method for the case of the transition from the lower-symmetry P 1a1 space group to the highersymmetry P 2an space group with increasing epitaxial lattice constant in the tensile-strain region of PTO/PZO, but similar considerations apply to the transition from P 1a1 to P 4mm in the compressive-strain region.

First, we chose a lattice constant near the transition, but far enough away that the distortion leading to the low-symmetry phase was well developed. In our example, we chose a=4.125 Å, in the P 1a1 phase but near the transition. We symmetrized this structure using the symmetry operations of the higher-symmetry phase, constructed the eigenvector in displacement space pointing from the symmetrized to the low-symmetry structure, and used this to define a symmetry-breaking mode amplitude. Next, for each lattice constant near the transition (here for a=4.100 to 4.250 Å in increments of 0.025 Å), we started from the relaxed high-symmetry structure and added distortions corresponding to this eigenvector with several amplitudes (and no further relaxation), fitting to an expression of the form

ΔE(u) = \frac{1}{2}κu^{2} + \frac{1}{4}αu^{4}  (1)

to the energy vs. mode amplitude u (relative to u=0) as computed from the first-principles calculations. Not surprisingly, we found that κ is strongly dependent on inplane lattice constant a, while α is not. We then plotted κ vs. a (typically obtaining a nearly linear behavior) and identified the critical lattice constant a_{crit} as the one at which this curve crosses through zero. As we shall see later, we found the transition to occur at a_{crit}=4.171 Å for the transition in this example. As we shall also see, using the same analysis on the compressive train side we found the transition from P 1a1 to P 4bm to occur at a_{crit}=3.890 Å.

D. Treatment of c-axis tilting

As mentioned above, some of the space groups considered in this work, specifically the monoclinic ones, allow by symmetry for the c axis to tilt away from the [001] direction, typically by developing a small [110] component. To study the effect of such a shear strain, especially on the location of second-order phase boundaries as described in the previous subsection, we again made use of the distortion pattern leading from the higher-symmetry tetragonal or orthorhombic state to the lower-symmetry monoclinic one as obtained in Sec. II C. Again denoting the amplitude of this normalized distortion by u, and starting from the higher-symmetry structure at the lattice constant of interest, we carried out calculations of changes in total energy and stress as a function of u and shear strain (axis tilt) η. We then parametrized the results to the first-principles calculations using a generalized version of Eq. (1) taking the form

ΔE(u, η) = \frac{1}{2}κu^{2} + \frac{1}{4}αu^{4} + βuη + \frac{1}{2}γη^{2}.  (2)

Again κ was the only parameter that was found to depend strongly on a. Setting ∂Ẽ(u, η)/∂η = 0 and ∂Ẽ(u, η)/∂u = 0 within the model, we can estimate both the shift of the phase boundary and the amount of c-axis tilt. For the phase boundary when the PTO/PZO superlattice system goes from the P 1a1 space group to the P 2an space group, for example, this resulted in a shift of the critical in-plane lattice constant by only 0.0016 Å, to 4.173 Å, leading to small shear strains near the phase boundary (∼0.2%). Similar results were found for other phase boundaries between monoclinic and tetragonal or orthorhombic phases.

The impact of the tilting of the c axis appears to be less profound in our study than in that of Ref. 4, in which octahedral rotations were not included. Indeed, we find it to be small enough that it will not change any of the results significantly, and so we have not included it in the results to be presented in the next section.

III. RESULTS

A. Parent phases

The results of the DFT calculations on PTO are given in Fig. 2. For in-plane lattice parameters less than about a = 3.90 Å, we find the stable structure to be the tetragonal P 4mm structure (space group 99) in which only P_{001} is present (ferroelectric “c phase,” 5-atom primitive cell), as in Fig. 1(a). For larger lattice parameters, we see a transition at a = 3.91 Å to an orthorhombic Ima2 structure (space group 46, 10-atom primitive cell) with nonzero P_{110} accompanied by an R_{110} rotation of ∼2^{∘}, combining the distortions of Figs. 1(b) and (d). Since there is no group-subgroup relationship between these

FIG. 2: (Color online) Calculated properties of PTO as a function of epitaxially constrained in-plane lattice constant.

(a) Minimum-energy structure: octahedral rotations (black) of type R_{110} (∗); polarizations (red) of type P_{001} (∘) and P_{110} (×). (b) Total energy (black squares) and c/a ratio (red triangles).

two phases, the transition must be first order. These results are in rather good agreement with the previous first-principles calculations of Refs. 4 and 18 (assuming mixed domain phases are not considered), even though those works did not include the octahedral rotations. This suggests that the small rotations that we predict on the higher-a side of the transition do not have a very profound effect on the transition in PTO. We also note that all of the first-principles results disagree with the zerotemperature prediction of a phenomenological analysis^{19} where a window of r phase (P_{001} and P_{110} both present) was predicted at T =0.

As is well known, bulk PZO has a complicated 40-atom unit cell built from a 2\sqrt{2} × \sqrt{2} × 2 enlargement of the parent perovskite unit cell, incorporating antiferroelectric distortions as well as octahedral rotations. While our supercell clearly cannot accommodate this structure, we think it unlikely that this kind of complex pattern would survive in the superlattice geometry. Since our immediate goal here is to provide a reference for the PTO/PZO superlattice calculations, we continue to restrict our calculations to structures realizable within our 20-atom simulation cell. After considering a wide variety of structures consisting of different combinations of polar and rotational distortions, we find that the monoclinic A1a1 phase (space group 9, C1c1 in the conventional set-

FIG. 3: (Color online) Calculated properties of PZO as a function of epitaxially constrained in-plane lattice constant.

(a) Minimum-energy structure: octahedral rotations (black) of type R_{110} (∗) and R_{001} (+); polarizations (red) of type P_{001} (∘) and P_{110} (×). (b) Total energy (black squares) and c/a ratio (red triangles).

ting, 10 atoms/cell), is lowest in energy throughout the range of in-plane lattice constants considered, as shown in Fig. 3. All four of the distortions P_{001}, P_{110}, R_{001}, and R_{110} shown in Figs. 1(a-d) are present in this phase.

B. PTO/PZO superlattice

We now turn to the case of the PTO/PZO superlattice. Compared to the relatively simple behavior for pure PTO and pure PZO, the phase sequence for the PTO/PZO superlattices presented in Fig. 4 is quite rich.

We first remind the reader of the previous results of Bungaro and Rabe,^{4} who confined their considerations to structures that do not include octahedral rotations. Within this constraint, these authors found a sequence of two phase transitions, first from a tetragonal P 4mm structure (“c phase”) for compressive strains a < 3.96 Åto a monoclinic Cm structure (“r phase”) at intermediate strains, and finally to an orthorhombic Amm2 structure (“aa phase”) for tensile strains a > 4.12 Å. Our results, presented in Fig. 4, show a qualitatively similar behavior for the polarization, although all of the phases we observe include octahedral rotations, and the phase boundaries (vertical dashed lines in Fig. 4) have shifted to make the region of monoclinic phase wider. We give a detailed presentation of our results in the following, first in the region of compressive strain, and then in the tensile region, where the behavior is especially subtle.

1. Compressive strain region

Starting at the smaller lattice constants a < 3.890 Å, we find a tetragonal P 4bm structure (space group 100) in which P_{001} and R_{001} are both present (20-atom primitive cell). The octahedral rotation angles for R_{001} in the PTO and PZO planes are quite different, with the PTO rotations remaining rather small, ∼1^{∘}. (Throughout the figure, only the PZO rotations are plotted.) When the lattice constant exceeds 3.890 Å, we find that P_{110} and R_{110} turn on simultaneously in a second-order transition, leading to a monoclinic P 1a1 structure (space group 7, P 1c1 in the conventional setting, 20-atom primitive cell). This structure is very similar to the A1a1 monoclinic structure obtained for PZO in Fig. 3, although the inequivalence of the PTO and PZO layers makes the P 1a1 structure more complicated. The quoted value of a_{crit}=3.890 Å was obtained using the approach described in Sec. II C, but applied in the compressive-strain region, and using the distorted structure at a=3.925 Å to define the distortion eigenvector.

2. Tensile strain region

Next, as the lattice constant increases and approaches a value of about 4.19 Å, we find that the amplitudes of P_{001} and R_{001} diminish and approach zero together. Again using the method of Sec. II C, we find that these mode amplitudes disappear in a second-order transition at a=4.171 Å, leading to an orthorhombic P 2an phase (space group 30, P nc2 in the conventional setting) in which only P_{110} and R_{110} are present (20-atom primitive cell). This structure is similar to the Ima2 phase seen at larger lattice constants for PTO in Fig. 2, except that the PTO and PZO planes are inequivalent.

However, at a marginally larger lattice constant (a=4.193 Å), a first-order transition occurs to an orthorhombic P 2_{1}am phase (space group 26, P mc2_{1} in the conventional setting, 20-atom primitive cell). The value of a_{crit}=4.193 Å was obtained from the crossing of curves of computed energy vs. lattice constant as computed in the P 1a1 and P 2_{1}am symmetries. The P 2_{1}am phase is one in which R_{110} disappears, P_{110} is approximately unchanged, and A_{¯110} and R_{001} distortions appear. A_{¯110}, denoted by the solid red dots in Fig. 4, consists of displacements (largest for Pb, then Ti, then Zr, but also involving O atoms) that are parallel and antiparallel to (¯110), patterned in such a way that no net polarization develops along this axis, as shown in Fig. 1(e). (Such antipolar distortions are also symmetry-allowed in the monoclinic P 1a1 phase, but we find them to be very small, of order 0.02 Å, suggesting that they are induced there

FIG. 4: (Color online) Calculated properties of PTO/PZO superlattice as a function of epitaxially constrained in-plane lattice constant. (a) Minimum-energy structure: octahedral rotations (black) of type R_{110} (∗) and R_{001} (+); polarizations (red) of type P_{001} (∘) and P_{110} (×); antipolar distortions of type A_{¯110} (filled red circles). (b) Total energy (black squares) and c/a ratio (red triangles). Heavy and light vertical dotted lines indicate phase boundaries determined from our work, and from Ref.4, respectively.

by intermode couplings as will be discussed in Sec. IV B.) Contrary to the usual case, the signs of R_{001} are the same in the PTO and PZO planes, with the rotations in both layers being small (∼1^{∘}). The fact that the displacements associated with A_{¯110} are much larger than those associated with R_{001} suggests that A_{¯110} and P_{110} constitute the primary instabilities; when these modes freeze in, they induce an R_{001} distortion, as will become clear from the symmetry analysis to be given in Sec. IV.

The fact that the second-order and first-order transitions occur at almost the same lattice constant has to be regarded as a coincidence. The window of P 2an phase found in our calculation, 4.171 Å< a <4.193 Å, is in fact so narrow that we cannot be confident of its existence at all. Small changes in the technical details of the calculations, such as the exchange-correlation potential or choice of pseudopotentials, could erase or widen this window of stability. The inclusion of c-axis tilting, as discussed in Sec. II D, also tends to reduce or eliminate the window of P 2an phase. We can only conclude that the P 2an phase is competitive with the other two phases near the transition from the P 4bm to the P 2_{1}am phase, and may appear in a narrow window between them. If the P 2an phase is not in fact present, then there would be a direct transition from P 4bm to P 2_{1}am. (Despite the presence of a group-subgroup relation, the transition would be first order because of the discontinuity in R_{110}.)

IV. DISCUSSION

This system illustrates several interesting aspects of complex perovskite oxides. Perovskite oxides generally are characterized by multiple unstable modes in the cubic structure, which couple and compete to produce lowenergy states. The instabilities also each separately couple to strain. Changing strain not only produces or eliminates instabilities, but it can also change the balance of competition and produce first-order phase transitions between distinct phases without a group-subgroup relation. Here, we see that combining two different perovskite constituents, each with characteristic unstable modes and strain dependence, can lead to a great increase in richness of the observed phases.

A. Parent phases

In the present case, the computed phases in the two constituent materials are relatively simple. The dominant instability in cubic PTO is the zone-center polar mode, which couples to strain so that elongation is favored along the polar direction. Compressive epitaxial strain thus favors P_{001}, while tensile epitaxial strain favors P_{110}. The absence of a monoclinic phase with rotating polarization, or equivalently the first-order character of the transition, can be understood as arising from the strong coupling of the polarization to tetragonal strain, as demonstrated by the quartic anisotropy of the polarization energy given in Table III of Ref. 18. Additional symmetry breaking for tensile strain results from the rotational instability R_{110}. While compressive strain is expected to favor the rotational distortion R_{001}, no instability occurs in the range of strain considered here.

In contrast, cubic PZO has multiple instabilities that occur together, yielding a monoclinic phase under the constraint of the supercell used for the first principles calculations. This phase is strongly favored and appears for all values of strain considered, with compressive strain favoring the modes P_{001} and R_{001} and tensile strain favoring the modes P_{110} and R_{110}, but not with such a strong dependence as to eliminate any of the modes at the highest strains considered. Also in contrast to PTO, the polarization freely rotates in the (1¯10) plane with changing epitaxial strain; this was already apparent in the first-principles results for strained PZO neglecting rotations, which included an intermediate monoclinic phase between the phases with normal and in-plane polarization.^{4,18}

B. PTO/PZO superlattice

We now turn to the superlattice phase sequence. The point of closest resemblance to the phases of the constituents is the monoclinic phase for intermediate strain, which shows the rotating polarization and accompanying rotations observed in PZO. However, the strain dependence of the mode amplitudes is much stronger. With increasing compressive strain, P_{110} and R_{110} disappear, yielding a phase with only P_{001} and R_{001}. With increasing tensile strain P_{001} and R_{001} disappear, yielding a phase with only P_{110} and R_{110}, which at slightly larger in-plane lattice constant undergoes a first-order transition to a phase dominated by P_{110}, A_{¯110} and R_{001}.

A symmetry analysis of the Taylor expansion of the energy around a high-symmetry reference structure in the amplitudes of the relevant distortion modes is useful for understanding the mode contents of the various phases observed. For the pure compounds, the highsymmetry reference structure is taken as the tetragonal P4/mmm structure with 5 atoms per cell and all atoms at high-symmetry positions. For the superlattice, the high-symmetry reference structure is the 10-atom P4/mmm structure. We consider the distortion modes shown in Figure 1; for the pure compounds there is an extra translation symmetry relating the two B cations and the modes shown are either even (P_{001}, P_{110}, A_{¯110}) or odd (R_{001}, R_{110}) under this translation.

For pure PTO, the analysis is simple. For compressive strain, the P_{001} mode is unstable and is the only mode included in the lowest-energy phase. For tensile strain, the lowest-energy phase includes both P_{110} and R_{110} distortions. In previous work^{4} we have seen that the P_{110} mode is unstable for tensile strain. R_{110} is thus nonzero either because it is independently unstable, as would be favored by tensile strain, or because the biquadratic coupling between P_{110} and R_{110} is cooperative and large enough to destabilize R_{110}.

For pure PZO, the four modes P_{001}, P_{110}, R_{001}and R_{110}are included in the lowest-energy phase. It is reasonable to assume that the polar distortions P_{001} and P_{110} are unstable throughout this range, with instability in R_{001} promoted by compressive strain and the instability in R_{110} promoted by tensile strain. However, the rotational modes need not be independently unstable. Biquadratic coupling to the polar modes could produce the effective instability. Furthermore, there is a fourth order invariant R_{110}P_{110}R_{001}P_{001} in the energy expansion that ensures that if three of the four distortions are present, the fourth will be induced through an improper mechanism even if it is not unstable.

The effect of this fourth-order coupling is particularly striking in the superlattice, as it is responsible for the coupled appearance/disappearance of P_{110} and R_{110} at a=3.890 Å and the coupled appearance/disappearance of P_{001} and R_{001} at a=4.171 Å. The fact that these modes appear and disappear together might at first seem surprising, since in the high-symmetry reference structure, these two distortions have different symmetry and the introduction of one mode into the reference structure would not induce the other. However, with the lowered symmetry due to the presence of P_{001} and R_{001}, the two modes P_{110} and R_{110} belong to the same irreducible representation and therefore are expected to occur together. This is formally equivalent to the presence of the fourth-order invariant R_{110}P_{110}R_{001}P_{001}, discussed above. For example, in the compressive-strain phase with P_{001} and R_{001}, if an instability with respect to P_{110} or R_{110} develops, R_{110} or P_{110} will be linearly induced through this fourth-order term, producing the intermediate monoclinic phase and ensuring that the two modes appear/disappear together.

It is in the extension of this symmetry analysis to include the A_{¯110} distortion that differences between the pure compounds and the superlattice arise. In the pure compounds, the A_{¯110} distortion is a symmetry-breaking distortion that would appear only if it were unstable, which is not the case for the strain range considered. In the superlattice, there are two third-order couplings R_{001}P_{110}A_{¯110} and R_{110}P_{001}A_{¯110}, not allowed in the pure compounds due to the additional translation symmetry. In the intermediate-strain monoclinic phase, both of these terms would linearly induce an A_{¯110} distortion, which is indeed present at small amplitude in the computed structure for this phase.

In contrast to its minor role in the intermediate-strain monoclinic phase, A_{¯110} appears to be a strong independent instability in the superlattice at higher tensile strain. This can explain the first-order transition at a=4.193 Å. The approximate continuity of P_{110} through the firstorder transition suggests that P_{110} is by far the dominant instability at higher tensile strain. Energy lowering from this polar phase can be produced in two ways: through freezing in an R_{110} instability to obtain the P nc2 phase, or by freezing in an A_{¯110} instability, which induces a small R_{001} distortion through the third-order invariant A_{¯110}R_{001}P_{110} and yields the P mc2_{1} phase. The key to the first-order transition is the competition of these two instabilities in the presence of P_{110}.

C. Growth of PTO/PZO and related superlattices

We have attempted the growth of [1:1] PTO/PZO superlattices using pulsed laser deposition (PLD), but without notable success. Typically, reflection high-energy electron diffraction (RHEED) indicates that an initial layer-by-layer growth rapidly gives way to a 3D growth mode. Our best result to date is for growth on a NdScO_{3} (NSO) substrate, where we observed five successful coherent repetitions of the [1:1] PTO/PZO unit before the 3D growth mode took over. Previous work on SrTiO_{3}/SrZrO_{3}^{20} and KTaO_{3}/KNbO_{3}^{21} systems shows the synthesis of B-site superlattices is possible. We speculate that the development of new substrates (or new substrate treatments) might help extend the layer-by-layer growth of PTO/PZO, especially for larger lattice parameters (a=4.0-4.2 Å). One option could be the new treatment developed for scandates (including NSO) by Kleibeuker et al.^{22}

Even if the extended growth of PTO/PZO superlattices remains elusive, this work may still provide useful guidelines for designing new ferroelectric materials based on B-site superlattices. This can be done either by creating new low-symmetry phases by making use of fourthorder coupling terms such as P_{001}R_{001}P_{110}R_{110}, or by inducing improper ferroelectricity via a third order term such as R_{001}P_{110}A_{¯110} or R_{110}P_{001}A_{¯110}. For example, in the latter scenario, a system with an A_{¯110} instability and an R_{001} or R_{110} rotational instability should be an improper ferroelectric. For A-site superlattices, recent work on the PTO/SrTiO_{3} system^{8} demonstrates that it is indeed possible to synthesize superlattices that induce improper ferroelectricity, even if the details of the symmetry analysis are somewhat different in that case.

V. CONCLUSIONS

In summary, we have presented the results of firstprinciples calculations of the structural phase transitions in [1:1] PTO/PZO superlattices, as well as in the parent compounds, as a function of in-plane epitaxial strain. A symmetry analysis was used to clarify the kinds of structural distortions occurring in each phase and to illuminate the nature of the transitions between phases.

While the PZO system remains in a single monoclinic phase over the whole region of lattice constants studied, and the PTO system has a single transition from a tetragonal to an orthorhombic phase, we find a significantly more complex phase sequence for the PTO/PZO superlattice. Specifically, with increasing lattice constant, we find a sequence of three transitions, first from a tetragonal to a monoclinic phase, then to a first orthorhombic phase, and finally to a second orthorhombic phase. The first orthorhombic phase occurs over such a narrow range of in-plane lattice constant that we cannot be confident of its existence in real PTO/PZO superlattices. The first three of these phases all contain strong polar and octahedral-rotation distortions, while the fourth has strong in-plane polar order, substantial antipolar distortions, and weak rotations.

Our results indicate that polarizations and rotational distortions are strongly coupled in the PTO/PZO superlattice. For example, for the monoclinic phase, in which four kinds of mode distortions are present simultaneously, we find these modes to be coupled such that the presence of any three of them will automatically induce the fourth, giving rise to a surprisingly rich behavior. Thus, we find that the inclusion of the octahedral rotations is crucial for an accurate determination of the phases and phase boundaries.

While the synthesis of short-period perovskite superlattices has the potential to lead to a wealth of new functional materials with improved properties, a thorough understanding of the interplay between epitaxial strain and spontaneous structural distortions will be needed to guide the search for such materials. The present work is a step along the way to this goal, and we hope that it might help identify new materials with enhanced ferroelectric and piezoelectric properties for next-generation applications.

Acknowledgments

The work was supported by ONR grants N00014-05-1-0054 and N00014-09-1-0300. This work is part of the research program of the Foundation for Fundamental Research on Matter [FOM, financially supported by the Netherlands Organization for Scientific Research (NWO)]. The technical assistance of V.R. Cooper in the early states of the work is gratefully acknowledged.

∗ Electronic address: j.l.blok@utwente.nl

1 B. Noheda, J. A. Gonzalo, L. E. Cross, R. Guo, S.-E. Park,

D. E. Cox, and G. Shirane, Phys. Rev. B 61, 8687 (2000). 2 G. Rijnders and D. Blank, Nature 433, 369 (2005). 3 J. Junquera and P. Ghosez, Journal of Computational and

Theoretical Nanoscience 5, 2071 (2008). 4 C. Bungaro and K. M. Rabe, Phys. Rev. B 69, 184101

(2004). 5 W. Zhong and D. Vanderbilt, Phys. Rev. Lett. 74, 2587

(1995). 6 P. Ghosez, E. Cockayne, U. Waghmare, and K. Rabe,

Phys. Rev. B 60, 836 (1999). 7 M. Fornari and D. Singh, Phys. Rev. B 63, 092101 (2001). 8 E. Bousquet, M. Dawber, N. Stucki, C. Lichtensteiger,

P. Hermet, S. Gariglio, J.-M. Triscone, and P. Ghosez, Nature 452, 732 (2008).

9 N. A. Benedek and C. J. Fennie, Phys. Rev. Lett. 106,

107204 (2011). 10 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 11 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566

(1980). 12 K. G. and H. J., Phys. Rev. B 47, R558 (1993).

13 G. Kress and J. Fürthmuller, Comput. Mater. Sci. 6, 15

(1996). 14 H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188

(1976). 15 B. P., Phys. Rev. B 50, 17953 (1999). 16 K. G. and J. D., Phys. Rev. B 59, 1758 (1999). 17 R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47,

1651 (1993). 18 O. Dieguez, K. M. Rabe, and D. Vanderbilt, Phys. Rev. B

72, 144101 (2005). 19 N. A. Pertsev, A. G. Zembilgotov, and A. K. Tagantsev,

Phys. Rev. Lett. 80, 1988 (1998). 20 H. M. Christen, L. A. Knauss, and K. S. Harshavardhan, Materials Science and Engineering B 56, 200 (1998). 21 E. D. Specht, H.-M. Christen, D. P. Norton, and L. A.

Boatner, Phys. Rev. Lett. 80, 4317 (1998). 22 J. Kleibeuker, G. Koster, W. Siemons, B. Dubbink, D. Kuiper, J. Blok, C. Yang, J. Ravichandran, R. Ramesh, J. ten Elshof, D. Blank, et al., Advanced Functional Materials 20, 3490 (2010).