Effects of anisotropic elasticity in the problem of domain formation and stability of monodomain state in ferroelectric films
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J a n Effects of anisotropic elasticity in the problem of domain formation and stability ofmonodomain state in ferroelectric films
A.M. Bratkovsky and A.P. Levanyuk
1, 2, 3 Hewlett-Packard Laboratories, 1501 Page Mill Road, Palo Alto, California 94304 Dept. Fiz. Mat. Cond., Universidad Autonoma de Madrid, Madrid 28049, Spain Moscow Institute of Radioengineering, Electronics and Automation, Moscow 117454, Russia (Dated: January 1, 2011)We study cubic ferroelectrics films that become uniaxial with a polar axis perpendicular to thefilm because of a misfit strain due to a substrate. The main present result is the analytical accountfor the elastic anisotropy as well as the anisotropy of the electrostriction. They define, in particular,an orientation of the domain boundaries and stabilizing or destabilizing effect of inhomogeneouselastic strains on the single domain state. We apply the general results to perovskite systems likeBaTiO /SrRuO /SrTiO films and find that at least not far from the ferroelectric phase transitionthe equilibrium domain structure consists of the stripes along the cubic axes or at 45 degrees tothem. We have also showed that in this system the inhomogeneous strains increase stability withregards to the small fluctuations of the metastable single domain state, which may exist not veryclose to the ferroelectric transition. The latter analytical result is in qualitative agreement withthe numerical result by Pertsev and Kohlstedt [Phys. Rev. Lett. , 257603 (2007)], but we showthat the effect is much smaller than those authors claim. We have found also that under certainconditions on the material constants, which are not satisfied in the perovskites but are not forbiddeneither, a checkerboard domain structure can be realized instead of the stripe-like one and that thepolarization-strain coupling decreases stability of a single domain state instead of increasing it.The single domain state is metastable at certain large thicknesses and becomes suitable for memoryapplications at even larger thicknesses when the lifetime of the metastable state becomes sufficientlylarge. PACS numbers: 77.80.Dj, 77.80.bn, 77.55.Px
I. INTRODUCTION
Properties of domain structures in thin ferroelectric films is currently a focus of extensive research. It is expected,quite naturally, that an understanding and an ability to control these properties will determine the prospects ofapplications of nanometer-size ferroelectrics. It depends critically on the external conditions like presence or absenceof electrodes. In this paper, we discuss domain structures in a system, which is, perhaps, the most important forapplications: a ferroelectric film with electrodes. The polar axis of the material is perpendicular to the film planeand the electrodes are ‘real’ meaning that the electric field penetrates into them, although only over tiny depths < A . This is an adequate model for the perovskite ferroelectric films on a substrate with compressive strain, likeBaTiO /SrRuO /SrTiO (BTO/SRO/STO) where the misfit strain drives the FE film into a uniaxial state. Wesupplement our analytical results with the relevant numerical estimates for BaTiO (BTO), PbTiO (PTO), andPb(Zr . Ti . )O (PZT) using the material constants available in the literatureIncomplete screening of the depolarizing field by SrRuO electrode leads to an absolute instability of a single domainstate and formation of a sinusoidal domain structure when thickness of BaTiO film is slightly above the minimalthickness compatible with ferroelectricity in this system . It seems that this situation is typical of real electrodesand we shall consider this case only. To find the minimal thickness, one does not need to take into account higherorder terms in the Landau-Ginzburg-Devonshire (LGD) free energy including the terms describing the electrostrictionsince the problem of stability of the paraelectric phase is linear . But in order to reveal the characteristics of thisstructure, i.e. to find out if the equilibrium structure is stripe-like or checkerboard and how the domain boundariesare oriented one has to take into account the anizotropy of elastic and electrostrictive properties of the ferroelectric.This is the main goal of the present paper. Specifically, we consider the case of cubic crystal anisotropy of elasticand electrostrictive properties only. This is relevant for films of cubic perovskites which become tetragonal because ofin-plane misfit strains due to cubic substrates like in the above-mentioned system. The change of cubic anisotropy totetragonal affects most strongly the dielectric properties since the crystals are “soft” dielectrically. They have muchsmaller effect on the elastic and electrostrictive properties, which can be considered to be the same as in cubic parentcrystals.Explicit account for the electrostriction and anisotropic elasticity is relevant also for study of stability of singledomain state. This has been correctly pointed out by Pertsev and Kohlstedt although these authors have missedseveral important points. Importantly, however , they overlooked that the state whose stability they were studyingwas actually metastable . Therefore, its stability with respect to small fluctuations did not mean that this statecan be used in memory applications. Indeed, its lifetime is very short if the film thickness is not sufficiently largerthan that calculated by Pertsev and Kohlstedt as the limit of single domain stability. We shall also discuss thisstability among other questions This makes sense because of several reasons. First , Pertsev and Kohlstedt performednumerical calculations using material constants for BaTiO and Pb(Zr . Ti . )O and the electrode parameters ofSrRuO while our results are analytical and apply to any material of the same symmetry. Moreover, our methodapplies to other symmetries as well. Second , Pertsev and Kohlstedt studied stability of the single domain state withrespect to “polarization wave”-like fluctuations with a single specific direction of the wave vector, while we considerwaves with k -vectors in arbitrary direction. Third , Pertsev and Kohlstedt apparently misinterpreted their own resultsby mixing together the well-known effect of homogeneous misfit strains and the effect of strains due to inhomogeneous polarization. In fact, the misfit strain simply results in renormalization of the materials constants and was effectivelytaken into account by all the previous authors. Only the account of the inhomogeneous polarization was pioneered inRef. . We show that this effect was vastly overestimated by Pertsev and Kohlstedt. In fact, the formal difference bymore than an order of magnitude between the results with and without account for electrostriction that they claimedstems from improper comparison of the compressed film with the materials constants renormalized by the misfit strainto one without any such renormalization at all, and not from the effect of the inhomogeneous strains on stability ofsingle domain state.Our analytical calculations provide a general view on the role of the inhomogeneous strains in stability of singledomain state. In particular, they reveal a possibility which seems academic at the moment but no reason is seento exclude it altogether. We mean a specific state where elasticity provokes domain formation of non-ferroelastic ferroelectric 180 ◦ domains. Such a state is realized if a certain condition on the electrostrictive and elastic constantsis met. We are not aware of an experimental realization of these conditions but we cannot find arguments prohibitingthem. It is worth mentioning that a qualitative conclusion about possibility of both stabilizing and destabilizingrole of inhomogeneous elastic strains for single domain state in ferroelectric films on substrates has been made inour previous paper where we considered an academic case of a single electrostriction constant and assumed isotropicelasticity . A surprising result of the present work is that the destabilizing effect of the inhomogeneous strains maybe very large contrary to the stabilizing one. Another unexpected result is the possibility of a checkerboard domainstate if some conditions on the material constants are met. Let us mention that without account for the anisotropicpolarization-strain coupling one comes to the conclusion about impossibility of such a state. For the perovskites thisconclusion remains valid but not in the general case.Studying the sinusoidal domain structure in BTO, PTO and PZT films on SrTiO , we find that the equilibriumorientation of the ”domain walls” is parallel (perpendicular) to the cubic axes in the film plane for BTO and PTO andis at 45 ◦ to these axes for PZT. In all cases the free energy of the sinusoidal domain structure depends very weakly onthe domain wall orientation. This is mainly due to both systems being nearly isotropic elastically and, additionally, therelevant electrostriction constant is relatively small. This observation may be important for understanding domaincreation at smallest thicknesses of the ferroelectric films. For BTO/SRO/STO system, our analytical calculationsprovide confirmation of the qualitative result of Pertsev and Kohlstedt about stabilizing effect of inhomogeneouselastic strains for single domain state in this system but with the above mentioned strong disagreement with theirstatement about importance of this effect.Having mentioned advantages and new possibilities provided by analytical calculations we should mention also theirinherent shortcomings. Our analytical method is feasible within a certain approximation only. This approximationimplies that the domain period is less than the film thickness. This condition is fulfilled for thick enough films butin very thin films the two quantities are in fact comparable. Therefore, the accuracy of our calculations should beinvestigated for these films, so the new numerical studies are desirable. We do not expect, however, that the differencebetween the results of approximated and more exact calculations either within a continuous medium theory or withinmicroscopic theories will be very large given close results of continuous and first principles theories even for films thatare just several unit cell thick (see, e.g., ).The paper is organized as follows: we describe the approximations used and define the terms in the LGD free energythat can be neglected within our approximation in Sec.II. This let us avoid unnecessary lengthy formulas in the rest ofthe paper. We spell out the constituent equations in Sec.III and then solve the general problem for the ‘polarizationwave’ (embryonic stripe domains) in the FE film with full account for elastic coupling. This is further used in Sec.Vto determine that the domain walls align with crystallographic cubic axes. Then, we find the conditions when themonodomain state first loses its stability with regards to the stripe domain structure in Sec.VI. One previouslyunexplored possibility is that the system can lose stability with regards to checkerboard domain structure, but ourresults in Sec.VII show that such a structure is absolutely unstablein perovskites although it is not necessarily so inthe general case. We summarize the present results in the Conclusions. dead/screening layer P x l P z dead/screening layer FE SUBSTRATE
FIG. 1: Schematic of the (perovskite) ferroelectric film with thickness l and metal electrodes (with screening length λ ) on amisfit substrate. The misfit makes the film a uniaxial ferroelectric with a spontaneous polarization along z − axis. II. OUTLINE OF THE METHOD AND THE APPROXIMATIONS USED
The main conclusions of this paper are made by analyzing the formula for free energy of the total system as afunction of the amplitude a of the ferroelectric ”polarization wave” presenting the sinusoidal domain structure andthe homogeneous part of the ferroelectric polarization, p . For the electrode and the film parameters of a system likeBTO/SRO/STO, the ferroelectric polarization that is perpendicular to the film plane, schematic of which is shownin Fig.1, has the form P z ( x, y, z ) = p + a cos kr cos qz, (1)where the orientation of k in the x, y plane is not fixed, 2 π/k is the period of the sinusoidal domain structure, q = π/l ,and l is the film thickness. To find the desired free energy, F ( a, p ), one has to find the elastic strains and the non-ferroelectric polarization P ⊥ = ( P x , P y ) as functions of a and p to present the total free energy as a function of a and p only. The total free energy contains contributions of the ferroelectric film, of the substrate and of the electrode.In principle, it should also contain a contribution of the voltage source but we consider here a short-circuited systemand are not concerned with the latter contribution.When calculating elastic strains in the ferroelectric which accompany the inhomogeneous polarization forming thesinusoidal domain structure, we follow the same philosophy as in our previous work. . In principle, when calculatingthese strains the inhomogeneous strains in the substrate should be taken into account. However, it is well knownthat they propagate into the substrate for about the same distances as the scale of inhomogeneity in the film ( x, y )plane. In our case, these inhomogeneities are due to the domain structure, i.e. this scale is the period of the domainstructure. Then, it is convenient to consider relatively thick films since the period of the domain structure is relativelysmall, specifically, it is much less than the film thickness , Fig. 2. The contribution of the substrate is its elasticenergy which, as we have mentioned above, is concentrated within a volume which is much smaller than the filmvolume, as defined by a small factor q/k = π/kl ≪
1. Another convenience of the thick film limit is that it is possibleto disregard the boundary conditions for the inhomogeneous parts of elastic strains and stresses at the surfaces of the
FIG. 2: (color online) Schematic of the ferroelectric film on the misfit substrate at the onset of sinusoidal polarization wave.The elastic coupling to the substrate allows inhomogeneous deformations but prohibits homogeneous strains in plane of thefilm. ferroelectric. Indeed, if we obtain a solution, which does not satisfy the boundary conditions, we can find correctionsto such a solution in a way that is conventionally used in the elasticity theory. First, we apply the external forcesto the surfaces, which are necessary to meet the boundary conditions with the strains corresponding to our solutionmaking this solution correct. Second, we apply forces opposite to the previous ones and find the strains producedby the new forces. These strains provide the correction to the original solution we were looking for. Once more, itis sufficient to observe that in our case the external forces have the period of the domain structure to understandthat the elastic energy associated with the corrections necessary to satisfy the boundary conditions can be neglectedquite similarly to the elastic energy of the substrate. Another convenience of thick film approximation is given by thepossibility to neglect those terms in the LGD free energy, which describe the electrostriction but contain componentsof non-ferroelectric polarization. According to Refs. , P ⊥ ( x, y, z ) = ( k /k ) a ⊥ sin kr sin qz, (2)where a ⊥ ≈ aq/k , q/k = π/kl ≪ . The electrostriction terms in the LGD free energy with non-ferroelectric compo-nents of polarization may contain the ferroelectric component, like P x P z u xz , or may not contain them, like in theterm P x P y u xy . In both cases, they contribute to a and p a terms in the free energy depending on a and p . In thefirst case, this contribution is proportional to ( q/k ) and in the second to ( q/k ) . Since there are also the terms a ,p a that do not contain the small factor q/k, the contribution of these terms can be neglected.Taking this into account, we write down the LGD free energy in the form: F ( P , u ik ) = F ( P ) + F ( u ik ) + F ( P , u ik ) , (3)where F ( P ) = A P z + B P z + 12 G ( ∇ ⊥ P z ) + 12 κP bz + A ⊥ P ⊥ , (4) F ( u ik ) = 12 λ (cid:0) u xx + u yy + u zz (cid:1) + λ ( u xx u yy + u xx u zz + u zz u yy ) + 2 µ (cid:0) u xy + u zy + u xz (cid:1) , (5) F ( P , u ik ) = q u zz P z + q ( u xx + u yy ) P z . (6)Here, q are the standard piezo-electric coefficients that should not be confused with the parameter q defining thetransversal profile of the polarization wave (1). In Eq.(4), A = γ ( T − T c ) , B, G = const , −→ ∇ ⊥ = ( ∂ x , ∂ y ) the gradientin the plane of the film, P bz is the non-ferroelectric (‘base’) part of the polarization perpendicular to the electrodes , A ⊥ > . Following Refs , we have neglected a term with the gradient in z − direction since it is much smaller thanthe one in plane of the film, ∂ z ≪ −→ ∇ ⊥ . It is worth mentioning that we have not included the energy of the electricfield.into the LGD free energy. The reason is that we shall use it to write down the constituent equations only.We shall eliminate u ik , P ⊥ , P bz as well the electric field components from the system of constituent electrostaticsequations to obtain two coupled equations of state for a and p . We shall obtain F ( p, a ) from the resulting equations.This is possible because of the thick films approximation. The most straightforward method to obtain F ( p, a ) wouldbe to substitute Eq.(1) into Eq.(3) supplemented by the electric field energy and to integrate over the film volume. Ingeneral, the result would not be the same as the one obtained from the constituent equations because of approximatecharacter of Eq.(1). However, for q = π/l the two results coincide and that makes it possible to use a more convenientmethod of the constituent equations. III. CONSTITUENT EQUATIONS
For the polarization components one has: AP z + BP z − G ▽ ⊥ P z + 2 q P z u zz + 2 q P z ( u xx + u yy ) = E z , (7) P bz = κE z , (8) P ⊥ = A ⊥ E ⊥ . (9)Before writing down the equations for the strain, we shall eliminate the electric field from the above three equations.Assuming Eq.(1) for P z , Eq.(2) for P ⊥ and putting : E z = E + E kz cos kr cos qz, E ⊥ = ( k /k ) E k ⊥ sin kr sin qz (10)we can replace Eqs.(7),(9) with Ap + (cid:2) BP z + 2 q P z u zz + 2 q P z u ⊥⊥ (cid:3) hom = E , (11) (cid:0) A + Gk (cid:1) a + (cid:2) BP z + 2 q P z u zz + 2 q P z u ⊥⊥ (cid:3) cc = E kz , (12) A ⊥ a ⊥ = E k ⊥ , (13)where u ⊥⊥ = u xx + u yy , [ . . . ] hom and [ . . . ] cc denote the homogeneous part ( k = 0) and the part proportional tocos kr cos qz of the expression in the brackets, correspondingly. Of course, as a result of this replacement, a part ofthe l.h.s. of Eq.(7) is lost but it corresponds to the higher harmonics of the sinusoidal distribution of the polarizationand these harmonics can be neglected close to the transition .The homogeneous part of the electric field E z can be calculated as, e.g., in Ref. yielding for the short-circuitedcase E z = − πdε b d + ǫ e l p, (14)where d is the thickness of the dead layer and ǫ e its dielectric constant. Recall that real electrodes have finite albeitsmall (Thomas-Fermi) screening length λ, which is completely analogous to a presence of the ‘dead’ non-ferroelectriclayers at the interface with thickness d/ λ. Using Eq. (14), the equation (11) gets the form: A p + (cid:2) BP z + 2 q P z u zz + 2 q P z ( u xx + u yy ) (cid:3) hom = 0 . (15)where A = A + 4 πdε b d + ǫ e l ≈ A + 4 πdǫ e l , (16)since usually the dead layer is very thin, ε b d ≪ ǫ e l. To transform Eqs. (12) we use the electrostatics equation,div D = 0 , (17)where D is the dielectric displacement, firstly for the ferroelectric material, taking into account that D =( ε ⊥ E ⊥ , ε b E z + 4 πP z ), where ε ⊥ = 1 + 4 π/A ⊥ , and ε b = 1 + 4 π/κ is the ‘ base ’ non-critical dielectric constant ,and together with the equation curl E = 0 , we find that E kz = − πq ε ⊥ k a. (18)Substituting (18) into the equation for the amplitude of the ‘polarization wave’ a, Eq.(12), we rewrite the latter as: (cid:18) A + Gk + 4 πq ε ⊥ k (cid:19) a + (cid:2) BP z + 2 q P z u zz + 2 q P z u ⊥⊥ (cid:3) cc = 0 . (19)The nontrivial solution of the above equation first appears when the coefficient in the first term in round bracketsbefore the amplitude a in the above equation first crosses zero, i.e. when A = (cid:2) − Gk − πq / (cid:0) ε ⊥ k (cid:1)(cid:3) max . Thattakes place at some A <
0, so upon lowering temperature at constant thickness the transition into domain stateoccurs somewhat below the bulk critical temperature T c , in other words. The transition for varying thickness of thefilm at some constant temperature T < T c takes place when the thickness exceeds some critical value. Therefore,the first nontrivial solution appears for the ‘polarization wave’ with the wave number k that minimizes the sum Gk + 4 πq / (cid:0) ε ⊥ k (cid:1) , so that (recall that q = π/l )4 πq ε ⊥ k = Gk , k = (cid:18) πq ε ⊥ G (cid:19) / = (cid:18) π ε ⊥ Gl (cid:19) / . (20)We can now rewrite Eq.(12) as the homogeneous one: A a + (cid:2) BP z + 2 q P z u zz + 2 q P z u ⊥⊥ (cid:3) cc = 0 , (21)where A = A + 2 Gk . (22)It is seen from Eq. (6) that the only source of elastic stresses and strains is P z ( x, y, z ) in our approximation. Since P z = p + 2 pa cos kr cos qz + a qz + cos 2 kr + cos 2 kr cos 2 qz ) , (23)we should expect that u zz = u (0) zz + u (1) zz cos 2 qz + u (2) zz cos kr cos qz + u (3) zz cos 2 kr + u (4) zz cos 2 kr cos 2 qz, (24)while for u xx , u yy , and for u ⊥⊥ we shall have similar formulas with the homogeneous part (first term in the aboveexpression) absent because of the substrate. The superscripts (0) − (4) denote contributions with different types ofthe coordinate dependencies as defined by Eq.(24). Below, we use the same superscripts for both the coefficients andthe functions.Substituting Eqs.(23),(24) and analogous equations for u xx and u yy into Eqs.(15),(21), we find: A p + (cid:2) BP z (cid:3) hom + 2 q pu (0) zz + au (2) zz ! + q a u (2) ⊥⊥ = 0 , (25) A a + (cid:2) BP z (cid:3) cc + 2 q " pu (2) zz + a u (0) zz + u (1) zz + u (3) zz u (4) zz ! + 2 q " pu (2) ⊥⊥ + a u (1) ⊥⊥ + u (3) ⊥⊥ u (4) ⊥⊥ ! = 0 . (26)We shall calculate the values u ( j ) ik in the next Section by solving the elastic problem explicitly.Importantly, the above equation of state (25) suggests that the film would tend to transform into a single domain(SD) state with p = 0 and a = 0 at temperature T SDc such that A = 0 or, in other words, A (cid:0) T SDc (cid:1) = − πd/ ( ǫ e l ) . (27)The second equation of state (26) yields a transition into a domain state ( p = 0 and a = 0) at the temperature T d such that A = 0 . or A ( T d ) = − Gk = − (cid:18) π Gε ⊥ (cid:19) / l . (28)Recall that in the present case, corresponding to BaTiO /SrRuO /SrTiO ,4 πdǫ e l > Gk ∼ πd at ε / ⊥ l , (29)where d at = √ πG ≈ A is the small ‘atomic’ length scale ( G = 0 . A for BaTiO ). The above relation means thatthe paraphase gives way to the domain phase , with a = 0 , thus preventing it from reaching the temperature T d whereit could have transformed into a single domain state. Obviously, same is true of the phase transformations in the filmas a function of thickness at constant temperature. There, one can introduce the critical thickness for domains, l d , where A = − Gk = − (cid:18) π Gε ⊥ (cid:19) / l d , (30)and the ‘critical thickness for the single domain state’ l SDc , such that A = − πd/ (cid:0) ǫ e l SDc (cid:1) . (31)These introduced critical thicknesses and temperatures are discussed in detail below in Sec.VI. IV. ELASTIC PROBLEM
Using Eqs. (5),(6), we obtain for the diagonal components of the elastic stress tensor: σ xx = λ u xx + λ ( u yy + u zz ) + q P z , (32) σ yy = λ u yy + λ ( u xx + u zz ) + q P z , (33) σ zz = λ u zz + λ ( u xx + u yy ) + q P z , (34)and formulas of the type σ xy = 2 µu xy , (35)for the off-diagonal components.We have already mentioned that the only u ik component which has a non-zero homogeneous part is u zz . This partis easily found from the condition at the free surface: σ zz = 0 at z = l/
2. From Eq.(34), one finds: u (0) zz = − q λ (cid:2) P z (cid:3) hom = − q λ (cid:18) p + a (cid:19) . (36)For the parts depending on z only , the equations of elastic equilibrium take the form: ∂σ (1) iz /∂z = 0 , (37)i.e. σ (1) iz = const = 0 since it should vanish at the free surface ( z = l/ u (1) zz = − q a / (4 λ ) , (38)and u (1) xx = u (1) yy = 0 , (39)because of the Saint-Venant’s elastic compatibility conditions for z − only dependent strains.When solving the rest of the elastic problem, we shall use the small parameter q/k ≪
1. This allows us to neglectthe derivatives with respect to z : formally, ∂/∂z ≪ ∂/∂x, ∂/∂y. As a result, the equations of the elastic equilibriumacquire the form: ∂σ (2 − xx ∂x + ∂σ (2 − xy ∂y = ∂σ (2 − yz ∂y + ∂σ (2 − zx ∂x = ∂σ (2 − yy ∂y + ∂σ (2 − yx ∂x = 0 , (40)where the superscripts (2 −
4) denote the part of stresses that is due to the three last terms in Eq.(23), which wedenote as P − z . Explicitly, λ ∂ u (2 − x ∂x + λ ∂ u (2 − y ∂y∂x + µ ∂ u (2 − x ∂y + µ ∂ u (2 − y ∂y∂x + q ∂P − z ∂x = 0 , (41) µ ∂ u (2 − z ∂y + µ ∂ u (2 − z ∂x = 0 , (42) λ ∂ u (2 − y ∂y + λ ∂ u (2 − x ∂y∂x + µ ∂u (2 − x ∂x∂y + µ ∂ u (2 − y ∂x + q ∂P − z ∂y = 0 . (43)Analogously to the isotropic case , we shall put the conditions u (2 − z = 0 that satisfy Eq.(42) but not, of course, theboundary conditions. The latter is not important in our approximation, as we argued above. Therefore, we concludethat u (2) zz = u (3) zz = u (4) zz = 0 , (44)and we are left with only two equations to solve. It is convenient to solve them separately for (2) and (3 ,
4) parts,since they correspond to different spatial harmonics.Simplifying the remaining equations (41),(43), we obtain: λ ∂ u (2 − x ∂x + ( λ + µ ) ∂ u (2 − y ∂y∂x + µ ∂ u (2 − x ∂y + q ∂P − z ∂x = 0 , (45) λ ∂ u (2 − y ∂y + ( λ + µ ) ∂ u (2 − x ∂y∂x + µ ∂ u (2 − y ∂x + q ∂P − z ∂y = 0 . (46)For the terms u (2) , we have P z ∝ cos kr , ∂ x ( y ) P z ∝ − k x ( y ) sin kr , meaning that u x ( y ) ∝ sin kr , ∂ u i /∂x = − k x u i ,etc. Then, λ k x u (2) x + ( λ + µ ) k x k y u (2) y + µk y u (2) x + q k x pa = 0 , (47) λ k y u (2) y + ( λ + µ ) k x k y u (2) x + µk x u (2) y + q k y pa = 0 . (48)The terms u (3) , (4) correspond to higher spatial harmonics in Eq. (23), but they should be taken into account since theycontribute to the terms with the main harmonic in the constituent equation (26). Since for this part P z ∝ cos 2 kr , we obtain a slightly different set of equations: λ k x u (3 , x + ( λ + µ ) k x k y u (3 , y + µk y u (3 , x + q k x a , (49) λ k y u (3 , y + ( λ + µ ) k x k y u (3 , x + µk x u (3 , y + q k y a . (50)Note that for the constituent equations we need the combinations: u (2) ⊥⊥ = u (2) xx + u (2) yy = k x u (2) x + k y u (2) y ,u (3 , ⊥⊥ = u (3 , xx + u (3 , yy = 2 k x u (3 , x + 2 k y u (3 , y . We find: u (2) ⊥⊥ = − q apf ( θ ) , (51)where θ is defined by k x = k cos θ , k y = k sin θ and the function f ( θ ) is f ( θ ) = 2 ( λ − λ ) sin θ + 2 µ cos θ ( λ + λ + 2 µ ) ( λ − λ ) sin θ + 4 λ µ cos θ . (52)For terms, corresponding to cos 2 kr and cos 2 qz, we obtain u (3) ⊥⊥ = u (4) ⊥⊥ = − q a f ( θ ) . (53)Using the results of solution of the elastic problem, Eqs.(36),(44),(51),(53), we can write Eqs.(25),(26) as A p + (cid:18) B − q λ (cid:19) p + pa (cid:20) B − (cid:18) q λ + q f ( θ ) (cid:19)(cid:21) = 0 (54) A a + a (cid:18) B − (cid:20) q λ + q f ( θ ) (cid:21)(cid:19) + ap (cid:18) B − (cid:20) q λ + q f ( θ ) (cid:21)(cid:19) = 0 . (55)From these two constituent equations corresponding to an extremum of the free energy, one can easily reconstructthe free energy: V − ˜ F ( p, a ) = A p + A a + e B p + 3 B ( θ )8 a p + 9 B ( θ )256 a , (56)where e B = B − q λ , B ( θ ) = e B + 43 (cid:20) q λ − q f ( θ ) (cid:21) , B ( θ ) = e B + 23 (cid:20) q λ − q f ( θ ) (cid:21) , (57)are the Landau coefficients before the quartic terms renormalized by the strain. The form of the free energy is thesame as in the isotropic case , but, importantly, the coefficients B and B depend on the orientation of the‘polarization wave’ given by the angle θ . V. ORIENTATION OF THE DOMAIN STRUCTURE
Consider the domain structure formed close to the paraelectric-ferroelectric transition. Although stability of theparaelectric phase is lost with respect to the polarization waves with the value of the k vector given by Eq.(20) andarbitrary orientation in the x − y plane, i.e. for any θ, the energy of sinusoidal domain structure depends on θ and onehas to find the ones corresponding to the equilibrium domain structure(s). Since we consider here one ‘polarizationwave’ only, we study the competition between the stripe-type structures. In Sec.VII we shall show that the square(checkerboard) domain structure is unstable in perovskite crystals that we study here. The checkerboard structurecould, in principle, be stable or metastable under some conditions on the material constants, but we are not aware ofany experimental example of this type, so it will be premature to study such a hypothetical case.Recall that we discuss the ferroelectric phase transition in a sample with short-circuited electrodes. Then, p = 0in the ferroelectric phase at least not far from the phase transition, and the phase transition into the inhomogeneousdomain phase occurs at A = 0. The free energy is: V − ˜ F ( p, a ) = A a + 9 B ( θ )256 a . (58)At a fixed θ, the minimum of this free energy is realized for a = − A B ( θ ) , (59)with the corresponding free energy V − ˜ F min ( p, a ) = − A B ( θ ) . (60)0We see that the equilibrium domain structure is realized for the angles θ which minimize the function B ( θ ) or,according to Eq.(57), maximize the function f ( θ ). Let us find maxima of this function. It can be written in the form: f ( θ ) = 2 λ + λ + 2 µ (cid:18) r − c tan θ + c (cid:19) , (61)where r = 2 µλ − λ , c = 4 λ µ ( λ + λ + 2 µ ) ( λ − λ ) . (62)One sees that for r > c or λ + 2 µ > λ , (63)the equilibrium domain structure corresponds to θ eq = 0 , π/ f ( θ ) max = f (0) = 1 /λ , while in the opposite case θ eq = π/ , π/ f ( θ ) max = f ( π/
4) = 2 λ + λ + 2 µ . (64)Throughout the present paper, we use the data for the material constants of BaTiO and PbTiO from Refs and of Pb(Zr . Ti . )O from Ref. . We see that for BTO and PTO the condition of Eq.(63) is met and, therefore,the equilibrium 180 ◦ domain structure consists of stripes parallel (perpendicular) to the cubic axes, while the oppositeinequality applies to PZT and the stripes make 45 ◦ with the cubic axes there. For BTO, our conclusion coincideswith that of Dvorak and Janovec who defined the equilibrium orientation of the 180 ◦ domain walls in BTO far fromthe phase transition. These authors were surprised by their conclusion about a very weak orientational dependence ofthe domain wall energy given that the experimental observations showed a clearly preferable orientation, the sameas suggested by the theory.It follows from our results that the weak orientational dependence of the domain structure energy takes place inthe sinusoidal regime too, and not only for BTO, but for all three perovskites we have made the numerical estimatesfor. Indeed, from Eq. (60) one sees that the orientational dependence of the domain structure energy comes from thefunction B ( θ ) . The maximum difference of values Eq.(57) for B ( θ ) can be used to characterize the anizotropy ofthe domain wall energy:∆ B = B − B = 23 q | f (0) − f ( π/ | = 43 q λ | λ − λ − µ | λ + λ + 2 µ . (65)We found ∆ B /B ∼ × − for BTO, much smaller anizotropy ∼ × − for PTO, and even smaller one for PZTwhere ∆ B /B ∼ × − (we have used the parameters listed in the footnote ). Such a weak angular dependenceof the domain structure energy is in accordance with phase field results of Ref. that shows domain walls mainlywith thermodynamically favorable orientations but also with strong deviations from them. VI. LOSS OF STABILITY OF A SINGLE DOMAIN STATE
It is convenient to study the loss of stability of the single domain state with respect to formation of a domainstructure using Eq. (56). With this, we mean the loss of stability with respect to an arbitrarily small ‘polarizationwaves’ so that the original single domain state may be, in principle, either stable or metastable. Specifically, in ourcase, when Eq. (29) is valid, this state is metastable .A solution of the equations, ∂ ˜ F /∂p = 0 , ∂ ˜ F /∂a = 0 , (66)corresponding to a single domain state ( p = 0, a = 0) is possible only if A < p extr = − A / e B , where thesubscript stands for the ‘extremum’. This extremum is a minimum (which is relative in our case) if ∂ ˜ F /∂p > ∂ ˜ F /∂a > , (67)1 FIG. 3: Schematic of the sinusoidal domain structure in (a) BTO and PTO, and (b) PZT. While the stripes are oriented alongthe crystal axes in case (a), in case (b) the stripes are at 45 degrees with respect to cubic axes, the difference being due to anopposite sign of elastic anizotropy in those two cases. p = p extr , a = 0 given that ∂ ˜ F /∂p∂a is evidently zero at this point. The first inequality in Eq.(67) isobviously valid for A <
0, while the validity of the second is not immediately evident.We find from (56):4 ∂ ˜ F∂a ! a =0 ,p = p extr = A + 3 B ( θ ) p extr = A − A − A (cid:16) B ( θ ) − e B (cid:17) / e B, (68)From the condition (cid:16) ∂ ˜ F /∂a (cid:17) a =0 ,p = p extr = 0 , we obtain the value of A corresponding to a loss of stability of the singledomain state with respect to appearance of a polarization wave with a given orientation, A pw ( θ ). It is convenient topresent it in the form: A pw ( θ ) = − πdε b d + ǫ e l + Gk + (cid:18) πdε b d + ǫ e l − Gk (cid:19) β ( θ ) , (69)where β ( θ ) = q − q f ( θ ) λ e Bλ + 2 [ q − q f ( θ ) λ ] . (70)The last term in Eq.(69) is the result of the polarization-strain coupling while the first two present the prior casewithout this coupling . According to Eq.(67), the corresponding single domain state will be (meta)stable at lowtemperatures such that A < min A pw ( θ ) . The actual loss of stability of the single domain state corresponds to the minimum of A pw ( θ ) . We have seen inSec.V that for perovskites the angular dependencies are very weak and we can neglect it putting f λ = 1 . Then β = q − q e Bλ + 2 ( q − q ) (71)Since in BTO, PTO, and PZT q > q , the last term in Eq.(69) is positive and, therefore, the region of meta stability of the single domain state in these systems is broader than according to and , in apparent accordancewith . However, there are serious reservations. First of all, the effect is not very spectacular. Indeed, the factor β inthe last term of Eq.(69) is always less than one half, β < / , approaching that value when q − q → ∞ . Therefore, A pw < − πdε b d + ǫ e l , (72)where the r.h.s. corresponds to a very strong strain coupling. Recall that A = − πd/ ( ε b d + ǫ e l ), or A = 0 , corresponds to what was calculated in several papers as a ‘ critical thickness of single-domain ferroelectricity ’ l SDc ,Eq.(31).To get the opposite limit of a weak strain coupling for A pw (0) , we put q − q = 0 and neglect Gk , as Pertsevand Kohlstedt did. We see that − πdε b d + ǫ e l < A pw (0) < − πdε b d + ǫ e l , (73)i.e. because of account for the polarization-strain coupling the value of A pw (0) changes always by less than 1 . ε b d < ǫ e l this is also the interval of change of the thickness corresponding to absolute loss ofstability of the single domain state at a fixed temperature. The above moderate, less than 50%, range of change is in striking disagreement with a statement by Pertsev andKohlstedt who claimed a more than an order of magnitude change by nullifying the electrostrictive constants. Theydo not report the details of their procedure, but it is clear from the rest of the paper that their suggestion of puttingthe electrostrictive constants to zero implied changes in the coefficients of the LGD free energy that should have beenrenormalized by the misfit strains. Evidently, it has nothing to do with the effects of the polarization-strain couplingomitted in Refs , since this renormalization is automatically taken into account there, while the effect of misfit strainon LGD coefficients was apparently neglected in a gedanken exercise performed in Ref. .Specifically, we find that for the perovskites BTO and PTO β = 0 .
4, while in PZT this parameter is 0 . , i.e. fourtimes smaller. We see that BTO and PTO are similar and closer to the limit q − q → ∞ , β = 0 .
5, i.e. the pointof loss of stability of single domain state is quite close in these materials to the ’critical thickness of single-domainferroelectricity’ studied in Ref. and elsewhere. However, the latter does not have any practical importance because3 FIG. 4: (color online) Regions of (meta)stability of the single domain and polydomain states in the ferroelectric film as afunction of temperature T at fixed thickness (top) and as a function of the film thickness l at fixed T (bottom). Upon loweringthe temperature at a fixed thickness l (top panel), the paraphase gives way to domains that are stable at all temperatures T < T d , where T d is below the critical temperature of the bulk ferroelectric transition T c . The single domain (SD) state is metastable at low temperatures
T < T
SDms , when a strain coupling is neglected. The strain coupling shifts the boundary ofmetastability towards so-called critical temperature for a single domain state, T SDc , as marked by the vertical arrows for theperovskites in question. The phase behavior of the films as a function of their thickness l at fixed temperature (bottom) isqualitatively similar. Very thin films are in a paraelectric phase that is replaced by domains at larger thicknesses l > l d . Thesingle domain state is metastable at thicknesses l > l
SDms and becomes suitable for memory applications at even larger (yet tobe determined) thicknesses when the life time of the metastable state becomes sufficiently large. Strain coupling may extendthe boundary of the metastability down to the so-called ‘critical thickness for SD ferroelectric state’ l SDc . The films with FEmemory would be somewhere at larger thicknesses in the metastable region, as marked on the diagrams. if one fixes the temperature and reduces the film thickness starting with a monodomain ferroelectric state at lowtemperatures or large thicknesses, that state will give way to domains before the thickness determined by the limit ofthe single domain state stability is reached. The fact of the matter is that the single domain state is metastable, itmay have a large lifetime at low temperatures and large film thicknesses but this lifetime goes essentially to zero (to‘atomically’ short times) when the above mentioned temperature or thickness are reached.Importantly, it follows from Eq.(69) that if q < q , the account for the inhomogeneous strains shrinks theregion of metastability of the single domain state. This shows that contrary to the claim by Pertsev and Kohlstedtthere is no general physical phenomenon such as stabilization of a single domain state because of inhomogeneousstrains accompanying formation of domains This may seem surprising, because solids are known to ‘dislike’ theinhomogeneous strains (free energy usually goes up). Moreover, the expectation of Pertsev and Kohlstedt is justified fora free-standing film, at least for elastically isotropic solid . But it is not certain for a film on substrate considered bothby them and in the present work. To explain the physical reason, we recall that the coupling with strain renormalizesthe coefficients before fourth-order terms in the LGD free energy, in our case we mean the coefficients before p and p a terms. Then, one has to take into account that the homogeneous strains in the plane of the substrate arenot possible while inhomogeneous ones are. Both homogeneous and inhomogeneous polarization create homogeneousstrain but to a different extent see Eqs.(36) and (38).while inhomogeneous strains are created, of course, by the4 -250-200-150-100-50050100 lmslcld FE film thickness l (nm) T ( o C ) SD metastable paraphase BTO/SRO/STO domains M e m o r y ? FIG. 5: (color online) The phase diagram for BaTiO /SrRuO /SrTiO films in the plane temperature-thickness. The linemarked l d delineates the para- and domain phases, while the one marked l ms marks the boundary of the metastability regions ofthe single domain state calculated for BTO accounting for strain coupling, Eq. (69). The broken line at larger thicknesses from l ms marks the metastable region calculated without accounting for the strain coupling, Eq.(76). The films with FE memorywould be somewhere at larger thicknesses in the metastable region. inhomogeneous polarization only. The out-of-plane and in-plane strains couple with the ferroelectric polarization P z by electrostriction terms with different coefficients and the final result of renormalization of the coefficient of p a term is due to several contributions and it is not clear upfront. It should be obtained by a consistent analysis, asit has been done above. No reason is seen to discard the possibility that the inequality q < q can be realized insome systems, and one cannot exclude, at least for the moment, the possibility of favoring the multidomain state bythe polarization-strain coupling. Interestingly enough, this favoring may be very strong: according to Eq.(69), theincrease of the region of absolute instability of single domain state becomes infinite when q − q tends to e Bλ / /SrRuO /SrTiO system are shown inthe temperature-film thickness ( T − l ) plane in Fig.5. They are found from the conditions that we discussed aboveand write down here for a reference: A = − Gk = − (cid:18) π Gε ⊥ (cid:19) / l d , (74) A = − πdǫ e l SDc , (75) A = − πdǫ e l SDms , (76)where the Landau coefficient A is evaluated for a given temperature T of interest, d/ λ = 0 . A, ǫ e = 8 .
45 for5SrRuO electrode, G = 0 . , and ε ⊥ the dielectric constant [see its definition below Eq.(17)] in the plane of theFE film has been found from the Landau coefficients . The last condition corresponds to l = l SDms found without accounting for the strain coupling. The critical line T − l SDms in the phase diagram, Fig.5 for BTO with an account forstrain coupling has been found from Eq.(69). The arrows on Fig.5 show the evolution of the state at either T = constor l = const.One should understand that in both illustrations it is implied that the corresponding critical points have physicalvalues as solutions to the conditions ( ?? ), or, equivalently, Eqs.(28),(27),(30),(31). Consider first the lowering of thetemperature at a fixed thickness l (Fig. 4, top panel), where the paraphase transforms into domain state below thetemperature T d that is smaller than the critical temperature of the bulk ferroelectric transition T c . We see that thesingle domain (SD) state would be metastable at low temperatures
T < T
SDms in the region overlapping with thedomain state. Note that the T SDms plotted in Fig.5 is found without accounting for the strain coupling. The straincoupling then shifts the boundary of metastability towards the so-called critical temperature for a single domain state T SDc thus broadening the range of metastability of the SD state, as shown in Fig.4. The phase behavior of the filmsas a function of their thickness l at fixed temperature (Fig. 4, bottom panel) is qualitatively similar. Very thin filmsremain in a paraelectric phase that is replaced by the domains at larger thicknesses l > l d . We see that both T SDc and l SDc are actually unreachable in the present case, since the system may get to those points only by moving fromthe paraphase down (right to left on the phase diagram, Fig.5), but such transitions are preempted by the domaininstability that sets in first. The single domain state is metastable at thicknesses l > l
SDms , and becomes suitable for memory applications at even larger (yet to be determined) thicknesses where its life time becomes sufficiently long. VII. INSTABILITY OF THE CHECKERBOARD DOMAIN STRUCTURE
In the previous Sections, we have assumed that the domain structure is stripe-like by taking into account onlyone ‘polarization wave’. This and other possibilities have been studied by Chensky and Tarasenko who consideredthe uniaxial ferroelectric isotropic in the x − y plane. Along with the stripe-like structure they discussed also thecheckerboard and the hexagonal domain structures. The latter can be realized in the presence of an external field onlywhich is not discussed in this paper. However, a checkerboard structure should be analyzed as an alternative to thestripe structure. In Ref. , the authors stated that the checkerboard structure never realizes, although, surprisingly,there is no proof of this statement. In this Section, we shall show that this structure is indeed unstable for theisotropic case treated in Ref. and then show that this conclusion holds also when one explicitly takes into accountthe polarization-strain interaction , apart from mere renormalization of the LGD coefficients by the misfit strains.Once again, we consider a short-circuited sample, i.e.the ferroelectric polarization is described by: P z = a cos k r cos qz + a cos k r cos qz, (77)where k and k are two noncollinear vectors whose modulus is given by Eq.(20) and whose (mutually orthogonal)directions remain unspecified for a moment, Fig. 6. A. Checkerboard domains without elastic coupling ( q = q = 0) In this case the solution for the fields is [cf. Eq.(10)] E z = E k z cos k r cos qz + E k z cos k r cos qz, (78) E x,y = E k x,y sin k r sin qz + E k x,y sin k r sin qz, (79)with Eq.(18) still applicable to spatial harmonics (as follows from linearity of Maxwell equations) and the equationof state for the fundamental harmonics is the same as Eq.(12): A a + (cid:2) BP z (cid:3) cc = E k z , (80)where one retains the terms (cid:2) BP z (cid:3) cc ∝ cos k r cos qz (symmetry dictates the analogous expressions for a ). In theabove equation, P z = ( a cos k r + a cos k r ) cos qz = 916 (cid:0) a + 2 a a (cid:1) cos k r cos qz + 916 (cid:0) a + 2 a a (cid:1) cos k r cos qz + . . . , (81)so that we obtain for the fundamental harmonic the following equations of state: A a + 9 B (cid:0) a + 2 a a (cid:1) = 0 (82)6 FIG. 6: Schematic of the checkerboard domain structure. It is absolutely unstable in case of BaTiO , PbTiO , andPb(Zr . Ti . )O typical perovskite ferroelectrics. and the analogous equation for a . Since these equations are obtained from extremum of the free energy, ∂ e F /∂a =0 , we again restore the full free energy, accounting for the symmetric contribution by the a harmonic: V − e F = A a + a ) + 9256 B (cid:0) a + a (cid:1) + 964 Ba a . (83)The equations of state have the checkerboard solution: a = a = − A / B. (84)Checking what type of extremum for the free energy is this solution, ∂ F∂a × ∂ F∂a − (cid:18) ∂ F∂a ∂a (cid:19) = − A < , we see that the checkerboard solution is the maximum of the free energy for some directions in the a , a plane andis absolutely unstable . B. Checkerboard domains with elastic coupling
We have seen above that the elastic coupling renormalizes the fourth order coefficients in formulas like Eq. (83)reducing them by some amounts which are different for different coefficients. Thus, instead of Eq.(83), we will have: V − e F = A a + a ) + 9256 (cid:0) B a + B a (cid:1) + 964 B a a , (85)where B and B are given by Eq.(57) for the corresponding angles and B is a new coefficient which depends on theboth angles and which can be, in principle, both positive and negative. Both from Eq.(57) and the cubic symmetry7one realizes that B = B = B . For what follows, it is important to mention that when B is negative it cannot beof large absolute value, otherwise there will be directions in the ( a , a ) plane along which the free energy diminisheswithout limits at large values of a , a what means a global instability of the system. By putting a = a , one seesfrom Eq. (85) that to avoid this instability the condition B + 2 B > B > a = a = − A B + 2 B ) . (87)To analyze stability of this solution, we calculate the second derivatives ∂ F∂a = ∂ F∂a = A B a + 932 B a = − A B B + 2 B ,∂ F∂a ∂a = 916 B a a = ± A B B + 2 B we then find the discriminant Z ≡ (cid:18) ∂ F∂a (cid:19) (cid:18) ∂ F∂a (cid:19) − (cid:18) ∂ F∂a ∂a (cid:19) = 14 A B − B B + 2 B . (88)Our further study is aimed at finding out if and when the condition of positiveness of Z , i.e. B > B , is compatiblewith the two conditions of the global stability mentioned above. Thus we need formulas for the coefficients B and B .Turning to taking into account explicitly the polarization-strain coupling, we recall that in our approximation ofsufficiently thick film the only source of the elastic strains is P z . This function contains now a cross term stemmingfrom P z = ( a cos k r cos qz + a cos k r cos qz ) (89)= (cid:2) a + a + 2 a a (cid:0) cos p + r + cos p − r (cid:1) + a cos 2 k r + a cos 2 k r (cid:3) qz , (90)where p ± = k ± k . Naturally, the components of the strain tensor will have terms depending on cos p + r , cos p − r :.We will have for u xx u xx = u (0) xx + u (1) xx cos 2 qz + u p + xx cos p + r + u p − xx cos p − r + u q + xx cos p + r cos 2 qz + u q − xx cos p − r cos 2 qz + u (3)1 ,xx cos 2 k r + u (3)2 ,xx cos 2 k r + (cid:16) u (4)1 ,xx cos 2 k r + u (4)2 ,xx cos 2 k r (cid:17) cos 2 qz, and the analogous equations for u yy and u zz . From our previous experience, it becomes immediately clear that u p ± xx ( yy ) = u q ± xx ( yy ) , since the equations for thosecomponents do not depend on z and we can now drop the indices p and q from the corresponding terms, leaving only u + , − xx ( yy ) . Also, due to the same reason as above, u (0) xx ( yy ) = u (1) xx ( yy ) = 0. Then, one can write: u xx = (cid:0) u + xx cos p + r + u − xx cos p − r (cid:1) (1 + cos 2 qz )+ u (3)1 ,xx cos 2 k r + u (3)2 ,xx cos 2 k r + ( u (4)1 ,xx cos 2 k r + u (4)2 ,xx cos 2 k r ) cos 2 qz. and a similar equation for u yy . The ‘diagonal’ terms for the first (second) k harmonic are u (3)1(2) ,xx + u (3)1(2) ,yy = u (4)1(2) ,xx + u (4)1(2) ,yy = − q a f (cid:0) θ (cid:1) . λ ∂ u ± x ∂x + ( λ + µ ) ∂ u ± y ∂y∂x + µ ∂ u ± x ∂y + q ∂P ± z ∂x = 0 , (91) λ ∂ u ± y ∂y + ( λ + µ ) ∂ u ± x ∂y∂x + µ ∂ u ± y ∂x + q ∂P ± z ∂y = 0 , (92)where P ± z = a a cos p ± r , in components u ± x ( y ) ∝ sin p ± r , ∂ u ± i /∂x = − ( p ± x ) u ± i etc. λ p ± x u ± x + ( λ + µ ) p ± x p ± y u ± y + µp ± y u ± x + q p ± x a a = 0 , (93) λ p ± y u ± y + ( λ + µ ) p ± x p ± y u ± x + µp ± x u ± y + q p ± y a a = 0 . (94)Then, u ± xx + u ± yy = p ± x u ± x + p ± y u ± y = − q a a f ( θ ± ) . (95)For u zz we conclude, as above, that only u (0) zz and u (1) zz are non-zero and following the same reasoning as for the stripephase, we obtain u (0) zz = u (1) zz = − q λ ( a + a ) . Having solved the elastic problem, we are now in a position to write down the constituent equations containing a and a only. To this end, we write two equations for a and a analogous to Eq.(21) but this time [ . . . ] cc would meanthe proportionality to cos k r cos qz or cos k r cos qz, respectively. Since both equations have the same structure, wewill discuss that for a only and, for the sake of brevity, we will mention only the terms containing a , the other termsare the same as for the one-sinusoid case discussed above.For clarity sake, we repeat Eq.(21) with a minor change for the present case: A a + (cid:2) BP z + 2 q P z u zz + 2 q P z ( u xx + u yy ) (cid:3) cc = 0 . (96)It is straightforward to find that the a -containing term stemming from (cid:2) P z (cid:3) cc is 9 a a / . From Eq.(26), one seesthat [ P z u zz ] cc = a (cid:16) u (0) zz + u (1) zz (cid:17) , recall that now we consider the case p = 0, and the contribution of this term is − q λ a a . Now [ P z ( u xx + u yy )] cc = a (cid:16) u (3)1 ,xx + u (3)1 ,yy (cid:17) + a (cid:16) u (4)1 ,xx + u (4)1 ,yy (cid:17) + 34 a ( u + xx + u + yy + u − xx + u − yy )and contribution of this term to the equation of state is − q a a (cid:2) f ( θ + ) + f ( θ − ) (cid:3) . Finally, the constituent equation for a takes the form: A a + 916 a B ( θ ) + 98 a a B (cid:0) θ + , θ − (cid:1) = 0 , where we have introduced B (cid:0) θ + , θ − (cid:1) = B − q λ − q (cid:2) f ( θ + ) + f ( θ − ) (cid:3) (97)9Similarly to the case of one sinusoid, we recover the free energy: F = A (cid:0) a . + a (cid:1) + 9256 B ( θ ) a + 9256 B ( θ ) a + 964 B (cid:0) θ + , θ − (cid:1) a a , (98)Recall that the square symmetry suggests that B ( θ ) = B ( θ ). and B ( θ ) is given by Eq. (57)Turning to examining the sign of the discriminant Z , we should mention that according to Eqs.(57) and (97): B − B = − B + 43 q (cid:20) f ( θ + ) + f ( θ − ) − f ( θ ) (cid:21) . One sees that the maxima of Z correspond to maxima of f ( θ + ) and f ( θ − ) [note that f ( θ +max ) = f ( θ − max ) because ofthe cubic symmetry], which are, automatically, the minima of f ( θ ) as we have seen in Sec.V. Then,[ B − B ] max = − B + 43 q (cid:20) f ( θ max ) − f ( θ min ) (cid:21) . Using the values of f ( θ max ) and f ( θ min ) found in Sec.V, we find that if λ + 2 µ > λ , [ B − B ] max = − B + 43 q (cid:18) λ − λ + λ + 2 µ (cid:19) = − B + 43 q λ + 2 µ ) + λ λ ( λ + λ + 2 µ ) (99)and if λ + 2 µ < λ , [ B − B ] max = − B + 43 q (cid:18) λ + λ + 2 µ − λ (cid:19) = − B + 43 q λ − λ − µ λ ( λ + λ + 2 µ ) (100)To prove that the checkerboard structure can be stable, at least in principle, with respect to small fluctua-tions, we should demonstrate that the positiveness of [ B − B ] max is compatible with the conditions B > B − B ] min > λ + 2 µ ≃ λ , and q = 0 . Both Eq.(99) and Eq.(100) then give[ B − B ] max ≃ − B + 2 q λ , (101)and give for the positiveness of [ B − B ] max the same condition: q > Bλ / B > q < Bλ / . One sees that for a nearly elastically isotropic ferroelectric with q = 0 , the checkerboard structure isat least metastable if 3 Bλ / > q > Bλ / B − B ] max ≃ − e B − λ (cid:0) q − q (cid:1) < . Indeed, we have already mentioned above that for the perovskites q > q . Also, e B > absolutely unstable
VIII. CONCLUSIONS
With the use of the Landau-Ginzburg-Devonshire theory, we have studied the effects of polarization-strain couplingwhen defining the character of equilibrium domain structures and the limits of absolute instability of a single domainstate in thin films of cubic ferroelectric films on a misfit substrate. On the compressive substrate, the cubic ferroelectric0behaves substantially as a uniaxial ferroelectric with the polar axis perpendicular to the film. The film is sandwichedbetween the electrodes that do not provide a perfect screening of the depolarizing field because of the finite Thomas-Fermi screening.length. Such a system is exemplified by (100) BaTiO /SrRuO /SrTiO film and similar perovskitestructures. Quantitative results have been obtained for BaTiO , PbTiO , and Pb(Zr . Ti . )O . We have foundthat close to the paraelectric-ferroelectric phase transition or at the film thicknesses close to the minimal thicknesscompatible with the ferroelectricity, the equilibrium domain structure in perovskites is the stripe-wise one withthe stripes parallel (perpendicular) to the cubic axes in BaTiO , PbTiO , while running at 45 ◦ to cubic axes inPb(Zr . Ti . )O . The energy of the domain structure depends very weakly on the stripe orientation, the maximumchange proves to be well below 1% in all three cases. We found that because of the polarization-strain couplinga competing checkerboard domain structure may, at least in principle, be equilibrium or metastable when certainconditions on the material constants are met, but we are not aware of any material system with such conditions.The limit of absolute instability of the single domain state changes due to the polarization-strain coupling. Thus, theinterval where the absolute instability is absent, meaning a metastability in the cases at hand, widens in perovskitesin agreement with the earlier conclusion by Pertsev and Kohlstedt . However, this effect is much smaller thanthat claimed by them. Increase of the metastability range is substantial in BaTiO and PbTiO where the absoluteinstability limit becomes close to what is often called the ”critical thickness for ferroelectricity” l SDc , Fig.5, butwithout accounting for the domain formation, which sets in first and prevents the system from ever reaching thispoint. The effect is much smaller in Pb(Zr . Ti . )O . We have found also that the polarization-strain coupling canlead to narrowing of the region of relative stability of the single domain state under certain conditions on the materialconstants, but we are not aware of an experimental.realization of these conditions.APL has been partially supported by Ministry of Science and Education of Russian Federation (State Contract D.J. Kim, J.Y. Jo, Y.S.Kim, Y.J. Chang, J.S. Lee, J.-G. Yoon, T.K. Song, and T.W. Noh, Phys. Rev. Lett. , 237602(2005). Y.S.Kim, D.H.Kim, J.D.Kim, Y.J.Chang, T.W.Noh, J.H.Kong, K.Char, Y.D.Park, S.D.Bu, J.-G.Yoon, J.-S Chung, Appl.Phys. Lett. , 102907 (2005). Y.S. Kim, J.Y. Jo, D.J. Kim, Y.J. Chang, J.H. Lee, and T.W. Noh, T.K. Song, J.-G. Yoon, J.-S. Chung, S. I. Baik andY.-W. Kim, C.U. Jung, Appl. Phys. Lett. , 072909 (2006). J. Junquera and P. Ghosez, Nature , 506 (2003). A.M. Bratkovsky and A.P. Levanyuk, Appl. Phys. Lett. , 253108 (2006). A.M. Bratkovsky and A.P. Levanyuk, J. Comput. Theor. Nanosci. , 465 (2009); arXiv:0801.1669v4 [cond-mat.mtrl-sci]. P.Aguado-Puente and J.Junquera, Phys. Rev. Lett. , 177601 (2008). E.V. Chensky and V.V. Tarasenko, Sov. Phys. JETP , 618 (1982) [Zh. Eksp. Teor. Fiz. , 1089 (1982)]. N.A. Pertsev and H. Kohlstedt, cond-mat/0603762. N.A. Pertsev and H. Kohlstedt, Phys. Rev. Lett. , 257603 (2007). N.A. Pertsev and H. Kohlstedt, Phys. Rev. Lett. , 149702 (2008). A.M. Bratkovsky and A.P. Levanyuk, Phys. Rev. Lett. , 149701 (2008). A.M. Bratkovsky and A.P. Levanyuk, Phil. Mag. , 113 (2010). A.K. Tagantsev and G. Gerra, J. Appl. Phys. , 051607 (2006). N.A.Pertsev, A.G.Zembilgotov, and A.K.Tagantsev, Phys. Rev. Lett. , 1988 (1998). Y.Li, S.Y. Hu, Z.K. Liu, and L.Q. Chen, Acta Mater. , 395 (2002). G. Sheng, J.X. Zhang, Y.L. Li, S. Choudhury, Q.X. Jia, Z.K. Liu, and L.Q. Chen, J. Appl. Phys. , 054105 (2008). J. Hlinka, Ferroelectrics , 132 (2008); J. Hlinka and P. M´arton, Phys. Rev. B , 104104 (2006); P. Marton, I. Rychetsky,and J. Hlinka, Phys. Rev. B , 144125 (2010). V. Dvorak and V. Janovec, Jpn. J. Appl. Phys. , 400 (1965). J.Fousek and M.Safrankova, Jpn. J. Appl. Phys. , 403 (1965). In cubic crystals the usual notations for the elastic constants are λ = c , λ = c , µ = c . We have used the followingvalues of the parameters in the SI units:(i) for BaTiO c = 1 . × Nm − , c = 8 . × Nm − , c = 1 . × Nm − , q = 1 . × JmC − , q = − . × JmC − , ˜ B = 3 . × Jm C − , G = 0 . A , Ref. (cf. );(ii) for PbTiO : c = 1 . × Nm − , c = 7 . × N m − , c = 1 . × Nm − , q = 1 . × JmC − , q = 4 . × JmC − , ˜ B = 2 . × Jm C − ;(iii) for Pb(Zr . Ti . )O : c = 1 . × Nm − , c = 8 . × Nm − , c = 3 . × Nm − , q = 7 . × JmC − , q = − . × JmC − , ˜ B = 1 . × Jm C −4