Mechanism of nucleation in ferroelastic domain switching

The mechanism of domain nucleation during ferroelastic domain switching is revealed by utilizing phase-field simulations. A complete compression-tension hysteresis loop is simulated, and the results show the nucleation of new domains at domain walls (DW). Due to relaxation of strain energy, the coherent DWs of the surviving domains often form features like steps, crevices and corners, and nucleation of new domains primarily takes place at these sites in a repeated fashion. The change in elastic strain energy around these sites is sufficient to provide the energy necessary for ferroelastic domain nucleation (FDN). The present simulations reveal that the apparent ellipsoidal shape of a nucleus might be an oblate spheroid when the nucleation takes place at a crevice. The present work also reveals that nucleation of triangular shaped domains is equally probable, and they take place at the steps or corners of DWs.

The mechanism of domain nucleation during ferroelastic domain switching is revealed by utilizing phase-field simulations.A complete compression-tension hysteresis loop is simulated, and the results show the nucleation of new domains at domain walls (DW).Due to relaxation of strain energy, the coherent DWs of the surviving domains often form features like steps, crevices and corners, and nucleation of new domains primarily takes place at these sites in a repeated fashion.The change in elastic strain energy around these sites is sufficient to provide the energy necessary for ferroelastic domain nucleation (FDN).The present simulations reveal that the apparent ellipsoidal shape of a nucleus might be an oblate spheroid when the nucleation takes place at a crevice.The present work also reveals that nucleation of triangular shaped domains is equally probable, and they take place at the steps or corners of DWs.Domain switching is the cause of several attractive functional properties in ferroic materials.Under the application of appropriate external stimulus such as an electric field in case of ferroelectric or a stress field in ferroelastic materials, energetically favorable domains nucleate and subsequently grow resulting in switching of domain/s inducing the reversal or reorientation of spontaneous polarization or strains, respectively.Several studies on the mechanism of domain switching have been conducted particularly in ferroelectric materials, in which new domains were observed to nucleate primarily at twodimensional defects such as domain walls (DWs) [1], free surfaces [2], interfaces between two phases [3], and grain boundaries [4].In contrast, in situ study of ferroelastic domain nucleation (FDN) in purely ferroelastic materials has not received much attention.The ex-situ diffraction and electron microscopic studies [5][6][7][8], and in-situ experimental studies of ferroelastic deformation [9][10][11] so far have not been able to reveal the mechanism of FDN.
Early attempt of modelling the domain switching in ferroelectric materials assumed various shapes of nucleus.For instance, Merz assumed needle shaped domain [12] while Miller theorized triangular shape of the nucleus [13].Some ab initio and MD calculations have shown that the shape of nucleus should be square or diamond [14,15].Several works have assumed spheroid shaped nucleus [16][17][18].Some of these shapes are also observed in experiments.For instance in LiNbO 3 [19] and LiTaO 3 [20], triangle shaped domains are observed while in situ STEM observation of BiFeO 3 film has suggested that the nucleus should rather be an oblate spheroid [2].To date, it is not understood exactly why different shapes of nuclei form during nucleation despite the fact that FDN in ferroelectric materials was reported to influence the performance of sensors and actuators [21,22].
Phase-field model (PFM) has significantly contributed to the study of ferroic behavior by simulating realistic domain switching in several ferroelectric materials [1,[23][24][25][26][27].However, in this model, DW energy is attributed to gradient of polarization and the contribution of interfacial energy has not been given due consideration.In the present work, we study FDN under compression-tension cyclic loading in ferroelastic t'phase of yttria-stabilized-zirconia (YSZ) using a quantitative phase-field model which exclusively considers twin boundary as its DW.Earlier, FDN has been proposed to be responsible for toughening behavior of t'phase [28].Using this model, we unravel the mechanism of FDN which has great scientific and technological implications.
In the present PFM, three crystallographic variants of t'phase having c-axis of tetragonal unit cell directed along [001], [010] and [100], are characterized by three order parameters p i (where i = 1, 2 & 3) and each order parameter varies from 0 to 1.The total free energy functional F(p i , ∇p i , ε) of a microstructure consisting of homogeneous domains and domain walls among them, is the sum of free energy of the bulk region consisting of single variant f bulk (p i ), free energy f grad (∇p i ) due to DW and that for elastic deformation f elast (ε) where ε is the strain tensor.The following form of f bulk and f grad are adopted in this work: k p are the gradient energy coefficients, and A, B and C are the constants.We note that the total volume of all the three domains remains conserved during domain switching.Therefore, the present model ensures the conservation of order parameters by setting a very high value for the coefficient A. In our previous work, we have shown that the symmetric phase field profiles of each order parameter defined by p(x) = 1/2[1 ± tanh(x/W)] automatically satisfies the conservation of order parameters [29].The DW thickness (W) and energy (γ) are related to coefficient B and k p by the following relations W = 2k p /√B and γ = /3 and the present model treats the DW characteristics as its input.Finally, the free energy of elastic deformation is given as where C is the fourth order elastic modulus tensor.The time dependent Ginzburg Landau equation for evolution of every order parameter is the following: where M p is the mobility.In the present model, the elastic strain is given by ε = ε tot − ε * = ε hom + δε − ε * where each variant is characterized by stress-free strain ε * which is measured with respect to a cubic lattice having the same volume of the tetragonal unit cell of t'phase.The total strain ε tot is the sum of the spatially independent homogeneous strain ε hom and inhomogeneous strain δε.Linear elasticity is adopted in the present model with an assumption of small deformation and mechanical equilibrium is allowed to prevail i.e., ∇.σ = 0 where σ is the Cauchy stress tensor.The details of 2D numerical setup and the simulation parameters can be found in our previous work where we have shown that the present model simulates realistic ferroelastic deformation in YSZ single crystal [29].In the present work, we have used the same material data for t' phase of YSZ at 500 • C [29], including the starting microstructures observed in experiments [30].The simulation domain in the present work has the dimensions of 720 × 720 nm 2 .In the present work, we have chosen W = Δx = 3 nm where Δx is the spatial resolution.In the Supplementary Materials, we have shown that such a choice of meshing does not affect the nucleation behavior simulated.Due to lack of reliable data, DW energy (γ) of 0.2 J.m − 2 at room temperature was considered in previous works.But in the present work, we have applied γ = 0.1 and 0.2 J.m − 2 which covers the entire range of coherent interfacial energy [31].It should be noted that below γ = 0.1 J.m − 2 , spurious domain switching is observed and hence lower γ is not studied.The strain rate for the applied deformation is 10 − 5 s − 1 .The 2D simulations are consistent with the applications of ferroic behavior in thin coating [32].The chosen microstructures consist of c (c-axis lying normal to x 1 x 2 plane), a 1 (c-axis along x 1 ) and a 2 (c-axis along x 2 as illustrated in Fig. 4a) domains which evolve in response to applied loading.
As represented by red colored plots in Fig. 1, due to application of compression along x 2 on both microstructures, a 2 domain switches to c and a 1 domains.The end of initial compression is the starting point of the hysteresis loop, marked by stage 1 in Fig 1, where the direction of loading is reversed.Between stages 1 and 2, the nucleation and growth of a 2 domains occur, and the nucleation sites are marked by red colored circles in both microstructures.Due to tensile loading (represented by green arrows), until stage 6 in Fig. 1a and stage 5 in Fig. 1b, c domains disappear while some a 1 domains remain unswitched.After these stages, the load is again reversed, and c domains reappear at stage 8 in Microstructure 1 and at stage 7 in Microstructure 2. By further compression until stage 11 (represented by red arrows), c and a 1 domains dominate over a 2 domains.The complete compression-tension hysteresis loop in each initial microstructure generates two different microstructures at the end of each loading path: one consisting of c and a 1 domains at the end of compression while the other one contains predominantly a 2 with some unswitched a 1 domains on completion of tensile loading.During compression of initial microstructures, domain switching occurs by movement of DWs and no nucleation takes place.The nucleation of a 2 domain is observed to occur at DW between c and a 1 domains along the tension path in the hysteresis loop, and in the present work we investigate this phenomenon.
It is well known that alternative arrangements of domains of two different tetragonal variants form DWs or twin boundaries along {110} plane [33,34].At a DW, spontaneous strains of two adjacent domains must satisfy the compatibility condition and this results in the formation of non-smooth features in very fine length scales which are otherwise smooth in microscopic scale [35].High resolution microscopic observations have indeed reported zigzag features of DWs [36,37].The microstructures shown in Fig. 2   A. Bhattacharya and M. Asle Zaeem on DWs such as steps, crevices and corners, and formation of these features as well as the nucleation of a 2 domain on these features result from the minimization of free energy ensuring mechanical equilibrium so that the domains always satisfy the compatibility condition [38,39].On application of tensile load to a microstructure consisting of surviving c and a 1 domains (i.e., between stages 1 and 2 in Fig. 1), the average stress gradually decreases in the entire microstructure.But some regions of microstructure may still retain a high concentration of compressive σ 11 .Since a 1 domain experiences compressive stress along its c-axis of the unit cell, the domain becomes unstable.Therefore, a 1 domain switches to a 2 domain which is stable under tension along x 2 direction as its c-axis is aligned favorably.The nucleation of a 2 domain can take place in certain localized regions where σ 11 is highly compressive.As shown in Fig. 2, the corners, crevices, and steps promote higher concentration of compressive σ 11 stress which results in nucleation of a 2 domain.The initiation of nucleation usually takes place with the formation of embryos at multiple sites and eventually these embryos grow and coalesce to form a nucleus.We note that the twin boundary between domains during nucleation is often curved because of the dominance of domain wall energy over the elastic energy, but later the situation reverses, and straight twin boundaries propagate during switching.

consist of various fine non-smooth features
In Fig. 2a, at the location marked by A in Microstructure-1, the elliptical shape of nucleus may be assumed to be the cross section of an oblate spheroid in three dimensions, while at B and C locations in Fig. 2b and Fig. 2c respectively, we observe triangular shaped nuclei.In Fig. 2c, it can be noted that the embryo of an apparently oblate shaped nucleus at location D quickly merges with the triangular nucleus at C. In the supplementary section, we have reported additional triangular and oblate spheroid shaped nuclei at E and F respectively (see Fig. S1a and  S1b).It can be noted that the nucleation of oblate spheroid shaped a 2 domain involves both a 1 →a 2 and c→a 2 switching reactions, but triangle shaped nucleation essentially involves a 1 →a 2 reaction.Similar nucleation phenomenon is also reported to occur at steplike feature 'ledges' on grain boundaries during creep fracture [40].The 'avalanche dynamics' of ferroelastic domain switching in LaAlO 3 is considered to originate from kink which is a crevice like feature on DW [41,42].
In Fig. 3, we have plotted the volume fractions of a 1 and a 2 domains and simultaneous variation of σ 11 at A, B and C locations in Fig. 2 during FDN.We observe that during nucleation, σ 11 relaxes very rapidly.After completion of FDN, the relaxation of σ 11 immediately stops or drastically slows down and this phenomenon allows identifying the completion of FDN.It is evident that the energy required for nucleation of new domain must be supplied by the elastic deformation energy.The change in elastic energy (ψ) must compensate the energy required for lattice translation (ΔG ε ) resulting in reorientation of spontaneous strains in domains as well as formation of new DWs.The change in free energy during FDN can therefore be expressed in the following way [43]: In Fig. 4a, we illustrate the change in crystallography during switching of a 2 →a 1 domains and based on this figure, we calculate the where r is the tetragonality of the unit cell.The elastic strain energy for this domain switching is then given by the following: where the elasticity moduli Substituting volume (V) and surface area (A) for these two shapes in Eq. ( 5), we determine dimensions of critical nuclei (r * for oblate spheroid and b * for triangle) as the following: where a = k t b for triangle and c = k s r for oblate spheroid shaped nuclei and the driving force for nucleation is given by (ψ − ΔG ε ).We have determined the values of ψ from Fig. 3 at these points to calculate the critical dimensions of triangular and oblate spheroidal shaped nuclei as well as the nucleation barrier (ΔG * ) and compared them with the simulations in Table 1.In these calculations, the constants k t and k s are obtained from the values reported in the second and third columns in Table 1 based on the simulations.It should be noted that the change in elastic energy due to change in applied load during FDN is negligible compared to ψ as it mostly involves the change in applied strain up to 0.5 %.
Miller considered thickness t of triangular shaped nucleus to be equal to one lattice spacing, and in this work we assumed the thickness t = 1 nm [13].The critical dimensions of triangular shaped nuclei from simulations mostly agree with those predicted by the nucleation model.We have also assessed the sensitivity of ΔG * to thickness of triangular shaped nucleus as shown in Fig. 4c, and it appears that nucleation of very thin triangular shaped domain can very easily take place.On the other hand, based on Table 1, one may speculate that the current 2D simulations have failed to capture the 3D features of nucleation of oblate spheroid shaped domains at A and F as the critical dimensions predicted by the simulations do not agree with the nucleation model.The nucleation model assumes that during FDN, the shape of nucleus remains preserved (i.e.c = k s r and k s is constant), but by Fig. 2a and Fig. S1, this assumption is clearly not valid for nucleation of spheroid domains.Furthermore, it should be noted that the nucleation of spheroid involves two simultaneous switching reactions which the current nucleation model does not consider.We attribute the complex dissipation of free energy resulting from simultaneous double switching reactions to the deviation of the present simulations from the nucleation model.Despite these deviations, the nucleation of proposed oblate spheroid shaped domain at location F (Fig. S1b) has two orders smaller ΔG * compared to that at location A (Fig. 2a), and it is comparable to ΔG * for triangle shaped nuclei.We observe that the nucleation at the neighboring C and D in Fig. 2c and similarly at G and F in Fig. S1b take place simultaneously indicating similar ΔG * values of both the nucleation events.Therefore, ΔG * for nucleation of a 2 domain is expected to be in the range of 10 − 18 − 10 − 17 J which is several orders lower than the

Table 1
The dimensions of critical nuclei, nucleation barrier and change in elastic energy as evaluated from the simulations and model calculations.calculations by Landauer [44], but 1-2 orders higher than that predicted by Miller and Weinreich [13].
The present work describes FDN in two microstructures, and the same mechanism of FDN is expected to occur in other microstructures of t'phase.Based on the work presented here and in Supplementary Materials, the described mechanism of FDN appears to be independent of DW energy.Considering twin boundary energy in the PFM, for the first time, the different shapes of nucleation are explained in the present work in a purely ferroelastic material.Some studies of domain switching in ferroelectric material have neglected the contribution of interfacial energy [15,23], but the present work emphasized its role in FDN.As elucidated in the present work for purely ferroelastic YSZ, similar strain compatibility conditions in other ferroelastic materials may also result in similar features on their domain walls where the nucleation is energetically more favorable.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Stress-strain hysteresis loop at a strain rate of 10 − 5 s − 1 in (a) Microstructure-1 and in (b) Microstructure-2 for γ = 0.2J.m − 2 .The red colored curve represents the stress-strain response under initial compressive loading on the starting microstructure which is included at the right bottom corner of each plot.The evolution of microstructure in the hysteresis loop for γ = 0.1J.m − 2 is similar to that of γ = 0.2J.m − 2 .

Fig. 2 .
Fig. 2. The nucleation of a 2 domains at DWs between c and a 1 domains at (a) crevice in Microstructure-1 for γ = 0.2J.m − 2 , (b) corner in Microstructure-2 for γ = 0.2J.m − 2 , (c) step and crevice in Microstructure-1 for γ = 0.1J.m − 2 .Below, the stress σ 11 acting along the direction x 1 is plotted for the same microstructures.The green colored contour lines representing the DWs of nucleated a 2 domains are observed in the locations marked as A, B, C, D. The blue colored arrows indicate the simultaneous formation of embryos elsewhere.

Fig. 3 .
Fig. 3.The variations of volume fractions of a 1 and a 2 domains during nucleation at (a) A, (b) B and (c) C. The variation of stress σ 11 in these locations are also included.The plots of elastic strain energy during nucleation of domain suggest major reduction of strain energy during FDN.

Fig. 4 .
Fig. 4. (a) The schematic change in crystallography during domain switching.The triangular and oblate spheroid shaped nuclei assumed are shown in schematic figures in (b).The change in free energy during nucleation of (c) triangular domain and (d) oblate spheroid shaped nuclei as function of dimensions are plotted for points C, E in (c) and A in (d).