The morphologies of twins in tetragonal ferroelectrics

Twins, as a special structure, have been observed in a small number of experiments on ferroelectrics. A good understanding of the morphologies of twins is very important for twin engineering applications in ferroelectric materials. In this work, the morphologies of twins in tetragonal ferroelectrics are investigated using the compatibility analysis of the transformation strains and spontaneous polarization and the energetic analysis. The tetragonal BaTiO3 single crystal is chosen as an example for the material system. The results show that the lamellar twin structures with 111 as the twin plane have been identified by both compatibility analysis and energetic analysis, which is consistent with the experimental observations. In addition to 111 twin structures, 1 2 ¯ 1 and 2 1 ¯ 5 twin structures can also appear in tetragonal ferroelectrics. Furthermore, stable state uncharged twin boundaries and metastable charged twin boundaries are distinguished by energetic analysis, where the spontaneous polarization is compatible across the uncharged twin boundaries, while it is not compatible across the charged twin boundaries. The results predicted by the compatibility analysis and energetic analysis are consistent, thus providing a pathway for understanding the morphologies of twins.

across which the strain compatibility condition is satisfied, but the polarization compatibility condition is not satisfied, have been reported in the theoretical and experimental aspects (Liu et al., 2007;Liu and Li, 2009;Li et al., 2016). Therefore, the formation of domain structure is the basis for understanding the functional properties of ferroelectric materials. It is also the technological application of underpinned ferroelectric materials, making domain engineering being an effective way to enhance the performance of ferroelectrics (Shelke et al., 2011;Li et al., 2017;Geng et al., 2020;Qiu et al., 2020;Chen et al., 2022;Lei and Liu, 2022;Liu et al., 2022;Xu et al., 2022).
Domain engineering for ferroelectric materials has been extensively studied from both experimental and theoretical perspectives (Shelke et al., 2011;Li et al., 2017;Geng et al., 2020;Qiu et al., 2020;Chen et al., 2022;Lei and Liu, 2022;Liu et al., 2022;Xu et al., 2022). It is well known that twin structures are often observed in metals and alloys (Li et al., 2018b;Song et al., 2020;Li et al., 2023). Based on the formation of twins, twin engineering is widely used to enhance the mechanical properties of materials (Cheng et al., 2018;Li et al., 2022;Chen et al., 2023), thermal properties (Gao et al., 2021), electrical properties (Lei et al., 2018). Although domain engineering in ferroelectrics has received extensive attention, twins in ferroelectrics have rarely been reported. It is noticed that the 90°ferroelectric domain in some literature is called "twin", but it is not a twin structure whose lattice should be symmetrical about the twin planes.
Twin structures have been experimentally observed in BaTiO 3 powders (Eibl et al., 1987;Qin et al., 2010). For example, Qin et al. synthesized twinned BaTiO 3 powders using BaCl 2 and TiO 2 at low temperatures, and their twin plane is (111) plane (Qin et al., 2010). By calcining a powder mixture of BaCO 3 and TiO 2 at 1100°C, Eibl et al. observed a (111) twin structures in the BaTiO 3 ceramics (Eibl et al., 1987). Also, Cao et al. fabricated twinned structures in BaTiO 3 films by a thermomechanical treatment technique, and BaTiO 3 films with twin structures exhibit a smaller coercive field (Cao et al., 2015). This suggests that twin engineering may be used to enhance the properties of ferroelectrics. However, there are currently few reports on ferroelectric twinning, while the morphologies of ferroelectric twin, including its grain shape and orientation remain unclear. In addition to the (111) twin plane, it is still unclear whether ferroelectric materials can form other twin planes, while the revealing of morphologies and twin planes for ferroelectric twins is very important for utilizing twin engineering to enhance the performance of ferroelectrics. It is well known that the stable structures in ferroelectrics are energetically stable states, which should satisfy compatibility conditions on transformation strains and spontaneous polarization. Our previous work developed an energetic analysis method for ferroelectric domain patterns based on the equivalent inclusion method (Liu et al., 2007;Liu and Li, 2009). Thus, in this article, we analyze the morphologies of twins in tetragonal ferroelectrics by using compatibility analysis and energetic analysis, and the energy minimizing twin morphologies, including orientation and shape, are identified, as well as the determined twin planes.

Tetragonal twin variants
In experiments, (111) twins have been observed in BaTiO 3 powders and films at room temperature (Eibl et al., 1987;Qin et al., 2010;Cao et al., 2015;Cao et al., 2017). It is well known that BaTiO 3 is tetragonal at room temperature. Guided by the experimental results (Eibl et al., 1987;Qin et al., 2010;Cao et al., 2015), the schematic of a tetragonal twin structure is plotted in Figure 1A, where the tetragonal unit cells TW1 and TW2 are symmetric about the twin plane. Each tetragonal unit cell has six possible polarization directions, corresponding to the six ferroelectric variants. Thus, there are a total of 12 ferroelectric variants in tetragonal twinned ferroelectrics. In their respective cubic crystallographic axes, such as TW1 in its local coordinate system x TW1 1 x TW1 2 x TW1 3 and TW2 in its local coordinate system x TW2 1 x TW2 2 x TW2 3 , the transformation strain ε i (i = 1,. . .,6) and the spontaneous polarization P i of each variant are given as: where α and β are lattice constants of the materials, and P t is the module of the spontaneous polarization. It is noticed that the variants with opposite spontaneous polarizations have the same transformation strains. Considering that the tetragonal unit cells TW1 and TW2 have different local crystallographic axes, they have different transformation strains and spontaneous polarizations, which can be obtained by the tensor rotation of Eq. 1, which are denoted by (ε S−TW1

Compatibility condition
In ferroelectrics, the stable interface between various ferroelectric variants is a consequence of energy minimization, across which the transformation strains and spontaneous polarizations of the variants satisfy the compatibility conditions (Shu and Bhattacharya, 2001;Li and Liu, 2004;Liu and Li, 2009;O'Reilly et al., 2022). For twin structures, the variants quantified by (ε S−TW1 form a twin boundary, as shown in Figure 1B. Thus, the compatibility conditions on the transformation strains and the spontaneous polarizations across the twin boundary are given by: where n represents the normal to the twin boundary, a represents the shear vector along the TB, and ⊗ refers to the standard tensor product of two vectors. Thus, the compatibility conditions Eqs 2, 3 allow us to determine the variants (ε S−TW1 i , P S−TW1 i ) and (ε S−TW2 j , P S−TW2 j ) that form an energetically stable twin boundary plane.

Energetic of morphology
While the compatibility analysis is capable of identifying the orientation of twin structures, it cannot identify the grain shapes of the twins, as well as the charged interface, which is energetically metastable due to polarizations incompatibility (Liu et al., 2007;Liu and Li, 2009). These deficiencies can be overcome by energy analysis using the equivalent inclusion method, which has been successfully applied to analyze the morphologies of ferroelectric domain patterns (Liu et al., 2007) and the precipitates in thermoelectric compounds (Liu et al., 2014).
For a ferroelectric crystal possessing a transformation strain ε s and spontaneous polarization P s , the electromechanical behavior is governed by the following constitutive equations: where ε and σ are the strain and stress, respectively. D and E are the electric displacement and field, respectively. F E is the elastic compliance tensor under constant electric field, k σ is the dielectric constant under constant stress, d is the piezoelectric coefficient. The superscript t denotes the transpose of the tensor. By change of variables, the ferroelectric constitutive Eqs 4, 5 can be transformed into: the two types of constitutive Eqs. 4-7 can be rewritten as the following matrix forms: with the superscript −1 denoting matrix inverse. Since the electromechanical moduli M and L are positive definite, it is convenient to analyze the energy of twin structures using the above constitutive Eqs. 4-7.
In order to analyze the stable morphologies of stable twin structures in tetragonal ferroelectrics, the equivalent inclusion method is adopted (Dunn and Taya, 1993;Liu et al., 2007;Liu and Li, 2009). The twin structure is treated as a matrix region D Frontiers in Materials frontiersin.org and an inhomogeneous inclusion Ω, where the matrix is composed of TW1 with electromechanical moduli denoted by L TW1 and transformation strain and spontaneous polarization given by Y S−TW1 , while the inhomogeneous inclusion is composed of TW2 with L TW2 and Y S−TW2 , as shown in Figure 2. Using the matrix TW1 as a reference, the effective transformation strain and spontaneous polarization of the inhomogeneous inclusion TW2 can be written as: the electromechanical fields can be obtained by Eshleby's classical solution (Liu et al., 2007;Liu and Li, 2009). Consider a uniform external field X 0 is applied to the twin structure in Figure 2. If an inhomogeneous inclusion does not exist, a uniform Y 0 will be induced in matrix. But due to the emergence of inhomogeneous inclusion TW2, the additional disturbance fields X′ and Y′ are induced. Thus, the electromechanical field for the inhomogeneous inclusion TW2 can be expressed by (Liu and Li, 2009): which can be computed by Eshelby's equivalent inclusion method (Liu et al., 2007;Liu and Li, 2009), such as: where the inhomogeneous inclusion TW2 is replaced by an equivalent inclusion, which has the same electromechanical moduli L TW1 of matrix and additional transformation strain and spontaneous polarization Y**, which ensures the equivalency of electromechanical fields between TW2 and the equivalent inclusion. For ellipsoidal inclusion, it is well known that the disturbance electromechanical field in the inclusion can be calculated by (Wang, 1992;Dunn and Taya, 1993;Liu et al., 2007;Liu and Li, 2009): where in the above expressions, I is the unit four-rank tensor, S is the piezoelectric Eshelby's tensors, which is a function of the electromechanical modulus L TW1 of the matrix TW1, as well as a function of the shape and orientation of the inclusion TW2 (Dunn and Taya, 1993). Combining Eqs 12, 13, the fictitious transformation strain and spontaneous polarization Y ** can be expressed as from which the electromechanical field in the twin structure can be determined. The potential energy of the twin structure system can then be expressed by:

FIGURE 2
Schematic of the inhomogeneous inclusion, where TW1 and TW2 represent the unit cells in the twin structure.
where the first term represents the elastic and electric energies, and the second term represents the work done by external loading. In the absence of the inhomogeneous inclusion TW2, the effective transformation strain and spontaneous polarization Y S and the disturbance fields vanish, leading to the potential energy expressed by: Thus, due to the existence of the inhomogeneous inclusion TW2, which infers the formation of a twin structure, the resulting energy variation is where the Hill's condition is adopted for simplification (Dunn and Taya, 1993;Li and Dunn, 1999). In this work, we consider the scenario of no external field, and the energy variation is further reduced to: Eq. 18 allows us to analyze the energy variation due to the emergence of twin structures, including its orientation and shape of the inhomogeneous inclusion TW2. The equilibrium morphologies of the twins can be identified by minimizing the energy variation with respect to the orientation and shape of the inhomogeneous inclusion TW2. For the inhomogeneous inclusion TW2, its shape can be assumed to be ellipsoidal, and the orientation can be described using Euler angles in the global coordinate system.

Results and discussions
In order to determine the morphologies of twins in tetragonal ferroelectrics, we employ the above theory to implement the compatibility analysis and energetic analysis of twin structures. We choose BaTiO 3 as the material system, whose electromechanical moduli associated with the single domain with the polarization along [001] in its cubic crystallographic axes is M. The transformation strain and spontaneous polarization are listed in Table 1. Notice that the physical quantities of other variants in different coordinate systems can be obtained by tensor rotations.

Compatibility analysis
Guided by the experimental results (Eibl et al., 1987;Qin et al., 2010;Cao et al., 2015), the tetragonal twin structure is schematically illustrated in Figure 1A x TW1 3 coordinate is then given by: In the x TW1  in its cubic crystallographic axes (F E : 10 -12 m 2 /N; d: 10 -12 C/N; κ σ /κ 0 ; where κ 0 = 8.85 × 10 -12 C 2 /Nm 2 is the permittivity of free space), transformation strain ε; and spontaneous polarization P(C/m 2 ) of tetragonal BaTiO 3 single crystal (Zgonik et al., 1994;Liu and Li, 2009 resulting in: For the twin structure composed of TW1 (ε S−TW1 i , P S−TW1 i ) and TW2 (ε S−TW2 j , P S−TW2 j ), the normal to the twin plane can be computed by solving Eqs. 2 and 3. Based on this compatibility analysis of the transformation strains and spontaneous polarizations, the twin planes between different tetragonal variants from TW1 and TW2 that may form twin structures are summarized in Table 2. It is found that the following pairs of variants can form (111)-twin in tetragonal ferroelectric: variants ), which have also been observed in experiments (Eibl et al., 1987;Cheng et al., 2006;Wu et al., 2006;Qin et al., 2010;Cao et al., 2015). Three representatives of these (111)-twins are schematically shown in Figure 3. Additionally, Table 2 shows that 121 and 215 twinned structures may also exist in tetragonal ferroelectrics, which have not been observed in currently reported experiments.     indicate the local coordinate system along cubic crystallographic axes of TW1 and TW2, respectively. ) as the matrix; (C) the potential energy variation as a function of the shape aspect ratio a 3 /a 1 and orientation angle θ of the inhomogeneous inclusion TW2; and (D) the corresponding potential energy variation at a 3 /a 1 = 0.0001; The schematics of the uncharged (E) and charged (F) twin boundaries corresponding to the two wells in the energy landscape. ), for the convenience of calculation, we choose the global coordinate system x 1 x 2 x 3 as shown in Figure  5A, where

Energetic analysis
TW1 and x 3 0 0 1 G 0 0 1 TW1 , which makes the polarization vectors of the two variants lie on the same plane denoted by HJOQ, as shown in Figure 5B. The orientation of inhomogeneous inclusion TW2 can then be described by only one Euler angle θ, which is the angle between the x 3 axis and the rotational axis X 3 of the spheroid, as shown in Figure 5B. Notice that x 2 X 2 . The transformation matrix that rotates from local coordinates to the global x 1 x 2 x 3 coordinate are respectively: The shape aspect ratio of the inhomogeneous inclusion TW2 can be defined by a 3 /a 1 a 3 /a 2 , where a 1 , a 2 , a 3 are the dimensions along the principal axes of the spheroidal. Based on Eq. 18, we calculate the potential energy variation as a function of the shape aspect ratio a 3 /a 1 and the orientation angle θ of inhomogeneous inclusion TW2, as shown in Figure 5C. It is found that, when the aspect ratio a 3 /a 1 approaches zero, two energy wells exist in the energy landscape, indicating that the stable twin structure is lamellar. Indeed, lamellae twins have been observed in BaTiO 3 ceramics (Oppolzer and Schmelz, 1983). Furthermore, one of the energy wells occurs at θ 125.3°, while another occurs at θ 36.6°. The well at θ 125.3°refers to the uncharged twin boundary with (111) as the twin plane, across which the transformation strains and spontaneous polarizations are compatible, as shown in Figure 5E, which is consistent with those predicted by the compatibility analysis [as listed in Table 2]. On the other hand, the well at θ 36.6°refers to the charged twin boundary with (112) as the twin plane, across which only the transformation strain is compatible, but the spontaneous polarization is not compatible, as shown in Figure 5F. Figure 5D reveals that the energy well at θ 125.3°i s lower than that at θ 36.6°, suggesting that the uncharged twin boundary is in a stable state, while the charged twin boundary is in a metastable state. A steep energy barrier exists between the two energy wells, making it hard to convert from the charged twin boundary to the uncharged twin boundary. The charged domain walls have been observed in many experiments (Sluka et al., 2012;Li et al., 2016;Wu et al., 2021), which are similar to the charged twin boundaries. Both of them are in ) as the matrix; (C) the potential energy variation as a function of the shape aspect ratio a 3 /a 1 and orientation angle θ of the inhomogeneous inclusion TW2; and (D) the corresponding potential energy variation at a 3 /a 1 = 0.0001; The schematic of the uncharged (E) twin boundary corresponding to the well in the energy landscape.
Frontiers in Materials frontiersin.org 08 metastable states, across which the spontaneous polarizations are also not compatible. ), for the convenience of calculation, we choose the global coordinate system x 1 x 2 x 3 as shown in Figure  6A, where

Twinned structure with variants
, and x 3 0 0 1 G 0 0 1 TW1 , such that polarization vectors of the two variants lie on the same plane HJOQ, as shown in Figure 6B. For this case, the potential energy variation due to the emergence of the inhomogeneous inclusion TW2 as a function of its shape aspect ratio a 3 /a 1 and the orientation angle θ is plotted in Figure 6C. It is observed that one energy well occurs at θ 114.1°in the energy landscape with a 3 /a 1 approaching zero, which corresponds to the uncharged twin boundary with (121) as the twin plane, across which the transformation strains and spontaneous polarizations are compatible, as shown in Figure 6E, consistent with the prediction by the compatibility analysis [as listed in Table 2]. ), we choose the global coordinate x 1 x 2 x 3 , as shown in Figure 7A, where x 1 1 0 0 G − 2

Twinned structure with variants
, and x 3 0 0 1 G 0 0 1 TW1 . In the global coordinate x 1 x 2 x 3 , the polarization vectors of the two variants lie on the same plane HJOQ, as shown in Figure 7B. We calculate the potential energy variation due to the emergence of the inhomogeneous inclusion TW2 as a function of its shape aspect ratio a 3 /a 1 and orientation angle θ in Figure 7C. It can be seen that, when the aspect ratio a 3 /a 1 approaches zero, there exists two energy wells, at θ 155.9°F ) as the matrix; (C) the potential energy variation as a function of the shape aspect ratio a 3 /a 1 and orientation angle θ of the inhomogeneous inclusion TW2; and (D) the corresponding potential energy variation at a 3 /a 1 = 0.0001; The schematics of the uncharged (E) and charged (F) twin boundaries corresponding to the two wells in the energy landscape. and θ 69.3°respectively. The energy well at θ 155.9°c orresponds to the uncharged twin boundary with (215) as the twin plane. The transformation strains and spontaneous polarizations across this twin plane are compatible, as shown in Figure 7E, which is consistent with the prediction by the compatibility analysis [as listed in Table 2]. The other energy well located at θ 69.3°is metastable, resulting in the charged twin boundary with (211) as the twin plane, across which only the transformation strains are compatible, while the spontaneous polarizations are not compatible, as shown in Figure 7F. ), the predicted results of the compatibility analysis and energetic analysis are summarized in Table 3. It can be found that the compatibility analysis on transformation strains can identify twin planes, but it cannot distinguish whether the twin plane is uncharged or charged. The energetic analysis not only can identify the twin plane, but also can further determine whether the twin plane is uncharged in a steady state or charged in a metastable state. The compatibility analysis on both transformation strains and spontaneous polarizations can only determine uncharged twin planes, as charged twin planes do not satisfy the compatibility in spontaneous polarizations. Overall, the results from compatibility analysis and energetic analysis are consistent. It has been demonstrated that twin engineering can reduce the coercive field of ferroelectric materials (Cao et al., 2015) due to the increased ferroelectric variants, suggesting that twin engineering may be used to enhance the other properties of ferroelectrics.

Conclusion
We have investigated the morphologies of twins in tetragonal BaTiO 3 single crystals by using the compatibility analysis of the transformation strains and spontaneous polarizations and the energetic analysis. It is found that the lamellar twin structures with 111 { } as the twin plane have been identified by both compatibility analysis and energetic analysis, which agrees well with experimental observations. The results also show that the 121 and 215 twin structures can exist in tetragonal ferroelectrics. In addition, the stable state uncharged twin boundaries and metastable charged twin boundaries are distinguished by energetic analysis, since the spontaneous polarizations are compatible across the uncharged twin boundaries, while the spontaneous polarizations are not compatible across the charged twin boundaries. The results predicted by the compatibility analysis and energetic analysis agree well with each other, offering a pathway for understanding the morphologies of twins in ferroelectrics.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, upon reasonable request.