Mesoscopic theory of defect ordering–disordering transitions in thin oxide films

Ordering of mobile defects in functional materials can give rise to fundamentally new phases possessing ferroic and multiferroic functionalities. Here we develop the Landau theory for strain induced ordering of defects (e.g. oxygen vacancies) in thin oxide films, considering both the ordering and wavelength of possible instabilities. Using derived analytical expressions for the energies of various defect-ordered states, we calculated and analyzed phase diagrams dependence on the film-substrate mismatch strain, concentration of defects, and Vegard coefficients. Obtained results open possibilities to create and control superstructures of ordered defects in thin oxide films by selecting the appropriate substrate and defect concentration.

(e.g. see Fig. 3 in Ref. 46 ). Following the theory by Khachaturyan 47 , we assume that concentration c is coordinate-independent.
The defect ordering direction can be different with respect to the film surfaces. Below we will distinguish perpendicular [see the left side (a) in Fig. 1] and parallel [see the central part (b) in Fig. 1] orientations of the planes of a constant order parameter and, thus, of equal sublattice occupations.
The η-dependent free energy of a thin film is the sum of the surface and bulk energies: where h is the film thickness, and the first term is the surface/interface energy that is not negative under the condition α S ≥ 0 . Hereinafter we will consider the case α S = 0 , corresponding to the so-called natural boundary conditions. The bulk density, f, of the Helmholtz free energy dependence on the order parameter η has the form: The first term in Eq. (1b), proportional to η 2 , originates from the interaction between defects from the different sublattices. The term can be readily derived from the term δn a δn b where δn a and δn b are the relative occupation numbers for the two sub-lattices. Hereinafter we will consider both cases α < 0 and α > 0 . For the former, defects may order spontaneously in a bulk material, while in the latter case the ordering in bulk is unfavorable, but can be induced by external factors (such as misfit strain).
The second term in Eq. (1b) is proportional to the entropy of the system with a minus sign, −TS 42,44 , since F = U − TS . The relative dimensionless entropy, defined as −(1 − η)ln(1 − η) − (1 + η)ln(1 + η) , tends to the minimum of −2ln2 under the condition η 2 → 1 , corresponding to the complete occupation of one of the sublattices ( δn b = 1 and δn a = 0 , or vice versa δn a = 1 and δn b = 0).We note that the entropy contribution stabilizes the system thermodynamically since its expansion contains all positive even powers of the order parameter.
The third and the fourth terms are the series expansion of generalized gradient energy in the derivatives ∂η/∂x i of the order parameter. The gradient energy was introduced by Cahn and Hillard in 1958 (see e.g. 48 ), since its inclusion becomes indispensable for correct description of inhomogeneous electrochemical systems with defect species. Below we consider the tensors g ij and w ijkl in isotropic approximation, g ij = gδ ij and w ijkl = wδ ij δ kl , and explore the case of w > 0 . This choice is because we study the case of the defects ordering for positive correlation energy.
Note that very often the higher gradient term is neglected ( w = 0 ). Rigorously it is justified only if g > 0 , assuming that its renormalization by the strain gradient energy is either positive or too small to change the positive sign of resulting effective gradient coefficient. Below we will show that w-term is mandatory to determine the threshold for emergence of the order parameter spatial modulation.
The elastic energy q η, σ ij includes the coupling between the order parameter and the stresses σ ij . Following the theoretical formalism of Vegard strains, proposed by Freedman 49 , and using the strain energy by Levanyuk et al. 50 , the coupling term between the order parameter and the stresses could be also expanded as a series on the order parameter: here s ijkl is the elastic compliances tensor, and u ij is the tensor of elastic strains. To describe the mixed mechanical state of the epitaxial binary oxide film on the substrate, we perform the Legendre transform by adding the term u kl σ kl . The third term in Eq. (1c) is the chemical expansion due to the appearance of elastic defects, ij characterize the Vegard strains for defects located within the sublattices (a) and (b), respectively. The symmetry of V (a,b) ij is determined by the local site symmetry of defect with respect to the lattice, which could be different from the lattice symmetry (see e.g. Ref. 49 considering different vacancies in SrTiO 3 ). Using the expressions of the partial sublattice occupation numbers through the order parameter η the third term in Eq. (1c) can be further transformed to c V c ij + ηV η ij with the mean value /2. It will be shown that the last term in Eq. (1c) plays a central role in the mechanism of defect ordering. It is a gradient-type striction due to the defect ordering that is characterized by a fourth rank tensor B ijkl .
In the continuous media approximation, the thermodynamically stable state of the film can be derived from the variation of the free energy (1) on η and σ ij , leading to the Euler-Lagrange differential equations: here we used the identity, arctanh(η) ≡ 1 2 ln 1+η 1−η . The boundary conditions to Eq. (2a) are the following: Scientific Reports | (2020) 10:22377 | https://doi.org/10.1038/s41598-020-79482-w www.nature.com/scientificreports/ The first term in Eq. (2c) originated from the variation of surface energy in Eq. (1a), and the so-called natural boundary conditions correspond to the case α S = 0 (surface energy is absent in this particular case). The second term in Eq. (2c) originated from the variation of the gradient energy in Eq. (1b), and the third term originated from the variation of the elastic energy Eq. (1c). Equation (2b) should be considered along with the conditions of mechanical equilibrium ∂σ ij /∂x j = 0 . The elastic boundary conditions to Eq. (2b) are the following: The first condition, σ i3 = 0 , means that the top surface of the film ( x 3 = 0 ) is mechanically free; and the second condition, u 11 = u 22 = u m , means that the bottom surface of the film ( x 3 = h ) is clamped to a rigid substrate, where u m is a misfit strain induced by the film-substrate lattices mismatch.
Below we consider thicknesses h smaller than the critical thickness of misfit dislocation appearance 43 , typically this means that h ≤ 10 nm.
The nonlinear boundary problem (2) contains a number of material-dependent constants, the majority of which are poorly known even for simple binary oxides, such as ZnO, MgO, SnO 2 , CeO 2 , HfO 2 , and even more so for complex ternary oxides, such as manganites (La,Sr)MnO 3 , paraelectric perovskites SrTiO 3 , EuTiO 3 , KTaO 3 , ferroelectric perovskites (Ba,Sr)TiO 3 , (Pb,Zr)TiO 3 , and orthoferrites PbFeO 3 , pristine and rare-earth doped BiFeO 3 , etc., which all can be deficient in oxygen. Hence, prior to solving the boundary problem (2) by e.g. finite element modeling (FEM), requiring all tabulated parameters, these should be taken from the experimental or density functional studies 5,6,8,19 . However, to facilitate the search in the multi-parameter space and open the way for further FEM, here we elaborate the analytical theory.

Analytical solution in harmonic approximation
To explore the phase evolution in the system described by the free energy, Eq. (1), we consider the three-dimensional Fourier series of the order parameter, where the wave vector components k lmn = 2π L l, 2π L m, 2π h n, are determined by the sizes L × L × h of the considered film, and m, n, and l are integer numbers. For the order parameter η(x) to be a real (i.e. observable) value, the equality η * l,m,n =η −l,−m,−n should be valid. Following Landau theory, below we assume that only one principal mode k dominates near the order-disorder phase transition, where δ is the phase shift and η 0 is the offset term. Next, we assume that the approximate equality is valid for the order parameter derivative in Eq. (1c), namely: Note that k is constant in the first harmonic approximation, but should be found in the self-consistent way in what follows.
Substitution of Eq. (3) in Eqs. (1b) and (1c) yields: where ∼ g (k) = g ij k i k j , ∼ Bij(k) = B ijmn k m k n and ∼ w (k) = w ijmn k i k j k m k n . Symmetry of the gradient-striction tensor B ijkl is determined by the local site symmetry of defect with respect to the lattice. Below we assume the tetragonal symmetry of the tensor ∼ Bij(k) , with a tetragonal axis along x 3 axis of the film.
In the continuum approximation, the thermodynamically stable state of the film can be analyzed by the variation of the free energy (1) that acquires the form: Scientific Reports | (2020) 10:22377 | https://doi.org/10.1038/s41598-020-79482-w www.nature.com/scientificreports/ As a next step, one can use a Galerkin procedure in the free energy (5a) with respect to the second harmonics assuming statistical averaging, equivalent to the spatial averaging in ergodic case. Such averaging "excludes" the contributions from the second and higher harmonics to cos[2(kx + δ)] 2n+1 due to a periodicity of trigonometric functions. Assuming the absence of correlations between the offset term η 0 and cos[2(kx + δ)] , all cross-terms proportional to η 2m+1 0 cos[2(kx + δ)] 2n+1 vanish after the averaging. Next, the terms proportional to η 2m+1 0 also vanish, while the even terms η 2m 0 do not contain any useful information about defect ordering. Thus, one can regard that �η 2 � ∼ ∼ η 2 ∼ η 2 , where the brackets ... designate the averaging, and obtain: (2), yields the algebraic equations of state: Near the phase transition, arctanh(η) ≈ η + η 3 3 in Eq. (6a). Assuming that the anisotropic Vegard tensor is diagonal (or at least can be diagonalized), that is true for many cases 49 Using this approximation, elastic solution for a thin oxide film on a rigid substrate is derived for cubic symmetry far from the structural domain walls. Zero components are σ 33 = σ 13 = σ 23 = σ 12 = 0 , u 12 = u 13 = u 23 = 0 , and nonzero components are: (7), and the Voight notations for B 1111 ≡ B 11 and B 1122 ≡ B 12 are used hereinafter.
Note, that in a freestanding film the stresses are zero σ ij = 0 , and elastic strains are Below we will suppose that η 0 = 0 and introduce the designations.
where the sum V m has the sense of partial molar volume, and the difference V n reflects the anisotropy impact.

Defect order-disorder transition
Long-range ordering parallel to the film surfaces. Here we consider the case when the only nonzero component of wave vector is k 3 , meaning that the harmonic modulation of the long-range order parameter η ∼ ∼ η cos[k 3 x 3 + δ] looks like planes parallel to the film surfaces x 3 = 0, h (see Fig. 3a, left). Since η is proportional to the difference of defect sub-lattice occupations δn a − δn b , this means the modulation of the sub-lattice occupation perpendicular to the film surfaces.
For the case the free energy density Eq. (1) has the following form (see Supplementary Materials):  4 3 . We refer the case with subscript "pr". Omitting the η-independent term α p , the free energy density (8a) can be expanded in power series in η as The derivation of Eq. (8) is given in Supplementary Materials 51 . Note that the renormalized coefficient before η 4 should be positive to ensure the stability of the phase described by the free energy (8b), otherwise higher terms should be included in the expansion. Since s 11 > |s 12 | for all elastically stable solids, the condition γ pr ≥ 0 is always valid.
The equilibrium values of the order parameter obtained from minimization of the energy (8b) and the corresponding free energy have the form: Long-range order exists if β pr c + k B T < 0.
Long-range ordering perpendicular to the film surface. We further consider the case when the only nonzero component of the wave vector is k 1 , , meaning that the harmonic modulation of the long-range order parameter η ∼ ∼ η cos[k 1 x 1 + δ] looks like planes perpendicular to the film surfaces x 3 = 0, h (see Fig. 3a, middle). Since η is proportional to the difference of defect sub-lattice occupations δn a − δn b , this means the modulation of the sub-lattice occupations parallel to the film surfaces.
The free energy density f in Eq. (1b) of the oxide has the following form in this case: We denote this case with subscript "pp". Omitting the η-independent term α p , the free energy density (9a) can be expanded in powers of η as Note that the renormalized coefficient before η 4 should be positive, otherwise higher terms should be included in the expansion. Thus, the inequality γ pp c 4 + ck B T > 0 should be verified. Since s 11 > |s 12 | for all elastically stable solids, and B22 the condition γ pp > 0 is always valid, and so γ pp c 4 + ck B T > 0 at finite temperatures.
The equilibrium values of the order parameter obtained from minimization of the energy (9b) and corresponding energy have the form similar to Eq. (8c): Long-range order exists if β pp c + k B T < 0 . The derivation of Eq. (9) is given in Supplementary Materials 51 . Note that, in the ordered phases, elastic dipoles may become polar due to the surface-induced piezoelectric coupling that, in turn, originates from the inversion symmetry breaking near the film surfaces. At the same time, the out-of-plane polar phase is affected by the strong depolarization field originated from the sharp gradient of polarization decay away from the surfaces.
Structural phase diagrams. The temperatures of transitions between defect-disordered and defectordered phases for parallel and perpendicular orientations of the planes of a constant order parameter can be determined from equations |g|(s11+s12) are dimensionless gradient-related parameters. Below we will analyze the positive roots k 1,3 > 0 only, since the negative ones describe the same physical states. Here we also introduced anisotropy term g a = (B 11 − B 12 ) cV n 2(s 11 −s 12) . Figure 2 illustrates the square-root like dependences of the dimensionless wave-vector components, k 1 /k 0 and k 3 /k 0 , on the effective strain u for g > 0 (red and magenta curves) and g < 0 (black and blue curves). Note that positive w determines the existence of the modulation and its characteristic wave number.
The transition temperatures, corresponding to the η-modulation vectors given by Eq. (11), are The conditions of the equilibrium between the ordered phases can be obtained from the equality of the corresponding free energy densities, Eqs. (8c) and (9c), namely.
allowing for the conditions of the long-range orders existence www.nature.com/scientificreports/ Note that the expressions (11) for k-vector components should be substituted in Eqs. (13). Equation (13a) can be solved with respect to temperature in the following form: The signs " ± " in Eq. (14a) correspond to two different roots, which have physical sense if correspond to T ≥ 0. These roots should further satisfy two conditions (13b). Depending on the parameters it appeared possible at least for the largest root of Eq. (14a). Since the sign of denominator in Eq. (14a) is not fixed, both signs in the numerator, +Det β pp , β pr or −Det β pp , β pr , are possible depending on the parameters. Hence, the number  It can be shown that, at the chosen parameter values, one of the roots (14a) is close to the expression T pp (u m , c) + T pr (u m , c) /2 . Hence, for β pp ≈ β pr expression (14a) can be expanded in series on the small difference β pp − β pr powers, namely When a strong inequality γ pp − γ pr << γ pp takes place, the second term in parenthesis in Eq. (14c) could be either close to unity (sign " − ") or much higher than unity (sign " + "). It is seen that in the case of sign " − " the solutions (14b)-(14c) will not satisfy both of relations (13b) simultaneously. The other root in Eq. (14a) deviates more significantly from the value T pp (u m , c) + T pr (u m , c) /2 , thus, in this case, conditions (13b) will be satisfied.
Using Eqs. (10)- (14) we plotted phase diagrams of the system and corresponding modulation amplitudes and wave vectors, shown in Figs. 3, 4 and 5, respectively. Note that the ordered phases OP ⊥ and OP || have different orientation of the order parameter modulation amplitude, η(x 1 ) and η(x 3 ) , with respect to the film surfaces. Figure 3a represents the scheme of the η(x 1 , x 3 ) distribution in the ordered (OP ⊥ and OP || ) and disordered (DP) phases. The ordered phases OP ⊥ and OP || are characterized by periodic changes of η = c c 0 (δn b − δn a ) , which are proportional to the difference of the sublattices "a" and "b" occupation numbers, δn b − δn a . The total concentration "c" remains constant, since the sum δn a + δn b is independent on η . The OP || phase is characterized by a periodic change η ∼ ∼ η cos[k 3 x 3 + δ] , corresponding to the long-range modulation of the defect sub-lattices "a" and "b" occupation perpendicular to the substrate plane x 3 = h . The OP ⊥ phase is characterized by a periodic change η ∼ ∼ η cos[k 1 x 1 + δ] , corresponding to the long-range modulation of the defect sub-lattices "a" and "b" occupation parallel to the substrate plane x 3 = h . The disordered DP phase is characterized by η = 0 , corresponding to the equal filling of both sub-lattices. Note that the modulation periods, 2π k 1 and 2π k 3 , should be significantly larger than the lattice constant for the validity of the continuum approach, and so Fig. 3a illustrates only these long-range modulations of the order parameter η(x) , but not the distance between the defect sub-lattices planes. Figures 3 and 4 are phase diagrams of an ultra-thin film in dependence on the misfit strain u m and the dimensionless defect concentration c/c 0 calculated for negative and positive α, respectively. Corresponding modulation wave vectors, k ⊥ and k || , are shown in Fig. 5a,b, respectively. Figure 3b-e, calculated for negative α, demonstrate the presence of OP ⊥ , OP || and DP phases, whose boundaries depend on the gradient g and striction B ij coefficients, which are different for Fig. 3b-e. The ultrathin almost vertical dark green stripe between OP || and OP ⊥ phases located at small |u m | < 0.1% in Fig. 3b,c is the coexistence www.nature.com/scientificreports/ region of the ordered phases, where the wave vector becomes very small, k → 0 . This phase is absent for the diagrams in Fig. 3d,e, which have a tricritical point, where OP ⊥ , OP || and DP phases coexist. The tricritical point {u m = 0, c/c 0 = 0.63} is almost independent on g and B ij signs. The OP || -DP and OP ⊥ -DP boundaries are curved for all diagrams. DP phase occupies mountain-shape region with the "top" at {u m = 0, c/c 0 = 0.63}, which hillsides start at small misfits. The region of its stability gradually decreases with |u m | increase. At that the noticeable asymmetry between compressive ( u m < 0 ) and tensile ( u m > 0 ) misfit strains is evident and related with positive Vegard effect ( V m > 0 ). It can be deduced from Fig. 3b-e, that the change of the B ij sign leads to an interchange between the OP || and OP ⊥ phases with different orientation of η-ordering, stemming from the chemo-strictive coupling between the strain u m and defect concentration c, expressed by the coupling terms B 11 +B 12 s 11 +s 12 (u m − cV m /2) and 2B 12 s 11 +s 12 (u m − cV m /2). The OP || -OP ⊥ boundary is almost vertical, and its small slope weakly depends on the V m value. The OP ⊥ phase exists for compressive misfit strain u m < 0 at B 11 + B 12 < 0 , and for tensile misfit strain u m > 0 at B 11 + B 12 > 0 . The OP || phase exists for compressive misfit strain u m < 0 at B 11 + B 12 > 0 , and for tensile misfit strain u m > 0 at B 11 + B 12 < 0 . The color maps of the order parameter amplitudes [ η(x 1 ) and η(x 3 ) ] and wave vectors (k ⊥ and k || ) calculated for α < 0 , g > 0 , B 11 + B 12 > 0 are shown in Fig. 5a,c, respectively. The continuous transition from OP || to OP ⊥ phase occurs at u m = 0 . Notably that k ⊥ = k || = 0 in the region u m ≈ 0 . The values of η and k gradually increase with |u m | increase.
The change of α sign critically affects the phase diagrams, as is shown in Fig. 4. The most pronounced effect is the appearance of the wide DP region between the ordered OP ⊥ and OP || phases, so that the tricritical point is absent for all concentrations "c". Note that the DP exists in the regions of misfits |u m | ≤ 2% and this region grows with defect concentration decrease. The width of the DP region slightly increases with |u m | increase. The change of the B ij sign leads to the interchange between the OP || and OP ⊥ phases, while the effect of the coefficient g is negligibly small in this case [compare Fig. 4a,b with 4c-d]. The OP ⊥ phase exists for high compressive misfit strain u m < −2.5% at B 11 + B 12 < 0 , and for tensile misfit strain u m > +2.5% at B 11 + B 12 > 0 . The OP || phase www.nature.com/scientificreports/ exists for compressive misfit strain u m < −1.9% at B 11 + B 12 > 0 , and for tensile misfit strain u m > +1.9% at B 11 + B 12 < 0 . The color maps of the order parameter amplitudes [ η(x 1 ) and η(x 3 ) ] and wave vectors (k ⊥ and k || ) calculated for α > 0 , g > 0 , B 11 + B 12 > 0 are shown in Fig. 5b,d, respectively. The OP || to OP ⊥ phases are separated by a wide region of the DP phase located at |u m | ≤ 2% . The values of k are nonzero at the DP-OP || and DP-OP ⊥ boundaries. The values of η and k gradually increase with |u m | increase. Note that the values of the order parameter η and the components of the wave vector, k 1 and k 3 , were calculated on the basis of numerical minimization of the more complex free energy functional [see Eqs. (8a) and (9a)], for which the series expansion was not used. It turned out that the differences between the approximate expressions (11) and the results of numerical calculations are insignificant for the wave vector. Also, to describe the phase transitions between the ordered and disordered phases, expressions (12) exactly correspond to numerical calculations, since these transitions are the second-order phase transitions and, the order parameter can be considered small near the transition points.
Note that phenomenological parameters listed in the caption to Fig. 3 were selected in such a way so as to satisfy the physical conditions αc 0 ∼ = k B T room , |α|/g ≥ k 0 , g /w < k 0 , while s ij and c 0 values are typical for oxides. As for the striction coefficients B ij , these are chosen so that the combinations of parameters, B 11 +B 12 s 11 +s 12 (u m − cV m /2) and 2B 12 s 11 +s 12 (u m − cV m /2) , are of the same order as g . To summarize the section, results shown in Fig. 3, 4 and 5 indicate that one can control the defect ordering-disordering by changing their concentration c/c 0 and misfit strain u m at fixed values of the other parameters.
Finally it may be interesting to compare the defect ordering with the cycloid ordering in a strained multiferroic BiFeO 3 thin films (see e.g. Ref. 11,21,22 and refs. to original papers therein). Actually, the strain-induced singularity of the cycloid period and phase diagram with various spin cycloid orientations in BiFeO 3 can be topologically analogous to the one considered in Figs. 3, 4 and 5 proving that chiral incommensurate phases of defect ordering (allowed by the theoretical formalism) can exist in the strained film.

Conclusion
We have analyzed the ordering of defects (e.g. oxygen vacancies) in thin oxide films in the framework of the continuum Landau-type theory. We derived analytical expressions for the energies of various defect-ordered states and calculated and analyzed phase diagrams dependence on the film-substrate misfit strain and concentration of defects for different gradient, striction and Vegard coefficients.
We have found that two defect-ordered phases, which are characterized by either parallel or perpendicular defect ordering in planes and corresponding wave vectors, can be stable. The stability conditions are determined by the misfit strain and the defect concentration at fixed values of the other parameters. The ordered phases border with the defect-disordered phase. Hence, we have shown that it is possible to control the defect orderingdisordering by changing their concentration and the film-substrate misfit strain (compressive or tensile). Thus, the obtained results open possibilities to create and control superstructures of ordered defects in thin oxide films by selecting the appropriate substrate and defect concentration.