Multiscale modeling of ferroelectrics with stochastic grain size distribution

Macroscopic properties of ferroelectrics are controlled by processes on the microscale, in particular the switching of crystal unit cells and the movement of domain walls, respectively. Besides these microscopic levels, the grains of a polycrystalline material constitute the mesoscopic scale. Interactions of grains with statistically distributed orientations, as a consequence of mechanical and electrostatic mismatch, give rise to for example, residual stress which in turn affects domain switching. A multiscale modeling thus has to incorporate at least three interacting scales. In this context, the condensed method has recently been elaborated as an efficient tool with low computational cost and effort of implementation. It is extended toward statistical distributions of grain sizes in a representative material volume element and amended with regard to the modeling of domain evolution. Each of the few parameters of the constitutive approach has a unique physical meaning and is adapted to available experimental values of macroscopic quantities of barium titanate taken from various sources.


Introduction
The class of so-called smart materials nowadays ranges from magnetostrictives and shape memory alloys to electro-or magnetoactive polymers and ferroelectrics. The latter belonging to the wider class of piezoelectrics and being characterized by the capability of switching spontaneous polarization, still play a crucial role herein, combining large stiffness and actuation power with distinguished responsivity. The nonlinear constitutive behavior of ferroelectrics, often illustrated by so-called butterfly hysteresis loops or those of electric displacement versus electric field, has extensively been investigated in the past both theoretically (Chen and Lynch, 1998;Cocks and McMeeking, 1999;Huber and Fleck, 2001;Huber et al., 1999;Hwang et al., 1995Kessler and Balke, 2001) and experimentally (Fo¨rderreuther, 2003;Franzbach et al., 2014;Huan et al., 2014;Mauck and Lynch, 2003;Wang and Li, 2020) to name only a few of the outstanding works. To validate constitutive models or identify model parameters, the maximum of strain and polarization versus electric field and the respective remanent values are essentially taken as a basis.
Appropriate constitutive models of ferroelectrics have to account for processes on different scales. In particular microphysically motivated models, which largely abstain from phenomenological parameters recorded on an engineering scale, have to account for the switching of crystal unit cells and the resulting movement of domain walls. Since only the latter can be observed under an optical microscope, however, both processes are assigned to different scales, that is, microand mesoscale. Another feature of the mesoscale in a polycrystalline material is its grain structure. Exhibiting statistically distributed orientations of crystalline axes and exposed to electromechanical loading, a mismatch of strain and polarization is observed at the grain boundaries. The resulting interaction gives rise to i. a. residual stress which, in turn, has an impact on domain switching. On the macroscopic scale, a representative volume element (RVE), containing a sufficiently large number of grains, yields the quantities being of interest for engineering analyses.
The first contributions to homogenization techniques probably trace back to the works of Voigt (1889) and Reuss (1929). More advanced analytical homogenization methods are, for example, the Mori-Tanaka method (Mori and Tanaka, 1973), the differential scheme (Norris, 1985), the self-consistent method (Hill, 1965;Kro¨ner, 1958) or Hashin-Shtrikman type formulations Shtrikman, 1962a, 1962b). However, these analytical approaches have their focus on the calculation of effective moduli of nonhomogeneous materials, rather than providing a coupled multiscale analysis. Numerical homogenization approaches are, for example, the FE2-method, introduced by Smit et al. (1998) and extended toward multiphysical material behavior, for example, by Schro¨der and Keip (2012), Labusch et al. (2014Labusch et al. ( , 2019 or Uetsuji et al. (2008Uetsuji et al. ( , 2012Uetsuji et al. ( , 2019, or FE-FFT based methods (Kochmann et al., 2016;Moulinec and Suquet, 1998). Semi-analytical methods have recently been applied to multiscale problems, for example, by Wulfinghoff et al. (2018) or Jaworek et al. (2020), where an analytical homogenization procedure for the microscopic boundary value problem and thus for each Gauss point is considered within a FE framework. Figure 1 illustrates multiscale modeling features of ferroic functional materials, where the macroscopic mixed boundary value problem is depicted on the left hand side, whereas the blue box on the right comprises aspects of different scales. An RVE represents a macroscopic material point and can involve grains of different compositions and for example, microcracks. Ferroelectric grains exhibit mesostructures with for example, 908 and 1808 domain walls in a tetragonal phase. Unit cells, in turn, constitute the domains, whereat different types may basically coexist in one grain, for example, for morphotropic compositions of PZT ceramics.
In the recent past, all these aspects have been included in modeling approaches of ferroelectric, ferromagnetic, and multiferroic materials, based on the so-called condensed method (CM) (Behlen et al., 2021;Ricoeur, 2015, 2016;Ricoeur and Lange, 2019;Uckermann et al., 2018;Warkentin and Ricoeur, 2020). This semi-analytical approach comprises homogenization within a multiphysical framework and scale-bridging interactions in a simple but robust and thermodynamically and electromechanically consistent manner. Its implementation is straightforward and results of RVEs are obtained with comparatively low computational cost. In comparison to most of the above mentioned multiscale approaches, the CM is able to calculate macroscopic as well as microscopic quantities without any kind of discretization scheme.
This work pursues two objectives. The model of tetragonal ferroelectrics developed in the context of the CM is extended at two points. Rhombohedral and ferromagnetic systems, as indicated in Figure 1, are not considered here, however, could straightforwardly be included in the extensions. In the averaging approach grains have hitherto been considered as equally sized, which is now replaced by an assumption of Gaussian distribution. The stochastic nature of the problem is thus incorporated twofold now, that is, in terms of grain orientation and relative size. The absolute sizes distinguishing fine and coarse grained behaviors in terms of different hierarchical domain levels, effects on lattice parameters (Huan et al., 2014;Li and Wang, 2017) or influences on material properties, see for example, Tan et al. (2015), are disregarded at this point. Another extension of the model introduces a lower limit for magnitudes of internal variables representing volume fractions of domains in a grain. Physically, it ensures that moving domain walls are not allowed to vanish due to external loading. Introducing an additional parameter on the one hand side, it replaces two tuning parameters on the other, which have hitherto been included in the constitutive equations to reduce switching strain and polarization.
The second objective of the work is to adapt essential parameters of the model to experimental findings of barium titanate (BT) as an exemplary material and to demonstrate the quantitative appropriateness of the modeling approach. In this context, experiments turn out to yield a considerable range of relevant quantities which is recorded compiling available hystereses of strain and electric displacement versus electric field. Predicted remanent and maximum values of strain and electric displacement are compared to the respective experimental ranges. The studies target accuracy just as computational costs, partly constituting competing factors. Residual stress is finally investigated due to its relevance for strength and reliability of the brittle ceramics.

Constitutive framework of a ferroelectric grain
The thermodynamic potential of the considered nonlinear ferroelectric material behavior reads, see for  Arlt (1990b), Egelkamp and Reimer (1990), and Scholehwar (2010). example, Wingen and Ricoeur (2019) or Ricoeur and Lange (2019): In equation (1) C ijkl , e lij , and k il describe the elastic, piezoelectric and dielectric properties and e irr kl and P irr i denote irreversible contributions due to domain wall motion. The independent variables within the constitutive framework are e kl and E l . The partial derivatives of equation (1), with respect to the independent variables, lead to the constitutive equations, see for example, Lange and Ricoeur (2015), where s ij and D i denote the associated variables mechanical stress and electric displacement. In equations (2) and (3) the material properties of a grain have been assumed constant in incremental changes of state and are determined as weighted averages, see for example, Avakian and Ricoeur (2016) or Huber et al. (1999), where domains with homogeneous material properties are assumed and p is associated with the number of domain variants per grain. In equation ( represent the elastic, piezoelectric, and dielectric coefficients of the domain n. Ferroelectric materials with tetragonal unit cells exhibit a domain structure with 908 and 1808 domain walls, which is depicted in Figure 2. On the right hand side the motivations for the micromechanical model and the internal variables n (n) are given. The 3D domain structure of a grain is represented by a microscopic material point (MiMP) with six possible polarization directions, that is, p = 6. Each direction is weighted by n (n) , where n = 1, . . . , p.
Following the work of Huber et al. (1999), the internal variables are interpreted as volume fractions of the domain species and have to satisfy the following conditions: The evolution of the irreversible contributions is obtained as follows (Ricoeur and Lange, 2019): The spontaneous strain e sp(n) kl is induced by 908 switching domains, whereas the change of spontaneous polarization DP sp(n) i depends on the kind of switching, see Appendix. The HEAVISIDE-functions H Àdn (n) À Á take the value 1 for domain species n which are switching, thus reducing the associated internal variables, and 0 for those being enriched by the switching. In equation (6) dn (n) describes the change of the volume fraction n (n) , see Figure 2, for a switching process from domain variant n to domain speciesk, where the latter is associated with the maximum energy dissipation which satisfies the necessary condition of exceeding a threshold w crit (n!k) . The dissipative work w diss (n!k) of a switching process from domain species n to k is equal to the negative thermodynamic driving force q em (n!k) and is thus obtained from the potential equation (1) by differentiation with respect to the internal variable (Wingen and Ricoeur, 2019): It should be mentioned, that equation (8) neglects the dependence of the material coefficients on the internal variables. An extension for the dissipative work was suggested by Kessler and Balke (2001) taking into account higher order effects. The dissipative work, outlined in equation (8), was originally introduced in the switching criterion by Hwang et al. (1995). A more detailed derivation of w diss (n!k) is found in Ricoeur and Lange (2019). The evolution law for the internal variables n (n) is formulated as follows: where dn 0 is a model parameter and the HEAVISIDEfunctions as usual take the value 1 for a positive argument including 0, and 0 for a negative value. The first HEAVISIDE-function in equation (9) verifies, if switching from domain n to k is feasible, whereas the second ensures that the one switching process from domain n to domaink associated with the maximum energy dissipation takes place. The energy threshold w crit (n!k) depends on the type of switching. For ferroelectric materials with tetragonal unit cells, motions of 908 and 1808 domain walls are possible. Hence, the energy threshold reads, see for example, Hwang et al. (1995) and Huber et al. (1999): where P 0 and E C denote the spontaneous polarization and the coercive field, respectively, see Table A2. The range of the internal variables outlined in equation (5) implies, that domains can vanish in favor of other domains. From a physical point of view, this behavior is inappropriate since domain walls are indispensable for the reduction of potential energy in a grain. Allowing for unconstrained switching within in the limits of equation (5) gives rise to an overestimation of the irreversible contributions e irr kl and P irr i . In Gellmann and Ricoeur (2016) or Lange and Ricoeur (2015) two additional parameters have been introduced in the constitutive equations to face this deficiency and eventually meet experimental results. Within FE-frameworks it was further required for numerical stability. A more physically motivated approach is the implementation of a lower limit n min for the volume fractions, thus adapting the range of the internal variables according to Since n min .0 it is guaranteed that domains of a MiMP and grain, respectively, cannot vanish. This improved approach introduces just one additional parameter which helps adapting numerical to experimental results. The second model parameter dn 0 of equation (9) just has to be chosen sufficiently small to attain convergence.

Micro-macro-transition with variable grain size
On the macroscopic scale, represented by an RVE, quantities, for example, mechanical stress s ij and electric displacement D i , are microscopic volume averages, see for example, Hori and Nemat-Nasser (1998) or Kessler and Balke (2001): In equation (12) and in the following, macroscopic quantities, obtained by homogenization, are specified by angled brackets. Assuming, homogeneous fields in a MiMP m with the volume V (m) , that is, , the volume averages are given as where M denotes the number of grains in the RVE. In contrast to Lange and Ricoeur (2015) and Ricoeur and Lange (2019), where V RVE = MV (m) has been considered, MiMPs and thus grains of the RVE do not hold equal sizes now. Macroscopic quantities, outlined in equations (12) and (13), finally read where x (m) = V (m) =V RVE describes the volume fraction of the MiMP m. Similar to the volume fraction of a domain n (n) , all x (m) have to satisfy the following conditions: Inserting equations (2) and (3) into equation (14), macroscopic stress and electric displacement are obtained as follows: It should be kept in mind that the constitutive equations of the previous section are interpreted as constitutive equations of a grain m. The material properties C (m) ijkl , e (m) ikl , and k (m) il as well as the irreversible contributions e irr(m) kl and P irr(m) i depend on the internal variables n (n) , see equations (4) and (6). A generalized Voigtapproximation, that is is applied for the sake of a scale transition. The macroscopic constitutive equations (16) and (17) finally read: The macroscopic material properties C ijkl , e ikl h i, k il h i and the inelastic contributions C ijkl e irr kl and P irr i depend on the volume fractions x (m) : Equations (19) and (20) contain known and unknown quantities of an RVE, depending on the given problem.
In particular, these are the macroscopic mechanical stress s ij , electric displacement D i h i, strain e kl , and electric field E l . Prescribing stress and electric field as external loads, that is the macroscopic strain results from equation (19): Inserting equation (23) in equation (20), the macroscopic electric displacement reads Residual stresses s (m) ij and electric displacement D (m) i of a grain m are obtained inserting equation (23) in equations (2) and (3): An inhomogeneous distribution of local stress and electric displacement, as typically observed in polycrystalline materials as a result of grain interactions, is thus represented in the model. Interactions of domains within a grain, however, are neglected at this point. The influence of grain size distributions on macroscopic and microscopic quantities will be investigated based on the following statistical considerations.

Gaussian distribution of grain sizes
The left hand side of Figure 3 exhibits a micrograph of a BT grain. On the right hand side a truncated octahedron, as a model of a three-dimensional grain, is illustrated. The volume of a truncated octahedron is given by (Mendelson, 1969;Weisstein, 2003) where a represents the edge length according to Figure  3. The grain size j being interpreted as the diameter 2R of the circumsphere, the following relation is obtained: It should be mentioned at this point that the volumes of the truncated octahedron and the circumsphere are not identical. The main characteristic of the circumsphere is that all vertices of the truncated octahedron are located on its surface. Inserting equation (28) into equation (27), the volume V (m) is finally given as a function of the grain size j (m) : The volume fraction of a grain m is thus obtained as In the model, the grain sizes are assumed normally distributed. In Figure 4(a) and (b) the probability density function and the volume fraction x (m) according to equation (30), respectively, are plotted versus the grain size j (m) for an RVE with M = 40 grains. As an example, an average grain size j = 5:9mm and a standard deviation s j = 2:95mm are chosen. Figure 4(a) reveals the well known Gaussian distribution, where around 68% of the grain sizes are located in a range of 6s j with respect to the averaged size j. In case of constant grain sizes, that is j (m) = j, equation (30) yields From equations (30) and (32) it is obvious that the volume fraction x (m) does not depend on the volume of the RVE for a given number of grains M, which uniquely determines the volume fractions in the case of uniform grains. For an RVE with 40 grains x = 2:5 % is obtained, which is indicated in Figure 4(b). Since the grain size is Gaussian distributed, around 68% of the volume fractions are in a range of 10 À1 %\x (m) \3 %, whereas for the whole RVE the range is 10 À3 %\x (m) \20%. The MiMP with a grain size j'16:5 mm and a corresponding volume fraction x'17% can be mentioned as one example.

Parametric studies of ferroelectric RVEs
As an example for a ferroelectric material, BT is employed. The set of material parameters is found in the Appendix. For the numerical simulations a pure electric loading into the e 2 -direction, that is is considered. The applied electric load is realized by a trilinear function, where E max 2 = 62 Á 10 6 Vm À1 , see Figure 3. Micrograph of a BT grain taken from Arlt (1990b) and truncated octahedron as model of a three-dimensional grain (in the style of Pearce, 1978).  Before keeping the focus on characteristic magnitudes such as remanent or maximum values for quantitative investigations, Figure 6 shows hysteresis loops of strain and electric displacement taken from the first two load cycles. M = 40 grains are employed in the RVE with the same arbitrary initial domain orientations as a basis for all three calculations. Whereas the butterfly loop, Figure 6(a), exhibits a distinct influence of statistical grain size distribution, represented by relative standard deviations of 30% and 40%, the electric displacement, Figure 6(b), is scarcely affected by nonuniform grain sizes. This effect can be explained by the dominant 1808 domain wall motion, which is, with regard to equation (8), purely electrically driven and occurs instantaneously at 6E C . Due to the generalized VOIGT-approximation, s. equation (18), the electric field E i is assumed homogeneous in the RVE and is prescribed, that is, E i = E ext i . Therefore, 1808 domain wall motion is independent of the relative standard deviation.

Magnitudes of grains and lower limits of domain volume fractions
In this subsection the influence of the number of grains (MiMPs) M and the lower limit of the volume fractions n min on remanent strain and polarization as well as maximum strain and computing time are investigated assuming uniform grain size. For this purpose, RVEs with 20, 30, 40, 50, 75, 100, and 200 grains are considered. Ten calculations with arbitrary initial domain orientations were performed in each case, requiring 70 numerical calculations in total. It has to be kept in mind that since absolute grain sizes are not relevant in the model, the parameter M does not allow conclusions to be drawn about the size of the RVE. The ''representative'' aspect of the RVE is rather given by the number of arbitrary grain orientations increasing with the parameter M.
In Figure 7 error bars of the macroscopic remanent strain e r 22 , maximum strain e max 22 , remanent polarization P r

2
, and the computing time t sim are plotted versus the number of grains M. At this point n min = 0 has been chosen for all calculations. Error bars depict the minimum, maximum, and mean values for each quantity. The red filled areas in Figure 7(a) to (c) represent the range of experimental data, taken from Enderlein (2007), Fo¨rderreuther (2003, and Wang and Li (2020). Having a look at the results of the remanent strain in Figure 7(a) the simulations are in a good accordance with experimental results independent of the number of MiMPs. Here, the mean values as well as the minima  and maxima are within the area of experimental findings. For the maximum strain, see Figure 7(b), the mean values of the numerical simulations are located at the lower limit of experimental results. For 50, 100, and 200 MiMPs, the mean values are even slightly below the experimental range. Concerning the variances, two issues should be highlighted: an increasing number of MiMPs is basically attended by a decrease of the variance and the scattering of the maximum strain is larger than of the remanent strain. Both effects can be explained by the arbitrary orientations of MiMPs. The influence of a specific orientation in an RVE with 20 or 30 MiMPs is much larger than in an RVE with 100 or 200 MiMPs, since the strain e kl is an average quantity, see equation (23).
The number of MiMPs does not have a significant influence on the remanent polarization, Figure 7(c). In comparison to experimental results, however, the latter is larger, exceeding the upper value of the experimental range by about 25%. Since n min = 0 and interactions between domains of a grain are not considered at this point, domains vanish in favor of others, finally overestimating predominantly the polarization. In Figures 8 and 9 the effect of a lower limit of domain volume fractions n min is investigated. For each simulation discussed in Figure 7, five more with n min = 2%, 4%, 6%, 8%, and 10% are considered. Requiring a total of 420 simulations, only the mean values are presented. Figure 8 Values within the red colored areas match the range of experimental data outlined above. The impact of n min on the remanent strain e r 22 is negligible, see Figure 8(a). A significant influence of n min , however, is observed at the maximum strain e max 22 , see Figure 8(b). An increasing n min reducing the amount of domain wall motion, a decrease of the maximum strain is one consequence. A n min ) 0 thus underestimates the maximum strain and should be excluded. Figure 9(a) and (b) illustrates the remanent polarization P r 2 and the computing time t sim , respectively, versus the number of MiMPs and n min . As expected an increasing n min leads to a decreasing remanent polarization and for n min .2 % it is in a good accordance with the experimental data. Compared to the remanent strain of Figure 8(a), n min has a significant influence on the remanent polarization. This effect is caused by 1808 switching processes, having no impact on the irreversible strain. Reducing the magnitude of domain switching, an increasing n min leads to a decreasing computing time, see Figure 9(b). For M = 40 and n min = 4 %, for example, a time saving of 30% is achieved. On the basis of the above investigations, the parameters M and n min can now be chosen reasonably with regard to experimental results and computational costs. In this context, M = 40 and n min = 4 % seem to be a good choice for all following investigations, whereupon e r 22 and P r 2 are in a very good accordance with experimental data, while the mean value of e max 22 = 6:8073 Á 10 À4 is just slightly below the experimental range. The corresponding average computing time is approximately 7 s per load cycle.

Stochastic distribution of grain size
In this subsection, the influence of normally distributed grain sizes on the remanent strain and polarization as well as maximum strain and computing time are investigated. Based on the study outlined in the previous subsection, an RVE with M = 40 MiMPs and n min = 4 % will be considered in the following. Concerning a reasonable range for the standard deviation of ceramic grain structures, relative values ranging from 4.25%, see Chinn (1994), to 35%, see Manosso et al. (2010), are found in the literature, where the standard deviation is normalized with respect to the average grain size. For the numerical simulations an average grain size of j = 5:9 mm is considered, representing a reasonable value for fine-grained BT. Based on the literature, standard deviations between s j = 0:295mm and s j = 2:36mm should be an appropriate range. However, for the sake of a convergence study, standard deviations from 5:9 Á 10 À4 mm to 2:95mm are considered, the former representing the limiting case of a uniform grain size. Again, ten different RVE s are taken into account for each value of s j , whereat the same sets of domain orientations as used in the previous subsection were chosen.
In Figure 10e r 22 , e max 22 , P r 2 and t sim are plotted versus the normalized standard deviation ranging from 0.01% to 50%. According to Figure 10 (Enderlein, 2007;Fö rderreuther, 2003;Wang and Li, 2020). and (b) computing time t sim taken from 10 simulations versus number of grains M and the lower limit of the domain volume fractions n min . Red area in (a) matches the scope of experimental data (Enderlein, 2007;Fö rderreuther, 2003;Wang and Li, 2020). standard deviations above 0.885 mm, that is, s j =j ø 15%, have a noticeable impact on the remanent and maximum strains in terms of the interval widths. The grain size distribution further has no essential influence on the remanent polarization, see Figure  10(c), where the intervals are small anyway. For the computing time t sim , Figure 10(d), no substantial impact is observed either. Figure 11 presents the results of a convergence study. Here, remanent strain and polarization, normalized with respect to the values of the same RVE with identical initial domain orientations, however, a uniform grain size e r, 0 , P r, 0 , are plotted versus the normalized standard deviation in a range from 0.01% to 50%. All asterisks here and in the following represent the same set of grain orientations. The limiting values of e r 22 and  (Enderlein, 2007;Fö rderreuther, 2003;Wang and Li, 2020). for s j ! 0, as expected, correspond to the magnitudes of uniform grain size, that is, lim s j !0 e r 22 s j ð Þ= e r, 0 22 , lim which is confirmed by Figure 11 and emphasizes the consistency of the model in this respect. Figure 11 also illustrates the sudden increase of scattering of strain and polarization for relative standard deviations above 10%. Figure 12 finally presents principal residual stresses at E 2 = E max 2 and E 2 = 0 after unloading. The black circles represent the results of uniform grain size, while the blue asterisks denote those of grain size distribution, where s j = 2:36mm and s j =j = 40%, respectively.
Since s I .s II .s III by convention, the stresses are restricted to the area below the red dash-dotted lines. Independent of the applied electric field, the black circles are concentrated to restricted areas in the space of principal stresses. This issue has already been observed in Lange and Ricoeur (2015), where the principal stresses of the CM were compared to those from FE simulations based on the same micro-mechanical model of domain switching. It should be noted, that in Lange and Ricoeur (2015) a 2D approach was considered, while 3D is adopted in this work, however, not showing a considerable impact on this aspect. For E 2 = 0, Figure 12(b), (d), and (f), there are no significant differences between RVE s with distributed or uniform grain sizes, whereas at E 2 = E max 2 , Figure 12 an influence is observed, although all values are basically arranged in the same areas.

Conclusion
A basically approved microphysically motivated scalebridging polycrystalline constitutive approach has been augmented in two different aspects. Domains are not allowed to vanish in favor of others by introducing constraints to domain volume fractions, and grain sizes are non-uniform, following a Gaussian distribution. Barium titanate, as extensively investigated ferroelectric with unique chemical composition, has been chosen for the sake of experimental validation. Data have been compiled in this regard from various references to compare key quantities of hysteresis loops to results of simulations. Arbitrary orientations of grains in an RVE require a statistical analysis, finally providing a very good agreement with experiments, adapting not more than two model parameters, that is, number of grains in an RVE and lower limit of domain volume fractions. The grain size distribution turns out to have a noticeable impact for relative standard deviations above 10%-15%, leading to an increasing scatter of predominantly strain, whereas polarization and residual stresses are scarcely influenced in their magnitudes. The computing time of one electric load cycle on a MacBook Pro is just a few seconds, emphasizing the efficiency of the modeling approach.

Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Material properties and switching quantities
In compressed VOIGT notation and for a poling into the positive e 1 -direction, see Figure 2, the material properties of BT are introduced as follows: where the elastic, piezoelectric, and dielectric coefficients are listed in Table A1.
A cubic cell is introduced as an interim configuration of a domain switching process to derive the spontaneous strain, related to a local Cartesian coordinate system e i , with e 1 as the dipole axis, see for example, Li and Rajapakse (2007): The coefficient e D describes the elongation of the unit cell in the dipole direction with respect to the reference cubic unit cell of the same volume. For a tetragonal unit cell, e D is defined as (Hwang et al., 1995): Here, c and a 0 denote tetragonal and cubic lattice parameters, where e D , c and a 0 are listed in Table A2. The differential irreversible strain going along with domain wall motion, reducing species n in favor of speciesk, with equation (9) bearing in mind that dn 0 = dn (k) =À dn (n) is a model parameter, see Table A2, and the internal variables for all other species, except n andk, remain unchanged. As an example, applying equation (A4) to a 908-switching in the e 1 -e 3 plane, the corresponding differential irreversible strain reads: Here, the domain n = 1 points into the positive e 1direction and the domaink = 5 into the positive e 3direction, see Figure 2.
The polarization vector of a tetragonal unit cell with e 1 as dipole axis is given by with the spontaneous polarization P 0 , see Table A2. The differential change of the spontaneous polarization from domain n to domaink is given as follows:   Table A2. Physical parameters of BT taken from Jaffe et al. (1971) and model parameter of equation (9).