A Simulation Model for Narrow Band Gap Ferroelectric Materials

Various ferroelectric simulation models have been developed in recent decades in order to study the mechanisms and predict the behaviors of ferroelectrics by simulating their hysteresis loops. Conventional ferroelectrics have wide optical band gaps (>2.7 eV), making them difficult to respond to visible light. Although their ferroelectricity can be affected by higher‐energy radiation like ultraviolet, little attention is paid to the development of their models incorporating light‐dependent factors. However, in recent years, narrow band gap (<2.5 eV) ferroelectrics have been discovered and are increasingly researched. These special ferroelectrics effectively absorb visible light and hence exhibit strongly light‐dependent ferroelectricity, triggering potentially a broad range of applications. Therefore, there is a need to develop a model for these ferroelectrics in order to predict their behavior under visible light. Such a model is also needed to improve the understanding of the interaction mechanisms between light and domains. In this paper, a ferroelectric simulation model based on the Jiles–Atherton theory considering light dependence is developed for the first time, and its accuracy is validated by experiments. The model shows an average error of 7.5% on polarization values compared to experimental results and thus can be employed to reliably predict photo‐induced ferroelectricity.


Introduction
Photoferroelectrics are ferroelectric materials exhibiting photovoltaic effect at the same time. [1] Conventional and widely known ferroelectrics like PZT (Pb(Ti, Zr)O 3 ) and BaTiO 3 are the most studied photoferroelectrics. [2] These materials have been extensively used as kinetic and thermal sensors, actuators, transducers, and energy harvesters thanks to their excellent DOI: 10.1002/adts.202000052 ferroelectricity thus strongly induced piezoelectricity and pyroelectricity. Nevertheless, there has not been any practical case for photovoltaic or optoelectric purposes due to the fact that conventional ferroelectrics yield wide optical band gaps, usually >2.7-3 eV (Si band gap 1.1 eV and narrowgap semiconductors 1.4-1.8 eV). [3,4] When energy of the incident radiation is high enough, i.e., >3 eV, these ferroelectrics can show light-dependent ferroelectric hysteresis loops (or the so called photo-induced ferroelectricity). [5] However, since the majority of the solar radiation is visible light, such materials can hardly trigger practical applications on this regime and therefore they have not become a subject of deep theoretical investigations in terms of simulation models for the photo-induced ferroelectricity.
The situation has changed since the discovery of narrow band gap photoferroelectric materials, for instance, (K, Ba)(Nb, Ni)O 3-(band gap 1.1 eV), [6] Bi(Fe 0.5 Cr 0.5 )O 3 (band gap 1.4 eV), [7] and (K, Na, Ba)(Nb, Ni)O 3-(KNBNNO, band gap 1.6 eV). [8] These narrow-gap ferroelectrics show a strong visible light absorption as well as lightdependent ferroelectric hysteresis loops. [9] The interaction between domain wall motion and incident light has greatly allured research interest for potential applications as the next generation opto-ferroelectric devices like light-effect transistors, lightwriting memories, multi-functional sensors, and multi-source energy harvesters. [1,10] However, there has not been any established model for the simulation of photo-induced ferroelectricity despite there is an urgent need to predict the ferroelectric behavior under visible light and to theoretically study the mechanisms of domain-light interaction in these narrow-gap photoferroelectrics. Therefore, this paper is dedicated to address this issue.
It is well known that ferroelectrics are generally sensitive to electric, mechanical, and thermal stimuli. The polarization state of ferroelectric materials is dependent on ambient physical excitation, such as electric field (E), pressure ( ), and temperature (T). [11] The principle of minimum energy justifies the domain presence and the highly nonlinear hysteretic behavior, linking the polarization field to all potential external excitations. Simulating variations of the polarization state is crucial for understanding all fundamental physical behaviors of the materials. It is also a mandatory route toward the simulation of complex systems based on piezoelectric/pyroelectric energy conversions.
www.advancedsciencenews.com www.advtheorysimul.com Photoferroelectric KNBNNO, as mentioned above, is a perovskite-structured, narrow band gap ferroelectric material. [8] Its polarization state, like all counterparts, depends on E, , and T. However, KNBNNO's particularly narrow band gap (1.6 eV, responding to all visible wavelengths) makes the polarization state highly sensitive to visible light. This largely differs from the situations of its wide-gap counterparts where only ultraviolet can cause a noticeable polarization variation. In this paper, the E-dependent ferroelectric version of the classic ferromagnetic Jiles-Atherton model [12] is developed to T-dependent and light-dependent models for the simulation of the ferroelectric behavior of KNBNNO ceramics, when subject to temperature and lighting condition changes. The T-dependent ferroelectric model for KNBNNO has been preliminarily discussed in our previous research. [13] However, due to unrecognizable change of spontaneous and remanent polarizations of KNBNNO at relatively low temperatures (up to 150°C) and with electric fields (up to 6 × 10 6 V m −1 ) well beyond the coercive field (1.5-2 × 10 6 V m −1 ), the preliminary model could not be validated at higher temperatures (i.e., across the phase transition range). In this paper, smaller electric fields (up to 2.5 × 10 6 V m −1 , but still larger than the coercive field) are applied and the temperature range is extended to beyond the orthorhombic-tetragonal ferroelectric-ferroelectric phase transition temperature (170°C). With these external excitations, the changes of polarization status become recognizable and thus the T-dependent model is completed. Meanwhile, a ferroelectric model of light-dependent polarization status for KNBNNO is proposed for the first time in this paper, paving the way toward models for the prediction of photo-ferroelectric behaviors. The light-dependent model shows an average error of 7.4% on polarization values compared to experimental results, especially in the saturation range (nearly zero error). This qualifies the model to be used for reliably predicting the photoinduced ferroelectricity in narrow band gap photoferroelectrics like KNBNNO.

Simulation Method
It should be noted that even if E, T, and polarization (P) are vectors, the experimental validations used in this paper can be assumed to have perfect collinearity. As a consequence, all these physical quantities can be expressed directly by their magnitudes. Hysteresis in ferroelectrics can be considered as a nonlinear delay separating external stimulus to the polarization response. The polarization variation of a ferroelectric material becomes insensitive to the excitation dynamics when subject to certain frequency which is the quasi-static threshold. The state below the quasistatic threshold is called the quasi-static state. Similar phenomena can be observed with ferromagnetic materials.
In 1980s, Jiles and Atherton established a new theory for the simulation of the quasi-static scalar hysteresis behavior of ferromagnetic materials. [12,14,15] This theory is particularly elegant as the proposed model can be simulated with only five parameters and it consists of solving a first-order differential equation for each timestep of the simulation. A large number of extensions and improvements have been proposed since the establishment of the Jiles-Atherton model, suggesting to take into account newly discovered behaviors in the simulation. These considerations include inversion of the model, [16] vectorization of the model, [17,18] and frequency dependence of magnetization variations beyond the quasi-static threshold. [19] Following the ferromagnetic simulation, the Jiles-Atherton model has also been developed to satisfy the simulation of ferroelectric bulk ceramics. [20] Those five major parameters have been kept the same while M (magnetization), M rev (reversible magnetization), and M irr (irreversible magnetization) have been replaced by P, P rev (reversible polarization), and P irr (irreversible polarization), respectively.
Although the Jiles-Atherton model as well as its ferroelectric modification is not the most accurate among all simulation schemes (e.g., it suffers from accommodation issue [21] ), it is a versatile and flexible model which can be well adapted on multiple couplings between polarization and external stimuli as observed in KNBNNO. [13] Therefore, this model is chosen for studying the functional simulation methods of KNBNNO and for the interpretation of experimental observations regarding to the temperature-and light-induced polarization variations. The improvements added to the original model (mentioned above) are also considered.

Ferroelectric Jiles-Atherton Quasi-Static Hysteresis Model
The ferroelectric version of the Jiles-Atherton model, as proposed in reference, [20] provides the evolution of P as a function of E. It gives relatively accurate results but constrains experimental conditions. It requires (1) low frequency (f < 1 Hz, quasi-static state) and (2) collinearity of P and E. In this model, P is composed from both its reversible and irreversible constitutes, defined by Equation (1) where c is the polarization reversibility coefficient and P an is the anhysteretic polarization. P an is a sigmoid function of E which can be expressed from a Langevin-type function shown as Equation (2). P an can also be derived from the Ising spin model, shown in Equation (3). In Equations (2) and (3), P s is the saturation polarization, a is the density of domain walls in the ferroelectric material, E e is the effective electric field which is equivalent to the Weiss mean field acting on atomic moments within a ferroelectric domain, and is the inter-domain coupling. E e is described by Equation (4).
Equation (2) is valid only for isotropic materials. In this work, the initial states of KNBNNO samples were considered macroscopically isotropic because they were ceramics with randomly distributed/oriented grains and domains after fabrication and before any electrical measurement. In principle the anhysteretic polarization curve can be observed experimentally, and in practice this signature should be obtained when a ferroelectric ma- terial is fully depolarized under the influence of a constant bias electric field. However, such a sophisticated experimental process means accurate measurements of the current density variations and thus exhibits significant uncertainty due to various sources of drifts that create divergences during the integration steps, leading to extremely complex verification. Therefore, two assumptions were made in this work: (1) As no anisotropy was intentionally introduced into the KNBNNO samples before measurements (apart from their inherent anisotropic but randomly oriented unit cells as explained above), the polarization behavior was treated as isotropic and the anhysteretic properties could be resumed to an isotropic contribution (Equations (2) and (3)); (2) The anhysteretic polarization curve could be calculated from a full ferroelectric hysteresis loop as illustrated below.
For KNBNNO, a better accuracy for the simulation of the anhysteretic polarization can be achieved by introducing a second term into Equation (4). This second term behaves as an equivalent parallel capacitance, thus, its influence is significant during the saturation stage where it creates a small positive drift. Equation (5) expresses the situation considering the second term.
In Equation (5), 0 and r are vacuum and relative permittivity, respectively. The practical anhysteretic curves are complicated to obtain in experiments due to the need of the superimposition of exotic waveform excitation contributions. The experimental acquisition processes are considerably time-consuming and the results are hardly trustworthy. Reasonable approximation can be made by averaging an experimental ferroelectric hysteresis loop. For a given E (E i ), the correspondingly estimated P an (P an-i est ) can be calculated from Equation (6) where P i inc and P i dec are the increasing and decreasing parts of the polarizations in a hysteresis loop, respectively. Figure 1 shows the comparison of an experimental hysteresis loop (measured at 25°C with 0.1 Hz excitation frequency in dark), the reconstituted anhysteretic curve (according to Equation (6)) and the simulated anhysteretic curve (according to Equation (5)). The parameters used here are listed in the figure.
Equation (7) gives the relation between P irr and P an , where k is the pinning energy factor and is a directional By combining Equations (1), (4), (5), and (7), the Jiles-Atherton model's main differential equation can be obtained in Equation (8). Equation (8) and its constituents are calculated and solved for each timestep of the simulation. dP/dE represents the differential permittivity.
In this study, the extension of the Jiles-Atherton model required an inverse hysteresis quasi-static contribution to consider the frequency effect, [13], i.e., an evolution of E as a function of P (E(P)). The inverse version of the ferroelectric modification of the Jiles-Atherton model in this study was adapted from the www.advancedsciencenews.com www.advtheorysimul.com ferromagnetic counterpart, [16] shown by Equations (9)- (11). P e is the effective polarization. Figure 2 illustrates the accuracy of the model by comparing the simulated and measured ferroelectric hysteresis loops. Parameters which are not given in Figure 1 are provided in Figure 2.
In the following sections of the simulation description, the influence of frequency, temperature, and lighting conditions will be considered. The variations of the dielectric properties directly on the evolution of ferroelectric hysteresis loops will be studied. The studies are compared to the reference state shown in Figures 1  and 2

Ferroelectric Dynamic Behavior
The hysteretic evolution of polarization state as a function of external excitation in a ferroelectric ceramic material is always frequency dependent. In the range of an ultralow frequency up to the above-mentioned quasi-static threshold, the shapes of the hysteresis loops obtained in this frequency range barely differ from each other. With an increase of frequency, the area enclosed by the hysteresis loop increases first toward the maximum and then continuously decreases, until the frequency is too high to be able to induce any domain wall motion thus no change can be observed anymore. The frequency dependence of hysteresis in a ferroelectric ceramic, in most cases, is simulated as a firstorder viscous damping which is the consequence of a constant ( ) multiplying the time derivation of polarization state. The validity of a model in comparison to experimental results across a large frequency range can be significantly improved by replacing this first-order time derivation with a fractional one. [22][23][24][25] The physical explanation of using fractional derivatives in the field of mechanics is known as the anomalous diffusions. In ferroelectricity, anomalous diffusions can also be envisaged, that is a mix of thermal and strain diffusions incorporated into the ferroelectric nonlinear domain wall motion. In this study, by solving Equation (12), the frequency dependence in the ferroelectric hysteresis model was derived. In Equation (12), n is fractional order, E dynamic is the dynamic electric field, and J.A −1 static represents the quasi-static contribution of the inverse of Jiles-Atherton model (ferroelectric modification). For a given polarization state, the difference between E dynamic and J.A −1 static is equal to the product of multiplying the time fractional derivative of polarization. Figure 3 compares ferroelectric hysteresis loop obtained from the simulation to those of measurements with frequencies of 1-100 Hz. The loops in Figure 3 do not show any significant change with frequency, implying that all the experimental measurements were performed below the quasi-static threshold. It can then be concluded that the quasi-static model is suitable for the simulation of KNBNNO ceramics with frequencies up to 100 Hz. Higher frequencies were not investigated due to the frequency limitation of the voltage amplifier (609B, Trek, USA) used to apply electric fields.

Influence of Temperature
The effect of temperature on ferroelectric properties of classic ferroelectric materials are, although complex, well established. Starting from the 0 K state, the ferroelectric properties will first remain unchanged, followed by a slight decrease of both saturation polarization and coercive field-Type I. Alternatively, it may show an increase of hysteresis indicators until a maximum and then a decrease until the hysteresis completely vanishes when approaching to Curie temperature-Type II. These changes can be expressed by the area enclosed in a hysteresis loop A as a function of temperature-A(T).
Due to the increase of leakage current with temperature, [26] the hysteresis loops of KNBNNO could be reliably measured up to only 200°C (i.e., the polarization at the maximum electric field was not smaller than the remanent polarization) (see Figures S1-S3, Supporting Information). In the measurement, an increase of the loop area with temperature was observed, therefore, it could be concluded that KNBNNO belongs to Type II ferroelectrics theoretically. In simulation, this means that the parameters set for the increasing part of the hysteresis needed to be changed accordingly from the decreasing one. The temperature effect on ferromagnetic hysteresis based on the Jiles-Atherton model has been proposed. [27] A modified approach for ferroelectric hysteresis has also been reported. [13] The same method was also used here for verifying the KNBNNO. Five parameters are given here to adjust the temperature dependence, and they are as follows: Here "0" corresponds to the 0 K state, for instance, P s (0) is the saturation polarization at 0 K. T C is Curie temperature. As experimental conditions at 0 K could not be achieved, in the simulation P s (0), k(0), a(0), (0), and c(0) are assimilated to their room temperature values. is a material dependent critical exponent. Figure 4 shows the comparison between the simulated and measured ferroelectric hysteresis loops at different temperatures. Relevant parameters are also given. The positive value confirms that the experimental measurements have all been done in the increasing part of the A-A(T) curve. The T C of 703.15 K is consistent with the reported experimental value (673.15 K), [8] proving the high accuracy and wide applicability of the model developed in this study. A complete set of data (hysteresis loops and polarization values) can be found in Figures S1-S6 (Supporting Information).

Influence of Incident Light
As discussed above, conventional wide-gap ferroelectrics are generally sensitive to electric, thermal, and mechanical stimuli. However, the effect of incident light on ferroelectric properties has rarely been simulated due to their insensitivity to visible light. KNBNNO, unlike those conventional wide-gap ferroelectrics, has a strong dependence of polarization on light due to its visiblerange band gap. Here in this study, KNBNNO ceramics have been exposed to five different illumination conditions including that in dark. Figure 5 illustrates the characteristics of the five light states. The conditions were set with the controlling software of   the light source using the "natural white light" mode and then by changing the intensity (brightness). Subsequently, the parameters describing the lighting conditions were measured (see the Experimental Section) and are shown in Figure 5.
As indicated by experimental data, [9] the shapes of the ferroelectric hysteresis loops show a strong dependence on the lighting condition. All typical hysteresis indicators including coercive field, saturation polarization, and remanent polarization increase with the average power density of the incident light. The change can be physically interpreted as the consequence of simultaneous generation of free charges and surface dielectric macro-dipoles. Incident light will affect the interfacial polarization and thus the related domain wall motion. [9,10] For KNBNNO, the shape variation of the hysteresis loops with lighting condition shows a similar trend to that of temperature. Therefore, a similar mechanism will be employed here to introduce the factor of light into the model, with five light-dependent parameters added as follows: Here in the parameters, LPD stands for light power density and LPD max corresponds to the maximum of the power density. In this study, it was a fictitious value obtained by comparing the simulation and measurement results and thus was set when the optimum results were observed. The subscription "D" represents the condition in dark. is a material dependent critical exponent. Figure 6 compares the simulated and measured ferroelectric hysteresis loops in different light states. Relevant parameters are also given in the figure. The LPD max value shown in the figure implies that the model works with incident light that has an average power density of up to 70 mW cm −2 . This is close to the widely used standard of AM 1.5 G (100 mW cm −2 ) for characterization of solar cells. Therefore, this model can also be used to roughly estimate the ferroelectric behavior under sunlight.
It should be noted that, before conducting each measurement the KNBNNO samples were exposed to the incident light for 30 min until the temperature change induced by light was eliminated in order to rule out any effect caused by temperature variation. The sample temperature was monitored by a thermometer. At 25°C, the incident light could lead to a temperature increase of only 1-2°C, which is negligible for the temperature effect on the sample. Therefore, the light-dependent change of hysteresis loops was believed to be induced by incident light rather than temperature change. This issue has also been discussed in literature. [9,10] A further improvement of the model can be done by considering the leakage current which was caused by the decrease of resistivity under light for KNBNNO, as shown in Figure 7. The leakage current was taken into account by introducing an electric field-dependent contribution as follows: R is the measured resistivity and is a constant dependent on sample dimensions. This improvement of the simulation is particularly suitable for the saturation area of the ferroelectric Adv. Theory Simul. 2020, 3, 2000052 www.advancedsciencenews.com www.advtheorysimul.com  (5) and (6) All parameters used in conventional ferroelectric J-A models P s = 0.1 C m −2 a = 1.8 × 10 6 r = 2 × 10 −9 2 P at room temperature in dark  hysteresis loop, as shown in Figure 8. Relevant parameters are also given in Figure 8. Table 1 summarizes the simulations carried out in this paper with the developed model. The target simulation factor, equations used and parameters kept constant and adjusted for each step are listed.

Quantitative Indicators and Errors
In this work, three quantitative indicators were used to assess the accuracy of the developed simulation method. They are: 1) Accuracy of the hysteresis area, 2) Accuracy of the coercive field, 3) Accuracy of the remanent polarization, where A, E c , and P r are area of hysteresis loop, coercive field, and remanent polarization, respectively and the subscripts "sim" and "mes" represent simulation and measurement results, respectively. Table 2 lists the values of the above quantitative indicators for the results shown in Figures 4 and 6.
According to Table 2, corresponding errors can be calculated by 100-〈A〉, 100-〈E c 〉, or 100-〈P r 〉. The most important parameter is P r as its variation changes the vital functions like piezoelectric, www.advancedsciencenews.com www.advtheorysimul.com pyroelectric, and photovoltaic properties. The average error of P r in this study was calculated to be 7.5%. In addition, the values of the quantitative indicators did not show any trend with changing temperature or lighting condition, affirming that the accuracy of the model is independent of external stimuli thus is reliable.

Conclusions
A ferroelectric simulation model adapted from the Jiles-Atherton theory for ferromagnetic materials has been developed. Dynamic, thermal, and optical factors of external electric field frequency, temperature, and lighting conditions, respectively, have been considered in the model. It is the first time that the incident lighting condition has been taken into account for a ferroelectric model. A recently discovered narrow band gap photoferroelectric ceramic, KNBNNO, has been used to validate the model. The simulation results have matched the experimental ones well, with average errors of 7.7%, 9.5%, and 7.5% for area of the hysteresis loop, coercive field, and remanent polarization, respectively (compared to experimental results). This model can be widely used to predict the behaviors of narrow band gap ferroelectrics, especially when exposed to visible light. Examples of practical applications of the model include rapid and non-experimental prediction, assessment and simulation of opto-thermal-electric coupling phenomena of narrow band gap ferroelectric materials. The change of polarizations significantly affects the pyroelectric and photovoltaic behaviors, but carrying out comprehensive experiments to study such influence is time-consuming while expected results cannot be guaranteed. Therefore, using the model to simulate changes of polarizations under different thermal and optical conditions thus to correlate them to the coupled pyroelectricphotovoltaic effects becomes essential in experimental studies of narrow band gap ferroelectrics. Along with potential practical applications, this model also stimulates the theoretical understanding of the mechanisms of domain-light interaction in photoferroelectrics under different thermal and optical conditions by quantitatively correlating the change of polarizations with incident light intensity.
Experimental Measurement: The fabricated KNBNNO ceramic samples were prepared with two electrode configurations. Configuration 1: Approximately 10 µm silver electrode was screen-printed on both sides, followed by firing at 600°C for 20 min. This configuration made at a high temperature allowed the samples to be measured with elevated temperatures. Configuration 2: 200 nm thick indium tin oxide (ITO) electrode was coated on both sides with sputtering. Visible light could pass through the transparent ITO, and thus this configuration was used for measurements under light. The ferroelectric hysteresis loops and resistivity data were obtained with a ferroelectric test system (Precision LCII, Radiant Technologies, Inc., USA). The measurement temperature was controlled by a sample stage (LTS 350, Linkam Scientific Instruments, UK). A smart light-emitting diode bulb (11586, 806 lm, AwoX, France) of which the illumination spectrum could be changed was used as the light source. The conditions were set with the controlling software of the bulb using its "natural white light" mode and then by changing the intensity (brightness). The powers and power densities at different wavelengths of the incident light were characterized with a PM100D optical power and energy meter integrated with an S120C silicon photodiode detector (Thorlabs, Germany). The color temperature and illuminance were measured by an optical sensor (RGBW200, Vishay, USA) integrated with a color spectrometer (ELV, Germany). Before conducting any measurement, the KNBNNO samples were exposed to the incident light for 30 min in order to wait for the temperature to stabilize and to rule out any effect cause by temperature variation. The sample temperature was monitored by a thermometer. At 25°C, the incident light could lead to a temperature increase of only 1-2°C, which was negligible for the temperature effect on the sample. Therefore, the light-dependent change of hysteresis loops was believed to be induced by incident light rather than temperature change.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.