Modeling structure–properties relations in compositionally disordered relaxor dielectrics at the nanoscale

Relaxor ferroelectrics are compositionally disordered perovskite oxides that, unlike ordinary ferroelectric perovskites, exhibit broad or diffuse peaks in the temperature dependence of their dielectric permittivity, with the location of the peak maximum also being strongly dependent on the frequency of the applied AC electric field.^1,2 Due to their excellent electroactive properties, these materials are a subject of research for high efficiency energy storage and conversion applications,^3,4 as well as for integration into tunable microwave devices (e.g., for telecommunications)⁵ or transducers with enhanced piezoelectric properties.^6,7 In particular, the relaxor barium strontium titanate, Ba $1 − x$Sr $x$TiO $3$ (BSTO), system is being investigated as an efficient, lead-free, and inexpensive component of microwave devices with characteristic low loss,⁸ that performs better than other semiconductor candidates in the GHz range.^8–11 

Although relaxor ferroelectrics have been discovered more than 50 years ago, they remain at the focus of intense interest, as possible connections between their local structure (including crystallographic defects, polar domain orientations, polycrystallinity, and other factors) and dielectric properties still require better qualitative and quantitative understanding. A number of different theories attempting to connect relaxor behavior with specific features of the system geometry and composition have been put forward during the past few decades. Earlier models represent relaxors as a linear-dielectric matrix into which nanoscale polar regions are embedded.¹² Another group of models views these materials as comprised only of quasi-randomly oriented nanoscale polar domains oriented, i.e., as a form of a polar glass.¹³ The work by Geneste et al.¹⁴ used a coarse-grained Hamiltonian approach to study the dielectric relaxation in Ba(Zr $0.5$Ti $0.5$⁠)O $3$ across the radio frequency range. Their findings demonstrate that random local electric or elastic fields can lead to collective relaxor-like dielectric relaxation in that frequency range.

In the recent work of Takenaka et al.,^15,16 a so-called polar slush model has been proposed for the canonical (PbMg $1 / 3$Nb $2 / 3$O $3$⁠) $1 − y$–(PbTiO $3$⁠) $y$ relaxor system. In this model, local polarization is arranged in $∼$2–10 nm size domains with no linear-dielectric matrix in between, while mutual orientations of the neighboring domains are observed to form high- or low-angle domain walls.¹⁷ As proposed in Ref. 18, it is the slush-like texture of the polarization field that gives rise to the relaxor behavior in the dielectric properties.

On the other hand, the structural inhomogeneities present in complex ferroelectric systems—with relaxors certainly qualifying as such—that include local variations of the composition and/or polarization (i.e., domains or polar nanoregions), polycrystallinity, defects, etc., are well known to produce deviations from the Debye-like behavior at frequencies in the GHz range.² Such relaxation processes are usually described by a variety of other empirical expressions, e.g., Cole–Cole, Havriliak–Negami, or fractional dynamics,^22–26 some of which are discussed and utilized in the sections below.

In general, these expressions may provide evidence that the system dispersion can no longer be described by a global parameter, such as, e.g., $τ$ in Eq. (1) or a similar one. Instead, a spatial distribution of multiple processes having different characteristic relaxation times could be occurring simultaneously, which makes establishing any empirical quantitative connections more difficult. However, although these formulas are flexible enough to provide good fits of experimental $χ ′(ω), χ ′ ′(ω)$⁠, or $χ ′ ′( χ ′)$ dependencies (the latter one usually being in the form of a Nyquist plot), any relations between their empirical parameters and the specifics of the relaxor system structure and composition are largely unknown. We should point out that at high frequencies (in the THz range), or in the paraelectric state, i.e., at temperatures well above the transition temperature $T C$⁠, Debye-like behavior illustrated by Eq. (1) is recovered.

In this investigation, we utilized time-dependent Landau–Ginzburg (TD-LG) theory combined with finite-element method (FEM)-based simulations to establish and quantify direct connections between variations of local composition in a relaxor system and its dielectric properties. Monocrystalline (no grain boundaries) and monodomain-polarized (no domain walls) BSTO was used as a representative model of a relaxor ferroelectric compound due to availability of a thermodynamic potential that is parameterized for a wide range of temperatures and average Sr concentrations $x$⁠.²⁷ The variation of local composition within the system was introduced by means of a correlated random field $x( r)$ that is characterized by a correlation length $L$⁠, corresponding to an approximate size of a region with the same Sr concentration, and the allowed deviation $δx$ of the local concentration away from the average.

The dielectric response of this representative system, including both susceptibility $χ ′$ and loss $χ ′ ′$⁠, was probed as a function of the driving AC field frequency $ω$⁠, temperature $T$⁠, and composition-field parameters ${x,δx,L}$⁠. We show that the presence of local composition variations in the system—that does not possess any other spatial inhomogeneities, such as crystal grains and/or polar domains—results in the emergence of all the major attributes of relaxor behavior, such as broadening of the peak in the $χ ′(T)$ curve, the dependence of the peak location on $ω$ and a prominent non-Debye character of the system relaxation near the transition temperature $T C$⁠.

Furthermore, the adopted modeling approach allowed us to sample a large number of different $x( r)$ configurations for the chosen $T$ and, at the same time, evaluate the response and relaxation of the entire polar dipole population (in the continuum approximation) within each configuration. The availability of such “fine-grained” information, which is usually not accessible through experimental measurements of dielectric response, enabled us to study the influence of different data averaging schemes on the shape of the obtained system relaxation curves and the fitted dependencies of their parameters on $T$ and $L$⁠. We uncovered that, unexpectedly, the dipole population relaxation behavior within each configuration remains mostly Debye-like near $T C$⁠, and that it is the averaging over multiple configurations with distinct (Debye-like) dielectric responses that recovers the non-Debye relaxation features observed in experiments.

The rest of this paper is organized as follows. In Sec. II, we present the details of the adopted simulation methodologies, including the techniques employed to compute the dielectric permittivity and loss with TD-LG approach, and the statistical schemes used for averaging the obtained data. In Sec. III, we discuss the results of our computational studies. Finally, in Sec. IV, we provide some conclusions and take-home ideas, as well as outline possible future directions for this work.

All of the FEM structural models for the compositionally disordered ferroelectric relaxor system utilized in this investigation were generated as a single quasi-two-dimensional square (periodic) material block $Ω$ that is $N×N$ nm in size. The $X,Y$ axes of the global Cartesian reference system were aligned with the edges of the square material block. Crystallographic axes $a,b$ describing the orientation of perovskite unit cells within the block were chosen to be collinear with the $X,Y$ axes. The crystallographic $c$ axis was considered to be along the homogeneous (⁠ $Z$⁠) direction, which implies that both the composition fields $x( r)$ and the obtained solutions are homogeneous or “perfectly correlated” in that direction. No polycrystallinity, i.e., existence of separate regions with different orientations of the $a,b,c$ axes within the material block, was allowed. Typical material block sizes utilized in computer simulations were $8×8$ or $16×16$ nm, with the square FEM mesh size being 0.5 nm. In order to study the system behavior, these structural models were processed with a three dimensional computational approach (described in Sec. II B), but on a two dimensional (⁠ $X,Y$⁠) grid, which resulted in a considerable reduction of computational efforts and allowed us to sample a large pool of different composition-field configurations. A number of models with much larger block sizes were also generated for illustrative purposes, as shown, e.g., in Fig. 1.

For the stoichiometric Ba $1 − x$Sr $x$TiO $3$ sample with the average Sr concentration $x$⁠, the presence of inhomogeneities in local composition (due to the disordering of Ba $2 +$ and Sr $2 +$ ions, that have relatively close ionic radii, on the A-site of the perovskite structure) was represented by a (position dependent) correlated random field $x( r,L,δx)$⁠. Here, the correlation length $L$ is a parameter controlling the desired size of Sr- or Ba-rich regions, while $δx$ governs the spread of the field values, with the local Sr concentration varying in between $x−δx$ and $x+δx$⁠. All the generated composition fields were subject to periodic boundary conditions (PBCs) in two dimensions. Additional details of the computational approach utilized for the construction of the composition field $x( r,L,δx)$ are presented in the supplementary material.

Since the typical sizes of the FEM structural models utilized in this study were relatively small, for each system state—characterized by the set of parameters ${x,δx,L,T}$—multiple models, referred in what follows as seeds, with different randomly generated composition-field patterns were constructed. Usually, 100+ seeds were used for each temperature point away from $T C$⁠, while 200+ seeds were employed for any temperature points close to $T C$⁠, with simulation results obtained for individual seeds averaged as described in Sec. II C. Figure 1(a) provides an illustration of a periodic composition field $x( r)$ in a $60×60$ nm simulation box, generated for $x=0.4$⁠, $δx=0.1$⁠, and $L=24$ nm. Additional images of representative composition field patterns are shown in Figs. 1 and 2 in the supplementary material.

The values of the composition-dependent parameters for the $f bulk$ term were taken from the work of Huang and collaborators,²⁷ who mapped that dependence to the one on pressure. The A-site substitution of Sr for Ba leads to a lattice contraction,²⁸ which is equivalent to applying hydrostatic pressure to the BaTiO $3$ (BTO) unit cell.²⁹ Huang and co-workers demonstrated that such an approach reproduces the appropriate phase diagrams for the BSTO system. Since composition-dependent gradient energy coefficients are not known for BSTO, we utilized appropriate coefficients for BTO instead.³⁰ We note that it is now possible to determine the gradient coefficients in $f wall$ directly from first-principles calculations.³¹ However, in the case of BSTO, this requires careful treatment of the A-site compositional disorder for a large number of different domain wall configurations (possibly even separating different phases). Since, here, we are considering small isotropic gradients in $P$ due to the compositional field variation in a monodomain sample, we do not expect the specific choice of the gradient terms to influence the results appreciably. Detailed expressions for all of the free-energy density terms in Eq. (2) as well as values of all parameters used in these expressions are provided in the supplementary material.

The system of coupled Eqs. (3) and (4) was solved in two dimensions with the help of a general real-space finite-element approach, implemented in the open-source Ferret package,³⁶ which is based on the Multiphysics Object-Oriented Simulation Environment (MOOSE) framework.³⁷ Both equations were cast into weak-form representations suitable for the application of Galerkin’s FEM. Automatically adjustable time step values were used in the gradient-flow energy-minimization algorithm. The total free energy of the system was considered converged when its relative change between consecutive time steps during the evolution of Eq. (3) was below 0.1 $%$ and the relative residual norm was less than $10 − 8$⁠.

Figure 1(b) presents an example of an equilibrium state of the $P( r)$ field (at zero applied field $E$⁠), corresponding to the composition field $x( r)$ shown in panel (a). The $P( r)$ field has variations in the magnitude of the polarization vector, with the Ba-rich areas (red) being polarized more strongly compared to the Sr-rich areas (blue). Slight bending of the field vectors (accented by the white arrows) is also evident at the boundaries between the regions. This behavior is similar to the features exhibited by the polar slush model introduced by Takenaka and co-workers^15,16 as well as to the experimental observations conducted in a relaxor system by local piezo-force microscopy.¹⁷ 

In order to evaluate the system dielectric response, the TD-LG Eq. (3) was solved for each applied electric field direction $j$ and frequency $ω≡2πf$ to obtain the time-dependent change in the system polarization $( Δ P ) j$⁠. Hence, fixed time step values in the energy-minimization algorithm were scaled with the applied electric field frequency as $∼0.1 f − 1$⁠. A frequency interval of $10 8$– $10 14$ rad/s (⁠ $∼ 10 7$– $10 13$ Hz) was sampled.

We also need to point out that, within the developed methodology presented here, the response of the system polarization to the applied electric field can be easily evaluated on multiple levels. For example, it is possible to fit the polarization variations described by Eq. (6) at individual nodes of the FEM mesh. In this case, the local polarization is $P s( r)+Δ P( r)$⁠, where $P s( r)$ is the spontaneous (zero applied field) value of polarization on that node, which provides a nodal-level estimate for $χ( r)$⁠. On the other hand, both $P s$ and $Δ P$ can be averaged over all the nodes in a particular seed configuration (losing their dependence on the position vector $r$⁠), and then fitted to Eq. (6), thus producing a seed-level estimate for $χ$⁠. As discussed in more detail in Sec. II D, both of these options were utilized in various approaches for averaging the system dielectric response.

Three possible approaches to averaging the computed data to obtain an “integral” dielectric response of the system state ${x,δx,L,T}$ were explored, which, in the order of decreasing data granularity, are referred to as (i) nodal, (ii) seed, and (iii) global averaging schemes.

In the nodal averaging scheme, for a given seed configuration, nodal-level dielectric responses $χ( r)$ were computed and fitted to the chosen empirical formula at each FEM mesh point of the simulation box. To evaluate the “integral” dielectric response of the system, the parameters of the fitted individual $χ( r,ω)$ curves were then averaged over all nodes within each seed configuration and over all seeds generated for the same temperature.

In the seed averaging scheme, for each seed configuration, variations of $P( r)$ were averaged over all the nodes to produce the seed-level estimate for $χ(ω)$ that was then fitted to the chosen empirical formula. On the next step, the parameters of the fitted individual $χ(ω)$ curves were averaged over all seeds at the same temperature to obtain the “integral” dielectric response of the system.

Finally, in the global averaging scheme, the seed-level estimates for $χ(ω)$ were averaged over all seeds at the same temperature to obtain an “integral” $χ(ω)$ dependence, which was then fitted to the chosen empirical formula. Since this scheme simply averages over the dielectric responses of multiple instances of the same system state and then fits the final result to an empirical expression, its output should closely resemble that of actual experiments that measure dielectric response in macroscopic samples. Naturally, the global averaging scheme requires the least amount of computational effort in data processing, whereas the seed and especially the nodal averaging schemes demand much more extensive processing efforts.

A nonlinear least-squares curve-fitting approach was used for curve fitting in all of the mentioned averaging schemes. Specifically, a trust-region reflective method, available in Python 3.9.13,³⁹ was employed, being a robust algorithm for fitting large sparse problems with bounds.^40–43 

Here, we briefly summarize the main approximations adopted in this study. On the structural level, the study is focused on evaluating the influence of compositional disorder on the dielectric response, while ignoring the effects of microstructure (e.g., polycrystallinity), as well as domain wall motion and rearrangements. Adhering to the progression of discussion in Sec. II, the approximations are separated into the three following groups:

• Structural models are realized as periodic monocrystalline material blocks, with variations in local composition represented by the correlated random field (which is one of many possible approaches to model compositional disorder). Before the application of the oscillating electric field, these models are pre-polarized into a monodomain configuration, i.e., they contain no discernible domain wall patterns where the polarization field has large changes in its orientation.

• System energy parameterization utilizes gradient energy coefficients for BTO instead of those for BSTO since the compositional dependence of the latter is unknown. The elastic interactions are not considered to speed up the calculations. Both of these are acceptable approximations for a system with no pronounced domain walls in the limit of weak perturbing electric fields. The domain wall mobility parameter $γ$ in Eq. (3) is, in principle, frequency dependent. However, since this dependence is also unknown, a fixed value of this parameter is used for all frequencies.

• Dielectric response evaluation is conducted in the limit of weak—i.e., being well below the coercive field in magnitude—harmonic applied fields that in turn produce a harmonic response of the polarization field at the same frequency. This approximation breaks down for strong applied fields, as well as at high (THz and above) field frequencies, where the response of the polarization field can become non-harmonic and include multiple components oscillating at different frequencies. The description of such processes within the adopted model can be improved by adding the second derivative of polarization with respect to time to the left hand side of Eq. (3) and by the harmonic analysis of the $Δ P(t)$ profile (efforts to implement these improvements are currently underway). The $χ ′(f)$ and $χ ′ ′(f)$ dependencies presented below include data for $f> 10 12$ Hz, which should be considered as a guide to the eye; both $χ ′$ and $χ ′ ′$ components become very small at these frequencies, compared to their sub-THz values.

In this section, we present and discuss the results of our computational studies of the frequency-dependent dielectric response of the compositionally disordered BSTO system. The evaluation of the system properties was focused on the temperature interval corresponding to the nonpolar cubic to polar tetragonal phase transition, while any subsequent polar-to-polar transformations (usually occurring at temperatures below 200 K) were not considered. Unless otherwise stated, all the information provided here is for a representative system with the $x=0.4$ average Sr concentration (a popular composition in which many BSTO samples are grown^8,44,45) and a local composition spread $δx=0.1$⁠. Results obtained for the longitudinal mode polarization-field excitations, demonstrating a much more pronounced non-Debye behavior, are discussed here, while results for the transverse mode excitations are provided in the supplementary material. Fits of the dielectric response of the system utilizing the HN expression are shown, while some representative results using other formulas are included in the supplementary material.

In Fig. 2, we present a temperature dependence of the system polarization for different Sr concentrations $x$ in the vicinity of the nonpolar cubic to polar tetragonal phase transition in BSTO. Results for both systems with homogeneous composition and for the ${x,δx,L}$ systems with varying local composition are shown, represented by bold and dashed/dotted lines, respectively. For the latter, the system polarization was obtained by computing a seed-averaged value of polarization magnitude $| P|$ and then averaging that value over all available seeds. Usually, at least 25 seeds were used for each $T$ point.

For the systems with homogeneous composition, sharp phase transitions are observed in agreement with the findings of Ref. 27. On the other hand, in the systems with composition variations, dissimilar $P( r)$ textures generated by different seeds possess distinct transition temperatures $T C$⁠. After averaging, this results in a diffuse character of the phase transition, which implies coexistence of tetragonal and cubic phases at the same temperature. The extent of the diffuse region around the sharp-transition temperature $T C$ can be expanded by increasing both $L$ (see yellow lines) and $δx$ (see red lines) at a given concentration $x$⁠.

In Fig. 3, we present the $χ ′(f), χ ′ ′(f)$⁠, and $χ ′ ′( χ ′)$ dependencies for the ${x,δx,L}={0.4,0.1,3}$ system. Curves for multiple temperatures in the vicinity of the nonpolar cubic to polar tetragonal phase transition for $x=0.4$ [panels (a) and (b)] and for $T≫ T C$ [panels (c) and (d)] are shown. The global averaging scheme was utilized to evaluate the dielectric response.

As depicted in Fig. 2, for this system (yellow dotted line) the polarization magnitude changes continuously between approximately 255 and 290 K, becoming zero at higher temperatures. From the data shown in Fig. 3, it is clear that the magnitude of the dielectric response $χ ′(f)$ in the sub-GHz region becomes very large near $T C≃$ 270–275 K, accompanied by the similar increase in the dielectric loss signal $χ ′ ′(f)$⁠. The respective magnitudes of both $χ ′(f)$ and $χ ′ ′(f)$ decrease substantially for $T< T C$ and $T≫ T C$⁠. Furthermore, the shape of the $χ ′ ′( χ ′)$ plot is distorted away from an ideal semicircle in the temperature interval near $T C$ [see panel (b)], indicating a non-Debye relaxation character for the system in the vicinity of the nonpolar-to-polar phase transition. On the other hand, the same plot is almost perfectly semi-circular for $T≫ T C$ [see panel (d)], suggesting Debye relaxation behavior in the system far away from the phase transition. The presented results for the dielectric response strength as a function of frequency and temperature look quantitatively similar to the experimental measurements of Ostapchuk et al.^44,45 that were obtained for thick polycrystalline samples of Ba $0.6$Sr $0.4$TiO $3$⁠.

In Fig. 4, the same dielectric response data is replotted as a function of temperature for a range of different frequencies, stretching approximately from 0.1 to 11 GHz. Results for the ${0.4,0,0}$ system with homogeneous composition and the ${0.4,0.1,3}$ system with varying local composition are shown in panels (a) and (b), respectively. In the vicinity of $T C$⁠, for the same probing-field frequency, sharper peaks with larger values of $χ ′$ are observed for the fixed-composition system, in comparison to the varying-composition one. Furthermore, in the latter system, when the probing-field frequency is increased, the corresponding positions of the peaks of the $χ ′(T)$ curves display a pronounced shift towards higher temperatures. The insert to panel (b) shows the relationship between the field frequency $f$ and the temperature $T m$⁠, marking the $χ ′(T)$ curve peak position. The data extracted from the curves in panel (b) (filled circles) are in excellent agreement with the fit to the Vogel–Fulcher expression $f= f 0exp{− E a/( T m− T f)}$ (dashed line), where $f 0$⁠, $E a$⁠, and $T f$ are constants.⁴⁶ This behavior suggests that the introduction of compositional inhomogeneities into the (otherwise monocrystalline and monodomain polarized) BSTO model system can reproduce some of the relaxor-like features exhibited by real materials.

The dielectric response curves presented in Fig. 3 for the ${0.4,0.1,3}$ system at different temperatures were fitted to the HN expression given in Eq. (9). The temperature dependence of the parameters in the HN expression, computed utilizing the nodal, seed, and global data averaging schemes, is presented in Fig. 5. As shown in panels (a) and (d), respectively, both the dielectric strength $χ HN 0$ and relaxation time $τ HN$ parameters undergo significant changes in the temperature region around $T C$⁠. The evolution of $χ HN 0$ as a function of $T$ is the same for all the three data averaging schemes. In the case of $τ HN$⁠, similar temperature profiles are obtained, with the peak value near $T C$ displaying some dependence on the scheme chosen—with it being the highest (slowest relaxation) in the global scheme and the lowest (fastest relaxation) in the nodal scheme.

On the other hand, the variation of the of asymmetry $σ HN$ and skewness $α HN$ parameters, shown respectively in panels (b) and (c), exhibits strong changes in behavior depending on the chosen data averaging scheme. Specifically, the values of both parameters have either minor or no deviations from unity—indicating a Debye-like relaxation—near $T C$ in the nodal and seed schemes, while in the global averaging scheme, the values of these parameters drop below 0.8, suggesting a non-Debye relaxation behavior.

This surprising discrepancy can be explained by comparing the dielectric responses—in the form of $χ ′ ′( χ ′)$ Nyquist plots—of multiple individual nodes within the same seed configuration, or multiple individual seed configurations at the same temperature point. It turns out that the majority of these responses, whether computed at temperatures near or far from $T C$⁠, display the Debye-like behavior—as manifested by mostly semicircular shapes of the respective Nyquist plots. However, in the vicinity of $T C$⁠, the values of the radii of the individual semicircular plots exhibit large variations. In the global scheme, averaging over a set of semicircular shapes with vastly dissimilar radii produces an integral shape that is distorted away substantially from the perfect semicircle and the consecutive fitting of this shape to the HN (or other) expression yields parameters that point to a deviation from the Debye-like behavior. On the other hand, fitting individual semicircular shapes first and then averaging their fitted parameters—following the nodal or seed data averaging schemes—results an integral shape that is a nearly perfect semicircle (as the averaged $σ HN$ and $α HN$ still remain close to unity).

With the system temperature moving away from $T C$⁠, the spread in the values of the individual response curves radii becomes much more narrow and the global averaging scheme produces a nearly semicircular integral shape as a result. A comparison of individual seed-level dielectric response curves with their integral global-average counterpart in the ${0.4,0.1,3}$ system for a variety of different temperatures near and far away from $T C$⁠, illustrating the discussion above, is shown in Fig. 6 in the supplementary material.

In Fig. 6, we provide a direct comparison between the temperature evolution of HN-fitted parameters of the systems with different average concentrations $x$ and the same correlation length $L$⁠: the ${0.4,0.1,3}$ system (i.e., the same data as the one shown in Fig. 5) and a ${0.5,0.15,3}$ system. In both cases, the global data averaging scheme was used to obtain the fittings. As can be clearly seen from the figure, the thermal evolution curves for all of the fitted parameters look similar for the both systems. The locations of the curve peaks in the ${0.5,0.15,3}$ system are shifted to a temperature that is an appropriate $T C$ for the average Sr concentration $x=0.5$ (cf. Fig. 2). Furthermore, since the latter system has an increased Sr content and, therefore, possesses lower average polarization, the magnitudes of the curve peaks for all the parameters—including those for $σ HN$ and $α HN$ that indicate deviations from the Debye-like behavior—are somewhat less pronounced in comparison to those associated with the more polar ${0.4,0.1,3}$ system.

Finally, we note that the analysis presented in this section could be conducted utilizing an empirical expression other than HN, producing thermal evolution curves for a different set of fitting parameters (which may in turn carry a somewhat different physical meaning). In Fig. 5 in supplementary material, we provide an alternative example of such analysis for the ${0.4,0.1,3}$ system utilizing a combination of the Debye and CC expressions.

In this section, we report the results of a different investigation of the sensitivity of the system dielectric response on its model parameters and applied conditions. The same ${0.4,0.1,L,T}$ system as before was utilized in the following study, however, unlike the setup discussed in Secs. III B and III C, the system temperature was fixed at 265 K (i.e., near but not too close to the $T C$ for the average Sr concentration $x=0.4$⁠), while the correlation length $L$ was allowed to vary between 0.1 and 8 nm. The model with a very small $L$ (0.1 nm) was chosen as a representative for a perfectly intermixed system configuration. In order to accommodate much larger $L$⁠, compared to the previously studied cases, all simulations described here were conducted for material block sizes of $16×16$ nm.

In panels (a) and (b), respectively, of Fig. 7, we present the $χ ′(f), χ ′ ′(f)$⁠, and $χ ′ ′( χ ′)$ dependencies for the ${0.4,0.1,L}$ system at $T=265$ K. Curves for multiple values of the correlation length $L$ are shown in color. The global averaging scheme was utilized to evaluate the dielectric response.

As shown in panel (a) of the figure, the shape of the $χ ′(f)$ and $χ ′ ′(f)$ displays a strong dependence on the value of $L$⁠. In particular, there is a crossover in behavior that occurs between $L=3$ and 4 nm: For $L≤3$ nm, the maximum value of $χ ′$ remains below 2000 and the peak in the $χ ′ ′(f)$ dependence is located around 10 GHz. On the other hand, for $L≥4$ nm, the maximum value of $χ ′$ keeps increasing with $L$⁠, saturating between 4000 and 5000, and the peak in the $χ ′ ′(f)$ dependence gradually shifts to lower frequencies, settling near 0.8–1 GHz. Similar changes are observed in the shape of the $χ ′ ′( χ ′)$ dependence, shown in panel (b): For $L≤3$ nm, these curves look mostly semicircular, while for $L≥4$ nm, they develop strong distortions away from that shape, indicating a switch from Debye-like to non-Debye relaxation behavior in the system.

The data presented in Fig. 7 demonstrate that concentration variations in BSTO that correspond to the region of “non-Debyeness” with $L≥4$ nm could be interesting for possible dielectric applications in the GHz range. In that region, compared to the “well intermixed” one where $L≤3$ nm, we predict a large improvement in dielectric strength—approximately by 3 $×$ or from $∼1500$ to $∼4500$—that is accompanied by a disproportional increase in the dielectric loss—less than 2 $×$ or from $∼800$ to $∼1400$⁠. [Also, cf. the results shown in Figs. 3(a) and 3(b) for the non-Debye transitional behavior of the ${0.4,0.1,3}$ system, where the level of dielectric loss remains relatively high compared to the dielectric strength magnitude.]

In the next step, following the same approach as the one discussed in Sec. III C, the dielectric response curves presented in Fig. 7 were fitted to the HN expression of Eq. (9). The evolution of the HN expression parameters under changing correlation length $L$⁠, computed utilizing the seed and global data averaging schemes, is displayed in Fig. 8. We observe that for $L≤3$ nm, both averaging schemes produce the same results, with $χ HN 0≃1500$⁠, $τ HN$ being slightly above $0.1× 10 − 10$ s, and $σ HN, α HN$ remaining close to unity, which indicates the Debye-like system relaxation behavior.

As $L$ grows beyond 3 nm, we still find reasonable agreement between the results obtained by both averaging schemes for the evolution of $χ HN 0$⁠, which gradually increases to $∼4500$⁠, $τ HN$⁠, which increases to $∼1.0$–1.2 $× 10 − 10$ s, and $α HN$⁠, which drops to $∼0.9$—with the seed averaging scheme predicting an earlier and more aggressive decrease. The main discrepancy between the averaging approaches is in their description of the evolution of the $σ HN$ parameter, which remains close to unity for all values of $L$ in the seed averaging scheme, while in the global averaging scheme, it drops to 0.7 for $L≥4$ nm—suggesting a crossover from Debye-like to non-Debye behavior, in agreement with the data presented in Fig. 7.

In closing, we note that we did not attempt to extend $L$ beyond $∼$10 nm, since such system configurations—on the microstructural level and depending on the specifics of allowed compositional variations between, e.g., Ba-rich and Ba-poor regions—can be considered, more appropriately, as ferroelectric–dielectric composites. We have already investigated (equilibrium) polar properties of such systems with the same approach as the one described in Sec. II B,^36,47–50 while their dielectric response can be probed straightforwardly by utilization of the techniques discussed in Sec. II C.

In this work, we have conducted numerical simulations of the dielectric properties of the quasi-two-dimensional monocrystalline and monodomain-polarized model of the BSTO relaxor-ferroelectric system. The local composition variations within the system were described by a correlated random field with parameters $x$⁠, $L$⁠, and $δx$⁠, representing, respectively, the average Sr concentration, the approximate size of a region with the same concentration, and the deviation of the local concentration away from the average. The system ground-state properties were evaluated with the TD-LG methodology combined with FEM-based simulations. An original technique was developed and implemented in the open-source Ferret package³⁶ within the MOOSE framework³⁷ to probe the frequency-dependent dielectric response of the system under the influence of the driving AC electric field, producing results that are in good quantitative agreement with some of the experimental measurements available for BSTO.⁴⁵ Extensive statistical sampling of random system configurations (or seeds) with different patterns of the local composition field was carried out for the most important simulation parameters, such as $T$ and $L$⁠, with three distinct data averaging schemes considered for processing the results of these simulations to obtain the integral dielectric response of the system. The extracted dielectric response curves were then fitted to empirical expressions describing relaxation processes in the system and the fitting parameter dependencies on $T$ and $L$ were studied.

As a result of this investigation, we have established connections between the system composition-field descriptors $x$⁠, $L$⁠, and $δx$ and properties, such as its average polarization, frequency-dependent dielectric strength $χ ′$ and loss $χ ′ ′$⁠, for a wide range of temperatures, including those near and far from the $T C$⁠. Specifically, it was determined that the presence of position-dependent compositional variations in the (otherwise monocrystalline and monodomain-polarized) structural model reproduces a number of characteristic features of relaxor systems, such as (i) the diffuse nature of the nonpolar to polar phase transition, (ii) broadening of the peak in the temperature dependence of the dielectric permittivity, and (iii) shift of the peak location under changing frequency of the applied electric field (in the GHz range).

Furthermore, since the accurate representation of the randomness of the local composition-field patterns in small unit cells requires statistical sampling, i.e., generation and processing of multiple instances (seeds) of the same system, with the following averaging of the obtained results, we found that, in the case of the dielectric response, these results depend strongly on the details of the chosen averaging approach. In particular, we observed that the local dipole population relaxation processes captured on the level of individual nodes or seed configurations exhibit mostly Debye-like behavior and, therefore, their averaging with nodal and seed schemes preserves this Debye-like character. Only in the global averaging scheme, where the “raw” seed-level dielectric response data is averaged first and then fitted to an empirical expression, such as, e.g., HN formula, is the non-Debye relaxation behavior detected (e.g., by the deviation of the asymmetry $σ HN$ and/or skewness $α HN$ parameters away from unity) under certain conditions. For the BSTO model system studied here, we have been able to identify two such regions of “non-Debyeness:” (i) close to the nonpolar to polar phase transition in the thermal evolution of the system, and (ii) in the general vicinity/not too far away of the system $T C$⁠, when the correlation length $L$ grows above a certain threshold value. The obtained insights could inform the design of novel relaxor compounds where the composition-field variations can be custom-tailored to produce the desired strength of the dielectric response, while minimizing losses.

The developed computational approach can be straightforwardly applied to elucidating the influence of other morphological features present in ferroelectric ceramic compounds and composites—polycrystallinity, structural defects, presence of polar domains, or other polar textures, etc.—on their dielectric response, and such work is currently underway. Furthermore, the random correlated field approximation—through the manipulation of the compositional and/or other dependencies of the materials coefficients of the system thermodynamic potential—would be useful for the continuum-level modeling of a variety of different systems, such as, e.g., micromagnetics (in the case of frustrated or non-stoichiometric magnets) or disordered thermoelectric materials.

In closing, we would like to comment on the possibilities for engineering composition variations, such as the ones considered in this study, into experimentally grown samples of BSTO or similar perovskite solid solutions. One attractive option that, in our opinion, could be investigated for this task is spontaneous phase separation by a spinodal nanodecomposition. This process and its growth conditions are a subject of active current research for semiconducting materials⁵¹ and for some perovskite oxide systems.^52–54 

The supplementary materials accompanying this manuscript include further details on the utilized computational approach, parameters for the employed total energy functional, as well as information on simulation results, such as the dielectric response in the transverse excitation mode, additional fittings of the dielectric response to empirical formulas, additional comparison of seed- and global-level averaged dielectric-response curves, and the dielectric response curves for the system with fixed Sr-composition.

A.G., S.P.A., and S.N. gratefully acknowledge the Air Force Research Laboratory, Materials and Manufacturing Directorate (AFRL/RXMS) for support via Contract No. FA8650–21–C5711. N.D.O. acknowledges support from the Innovations in Measurement Science Program at NIST. Some of this research was performed, while A.M.H. held an NRC Research Associateship award at NIST. The authors are grateful to Professor Menka Jain (UConn) and Dr. David Furrer (Pratt & Whitney) for useful suggestions, as well as to the developer team at Idaho National Laboratory for implementing the random correlated field in MOOSE. Distribution A. Approved for public release: distribution unlimited (No. AFRL-2023-1203). Date Approved 03-10-2023.

The authors have no conflicts to disclose.

Ashok Gurung: Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). John Mangeri: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Aaron M. Hagerstrom: Conceptualization (equal). Nathan D. Orloff: Conceptualization (equal). S. Pamir Alpay: Conceptualization (equal); Funding acquisition (supporting); Investigation (equal); Project administration (supporting); Resources (equal). Serge Nakhmanson: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (lead); Methodology (equal); Project administration (lead); Resources (lead); Software (equal); Supervision (lead); Validation (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.