Abstract. The influence of the flexoelectric and rotostriction coupling on the phase diagrams of ferroelastic- quantum paraelectric SrTiO3 films was studied using ...
Submitted to the Journal of Applied Physics
Roto-flexoelectric coupling impact on the phase diagrams and pyroelectricity of thin SrTiO3 films A.N. Morozovska1,2, E.A. Eliseev3, S.L. Bravina1, A.Y. Borisevich4, and S.V. Kalinin5 1
Institute of Physics, National Academy of Science of Ukraine, 46, pr. Nauki, 03028 Kiev, Ukraine
2
Institute of Semiconductor Physics, National Academy of Science of Ukraine, 45, pr. Nauki, 03028 Kiev, Ukraine
3
Institute for Problems of Materials Science, National Academy of Science of Ukraine, 3, Krjijanovskogo, 03142 Kiev, Ukraine 4
Materials Sciences and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA 5
Center for Nanophase Materials Science, Oak Ridge National Laboratory, Oak Ridge, TN, 37831, USA
Abstract The influence of the flexoelectric and rotostriction coupling on the phase diagrams of ferroelasticquantum paraelectric SrTiO3 films was studied using Landau-Ginzburg-Devonshire (LGD) theory. The phase diagrams in coordinates temperature - film thickness were calculated for different epitaxial misfit strains. Tensile misfit strains stimulate appearance of the spontaneous out-of-plane structural order parameter (displacement vector of an appropriate oxygen atom from its cubic position) in the structural phase. Compressive misfit strains stimulate appearance of the spontaneous in-plane structural order parameter. Gradients of the structural order parameter components, which inevitably exist in the vicinity of film surfaces due to the termination and symmetry breaking, induce improper polarization and pyroelectric response via the flexoelectric and rotostriction coupling mechanism. Flexoelectric and rotostriction coupling results in the rotoflexoelectric field that is antisymmetric inside the film, small in the central part of the film, where the gradients of the structural parameter are small, and maximal near the surfaces, where the gradients of the structural parameter are highest. The field induces improper polarization and pyroelectric response. Penetration depths of the improper phases (both polar and structural) can reach several nm from the film surfaces. An improper pyroelectric response of thin films is high enough to be registered with planar-type electrode configurations by conventional pyroelectric methods. 1
Introduction Thin ferroic films, their multilayers and heterostructures attract permanent fundamental scientific interest and contain the great potential for multiple conventional and novel applications in nanoelectronics, sensorics and optics [1, 2, 3]. Unique functionalities of thin films are believed to originate primary from the surfaces and interfaces influence on the physical properties via great versatility of the coupling mechanisms enhanced or induced by the gradient effects and symmetry breaking [ 4 , 5 ]. Since purely empirical way of the film properties study is extremely timeconsuming to explore fundamental mechanisms and novel applications, theoretical studies of the surface and interface induced couplings are necessary for effective search of thin films with improved functional properties. Elastic field gradients can induce polarization near the surfaces and interfaces via the flexoelectric coupling [6, 7, 8]. Note, that all materials are flexoelectrics [9], and all perovskite oxides with static rotations (such as octahedral rotations also called antiferrodistortions) possess rotostriction and rotosymmetry [ 10 ]. The joint action of flexoelectric effect and rotostriction (previously named roto-flexo-effect) can lead to spontaneous polarization at an interface across which the octahedral rotations changes [11, 12]. In particular, it was predicted that the roto-flexo field at oxide interfaces gives rise to improper ferroelectricity and pyroelectricity in the vicinity of antiferrodistortive boundaries and elastic twins in bulk SrTiO3 [11] as well as in the vicinity of semi-infinite SrTiO3 surface below 105 K [12]. In the work we theoretically study the influence of the roto-flexoelectric coupling on the phase diagrams of thin SrTiO3 films. Several mechanisms of the surface influence on the phase transitions in ferroelastic-quantum paraelectric films were included in our calculations, namely misfit strain between the epitaxial film and substrate, surface energy, rotostriction and electrostriction effects, and flexoelectric coupling. All the mechanics can essentially influence on the phase diagrams of thin SrTiO3 films. It is well-known that misfit strain can induce ferroelectric polarization in quantum paraelectric SrTiO3 [13, 14, 15, 16]. Within LGD phenomenology the surface energy coefficients determine so-called extrapolation lengths [17, 18]. Extrapolation lengths in dependence on their sign and value, which in turn are interface- and termination-dependent [19, 20 ], can enhance or suppress polarization via the spatial confinement of the film interfaces. Rotostriction and electrostriction effects lead to renormalization (in particular enhancement) of the biquadratic coupling between the structural order parameter (displacement vector of an appropriate oxygen atom from its cubic position in the oxygen octahedra) [13]. We demonstrated that the flexoelectric and rotostriction coupling, which are the strongest in the gradient regions of thickness
2
~5 nm under the film surfaces, induces improper polarization and pyroelectric response of thin SrTiO3 films.
2. Problem statement and basic equations Here we consider a thin epitaxial ferroelastic-quantum paraelectric film deposited on a rigid substrate. The film structure with capacitor geometry shown in Fig.1a corresponds to the shortcircuited electrical boundary conditions. The structure shown in Fig.1b corresponds to the opencircuited electrical boundary conditions. The misfit strain u11 = u22
= um originates from the
x 3 =0
lattice mismatch between the epitaxial film and substrate.
(a)
top electrode
x2
P2 A1
film
x1
substrate h/2 electrode
x2
A2
P1
-h/2 x1
(b)
P3 A3
substrate x3
x3
Figure 1. Ferroelastic film clamped on a rigid substrate. The film structure (a) corresponds to the short-circuited electrical boundary conditions; (b) corresponds to the open-circuited electrical boundary conditions. In the parent high temperature phase above the structural phase transition the spontaneous order parameters are absent and the free energy has the form:
F = ∫ d 3r f b (r ) + ∫ d 3r f S (r ) V
(1a)
S
The bulk energy density has the form [11]:
3
2 ⎛ ⎛ ⎞ ⎜ β (T )P 2 + β P 2 P 2 + ... + 1 g ⎜ ∂Pi ∂Pk ⎟ − P E − ε 0 ε b E i ij i j ijkl ⎜ ∂x ∂x ⎟ i i ⎜ i 2 2 ⎝ j l ⎠ ⎜ ⎜ ⎛ ∂P cijkl ∂u ij ⎞ 1 ⎟ + ξ ijkl Pi Pj Ak Al u ij u kl + f ijkl ⎜ u ij k − Pk f b = ⎜ − qijkl u ij Pk Pl + ⎜ ⎟ x x 2 2 ∂ ∂ ⎜ l l ⎠ ⎝ ⎜ vijkl ⎛ ∂Ai ∂Ak ⎞ ⎜ 2 2 2 ⎟−r u A A ⎜ ⎜⎜ + α i (T )Ai + α ij Ai A j + ⎜ ∂x ∂x ⎟ ijkl ij k l 2 l ⎠ ⎝ j ⎝
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎟ ⎠
(1b)
Pi is the polarization vector, uij (x ) is the strain tensor (i,j = 1 − 3); gijkl and vijkl are the gradient coefficients; f ijkl is the forth-rank tensor of flexoelectric coupling experimentally determined for SrTiO3 by Zubko et al [21], qijkl is the forth-rank electrostriction tensor, rijkl is the rotostriction tensor, cijkl is the elastic stiffness; ξijkl is the biquadratic coupling term [ 22 , 23 ]. Ai are the components of the structural order parameter (OP) (see e.g. Ref. [13]). Note, that for the case of the oxygen octahedra rotation transition as an order parameter could be chosen either the rotation angle or the displacement of an appropriate oxygen atom from its cubic position [24]. The temperature dependence of the coefficients α i and β i can be fitted with Barrett law [ 25 ],
(
(
)
(
αi (T ) = αT Tq( E ) coth Tq( E ) T − coth Tq( E ) T0( E )
)) ,
(
(
)
(
βi (T ) = βT Tq( Φ ) coth Tq( Φ ) T − coth Tq( Φ ) TS
))
for
quantum paraelectrics. The higher expansion coefficients βij and αij are weakly temperature dependent. In the dielectric limit (i.e. under the absence of screening by internal carriers) the depolarization electric field components E id are determined self-consistently from Maxwell equation ε 0 ε b
∂Eid ∂xi
=−
∂Pi ∂xi
, where ε0=8.85×10−12 F/m is the universal dielectric constant, ε b is the
“base” isotropic lattice permittivity, different from the ferroelectric soft mode permittivity [26, 27, 28]. It should be noted that the flexoelectric coupling could be rewritten in the form of Lifshitz
(
)
invariant, fijkl uij ∂Pk ∂xl − Pk ∂uij ∂xl 2 [7, 11-12, 29], or in the short form f ijkl uij both forms give the same equations of state, the choice between f ijkl uij
∂Pk . Though ∂xl
∂Pk and Lifshitz invariant ∂xl
has crucial implications on the boundary conditions [12, 30]. The surface energy density has the form:
f S = aijS Ai A j + bijS Pi Pj
(1c)
4
Note, that the surface energy coefficients aijS and bijS should be termination dependent. This immediately leads to the termination dependent extrapolation lengths for the polarization and structural order parameters. It is in agreement with numerous DFT studies for termination effects in ferroelectrics [31, 32, 33, 34]. After appropriate variation, surface energy (1c) contributes to the
⎛ ⎛ ⎞ ∂P ∂A ⎞ boundary conditions ⎜⎜ 2aijS A j + n j vijkl k ⎟⎟ = 0 and ⎜⎜ 2bijS Pj + n j g ijkl k + n j f klij u kl ⎟⎟ = 0 , ∂xl ∂xl ⎠ S ⎝ ⎝ ⎠S where n j is the outer normal components to the film surfaces. The equation of state for strain and stress (generalized Hooke’s law) is:
σ ij = cijkl u kl + f ijkl
∂Pk ∂xl
− qijkl Pk Pl − rijkl Ak Al
(2)
2.2. One-dimensional (1D) case Below we consider one-dimensional (1D) case, when the structural order parameter Ai ( x3 ) and polarization Pi ( x3 ) depend on the distance from the interface x3 . Elastic fields in the epitaxial thin film clamped on a rigid substrate have the view: σ33 = 0 ,
σ13 = 0 ,
σ12 = −q44 P1P2 − r44 A1 A2 ,
σ 23 = 0 ,
(3a)
σ 22 = (c11 + c12 )um + c12u33 + f12
∂P3 − q12 P32 + P12 − q11 P22 − r12 A32 + A12 − r11 A22 , ∂x3
(3b)
σ11 = (c11 + c12 )um + c12u33 + f12
∂P3 − q12 P32 + P22 − q11 P12 − r12 A32 + A22 − q11 A12 ∂x3
(3c)
u 33 = −2
(
)
(
)
(
)
(
)
c12 f ∂P 1 um + q11 P32 + q12 P22 + q12 P12 + r11 A32 + r12 A22 + r12 A12 − 11 3 , c11 c11 c11 ∂x3
(
)
u11 = u22 = um ,
u13 =
q 44 P1 P3 + r44 A1 A3 2c 44
−
f 44 ∂P1 , 2c 44 ∂x3
u12 = 0 , u23 =
q44 P2 P3 + r44 A2 A3 f ∂P2 − 44 . 2c44 2c44 ∂x3
(3d) (3e) (3f)
We regard that the bulk material in the parent high temperature phase has m3m symmetry in order to derive Eqs.(3). Substituting the elastic solution Eq.(3) into equations of state δf b δ Pi = 0 and δf b δ Ai = 0 , we derived the free energy functional with renormalized coefficients: ~ ~ ~ ~ f b = f hom + f grad + f electric
(4a)
5
(
) + β (P
(
)
(
)(
) P (A
~ * f hom ≅ β*3 T , um P32 + β*33 P34 + β13 P12 + P22 P32 + β1* T , um P12 + P22
) ( ) +A ) + ξ (P A + P A ) + ξ (P + P )A + ξ P P A A + ξ (P A + P A )P A + α (T , u )A + α A + α (A + A )A + α (T , u )(A + A ) + α A A + α (A * P12 P22 + β12 * 12
2 1
* 3
* 11
2 2
m
4 1
2 2
2 3
* P12 A12 + P22 A22 + ξ*33 P32 A32 + ξ*31 + P24 + ξ11
2 1
* 33
* 13
4 3
2 1
* 13
2 2
2 3
2 1
2 2
44 1 2
2 3
1
* 1
* 44
2
2 1
m
2 3
2 1
2 2
1 1
2
2
3
2 2
* 12
2 1
2 2
(4b)
3
* 11
4 1
+ A24
⎞ g ⎛ ∂P ⎞ 2 f q ⎟ + ⊥ ⎜ 3 ⎟ + 11 12 (P 2 + P 2 ) ∂P3 1 2 ⎟ 2 ⎜⎝ ∂x3 ⎟⎠ ∂x3 c11 ⎠ 2 2 2 ∂P2 ⎞ v 44 ⎛⎜ ⎛ ∂A1 ⎞ ⎛ ∂A2 ⎞ ⎞⎟ v11 ⎛ ∂A3 ⎞ f 44 q 44 ⎛ ∂P1 ⎟+ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ + P2 + P3 ⎜⎜ P1 ⎟ 2 ⎜ ⎜ ∂x ⎟ + ⎜ ∂x ⎟ ⎟ + 2 ⎜ ∂x ⎟ ∂ ∂ x x c 44 3 3 ⎠ ⎝ ⎝ 3⎠ ⎝⎝ 3 ⎠ ⎝ 3 ⎠ ⎠ ⎞ ∂P ⎛f r ⎛ ∂P ∂P ⎞ f r f r + ⎜⎜ 11 12 (A12 + A22 ) + 11 11 A32 ⎟⎟ 3 + 44 44 A3 ⎜⎜ A1 1 + A2 2 ⎟⎟ ∂x3 ⎠ c11 c 44 ⎝ ∂x3 ⎠ ∂x3 ⎝ c11
2 g || ⎛⎜ ⎛ ∂P1 ⎞ ⎛ ∂P2 ~ ⎜ ⎟ +⎜ f grad ≅ 2 ⎜ ⎜⎝ ∂x3 ⎟⎠ ⎜⎝ ∂x3 ⎝
⎞ ⎟ ⎟ ⎠
)
2
(4b)
Note, that the form of the energy (4) coincides with the one obtained by the substitution of Eq.(3) ~ into Eq.(1b) and making a Legendre transformation f b = f b − σ33u33 − 2σ13u13 − 2σ 23u23 . The expansion coefficients of the structural order parameter in Eq.(4b) are renormalized by rotostriction coupling ( rij ) and misfit strain u m that are introduced as:
⎛ ⎛ c ⎞ c ⎞ α1* T , u m = α1 (T ) − ⎜⎜ r12 + r11 − 2r12 12 ⎟⎟u m , α *3 T , u m = α 1 (T ) − 2⎜⎜ r12 − r11 12 ⎟⎟u m c11 ⎠ c11 ⎠ ⎝ ⎝
(5a)
r122 r122 r11r12 r442 r112 * * * α = α11 − , α12 = α12 − , α13 = α12 − − , α 33 = α11 − 2c11 c11 c11 2c44 2c11
(5b)
(
)
(
)
* 11
The biquadratic coupling terms ξ ij are renormalized by rotostriction ( rij ) and electrostriction ( qij ) coupling: * ξ11 = ξ11 −
r12 q12 , c11
* ξ12 = ξ *21 = ξ12 −
r12 q12 , c11
ξ *44 = ξ 44 − * ξ13 = ξ12 −
r44 q 44 , c 44
r11q12 , c11
ξ*33 = ξ11 − ξ*31 = ξ12 −
r11q11 , c11
r12 q11 . c11
(6a) (6b)
The expansion coefficients of the polarization are renormalized by electrostriction coupling ( q ij ) and misfit strain u m :
⎛ c ⎞ β1* T , u m = β1 (T ) − ⎜⎜ q12 + q11 − 2q12 12 ⎟⎟u m , c11 ⎠ ⎝
(
⎛ c ⎞ β*3 T , u m = β1 (T ) − 2⎜⎜ q12 − q11 12 ⎟⎟u m , c11 ⎠ ⎝
)
q122 β = β11 − , 2c11 * 11
q122 β = β12 − , c11 * 12
(
)
2 q11q12 q44 β = β12 − − , c11 2c44 * 13
q112 β = β11 − . 2c11 * 33
(7a)
(7b)
The gradient coefficients ( g ij ) are renormalized by the flexoelectric coupling ( f ij ) as: [7]
6
⎛ f2⎞ g || = ⎜⎜ g 44 − 44 ⎟⎟ , c 44 ⎠ ⎝
⎛ f2⎞ g ⊥ = ⎜⎜ g11 − 11 ⎟⎟ . c11 ⎠ ⎝
(7c)
For the short-circuited boundary conditions (i.e. when electric potential is zero at the conducting
(
electrodes), the depolarization field is E = − P3 ( x3 ) − P3 d 3
)
1 +h / 2 dx3 P3 ( x3 ) . For ε 0 ε b , where P3 = h − h∫/ 2
the open-circuited boundary conditions (i.e. when electric displacement is zero outside the film) depolarization field is E 3d = − P3 ( x3 ) ε 0 ε b . For the open-circuited boundary conditions the electric ε ε E2 ~ energy f electric = − Pi Ei − 0 b should be modified by adding the term Di E i [35] (see also [36]). 2
Variation of the bulk free energy (1b) and surface energy (1c) densities with respect to the structural order parameter Ai ( x3 ) and polarization Pi ( x3 ) leads to equations of state supplemented by the mixed boundary conditions
⎛ S ∂A ⎞ ⎜ 2a1 A1 + n3 v 44 1 ⎟ = 0, ⎜ ∂x3 ⎟⎠ S ⎝
⎛ S ∂A ⎞ ⎜ 2a1 A2 + n3 v 44 2 ⎟ = 0, ⎜ ∂x3 ⎟⎠ S ⎝
⎛ S ∂A ⎞ ⎜ 2a3 A3 + n3 v11 3 ⎟ = 0 (8a) ⎜ ∂x3 ⎟⎠ S ⎝
⎛ S ⎞ ∂P f ⎜ 2b1 P1 + n3 g || 1 + n3 44 (q 44 P1 P3 + r44 A1 A3 )⎟ = 0 ⎜ ⎟ ∂x3 c 44 ⎝ ⎠S
(8b)
⎛ S ⎞ ∂P f ⎜ 2b1 P2 + n3 g || 2 + n3 44 (q 44 P2 P3 + r44 A2 A3 )⎟ = 0 ⎜ ⎟ ∂x3 c 44 ⎝ ⎠S
(8c)
⎛ S ⎛ ⎞ ∂P c ⎜ 2b3 P3 + n3 g ⊥ 3 + 2n3 ⎜ f12 − 12 f11 ⎟u m + ⎜ ⎟ ⎜ c11 ∂x3 ⎝ ⎠ ⎜ ⎜ n f 11 q P 2 + q P 2 + q P 2 + r A 2 + r A 2 + r A 2 11 3 12 2 12 1 11 3 12 2 12 1 ⎜ 3c 11 ⎝
(
⎞ ⎟ ⎟ ⎟ =0 ⎟ ⎟ ⎠S
)
(8d)
The components of the outer normal to the film surfaces x3 = m h 2 are n3 = ±1 correspondingly. Hereinafter we assumed that the tensors aijS ≡ δij aiS and bijS ≡ δij biS are diagonal and, more important, chosen the same for both film surfaces for simplicity only. The assumption leads to the symmetric (with respect to the interfaces x3 = m h 2 ) boundary conditions (8) and thus to the symmetric solutions. In general case different interface energy should lead to the asymmetric boundary conditions. It may seem that one can neglect all nonlinear terms in the boundary conditions (8b)-(8d) for estimations, but exactly some of the terms couple polarization with structural order parameter and thus can strongly affect phase diagrams and finite size effects in thin multiferroic films.
7
3. Approximate analytical solution in decoupling approximation 3.1. Structural order parameter components
Numerical analyses showed that in the structural phase (e.g. octahedrally tilted phase α 1* < 0 ), polarization value, orientation and corresponding derivatives weakly influence on the structural order parameters Ai in such materials as perovskites SrTiO3, CaTiO3, EuTiO3. For the case all polarization-dependent terms can be omitted in equations for the structural order parameters, and, in this decoupling approximation, approximate analytical solution acquires the form:37
⎛ ⎞ x3 2m + c0 m ⎟ , sn⎜ ⎟ 1 + m ⎜⎝ L1A T , u m 1 + m ⎠
A1 ( x3 ) = Ab
(
A2 ( x3 ) = 0 .
)
(9)
Here sn(u|m) is the elliptical sine function. For the case of identically zero or negligibly small A3 , * and characteristic length L1A (T , u m ) = the amplitude Ab = − α 1* 2α 11
− v 44
(
2α 1* T , u m
)
, which should
be positive in the structural phase. For the case A3 ≠ 0 , one can use the approximations Ab ≈ −
* 2α1* (T , u m ) − (α *3 α13 α *33 )
4α − (α * 11
*2 13
α
* 33
)
and L1A (T , u m ) ≈
(
− v 44
) (
* 2α T , u m − α *3 α 13 α *33 * 1
)
, which are self-
* consistent with Eq.(12) under the condition α 13 A12 + α *3 < 0 .
The constants m and c0 should be determined from the symmetric boundary condition (8) at the film surfaces x3 = ± h 2 . Symmetry of the boundary conditions gives c0 = K (m ) 2 , where K(m) is the complete elliptic integrals of the first kind [38]. Similar solution is obtained with the substitution of subscripts: 1↔2. At m→0 the term A1 (x3 ) vanishes, corresponding to the onset of the structural phase transition into either cubic phase or another structural phase. The condition of A1 ( x3 ) vanishing is given by equation ⎛ h cos⎜⎜ A ⎝ 2 L1 T , u m
(
)
⎞ ⎛ λ S1 h ⎟= ⎜ sin A A ⎟ L T,u ⎜ 2L T , u 1 m m ⎠ ⎝ 1
(
)
where the structural order parameter extrapolation length λ S 1 =
(
v 44 2a1S
)
⎞ ⎟, ⎟ ⎠
(10)
is introduced. The way of
derivation of expression (10) is explained in Ref.[ 39 ]. The critical film thickness
(
h = 2L T , u m A cr
A 1
)
(
⎛ L1A T , u m arctan⎜⎜ ⎝ λ S1
) ⎞⎟ below which the structural order parameters ⎟ ⎠
A1, 2 disappear can
be determined from Eq.(10) at given temperature T, misfit strain u m and extrapolation length λ S 1 8
The solution different from (9) is A1 ( x3 ) = A2 (x3 ) = A( x3 ) ≈ Ac
⎛ 2m x3 K (m ) ⎞⎟ + sn ⎜ A m , ⎟ 1 + m ⎜⎝ L1 (T , um ) 1 + m 2 ⎠
(
(11)
)
* * * + α12 + α13 where the amplitude Ac = − α1* 2α11 for identically zero or negligibly small A3 . For
A3 ≠ 0 the amplitude Ab ≈ −
(4α
* 2α 1* (T , u m ) − (α *3 α 13 α *33 )
* * ) − 2(α13*2 α *33 ) + 2α12 + 2α13
* 11
, which is self-consistent with
* Eq.(12) under the condition 2α 13 A12 + α *3 < 0 . Length L1A (T , u m ) is the same as in Eq.(9). At m→0
A1 (x3 ) vanishes, which means the structural phase transition. The boundary of the structural phase A1 ( x3 ) = A2 (x3 ) = A( x3 ) is also given by Eq.(10). Typically the gradient coefficient v44 >> v11 [ 40]. As a result, the characteristic lengths
(
)
(
)
L1A T , u m = L2A T , u m , which determine the gradient scale of the structural order parameter
components A1 and A2, are typically higher than lattice constant(s), but the length L (T , u m ) = A 3
v11
12
2α *3 (T , u m )
that determines the gradient scale of the structural order parameter
component A3, is (much) smaller than the lattice constant, thus approximate analytical solution for A3 is:
(
)
⎧ α * A 2 (x ) + A 2 (x ) + α * * 3 ⎪ − 13 1 3 * 2 3 A12 + A22 + α *3 < 0 , α13 A3 ( x3 ) ≈ ⎨ 2α 33 ⎪ * 2 2 * ⎩0, α 13 A1 + A2 + α 3 > 0
(
(
)
(12)
)
Note, that the accuracy of approximate solution (12) is the best for the case a 3S = 0 , when
∂A3 ∂x3
=0
at the film surfaces x3 = m h 2 . To be self-consistent hereinafter we regard a 3S = 0 . The boundary * (A12 (x3 ) + A22 (x3 )) + α *3 = 0 . of the structural phase with A3 ( x3 ) ≠ 0 is given by equation α 13
3.2. Polarization components
Relatively simple results for polarization can be derived for the open-circuited electrical boundary conditions, since P3 ( x3 ) appeared negligibly small due to the strong depolarization field.
Moreover, corresponding spatial scale, L (T , u m ) = P 3
g⊥
2β *3 (T , u m ) + 1 ε 0 ε b
12
, appeared ~0.1 nm, i.e.
significantly smaller than the lattice constant (~0.4 nm) of considered perovskites. Thus putting P3 ( x3 ) ≈ 0 , we should solve two coupled equations for polarization components P1 and P2 only. 9
In the tetragonal-like structural, A1 ( x3 ) = A2 (x3 ) ≡ A( x3 ) , P1 (x3 ) = P2 ( x3 ) ≡ P(x3 ) , where P (x3 ) obeys equation:
(
)
* * γ ( x3 , A, A3 )P + 4β11 + 2β12 P 3 − g ||
((
∂2P ∂x
2 3
= E||FR ( A, A3 )
)
)
* * * + ξ12 A 2 ( x3 ) + ξ13 A32 ( x3 ) The function γ ( x3 , A, A3 ) = 2β1* + ξ 44 A 2 + 2 ξ11
E||FR ( A, A3 ) =
(13) and flexo-roto field
f 44 r44 ∂( A3 A) depends on the structural parameter components. The structural c 44 ∂x3
parameter components A1 = A2 ≡ A( x3 ) and A3 (x3 ) are given by Eq.(9) and (12) correspondingly. To be consistent with approximation (12) we should put ⎛ ⎞ f r ∂P conditions to Eq.(13) is ⎜⎜ P ± λ P1 ± 44S 44 AA3 ⎟⎟ ∂x3 2b1 c 44 ⎝ ⎠
∂ ( A3 A) ∂x3
≈ A3
= 0 , where λ P1 = x3 = ± h 2
∂A . The boundary ∂x3 g || 2b1S
is extrapolation
length of the polarization. In the monoclinic-like structural phase, A2 ( x3 ) = 0 and A1 ( x3 ) ≠ 0 (or A1 ( x3 ) = 0 and A2 ( x3 ) ≠ 0 ), P1 (x3 ) ≠ 0 and P2 (x3 ) = 0 (or P1 ( x3 ) = 0 and P2 ( x3 ) ≠ 0 ), where P1 ( x3 ) obeys equation: γ ( x3 , A1 , A3 )P1 + 4β P − g || * 3 11 1
The
function
E1FR ( A1 , A3 ) =
∂ 2 P1 ∂x
(
2 3
= E1FR ( A1 , A3 )
)
* * γ (x3 , A1 , A3 ) = 2 β1* + ξ11 A12 ( x3 ) + ξ13 A32 ( x3 )
and
the
(14) flexo-roto
field
f 44 r44 ∂ ( A3 A1 ) are introduced. The structural parameter components A1 and A3 are c 44 ∂x3
given by Eq.(10) and (12) correspondingly. To be consistent with approximation (12) we should put
∂( A3 A1 ) ∂x3
≈ A3
⎛ ⎞ f r ∂P . The boundary conditions to Eq.(14) are ⎜⎜ P1 ± λ P1 ± 44S 44 A1 A3 ⎟⎟ ∂x3 ∂x3 2b1 c 44 ⎝ ⎠
∂A1
The characteristic lengths L (T , u m ) = L (T , u m ) = P 1
P 2
g ||
2β1* (T , u m )
= 0. x3 = ± h 2
12
, which determine the
gradient scale of the polarization components P1 and P2, are typically higher than lattice constant(s). The fact allows us to solve Eqs.(13), (14) along with extrapolation length-dependent boundary conditions. Solution of Eqs.(14) can be found numerically, but approximate analytics is possible for determination of the ferroelectric and polar phase boundaries [as argued in the Supplementary Materials ]. The following transcendental equation determines the in-plane ferroelectric phase
boundary in dependence on film thickness h, misfit strain um and temperature T: 10
h 2 ⎛ h 2 γ ( x, Ai (x )) ⎞ γ (h 2, Ai (h 2)) ⎛⎜ γ ( x, Ai ( x )) ⎞⎟ ⎜ ⎟ cos ∫ sin ∫ dx dx = λ P1 (15) ⎜0 ⎟ ⎜0 ⎟ − − g g − g || || || ⎝ ⎠ ⎝ ⎠ The boundary (15) depends on the extrapolation length λ P . The function γ ( x, Ai ( x )) is different for
the structural phases A1 ( x3 ) = A2 (x3 ) ≡ A( x3 ) and A2 ( x3 ) = 0 , A1 ( x3 ) ≠ 0 (compare Eq.(13) and (14) correspondingly). The critical thickness hcrP , below which the switchable ferroelectric polarization components disappear can be determined from Eq.(15) at given temperature T, misfit strain u m and extrapolation length λ P . The boundary of the local polar phase induced by the flexo-roto field is diffuse and can be estimated from the condition E FR ( A1 , A2 , A3 ) ≈ 0 indicating the flexo-roto field vanishing far from the film surfaces. For the short-circuited electric boundary conditions out-of-plane polarization component appears and it is coupled with corresponding flexo-roto field component
E3FR =
f11 r11 2 ⎞ ∂ ⎛ f11 r12 2 2 ⎜ A A A3 ⎟⎟ (see Eq.(A.2g) in Suppl. Mat. ). For the case of short+ + 1 2 c11 ∂x3 ⎜⎝ c11 ⎠
(
)
circuited conditions the critical thickness of the out-of-plane switchable ferroelectric polarization appearance can be determined from equation: ⎛ h 2 γ(x ) ⎞ ⎛ h 2 γ(x ) ⎞ γ (h 2) ⎜ ⎟ dx + λ P 3 cosh ∫ sinh⎜ ∫ dx ⎟ = 0 , ⎜ ⎟ ⎜ ⎟ g g g ⊥ ⊥ ⊥ ⎝0 ⎠ ⎝0 ⎠ 1 1 where γ ( x ) = 2β*3 + 2(ξ*33 A32 + ξ*31 (A12 + A22 )) + and typically γ ( x ) ≈ . ε 0εb ε 0εb
(16)
4. Phase diagrams 4.1. Approximate expressions for the phase boundaries
In order to calculate phase diagrams of thin films, of interest are the surface energy coefficients a iS ≥ 0 and biS ≥ 0 . These parameters uniquely determine corresponding extrapolation lengths
λ P1 =
g|| S 1
2b
,
λ P3 =
g⊥ S 3
2b
,
λ S1 =
v 44 S 1
2a
,
λ S3 =
v11 2a3S
,
(17)
The physical meaning of these parameters is discussed in detail in Refs.[17, 18]. Since the value and temperature dependence of the surface energy coefficients (and thus extrapolation lengths) are priory unknown for a given ferroic material, we can vary it in the actual range of 1 – 100 lattice constants. Analytical calculations, based on the free energy (1) minimization by direct variational method with trial functions taken in the form of elliptic (when depolarization field is absent) or hyperbolic (when depolarization field is present) functions [39], after the integration on x3 over the 11
film thickness h, leads to the renormalization of the expansion coefficients for the structural order parameter components:
(
)
(
π 2 ν 11
)
α T , um , h ≈ α T , um + R 3
(
* 3
)
(
)
α T , um , h ≈ α T , um + R 1
* 1
π 2 hλ S 3 + 2h 2 π 2 ν 44
π hλ S 1 + 2h 2
2
≡ α(T ) +
≡ α(T ) +
⎛ c12 ⎞ ⎜ ⎟u m , − − 2 r r 12 11 ⎜ ⎟ c + 2h 2 11 ⎠ ⎝
π 2 ν11 π 2 hλ S 3
(18a)
⎛ c ⎞ − ⎜⎜ r12 + r11 − 2r12 12 ⎟⎟u m , (18b) c11 ⎠ π hλ S 1 + 2h ⎝ π 2 ν 44
2
2
The expansion coefficients for the in-plane and out-of-plane polarization components are:
(
)
(
)
β1R T , u m , h ≈ β1* T , u m +
(
)
(
π 2 g || π hλ P1 + 2h 2
)
β 3R T , u m , h ≈ β *3 T , u m +
2
≡ β(T ) +
π 2 g ||
⎛ c ⎞ − ⎜⎜ q12 + q11 − 2q12 12 ⎟⎟u m , c11 ⎠ π hλ P1 + 2h ⎝ 2
2
⎛ g⊥ g⊥ c ⎞ ≡ β(T ) + − 2⎜⎜ q12 − q11 12 ⎟⎟u m . (λ P3 + L⊥ )h ⎝ (λ P 3 + L⊥ )h c11 ⎠
(18c)
(18d)
Note that Eq.(18b) is valid for the short-circuited boundary conditions and L⊥ = ε 0 ε b g ⊥ ~0.2 nm. For the open-circuited boundary conditions the renormalized coefficient β3R (T , um , h ) ≈ β*3 (T , um ) +
1 π2 g ⊥ + 2 2ε 0εb π hλ P 3 + 2h 2
(18e)
is typically positive due to the big positive additive (2ε0εb ) originated from the depolarization −1
field. Using Eqs.(18) allows the boundary between the structural or polar phases and “high temperature” phase with zero A1,2, A3, P1,2 and P3 to be estimated from conditions α 1R (T , u m , h ) = 0 ,
(
)
(
)
α 3R T , u m , h = 0 , β1R T , u m , h = 0
and β 3R (T , u m , h ) = 0 correspondingly.
The
equilibrium
boundaries between different structural and/or polar phases should be determined numerically from the condition of the corresponding energies equality for a given misfit strain u m .
4.2. Simulation results 4.2.1. Phase diagrams in coordinates “temperature –SrTiO3 film thickness”
The phase diagrams in coordinates “temperature –SrTiO3 film thickness” are shown in Figs 2-5. All diagrams are calculated for the short-circuited electric boundary conditions, which are the
most favorable for P3 appearance and most typical for thin film applications as e.g. corresponding to the thin film placed between conducting electrodes (planar capacitor). Different plots (a)-(d) on each diagram are calculated for different misfit strains u m = 0, -0.5, -2, 1 %. The phase labels indicate nonzero components of the structural order parameter and ferroelectric polarization components. All polar phases correspond to switchable ferroelectric polarization, while the flexoroto fields (proportional to the gradients of the structural parameter components) can induce 12
improper polarization components, which will be demonstrated in the next section. For instance, designation P1P2A3-phase means the phase with nonzero ferroelectric polarization components P1, P2 and the structural order parameter A3; A3-phase means the phase with nonzero structural order parameter A3; HT-phase means the parent phase with identically zero Ai and Pi entire the film. Solid lines denote the phase transitions between the ferroelectric, structural and HT phases. Phases with improper (and thus non-switchable) polarization components induced by the flexo-roto fields are designed as IPi.. For instance, an improper IP3-phase can be induced by the flexo-roto field component E3FR =
f11 r11 2 ⎞ ∂ ⎛ f11 r12 2 2 ⎜ A A A3 ⎟⎟ (see Eq.(A.2g) in Suppl. Mat. ) as well as by + + 1 2 c11 ∂x3 ⎜⎝ c11 ⎠
(
)
⎞ ⎛ c inhomogeneity in the boundary conditions 2⎜⎜ f12 − 12 f11 ⎟⎟um (see Eq.(A.3c) in Suppl. Mat. ) and c11 ⎠ ⎝ likely to be observed experimentally when the flexoelectric coupling is strong. The term
⎞ ⎛ c 2⎜⎜ f12 − 12 f11 ⎟⎟u m leads to the transformation of the HT-phase into the IP3-phase for nonzero um . c11 ⎠ ⎝ Similarly, phases with surface-induced structural order parameter components are designed as IAi.. For instance, designation A3P3-IA1P2 phase means the phase with nonzero ferroelectric polarization P3 and spontaneous structural parameter A3, and with components A1 and improper polarization P2 induced by the surface influence. Let us underline that improper (surface) I-phases primary appear due to the flexo-roto coupling. For the case of zero misfit strain ( u m = 0), which corresponds to the film on a matched substrate, phase diagram contains the maximal amount of phases including the ferroelectricstructural P1P2A3-phase (see plots (a) in Figs 2-5). Compressive (i.e. negative) misfit strain ( um 0.05%) leads to the emergence of in-plane A1,2-phases and P1,2-phases as well as to the simultaneous disappearance of out-of-plane A3-phases and P3-phases (see plots (d) in Figs 2-5). The conclusion regarding the possibility to trigger between in-plane and out-of-plane components of spontaneous polarization and structural order parameter by the sign of misfit strain is in agreement with the phase diagram calculated by Pertsev et al. [13]. Note, that the phase diagram [13] is independent on the film thickness, since it was calculated without taking into account the gradient effects (and so without flexoelectric coupling). The gradient effects can be minimized in the limiting case λ P, S → ∞ (i.e. aiS = biS = 0 ), that corresponds to the so-called “natural” boundary conditions, when the surface energy is independent on the structural order parameter and polarization and all gradients are zero inside a homogeneous film. The natural 13
boundary conditions correspond to the constant solution, Pi ( x3 ) = Pi 0 and Ai ( x3 ) = Ai 0 , independent of the film thickness. For the case of six constants, Pi 0 and Ai 0 are determined from the minimum of the bulk density (4b), and we checked that our results reproduce those in [13]. The gradient and flexoelectric effects contribution to the phase diagrams inevitably appears for finite extrapolation lengths, at that the values of critical thicknesses of the structural order parameter ( hcrA ) and polarization ( hcrP ), phase boundaries position and amount of phases strongly depend on the value of polar ( λ P ) and structural ( λ S ) extrapolation lengths and relations between them. In particular, the critical thicknesses hcrA and hcrP decrease with extrapolation length increase. To demonstrate the fact, Figs 2-5 are calculated for different extrapolation lengths: small λ P, S ~lattice constant, high λ P, S >>lattice constant, equal λ S = λ p and different λ S ≠ λ p ).
The phase diagrams in Figure 2 are calculated for equal small extrapolation lengths
λ P = λ S =0.4 nm. Phase boundaries are strongly dependent on the film thickness and temperature for the case. The regions of the structural and polar phases increase with the film thickness increase. The boundary of parent HT-phase is almost independent on the film thickness. The critical thicknesses hcrA and hcrP are maximal for the case of small extrapolation lengths, but strongly dependent on misfit strain. For most cases hcrA ≤ hcrP . The phase diagrams in Figure 3 are calculated for high enough equal extrapolation lengths
λ P = λ S =40 nm. Some phase boundaries are close to the horizontal lines for the case. In the sense the results tends to the limiting case λ P, S → ∞ . However size effect is still pronounced at least for small misfits (see plot (a)). The phase diagrams in Figure 4 and 5 are calculated for different extrapolation lengths λ P ≠ λ S . The cases λ P >> λ S and λ P > λ S ferroelectric phases exist in the essentially wider temperature range than in the case λ P >h can be registered under quasi-uniform temperature distribution entire the film thickness and gives information about the pyroelectric activity integrated over the film thickness. In the case of thin films (shown in Figs.9a,b), which are inversely polarized ( P3 ( z ) = − P3 (− z ) ), the pyroelectric response measured in
a sandwich electrode configuration should be negligibly small. For thicker films when λT < h/2 at f < 100 MHz it is possible to register a pyroelectric response from near-surface unipolarized region connected with P3 component (see Figs. 9c,d). Comparison of Fig.9c and 9d shows that a strong pyroelectric asymmetry can be observed in thick films. At that, the asymmetry is not related with different extrapolation lengths (they are set equal for z = ± h 2 , i.e. λ P1 (h 2) = λ P1 (− h 2) , λ P 3 (h 2) = λ P 3 (− h 2) ), but originated from misfit strain relaxation in thick enough films. For a film thickness more than several hundreds of nanometres misfit strain exists only in the vicinity of the clamped interface z = + h 2 , and it vanishes when approaching mechanically free surface z = − h 2 . The asymmetry in pyroactivity along film thickness allows registration of the total pyroelectric response by integral pyroelectric methods and pyroactivity probing by modulation methods. The pyroelectric coefficients Π 1, 2 caused by in-plane polarization components P1,2 can be measured in planar electrode configurations [46], which give the information about pyroactivity arisen due to existence and interaction between the components P1 and P2. Such type of measurements were considered theoretically in [47] along with analysis of primary, secondary and tertiary pyroelectric effect contributions. Feng et al. [47] studied pyroelectric response of LiTaO3 type II sensor with planar electrode configuration.
23
1
0.6
Πi (10-4 C/Km2)
0
Π3
-1
Πi (10-4 C/Km2)
(b) um=0.05%, Π1 (a) um=0.05%,
T=20 K, h=20 nm -2
A1P2-IA3P1P3-phase
T=60 K, h=20 nm
0.3
Π3
0
Π2 Π1
-0.3
A1-IA3P1P3-phase
Π2
-3 -10
0
-0.6
10
coordinate (nm)
0.6
0
-10
2
Π1
0 -500
Πi (10-4 C/Km2)
Πi (10-4 C/Km2)
(c) um=0.05%, T=60 K, h=1000 nm 0.3
10
coordinate (nm)
Π1
Π3 -490
coordinate (nm)
0
Π3
-2
(d)
-4 -480
480
490
coordinate (nm)
500
Figure 9. Distributions of pyroelectric coefficient components in a short-circuited SrTiO3 film of
thickness h=20 nm (a,b) and 1000 nm ((c) –near the interface z = − h 2 , (d) –near the interface
z = + h 2 ). Misfit strain um=0.05 % and temperature T=20 K (a, c, e) and 60 K (b, d, f). Extrapolation lengths are λP1=40 nm, λP3=40 nm, λS1=0.4 nm λS3=40 nm
5. Summary
Using Landau-Ginzburg-Devonshire (LGD) theory we study the influence of the flexoelectric and rotostriction coupling on the phase diagrams of ferroelastic-quantum paraelectric SrTiO3 films. We determined phase diagrams in coordinates temperature - SrTiO3 film thickness for different epitaxial misfit strains. Tensile misfit strains stimulate appearance of the spontaneous out-of-plane structural order parameter in the structural phase. Compressive misfit strains stimulate appearance of the spontaneous in-plane structural order parameter. Improper structural parameter components appear the vicinity of film surfaces as induced by the boundary conditions (surface energy). Gradients of the structural order parameter components, which inevitably exist in the vicinity of film surfaces, can induce improper polarization and pyroelectric response via the flexoelectric and 24
rotostriction coupling mechanism. Flexoelectric and rotostriction coupling results into the rotoflexoelectric field that is antisymmetric inside the film, small in the central part of the film, where the gradients of the structural parameter are small, and maximal near the surfaces, where the gradients of the structural parameter are the highest. The field induces antisymmetric improper polarization and pyroelectric response, which are zero in the central part of the film and maximal near the surfaces. Penetration depths of the improper phases (both polar and structural) can reach several nm from the film surfaces. Results obtained for SrTiO3 films evidence in favor of the rotoflexo coupling importance in such oxide perovskites as CaTiO3 and EuTiO3 films. The temperature dependence of improper and spontaneous polarization (via the temperature dependence of Ai and β iR (T , u m , h ) ) determines existence of pyroelectric response. Theoretically estimated Пi values ~ 10−5 - 10−4 C/m2K are comparable with those of known pyroactive materials PVDF and LiNbO3. The pyroelectric response originated from out-of-plane polarization component with symmetrical inversely polarized thickness distribution for thin films is negligible under thermal probing at frequencies f > v11 [ 1 ]. As a result, characteristic lengths
(
)
(
)
L1A T , u m = L2A T , u m =
− v 44
(
2α 1* T , u m
)
, which determine the gradient scale of the structural order
parameter components A1 and A2, are typically higher than lattice constant(s), but the length
(
)
L T , um = A 3
v11
(
2α *3 T , u m
12
)
that determines the gradient scale of the structural order parameter
component A3 is (much) smaller than the lattice constant, thus
(
)
* 2α *3 A3 + 4α *33 A33 + 2α 13 A12 + A22 A3 = v11
∂ 2 A3 ∂x32
≈0
(A.5c)
Approximate analytical solution of Eqs.(A.5a,b) in the structural phase (i.e. at α1* < 0 ) is: A1 ( x3 ) = Ab
⎛ ⎞ x3 2m + c0 m ⎟ , sn ⎜ ⎟ 1 + m ⎜⎝ L1A T , u m 1 + m ⎠
(
)
A2 ( x3 ) = 0 .
(A.6a)
* and constants m and Here sn(u|m) is the elliptical sine function. The amplitude Ab = − α1* 2α11
c0 should be determined from the symmetric boundary condition (A.4) at film surfaces x3 = ± h 2 .
3
Symmetry of the conditions gives c0 =
K (m ) , where K(m) is complete elliptic integrals of the first 2
kind [2]. Similar solution (1↔2) is A1 ( x3 ) = 0 ,
A2 ( x3 )Ab
⎛ x3 2m K (m ) ⎞⎟ + m . sn ⎜ ⎟ 1 + m ⎜⎝ L1A T , u m 1 + m 2 ⎠
(
(A.6b)
)
Solution, different from (A.6a-b) is A1 ( x3 ) = A2 ( x3 ) = Ac
⎛ x3 K (m ) ⎞⎟ 2m m . + sn ⎜ ⎟ 2 1 + m ⎜⎝ L1A T , u m 1 + m ⎠
(
(A.6c)
)
* * * Where the amplitude Ac = − α 1* (2α 11 ). + α 12 + α 13
Approximate analytical solution of Eqs.(A.5c) is:
(
)
⎧ α * A 2 ( x ) + A 2 ( x ) + α *3 * ⎪ − 13 1 3 * 2 3 A12 + A22 + α *3 < 0 , α 13 A3 ( x3 ) ≈ ⎨ 2α 33 ⎪ * 2 2 * ⎩0, α 13 A1 + A2 + α 3 > 0
(
(
)
(A.7)
)
Free energy (4) minimum determines which of the solutions (A.6) or (A.7) is more energetically preferable. Then one should substitute expressions for Ai into Eqs.(A.2d-g) for polarization. Relatively simple results for polarization can be derived open-circuited electrical boundary conditions, since P3 ( x3 ) appeared negligibly small due to the strong depolarization field. Moreover,
corresponding spatial scale, L (T , u m ) = P 3
12
g⊥
2β *3 (T , u m ) + 1 ε 0 ε b
, appeared ~0.1 nm, i.e. significantly
smaller than the lattice constant (~0.4 nm) of considered perovskites. Thus, putting P2 ( x3 ) ≈ 0 , we simplify Eqs.(A.2d-e) as:
(
∂ 2 P1
)
* * * * * 2β1* P1 + 4β11 P13 + 2β12 P1 P22 + 2 ξ11 A12 + ξ12 A22 + ξ13 A32 P1 + ξ 44 A2 A1 P2 = g ||
∂x
2 3
+
f 44 r44 ∂ ( A3 A1 ) c 44 ∂x3 (A.8a)
(
)
* * * * * 2β1* P2 + 4β11 P23 + 2β12 P2 P12 + 2 ξ11 A22 + ξ12 A12 + ξ13 A32 P2 + ξ 44 A2 A1 P1 = g ||
∂ 2 P2 ∂x
2 3
+
f 44 r44 ∂ ( A3 A2 ) c 44 ∂x3 (A.8b)
For the considered case either A1 ( x3 ) = 0 , or A2 (x3 ) = 0 or A1 ( x3 ) = A2 ( x3 ) ≠ 0 , and so Eqs.(A.8)
(
)
(
)
allow evident simplifications. Characteristic lengths L T , u m = L T , u m = P 1
P 2
g ||
(
2β1* T , u m
12
)
, which
determine the gradient scale of the polarization components P1 and P2, are typically higher than
4
lattice constant(s). The fact allows us to solve Eqs.(A.8) along with substitution either (A.6) or (A.7) and boundary conditions (A.3a,b) ⎛ S ⎞ ∂P f ⎜ 2b1 P1 ± g || 1 ± 44 r44 A1 A3 ⎟ ⎜ ⎟ ∂x3 c 44 ⎝ ⎠ ⎛ S ⎞ ∂P f ⎜ 2b1 P2 ± g || 2 ± 44 r44 A2 A3 ⎟ ⎜ ⎟ ∂x3 c 44 ⎝ ⎠
= 0,
(A.9a)
= 0,
(A.9b)
x3 = ± h 2
x3 = ± h 2
Analyses of equations (A.8) and conditions (A.9) along with expressions (A.6)-(A.7) lead to the conclusion that several polar phases are possible. 1) For the case A1 ( x3 ) = A2 (x3 ) ≡ A( x3 ) and P1 (x3 ) = P2 ( x3 ) ≡ P(x3 ) , polarization P( x3 ) obeys equation:
(2β
* 1
((
)
)) (
)
* * * * * + ξ 44 A 2 + 2 ξ11 + ξ12 A 2 + ξ13 A32 P + 4β11 + 2β12 P 3 − g ||
where A3 ( x3 ) ≈ −
* 2α 13 A 2 ( x3 ) + α *3
2α
* 33
∂ 2 P f 44 r44 ∂ ( A3 A) , (A.10) = c 44 ∂x3 ∂x32
* * for 2α13 A2 + α *3 < 0 , and A3 ( x3 ) ≈ 0 for 2α13 A2 + α *3 < 0 .
⎛ ⎞ ∂P f 44 ± r44 AA3 ⎟⎟ Boundary condition is ⎜⎜ 2b1S P ± g || ∂x3 c 44 ⎝ ⎠
= 0 . Solution of Eq.(A.10) can be found x3 = ± h 2
numerically, but approximate analytics is possible for determination of the ferroelectric and polar phase boundaries. Ferroelectric phase boundary can be determined as the instability threshold of the homogeneous solution. Analytical expression for the homogeneous solution of Eq.(A.10) was obtained in the VKB-type approximation ⎛ x3 γ ( x ) ⎞ ⎛ x3 γ ( x ) ⎞ ⎛ x3 γ ( x ) ⎞ Phom (x3 ) ≈ C1 exp⎜ i ∫ dx ⎟ + C 2 exp⎜ − i ∫ dx ⎟ ≡ C cos⎜ ∫ dx ⎟ (A.11) ⎜ 0 − g || ⎟ ⎜ 0 − g || ⎟ ⎜ 0 − g || ⎟ ⎠ ⎠ ⎝ ⎠ ⎝ ⎝ 2 g f Extrapolation length λ P = ||S , gradient coefficient g || = g 44 − 44 , the even function 2b1 c 44
(
((
)
))
* * * γ ( x ) = 2β1* + ξ 44 A2 + 2 ξ11 + ξ12 A2 + ξ13 A32 . Substitution of Eq.(A.12) into the homogeneous
⎛ ∂P ⎞ = 0 gives the instability threshold as the condition of zero boundary conditions ⎜⎜ P1 ± λ P 1 ⎟⎟ ∂ x 3 ⎠ x =±h 2 ⎝ 3 determinant h 2 ⎛ h 2 γ(x ) ⎞ γ (h 2) ⎛⎜ γ ( x ) ⎞⎟ ⎟ ⎜ (A.12) dx − λ P cos ∫ sin ∫ dx = 0 ⎜ 0 − g|| ⎟ ⎜ 0 − g|| ⎟ g − || ⎠ ⎠ ⎝ ⎝ Transcendental Eq.(A.12) determines the ferroelectric phase boundary in coordinates film thickness
h, misfit strain um and temperature T. The boundary depends on the extrapolation length λ P . 5
The boundary of the polar phase induced by the flexo-roto fields (inhomogeneities) is diffuse and can be estimated from the condition
∂ ( A3 A) ∂x3
≈0.
2) For the case A2 ( x3 ) = 0 and A1 ( x3 ) ≠ 0 (or A1 ( x3 ) = 0 and A2 ( x3 ) ≠ 0 ), P1 (x3 ) ≠ 0 and P2 ( x3 ) = 0 (or P1 ( x3 ) = 0 and P2 (x3 ) ≠ 0 ), we obtained that P1 ( x3 ) obeys equation: * * * 2(β1* + ξ11 A12 + ξ13 A32 )P1 + 4β11 P13 − g ||
where A3 ( x3 ) ≈ −
* α13 A12 ( x3 ) + α *3
2α
* 33
∂ 2 P1 ∂x
2 3
=
f 44 r44 ∂ ( A3 A1 ) , c 44 ∂x3
(A.13)
* * for 2α13 A12 + α *3 < 0 , and A3 ( x3 ) ≈ 0 for 2α13 A12 + α *3 < 0 .
⎛ ⎞ ∂P f Boundary condition is ⎜⎜ 2b1S P1 ± g || 1 ± 44 r44 A1 A3 ⎟⎟ ∂x3 c 44 ⎝ ⎠
= 0 . The ferroelectric phase boundary x3 = ± h 2
in coordinates film thickness h, misfit strain um and temperature T can be determined from
(
)
* * A12 ( x ) + ξ13 A32 ( x ) . Eq.(A.12), where γ ( x ) = 2 β1* + ξ11
For the case of short-circuited electric boundary conditions the critical thickness of the outof-plane ferroelectric polarization appearance can be determined from equation:
⎛ h 2 γ(x ) ⎞ ⎛ h 2 γ(x ) ⎞ γ (h 2 ) cosh⎜ ∫ dx ⎟ + λ P sinh ⎜ ∫ dx ⎟ = 0 , ⎜ ⎟ ⎜ ⎟ g g g ⊥ ⊥ ⊥ ⎝0 ⎠ ⎝0 ⎠ 1 1 and typically γ ( x ) ≈ . where γ ( x ) = 2β*3 + 2 ξ*33 A32 + ξ*31 A12 + A22 + ε 0εb ε 0εb
(
(
(A.14)
))
A.4. Material parameters used in calculations Table 1. SrTiO3 material parameters [3, 4]. Parameter εb αT T0( E )
SI units dimensionless 106×m/(F K) K
Value
43 0.75 30
Source and notes [ ] [c, d] ibidem
Tq( E )
K
54
ibidem
aij
109×m5/(C2F)
ibidem calculated from aiju
qij gijkl
1010× m/F 10-11×V⋅m3/C
a11u =2.025, a12u =1.215, a11σ =1.724, a12σ =1.396 q11=1.251, q12= −0.108, q44=0.243 g11=g44=1, g12=0.5
βT TS Tq( Φ )
1026×J/(m5 K) K K
9.1 105 145
[c] Estimation based on Ref. [e] [c] [c] [c]
bij
1050×J/m7
b11u =1.94, b12u =3.96, b11σ =1.69, b12σ =3.88
[c] calculated from biju
a, b
6
rij ηijkl
1030×J/(m5) 1029 (F m)-1
vijkl cij sij fijkl
1010×J/m3 1011×J/m3 10-12×m3/J V
r11=1.3, r12= −2.5, r44= −2.3 u u η11 = −3.366, η12 = 0.135, ηu44 =6.3 σ σ η11 = −2.095, η12 = −0.849, ησ44 =5.860 v11=0.28, v12= −7.34, v44=7.11 c11=3.36, c12=1.07, c44=1.27 s11=3.52, s12= −0.85, s44=7.87 f 11e = − 3.24 , f12e = 1.44 , f 44e = 1.08
[c] [c] calculated from ηuij [d, f] [c, d] calculated from cij Recalculated from Ref.[g] at given stress
a
G. Rupprecht and R.O. Bell, Phys. Rev. 135, A748 (1964).
b
G.A. Smolenskii, V.A. Bokov, V.A. Isupov, N.N Krainik, R.E. Pasynkov, A.I. Sokolov, Ferroelectrics and
Related Materials (Gordon and Breach, New York, 1984). P. 421 c
N. A. Pertsev, A. K. Tagantsev, and N. Setter, Phys. Rev. B 61, R825 (2000).
d
A.K. Tagantsev, E. Courtens and L. Arzel, Phys. Rev. B, 64, 224107 (2001).
e
J. Hlinka and P. Marton, Phys. Rev. B 74, 104104 (2006).
f
W. Cao and R. Barsch, Phys. Rev. B, 41, 4334 (1990)
g
P. Zubko, G. Catalan, A. Buckley, P.R. L. Welche, J. F. Scott. Phys. Rev. Lett. 99, 167601 (2007).
References 1
W. Cao and R. Barsch, Phys. Rev. B, 41, 4334 (1990).
2
I.S. Gradshteyn and I.M. Ryzhik Table of Integrals, Series, and Products, 5th ed edited by A
Jeffrey. Academic, New York (1994) 3
A.N. Morozovska, E.A. Eliseev, M.D. Glinchuk, Long-Qing Chen, Venkatraman Gopalan.
Phys.Rev.B. 85, 094107 (2012) 4
A.N. Morozovska, E.A. Eliseev, S.V. Kalinin, Long-Qing Chen and Venkatraman Gopalan.
Surface polar states and pyroelectricity in ferroelastics induced by flexo-roto field. Accepted to APL (http://arxiv.org/abs/1201.5085)
7
