A phase-ﬁeld approach to studying the temperature-dependent ferroelectric response of bulk polycrystalline PZT

Ferroelectric ceramics are of interest for engineering applications because of their electro- mechanical coupling and the unique ability to permanently alter their atomic-level dipole structure (i.e., their polarization) and to induce large-strain actuation through applied elec- tric ﬁelds. Although the underlying multiscale coupling mechanisms have been investigated by modeling strategies reaching from the atomic level across the polycrystalline mesoscale to the macroscopic device level, most prior work has neglected the important inﬂuence of temperature on the ferroelectric behavior. Here, we present a phase-ﬁeld (diffuse-interface) constitutive model for ferroelectric ceramics, which is extended to account for the effects of ﬁnite temperature by considering thermal lattice vibrations based on statistical mechanics and by modifying the underlying Landau-Devonshire potential to depend on temperature. Results indicate that the chosen interpolation of the Landau en- ergy coeﬃcients is a suitable approach for predicting the temperature-dependent spontaneous polarization accurately over a broad temperature range. Lowering the energy bar- rier at ﬁnite temperature by the aforementioned methods also leads to better agreement with measurements of the bipolar hysteresis. Based on a numerical implementation via FFT spectral homogenization, we present simulation results of single- and polycrystals, which highlight the effect of temperature on the ferroelectric switching kinetics. We observe that thermal ﬂuctuations (at the phase-ﬁeld level realized by a thermalized stochastic noise term in the Allen-Cahn evolution equation) promote the nucleation of needle-like domains in regions of high heterogeneity or stress concentration such as grain boundaries. This, in turn, leads to a faster polarization reversal at low electric ﬁelds and a simulated domain pattern evolution comparable to experimental observations, stemming from the competition between nucleation and growth of domains. We discuss the development, implemen- tation, validation, and application of the temperature-dependent phase-ﬁeld framework for ferroelectric ceramics with a focus on tetragonal lead zirconate titanate (PZT), which we demonstrate to admit reasonable model predictions and comparison with experiments. RVE down atomic i.e., every pixel mimics exactly one tetragonal atomic-level unit cell exhibits temperature-dependent Brownian motion through the space-time random Instead, the utilization of pathways with lower energy barriers leads to switching predominantly by two consecutive 90 ◦ -rotations. Second, thermal noise leads to signiﬁcantly faster growth of a domain nucleus at elevated temperature in both the longitudinal and transverse directions – whose relation | v tip | (cid:20) | v wall | is responsible for the characteristic needle-like shape of ferroelectric domains. In addition, thermal ﬂuctuations promote the branching of existing domains and nucleation of new domains. While the nucleation spots are randomly distributed in a defect-free single-crystal, grain boundaries in a polycrystal (like any other location with stress or charge concentrations) act as natural sites for nucleation. The emerging simulated microstructures during polarization switching incorporate qualitatively various characteristic features known from experimental observation, including ﬁrst- and higher-order laminates, and wedge-like structures. A detailed comparison with experiments is unfortunately out of reach since in-situ measurements of ferroelectric microstructures, especially over a broad temperature range and under applied electric ﬁelds, are a rare ﬁnd. Our simulations capture general qualitative trends while a quantitative comparison will require further experimental data and may require a re-calibration of model parameters (speciﬁcally of the drag coeﬃcient μ , which may also be assumed temperature-dependent in general). Yet, our model demonstrated the salient features of ﬁnite-temperature ferroelectric switching in a promising fashion (based on energetic potentials obtained from ﬁrst principles). We have thus presented an approach to “thermalize” a 0K ﬁrst-principles-based model for ﬁnite-temperature phase-ﬁeld simulations.


a b s t r a c t
Ferroelectric ceramics are of interest for engineering applications because of their electromechanical coupling and the unique ability to permanently alter their atomic-level dipole structure (i.e., their polarization) and to induce large-strain actuation through applied electric fields. Although the underlying multiscale coupling mechanisms have been investigated by modeling strategies reaching from the atomic level across the polycrystalline mesoscale to the macroscopic device level, most prior work has neglected the important influence of temperature on the ferroelectric behavior. Here, we present a phase-field (diffuse-interface) constitutive model for ferroelectric ceramics, which is extended to account for the effects of finite temperature by considering thermal lattice vibrations based on statistical mechanics and by modifying the underlying Landau-Devonshire potential to depend on temperature. Results indicate that the chosen interpolation of the Landau energy coefficients is a suitable approach for predicting the temperature-dependent spontaneous polarization accurately over a broad temperature range. Lowering the energy barrier at finite temperature by the aforementioned methods also leads to better agreement with measurements of the bipolar hysteresis. Based on a numerical implementation via FFT spectral homogenization, we present simulation results of single-and polycrystals, which highlight the effect of temperature on the ferroelectric switching kinetics. We observe that thermal fluctuations (at the phase-field level realized by a thermalized stochastic noise term in the Allen-Cahn evolution equation) promote the nucleation of needle-like domains in regions of high heterogeneity or stress concentration such as grain boundaries. This, in turn, leads to a faster polarization reversal at low electric fields and a simulated domain pattern evolution comparable to experimental observations, stemming from the competition between nucleation and growth of domains. We discuss the development, implementation, validation, and application of the temperature-dependent phase-field framework for ferroelectric ceramics with a focus on tetragonal lead zirconate titanate (PZT), which we demonstrate to admit reasonable model predictions and comparison with experiments.
A key property of ferroelectric ceramics is the coupling of electric and mechanical fields. While primarily used in the linear regime ( Taylor, 1985;Yang, 2006 ), ferroelectric ceramics under sufficiently large electric or mechanical loading enter a nonlinear regime, where a remnant polarization remains after the load is removed ( Bhattacharya and Ravichandran,20 03;Chaplya and Carman,20 01 ). Such permanent changes in the atomic level dipole-structure offer avenues to adjust material properties ( le Graverend et al., 2015 ), induce significant shape changes ( Burcsu et al., 2004 ), or store information ( Buck, 1952 ).
Since the discovery of ferroelectricity in 1921 ( Valasek, 1921 ), numerous ferroelectric materials have been found, among which the family of perovskite oxides is the most technically relevant. At the Curie temperature, the crystallographic structure of perovskites exhibits a phase transition from a high-symmetry, cubic lattice to a lower-symmetry, tetrahedral or rhombohedral phase ( Jaffe et al., 1971;Jona et al., 1957;Shirane and Hoshino, 1954 ). As a result, individual ions shift from their centrosymmetric positions, which leads to spontaneous polarization and spontaneous strains below the Curie temperature ( Lines and Glass, 2001 ). Above the Curie temperature, crystals are centrosymmetric, hence the electric dipole vanishes in the absence of an applied electric field -a quality of the material referred to as paraelectricity. Below the Curie point, the atomic-level polarization of perovskite oxides is electrically alterable, and many are also mechanically alterable ( Chaplya and Carman, 2001 ) -these properties are referred to as ferroelectricity and ferroelasticity, respectively. In this regime, any microstructural rearrangement is accommodated by the nucleation and growth of an intricate network of ferroelectric domains, involving (for reasons of compatibility) primarily 90 • -and 180 • -domain walls, whose temperature-dependent kinetics have remained a challenge to model.
Early studies ( Abe, 1959;Drougard, 1960;Miller and Weinreich, 1960 ) suggested that the motion of a 180 • -domain wall in a defect-free single crystal is driven by the nucleation and growth of triangular shaped domains in a staggered manner. Although such nucleation-driven domain wall motion explains certain experimental observations, it unfortunately fails to predict the required activation fields for nucleation ( Paruch et al., 2006;Tybell et al., 2002 ) and the absolute wall velocity. To fill this gap, Hayashi (1972) proposed an analytical model to account for the kinetics of domain wall motion based on the theory of absolute reaction rates, whereas Logé and Suo (1996) described ferroelectric domain wall motion as a non-equilibrium thermodynamic process, deriving a kinetic model based on variational principles. While the quasistatic material response of idealized ferroelectric ceramics is generally well understood and captured by such models, the complex microstructural mechanisms in real materials -from oxygen vacancies on the atomic level to grain boundary (GB) mechanisms on the polycrystalline mesoscale to boundary conditions on the macroscopic device level -and their influence on the switching kinetics is far less established. As an example, consider the intricate effect of stress concentrations near defects, cracks, and GBs that promote switching and interfere with domain wall motion ( Lambeck and Jonker, 1986;Marincel et al., 2015;Rodriguez et al., 2008 ).
Another open challenge is the rate-and temperature-dependent kinetics of ferroelectric switching ( Arlt and Dederichs, 1980;Merz, 1956;Schultheiß et al., 2018;Wojnar et al., 2014;Zhou et al., 2001 ), which emerges on the atomic level but is strongly influenced by the mesoscale defect network through its impact on domain wall motion and nucleation. This broad range of length and time scales involved presents a challenge for both experimental observation and computational modeling. For better accessibility using TEM imaging, experimental research has focused on ferroelectric switching in thin films ( Chen, 2008;Lohse et al., 2001;Tagantsev et al., 2002 ). The thus observed response, however, does not necessarily capture the behavior of bulk ferroelectrics, since it involves both material and structural effects. Recent experiments by Schultheiß et al. (2018) studied bulk PZT using a fast high-voltage switch setup; the step response of polarization and strain was measured, providing insight into fast-switching kinetics and demonstrating not only clear rate dependence but also a dependence on grain size and texture.
The influence of temperature on the ferroelastic and ferroelectric material response has been assessed experimentally. For the ferroelectric case, the temperature dependence of the piezoelectric and dielectric coefficients of PZT has been measured by electric cycling from room temperature to (or close to) the athermal limit ( Zhang et al., 1994 ). Hooker (1998) performed high-temperature experiments reaching up to 500 K. More recent measurements with co-doped soft PZT covered an even broader temperature range that approached the Curie temperature ( Kaeswurm et al., 2018 ). The influence of temperature on the ferroelastic material properties has been measured in uniaxial compression experiments, e.g., by Ji and Kim (2013) ; Kaeswurm et al. (2018) ; Marsilius et al. (2010) ; Webber et al. (2009) .
When modeling ferroelectric ceramics, three approaches are popular: phenomenological macroscale models, sharpinterface models, and diffuse-interface phase-field models; see, e.g., Vidyasagar et al. (2017) for a discussion and examples. Since we are interested in domain evolution at the mesoscale, we here follow the phase-field approach of Zhang and Bhattacharya (2005) and Su and Landis (2007) , who first modeled ferroelectric domain structures at the mesoscale by solving the electro-mechanically coupled boundary value problem (BVP) based on a finite-element (FE) discretization. By contrast, we adopt Chen 's (2008) spectral strategy to solve the BVP efficiently in Fourier space and specifically adopt the formu-lation of Vidyasagar et al. (2017) . Our ferroelectric constitutive model is derived from the thermodynamic potentials of Völker et al. (2011) , who used first-principles data based on density functional theory (DFT) and atomistic simulations to calibrate the (zero-temperature) electric enthalpy density. While those studies all neglected thermal effects, Landis (2016, 2019) used phase-field methods to characterize the structure of ferroelectric-to-paraelectric phase boundaries near the Curie temperature and derived a thermodynamic framework that accounts for spatially heterogeneous temperature fields. Vopsaroiu et al. (2010) investigated thermally activated switching kinetics by using a non-equilibrium statistical model that describes the polarization switching of a nucleus. Liu et al. (2016) performed molecular dynamics simulations to investigate ferroelectric domain wall motion at finite temperature beyond Merz's law ( Merz, 1956 ). Finite-temperature effects in the continuum phase-field framework, however, have remained a rare find.
In our approach presented here, temperature enters the phase-field description of ferroelectric ceramics in two ways. First, the underlying polarization potential is adjusted to depend on temperature by interpolating between the firstprinciples-informed energy landscape at zero temperature ( Völker et al., 2011 ) and the convex energy potential at the Curie temperature -taking inspiration from the temperature dependence of the order parameter in continuous phase transitions close to the transition point being characterized by a power law and an associated critical exponent ( Toda et al., 1983 ). Such interpolation models, have previously been proposed in the context of, e.g., lambda phase transitions of liquid helium ( Ferrell et al., 1968 ), liquid-gas phase transitions in nuclear reactors ( Panagiotou et al., 1984 ), glass transitions of amorphous oxides ( Ojovan and Lee, 2006 ), paramagnetic-ferromagnetic phase transitions ( Mohan et al., 1998 ), and ferroelectric phase transitions in single-crystalline barium titanate ( Li et al., 2005;Wang et al., 2010;Woldman and Landis, 2016 ). Second, we account for thermal lattice vibrations by a statistical mechanics-based thermalization of the Allen-Cahn evolution equation through temperature-dependent random noise ( Funaki, 1995;Rolland et al., 2016;Shardlow, 20 0 0 ). Related stochastic phase-field models haven been employed to model, e.g., the microstructure evolution in magnetic materials ( Koyama, 2008 ), dendritic crystal growth ( Karma and Rappel, 1999;Shang et al., 2016 ), confined nanoferroelectrics ( Slutsker et al., 2008 ), solidification of austenitic nickel-chromium-based superalloys ( Radhakrishnan et al., 2019 ), plasticity in Ti-alloys ( Zhu et al., 2017 ), and GB motion ( Baruffi et al., 2019 ). Here, we introduce thermal fluctuations to affect the polarization evolution. The resulting finite-temperature phase-field model is validated against experimental measurements (in terms of the ferro-and piezoelectric properties), and we discuss the predicted influence of temperature on ferroelectric microstructures and the associated switching kinetics.

Ferroelectric constitutive model and RVE-problem at zero temperature
We adopt and extend the constitutive model of Vidyasagar et al. (2017) , which was introduced for zero-temperature simulations and which we here summarize briefly to present our modifications and extensions in the proper context. We consider tetragonal perovskite ceramics below their Curie temperature and use continuum mechanics to describe a body ⊂ R n in n -dimensional space. The small strains in brittle ceramics allow the use of linearized kinematics with an infinitesimal strain tensor ε = sym (∇ u ) derived from a (mechanical) displacement field u ( x , t ) : × R → R n , depending on position x ∈ and time t > 0. If inertial and body forces are neglected (as in our applications), the mechanical problem is govern by the balance of linear momentum, which requires with the infinitesimal Cauchy stress tensor σ (and appropriate Dirichlet and Neumann boundary conditions in terms of prescribed displacements and tractions, respectively). The governing equations for the electric problem are derived from Maxwell's equations. Gradients in the voltage potential φ : R n × R → R produce an electric field e = −∇φ, which is connected to the electrical displacement field d : × R → R n and the polarization field p : × R → R n through d = κ 0 e + p , where κ 0 is the permittivity in vacuum. By assuming that no free charges are present within the body , Gauss' law reduces again assuming appropriate Dirichlet and Neumann boundary conditions (in terms of prescribed voltages and surface charges, respectively).
In order to close the above system of equations, we require constitutive relations as well as a dissipative evolution equation for the polarization field. We derive all constitutive relations from the electric enthalpy density W , which for a ferroelectric perovskite decomposes as ( Su and Landis,20 07;Zhang and Bhattacharya,20 05 ) such that The first two terms in (3) arise from decomposing the linear elastic strain energy density into a purely mechanical strain energy mech and an electrostrictive coupling contribution coupl . Writing ε i j = ε e i j + ε r i j , where ε e i j denotes elastic strains and ε r i j remnant strains, we have elastic ( ε ) = 1 2 ε i j − ε r i j C i jkl ε kl − ε r kl with fourth-order elasticity tensor C i jkl , such that the resulting mechanical energy density mech reads whereas the anisotropic electro-mechanical coupling energy, according to Völker et al. (2011) , is expressed as coupl ( ε , p ) = ε i j B i jkl p k p l + p i p j A i jkl p k p l .

(6)
Higher-order coupling tensors F i jklmn and G i jklmn , introduced by Su and Landis (2007) , are not here considered to enforce stress-free conditions on average (as described below). Here and in the following, we use classical index notation with Einstein's summation convention.
The non-convex polarization energy for tetragonal PZT is described by a Landau-Devonshire polarization potential pol , which is assumed as a sixth-order polynomial. We here adopt the potential proposed by Völker et al. (2011) , who calibrated the polynomial coefficients using first-principles DFT data at the athermal limit and exploiting the known crystal symmetries; see Vidyasagar et al. (2017) for a discussion. Finally, the energy contained in the domain walls and the electric energy in free space are, respectively, where G 0 represents an interface energy. In polycrystals the energy density is rotated into the local coordinate system of each grain, transforming all vector-and tensor-valued fields according to the local rotation R ∈ SO( n ).
For a fixed polarization p , (1) and (2) -along with constitutive relations (4) -provide an equilibrium solution of the In reality, the polarization p ( x , t ) evolves over time in a dissipative manner and requires a kinetic evolution law. The latter is usually modeled by the Allen-Cahn equation of a linear gradient flow ( Su and Landis, 2007;Zhang and Bhattacharya, 2005 ): with an inverse mobility (or drag coefficient) μ > 0. The unknown fields u ( x , t ), φ( x , t ), and p ( x , t ) are now obtained from simultaneously solving linear momentum balance (1) , Gauss' law (2) , and the kinetic evolution law (8) , based on the enthalpy (3) and constitutive relations (4) . All material parameters are summarized in Appendix A . We follow Vidyasagar et al. (2017) and compute the effective material response by solving the above equations within a representative volume element (RVE), using spectral homogenization to impose volume-average strains ε and average electric fields e , see Appendix C for details. By assuming elastic homogeneity (which is the case in (an)isotropic singlecrystals as well as in polycrystals when assuming elastic isotropy), we avoid an iterative FFT-based solution scheme and can impose average stresses directly: the average stress in the RVE is σ i j = C i jkl ε kl + B i jkl p k p l , which for a spatially homogeneous C -tensor can be inverted for ε since C i jkl ε kl = C i jkl ε kl . Because samples in our experiments are unconstrained ( Tan et al., 2019;Vidyasagar et al., 2017 ), we assume a negligible average stress in the specimen, so that the average strain components to be imposed are obtained as This allows us to solve for the mechanical strains, displacement field, electric voltage, electric displacements, and the incrementally updated polarization field, using constant time steps t > 0, in a staggered fashion based on the FFT scheme summarized in Appendix C ; for numerical details see Vidyasagar et al. (2017) . In order to reduce computational costs and for ease of visualization, we conduct two-dimensional (2D) simulations in the e 1 − e 3 -plane in most numerical examples. Specifically, we assume vanishing out-of-plane components e 2 = 0 , p 2 = 0 , d 2 = 0 , and σ i 2 = 0 ( i = 1 , 2 , 3 ), whereas ε 22 = 0 in general. This allows us to use 3D material constants (as obtained from first principles) while simulating a planar RVE (thus allowing the out-of-plane strains to accommodate the remnant strains as in a bulk ferroelectric that is stress-free on average). The planar assumption allows for an inexpensive computation of the bulk material response, in which the mechanical and electric fields are restricted to the plane without thin-film effects (e.g., thickness-dependent material response, depolarization fields).

Extension to and effects of finite temperature
To account for temperature dependent material behavior, we modify the above constitutive model as follows: first, we render the polarization potential pol temperature-dependent; second, we append the Allen-Cahn evolution law (8) by a stochastic noise term to mimic the effects of thermally induced lattice vibrations. We acknowledge that this is a first-order approximation; i.e., we modify those terms which, in our view, show the strongest influence on the resulting predicted material response. In principle, we could also account for temperature-dependent elastic and coupling coefficients as well as mobility. Further, one could consider heat conduction, thermal expansion, and thermal heating due to the dissipative evolution kinetics. Here, we assume that all these contributions have a marginal impact on the ferroelectric hysteresis and the microstructural domain evolution compared to the two former aspects taken into account. We therefore assume a uniform known temperature across the RVE and simulate the material response at different temperatures. We note that, for more accurate predictions, the effects of thermal expansion should be included in the model in order to account for secondary pyroelectricity and for thermally induced stresses at GBs in the polycrystalline case. As we simulate unconstrained, elastically isotropic samples under isothermal conditions, we neglect thermal expansion.

Temperature-dependent polarization potential
We exploit our knowledge of the zero-temperature polarization potential of tetragonal PZT from Völker et al. (2011) as well as of the paraelectric phase implying a convex potential landscape at the Curie temperature θ C . In general, the polarization potential at a finite temperature θ is unknown and must be modeled properly between θ = 0 K and θ = θ C . Following Landau (1908Landau ( -1968 and Devonshire (1949) , we introduce a linear interpolation of the polarization potential with respect to temperature θ , such that the polarization enthalpy density in 3D becomes where p i = p · e i denotes the polarization component in the x i -direction, and the Cartesian unit vectors e i ( i = 1 , 2 , 3 in three dimensions) are chosen to align with the tetragonal crystal axes 100 , 010 , and 001 . α 1 through α 123 are material constants adopted from the DFT-based 0K potential of Völker et al. (2011) ; see Appendix A .
Due to symmetry of the tetragonal unit cell, the polarization potential pol ( p , θ ) has six minima (and 2 n minima in n dimensions in general). In those ground states, the polarization is aligned with one of the tetragonal crystal axes and p (θ ) = ±p 0 (θ ) e i , where p 0 ( θ ) > 0 denotes the (now temperature-dependent) spontaneous polarization. Consider now a ferroelectric single-crystal forming a single, homogeneous domain, whose polarization is aligned with one of the tetragonal crystal axes. In the absence of any external mechanical or electrical loading ( e = 0 ), minimizing the electric enthalpy (which is equivalent to minimizing (10) ) with respect to the polarization and considering only positive and real solutions p 0 ∈ R + identifies the temperature-dependent spontaneous polarization as We note that (11) is identical to the theory of Devonshire (1949) only in the limit | (θ C − θ ) /θ C | 1 , in which case a Taylor expansion of (11) results in the classical relation p 2 0 = β(θ C − θ ) /α 11 with a constant β > 0. By contrast, we here do not make this simplifying assumption since we aim to cover the full temperature range from 0K to the Curie point (and we will demonstrate that retaining the exact form (11) is important to arrive at accurate predictions). A further intrinsic ferroelectric property, which is predicted by Landau-Devonshire theory, is the coercive field e c , which refers to the electric field required in a single-crystal for complete 180 • polarization reversal. Considering a stress-free singledomain single-crystal with an applied electric field (aligned with a tetragonal crystal axis), we solve ∂ W pol ( p , e , θ ) /∂ p i = 0 using (10) with p = p e i , e = e e i , to find a relation between the electric field e and the equilibrium polarization p (at a given This relation is visualized in Fig. 1 for various temperatures. The coercive field corresponds to the local maximum in the electric field (illustrated as dashed lines in Fig. 1 ) and hence follows from solving ∂ 2 W pol /∂ p 2 i = 0 for p * = p * (θ ) and inserting the solution into (12) so e c (θ ) = | e (p * (θ ) , θ ) | . We omit the lengthy analytical solution here. For this 1D scenario, the temperature-dependent polarization potential and its corresponding electric field are plotted as functions of the polarization for various temperatures in Fig. 1 . The minima in the polarization potential are located at ± p 0 (θ ) , whereas the coercive field is identified as the points of bifurcation in the electric hysteresis.

Thermal fluctuations via stochastic noise
While the above temperature-dependent potential reflects variations in the spontaneous polarization and coercive field, it affects the kinetics of ferroelectric switching only through changes in the driving force (due to changes in the energy landscape). This, however, neglects another important effect of temperature. In any ferroelectric sample, the abundant network of defects (including point defects such as oxygen vacancies as well as higher-dimensional defects such as GBs and existing domain walls) serves as nucleation sites for the heterogeneous nucleation of new domains, while also impeding domain wall motion through pinning and drag effects ( Jo et al., 2009;Puchberger et al., 2017 ). Such mechanisms are generally Influence of temperature on the polarization potential pol (p, θ ) (left) and the corresponding electric field e ( p, θ) (right) for a single-domain singlecrystal and an electric field parallel to the polarization. The minima of pol (p, θ ) are at ± p 0 ( θ), whereas the local maxima/minima of the electric field represent the coercive field ±e c (θ ) (indicated as dashed lines). temperature-dependent, and one underlying causal mechanism are atomic lattice vibrations whose amplitude grows with temperature. Although generally being of small amplitude compared to the atomic unit cell rearrangements during ferroelectric switching, these small perturbations can be sufficient for promoting nucleation and growth of domains by helping the material to locally overcome the respective energy barriers. Simply put, not only do energy wells in the potential of Fig. 1 become shallower with increasing temperature, but also do atoms fluctuate at higher amplitude within those wells, which promotes switching to the respective other well and hence increases the escape rate, especially near lattice defects.
Because it is neither possible nor desirable to compute the motion of individual atoms inside the RVE, we here use a statistical mechanics-based approach to capture the influence of atomic vibrations by amending the polarization kinetics to include a term of Brownian motion at the RVE-/mesoscale. Specifically, we turn the Allen-Cahn Eq. (8) into the stochastic form in which η( x , θ ) represents a random noise term that mimics the effect of lattice vibrations. To comply with the second law of thermodynamics, we consider only conditions of constant uniform temperature within the RVE. We require that the noise term satisfies the following constraints (for a detailed derivation of the statistical mechanics considerations see Appendix D ): 1. For a truly stochastic noise that does not bias the evolution of the polarization field p in any direction, the random noise must average to zero over time at any point within the simulated RVE: for all x ∈ , for any sufficiently large time window τ > 0 .
2. The random noise must average to zero over the RVE at any given time: 3. The random noise is uncorrelated in space and time, and its variance σ 2 depends on temperature θ and time increment t according to where V char = a 2 tetr c tetr is the volume of the perovskite's atomic unit cell, which is used for normalization of the thermodynamic potential. We do not account for local variations in the lattice volume V char , which could be important for a more accurate representation of GBs and multiple, low-symmetry phases. As detailed in Appendix D , the correlation constraint (16) stems from a statistical mechanics consideration, interpreting the random noise term analogous to a random walk whose overall effect, over sufficiently long times, obeys a Boltzmann-type equilibrium probability distribution. By solving the associated Fokker-Planck equation in the equilibrium limit, the above condition (16) emerges.  Table A2 ). The predicted normalized (single-crystal, single-domain) coercive field e c (θ ) is also included (dashed blue line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) The modified Allen-Cahn Eq. (13) for the polarization field, along with the conditions (14) -(16) , effects a kinetic evolution of the polarization that depends on temperature in a stochastic sense -and the thermal fluctuations grow with increasing temperature according to (16) . To enforce the above conditions in practice, we pick real, uncorrelated random numbers out of a standard normal distribution N (μ, σ 2 ) with mean μ = 0 and variance σ 2 = 1 . This is achieved, e.g., by the Muller-Box sampling ( Muller, 1958 ) or the polar method of Marsaglia and Bray (1964) . For 3D simulations, noise is generated by picking random numbers { x 1 , x 2 , x 3 } at each time step and for each point inside the RVE, so that rescaling gives the sought random noise (at each discrete time step and at each point) as We point out that the stochastic Allen-Cahn equation is assumed ill-posed for dimensions n ≥ 2 (i.e., its continuum limit does not have a reasonable meaning), which may introduce mesh dependence ( Ryser et al., 2012 ). In our scenario, however, there exists a natural, finite length scale, since the electric dipole within the atomic unit cell is the smallest unit exposed to lattice vibrations acting on the surrounding ions. Hence, the size of the atomic-level unit cell (of volume V char ) provides a physical length scale that relates the random noise to the numerical discretization x used in simulations. Choosing x at the level of the atomic unit cell hence provides a reasonable solution. (While shrinking the mesh size below the atomic unit cell is physically questionable, coarser grids generally underestimate the number of possible nucleation sites and therefore slow down the switching kinetics.) An alternative would be to regularize the noise with a correlation length that depends on the length scale of the dipole-dipole interactions (see, e.g., Kohn et al. (2007) for an investigation of the stochastic Allen-Cahn equation at the sharp-interface limit by using large-deviation theory). Here, the interface energy introduces a characteristic length scale for dipole-dipole interactions, which acts as a natural regularization by limiting the impact of a unit cell's noise on its neighbors.

Influence of the temperature-dependent polarization potential
To assess the accuracy of the chosen linear interpolation of the polarization potential with temperature, Fig. 2 illustrates the spontaneous polarization p 0 vs. temperature -comparing computed results obtained from the linearly interpolated polarization potential (10) as well as from the approximation by Devonshire (1954) to experimental data for different types of PZT. Unfortunately no complete set of data for a single type of PZT across the full temperature range is available to our knowledge. Hence, for an accurate comparison the temperature and spontaneous polarization are normalized by, respectively, the Curie temperature θ C and the extrapolated polarization at 0K, p 0 (0) , for each material (see Table A2 for the exact reference values used for normalization). We note that the drop in the experimental data of Hooker (1998) at low temperatures is questionable in our view (one may question whether complete polarization reversal was achieved at those low temperatures, since all other data clearly report a different trend). If we ignore the low-temperature data of Hooker (1998) , the normalized spontaneous polarization measurements in Fig. 2 coincide reasonably well with the prediction by our finitetemperature model for all shown PZT compositions (also demonstrating the continuous, second-order phase transition expected for PZT).
Important characteristics of ferroelectrics are their electric hysteresis and butterfly curve, which we extract from singlecrystal RVE simulations at different temperature levels, using bipolar electric field cycling. To this end, a triangular-shaped average electric field in the x 3 -direction with amplitude e 3 and cycling period T is applied; simultaneously, the average electric displacement parallel and the average strain perpendicular to the electric field, d 3 and ε 11 , respectively, are recorded.
Numerical results of the bipolar switching hysteresis at different temperatures, computed with a single-crystal 2D RVE of grid resolution 256 × 256, are plotted in Fig. 3 (a). Similar to the predicted temperature dependence of the coercive field in Fig. 2 , we notice an approximately linear decrease of the coercive field coercive field e c = e | p=0 with increasing temperature. In Fig. 3 (c) we plot the normalized polarization (i.e., the polarization normalized by its value p 0 = p| e =0 for each temperature) vs. the normalized electric field (i.e., the applied electric fields normalized by the coercive field at each temperature). The normalized polarization at zero electric field is approximately 1 across the full temperature range tested, so the polarization converges to its equilibrium state, which implies that the simulation indeed captures the quasistatic material response at the chosen cycling rate. As the only exception, results for 600K reveal a polarization at zero electric field that is considerably higher than the spontaneous polarization at that temperature, so we observe a strong effect of temperature on the hysteresis. Also included in Fig. 3 (c) are experimentally measured data for (polycrystalline) PZT-5A at room temperature ( Tan et al., 2019 ), whose normalized curve agrees well with the simulated hystereses.
As a further characteristic of ferroelectric ceramics we compute the evolution of strain with electric field. Fig. 3 (b) plots the negative lateral strain ε 11 vs. the applied electric field as the classical butterfly curve. Analogous to the polarization hysteresis, an increase in temperature leads to a decrease of the electric field at maximum strain ( e ε c ), which is slightly higher than the corresponding field from the polarization hysteresis curve ( e c ). Furthermore, we observe a decrease of the strain magnitude from polarization reversal with increasing temperature. The corresponding normalized curves are shown in Fig. 3 (d), again indicating good qualitative agreement with experimental data.
From the bipolar switching hysteresis and the butterfly curves, small-signal properties such as the piezoelectric coefficients d 31 , d 33 and dielectric constants κ 11 , κ 33 can be determined as (no summation implied) The dielectric constant or relative permittivity κ 33 (the slope of the polarization hysteresis at zero electric field) is a measure of the capacitance of a medium. The piezoelectric coefficients d 31 and d 33 (the slopes of the strain perpendicular and parallel to the switching direction, ε 11 and ε 33 , respectively, at zero electric field) provide a relation between the induced strain and the applied electric field and can be interpreted as a force sensitivity (i.e. charge released per Newton force). The temperature dependence of the piezoelectric coefficient and of the dielectric constant as obtained from our phase-field model is shown in Fig. 4 in comparison with experimental data. Since the bipolar switching hysteresis and the butterfly curve depend strongly on a particular material's microstructure and composition (i.e., its grain size and texture, defect distribution, titanium concentration, dopants, etc.) which are not considered in our model, all reported small-signal properties are normalized with respect to their value at 300 K. The overall trends of the temperature-dependent piezoelectric coefficients d 31 and d 33 are captured reasonably well, independent of the specific ferroelectric ceramic (and unbiased by microstructural variations). We note that the dielectric constant κ 33 shows a stronger dependence on the particular material. Our model (based on the first-principles-informed 0K potential of Völker et al. (2011) ) comes closest to Hooker 's (1998) measurements of PZT-5A. However, effects at the polycrystalline mesoscale, such as domain wall motion and defect pinning, are known to have an impact on the large-signal and small-signal properties. Considering that we used a single-crystal in simulations, the agreement with measurements is reasonably good.

Influence of thermal fluctuations
To assess the impact of the thermalized random noise on the ferroelectric switching kinetics, we deliberately deactivate the temperature dependence of the polarization potential (discussed in the previous section) in order to isolate the effect of the stochastic noise (this ensures that varying the temperature does not alter the coercive field, so that a constant applied electric field is a legitimate test case for evaluating the influence of the introduced random noise for varying temperature levels). The combined effects of temperature-dependent potential and thermal noise will be investigated in the following Section 4.3 . Subsequent numerical examples use a 2D RVE with 1024 × 1024 grid points and resolve the ferroelectric microstructure down to the atomic level; i.e., as discussed before, every pixel mimics exactly one tetragonal atomic-level unit cell and exhibits temperature-dependent Brownian motion through the space-time random process.  As an instructive scenario, we use the well-defined environment of a single-crystal to study the kinetics of domain nucleation and growth under the influence of thermal noise. To initialize the nucleus, we seed an elliptically shaped a + -domain at the center of the RVE, as depicted in Fig. 5 (a), and -for its stabilization -apply a constant electric field e 3 = 10 8 V/m significantly below the coercive field (which in this case is e c = 5 · 10 8 V/m). At varying noise levels, we observe the isolated nucleus grow in two directions: in the longitudinal direction (spreading with the needle-tip speed v tip ) and in the transverse direction (accommodated by classical domain wall motion at a speed v wall ), see Fig. 5 (a). As shown in Fig. 5 (b) and (c), both velocities are strongly influenced by the thermal fluctuations, with the propagation speeds increasing approximately linearly with temperature and consistent with 0K results obtained at the athermal limit without thermal fluctuations. Independent of temperature, the needle-tip velocity v tip is considerably higher (about a factor of 33) than the domain wall velocity v wall , which consequently results in a slender, needle-like shape of the growing a + -domain. This predominant growth in the longitudinal direction is illustrated in Fig 5 (d), showing the computed polarization in the vertical direction p 3 ( x , t )/ p 0 after an elapsed time of t = 170 μ/ | α 1 | . Fig 5 (f) reveals that the bulk of the nucleus occupies an equilibrium polarization state ( a + -domain) with a low polarization energy, whereas the domain walls and the needle-tip are in a non-equilibrium polarization state with a locally high polarization energy. This high-energy polarization state makes the needle tip and walls prone to thermally-driven switching due to the lower energy barrier E 90 in the polarization energy pol ( p , θ ) that stands in competition with the thermal energy k B θ .

Combined effects of thermal fluctuations and temperature-dependent energetics
To understand the behavior observed when including both the temperature-dependent polarization potential (affecting the energetic switching barriers) and the thermal noise (causing fluctuations that help overcome those barriers), we illustrate in Fig. 6 (a) a typical landscape of the polarization enthalpy density W pol ( p , e , θ ) vs. the (normalized) polarization p = (p 1 , p 3 ) in 2D, at a fixed applied electric field e 3 = 8 · 10 7 V/m and temperature θ = 300 K . Consider as the initial state p = (0 , −p 0 ) T . Under the applied field, switching from p = (0 , −p 0 ) T to p = (0 , p 0 ) T is most easily accommodated by two subsequent 90 • -switching events. The minimum energy pathway (MEP) connecting those two polarization states is obtained by using the simplified string method ( Sheppard et al., 2008 ) and is indicated as a magenta curve in Fig. 6 (a). Plotting the polarization enthalpy density along this MEP reveals the energy barrier E 90 of a 90 • -domain wall, see Fig. 6 (b). (The barrier for 180 • -switching is significantly higher.) As summarized in Fig. 6 (c), the energy threshold E 90 , which separates two 90 • -adjacent polarization states, depends on the applied electric field e 3 as well as on temperature θ , the latter dependence enters through the polarization potential pol ( p , θ ) introduced in Section 3.1 . Data in Fig. 6 (c) indicates that increasing the temperature reduces the energy barrier for 90 • -switching, so that maintaining a constant applied electric field induces domain switching more readily with increasing temperature. This effect becomes apparent in Fig. 7 , which shows the same single-crystal example from Fig. 5 but this time at the three temperatures θ = 275 K, 325 K, and 375 K, while applying the same electric field e 3 = 8 · 10 7 V/m. We observe three distinct switching mechanisms: (a) growth of the nucleus predominantly as a needle in the longitudinal direction at 275 K, (b) branching of the existing a + -domain into multiple a − -domains at 325 K, and (c) nucleation of mainly a + -domains at randomly distributed locations inside the c − -domain at 375 K. For the given choice of temperature and electric field, classical domain wall motion perpendicular to the wall plays only a minor role 1 , which is in agreement with experimental observations and analytical considerations ( Ayoub et al., 2017;Hayashi, 1972;Meng et al., 2015;Merz, 1956 ). The temperature dependence of the polarization potential hence globally reduces the energy barrier at elevated temperatures and stimulates thermally-driven polarization reversal by branching of existing domains and nucleation of new ones.
The conditions of the three snapshots in Fig. 7 correspond to points A, B, C highlighted in Fig. 6 (c). By more broadly covering the space of electric fields and temperatures, numerical simulations were used to identify regions in Fig. 6 (c) in which polarization switching occurs primarily by growth only as in Fig. 7 (a) (blue shaded area), growth, branching, and nucleation as in Fig. 7 (c) (red area), or shrinkage and extinction of the nucleus (white area). This illustrates the competing microstructural mechanisms and the influence of temperature and electric field. Other fluctuation fields (not considered in this work, caused, e.g., by thermally-driven migration of oxygen vacancies or free charges) are expected to have a similar effect as lattice vibrations. On the other hand, microstructural imperfections such as GBs, lattice defects, cracks and voids result in localized high-energy spots, leading to heterogeneous nucleation instead of at random locations as seen in the single-crystalline RVE in Fig. 7 (c).
To probe the impact of heterogeneity, we simulate a polycrystalline RVE of PZT with randomly-oriented grains, whose orientations are assigned based on a Gaussian distribution with zero mean and 22 • standard deviation. The sample is poled initially in the negative vertical direction, resulting in a single c − -domain. After equilibration, the polarization adjusts slightly according to the preferred orientation of each grain. Finally, an electric field is applied and kept constant during the domain evolution, comparable to the experimental step-load procedure described by Schultheiß et al. (2018) (this step-response loading shows the system kinetics in a clean fashion without dependence on, e.g., the frequency during bipolar electric cycling). Fig. 8 shows various snapshots of the same simulated ferroelectric microstructure, showing the normalized polarization p 3 ( x , t )/ p 0 in the vertical direction, the polarization energy density pol ( p , θ ) , and the elastic energy density mech ( ε ) .
The grain orientations within the RVE are shown schematically in Fig. 9 . Analogous to the single-crystal example, we observe that nucleation of new a -domains inside a c − -domain lowers the polarization potential in the bulk, but the domain wall and the needle-tip of the nucleus remain in a non-equilibrium polarization state and are therefore energetically unfavorable. The red spots in the polarization energy map indicate locations where the energy has reached the threshold of a 90 • -domain wall (for a grain with zero misorientation); since these are unstable states, an immediate polarization switching can be expected. We also illustrate the elastic energy density, which highlights locations of high stress concentrations, such as grain triple junctions, mismatching domain interfaces caused by c / a -lattice distortion, and perpendicular branches   growing out of existing domains. The shown microstructure reveals primarily 90 • -domain patterns arranged in laminate structures, including more complex domain patterns such as second-order laminates or (in the magnified view of the polarization distribution) a wedge-like microstructure along the horizontal GB, reminiscent of ferroelectric domain patterns observed experimentally; see, e.g., the TEM images of Schmitt et al. (2007) ; Woodward et al. (2005) .
Figs. 9 and 10 visualize the influence of temperature on the evolution of such ferroelectric microstructures at the three temperatures θ = 275 K, 325 K, and 375 K. Under a constant electric field of e = 8 · 10 7 V/m, applied instantaneously at t = 0 s, the polarization evolves, whose average p 3 ( x , t )/ p 0 is shown in Fig. 9 , while the corresponding microstructures at the strain polarization levels indicated as A through C, at three different temperatures, are illustrated in Fig. 10 . That is, the shown microstructures within each row of Fig. 10 have the exact same average polarization p 3 ( x , t )/ p 0 (and the same applied electric field) but the underlying microstructures differ significantly due to the three distinct temperature levels.
The domain pattern evolution, also shown in Movies S1 and S2 (see supporting online material), shows two distinct switching mechanisms: (i) various nucleation events of needle-like domains at GBs and triple junctions, followed by (ii) subsequent growth -predominantly in the longitudinal (needle-tip) direction. The ratio of the speeds of the aforementioned mechanisms is an important factor that determines the appearance of ferroelectric microstructures. Increasing the temperature generally leads to more detailed and finer domain structures. This is traced back to the competition between nucleation and growth, yet we reiterate that two competing effects are at play here. On the one hand, with increasing temperature the polarization energy landscape becomes shallower and the coercive field e c is reduced, see Fig. 6 ; this increases the number of possible nucleation sites at a constant electric field with increasing temperature, so it becomes easier to overcome the energy barrier between adjacent spontaneous polarization states. On the other hand, the noise amplitude increases with temperature ( | η| ∝ √ θ ), so the larger step size of the random walk enables statistically more locations to escape from local energy minima to energetically lower polarization states. These two effects explain the temperature dependence of the domain nucleation sites, which becomes apparent in Fig. 10 : nucleation at low temperature occurs primarily at GBs, while at high temperature the formation of new domains is not restricted to locations with high stress concentrations. Instead, the shallow energy landscape ( Fig. 6 c) in combination with higher thermal fluctuations allows the random walk to overcome the energy barrier of the polarization potential, resulting in nucleation at random locations -similar to the single-crystal results in Fig. 7 c. (Note that the chosen periodic boundary conditions of the computational domain do not affect the nucleation sites directly, but the propagation of growing domains across the boundaries can trigger the nucleation of new domains indirectly.) For completeness, Fig. 9 also includes (as dashed lines) results obtained without the stochastic noise (so the temperature dependence stems solely from the polarization potential), which highlights the impact of the fluctuations: with thermal noise, we observe a considerably faster response time, which is explained by the increasing nucleation rate of a -domains and, as a consequence, polarization reversal being dominated by nucleation as opposed to domain growth.

Conclusions
We have presented a finite-temperature continuum model for ferroelectric ceramics, which is based on a temperaturedependent Landau-Devonshire potential and on a temperature-dependent stochastic Allen-Cahn equation for the evolution of the total polarization. The former was shown to provide an accurate prediction of the spontaneous polarization, the coercive field, and the piezoelectric and dielectric constants across a broad temperature range in agreement with experimental data for PZT (after normalization of the electric field by the coercive field, the computed butterfly curve also showed convincing agreement with room-temperature measurements). Because of the large spread among measured data for different PZT compositions and the fact that we do not account for dopants in the model, the piezoelectric and dielectric constants required normalization for comparison. However, when considering that we compare measurements from different types of polycrystalline PZT at the macroscale with numerical results computed with a single-crystalline 2D model at the mesoscale, the presented framework captures the salient macroscopic temperature effects reasonably well.
Based on statistical mechanics, we introduced a temperature-dependent Gaussian noise into the evolution equation for the polarization, which mimics atomic-level lattice vibrations at the continuum scale. Typical for diffusive processes such as ferroelectric domain wall motion, the noise amplitude is proportional to the square root of temperature and time increment. Simulations revealed that the thermal noise has a considerably effect on the ferroelectric switching kinetics. First, superimposing random small perturbations onto the deterministic gradient-flow kinetics breaks the symmetry of the singlecrystalline polarization energy, such that 180 • -switching becomes less probable. Instead, the utilization of pathways with lower energy barriers leads to switching predominantly by two consecutive 90 • -rotations. Second, thermal noise leads to significantly faster growth of a domain nucleus at elevated temperature in both the longitudinal and transverse directions -whose relation | v tip | | v wall | is responsible for the characteristic needle-like shape of ferroelectric domains. In addition, thermal fluctuations promote the branching of existing domains and nucleation of new domains. While the nucleation spots are randomly distributed in a defect-free single-crystal, grain boundaries in a polycrystal (like any other location with stress or charge concentrations) act as natural sites for nucleation. The emerging simulated microstructures during polarization switching incorporate qualitatively various characteristic features known from experimental observation, including first-and higher-order laminates, and wedge-like structures. A detailed comparison with experiments is unfortunately out of reach since in-situ measurements of ferroelectric microstructures, especially over a broad temperature range and under applied electric fields, are a rare find. Our simulations capture general qualitative trends while a quantitative comparison will require further experimental data and may require a re-calibration of model parameters (specifically of the drag coefficient μ, which may also be assumed temperature-dependent in general). Yet, our model demonstrated the salient features of finitetemperature ferroelectric switching in a promising fashion (based on energetic potentials obtained from first principles). We have thus presented an approach to "thermalize" a 0K first-principles-based model for finite-temperature phase-field simulations.
We close by pointing out that we linked the stochastic noise in our model to thermal lattice vibrations. By using as the normalization volume the primitive unit cell of the crystal lattice (known from DFT calculations), we ensure that the noise amplitude is intrinsically connected to material properties without any fitting parameter. As a downside, this restricts simulations to small length scales (effectively limiting the pixel or voxel size to that of an atomic unit cell), as demonstrated in the presented examples with RVEs at the nanoscale. This also results in realistic domain wall thicknesses in simulations, not achievable at considerably larger scales. Note that one may alternatively interpret the introduced fluctuation field at larger scales, e.g., as the joined impact of temperature and fluctuating point defects and charges on the mesoscale, in which case larger spatial simulation domains are feasible but at the cost of rendering the random noise phenomenological and its amplitude a fitting parameter. Irrespectively, we conclude that the presence of random fluctuations is key to achieving realistic predictions of ferroelectric microstructures not predictable in a perfect, noise-free system.

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.

Appendix B. Homogenization Problem
The constitutive material law introduced in Section 2 describes the behavior of a ferroelectric perovskite at the singlecrystal, single-domain level. The transition from that scale to the macroscale is accomplished by computing the effective response of a representative volume element (RVE), as is customary in classical first-order homogenization (see e.g. Miehe et al. (2002) ;Schröder (2009) ). For a sample with an approximately statistically homogeneous microstructure, we hence define an effective property as the volume average over the RVE, writing where | | denotes the volume of the RVE. We solve the balance of linear momentum (1) and Gauss' law (2) , while imposing periodic boundary conditions over the surface (or boundary in 2D) of the RVE. Technically, we decompose the RVE boundary into opposite parts such that ∂ = ∂ + ∪ ∂ − , and we enforce where x + and x − are pairs of opposing points on ∂ + and ∂ − , respectively. The volume-averaged strain is denoted by ε 0 = ε and volume-averaged electric field is e 0 = e . For homogenization to work, we assume a separation of scales and we postulate that body and inertial forces are negligible since we are interested in the quasistatic material behavior; hence the mechanical and electrical RVE problems are solved quasi-statically.
The only time-dependent governing equation, the modified Allen-Cahn Eq. (13) , is solved by assuming the periodic boundary conditions which does not impose an average but instead allows the polarization field to evolve freely (aside from periodicity on the RVE surfaces).
In our experiments ( le Graverend et al., 2015;Tan et al., 2019;Wojnar et al., 2014 ), a uniform electric field ē = φ/h is applied over the specimen thickness h , with φ denoting the corresponding voltage differential. In this setup, the electric field is applied at the macroscale, thus we assume a separation of scales. The measured electric charge Q can be directly linked to the average electric displacement d = Q/A, viz. ( Vidyasagar et al., 2017 ) (B.4) where A is the area of the electrodes, z the unit vector pointing through the sample thickness and q s the charge density. We emphasize that our study is dedicated to the bulk response of ferroelectric ceramics and not to thin films exhibiting considerable free-surface effects. As a consequence, the depolarization field plays only a marginal role, such that the above relations hold and we may assume that the applied voltage differential can be directly interpreted as the average applied electric field at the RVE-level. Finally, since both the electric field ē = e = e 0 and the electric displacement d = d are related to their measured counterparts, the average polarization p (based on an isotropic permittivity) is obtained as

Appendix C. Spectral Solution Scheme
We follow the approach of Vidyasagar et al. (2017) and solve all governing equations in Fourier space, encouraged by the periodic homogenization scheme. To this end, we discretize the RVE into N grid points in each dimension, such that the position vector over all grid points becomes x = { x 1 , . . . , x N } . For any function f ( x ), we define its inverse discrete Fourier transform as where k denotes the wave vector in the reciprocal lattice (the complete set being T ), and ˆ f ( k ) are the Fourier coefficients.
We solve Gauss' law in Fourier space to obtain the complex voltage potential and electric field as, respectively, Since we assume linearized kinematics and elastic isotropy, we can analogously solve the balance of linear momentum directly in Fourier space -without the need for stress perturbation terms as classically done for heterogeneous media ( Lebensohn et al., 2012;Moulinec and Suquet, 1998;. To avoid such an iterative solution scheme, we assume a homogeneous material and approximate an isotropic elastic material behavior with Voigt stiffness moduli C 11 , C 12 (as obtained from first principles by Völker et al. (2011) ) and define C 44 = (C 11 − C 12 ) / 2 , so that the components of the fourth-order elasticity tensor can be written as C i jkl = λ e δ i j δ kl + μ e (δ ik δ jl + δ il δ jk ) with Lamé moduli λ e and μ e . Of course, this is a simplifying assumption and we admit that the elastic anisotropy may have an impact on the ferroelectric response (especially when considering, e.g., the elastic mismatch near GBs). However, given the variations of reported (experimental and computed) elastic moduli, especially for PZT near the morphotropic phase boundary, we are not in a position to quantify the exact influence of elastic anisotropy and therefore limit our study to isotropy.
Consequently, the total strains ε = ε e + ε r decompose additively into the elastic strains ε e and the coupling strains ε r , where the latter stems from the coupling energy density. Defining the associated stress tensor as σ r = ∂ coupl /∂ ε and applying the Fourier transform to the balance of linear momentum leads to where A ik ( k ) = C i jkl k j k l is the acoustic tensor, and ˆ σ r nm ( k ) represents the Fourier-transformed coupling stresses (which are computed in real space from the constitutive law (4) ). Eq. (C.3a) can be solved directly in Fourier space without iterations, which enables fast and efficient simulations without the need for computing or storing a consistent tangent (hence enabling the presented high-resolution simulations).
Since samples in experiments were unconstrained, we assume a negligible average stress in the sample and hence inside the RVE: σ i j = C i jkl ε kl + B i jmn p m p n = 0 .

(C.4)
Exploiting the assumption of elastic isotropy, (C.4) allows us to compute the average strain in the RVE as ε 0 kl = ε kl = C −1 i jkl σ i j − B i jmn p m p n = −C −1 i jkl B i jmn p m p n = ε r kl . (C.5) We discretize the modified Allen-Cahn equation by an implicit backward-Euler finite-difference scheme, based on constant time increments t > 0 such that t n = n t. Thus, (8) is turned into , j (C.6) for every grid point inside the RVE, where p n +1 = p ( x , t n + t ) and p n = p ( x , t n ) denote the polarization at the current and previous time increment, respectively. For convenience, we define a thermodynamic driving force with components and the nonlocal gradient term, which is computed in Fourier space, is evaluated as The overall ferroelectric problem is solved in a time-incremental, staggered manner. We first solve for the electric field ˆ e n +1 and the strains ˆ ε n +1 in Fourier space, using (C.3a) and (C.2) based on the polarization p n from the previous time step. Applying the inverse Fourier transform yields the real-space quantities e n +1 and the strains ε n +1 at all RVE grid points. Next, using implicit Euler time integration to solve leads to the sought new polarization p n +1 . The time step size t was verified by numerical experiments to be sufficiently small to achieve convergence of this staggered scheme.

Appendix D. Derivation of the stochastic noise term
For a thorough discussion of the stochastic noise term introduced in the Allen-Cahn equation to mimic thermal fluctuations, we briefly revisit the random walk concept and Brownian motion, starting in 1D for simplicity. We start with the Langevin equation ( Langevin, 1908 ), considering only the overdamped solution (inertial terms are neglected) with a polarization p that is attached to its equilibrium position p 0 through a potential W ( p ) and has an inverse mobility μ. With the added random-walk term η( t ), the equation of motion becomes The double-well potential W ( p ) keeps the polarization close to the spontaneous polarization p 0 and prevents it from drifting over time. Therefore, we expect that the variance of the polarization does not diffuse to zero over time but that the distance from the equilibrium position remains bounded, so that states far from p 0 (of high energy) become unlikely (even more so then before). The viscous damping slows down the polarization motion and we expect that over long times the polarization may not assume an equilibrium position (the random noise prevents this) but will attain an equilibrium distribution with constant mean and variance. One may expect that over long times ( t → ∞ ) this process attains a thermal equilibrium, for which the probability of finding a polarization p is given by a Boltzmann distribution ( Boltzmann, 1868;Gibbs, 1902 ) with temperature θ , Boltzmann's constant k B and the energy V (p) = V char W (p) , where V char denotes a characteristic volume used for normalization (since W is an energy density, V char is required to arrive at an energy and may be interpreted as the volume of the material or grid point of interest). The question now is how to choose η( t ) such that we indeed attain thermal equilibrium in the long-term limit as t → ∞ . Let us discretize the governing Eq. (D.1) in time with a constant step size t > 0, so that applying a first-order forward-Euler finite-difference stencil leads to and η t (t) = t η(t) is a short notation for the random fluctuation (whose amplitude needs to be found, so multiplication by the time step does not affect the final result). If we look only at the stochastic contribution to the polarization change, then (D.5) As for a random walk, we require the random noise term to have zero mean and to be uncorrelated, i.e., respectively, with some constant diffusion coefficient D * ≥ 0 that captures the (yet to be determined) noise amplitude. After n time steps, the random noise has altered the solution by (D.7) By exploiting the uncorrelated nature of noise from distinct time steps, we conclude that the mean squared difference between the initial and final position over the above n steps is (D.8) Assuming an unbiased random walk, this implies that, with the total elapsed time n t , [ p(t + n t ) − p(t )] 2 = 2 D n t, (D.9) so that a comparison of (D.8) and (D.9) yields η 2 t (t) = 2 D t.
(D.10) Such a scenario is achieved by choosing a Gaussian noise of average 0 and amplitude 2 D t , whose probability distribution is ( Gauss, 1809;Laplace, 1774 ) ρ(η t ) = Next, we consider the full governing Eq. (D.3) , including the nonconvex potential, to identify the unknown constant D . We start with p(t + t ) = p(t ) + v (p) t + η t (t) . (D.12) Using a generalized version of the time evolution of the probability distribution, we may write ρ(p, t + t ) = ∞ −∞ P t (p, q ) ρ(q, t )d q, (D.13) where P t ( p, q ) is the probability that the particle moves from q at time t to the position p at time t + t, and we integrate over all possible positions q . (D.13) is known as the Chapman-Kolmogorov equation ( Kampen, 2007 ). Moving from q to p in our scenario implies that p = q + v (q ) t + η t (t) , cf. (D.12) . Simply speaking, the random noise term has the right magnitude to help move the particle from q to p (while the potential is also acting). The probability that the noise has exactly a magnitude of η t = p − q − v (q ) t is defined by the Gaussian distribution (D.11) : (D.14) A Kramers-Moyal expansion ( Kramers, 1940;Moyal, 1949 ) of the master equation ( Kampen, 2007 ) is used to finally derive the well-known Fokker-Planck equation ( Fokker, 1914;Kolmogoroff, 1931;Planck, 1917 ) where we inserted the definition of ν from (D.4) .

(D.17)
It is straightforward to verify by substitution that the following presents a solution: ρ eq (p) = Therefore, in the limit of long times, the probability distribution approaches the steady-state Boltzmann distribution (D.20) if we choose the random noise (using (D.10) ) according to η t (t) = 0 and η t (t) η t (t ) = 2 k B θ μV char t δ(t − t ) .

(D.21)
This defines the temperature-dependent random noise term in 1D. Since the random noise term must satisfy these relations in each direction (and at every point inside the RVE), the generalization to 3D leads directly to the relations presented in Section 3.2 .

(E.2)
Both constants are zero because of stress-free boundary conditions. If we assume that the domain wall moves at a constant speed v , then we may express all fields of interest as a traveling wave solutions depending on the coordinate ξ = x · e 1 − v t.
(E.6) Fig. E1. Schematic representation of a ferroelectric single-crystal with a 180 • domain wall being subject to an electric field in the e 2 -direction.
By substitution into (E.5) , we finally obtain where f represents the Eshelby traction acting on the domain wall and, for a constant domain wall shape, C = ∞ −∞ μp 2 2 ,ξ d ξ is a constant characterizing the dissipation in the moving domain wall and depending on the shape of the domain wall. For a constant applied electric field (at a constant temperature), we have f = W = const . for all times, so that we conclude C v = f = const. , (E.8) so we expect a constant domain wall speed v emerging from the linear gradient descent kinetics (as one may expect).
It is important to note that this holds true only for electric fields below the coercive field, so ferroelectric switching is accommodated purely by the motion of domain walls. We point out further that (E.7) hints at the influence of the evolution equation: if one replaces the Allen-Cahn evolution law by, e.g., a more complex dissipative mechanism (e.g., defining a power-law type dissipation potential as found in empirical macroscale models for ferroelectrics ( Miehe and Rosato, 2011 )), then this transforms the left-hand side of (E.7) accordingly and one can conclude the resulting relation between the domain wall speed v and the driving force f .