Guided search for desired functional responses via Bayesian optimization of generative model: Hysteresis loop shape engineering in ferroelectrics
11 Notice: This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).
Guided search for desired functional responses via Bayesian optimization of generative model: hysteresis loop shape engineering in ferroelectrics
Sergei V. Kalinin, Maxim Ziatdinov , and Rama K. Vasudevan The Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37831 The Computational Sciences and Engineering Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831 Advances in theoretical modelling across multiple disciplines have yielded generative models capable of high veracity in predicting macroscopic functional responses of materials emerging as a result of complex non-local interactions. Correspondingly, of interest is the inverse problem of finding the model parameter that will yield desired macroscopic responses, such as stress-strain curves, ferroelectric hysteresis loops, etc. Here we suggest and implement a Gaussian Process based methods that allows to effectively sample the degenerate parameter space of a complex non-local model to output regions of parameter space which yield desired functionalities. We discuss the specific adaptation of the acquisition function and sampling function to make the process efficient and balance the efficient exploration of parameter space for multiple possible minima and exploitation to densely sample the regions of interest where target behaviors are optimized. This approach is illustrated via the hysteresis loop engineering in ferroelectric materials, but can be adapted to other functionalities and generative models. The code is open sourced and available at [github.com/ramav87/Ferrosim]. [email protected]
Introduction
Since the discovery of BaTiO solid solutions in the late 1940s, ferroelectric materials have become one of the central topics in condensed matter physics and materials science alike. Devices such as SONARs and later piezoelectric and optoelectronic ceramics have attracted the strong interest from application community, while the unique physics of this materials have stimulated intensive fundamental research and theoretical studies of these material. Initially, fundamental research was focused on the properties of single crystals of simple ferroelectrics, whereas many applications were preponderantly based on the ceramic materials. With the development of the field, the attention of scientific community has been shifted towards the materials systems such as disordered ferroelectrics, morphotropic phase boundary materials, and ferroelectric relaxors.
These materials generally offer giant dielectric and electromechanical responses, both of interest to applications and necessitating understanding the fundamental origins of these phenomena. After multiple studies, it was established that a common feature of these systems is the presence of nanoscale polarization inhomogeneities emerging due to the interplay between polarization instabilities, mesoscale elastic and electric depolarization fields, and structural defects. Correspondingly, the materials response to even small perturbations can result in strong changes in the inhomogeneous ground state of the system, giving rise to a rich spectra of phenomena such as broad relaxation time distributions, strong frequency dispersion of responses, giant responses, smeared phase transitions, etc. These behaviors are not limited to disordered ferroelectrics and similar behaviors are observed in other materials including phase separated manganites, spin and cluster glasses, etc.
However, these developments have further brought forth the challenge of understanding these materials. In classical ferroelectrics that are homogeneous on the length scales between unit cell and domain size, the functionalities can be readily described using mesoscopic Ginzburg Landau Devonshire (GLD) type theories with the parameters readily available from macroscopic measurements.
More recently, these can also be obtained from the density functional theory (DFT) models. The gradient and interfacial terms necessary for any phase-field type model are often postulated, or can be derived from atomistic imaging.
Jointly, this combination of DFT and GLD allows comprehensive description of ferroelectric materials and microstructures on all length scales. With these, materials design and optimization on the unit cell level can be approached via a theory guided search through the chemical space of existing materials.
The situation is fundamentally different for the disordered ferroelectrics. Here, the functionalities emerge as a result of collective phenomena on multiple length scales, including processes such as polarization rotations and ferroelastic domain wall motion. Often, these mechanisms are defined only via macroscopic descriptors, whereas associated local mechanisms are virtually unknown. The existence of order parameter that can describe these behaviors is heavily debated, as is the case for other disordered systems. Finally, the applicability of DFT type models for the description of the stochastic mesoscopic systems is naturally limited. These considerations have stimulated the intensive development of the mesoscopic lattice type models, where the functionality of the material is represented via collection of spin-like variable on the periodic lattice. The properties (scalar, vectors) and corresponding dynamics are introduced, and time dynamics of the system for various field, temperature, and time histories is explored using classical Monte-Carlo methods.
While these models provided the certain insight into the evolution and properties of the ferroelectric and relaxor state, the nature of the relaxor states, appropriate reduced descriptors, and especially pathways to tailor macroscopic functionalities via microstructure have not been explored. Here, we explore the targeted optimization of the relaxor functionalities via microstructural engineering. As the descriptors for relaxor behavior, we choose the strength and concentration of the defects and the parameters of the embedding GLD functional. In this description, the global polarization behavior emerges as a collective response of polarization fields interacting with the defects. We show that the Gaussian Process optimization can be used to explore the phase space of this system and discover the functionalities of interest. This approach allows to both explore the general parameter space and localize target areas of interest (of which there can be more than one), as demonstrated here for the macroscopic observable of microscopic model. We note that while here this analysis is performed for random defect distribution, this analysis can be performed for microstructural models based on interacting defects or sampled from generative models derived from experimental observations, for example from direct imaging studies, or from atom probe tomography, pair distribution functions from scattering, etc.
Results and Discussion
As a model of ferroelectric material, we chose the lattice model introduced by Ricinschi et al. Here, the polarization at each lattice site is represented as a continuous scalar variable, p , as shown in Fig. 1 (a). The local energy as a function of this variable follow the standard GLD form, πΉ (cid:3039)(cid:3042)(cid:3030) = (πΌ 2β )π (cid:2870) + (π½ 4β )π (cid:2872) β πΈ (cid:3039)(cid:3042)(cid:3030) π , where Ξ± = πΌ (cid:2868) (π β π (cid:3004) ) and Ξ² are GLD coefο¬cients and π (cid:3004) is Curie temperature. The term βπΈ (cid:3039)(cid:3042)(cid:3030) π describes the usual coupling between polarization and the local electric field. The collective effects in this model are introduced via the gradient like interactions, πΉ = β (cid:3435)πΉ (cid:3039)(cid:3042)(cid:3030) + πΎ β (π (cid:3036),(cid:3037) β π (cid:3036)(cid:2878)(cid:3038),(cid:3037)(cid:2878)(cid:3039) ) (cid:2870)(cid:3038),(cid:3039) (cid:3439) (cid:3015)(cid:3036),(cid:3037) , (1) where i , j = 1, .. N are the lattice sites, K is the coupling constant, and the sum over k, l is taken over chosen neighborhood for each site. The disorder in the model can be introduced in multiple ways, e.g. via the position dependence of coupling constant K (bond disorder), or presence of the local electric fields πΈ (cid:3031) (π, π) (field disorder). The depolarization field effects are introduced via spatially uniform depolarization field πΈ (cid:3031)(cid:3032)(cid:3043) = βπΌ (cid:3031)(cid:3032)(cid:3043) β©πβͺ , where πΌ (cid:3031)(cid:3032)(cid:3043) is the depolarization factor and β©πβͺ is averaged polarization. In realistic systems, πΌ (cid:3031)(cid:3032)(cid:3043) is determined by the geometry of the system and efficiency of screening; here we consider it to be a control parameter. The local field acting on each lattice site is then the sum of external, depolarization, and disorder fields, πΈ (cid:3039)(cid:3042)(cid:3030) = πΈ (cid:3032)(cid:3051)(cid:3047) + πΈ (cid:3031)(cid:3032)(cid:3043) + πΈ (cid:3031) (π, π) . Note that the separation of the acting field into the local on-site component and averaged off-site component is the simplest possible mechanism accounting for the physics of the problem. Finally, the dynamics of the local spins is given by the classical Landau-Khalatnikov equation, πΎππ (cid:3036)(cid:3037) ππ‘β = βππΉ ππ (cid:3036)(cid:3037) β . The output of the models can be calculated in the usual fashion and can represent the average polarization β©πβͺ , its evolution as a function of external control stimuli, or some descriptors of microscopic states. Figure 1. (a) Schematic representation of the lattice GLD model and (b) free energy of each spin with (red) and without (blue) random field disorder. Here the model is realized on a square lattice of size N * N , with the starting polarization at each site equal to the nominal remnant state of the system, and given by the solution to (1) in the absence of field, i.e. π (cid:3045) = β(cid:3493)β πΌ π½β . An applied electric field, in this case a simple sine wave with prescribed frequency and amplitude is used. The simulation begins by selecting a lattice site, and then computing the gradient at time t n : (cid:3031)(cid:3043) (cid:3284)(cid:3285) (cid:3031)(cid:3047) = βπΎ (cid:2879)(cid:2869) (cid:3435)π½π (cid:3036)(cid:3037)(cid:2871) + πΌπ (cid:3036)(cid:3037) + πΎ β (π (cid:3036)(cid:3037) β π (cid:3038)(cid:3039) ) (cid:3038),(cid:3039) β πΈ (cid:3039)(cid:3042)(cid:3030) (cid:3439) (2) The polarization at the chosen site is updated in the direction of the gradient, i.e. π (cid:3036)(cid:3037) (π‘ (cid:3041)(cid:2878)(cid:2869) ) = π (cid:3036)(cid:3037) (π‘ (cid:3041) ) + βπ‘ β (cid:3031)(cid:3043) (cid:3284)(cid:3285) (cid:3031)(cid:3047) . Once the polarization is updated across all lattice sites for t n , the simulation progresses to the next time step where the same loop resumes. The output of the simulation is the polarization evolution of the system, discretized in space and time, i.e. a matrix of size ( t , n , n ) where t is the number of time steps of the simulation, and n is the simulation size (assuming a square lattice of size n x n ). It should be noted that this lattice model represents the intermediate case between the fixed spin lattice models such as Ising or Heisenberg models, and the classical phase-field approach. Compared to the Ising model, the spin variable now is continuous and follows the classical Khalatnikov equation for order parameter. At the same time, the depolarization fields are assumed uniform, unlike the PFM models. The latter approximations allows avoiding the use of semi-implicit spectral methods necessary to recalculate the depolarization field distributions at each step, severely reducing the calculation time. That said, the approach developed here is universal and can be applied for any generative model, as limited by available computational capabilities. Figure 2. (a) Sequence of the polarization states of the 10x10 lattice during the sinusoidal sweep of the electric fields. (b) Evolution of the hysteresis loops as a function of the depolarization factor. (c) Evolution of the hysteresis loops as a function of the defect concentration. To illustrate the behavior of the model, shown in Fig. 2(a) is the evolution of the microscopic states of the system during the polarization switching induced by external electric field. Here, the domain nucleation, formation of domain walls, and interaction of domain walls with the random field defects are clearly visible. The evolution of the average polarization as a function of the field is show in Fig. 2 (b) for different depolarization factors. Here, classical well-saturated hysteresis loops for zero depolarization are obtained. On increasing the depolarization factors, the loops become more constrained and ultimately collapse to a straight line, as expected for a ferroelectric material. Similarly, evolution of the hysteresis loop as a function of the random field strength is shown in Fig. 2 (c), exhibiting successful narrowing of the hysteresis loop and pinning to the preexisting defects. Similarly, effects of other control variables including the concentration of the defects, strength of local coupling, or GLD coefficients can be explored. With this model in hand, we aim to explore which combinations of parameters can be used to yield the desired functionalities or microstructures. For ferroelectric material, desired functionality can include certain characteristics of the hysteresis loop. These can include classical descriptors such as hysteresis loop area, remanent polarization, imprint, or asymmetry. Alternatively, ad hoc descriptors can be introduced to describe the deviation of the loop from square, maximum slope, etc. as required by specific applications. Similarly, descriptors can be constructed for the more complex field histories, for example specific Preisach density characteristics. Finally, descriptors can be constructed based on spatial variability of polarization field and its field history dependence, e.g. variability of polarization field for specific voltage and field history, critical bias for domain formation, etc. It is important to note that the system behavior and associated descriptors can be sensitive to the details of the microstructure, i.e. the exact positions of the defects in this case. It is generally expected that macroscopic descriptors such as loop shape will be robust with respect to exact microscopic defect configurations and be determined by the average concentrations and defect strength, whereas specific polarization field distributions are not. These considerations will be important for future discussion of materials design.
Figure 3.
Evolution of the hysteresis loop parameters including (a) total polarization (integral of the upper branch of the loop), (b) remanent polarization on the upper branch of the loop, (c) area under the loop, and (d) switchable polarization with the number of defects ( n d , horizontal axis) and random field strength ( E d , vertical axis), both in the (0,100) interval. In this setting, the material design for specific applications becomes the inverse problem of finding the model parameters that maximize the required functionality, i.e. specific hysteresis loop or robust microscopic descriptor. The direct approach for solving such problem can be based on direct grid search in the parameter space, as illustrated in Fig. 3 for dependence of hysteresis loop parameters on the defect strength and concentration. From the graphs, it is obvious that here with a good degree of approximation the required properties will be uniform along the lines n d E d = const , as can be surmised from the fact that in this case domain formation is impeded and hence total defect field emerges as a control variable. This observation, along with the available control variables such as dopant solubility and tendency for clustering, can be used as a basis for materials selection. Figure 4.
The Bayesian optimization is used to discover the region of the parameter space of the non-local model that gives rise to the target functional behavior. Here, (a) the numerical model is evaluated to give rise to (b) functionality of interest. Here, we choose the lattice ferroelectric model and ferroelectric hysteresis loop, but these can be arbitrary. (c) Based on the target functional, we evaluate a single scalar parameter that represents the βgoodnessβ of the system for applications. (d) Bayesian optimization is used to efficiently explore parameter space for the model balancing the exploration and exploitation, i.e. targeting the optimum value of the target function. However, exploration of the parameter space of the system targeting desired functionalities is a highly non-trivial problem. The dimensionality of parameter space even for the simplified model used here precludes simple grid search approach. At the same time, the methods based on the gradient descent tend to behave poorly when the basins of attraction have complex shape and for the case of the behaviors in Fig. 3 will yield a single combination of n d and E d , while the presence of other solutions tend to remain undiscovered. Finally, if the parameter space contains multiple optima corresponding to disparate degenerate low-dimensional manifolds in high dimensional parameter spaces, many of these will remain undiscovered. Here, we explore the use of Bayesian optimization (BO) for the exploration of these complex parameter spaces to search for target functionalities (Fig. 4). The BO is underpinned by the Gaussian Process (GP) regression methods. Below, we briefly introduce the salient aspects of GP and describe how to build the BO search on top of it. The GP refers to a general problem of learning a function, f , from given set of observations D = ( x , y ), . . .( x N , y N )}. It is assumed that observations are noise-corrupted values of the function, y = f ( x ) + ο₯ , where ο₯ is Gaussian observation noise, whereas the arguments are known exactly. The learning is performed via Bayesian inference in a function space. The key assumption in this analysis is that the function f has a prior distribution f ~ π’π«(0, πΎ (cid:3033) (π₯, π₯β²) , where πΎ (cid:3033) is a covariance function (kernel). The kernel function defines the strength and functional form of correlations between the values of the function across the parameter space. The functional form of the kernel is postulated in the beginning of the regression process. In the cases when physical guidance is absent, the typical choice can be the Gaussian (or Radial Basis Function, RBF) kernel. During the GP regression, the expected value of the function, corresponding uncertainties, and kernel hyperparameters are optimized simultaneously. The output of the GP process is then the predicted data set and uncertainty maps representing the quality of prediction. Additionally, kernel parameters can yield the information on the physics of the system. The GP can be further used for the targeted exploration of the parameter space. Here, the process starts from a small number of seed points that serve to define initial guesses on kernel hyperparameter, function values, and uncertainties. With these, the subsequent points in the parameter space can be selected. In a purely exploratory strategies, the points are selected such as to minimize the uncertainty in the system. In the exploration-exploitation strategies, the acquisition function that includes both the expectation values of the parameter and uncertainty is formed and the process is driven to minimize this acquisition function. The choice of the acquisition function for each specific problem requires tuning to avoid the process to be stuck in a certain region of the parameter space. Thus, GP allows to both explore the general parameter space, and localize target areas of interest (of which there can be more than one). Here we implemented the BO with GP based on the GPim package. We used the RBF or Matern kernels, defined as π (cid:3019)(cid:3003)(cid:3007) (π₯ (cid:2869) , π₯ (cid:2870) ) = π (cid:2870) exp (β0.5 Γ |(cid:3051) (cid:3117) (cid:2879)(cid:3051) (cid:3118) | (cid:3118) (cid:3039) (cid:3118) ) (2) and π (cid:3014)(cid:3028)(cid:3047)(cid:3032)(cid:3045)(cid:3041) (π₯ (cid:2869) , π₯ (cid:2870) ) = π (cid:2870) exp (cid:4672)ββ5 Γ |(cid:3051) (cid:3117) (cid:2879)(cid:3051) (cid:3118) |(cid:3039) (cid:4673) (cid:4672)1 + β5 Γ |(cid:3051) (cid:3117) (cid:2879)(cid:3051) (cid:3118) |(cid:3039) + (cid:2873)(cid:2871) Γ |(cid:3051) (cid:3117) (cid:2879)(cid:3051) (cid:3118) | (cid:3118) (cid:3039) (cid:3118) (cid:4673) , (3) where l and ο³ are kernel length scale and variance, respectively, which are learned from the data by maximizing the log-marginal likelihood. Here, we explored the effect of the defect concentration, defect strength, and depolarization field on the hysteresis loop shape. The 2D or 3D parameter space is defined as dense uniform grid of possible parameter values. Note that the function values are not evaluated initially and the grid is used only as discretization of the parameter space. We subsequently implement the Bayesian optimization strategy exploring the parameter space based on the maximum uncertainty (exploration), searching for target function values (exploitation), or both. To avoid the exploration of only narrow region of parameter space once the target function value is pursued in exploitation strategy, we have introduced the approach where the next measurement point cannot be chosen closer then given distance l to the previous point. The distance can be set as equal to the kernel length scale or defined ad hoc . Subsequent experimentation has suggested that optimal results can be obtained if the short-term memory is introduced, namely that location at time t is chosen no closer than l to the location at step t -1, al at step t -2, a l at step t -3 and so on. Here, we have chosen the length of the memory to be equal to 10 and memory decay coefficient a = 0.8; however other values are possible. Finally, to combine the exploration of the broad parameter space and maximization of required behavior, we have introduced combined strategy where the acquisition function was switched every ten steps from purely exploratory to the exploitation. In the exploratory regime, the points in the parameter space are chosen based on the maximum uncertainty. In the exploitation step, the algorithm maximized the target function, defined here as πΏπ = ππ₯π(β(ππππ β π‘πππππ‘) (cid:2870) ) , (4) where mean is the expected value of the function from GP regression, and target is the target value of specific descriptor. The effect of the target function is to focus exploration in the locations where the expected value of the function is closest to target. Figure 5.
GP-based optimization of the hysteresis loop asymmetry targeting the value 0.5. To illustrate the Bayesian optimization of the hysteresis loop shape, shown in Fig. 5 is the evolution of the search process in the two-dimensional parameter space for n d in the range (0,80) and E d in the range (0,80). This range covers the complete evolution of the hysteresis loop from completely open for ( n d , E d ) = (0,0) to collapsed at (80,80). As a target function, we have somewhat arbitrarily chosen the loop offset being 0.5, corresponding to the remanent polarization on the lower branch of the loop being zero. At the initial stages of search (10 steps), the measurement points are randomly distributed in the parameter space and the corresponding uncertainty map determined by the convolution of measurement points and kernel function shows large variability. Upon sequential evolution (step 50), the algorithm identifies the manifold on the parameter space where the function adopts the target value. From this moment, the exploitation step preponderantly focuses on the points in the vicinity of this manifold, whereas the exploration steps continue sampling the parameter space (steps 150 and 500). The kernel hyperparameters in this process converge relatively rapidly and remain stable during the optimization, reflecting the monotonic character of response function. Figure 6.
Bayesian optimization of the area under the loop as a function of random field and strength and defect concentration. Shown are exploration histories, prediction and uncertainty maps after 100 iteration steps for (a) target value of 0.7 and (b) 0.9. This process is further illustrated in Fig. 6 showing the exploration histories, predictions, and uncertainty maps for the hysteresis loop area for different target values. Here, when the target value is chosen to be 0.7, the search points are concentrated in the top right corner of the parameter space, leading to the reduced uncertainty. The change of the target value to 0.9 concentrates the search in the vicinity of the maximum. Remarkably, the prediction maps converge very rapidly and are very close to each other. The same optimization process can be extended to 3- and 4D space, including e.g. depolarization field, coupling strength, or field limits. However, visualization of the resulting dependencies in this case is more complicated, and hence is not shown. Conclusions
To summarize, here we suggest and implement a Gaussian Process based methods to effectively explore and exploit the multidimensional parameter space of the complex non-local models to discover the regions of interest, in which the desired response adopts maximal value, i.e. yields desired functionalities. We discuss the specific adaptation of the acquisition function and sampling function to make the process efficient and balance the efficient exploration of parameter space for multiple possible minima and exploitation to densely sample the regions of interest where target behaviors are optimized. This approach is illustrated via the hysteresis loop engineering in ferroelectric materials using continuous lattice model. However, it can be adapted in a straightforward manner to other functionalities and generative models, complementing previously developed stochastic optimization approaches.
The interactive notebook that allows reproducing the manuscript results is freely available at https://git.io/Jftuy.
Data Availability:
The data that support the findings of this study are openly available at [github.com/ramav87/Ferrosim].
Acknowledgements:
This effort (Gaussian Process) is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), Materials Sciences and Engineering Division (S.V.K., R.K.V.) and was performed and partially supported (M.Z.) at the Oak Ridge National Laboratoryβs Center for Nanophase Materials Sciences (CNMS), a U.S. Department of Energy, Office of Science User Facility. References
1 Setter, N. et al.
Ferroelectric thin films: Review of materials, properties, and applications.
J. Appl. Phys. , doi:051606 10.1063/1.2336999 (2006). 2 Muralt, P., Polcawich, R. G. & Trolier-McKinstry, S. Piezoelectric Thin Films for Sensors, Actuators, and Energy Harvesting.
Mrs Bulletin , 658-664, doi:10.1557/mrs2009.177 (2009). 3 Damjanovic, D. Ferroelectric, dielectric and piezoelectric properties of ferroelectric thin films and ceramics. Reports on Progress in Physics , 1267-1324, doi:10.1088/0034-4885/61/9/002 (1998). 4 Shaw, T. M., Trolier-McKinstry, S. & McIntyre, P. C. The properties of ferroelectric films at small dimensions. Annual Review of Materials Science , 263-298, doi:10.1146/annurev.matsci.30.1.263 (2000). 5 A.K. Tagantsev, L. E. C., and J. Fousek. Domains in Ferroic Crystals and Thin Films . (Springer, 2010). 6 Rodel, J. et al.
Perspective on the Development of Lead-free Piezoceramics.
Journal of the American Ceramic Society , 1153-1177, doi:10.1111/j.1551-2916.2009.03061.x (2009). 7 Viehland, D., Li, J. F., Dai, X. H. & Xu, Z. Structural and property studies of high Zr-content lead zirconate titanate. Journal of Physics and Chemistry of Solids , 1545-1554, doi:10.1016/0022-3697(96)00025-x (1996). 8 Dai, X. H., Xu, Z. & Viehland, D. Disordered oxygen octahedral rotations and glasslike polarization characteristics in rhombohedral lead zirconate titanate. J. Appl. Phys. , 1919-1921, doi:10.1063/1.363009 (1996). 9 Teslic, S., Egami, T. & Viehland, D. in Applications of Synchrotron Radiation Techniques to Materials Science Iii
Vol. 437
Materials Research Society Symposium Proceedings (eds L. J. Terminello, S. M. Mini, H. Ade, & D. L. Perry) 169-174 (Materials Research Society, 1996). 10 Grinberg, I., Suchomel, M. R., Davies, P. K. & Rappe, A. M. Predicting morphotropic phase boundary locations and transition temperatures in Pb- and Bi-based perovskite solid solutions from crystal chemical data and first-principles calculations.
J. Appl. Phys. , doi:094111 10.1063/1.2128049 (2005). 11 Woodward, D. I., Knudsen, J. & Reaney, I. M. Review of crystal and domain structures in the PbZrxTi1-xO3 solid solution. Phys. Rev. B , doi:104110 10.1103/PhysRevB.72.104110 (2005). 12 Xu, Z., Kim, M. C., Li, J. F. & Viehland, D. Observation of a sequence of domain-like states with increasing disorder in ferroelectrics. Philosophical Magazine a-Physics of Condensed Matter Structure Defects and Mechanical Properties , 395-406 (1996). 13 Dai, X. H., Xu, Z. & Viehland, D. Normal to relaxor ferroelectric transformations in lanthanum-modified tetragonal-structured lead zirconate titanate ceramics. J. Appl. Phys. , 1021-1026, doi:10.1063/1.360889 (1996). 14 Xu, Z., Dai, X. H., Li, J. F. & Viehland, D. Coexistence of incommensurate antiferroelectric and relaxorlike ferroelectric orderings in high Zr-content La-modified lead zirconate titanate ceramics. Appl. Phys. Lett. , 1628-1630, doi:10.1063/1.115673 (1996). 15 Dai, X. H., Xu, Z. & Viehland, D. Long-time relaxation from relaxor to normal ferroelectric states in Pb0.91La0.06(Zr0.65Ti0.35)O-3. Journal of the American Ceramic Society , 1957-1960, doi:10.1111/j.1151-2916.1996.tb08019.x (1996). 16 Glazounov, A. E., Tagantsev, A. K. & Bell, A. J. Evidence for domain-type dynamics in the ergodic phase of the PbMg1/3Nb2/3O3 relaxor ferroelectric. Phys. Rev. B , 11281-11284, doi:10.1103/PhysRevB.53.11281 (1996). 7 17 Glinchuk, M. D. & Stephanovich, V. A. Dynamic properties of relaxor ferroelectrics. J. Appl. Phys. , 1722-1726, doi:10.1063/1.369316 (1999). 18 Tagantsev, A. K. & Glazounov, A. E. Does freezing in PbMg1/3Nb2/3O3 relaxor manifest itself in nonlinear dielectric susceptibility? Appl. Phys. Lett. , 1910-1912, doi:10.1063/1.123710 (1999). 19 Glinchuk, M. D. & Stephanovich, V. A. Theory of the nonlinear susceptibility of relaxor ferroelectrics. Journal of Physics-Condensed Matter , 11081-11094, doi:10.1088/0953-8984/10/48/027 (1998). 20 Glazounov, A. E. & Tagantsev, A. K. Direct evidence for Vogel-Fulcher freezing in relaxor ferroelectrics. Appl. Phys. Lett. , 856-858, doi:10.1063/1.122024 (1998). 21 Bdikin, I. K., Shvartsman, V. V. & Kholkin, A. L. Nanoscale domains and local piezoelectric hysteresis in Pb(Zn1/3Nb2/3)O-3-4.5%PbTIO3 single crystals. Appl. Phys. Lett. , 4232-4234, doi:10.1063/1.1627476 (2003). 22 Levy, J., Hubert, C. & Trivelli, A. Ferroelectric polarization imaging using apertureless near-field scanning optical microscopy. Journal of Chemical Physics , 7848-7855, doi:10.1063/1.481389 (2000). 23 Pan, X. Q., Kaplan, W. D., Ruhle, M. & Newnham, R. E. Quantitative comparison of transmission electron microscopy techniques for the study of localized ordering on a nanoscale.
Journal of the American Ceramic Society , 597-605 (1998). 24 Iwata, M. et al. Domain wall observation and dielectric anisotropy in PZN-PT by SPM.
Materials Science and Engineering B-Solid State Materials for Advanced Technology , 88-90, doi:10.1016/j.mseb.2005.02.004 (2005). 25 Dagotto, E. Complexity in strongly correlated electronic systems.
Science , 257-262, doi:10.1126/science.1107559 (2005). 26 Dagotto, E., Hotta, T. & Moreo, A. Colossal magnetoresistant materials: The key role of phase separation.
Physics Reports-Review Section of Physics Letters , 1-153, doi:10.1016/s0370-1573(00)00121-6 (2001). 27 Tao, J. et al.
Direct Imaging of Nanoscale Phase Separation in La(0.55)Ca(0.45)MnO(3): Relationship to Colossal Magnetoresistance.
Physical Review Letters , doi:097202 10.1103/PhysRevLett.103.097202 (2009). 28 Fath, M. et al.
Spatially inhomogeneous metal-insulator transition in doped manganites.
Science , 1540-1542, doi:10.1126/science.285.5433.1540 (1999). 29 Tokura, Y. Critical features of colossal magnetoresistive manganites.
Reports on Progress in Physics , 797-851, doi:10.1088/0034-4885/69/3/r06 (2006). 30 Wei, J. & Natelson, D. Nanostructure studies of strongly correlated materials. Nanoscale , 3509-3521, doi:10.1039/c1nr10457h (2011). 31 Cao, W. W. & Barsch, G. R. LANDAU-GINZBURG MODEL OF INTERPHASE BOUNDARIES IN IMPROPER FERROELASTIC PEROVSKITES OF D(18)4H SYMMETRY. Phys. Rev. B , 4334-4348, doi:10.1103/PhysRevB.41.4334 (1990). 32 Haun, M. J., Furman, E., Jang, S. J. & Cross, L. E. THERMODYNAMIC THEORY OF THE LEAD ZIRCONATE-TITANATE SOLID-SOLUTION SYSTEM, .1. PHENOMENOLOGY. Ferroelectrics , 13-25, doi:10.1080/00150198908221436 (1989). 33 Haun, M. J. et al. THERMODYNAMIC THEORY OF PBZRO3.
J. Appl. Phys. , 3173-3180, doi:10.1063/1.342668 (1989). 34 Haun, M. J., Furman, E., McKinstry, H. A. & Cross, L. E. THERMODYNAMIC THEORY OF THE LEAD ZIRCONATE-TITANATE SOLID-SOLUTION SYSTEM, .2. TRICRITICAL BEHAVIOR. Ferroelectrics , 27-44, doi:10.1080/00150198908221437 (1989). 8 35 Haun, M. J., Zhuang, Z. Q., Furman, E., Jang, S. J. & Cross, L. E. THERMODYNAMIC THEORY OF THE LEAD ZIRCONATE-TITANATE SOLID-SOLUTION SYSTEM, .3. CURIE CONSTANT AND 6TH-ORDER POLARIZATION INTERACTION DIELECTRIC STIFFNESS COEFFICIENTS. Ferroelectrics , 45-54, doi:10.1080/00150198908221438 (1989). 36 Haun, M. J., Furman, E., Halemane, T. R. & Cross, L. E. THERMODYNAMIC THEORY OF THE LEAD ZIRCONATE-TITANATE SOLID-SOLUTION SYSTEM, .4. TILTING OF THE OXYGEN OCTAHEDRA. Ferroelectrics , 55-62, doi:10.1080/00150198908221439 (1989). 37 Borisevich, A. Y. et al. Exploring Mesoscopic Physics of Vacancy-Ordered Systems through Atomic Scale Observations of Topological Defects.
Physical Review Letters , doi:065702 10.1103/PhysRevLett.109.065702 (2012). 38 Li, Q. et al.
Quantification of flexoelectricity in PbTiO3/SrTiO3 superlattice polar vortices using machine learning and phase-field modeling.
Nat. Commun. , doi:1468 10.1038/s41467-017-01733-8 (2017). 39 Vugmeister, B. E. & Rabitz, H. Coexistence of the critical slowing down and glassy freezing in relaxor ferroelectrics. Phys. Rev. B , 14448-14453, doi:10.1103/PhysRevB.61.14448 (2000). 40 Juhas, P. et al. Correlations between the structure and dielectric properties of Pb(Sc-2/3(1/3)3(3)W()O)-Pb(Ti/Zr)O relaxors.
Phys. Rev. B , doi:214101 10.1103/PhysRevB.69.214101 (2004). 41 Vugmeister, B. E. Polarization dynamics and formation of polar nanoregions in relaxor ferroelectrics. Phys. Rev. B , doi:174117 10.1103/PhysRevB.73.174117 (2006). 42 Grinberg, I., Juhas, P., Davies, P. K. & Rappe, A. M. Relationship between local structure and relaxor behavior in perovskite oxides. Physical Review Letters , doi:267603 10.1103/PhysRevLett.99.267603 (2007). 43 Takenaka, H., Grinberg, I., Liu, S. & Rappe, A. M. Slush-like polar structures in single-crystal relaxors. Nature , 391-+, doi:10.1038/nature22068 (2017). 44 Ricinschi, D. et al.
Analysis of ferroelectric switching in finite media as a Landau-type phase transition.
Journal of Physics: Condensed Matter , 477-492, doi:10.1088/0953-8984/10/2/026 (1998). 45 Rasmussen, C. E. & Williams, C. K. I. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) . (The MIT Press, 2005). 46 Ziatdinov, M. GPim: Gaussian processes for images and hyperspectral data.
GitHub repository, https://git.io/JvK1s (2020). 47 Jayachandran, K. P., Guedes, J. M. & Rodrigues, H. C. Ferroelectric materials for piezoelectric actuators by optimal design.
Acta Mater. , 3770-3778, doi:10.1016/j.actamat.2011.02.005 (2011). 48 Jayachandran, K. P., Guedes, J. M. & Rodrigues, H. C. Solutions for maximum coupling in multiferroic magnetoelectric composites by material design. Sci Rep , 9, doi:10.1038/s41598-018-22964-9 (2018)., 9, doi:10.1038/s41598-018-22964-9 (2018).