A phase-field model for study of ferroelastic deformation behavior in yttria stabilized zirconia

A modified phase-field model coupled with linear elasticity is formulated to study the effect of ferroelastic domain switching on the deformation behavior of tetragonal-prime yttria-stabilized zirconia. Simulations are performed under uniaxial tension and compression using realistic domain microstructures. This work provides new insights into the mechanism of ferroelastic deformation by studying the evolution of domains under different loading directions and strain rates. The model predicts realistic stress-strain curves consisting of stress-plateaus with accurate domain switching behavior without setting a switching criterion based on a critical stress for nucleation of new domains. The results reveal that the origin of coercive stress lies in the balance between the energy absorbed by the domain switching process and the externally applied energy. Loading in different directions but parallel to the spontaneous strain results in different levels of coercive stress as domain switching absorbs different amounts of energy. For the same reason when the strain-rate increases, the coercive stress also increases. This work establishes the relation between the coercive stress and microstructural changes during ferroelastic domain switching.

In zirconia based ceramics, there exists three primary sources of inelastic deformation behavior [22].First is the martensitic transformation of tetragonal phase to monoclinic phase governing some interesting functional properties such as shape-memory effect and superelasticity [23][24][25][26].The second inelastic deformation behavior is dislocation-movement based plasticity, and the third is ferroelastic domain switching.Each deformation mechanism has a distinct stress-strain response.For instance, the stress-strain curve due to martensitic transformation shows serration under compression [27,28], while under tension it shows small change in slope during deformation [29].At room temperature, YSZ is brittle and does not show any plasticity, but at higher temperatures [30] 〈011〉 slip system is activated [27].A typical plastic flow behavior is observed when the material is deformed beyond the yield point.In tetragonal-prime (t′) phase of YSZ and after an almost linear stress-strain response, ferroelastic deformation mostly produces a plateau in the stress-strain curve, and the stress at this plateau is considered to be the coercive stress where the rate of ferroelastic domain switching is maximum [8,13,31,32].Some researchers have even treated this time-dependent deformation behavior due to domain switching as ferroelastic creep [33,34].Generally, it is believed that ferroelastic domain switching begins when the applied stress surpasses a critical limit, and the critical stress is the point where stress-strain curve deviates from the linear response [34][35][36].However, a recent attempt to measure the critical stress experimentally has failed, rather it has been observed that domain switching actually starts before reaching the critical stress [34].
Some researchers have considered the domain switching process as a first order phase transition consisting of three stages: the initial nucleation of domain in accordance with the classical nucleation theory, followed by bulk reorientation, and finally Ostwald ripening [37,38].Ishibashi has modified the classical Kolmogorov-Avrami model to describe the kinetics of domain switching more accurately by considering different modes of nucleation [39,40].According to the classical theory, the nucleation of domain is directly related to the critical stress [41].However, it is not clear how the nucleation takes place in a microstructure that already consists of many ferroelastic domains.
The t′ phase in YSZ is a metastable phase, and it is obtained by quenching the cubic phase when the paraelastic phase undergoes a displacive but non-martensitic phase transformation [42].This phase is considered non-transformable as it resists martensitic transformation when the domain size decreases [43,44].It has a fluorite structure with c/a ratio just greater than 1 and its ferroelastic deformation behavior is strongly influenced by its crystal structure [45,46].During ferroelastic domain switching in tꞌ phase, the interface parallel to {101} plane separating the domains migrate microscopically, while at the atomic scale the oxygen ion jumps into neighboring anion sublattice position due to shear distortion of cation sublattice [28].Since, the jump of ions requires to surpass the activation energy barrier, a critical amount of stress or heat or both must be provided externally for the jump to initiate [47].
Understanding the behavior of ferroelastic domain switching under external stress is key to the application of ferroelastic materials.Since, the materials are highly sensitive to inherent defects which are introduced during manufacturing [48], it is very difficult to study the effect of a wide range of stresses on the evolution of ferroelastic domains.Earlier, some numerical attempts have been made to simulate ferroelastic domain reorientation using the framework of finite-element method with an algorithm for domain switching [49][50][51].In those works, it was postulated that the domain switching takes place when the elastic energy surpasses a critical value which depends on the input coercive stress.This postulate is similar to the assumption of having a critical stress criterion for domain switching to start [28,35,36,47].Though, it has been possible to simulate stress-strain response in ferroelastic materials, these works have a major limitation is assuming that domains do not switch simultaneously.Improvised model has considered the appropriate kinetics of the evolving domains [51,52], but these models cannot reveal any mechanism related to the evolution of microstructures consisting of domains.
The phase-field modeling approach has become a reliable and powerful tool for computation of microstructures in materials.Since the interfaces need not to be tracked, it has been used in simulation of solidstate phase transformations [53][54][55][56].This technique has also been employed to investigate domain switching in multiferroics [57][58][59].Since ferroelectric material with perovskite structure suffers from low fracture toughness, some numerical investigations are carried out to study the effect of toughening by ferroelastic domain switching [60,61].In these works, the crystallographic variants were order-parameters and Landau-Devonshire type model were used.Phase-field model has been employed profusely to study the evolution of microstructures controlled by deformation such as martensites [53], however, these models cannot be employed to simulate domain switching as the conservation of the order parameters is not implemented [62][63][64]56,65].Consequently, there have been a very few attempts to simulate purely ferroelastic domain switching.Using multiphase-field model, Pi et al. demonstrated transitions of crystallographic domain variants under external loading [66].They also investigated the nucleation of ferroelastic twins and its effect on toughening mechanism using a phase-field model and Eshelby's theoretical approach in YSZ [67].Separately, Sun et al. used a simple phase-field model considering only two crystallographic variants and investigated the effect of domain switching on the propagation of crack [68].None of these previous works were able to concurrently simulate the stress-strain response of ferroelastic deformation and the evolution of microstructures, therefore the domain switching mechanisms were not fully understood.
In biferroic ferroelectric-ferroelastic materials, details of domain switching mechanism have been unraveled by in-situ transmission electron microscopy and piezoelectric force microscopy observations coupled with phase-field simulations [69][70][71].The understanding of the domain switching mechanism in these materials has resulted in various applications exploiting the functional properties such as design of logic gate [72], memory storage [73] and actuators [74].On the other hand, the purely ferroelastic domain switching has been mostly studied by ex-situ diffraction experiments [13,75,76] and therefore, the domain switching mechanism in purely ferroelastic material such as YSZ has not been completely understood.An understanding of the domain switching mechanism in these materials can pave the way for tailoring their compositions and microstructures to achieve desired functional properties.In the present work, we formulate a novel phase-field model to simulate the domain switching behavior of tꞌ-YSZ and its stress-strain response.In order to validate the numerical simulations, we have compared the results with a compression test of tꞌ-YSZ single crystal performed at a fixed temperature of 500 • C [77].Since tꞌ phase is not observed to experience martensitic transformation, the contribution of martensitic tetragonal to monoclinic phase transformation to the deformation behavior can be neglected.Furthermore, it has been observed that ferroelastic deformation precedes the dislocation-controlled plasticity [77], therefore, the simulated ferroelastic stress-strain response of tꞌ phase in YSZ can be directly compared with experimental observation without any ambiguity.We discuss the formulation of the model in Section 2. In Section 3, we give the details of simulation of stress-strain response of ferroelastic domain switching in YSZ.Finally, we present our concluding remarks in Section 4.

Phase-field model
In the present work, the three crystallographic variants of the t′ phase in YSZ are denoted by phase-field variables p 1 , p 2 and p 3 .These phasefield variables represent the variants which have their c-axes of the tetragonal unit cell aligned along [100], [010] and [001], respectively, as shown in Fig. 1(a).If variant p 1 is present at a point in space then p 1 = 1, and along the interface with other variants or phases p 1 smoothly varies from unity to zero.The other phase-field variables also similarly describe the other crystallographic variants and the interfaces between them.The definition of phase-field variables based on the orientation of the crystallographic variants adopted in the present work is similar to the order-parameters employed in earlier phase-field models of shapememory materials [56,65,78].But the present phase-field model for domain switching additionally applies the conservation of the phase-field variables, i.e., ∑ 3 m=1 p m = 1 without employing complicated Lagrange-multiplier based approach of multiphase-field model where the phase-field variable is actually volume fraction of the phases [66,79,80].It should be emphasized that when one domain switches into another, the total volume of t′ domains does not change, and since the density of variants are essentially the same, the conservation of mass holds during the switching.Therefore, phase-field variables representing domains must be conserved at every point of the numerical domain, i.e., ∑ 3 m=1 p m = 1.The total free energy F of the domain consists of the bulk free-energy density f bulk in homogeneous ferroelastic domains, the gradient energy density f grad in the domain-walls, and the elastic deformation energy f elast : where ε ij is the elastic strain tensor.The free energy of homogeneous f bulk in the present work is the following: A. Bhattacharya where the coefficients A, B and C are constants.The first term on the right-hand side of Eq. ( 2) penalizes the total free-energy if ∑ 3 m=1 p m ∕ = 1 or the conservation of variants is violated, and the coefficient A is set to be very high.A similar strategy was also adopted in an earlier work [81].It can be noted that second and third terms are the multiwell and multiobstacle potential functions to model two-domains and three-domains interactions [82].The gradient energy density function f grad is generally formulated in the following way: where k p1 , k p2 and k p3 are the gradient energy coefficients.According to the Ginzburg-Landau theory, the phenomenological minimization of total free energy for every phase-field variable leads to the following equations: where M pm is the mobility of p m th variant signifying the dynamics of relaxation.In Appendix-2, we have provided some details on the physical significance of the phase-field model parameters and how they can be determined.It can be noted that the phase-field mobilities may depend on temperature and strain rate.
The free energy due to elastic deformation is the following: where C ijkl is the fourth order elastic modulus tensor.The elastic strain is obtained by subtracting the stress-free transformation strain ε * ij from the total strain ε tot ij .The total strain is sum of the spatially independent homogeneous strain ε hom ij and inhomogeneous strain δε ij .
According to the theory of linear elasticity, σ ij = C ijkl ε kl where σ ij is the second order stress tensor.
are the stress-free transformation strain of the variants, and h(p) = 3p 2 − 2p 3 is an interpolation function.It is assumed that the relaxation of stress is fast enough that mechanical equilibrium prevails, therefore: In the present work, two-dimensional (2D) square domain as shown in Fig. 1(b) is considered and plane stress condition is applied.The following relations between strains and displacements are adopted in accordance with the theory of linear elasticity in 2D: , and δε 12 = 1 2 where u 1 and u 2 are the displacements along x 1 and x 2 axes respectively where σ 11 and σ 22 are the normal stresses along x 1 and x 2 respectively and σ 12 is the shear-stress.It can be noted that the model presented above considers the 3 variants in YSZ, but in Appendix-4, we have presented the generalized phase-field model considering arbitrary n number of crystallographic variants, and this model can be applied to study any ferroelastic material.

Equilibrium condition
We consider an interface between two variants and neglect the strain-energy of these two variants for the sake of simplicity.We also assume that the equilibrium between the two variants prevails so that the interface is static due to balance between interfacial energy γ and the bulk free energy of the variants.Therefore, we set p 3 = 0 and let p 1 = p so p 2 = 1 − p following the conservation rule.For convenience, we also assume that k p1 = k p2 = k p3 = k p .Therefore, Eq. ( 1) for two variant interface becomes the following: For the interface is in equilibrium, the following one dimensional relation can be derived [83]: where W = 2k p /√B is the measure of interface width as mathematically the interface is extended to infinity and also γ = k p ̅̅̅ B √ /3 (see Appendix-2 for details) [84].Thus, the parameters k pi and B can be evaluated from the values of interface width and interfacial energy or domain-wall energy.It should be noted that by definition the phase-field variables in the present model represent crystallographic orientation similar to the order parameters in the phase-field models for martensitic transformation [56,85].But additional requirement of conservation of phase-field variables for ferroelastic domain switching is ensured by choosing appropriate f bulk and f grad functionals so that the interface between two domains at equilibrium results in the phase-field profile expressed by Eq. (11).Such an interface profile has an advantage that at every point in the domain wall, ∑ 3 m=1 p m = 1 is naturally satisfied and essentially addresses the conservation in the present model particularly for symmetric profiles with k p1 = k p2 = k p3 .This originates from the fact that the three degenerate tetragonal variants are crystallographically equivalent, and therefore the width and energy of all the domain walls are the same.In the schematic Fig. 1(d) which depicts the profiles of phase-field variables of the current model in domain wall, at three points of A, B and C in the domain wall, we clearly observe that ∑ 3 m=1 p m = 1.In fact, for the phase-field profile in Eq. ( 11), the conservation of phase-field variables is maintained in all the domains and domain walls.

Numerical implementation
The partial differential equations in Eq. ( 4) and ( 9), for all the three variants, are solved using finite difference method by numerical implementation of discrete equations in C programming language.The Laplacian operator is approximated by a second-order accurate scheme while the explicit first-order scheme is used for time stepping.Eq. ( 9) for mechanical equilibrium is solved iteratively until the difference between the two successive displacements at every point is less than Δx.10 − 5 , where Δx is the spatial step in the discrete domain.In the present work, we have adopted W = Δx.In Appendix-3, the effects of W and Δx on simulation results are discussed in detail.Neumann boundary condition is applied while solving Eq. ( 4), and the mixed boundary conditions applied for Eq. ( 9) are depicted in Fig. 1(b).Tensile and compressive loading is implemented by varying boundary values of u 1 and u 2 with time in regular intervals along one side of the square domain as shown in Fig. 1(b) [56].In some compression tests, stress-relaxation tests are simulated by adjusting the displacement at the loading boundary so that the average stress in the domain does not exceed the prescribed stress limit σ R .This is done by setting the further increment of displacement at the loading boundary to zero whenever average stress is greater than σ R.
Baither et al. has performed compression tests on single crystal of YSZ at 500 • C at strain rate of 6.6 × 10 − 6 s − 1 [86].Therefore, the simulations in the present work are carried out for the aforementioned temperature and strain rate.It should be noted that such slow strain rate does not violate the assumption of mechanical equilibrium in the present model.
The experimental observation suggests that the tetragonality in yttria-zirconia system does not vary with the temperature noticeably up to 500 • C, therefore in the present work, the stress-free transformation strain at 500 • C is assumed to be the same as that of the room temperature for YSZ [66,87].On the other hand, the elastic modulus is considered to vary with temperature, and all the variants are assumed to have identical isotropic elasticity [88,89].It should be noted that the deformation behavior of tetragonal YSZ is known to be anisotropic due to a large degree of tetragonality of unit cell (around 1.45) [26,90].In contrast, tꞌ phase has a low tetragonality, and this phase is not expected to show a noticeable anisotropic deformation behavior.Consistent with the experimental observations, a coherent interface is considered in the present work [66].The details of the material properties and simulation parameters are given in Table 1.
The model size in the present work is 720 × 720 nm.We have adopted realistic numerical domains consisting of all the three variants.It has been experimentally observed that the variants are arranged in ordered colonies; in each colony, two variants of plates are arranged alternatively while one domain variant passes through to the next colony [77,91].In the present work, the simulation domain in the same way consists of two different sets of colonies as shown in Fig. 1(c).Based on the experimental observations, two microstructures, which are reported to be the main patterns, are adopted for the present study [77].In Microstructure-1, domain variants p 1 and p 2 are arranged alternatively in one colony, while the next colony is formed of variants p 3 and p 2 .Furthermore, both p 1 and p 2 are also inclined with both loading axes.In Microstructure-2, one colony consists of p 2 and p 3 , while the next colony consists of p 1 and p 3 , but p 1 and p 2 in this case lie along the either of the loading axes.Also, Microstructure-2 contains larger amount of p 3 domain as compared to Microstructure-1.The width of each domain is approximately 0.1 μm.In most of the simulations, the numerical model is subjected to compressive and tensile loads at a strain rate of 6.6 × 10 − 6 s − 1 (referred to as low strain rate in this paper).In Section 3.3 and to study the effect of strain rate, simulations with a higher compressive strain rate of 1.32 × 10 − 5 s − 1 are also completed (referred to as high strain rate in this paper).We assume that the phase-field mobility does not change noticeably in the range of strain rates applied in the present

Mechanisms of ferroelastic domain switching
We employ phase-field simulations to understand the mechanism of domain switching and its effect on stress-strain response.While studying the domain reorientation under compression and tension, we analyze the evolution of domains in order to unravel the mechanisms and identify the role of microstructures, loading direction and rate in ferroelastic deformation behavior of tetragonal prime YSZ.

Compressive loading
The morphological evolution of t′ phase and the stress-strain response under compressive loading are shown in Fig. 2 for both microstructures and loading directions.In order to understand the mechanisms of evolution of domains, we investigate the morphological evolution in Microstructure-2 under compression along direction 1.For other cases, some additional features of ferroelastic domain switching are presented in the Supplementary Material.
In Fig. 3, we plot the evolution of domains in Microstructure-2 and the corresponding evolution of strain along the loading direction, ε 11 is also included in the same figure.A. Bhattacharya and M. Asle Zaeem domain appears to be nearly completed.In order to understand the reorientation in the microstructure, we mark different regions with different colored circles and track the switching in these regions.We observe that in the region inside the red circles, domain p 1 switches into p 3 while in the region within green circles, p 1 reorients into p 2 .It appears that across the domain walls different switching response takes place in different locations.As a result of the switching, the twin boundary between p 1 and p 3 moves approximately along [110], while the same between p 1 and p 2 moves in . Since p 1 domain has c-axis along [100], on the application of compressive load along x 1 direction, it switches into p 2 and p 3 .On a closer look, we observe that at the strain of ~0.6 %, p 1 domain completely reorients within red circle while within green circle, some p 1 domain remains unswitched.Therefore, it appears that p 1 →p 3 is slightly faster than p 1 →p 2 .
We also study the evolution of ε 11 in Fig. 3, and the region with compressive ε 11 is clearly distinguished from the region with zero or small tensile strain by the dashed green lines.It is evident that most of the microstructures in the early stage of loading at 0.02 % strain have nearly zero ε 11 (which is represented by maroon color), but p 1 domains are evidently under compression which results in its reorientation .As the deformation proceeds, the region with compressive ε 11 expands.
However, we observe that within red and green circles where domain reorientation is in progress, the region with newly reoriented domains at 0.45 % strain have nearly zero ε 11 .By comparing the stress-strain curve in Fig. 2, we note that the relaxation of strain continues in the plateau of the stress-strain curve.
To gain insights into the kinetics of domain reorientation in the whole microstructure, we have plotted the variation of volume fractions of all the domains in Fig. 4(a).Furthermore, in Fig. 2, we mark the end of plateaus in the strass-strain curves with vertical dashed lines and use the same lines in Fig. 4 for the cases under compressive loading in Microstructure-2.We clearly distinguish stage-1 as the 'before stress plateau' from stage-2 as 'post stress plateau' segments in the stress-strain curves.Fig. 4(a) suggests that the domain reorientation of p 1 is complete at the end of stage-1 when the plateau is reached in the stress-strain curve.Fig. 4(b) depicts the variation of elastic energy of individual domains (which are distinguished by the phase-field variables p m ) as a function of applied strain.To simplify the calculation of volume, we have assumed a unit distance perpendicular to the 2D simulation domain.In the same plot, we have also included the homogeneous bulk free energy f bulk of the entire simulation domain consisting of three domains by using Eq. ( 2).Similar to Fig. 4(a), we also draw the vertical line marking the end of the plateau in the stress-strain curve.We observe that in the stage-1, the elastic energy of the domains does not increase substantially, but subsequently in stage-2, the elastic energy of the domains increases.The nearly constant level of elastic energy in stage-1 must be attributed to domain switching.In Fig. 3, we have observed the relaxation of strain in the switched region.The relaxation of strain indicates that the strain energy should not increase during domain reorientation despite the application of load.
The variation of volume fractions of domain variants shown in Fig.(a) is not entirely consistent with the observed domain-switching behavior explained in Figs. 2 and 3. Though the volume fraction of the domain p 1 steadily decreases to zero and that of p 3 increases, the volume fraction of p 2 decreases despite the fact that p 1 →p 2 is clearly identified in the region within the green circle.On a closer look in Fig.
(a), we observe that in the early stage of loading, the volume fraction of p 3 decreases and that of p 2 increases in exactly opposite manner.This is attributed to additional domain switching of p 3 →p 2 in the initial stage, as shown in the region within yellow circle up to the strain of 0.25 % in Fig. 3. Furthermore, in stage-2 after the strain of 0.8 %, we similarly observe the change in volume fractions of p 2 and p 3 in Fig. 4(a) which is marked by black arrows.In this case, the reverse switching f p 2 →p 3 takes place.This additional domain switching is also captured in Fig. 4(b) where the reverse switching results in dip in elastic energy as marked by the black arrow indicating that this additional domain switching also absorbs the elastic energy.In Fig. 3 at the strain level of 0.25 %, we observe the final stage of the additional domain switching of p 3 →p 2 where the small bridge of p 3 domain on p 2 disappears.The bridge of p 3 domain in the initial stage of loading has ε 11 ∼ 0 but the same location develops compressive ε 11 after the strain of 0.45 %.It appears that the compressive ε 11 does not develop in the p 3 domain in the initial stage of loading which force p 1 and p 2 domains to bear the applied compressive load.When the compressive load increases further, p 1 and p 2 domains, due to smaller volume fractions, can no longer support the applied load.Therefore, p 3 switches into p 2 .In the later stage when p 3 domains also develop compressive ε 11 , some p 2 reorients back to p 3 to relax ε 11 strain as the relaxation of ε 11 is evident after the strain of 0.45 % in the switched p 3 domain (within yellow circles in Fig. 3).The relaxation of strain results in dip in elastic energy as shown in Fig. 4(b).
Under compression, we have identified the usual p i →p j + p k switching, and we have also discussed additional switching mechanism which similarly absorbs the externally supplied energy.However, under compression, p i →p j + p k switching does not always take place.For instance, in case of loading in direction 1 in Microstructure-1, we have identified the mechanism of "partial switching" where switching between two domains takes place while the third domain does not A. Bhattacharya and M. Asle Zaeem participate in any switching and its volume fraction does not significantly change with applied strain, as shown in Figs.S1 and S4.The switching mechanism widely varies with the loading directions and microstructures.

Tensile loading
In Fig. 5, we have plotted the stress-strain response of both the microstructures under both the tensile loading directions.We elucidate the mechanism of domain switching in Fig. 6 in Microstructure-2 under tensile loading direction-2.In the same figure, we plot the strain in the loading direction ε 22 and distinguish the regions under compression and tension in the similar way.We mark two regions bounded by red and green circles to track the switching in the similar way and observe that within the red circle p 3 →p 2 and within the green circle p 1 →p 2 take place.Since p 1 and p 3 domains have their c-axes lying perpendicular to loading direction x 2 , both domains under tension switch into p 2 domain which has its c-axis along [010].Furthermore, we note that the movement of the domain wall is exactly opposite to what observed in Fig. 3 during compression in the same microstructure.
We similarly study the evolution of strain ε 22 in the direction of loading to relate the evolution of microstructure to loading.We observe that in the region bounded by the red circle within p 3 domain and near it in p 1 domain, small tensile ε 22 results in domain switching.In comparison, the regions with p 2 domains have in fact slightly compressive ε 22 until the applied strain of 0.7 %.We also observe that p 3 →p 2 is dominant over p 1 →p 2 , and after strain level of 1 % significant amount of p 3 →p 2 has taken place whereas p 1 →p 2 is not quite visible.We also observe that as the applied strain increases, tensile state of ε 22 strain significantly increases particularly in p 1 domains, but the region with p 2 domains have mostly remained around ε 22 ∼ 0 above the strain level of 0.7 %.This behavior is attributed to the relaxation of strain after domain reorientation as explained in the earlier subsection for compression cases.
In Fig. 7(a), we plot the volume fractions of all the three variants as a function of applied strain.In this plot, the end of plateau in the stressstrain curve in Fig. 5 is also marked by the dashed vertical line.It is evident that globally in Microstructure-2 under tensile loading, p 3 reorients into p 2 which results in the reduction of the volume fraction of p 3 , but p 1 initially slightly increases up to strain of 0.8 % and then slowly decreases.It should also be noted that according to Fig. 7(a), the switching is far from completion till the strain level of 1.2 %.Till the end of the plateau, though p 3 →p 2 takes place steadily but globally the volume fraction of p 1 decreases only in stage-2.This is consistent with our observation of domain switching in the region within the green circle in Fig. 6.The same figure also shows that in stage 2, the region with p 1 domains have high tensile ε 22 which promotes its switching behavior.
Apparently, the activation of switching of p 1 domain requires much higher load which is achieved in stage 2. In Fig. 7(b), we observe that the elastic energy of all the domains remain nearly constant in stage-1, but in stage-2 they all increase rapidly.
In general, the switching response p i + p j ↔ p k under tension as well as compression actually consists of two partial switching responses p i ↔ p k and p j ↔ p k and domain wall moves as a result of domain reorientation across it.These two partial switching responses may have fast identical kinetics resulting in near completion of switching response.The evolution of domains under compressive loading in both the directions in Microstructure-2 and in Microstructure-1 in loading direction 2 can be considered to belong this category (see Supplementary Material).In the other cases of tensile loading, these partial switching responses are slow and have different kinetics resulting in domain switching response which is far from completion at the end of the plateau in the stress-strain curve.The evolution of domains discussed here for all the cases are therefore essentially controlled by the relative kinetics of these two partial switching responses.Physically, different kinetics of partial switching response must be related to friction [92].It is expected that the stress-strain response will also depend on the kinetics of these two partial switching responses.

Domain switching and the formation of plateau in stress-strain curve
The homogeneous bulk free energy f bulk of the entire simulation domain decreases under all the compressive and tensile loading in both the microstructures as shown in Figs.4(b) and 7(b), signifying the fact that domain reorientation indeed results in minimization of energy.Furthermore, analyzing the result in Fig. 7(a) suggests that the switching process is slow initially due to very low driving force for switching but subsequently it accelerates as the driving force for switching increases with the applied strain, and finally it stabilizes when the domain reorientation ceases.In Fig. 4(a), the additional domain switching complicates the behavior of time dependent evolution of domains, but at the end of domain switching, the volume fractions of domains indeed do not change.It can be noted that previous works on domain switching have assumed that onset of reorientation takes place when the elastic energy surpasses a critical limit which is determined by the critical stress [35,36,49,50].In Fig. 7, we have observed that the switching of p 1 domains starts in stage-2 when the strain in the loading direction attains sufficiently high value.This observation might be consistent with the theory of critical stress/energy [34,36].It can be noted that the starting microstructure consists of all the three domains and no exclusive Fig. 5.The evolution of domains in (a) Microstructure-1 and (b) Microstructure-2 under tension and the corresponding stress-strain curves at strain rate 6.6 × 10 − 6 s − 1 .Orange-colored arrows mark the stress plateaus.Dashed vertical lines represent the end of plateau in the stress-strain curve for tensile loading conditions 1 and 2, respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)nucleation process is implemented in the current phase-field model, therefore this model can be construed as the second category of domain switching process of Ishibashi-Kolmogorof-Avrami's model which does not require any nucleation step [39].After the application of load without the need of nucleation of p 1 or p 2 domain in tension and p 1 or p 2 and p 3 domains in compression, the domain switching starts almost immediately, however it needs to overcome small resistance provided by the domain wall.The fact that the domain reorientation begins near the start of the loading is also corroborated by an experimental observation [34].Furthermore, we have not applied any switching criteria exclusively and yet the domain switching is simulated reasonably, and the stress-strain curves also show a stress plateau consistent with the experimental observations [13,86].It appears that the migration of domain wall requires to overcome a critical level of stress/energy rather than the nucleation of new ferroelastic domain, and domain switching can proceed without any further nucleation.
The stress-strain curves simulated in Figs. 2 and 5 have shown initial linearity followed by well-defined stress plateaus representing coercive stress which is a characteristic of ferroelastic deformation.Similar to the experiment, we have also simulated stress-relaxation test as shown in Fig. 2. In experiments, the selected values of σ R appear to be lower than the stress level of the plateau in the stress-strain curve [86].It can be observed in Fig. 2 that stress is relaxed whenever the loading is stopped, however the drop of the stress is different in different loading conditions and microstructures.Apparently the internal stress varies with microstructures and loading [93].Furthermore, the domain switching still proceeds at a nearly constant or a very close range of stresses and this is similar to creep-like behavior of ferroelastic domain switching [34,94].In Figs.4(b) and 7(b), we also plot the elastic energy of the alldomain variants averaged by their respective volume fractions and it is identified as the average elastic energy in these plots.Considering all the loading conditions, in tension and compression, and microstructures, the average elastic energy very slowly increases or remains nearly constant in stage 1, but rapidly increases with time in stage 2. The constancy of average elastic energy in stage 1 indicates that the energy absorbed by the domain reorientation per unit time is the same as that supplied externally.In the earlier experimental works, it was proposed that the reason for the plateau formation is the absorption of the externally supplied strain energy by the domain switching process and the present work for the first time confirms it [13,44,95].In the initial stage of loading, the major portion of the energy supplied externally remains unconsumed by the domain switching process and it results mainly in elastic deformation of the material.This stage results in a linear stress-strain curve.Subsequently, the extent of domain reorientation increases with time due to either acceleration or steady rate of domain switching, and more energy is absorbed by the process of domain reorientation.When the energy supplied by the external loading is almost fully consumed by the switching process, the stress-strain curve displays a plateau in the stress-strain curve over a very narrow range of stress, and average elastic energy does not increase significantly at this stage.The competition between externally supplied elastic deformation and its absorption by domain reorientation is captured in stage-1 in Figs.4(b) and 7(b), and it is also characterized by the maximum domain switching rate [32].Subsequently, the exhaustion of orientable domains or small remnant domain limits the absorption of additional elastic energy due to externally applied load, the linear elastic deformation returns in stage 2 and elastic energy rapidly increases.
The final configurations of the microstructures and the migration of the domain-wall simulated in the present work agree well with the experimental observations and an earlier phase-field simulations of domain switching [9,10,67,78,92].Though, the present work simulates the realistic stress-strain response of ferroelastic domain switching as opposed to arbitrarily applied constant load in the earlier work [66].Additionally, it also for the first time unravels details of the mechanism of domain switching identifying the role of microstructure and loading in ferroelastic deformation behavior.Though, we have understood the relation between microstructure and the plateau in the stress-strain curve qualitatively, a quantitative relationship between these two have not been investigated yet.

The relation between coercive stress and evolving ferroelastic microstructures
Ferroelastic domain switching is primarily characterized by the coercive stress which is a plateau in the stress-strain curve.We estimate the coercive stress ranges for all the stress-strain curves by adopting the procedure described in an earlier work and the method is explained in Appendix-1 [96].In Table 2, we have summarized the calculated coercive stresses for all the simulations.It can be immediately observed that coercive stress strongly depends on the loading type and direction as well as the microstructures.It can be noticed that the coercive stress obtained under loading-1 is consistently higher than that obtained for loading-2 for the same microstructure and the sense of load.
It is already realized that the plateau in stress-strain curve is related to the evolution of domains during switching.We note that a microstructure consisting of domains can be quantitatively described by the volume fractions of the three domains.We construct a 3D plot with the three axes representing the volume fractions of the three domains of p 1 , p 2 and p 3 , and a particular microstructure therefore represents a point in p 1 -p 2 -p 3 space.The evolution of microstructure is thus represented by trajectory of these points in p 1 -p 2 -p 3 space.In Fig. 8, we plot the evolution of microstructures for some simulations performed at low strain rate in p 1 -p 2 -p 3 space.For convenience, we have also included the projection of these trajectories on p 1 -p 2 plane in the same figure.Ideally, all the trajectories should start from two points corresponding to the two microstructures.However, due to the initial transient state during loading, the starting microstructures differ in all the cases.Fig. 8 clearly shows that the direction of trajectories in p 1 -p 2 -p 3 space under the same loading direction (loading-2) and state of stress (compression) remains nearly parallel.If the state of stress is reversed (i.e., under tension), the direction of trajectory also changes irrespective of the microstructure.Nevertheless, this plot allows us to measure the microstructural changes quantitatively.We calculate the length of these trajectories in p 1 -p 2 -p 3 space between the initial and arbitrarily defined final state by the following relation: where V pi is the volume fraction of p i domain variant.In Table 2, we have also included the calculated values of L for all the simulations.In these calculations, the initial state is defined by the start of the loading whereas the final state represents the microstructure corresponding to a point just above the plateau in the stress-strain curve.Table 2 clearly shows that L has consistently lower values under loading-1 in comparison to loading-2.We also observe that value of L calculated is fairly consistent with the observation of evolving microstructure discussed in Section 3.1 and in the supplementary material.In Microstructure-1 under tensile loading in direction 2, Fig. S7 and Fig. S9 show that both p 1 →p 2 and p 3 →p 2 partial switching reponses proceed steadily and its L value is the highest.On the other hand, additional domain switching ultimately result in a lower value of L in Microstructure-2 under

Table 2
The coercive stress-range and the parameter L evaluated for all the conditions.compressive loading in direction 1.Small values of L indicates small net change in microstructure and less absorption of energy leading to domain reorientation.Therefore, the unconsumed balance of externally supplied energy will result in a larger elastic strain energy in the numerical domain and hence a higher coercive stress.The parameter L is certainly related to the friction of the domain-wall motion [28,92].Larger friction should result in smaller values of parameter L. Different coercive stress values in two different loading directions parallel to stress-free strain suggest that the friction to domain motion is different even in the same microstructure.This is attributed to the difference in absorption of the elastic energy due to difference in orientations.
Lankford reported that tꞌ single crystal of YSZ showed different stress-strain responses when loaded along crystallographic direction [100] and [123] directions.Crystal loaded along crystallographic direction [100] showed flat stress-plateau whereas loading along [123] resulted in moderate change in the slope of stress-strain curve at much higher stress level than the stress-plateau formed when loaded in crystallographic direction [100] direction [28].Lankford considered the fact that the stress-plateau at crystallographic direction [100] represented a critical stress limit for domain switching.Therefore, the inclination of c-axis of the same domain with respect to [123] loading direction must require higher applied stress level for domain switching because of the inclination.Though, the effect of orientation of domains with respect to the loading direction is already recognised in the present work, but it still offers an alternative explanation of the effect of loading direction in terms of the microstructural change parameter L. Due to inclination of the c-axis in crystal oriented along [123], less elastic strain energy will be available for domain reorientation and therefore the microstructure should absorb less energy for domain switching resulting in smaller value of L. As mentioned, the plateau in the stress-strain curve appears when the energy absorbed by the process of domain switching and that supplied externally per unit time are comparable, but in order to reach this stage, sufficient extent of domain reorientation must take place and it is achieved at higher externally supplied strain energy resulting in increase in coercive stress.

The relation between evolving ferroelastic microstructures and strainrate sensitivity
We also discuss the strain-rate sensitivity of ferroelastic deformation.It was reported that the stress level of the plateau in the stress-strain curve increases with strain rate [28].Fig. 9(a) and (b) show the effect of strain rate on the stress-strain curve under compressive loading in both the microstructures for both the loading conditions and coercive stresses can be observed to increase when the strain rate has increased from 6.6 × 10 − 6 s − 1 to 1.32 × 10 − 5 s − 1 indicating the consistency with the experimental observation of ferroelastic deformation.In Fig. 9(c) and (d), we plot the variation of volume fractions of the domain-variants as functions of applied strain at higher strain rate of 1.32 × 10 − 5 s − 1 .Comparing Fig. 4(a), Fig. S4, Fig. S9, Fig. 9(c) and Fig. 9(d), it can be observed that the kinetics of evolution of domains at high strain rate is faster than that at low strain rate, but it does not reveal any remarkable A. Bhattacharya and M. Asle Zaeem change in microstructural evolution nor any difference in the mechanism of domain-evolution.Faster kinetics at a high strain rate is attributed to the increase in strain energy when strain rate is increased.
According to the theory of ferroelastic creep, the strain rate is directly proportional to the activation energy [34,94,98].When the applied strain rate increases, the additional activation energy must be provided by the application of stress and therefore, coercive stress increases.In the same Table 2, the microstructural change parameter L for high strain rate is also included.The values of L for higher strain rate are in general higher than that for low strain rate.This fact is also corroborated by Figs.4(a), S4, 9(c) and (d) where the final volume fractions of domains differ for two different strain rates.It indicates that the partial switching responses proceed faster at a higher strain rate.However, it may lead to an apparent contradiction as both the coercive stresses and the parameter L in Table 2 are higher for deformation at higher strain rate.However, it should be noted that the amount of deformation energy per unit time increases with strain rate therefore the kinetics of domain-switching at higher strain rate will be obviously higher.The comparison will be meaningful at the equal strain level which also indicates the equality of driving-force for switching.Since a higher strain rate requires less time to reach the same level of strain, the change in microstructure is actually less due to less transformation time available as we assume that the phase-field mobility does not change in the range of strain rates applied in the present work.Therefore, a higher coercive stress at a higher strain rate is again attributed to the low value of L. Finally, we summarize the relation between coercive stress and L in Fig. 10.We have considered average values of coercive stress ranges and plotted them as function of L in Fig. 10 for the deformations carried out in both the microstructures at both the strain rates.It is evident that the average coercive stress decreases with the microstructural change parameters L for all the cases discussed in this work.The present work therefore not only describes the mechanism of domains switching in detail, but also it identifies the parameters affecting the ferroelastic deformation behavior.

Conclusion
We have presented a novel phase-field model for simulating domainswitching and ferroelastic deformation.The model accurately predicts the switching behavior in response to external loading.The volume fractions of domains continuously vary with time under external load.The switching of the domains starts almost immediately after the application of load.The present work reveals that the domain reorientation requires the presence of favorable distribution of strain, but after switching, the strain is immediately relaxed.It has also been observed that orientation of domains in the microstructure plays a significant role in the kinetics of the domain-switching.Our simulations also capture the formation of plateaus in the stress-strain curves.During loading, a stage is attained when a substantial fraction of externally supplied energy is consumed [65]by the domain switching process.At this stage, strain increases while stress does not vary and the plateau in stress-strain curve forms resulting in a significant amount of domain reorientation in the microstructure.Till the completion of plateaus in the stress-strain curve, the average elastic energy of the microstructure does not increase due to domain-reorientation, but subsequently elastic deformations of the surviving domains increase the elastic energy.We have presented a scheme for quantification of the microstructural changes during the domain-reorientation in YSZ in terms of volume fractions of domains and established its relationship with the coercive stress.Our simulations suggest that coercive stress is inversely proportional to the net microstructural change irrespective of applied strain-rate and the loading directions.Though the model is developed for YSZ, it is certainly capable of simulating ferroelasticity in other materials as well, and the generalized phase-field model to simulate domain switching in any purely ferroelastic material is presented in Appendix-4.
In the present work, loading in two mutually perpendicular directions is considered and coercive stress is found to depend upon loading directions despite the fact that stress-free strain is parallel to loading direction.For better understanding of the influence of loading direction, arbitrary loading angles should be explored when the microstructures are oriented at different angles with respect to loading axis in the real polydomain single crystal.Since the starting microstructure consists of all the domains, the effect of domain nucleation cannot be assessed.But domain-nucleation can be studied with the same model if an appropriate microstructure is chosen.In the present work, we have studied only two different strain rates, more strain rates should be explored for better insight.At high strain rates, the temperature of the Fig. 10.The variation of coercive stress with the microstructural change parameter L for all the deformations in both the microstructures at both the strain rates.The red and blue colored plots represent slow strain rate compressive and tensile deformations, respectively, while the green-colored plot illustrate high strain rate compressive deformation.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)material increases, and this may activate the slip system for plastic deformation.Though the current model is not capable of capturing the effect of high strain rate, nor it is coupled with the plastic deformation, but in future it will be certainly an interesting numerical study with appropriate model if not possible in the real experiments to enhance our insight into the ferroelastic deformation behavior [65].

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.
where the driver of the domain wall motion is C ijkl ε kl ], and the elasticity tensor is identical in all the three variants.Considering 1D form of Eq. (A1) and multiplying dp/dx to both sides of this equation, we obtain: where v is the steady velocity of the domain wall.At equilibrium condition, the left-hand side and the last term of the right-hand side of Eq. (A2) become zero, and we obtain: By solving Eq. (A3), one can derive the equilibrium phase-field profile p(x) in Eq. (11).On the other hand, using the 1D form of Eq. ( 10) and considering the fact that the total free energy of a two-domain system at equilibrium is nothing but the interfacial or domain wall energy (γ), we obtain the following relation for equilibrium condition: Assuming that the deviation from equilibrium is very small so that the phase-field profile of the moving domain wall does not differ noticeably from the equilibrium phase-field profile, and by integrating Eq. (A2) from − ∝ to + ∝, we obtain the following relation for the mobility: where v = MΔf elast and M is the mobility of domain wall migrating in diffusionless manner [99].Therefore, the phase-field mobility M p can be directly evaluated from the domain wall mobility by Eq. (A5).However, due to lack of knowledge of this mobility, we have run several simulations with different M p values and the corresponding stress-strain responses are reported in Fig. A2.For a low value of M p , the kinetics of domain switching is slow, and the stress-strain curve shows small inelastic deformation.But, for a higher value of M p = 1.8 × 10 − 9 m 3 J − 1 s − 1 , a plateau forms in the stress-strain curve.In present simulations, we have applied a value of M p , beyond which the model does not converge.A4 all domains tend to switch to p 3 which contradict the experimental observations where under compression two domains should survive [86].Such unrealistic domain switching behavior also results in unrealistic stress-strain response.But for a higher value of C, the unrealistic domain switching behavior is restricted.Floquet reported that in BaTiO 3 , the 90 • ferroelastic domain wall is 4-6 nm wide [97], and a later study suggested that the width of domain walls varies between 1.5-3 nm [100].In the present study, we have chosen W=1.5 and 3 nm.By Eq. ( 11), the equilibrium phase-field profile is extended to infinity, but W = 2k p /√B is a finite measure of interface width or domain wall thickness.One can choose the spatial discretization step Δx same as W or W can be integral multiple of Δx.We have simulated ferroelastic domain switching under compression along x 2 direction (see Fig. 1(b)) for three conditions: (i) W = Δx = 3 nm, (ii) W = Δx = 1.5 nm, and (iii) W = 3 nm and Δx = 1.5 nm.In these simulations, the phase-field mobility (M p ) values are calculated by using the relation in Eq. (A5), and we have summarized the results in Fig. A5.The stress-strain response and the evolution of domains do not noticeably depend on the choice of W and Δx, particularly in the range of 1.5-3 nm.Therefore, in the present work all the simulations are performed using W = Δx = 3 nm.Using these larger values for W and Δx significantly expedites the simulations of the 720 nm × 720 nm study domain which is chosen based on the experimental observation [86].

Fig. 1 .
Fig. 1.(a) Crystallographic variants of t′-YSZ.(b) The schematic of boundary conditions for displacements and loading conditions.(c)The initial numerical domains consisting of the three phase-field variables p 1 , p 2 and p 3 and the microstructures created based on the experimental observation [77].Colony boundaries are marked by red colored lines.(d) The symmetric phase-field profiles of p 1 & p 2 satisfy p 1 + p 2 = 1 at A, B and C. In the present model, all along the domain wall p 1 +p 2 +p 3 = 1 ensuring the conservation of phase-field variables.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) kp 1 = kp 2 = kp 3 = kp 3.0 × 10 − 5 J 1/2 /m 1/2 Coeff.A 3 × 10 10 J/m 3 Coeff.B 4 × 10 8 J/m 3 Coeff.C 3 × 10 10 J/m 3 Phase-field mobility Mp 1 = Mp 2 = Mp 3 = Mp 1.8 × 10 − 9 m 3 /Js A. Bhattacharya and M. Asle Zaeem work.

Fig. 2 .Fig. 3 .
Fig. 2. The evolution of domains under compression and corresponding stress-strain curve at strain rate of 6.6 × 10 − 6 s − 1 in (a) Microstructure-1 and (b) Microstructure-2.Stress plateaus are marked by orange-colored arrows.Dashed vertical lines represent the end of plateaus in the stress-strain curves under compression for loading conditions 1 and 2, respectively.The dotted horizontal black lines represent the curves with stress-relaxation in respective loading directions.The simulation is compared with experiment [77].(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. (a) The variation of volume fractions of the three variants with applied compressive strain at the strain rate of 6.6 × 10 − 6 s − 1; (b) The corresponding variation of elastic energy of all the three domain variants and their average, and the bulk free energy.Dashed green vertical line represents the end of plateaus in the stressstrain curves in Fig. 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.The evolution of domains during tensile loading (top) and the corresponding distributions of ε 22 is plotted at the bottom.The arrows on the top row indicate the direction of migration of the domain wall.In the bottom row, the green-colored dashed line is the boundary of region having compressive ε 22 , and the blue colored line indicates the domain-wall.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. (a) The variation of volume fractions of the three variants with applied strain during tension at strain rate of 6.6 × 10 − 6 s − 1 .The vertical blue dashed lines represent the end of plateau in the stress-strain curves in Fig. 5.(b) The variation of elastic energy in three variants, average elastic energy and bulk free-energy of the microstructure with strain during tension.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Microstructure-2, Loading-1 & Compression 770-790 0.21 Microstructure-2, Loading-2 & Compression 380-410 0.35 Microstructure-2, Loading-1 & Tension 650-810 0.5 Microstructure-2, Loading-2 & Tension 640-780 0.53 High strain rate: 1.32 × 10 − 5 s − 1

Fig. 8 .
Fig. 8. Trajectories of evolution of domains during switching by low strain-rate deformation in p 1 -p 2 -p 3 volume fraction space represented by points.The projection of these trajectories onto p 1 -p 2 is represented by continuous lines of same colors.The arrows of the same colors indicate the direction of evolution in p 1 -p 2 plane.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Stress-strain curves at low (6.6 × 10 − 6 s − 1 ) and high (1.32 × 10 − 5 s − 1 ) strain rates in (a) Microstructure-1, and (b) Microstructure-2.The stress plateaus are marked by orange-colored arrows.Corresponding variation of volume fractions of the three variants with strain during compression at high (1.32 × 10 − 5 s − 1 ) strain rate in (c) Microstructure-1, and (d) Microstructure-2.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. A2 .
Fig. A2.The effect of mobility M p on stress-strain response in both the microstructures under compression along x 2 .A relativly larger mobility value results in formation of a plateau in the stress-strain curve.Free energy function parameters A, B, C and k p : The parameter A ensures the conservation of the phase-field variables.We ran several simulations with different values of A, and Fig. A3 summarizes the effect of A on the conservation of phase-field variables for three simulations.For a low value of A = 10 8 J.m − 3 , the summation of the phase-field variables p 1 + p 2 + p 3 significantly differ from unity in the domain wall.For A = 3 × 10 10 J.m − 3 , the conservation is perfectly maintained along the domain wall and this value is therefore used for the present simulation work.The parameter B and k p are related to domain wall width and domain wall energy by the following relations: W = 2k p /√B and γ = k p ̅̅̅ B √ /3.From the input of W and γ values, these two parameters are calculated.The parameter C defines the activation barrier for domain wall motion particularly at three domain junctions.For very low values of C, the extent of domain switching increases and as shown in Fig.A4all domains tend to switch to p 3 which contradict the experimental observations where under compression two domains should survive[86].Such unrealistic domain switching behavior also results in unrealistic stress-strain response.But for a higher value of C, the unrealistic domain switching behavior is restricted.

Fig. A3 .
Fig. A3.The effect of the free-energy parameter A (in Eq. (2)) on the conservation of the phase-field variables.Here, Microstructure-1 is subjected to compression along x 2 for different values of A.

Fig. A4 .
Fig. A4.The effect of the free-energy parameter C (in Eq. (2)) on the domain switching behavior.Microstructure-1 is subjected to compression along x 2 for different values of C. Appendix 3: The convergence behavior of the present model

Fig. A5 .
Fig. A5.The stress-strain curves in (a) and (b) and evolution of volume fractions of p 1 , p 2 and p 3 in (c) and (d) in both the microstructures for different values of Δx and W.

Table 1
Materials properties and model parameters.