Formation of polarization needle-like domain and its unusual switching in compositionally graded ferroelectric thin films: an improved phase field model

Abstract

Compositionally graded ferroelectrics (cgFEs), which possess a spatial variation in material composition, exhibit anomalous or even unprecedented properties. Despite several breakthroughs having been achieved in experiments, there surprisingly is a lack of an effective simulation approach for cgFEs, thereby greatly hindering a deep understanding about underlying mechanisms hidden behind the observed phenomena. In this study, an improved phase field model is proposed for a cgFE made of PbZr_(1−x)Ti_xO₃ based on the Ginzburg–Landau theory. The improved approach enables us to capture several key phenomena occurring in the cgFE PbZr_0.8Ti_0.2O₃ ⇔ PbZr_0.2Ti_0.8O₃ thin film that are observed from experiments, such as the formation of needle-like domains with curved domain walls, ferroelastic switching under an electric field, and the voltage offset of the hysteresis loop. Results obtained from the improved approach indicate that the high elastic energy near the needle-tip of an a-domain gives rise to a deflection of the domain wall from the regular a/c domain wall, while the high concentration of the depolarization field shrinks a part of the a-domain near the needle-tip and terminates the a-domain within the cgFE thin film. These facilitate the stabilization of needle-like domains with curved domain walls in the cgFE thin film. Furthermore, the flexoelectricity is proven to play an important role in the voltage offset of the hysteresis loop. On the other hand, the ferroelectric switching process in the cgFE thin film exhibits a vastly different response in comparison to that in homogeneous ferroelectric thin films, including local switching initiation and formation and annihilation of vortex–antivortex pairs during the switching. The present study, therefore, provides an incisive approach for investigations on cgFEs, which may bring new understanding and unique insights into these complex materials, as well as novel potential applications.

---

1 Introduction

Spontaneous breaking of spatial-inversion symmetry at the unit cell level results in the formation of a stable yet switchable electrical polarization in ferroelectrics, which, in turn, allows the materials to be applied across a broad range of advanced technologies, including nonvolatile random access memory (FeRAM) devices;¹ sensors, actuators, and transducers;² field-effect transistors;^3,4 tunnel junctions;⁵ magnetoelectric-coupling devices;^6,7 and switchable photovoltaics.⁸ Recently, an additional macroscopic broken spatial-inversion symmetry has been introduced into the system by using compositionally graded ferroelectric (cgFE) thin films, in which the chemical compositions are spatially varied throughout the thickness of the thin films. This macroscopic spatial-inversion symmetry breaking brings about novel physical properties, such as the presence of built-in electric fields,⁹ shifted hysteresis loops,^9–12 large susceptibilities,¹³ and a potential for geometric frustration.¹⁴ Therefore, cgFEs have emerged as a potential candidate for novel applications in electric devices,¹⁵ and attracted considerable experimental and theoretical attention in the last decade.

Experimentally, it was reported early on that an offset in the vertical (polarization) axis of the hysteresis loop occurs in cgFE heterostructures, and the polarization gradients in the material were primarily believed to be the origin of the observed offsets.^11,16 Later, an experimental investigation on the origin of hysteresis loop shifts in cgFEs⁹ indicated that the vertical offsets were explicitly dependent on the measurement circuit^16,17 and the offsets should be along the horizontal (voltage) axis because of a built-in electric field. The origin of horizontal offsets results from the macroscopic broken spatial-inversion symmetry, and thereby, is intrinsic to the cgFE materials. Recently, a large horizontal offset of the ferroelectric hysteresis loop has been further demonstrated to originate from strong internal long-range electric fields that arise from compositional variations, by using high quality cgFE thin films.¹⁸ On the other hand, a high dielectric response over a wide temperature range has been reported for BaTiO₃ cgFE thin films with altered Ba-site stoichiometry along the thickness.^19–21 In these studies, the interesting characteristic overlays the broad and flat profile of the dielectric constant versus temperature, i.e., there is an absence of anomaly dielectric that is commonly observed in homogeneous ferroelectrics. The temperature independent dielectric constant originates from a strong smearing of the paraelectric–ferroelectric phase transition of the entire heterostructure accompanied by a reduction of the global transition temperature of the entire layer, with the highest Curie temperature due to the internal (depolarizing) electric fields. More recently, unexpected crystal and domain structures, and thereby novel properties, have been observed in cgFE thin films.^22–24 These studies have demonstrated that large polarization gradients and desirable temperature-stable susceptibilities can be achieved through the control of composition and strain gradients in ferroelectrics.²³ More interestingly, strain gradients in compositionally graded heterostructures made of PbZr_(1−x)Ti_xO₃ can stabilize unusual needle-like ferroelastic domains that terminate inside the film.²⁴ The stabilization of such a needle-like domain structure in a cgFE originates from complex interactions among the electrical and mechanical fields, and rendering intuitive guesses about these systems becomes quite difficult. Despite a significant scientific interest and an immense potential for applications, the structure–property relationships in cgFE thin films remain unclear since the physical mechanisms behind the observed behaviors, which arise from a plethora of intrinsic and extrinsic factors, have been difficult to identify. In turn, this has limited the adoption and utilization of cgFEs in advanced ferroelectric-based devices.

Theoretically, an attempt to understand the exotic properties in cgFEs has been preliminarily made with the use of transverse Ising models.^25,26 Recent development in theoretical approach based on principles of non-linear thermodynamics can shed light on the polarization stabilities, phase transformation characteristics and dielectric properties of cgFEs.^27–29 These theoretical studies provide an explanation as to why cgFEs should have lower dielectric losses and lower leakage currents compared to homogeneous ferroelectrics. In addition, role of the flexoelectric effect in the enhancement and stabilization of polarization components in cgFEs has also been demonstrated by using an extended thermodynamic theory.^27,28,30 Although the use of theoretical approaches can explain some experimental observations, there still lacks a consideration of domain structures and their response to the mechanical and electric excitations, despite of their crucial roles in determination of macroscopic properties. This thus limits the use of theoretical methods for cgFEs, particularly with complex domain structures. On the other hand, the domain structure in cgFEs can be investigated using first-principles calculations.¹⁴ However, the heavy computational load imposes restrictions on the number of atoms and temperature in the simulation system that makes it critically difficult for further investigation on the behavior of polarization and domain structures in cgFEs under mechanical and electric fields. To circumvent these difficulties and to provide an effective approach for investigations on cgFEs, a distinct strategy is essential. As a potential alternative, phase field modeling based on the Ginzburg–Landau theory has been commonly used to study domain structures in ferroelectric materials, and polarization switching under a static electric field and/or mechanical loading.^31–37 Although the phase field model is, until now, only developed for homogeneous ferroelectric systems, an integration of a non-linear thermodynamic formalism based on Ginzburg–Landau theory for the cgFE system into phase field frameworks is quite promising.

In this work, a new phase field model is developed for cgFE thin films made of PbZr_(1−x)Ti_xO₃, which incorporates an extended non-linear thermodynamic formalism based on the Ginzburg–Landau theory. The developed phase field model is validated through an investigation on the formation of needle-like a-domains that have recently been experimentally observed in a cgFE PbZr_0.8Ti_0.2O₃ (top) ⇔ PbZr_0.2Ti_0.8O₃ (bottom) thin film.²⁴ Based on the simulation results, the mechanism hidden behind the formation of needle-like a-domains with curved domain walls is unveiled. In addition, the ferroelastic and ferroelectric switching of the needle-like domain in the cgFE thin film are also investigated. Furthermore, the underlying mechanism for the voltage offset of the hysteresis loop that was observed in experiment is explored and discussed.

2 Computational details

2.1 Phase field model of cgFEs based on the Ginzburg–Landau theory

Here, we carry out the groundwork for an extension of the thermodynamic formalism for cgFEs based on the Ginzburg–Landau theory. To study PbZr_(1−x)Ti_xO₃ (PZT), we have utilized a formalism expanded to the sixth order. The material parameters of cgFEs are functions of the mole fraction x of PbTiO₃ in PZT, e.g., a^(x) = f(x), where a^x is any material parameter, such as thermodynamic stiffness coefficients, elastic compliance or stiffness, electrostrictive coefficients, and Curie temperature of the constitutive materials. x = 0 corresponds to the pure PbZrO₃ material and x = 1 corresponds to the pure PbTiO₃ ferroelectric. In the present study, the compositional materials are assumed to linearly vary along the thickness of the thin film; that is the x₃ direction. Therefore, material parameters are also functions of the coordinate along thickness, a^(x) = f(x) = g(x₃). The schematic of the finite element model for the cgFE structure and the change of material properties along the x₃ direction are provided in Fig. 1. A set of material parameters corresponding to a particular mole fraction x are assigned to an element layer depending on their coordinates in the x₃ direction. The model, thus, strongly depends on the mesh sizes. To provide a good presentation for the cgFE model, small uniform meshing along the material gradation direction are required. In this study, uniform elements with the sizes of 0.4 nm in the thickness direction, which is as small as the PbTiO₃ lattice constant,³⁸ are employed.

Fig. 1 Sketch outlining the gradual variation of material properties along the x₃ coordinate as captured by the cgFE model.

The total free energy of the cgFE system, F_FE, is essentially the sum of the free energies of each element layer, and is described as:

where V_n is the volume of the n^th element layer, f_FE^(x) is the total free energy density of ferroelectric with the mole fraction x, and m is the number of element layers to describe the gradation of material constituents. f_FE^(x) is described by:^31–34,39

f_FE^(x) = f_L^(x)(p_i) + f_G(p_i,j) + f_Elas^(x)(ε_ij, p_i) + f_Flex(p_i, ε_ij) + f_E^(x)(p_i, E_i), (i, j = 1, 2, 3) (2)

where the superscript (x) indicates that the functions are dependent of the mole fraction x and spatially-dependent; p_i is a polarization vector; p_i,j = ∂p_i/∂x_j denotes the derivative of the ith component of the polarization vector p_i, with respect to the jth coordinate x_j; ε_ij and E_i denote the components of strain and electric field, respectively. f_L^(x)(p_i), f_G(p_i,j), f_Elas^(x)(ε_ij, p_i), f_Flex(p_i, ε_ij), and f_E^(x)(p_i, E_i) are the Landau energy density, the gradient energy density, the elastic energy density, the flexoelectric energy density, and the electrostatic energy density, respectively.

The Landau energy density is expanded in terms of polarization components up to the sixth order as:^40–43

where α₁^(x) is the dielectric stiffness, and α₁₁^(x), α₁₂^(x), α₁₁₁^(x), α₁₁₂^(x), and α₁₂₃^(x) are higher order dielectric stiffness constants. The corresponding Landau coefficients are given by:^43–45

α₁^(x) = (T − T₀)/2κ₀C₀,

α₁₁^(x) = (10.612 − 22.655x + 10.955x²) × 10¹³/C₀,

α₁₂^(x) = η₁/3 − α₁₁^x,

α₁₁₁^(x) = (12.026 − 17.296x + 9.179x²) × 10¹³/C₀,

α₁₁₂^(x) = (4.2904 − 3.3754x + 58.804e^−29.397x) × 10¹⁴/C₀,

α₁₂₃^(x) = η₂ − 3α₁₁₁^(x) − 6α₁₁₂^(x)

where T and T₀ denote the temperature and the Curie–Weiss temperature, respectively, C₀ denotes the Curie constant and κ₀ is the dielectric constant of a vacuum, which are given as:^43–45

κ₀ = 8.85 × 10⁻¹²,

η₁ = [2.6213 + 0.42743x − (9.6 + 0.012501x) × e^−12.6x] × 10¹⁴/C₀,

η₂ = [0.887 − 0.76973x + (16.255 − 0.088651x) × e^−21.255x] × 10¹⁵/C₀,

T₀ = 462.63 + 843.4x − 2105.5x² + 4041.8x³ − 3828.3x⁴ + 1337.8x⁵,

when 0.0 ≤ x ≤ 0.5,

when 0.5 ≤ x ≤ 1.0.

The gradient energy density for ferroelectric materials represents the domain wall energy density that gives the energy penalty for the spatially inhomogeneous polarization. The gradient energy density is described by:⁴⁶

f_G(p_i,j) = G_ijklp_i,jp_k,l, (k, l = 1, 2, 3) (4)

where G_ijkl are the gradient energy coefficients. In the present study, the gradient coefficients are assumed to be independent of the mole fraction x, and equal to the values for PbTiO₃.^35,39 The elastic energy density induced by mechanical strain is given by:^31,43where c_ijkl^(x) are the elastic coefficients and ε⁰_ij = Q_ijkl^(x)p_kp_l is spontaneous strain or eigenstrain and Q_ijkl^(x) is the electrostrictive coefficient tensor. The elastic constants for a number of compositions are adopted from a previous study,⁴⁷ and are given in Table 1. Note that, as the elastic coefficients slightly change with the change of mole fraction x, the elastic coefficients of material with 0.2 ≤ x ≤ 0.4 adopt the values for x = 0.4. For cubic symmetry, electrostrictive coefficients are given as:^43,45

Table 1 Elastic coefficients of PbZr_(1−x)Ti_xO₃ for a number of compositions

Ti content x	0.4	0.5	0.6	0.7	0.8	0.9	1.0
c11 (10^11 N m^−2)	1.681	1.545	1.696	1.712	1.728	1.704	1.746
c12 (10^11 N m^−2)	0.826	0.840	0.819	0.811	0.802	0.760	0.7937
c44 (10^11 N m^−2)	0.406	0.348	0.472	0.571	0.694	0.833	1.111

---

Flexoelectricity, the electromechanical coupling between the polarization and strain gradient (direct effect) or mechanical stress and electric field gradient (converse effect), is a property of all insulators.⁴⁸ The flexoelectric coupling energy density is expressed as:^48–51

where ε_ij,l = ∂ε_ij/∂x_l and f_ijkl are components of the flexocoupling tensor. The electrostatic energy density, which is obtained through Legendre transformation, is expressed as:^31,52where κ_c is the dielectric constant of the background material.

The temporal evolution of the polarization field towards its thermodynamic equilibrium state is described by the time-dependent Ginzburg–Landau (TDGL) equation:

where L is a kinetic coefficient, r = (x₁, x₂, x₃) is the spatial vector, and t denotes evolution time.

On the basis of Mindlin’s strain gradient theory⁵³ the mechanical equilibrium equation accounting for contributions of strain gradients is described as:

where b_k are the components of a body force per unit volume. The equilibrium of the electric field is governed by Maxwell’s (or Gauss’) equation as:where ζ is the volume charge density.

To solve the equilibrium equations, a numerical algorithm based on the finite element method (FEM) is employed in the present study. Using the variation or principle of virtual work, the weak form of the governing eqn (8)–(10) is expressed as:

where u_i is displacement; t_i and ω are the surface traction and the surface charge, respectively; ϕ is electrical potential; (∂f_FE^(x)/∂p_i,j)ξ_j denotes the surface gradient flux; and ζ_j is the component of the normal unit vector of the surface. For the space discretization we use hexahedral elements, of which each node is assigned with seven degrees of freedom, including three polarization components, one electrical potential variable, and three displacement components. In comparison with the phase field model that uses the Fourier spectral iterative perturbation algorithm,^32,54 the phase field model based on the finite element algorithm does not require uniform hexahedral (3D) and quadrilateral (2D) meshes and periodic boundary conditions for numerical implementation.⁵⁵ The present phase field model is able to simulate ferroelectrics with arbitrary boundary conditions and geometrical shapes. Therefore, the phase field model with an FEM-based algorithm is suitable and effective for the study of domain evolution in cgFEs.

In general, the finite element stiffness equations can be written as:

K^e × d^e = F^e (12)

where, d^e = (u^e, ϕ^e, p^e) are the degrees of freedom of an element in the cgFE system, K^e is the tangent matrix, and F^e denotes the load vector. The tangent matrix K^e has a block structure and can be described as:where,





here, Ω_e is the volume of an element. In the present work, the matrices [c(x₃)], [Q(x₃)], [Q′(x₃)], [κ_c(x₃)], and [α(x₃)] are functions of spatial coordinates along the x₃ direction, which are distinguished from the previous phase field model for homogeneous ferroelectrics where all the functions are independent of the coordinates. Details of these matrices are provided below, while details of the remaining matrices are provided in previous studies.^56,57 Therefore, tangent matrix K^e can be described for the spatial variation of material in cgFE systems.where,

q₁₁^(x) = c₁₁^(x)Q₁₁^(x) + 2c₁₂^(x)Q₁₂^(x),

q₁₂^(x) = c₁₁^(x)Q₁₂^(x) + c₁₂^(x)(Q₁₁^(x) + Q₁₂^(x)),

q₄₄^(x) = c₄₄^(x)Q₄₄^(x).

where,

AA₁₁^(x) = 2α₁^(x) + 4α₁₁^(x)p₁₂ + 2α₁₂^(x)(p₂₂ + p₃₂) + 6α₁₁₁^(x)p₁₄ + α₁₁₂^(x)[4p₁₂(p₂₂ + p₃₂) + 2(p₂₄ + p₃₄)] + 2α₁₂₃^(x)p₂₂p₃₂,

AA₂₂^(x) = 2α₁^(x) + 4α₁₁^(x)p₂₂ + 2α₁₂^(x)(p₁₂ + p₃₂) + 6α₁₁₁^(x)p₂₄ + α₁₁₂^(x)[4p₂₂(p₁₂ + p₃₂) + 2(p₁₄ + p₃₄)] + 2α₁₂₃^(x)p₁₂p₃₂,

AA₃₃^(x) = 2α₁^(x) + 4α₁₁^(x)p₃₂ + 2α₁₂^(x)(p₁₂ + p₂₂) + 6α₁₁₁^(x)p₃₄ + α₁₁₂^(x)[4p₃₂(p₁₂ + p₂₂) + 2(p₁₄ + p₂₄)] + 2α₁₂₃^(x)p₁₂p₂₂.

Note that in the present study the thermodynamic formalism and the material coefficients of PbZr_(1−x)Ti_xO₃ measured by M.J. Haun et al.^{41,44,45,58,59} are adopted into the phase field model for the cgFE thin film. Such thermodynamic formalism and material coefficients of PbZr_(1−x)Ti_xO₃ have been widely used in phase field simulations and have yielded highly consistent results to experiment data.^60,61 The thermodynamic formalism and the material coefficients of PbZr_(1−x)Ti_xO₃ measured by M. J. Haun et al. are reliable and can be used to provide physical meanings for experimental observations. Therefore, it is reasonable to further employ the thermodynamic formalism and the material coefficients of PbZr_(1−x)Ti_xO₃ in the phase field model of cgFEs, in order to reproduce the experimental data with support of the proper physics.

2.2 Simulation model and procedure

According to the previous experimental study,²⁴ here we consider 100 nm-thick, linear compositionally graded PbZr_0.8Ti_0.2O₃ (top) ⇔ PbZr_0.2Ti_0.8O₃ (bottom) heterostructures. The simulation model is constructed with dimensions of 100 × 10 × 100 nm. We introduce a Cartesian coordinate system with the x₁x₂ plane coincident with the thin film plane and the x₃ direction along the thin film thickness. The mechanical boundary condition on the bottom surface is: the node at the center is fixed and other nodes can only move in the x₁ − x₂ plane, such that the substrate strain can be relaxed by the applied field. Periodic boundary conditions are applied along the x₁ and x₂ axes. The spatial mesh of the simulation model is composed of 2.5 × 10⁵ hexahedral elements. The thin film consists of a linear cgFE with the [100], [010], and [001] crystallographic orientations along the x₁, x₂, and x₃ directions, respectively.

The thermodynamic equilibrium of the polarization field in the cgFE thin film without external fields is obtained by starting from an initial paraelectric state with small random perturbations. The evolution of the polarization field is numerically simulated by iteratively solving the TDGL equation. A backward Euler scheme and Newton iteration method are used in the FEM algorithm for the time integration and nonlinear iteration, respectively. The polarization configuration is stably formed when an insignificant change in polarization is observed.

The electrical potential at the top and bottom surfaces is firstly set to zero to allow the film to reach equilibrium. After equilibrium is reached, electrical switching is simulated by altering the electrical potential at the top surface of the thin film with a time-dependent sinusoidal electric field as E = E_max sin(2πft), where E_max and f are the amplitude and frequency of the applied field, respectively. The amplitude is set to be sufficiently high to ensure that saturation of the polarization switching could be reached. At each applied electric field, the simulated ferroelectric film is allowed to evolve one step with a simulation time step of Δt* = 0.01. For computational convenience, the used material parameters are normalized, as given in the previous study.³⁵

3 Results and discussion

3.1 Formation of a polarization domain structure in the cgFE thin film

Typical snapshots of temporal polarization distributions in the cgFE thin film at different normalized evolution time are illustrated in Fig. 2. At the beginning of the polarization evolution, a random distribution of polarization is assigned as an initial state (Fig. 2a). Then, polarizations in some local regions gradually align each other to form local rectilinear domains (Fig. 2b), which energetically favour four crystallographic orientations in the x₁x₃ plane. The formed domains mutually compete, in which some domains expand their sizes while the others shrink and then disappear, i.e., coarsening domain process (Fig. 2c). Specifically, at t* = 6.0, the c-domains cover a large portion of the domain structure in the cgFE thin film, while the a-domains become small and mostly cluster near the bottom surface (Fig. 2c). The c-domain continually expands during the evolution, while the a-domains merge into a single domain, then shrink (Fig. 2d). At the end of evolution, the cgFE thin film exhibits a c/a/c/a-like domain structure with an exclusively needle-like a-domain that terminates in the middle of the film. The system then achieves the equilibrium state. The averaged width of the needle-like domain at the stable state is as small as 10 nm. Compared with the recent experimental report,²⁴ the needle-like domain structure obtained from our developed phase field model achieves highly consistent characteristics, which clearly demonstrates the reliability of the proposed approach.

Fig. 2 Formation of the domain structure in the cgFE PbZr_0.8Ti_0.2O₃ (top) ⇔ PbZr_0.2Ti_0.8O₃ (bottom) thin film. Snapshots of temporal polarization distributions at different evolution times: (a) t* = 0.0, (b) t* = 2.0, (c) t* = 6.0, (d) t* = 10.0, (e) t* = 15.0, and (f) t* = 20.0.

In homogeneous ferroelectric thin films, a-domains generally cross through the thickness of the film and form nearly charge-neutral 90° domain walls to achieve the most stable domain configuration in a tetragonal ferroelectric structure.^57,62 On the other hand, a-domains that do not extend through the thickness are also observed in real samples,⁶³ where defects such as dislocations likely stabilize the partial a-domains.^64–66 In such systems, the a-domains adopt a needle-like shape with straight domain walls and tapered needle point.⁶³ However, the polarization map in Fig. 2f indicates that the domain walls are no longer abrupt but curved, particularly near the needle tip. Such curved domain walls are also clearly observed in experiment.²⁴ This suggests a distinguished characteristic of the partial a-domain in the cgFE thin film.

To provide a comprehensive viewpoint on the formation of the needle-like domain and curved domain walls in correlation with the variation of material compositions, we further investigate the distribution of energy densities in the cgFE thin film. Because of the compositional variation, the material properties are spatially dependent in the thickness direction. This bring about the spatial dependence of the polarization magnitude along the thickness, as shown in Fig. S1 in the ESI.† Since both material coefficients and polarization magnitude are dependent of the coordinate in the thickness direction, these significantly effect the distribution of energy densities, particularly the elastic energy density due to polarization-induced strain and the electrostatic energy density due to the internal depolarization field. The distribution of elastic energy density is presented in Fig. 3a. The elastic energy density is large at the a/c domain walls and in the a-domain. In particular, the elastic energy density highly concentrates at the middle of the cgFE thin film. This is because the electrostrictive coefficient (Q₁₁ + Q₁₂) achieves a highest value at the middle of the cgFE PbZr_0.8Ti_0.2O₃ ⇔ PbZr_0.2Ti_0.8O₃ thin film, where the Ti content x is 0.5 in the PbZr_(1−x)Ti_xO₃ ferroelectric, as shown in Fig. 3b. Such a high concentration of elastic energy further enhances interaction between mechanical and electric fields through the electrostrictive coupling. Since the interaction between mechanical and electric fields is inhomogeneous and particularly high at the middle of the cgFE thin film, the domain wall near the needle-tip departs from the regular a/c domain wall, and thereby, is curved. On the other hand, the inhomogeneous distribution in elastic energy brings about an inhomogeneous distribution in depolarization field through the non-linear electromechanical coupling. As shown in Fig. 3c, the depolarization field in the x₁ direction achieves the highest and lowest fields also near the departed domain wall. Interestingly, the highest and lowest values of the depolarization field are at different altitudes from the bottom surface, which originates from the departed domain wall near the needle-tip, as clearly illustrated in Fig. 3d. The high magnitude of the depolarization field, in turn, causes a shrinkage in the size of the a-domain at the middle of the cgFE thin film and a termination of the a-domain within the film. Therefore, the high concentration of elastic energy at the middle of the thin film gives rise to a deflection of the domain wall from the regular a/c domain wall in the cgFE, while the high concentration of depolarization filed shrinks the a-domain near the needle-tip and terminates the a-domain within the cgFE thin film. These originate from the material variation along the thickness of cgFE thin film that brings about a domain structure distinguished from that in homogeneous ferroelectric thin films.

Fig. 3 (a) Elastic energy landscape. (b) Dependence of the electrostrictive coefficient on the content of Ti in PbZr_(1−x)Ti_xO₃. (c) Depolarization field landscape. (d) Distribution of depolarization field along the x₁ direction at x₃ = 50 nm and x₃ = 60 nm, where the depolarization field achieves maximum and minimum values, respectively.

3.2 Ferroelastic switching of domain structure in the cgFE thin film under electric/mechanical fields

Because of the termination of polarization head-to-tail arrangement, the needle-tip of the a-domain is charged and relatively high in energy, as discussed above. This configuration makes the partial a-domain effectively switchable. To investigate behaviours of the needle-like domain, a small homogeneous electric field is applied along the thickness direction. As a negative electric field is imposed to the cgFE thin film, the upward polar vectors in the c-domain are unfavourable to the applied field. The polar vectors at the needle-tip, where the barrier to polarization rotation is lowest, start to rotate 90° to coincide with the orientation of polar vectors in the a-domain. The switching of polarization under a negative electric field thus simultaneously enlarges and elongates the a-domain, as shown in Fig. 4b. At E = −150 kV cm⁻¹, the needle-tip reaches the top surface of the cgFE thin film, and the domain structure achieves a stripe pattern, where all polarization vectors are in the head-to-tail arrangement, as presented in Fig. 4c. In contrast to the stripe domain structure with straight domain walls in homogeneous ferroelectric thin films, the domain configuration in the cgFE thin film is composed of curved domain walls. On the other hand, as a positive electric field is applied to the cgFE thin film along the thickness direction, the partial a-domain becomes shorter as the polar vectors rotate 90° to coincide with the electric field. The switching also begins from the needle-tip. With an increasing electric field, the partial a-domain gradually shrinks along the long axis (Fig. 4d). However, the partial a-domain still remains, even up to 500 kV cm⁻¹ (Fig. 4e). The needle-like domain, therefore, can be seen to be ferroelastic switchable between several metastable states under an electric field, including partial domains, stripe domains, or nearly single c-domain.

Fig. 4 Partial ferroelastic domain switching by electric fields: (a) E = 0.0 (kV cm⁻¹), (b) E = −100 (kV cm⁻¹), (c) E = −150 (kV cm⁻¹), (d) E = 350 (kV cm⁻¹) and (e) E = 500 (kV cm⁻¹).

We then investigate the response of the needle-like domain subjected to a mechanical strain along the thickness direction, as shown in Fig. 5. Under compression, the upward polar vectors in the c-domain are unfavourable. The polar vectors at the needle-tip start to rotate 90° to coincide with the orientation of polar vectors in the a-domain. The ferroelastic switching under compression enlarges the a-domain, as shown in Fig. 5b. At ε₃₃ = −2.2%, the needle-tip reaches the top surface of the cgFE thin film, and the domain structure achieves a stripe pattern, as presented in Fig. 5c. Similar to the case with the negative electric field, the stripe domain configuration in the cgFE thin film under compression is composed of curved domain walls. In contrast, under tension, the partial a-domain becomes shorter as the polar vectors rotate 90° to coincide with the applied strain direction. Ferroelastic switching also initiates from the needle tip. With increasing tension, the partial a-domain gradually shrinks along the long axis (Fig. 5d and e). The needle-like domain, thus, can also be switchable between several metastable states under a mechanical field.

Fig. 5 Partial ferroelastic domain switching by mechanical fields: (a) ε₃₃ = 0.0, (b) ε₃₃ = −0.004, (c) ε₃₃ = −0.022, (d) ε₃₃ = 0.01 and (e) ε₃₃ = 0.03.

3.3 Ferroelectric switching of domain structure in the cgFE thin film: a role of flexoelectricity

To explore how the variation in partial a-domain influences the material and domain wall responses, we consider the ferroelectric switching of the cgFE thin film under a time-dependent sinusoidal electric field as E = E_max sin(2πft). The electric field is imposed along the thickness direction of the cgFE thin film. From preliminary simulation, the amplitude of the applied field, E_max, is selected equal to 350 kV cm⁻¹ to guarantee for the ferroelectric switching, and the loading frequency is about 0.1 kHz. Note that the macroscopic hysteresis loop of the cgFE thin film observed in experiment is voltage asymmetric, implying that there is preference for the up- or down-poled states. Previous studies have suggested that the primary energy responsible for the asymmetry is electrostatic energy, however, the underlying mechanism for the asymmetry is still unknown. Since the cgFE thin film intrinsically possesses a steep gradient in polarization and strain magnitudes²³ (see Fig. S1 in the ESI†), we anticipate that the voltage offset or built-in potential results from the presence of the strain and polarization gradient through the flexoelectric effects. To make this point clear, we investigate the influence of the flexoelectric effect on the ferroelectric switching in the cgFE thin film.

We first investigate the influence of f₁₁ and f₄₄ on the polarization response of the cgFE thin film under external electric loading. Two flexocoupling types (f₁₁ ≠ 0, f₁₂ = 0, f₄₄ = 0) and (f₁₁ = 0, f₁₂ = 0, f₄₄ ≠ 0) are separately investigated for the effects of f₁₁ and f₄₄, respectively. The values of f₁₁ and f₄₄ to be considered here vary from 0 V to 11.77 V.⁶⁷ The hysteresis loops obtained from different f₁₁ and f₄₄ are presented in Fig. 6a and b, respectively. The macroscopic hysteresis loop exhibits a rectangular-like shape in the absence of the flexoelectric effect (f₁₁ = f₁₂ = f₄₄ = 0), in which the coercive electric field is of about 200 kV cm⁻¹. However, the hysteresis loop slightly shifts towards the positive side of the electric field with the increase in f₁₁ and f₄₄ values, even at a value as large as 11.77 V. This demonstrates that both the f₁₁ and f₄₄ coefficients slightly affect the ferroelectric switching in the cgFE thin film.

Fig. 6 Electrical hysteresis loops for: f₁₁ = 0 V, 7.58 V, 10.46 V, and 11.77 V under E = 350 sin(2πft) (kV cm⁻¹).

In a similar manner, we continue to investigate the influence of the f₁₂ coefficient on the polarization response of the cgFE thin film under an external electric loading. In this consideration, f₁₁ and f₄₄ are set to be zero, while f₁₂ is considered in a range from 0 V to 10.46 V.⁶⁷ The hysteresis loop gradually shifts towards the positive side of the electric field as the value of f₁₂ increases from 0 to 7.85 V. The shift becomes significant when f₁₂ is about 10.46 V. Such a voltage offset in the hysteresis loop clearly demonstrates the effect of flexoelectricity on the ferroelectric switching in the cgFE thin film, and is consistent with the observation in the previous experiment.²⁴ On the other hand, the flexoelectric coefficient f₁₂ in the homogeneous PZT is evaluated up to 25 V,^48,68 while that in the pure PbTiO₃ is about 5 V.^48,69 Since flexoelectricity is proportional to permittivity,^70–73 the high value of f₁₂ of about 10/12 V for the PbZr_0.8Ti_0.2O₃ ⇔ PbZr_0.2Ti_0.8O₃ thin film, which is roughly as large as the averaged value of the pure PbTiO₃ and PZT, is reasonably used in the present study. We further plot in Fig. 7b the dependent profile of the voltage offset, which is characterized by the difference between the coercive electric fields in the positive (E^P_C) and negative (E^N_C) electric sides of the hysteresis loop, e.g., ΔE = |E^P_C| − |E^N_C|, on the flexoelectric coefficient f₁₂. As f₁₂ is small (f₁₂ ≤ 6 V), ΔE slightly increases with a rise in f₁₂. However, ΔE increases significantly as f₁₂ further increases. This, on the other hand, verifies the unusual and complex evolution of built-in potential in cgFE thin films reported in the previous study.¹⁸

Fig. 7 (a) Electrical hysteresis loops for: f₁₂ = 0 V, 5.23 V, 7.58 V, and 10.46 V under E = 350 sin(2πft) (kV cm⁻¹). (b) The dependence of voltage offset on the flexoelectric coefficient f₁₂. (c) Snapshots of the polarization structures during the temporal evolution of polarization switching under the coercive electric field, from A1 to A6 for the case of f₁₂ = 10.46 V.

On the other hand, to clarify the behaviour of polarization vectors during the switching in the presence of flexoelectricity, typical snapshots of polarization domain structures are presented in Fig. 7c that correspond to points A1–A6 in Fig. 7a for the case of (f₁₂ = 10.46 V, f₁₁ = f₄₄ = 0). At the beginning of the switching (point A1), the polarizations change their orientations at a small local area near the curved domain wall. Such a local change in polarization field originates from an inhomogeneous distribution of electrostatic energy along the curved domain walls. The local change in polarization field interestingly brings about the formation of a local vortex–antivortex pair^74,75 within the a/c/a/c-domain structures. As time passes, the vortex and antivortex increase in size, while the vortex and antivortex cores repel each other, as shown at point A2. At point A3, there are several vortex–antivortex pairs formed in the domain structure at the original domain walls. The stripe domain structure is thus almost faded away, and the domain structure with multiple vortices and antivortices emerges. Note that the polarization domains with polarization vector along the electric field are also formed with the formation of vortex and antivortex. Such domains gradually increase in size during the evolution (point A4), then cover a large portion of domain structures (point A5). Simultaneously, vortices and antivortices are forced to shrink. Continuous distortions of polarization can bring vortices and antivortices together and cause them to annihilate. At the end of evolution, new needle-like domains emerge, and the polarization orientation of the c-domain is switched towards the direction of the applied electric field, as shown by point A6. It is noted that a reversed switching phenomenon can occur in a similar manner when the electric field is applied in the opposite direction, yet at a larger coercive field.

To provide a comprehensive viewpoint on the characteristic of polarization switching in the cgFE thin film, we further investigate the switching in a homogeneous PbTiO₃ ferroelectric thin film with the same dimensions and simulation set-ups. The obtained results are provided in Fig. S2 within the ESI.† A careful consideration of polarization evolution during the switching suggests that there is no vortex or antivortex formed during the switching since the switching initiates homogeneously from the 90° straight domain walls. At the end of the switching, reversed stripe domains with straight domain walls are obtained. Details of the polarization stripe domain switching are also provided in a previous study.⁵⁷ Therefore, in the cgFE thin film, the polarization switching process exhibits vastly different responses in comparison to that in the homogeneous ferroelectric thin film. The local switching initiation and the formation and annihilation of vortex–antivortex pairs are unusual phenomena that characterize the ferroelectric switching in the cgFE thin film.

4 Conclusion

In the present study, a new phase field model is developed for cgFE thin films made of PbZr_(1−x)Ti_xO₃ based on the Ginzburg–Landau theory. The developed approach is validated by comparing the simulation results with several key phenomena occurring in cgFE thin films that are observed from experiment, such as the formation of needle-like domains with curved domain walls and the voltage offset of the hysteresis loop. We demonstrate that the high concentration of elastic energy near the needle-tip of the a-domain gives rise to a deflection of the domain wall from the regular a/c domain wall in the cgFE thin film, while the high concentration of the depolarization field shrinks the a-domain near the needle-tip and terminates the domain within the thin film. As a result, such a high concentration of elastic energy and depolarization field facilitates the stabilization of needle-like domains with curved domain walls in the cgFE thin film. The needle-like domain can be extended upwards to the surface of the film through ferroelastic switching under the application of electric field or compression at the surface of the film. Our results also show that the ferroelastic switching initiates at the high-energy point of the needle, where the polarization rotation is initiated. Furthermore, flexoelectricity is proven to play an important role in the voltage offset of the hysteresis loop. On the other hand, the ferroelectric switching process in the cgFE thin film exhibits vastly different responses in comparison to that in homogeneous ferroelectric thin films, including local switching initiation, and the formation and annihilation of vortex–antivortex pairs during the switching. The present study, therefore, provides an incisive approach for an increasingly important class of materials formed by cgFEs, which may bring new understanding and unique insights of these complex cgFE materials, as well as novel potential applications.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research is funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.02-2018.06.

Notes and references

1. J. F. Scott and C. A. P. De Araujo, Science, 1989, 246, 1400–1405.
2. J. Scott, Science, 2007, 315, 954–959.
3. S. Miller and P. McWhorter, J. Appl. Phys., 1992, 72, 5999–6010.
4. S. Mathews, R. Ramesh, T. Venkatesan and J. Benedetto, Science, 1997, 276, 238–240.
5. E. Y. Tsymbal and H. Kohlstedt, Science, 2006, 313, 181–183.
6. W. Eerenstein, N. Mathur and J. F. Scott, Nature, 2006, 442, 759.
7. L. V. Lich, T. Shimada, K. Miyata, K. Nagano, J. Wang and T. Kitamura, Appl. Phys. Lett., 2015, 107, 232904.
8. T. Choi, S. Lee, Y. Choi, V. Kiryukhin and S.-W. Cheong, Science, 2009, 324, 63–66.
9. M. Brazier, M. McElfresh and S. Mansour, Appl. Phys. Lett., 1999, 74, 299–301.
10. N. W. Schubring, J. V. Mantese, A. L. Micheli, A. B. Catalan and R. J. Lopez, Phys. Rev. Lett., 1992, 68, 1778.
11. J. V. Mantese, N. W. Schubring, A. L. Micheli, A. B. Catalan, M. S. Mohammed, R. Naik and G. W. Auner, Appl. Phys. Lett., 1997, 71, 2047–2049.
12. M. Brazier, M. McElfresh and S. Mansour, Appl. Phys. Lett., 1998, 72, 1121–1123.
13. D. Bao, X. Yao and L. Zhang, Appl. Phys. Lett., 2000, 76, 2779–2781.
14. N. Choudhury, L. Walizer, S. Lisenkov and L. Bellaiche, Nature, 2011, 470, 513.
15. S. Ezhilvalavan and T.-Y. Tseng, Mater. Chem. Phys., 2000, 65, 227–248.
16. S. Zhong, S. Alpay, Z.-G. Ban and J. Mantese, Appl. Phys. Lett., 2005, 86, 092903.
17. G. Akcay, S. Zhong, B. Allimi, S. Alpay and J. Mantese, Appl. Phys. Lett., 2007, 91, 012904.
18. J. C. Agar, A. R. Damodaran, G. A. Velarde, S. Pandya, R. Mangalam and L. W. Martin, ACS Nano, 2015, 9, 7332–7342.
19. X. Zhu, N. Chong, H. L.-W. Chan, C.-L. Choy, K.-H. Wong, Z. Liu and N. Ming, Appl. Phys. Lett., 2002, 80, 3376–3378.
20. X. Zhu, J. Zhu, S. Zhou, Z. Liu, N. Ming, H. L.-W. Chan, C.-L. Choy, K.-H. Wong and D. Hesse, Mater. Sci. Eng., B, 2005, 118, 219–224.
21. D. Maurya, F.-C. Sun, S. P. Alpay and S. Priya, Sci. Rep., 2015, 5, 15144.
22. R. Mangalam, J. Karthik, A. R. Damodaran, J. C. Agar and L. W. Martin, Adv. Mater., 2013, 25, 1761–1767.
23. A. R. Damodaran, S. Pandya, Y. Qi, S.-L. Hsu, S. Liu, C. Nelson, A. Dasgupta, P. Ercius, C. Ophus and L. R. Dedon, et al., Nat. Commun., 2017, 8, 14961.
24. J. Agar, A. Damodaran, M. Okatan, J. Kacher, C. Gammer, R. Vasudevan, S. Pandya, L. Dedon, R. Mangalam and G. Velarde, et al., Nat. Mater., 2016, 15, 549.
25. X. Wang, C. Wang, W. Zhong and P. Zhang, Phys. Lett. A, 2001, 285, 212–216.
26. H.-X. Cao and Z.-Y. Li, J. Phys.: Condens. Matter, 2003, 15, 6301.
27. J. Karthik, R. V. K. Mangalam, J. C. Agar and L. W. Martin, Phys. Rev. B, 2013, 87, 024111.
28. J. Zhang, R. Xu, A. R. Damodaran, Z.-H. Chen and L. W. Martin, Phys. Rev. B, 2014, 89, 224101.
29. I. Misirlioglu and S. Alpay, Acta Mater., 2017, 122, 266–276.
30. Y. Qiu, H. Wu, J. Wang, J. Lou, Z. Zhang, A. Liu and G. Chai, J. Appl. Phys., 2018, 123, 084103.
31. H.-L. Hu and L.-Q. Chen, J. Am. Ceram. Soc., 1998, 81, 492–500.
32. Y. Li, S. Hu, Z. Liu and L. Chen, Appl. Phys. Lett., 2001, 78, 3878–3880.
33. S. Hu and L. Chen, Acta Mater., 2001, 49, 1879–1890.
34. L.-Q. Chen, Annu. Rev. Mater. Res., 2002, 32, 113–140.
35. L. V. Lich, T. Shimada, K. Nagano, Y. Hongjun, J. Wang, K. Huang and T. Kitamura, Acta Mater., 2015, 88, 147–155.
36. L. V. Lich, T. Shimada, S. Sepideh, J. Wang and T. Kitamura, Acta Mater., 2017, 125, 202–209.
37. L. V. Lich, T. Shimada, J. Wang, V.-H. Dinh, T. Q. Bui and T. Kitamura, Phys. Rev. B, 2017, 96, 134119.
38. G. Shirane, R. Pepinsky and B. C. Frazer, Acta Crystallogr., 1956, 9, 131–140.
39. J. Wang, S.-Q. Shi, L.-Q. Chen, Y. Li and T.-Y. Zhang, Acta Mater., 2004, 52, 749–764.
40. M. J. Haun, E. Furman, S. Jang, H. McKinstry and L. Cross, J. Appl. Phys., 1987, 62, 3331–3338.
41. M. J. Haun, E. Furman, S. J. Jang and L. E. Cross, Ferroelectrics, 1989, 99, 13–25.
42. S. Choudhury, Y. Li and L.-Q. Chen, J. Am. Ceram. Soc., 2005, 88, 1669–1672.
43. Y. Li, S. Hu and L. Chen, J. Appl. Phys., 2005, 97, 034112.
44. M. J. Haun, E. Furman, H. A. McKinstry and L. E. Cross, Ferroelectrics, 1989, 99, 27–44.
45. M. J. Haun, Z. Q. Zhuang, E. Furman, S. J. Jang and L. E. Cross, Ferroelectrics, 1989, 99, 45–54.
46. W. Cao and L. E. Cross, Phys. Rev. B, 1991, 44, 5.
47. N. Pertsev, V. Kukhar, H. Kohlstedt and R. Waser, Phys. Rev. B, 2003, 67, 054107.
48. P. Zubko, G. Catalan and A. K. Tagantsev, Annu. Rev. Mater. Res., 2013, 43, 387–421.
49. P. Yudin and A. Tagantsev, Nanotechnology, 2013, 24, 432001.
50. Q. Li, C. Nelson, S.-L. Hsu, A. Damodaran, L.-L. Li, A. Yadav, M. McCarter, L. Martin, R. Ramesh and S. Kalinin, Nat. Commun., 2017, 8, 1468.
51. Y. Cao, A. Morozovska and S. V. Kalinin, Phys. Rev. B, 2017, 96, 184109.
52. J. Wang, Y. Li, L.-Q. Chen and T.-Y. Zhang, Acta Mater., 2005, 53, 2495–2507.
53. R. D. Mindlin, Int. J. Solids Struct., 1965, 1, 417–438.
54. Y. L. Li, S. Y. Hu, Z. K. Liu and L. Q. Chen, Acta Mater., 2002, 50, 395–411.
55. P. Song, T. Yang, Y. Ji, Z. Wang, Z. Yang, L. Chen and L. Chen, Comput. Phys. Commun., 2017, 21, 1325–1349.
56. W. Shu, J. Wang and T.-Y. Zhang, J. Appl. Phys., 2012, 112, 064108.
57. L. V. Lich, T. Shimada, J. Wang and T. Kitamura, Acta Mater., 2016, 112, 1–10.
58. M. J. Haun, E. Furman, T. R. Halemane and L. E. Cross, Ferroelectrics, 1989, 99, 55–62.
59. M. J. Haun, E. Furman, S. J. Jang and L. E. Cross, Ferroelectrics, 1989, 99, 63–86.
60. P. Gao, J. Britson, J. R. Jokisaari, C. T. Nelson, S.-H. Baek, Y. Wang, C.-B. Eom, L.-Q. Chen and X. Pan, Nat. Commun., 2013, 4, 2791.
61. F. Xue, J. J. Wang, G. Sheng, E. Huang, Y. Cao, H. H. Huang, P. Munroe, R. Mahjoub, Y. L. Li and V. Nagarajan, et al., Acta Mater., 2013, 61, 2909–2918.
62. B. Meyer and D. Vanderbilt, Phys. Rev. B, 2002, 65, 104111.
63. P. Gao, J. Britson, C. T. Nelson, J. R. Jokisaari, C. Duan, M. Trassin, S.-H. Baek, H. Guo, L. Li and Y. Wang, et al., Nat. Commun., 2014, 5, 3801.
64. D. Su, Q. Meng, C. Vaz, M.-G. Han, Y. Segal, F. J. Walker, M. Sawicki, C. Broadbridge and C. H. Ahn, Appl. Phys. Lett., 2011, 99, 102902.
65. M.-W. Chu, I. Szafraniak, R. Scholz, C. Harnagea, D. Hesse, M. Alexe and U. Gösele, Nat. Mater., 2004, 3, 87.
66. T. Kiguchi, K. Aoyagi, Y. Ehara, H. Funakubo, T. Yamada, N. Usami and T. J. Konno, Sci. Technol. Adv. Mater., 2011, 12, 034413.
67. P. Yudin, R. Ahluwalia and A. Tagantsev, Appl. Phys. Lett., 2014, 104, 082913.
68. W. Ma and L. E. Cross, Appl. Phys. Lett., 2003, 82, 3293–3295.
69. G. Catalan, A. Lubk, A. Vlooswijk, E. Snoeck, C. Magen, A. Janssens, G. Rispens, G. Rijnders, D. Blank and B. Noheda, Nat. Mater., 2011, 10, 963.
70. S. M. Kogan, Soviet Physics - Solid State, 1964, 5, 2069–2070.
71. G. Catalan, L. Sinnamon and J. Gregg, J. Phys.: Condens. Matter, 2004, 16, 2253.
72. L. E. Cross, J. Mater. Sci., 2006, 41, 53–63.
73. P. Zubko, G. Catalan, A. Buckley, P. R. L. Welche and J. F. Scott, Phys. Rev. Lett., 2007, 99, 167601.
74. I. I. Naumov, L. Bellaiche and H. Fu, Nature, 2004, 432, 737–740.
75. L. V. Lich, T. Shimada, J. Wang and T. Kitamura, Nanoscale, 2017, 9, 15525–15533.

---

Footnote

† Electronic supplementary information (ESI) available: Supplement figures; brief description of the polarization switching in homogeneous ferroelectric under electric field. See DOI: 10.1039/c8ra10614b

---