{'fraction_non_alphanumeric': 0.03936344882751652, 'fraction_numerical': 0.01465078095398476, 'mean_word_length': 5.045813184016289, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling of dielectric hysteresis loops in ferroelectric semiconductors with charged defects

We have proposed the phenomenological description of dielectric hysteresis loops in ferroelectric semiconductors with charged defects and prevailing extrinsic conductivity. Exactly we have modified Landau-Ginsburg approach and shown that the macroscopic state of the aforementioned inhomogeneous system can be described by three coupled equations for three order parameters. Both the experimentally observed coercive field values well below the thermodynamic one and the various hysteresis loop deformations (constricted and double loops) have been obtained in the framework of our model. The obtained results quantitatively explain the ferroelectric switching in such ferroelectric materials as thick PZT films.


Introduction
The spontaneous electric displacement switching under the external field is one of the essential features of the ferroelectric materials [1], [2]. The main comprehensively studied characteristics of ferroelectric hysteresis loop are the spontaneous displacement and coercive field values [2], [3]. However, conventional theoretical approaches give significantly incomplete picture of the displacement switching, namely: • Nucleation theory [4] proved that the ferroelectrics inhomogeneities promote the domains nucleation, appearance of several nuclei in the ferroelectric film do not need to overcome the energetic barrier due to the long range interaction between them [5]. Calculated value of the coercive field is inversely proportional to the film thickness [6] and thus nucleation theory can explain the drastic increase of coercive field with film thickness decrease experimentally observed in some thin films [7], [8]. However, nucleation theory predicts that coercive field tends to zero at film thickness increasing, and this theory seems unsuitable for thick films and bulk materials. • The domain kinetics theory evolved in [9], [10] uses Kolmogorov-Avrami statistical model (see chapter 4 in [11]). It allows one to describe experimental data with high accuracy due to the great amount of fitting parameters without clear interpretation. The obtained results represent themselves only computer simulations. • Landau-Ginsburg theory (see [3], chapter 7 in [11]) evolved for the mono-domain perfect ferroelectrics describes switching without domains (so called homogeneous switching) and neither domain pinning, nor domain nucleation and domain movement. Calculated values of thermodynamic coercive field [2], [3] are from several times to several orders greater than the experimentally observed values for the real bulk ferroelectric materials (see e.g. [12], [13]). Renormalized Landau-Ginsburg free energy with coefficients depending on the film thickness and boundary conditions characteristics [14], [15] does not describe adequately ferroelectric hysteresis even in epitaxial film, because it predicts very sharp coercive field relaxation up to the aforementioned thermodynamic coercive field with film thickness increase [16]. • As for the imperfect doped ferroelectrics, here the hysteresis loops look much thinner [17]- [20] than Landau-Ginsburg ones, and its shape undergoes various deformations ("slim" loops [21], chapter 5 in [12]; minor loops [22], [23]; constricted loops [24], [25]). Therefore, the modification of Landau-Ginsburg free energy (and thus Landau-Khalatnikov equation [3] derived from it by variational method) suitable for the description of inhomogeneous ferroelectrics switching seems necessary. Incorporation of the modified free energy with inhomogeneities contribution could be a solution to this problem.
To our mind, the inhomogeneous electric fields caused by charged defects should be taken into account. However in contrast to the random field theory, developed for the relaxor insulator ferroelectrics [26], semiconductor properties of these materials [27] (at least extrinsic conductivity created by charged impurities [28], [29]) should be taken into consideration. Right in the way we modified Landau-Ginsburg approach for the ferroelectricssemiconductors with charged defects and found both the essential coercive field decrease and the hysteresis loops deformation experimentally observed in inhomogeneous systems.

The problem
The majority of imperfect ferroelectrics with non-isovalent impurities or some unavoidable imperfections would be considered rather as extrinsic semiconductors [27], [28] than perfect insulators.
We assume that almost immovable non-stechiometric defects or non-isovalent impurity centers are embedded into hypothetical perfect uniaxial ferroelectric (z direction coincides with the polar axis). We suppose that impurity centers or defects are ionized (e.g. after UV, photo-or thermal excitation) and even in the absence of proper conductivity, they provide a prevailing extrinsic conductivity in the bulk sample.
For the sake of simplicity, we regard that the sample as a whole is the electro-neutral n-type extrinsic semiconductor with positively charged defects with density ρ s (r). The microscopic spatial distribution of these defects charge density ) (r . The modulation δρ s fluctuates due to the great variety of misfit effects (lattice distortions, local shift, possible clusterization at high defect concentration). For the sake of simplicity, let us assume that charged defects spatial distribution is quasi-homogeneous, i.e. the average distance between defects 3 d n a d ≅ (a is lattice constant, n d is defect concentration). The size distribution function of charged defects is well localized near its sharp maximum at the average size 2r 0 . The movable screening clouds δn(r,t) surround each charged center (see Scheme 1a). The characteristic size of these screening clouds is of the same order as the Debye screening radius R D . When one applies the external field E 0 , screening clouds of free charges are deformed, and nano-system "defect center + screening cloud" becomes polarized (Scheme 1b). Polarized regions "δρ s (r)+δn(r,t)" cause the additional inner electric field fluctuations δE(r,t). According to the equations of state, the fluctuations of the inner electric field δE cause displacement fluctuations δD.
The charged defects with the charge density ρ s (dark circles with radius r 0 ) are screened by the free charges with density δn (gray circles or ellipses with screening radius R D ). The parts "a" and "b" show the system with the zero and nonzero external field E 0 respectively.

General equations
Hereinafter we assume that the period of external field changing is much greater than the free carriers relaxation time, and therefore the quasi-static approximation is valid. Thus, Maxwell's equations for the quasi-static electric field Е and displacement D have the form: They have to be supplemented by the material expressions for current and charge density: Here ρ m , µ m and κ m are the m-type movable charge volume density (m=n,p), mobility and diffusion coefficient respectively, j c is the macroscopic free-carriers current, ρ s (r) is the fluctuating charge density of static defects. Keeping in mind that the sample is the extrinsic semiconductor with prevailing n-type conductivity, hereinafter we neglect the proper conductivity, put ∑ ρ ≈ m m n and omit the subscript "m". So equations (1), (2) can be rewritten as: Here n<0, µ<0 and κ>0. In accordance with Einstein relation [28]. Hereafter we suppose that homogeneous external field E 0 (t) is applied along polar z -axis. The infinite in the transverse directions rather thick sample with equivalent boundaries occupies the region l l < < − z . We suppose that the potentials of electrodes are given, so that the inner field satisfies the condition: Hereinafter we introduce the averaging over sample volume V : Hereinafter the dash designates the averaging of functions . The spatial distribution of the deviations δf(r,t) consist of the part caused by spontaneous displacement screening [27] and localized in the ultrathin screening surface regions [30] and the one caused by microscopic modulation δρ s with quasi-homogeneous spatial distribution. For the µm-thick sample the contribution from the ultrathin screening region to the average functions is negligibly small, and the averaging over sample volume is equivalent to the statistical averaging with homogeneous distribution function. Using this distribution function properties, one can regard that: and the correlation between the different f δ -functions is equal to zero if the total power of the functions is an odd number. It follows from (3)-(5) that: Here e z s the unit vector directed along z-axis, E is the applied uniform field ) ( 0 t E and . Notice that the average values D E , are determined experimentally [1]- [3] most of the times. Having substituted (7) into (3) and averaged, one can obtain the expressions for the average quantities, namely: The absence of the space charge average density ρ follows from the sample electroneutrality. Using and obtain the following system of coupled equations: Here Г>0 is the kinetic coefficient, α<0, β>0, γ>0 are parameters of the hypothetical pure (free of defects) sample. Hereafter we denote D D z δ ≡ δ . We would like to underline, that the sum of equations (11)-(12) coincides with Landau-Khalatnikov equation [3] for D only under the condition δρ s =0. This condition corresponds to the absence of inhomogeneities or their seeding, and thus only the homogeneous polarization switching can take place when the external field exceeds thermodynamic coercive field. In such case the inner field in the bulk of the sample is the sum of the external field and predetermined depolarization field originated from displacement screening [27]. In contrast to this at δρ s ≠0 the inner field δE contains the random component dependent over δρ s (r) and δD in accordance with (9)-(10).
The system of equations (9)-(12) is complete, because the quantities E δ δ , n can be expressed via the fluctuations of displacement δD and δρ s (r) allowing for (9), (10). The spatial distribution and temporal evolution of the displacement ) , ( t r D in the bulk sample is determined by the non-linear system (9)-(12) supplemented by the initial distributions of all fluctuating variables (e.g. distribution of charged defects). However, in our opinion, only the adopted for the stochastic differential equations implicit numerical schemes can be used in order to obtain the numerical solutions of (9)- (12). Therefore, hereinafter we consider only the average characteristics of the inhomogeneous ferroelectric semiconductors. The study of the mechanisms of domain wall pinning by the given distribution of charged defects, domain nucleation during spontaneous displacement reversal is beyond the scope of this paper. Similar problems for the ferroelectrics -ideal insulators were considered in details earlier (see e.g. [32], [33]).

Coupled equations.
In order to simplify the nonlinear system (9)-(12) the following assumption has been used. The charged inhomogeneities are surrounded by screening clouds (see Fig. 1), therefore n s δ δρ~ and 2 s Small positive function η is determined by the ratio 0 r R D (see Fig. 1 and Appendix A).
When external field amplitude increases, η value slightly decreases due to the polarization of the system "charged fluctuation + screening cloud". We suppose that η can be approximated by effective constant value at not very high external field. After elementary transformations of (10), the equation for electric field fluctuation δE acquires the form: The equation (9) gives Using (15) and (14) the electric field fluctuations δE caused by charged defects can be expressed via δD and δρ s : The equations (11), (12), and (16) is the self-consistent system of the nonlinear integraldifferential equations for δD, δE and D . Its non-homogeneity is proportional to charge fluctuations δρ s and external field E 0 .
The approximate system of first-order differential equations for average displacement D , its mean-square fluctuation 2 D δ and correlation s D δρ δ can be derived after some elementary transformations of (11), (12) (see [34], Appendix B and (13), (16)). Thus, we obtain three coupled equations: Here the renormalized coefficients have been introduced. The renormalization of α takes into account the finite value of correlation length and classical renormalization of gradient term γ by Debye screening radius R D in the bulk of a sample [30]. The renormalization of Γ takes into account the finite value of maxwellian relaxation time in the right hand-side of (18b). It is originated from diffusion field s n grad µρ κ ) ( (see (10), (14) and definition µ π κ = n R D 4 ). This term has been neglected in the classical equations of Landau-Khalatnikov type as well as in our previous paper [35].
The system (17)  We would like to note that derived system of coupled equations (17) might possess chaotic regions, strange attractors as well as strongly non-ergodic behaviour and continuous relaxation time spectrum [11]. Any new system of such type demands a separate detailed mathematical study that was not the aim of this paper.

Ferroelectric hysteresis.
In this section, we consider the equilibrium solution of (17), which corresponds to the quasi-static external field changing. Let us demonstrate, how the ferroelectric hysteresis loop (i.e. the dependence of displacement D over the external field E 0 ) changes its shape under the presence of charged defects. First, let us rewrite equations (17) in dimensionless variables: Here, . Below we numerically analyze the solution of the system (18) for the following parameters . The equilibrium solution of (18) corresponds to the quasi-static external field changing. For the harmonic modulated applied field . At such low frequencies, the initial conditions determine only the initial curve. They do not play any important role in the loop shape (hereinafter we put them zero: and different small g~1 values are shown in Fig. 2. The increase of g value firstly leads to the coercive field decrease then to the constriction appearance, subsequent transformation to the double loop and finally to the loop disappearance. The broadening, slight tilting and smearing of the loop shape under applied field frequency increasing can be explained by dielectric losses increasing. It should be noted, that for R>1 and ξ→1, g=0, w→0 hysteresis loops are absent [34].  increasing. This effect is somewhat similar to the known "square to slim transition" of the hysteresis loops in some relaxor ferroelectrics [21]. Basic plot demonstrates the drastic decrease of coercive field for R=0. 1 . This means that coupled system (17) reveals inhomogeneous displacement switching, e.g. sample non-polarized as a whole splits into the oppositely polarized regions when external field reaches coercive value. In contrast to this, the Landau-Khalatnikov equation describes homogeneous switching with 0 ≡ δD independently on the external field. Therefore, system (17) is not equivalent to the equations of Landau-Khalatnikov type.  . In Appendix C we obtain the estimations for static remanent displacement ) 0 , , namely at 1 >> ξ : Note, that (19)  . We can conclude that the increasing of charged defect concentration (as well as its fluctuations) leads to the drastic decreasing of the coercive field, appearance of constricted and double hysteresis loops related to the inhomogeneous displacement switching. We demonstrate that this result much better agrees with experiments in thick films and bulk samples, than predictions of conventional Landau-Ginzburg and nucleation theories in the next section.

Comparison with experiment
Solid solutions Pb(Zr,Ti)O 3 (PZT) are widely used ferroelectrics, however much more work is needed to determine the physics of these materials [13]. A wide range of low level additives (0.0-5.0%) significantly influence on the dielectric properties of PZT at room temperatures. At liquid helium temperatures, all these "extrinsic" effects disappear and the experimental data agree with "single domain" intrinsic permittivity calculated from thermodynamic theory. Such donor additives as La, Nb, Nd or Ce at room temperature "soften" PZT dielectric properties, in particular hysteresis loops have rather high remanent displacement and low coercive field [13]. We suppose that aforementioned dopants as well as numerous unavoidable Pb vacancies (originated due to the high volatility of PbO, see chapter 10 in [12]) play a role of randomly distributed charged defects.
For  (17)- (19). However, up to 5µm films such interface phenomena as film-substrate mismatch and charged layers cause loops asymmetry [15] due to the induced electric field, and mismatch-induced homogeneous elastic stresses renormalize the free energy coefficients α , β [37].
In accordance with (18) ) as fitting parameters determined by unavoidable cations vacancies depending over molar fraction x and/or Nb 5+ concentration. Our fitting is presented in Figs.6-8. The discrepancy at high external field between calculated and measured the loops is regarded to the 4 D -approximation used in our model. In order to improve the fitting the higher terms ... , 6

D
should be taken into account. almost independent over external field amplitude E 0 (so called saturated "square" loops), whereas experiment data for PZT show its strong dependence over E 0 (see Fig.7b). In contrast to Landau-Khalatnikov model, our model shows quantitatively agreement with experiment (see Fig.7).  (plot a). Coercive field as a function of applied electric field amplitude (plot b). Squares are experimental data for 20µm-thick Nb-modified PZT film [19], solid line -our theory. Parameters are the same as in Fig.6.
The right-hand shift of the 1.9µm-thick PZT film is related to the laking layer near the boundary metal-semiconductor at the bottom electrode (see chapter 14 in [38] and Fig.8).
). It is follows from Tab. 1 and Figs.6-8 that our model can quantitatively describe typical ferroelectric hysteresis loops in thick "soft" PZT films, in particular it gives correct coercive field values and its dependence over external field amplitude.
The constricted and double loops were observed in PZT ceramics doped with Nd [18], La [24], [25]. Sometimes constricted loops disappear after baking at high pressures, annealing in special atmospheres or several hundreds of switching cycles [18], [24], [25]. In order to explain this effect, we assume that ceramics treatment or fatigue can cause the significant changing of charged defects spatial distribution (e.g. the characteristic distance d and concentration n d ). The latter leads to the changing of gR 2 and ξ values in accordance with (18). Thus either appearance of constricted and double hysteresis loops in aged materials [18], [25] or their disappearance depends on the values of parameters R, ξ and g after sample treatment.
It should be noted, that the origin of the constricted or double loops in aged ferroelectric BaTiO 3 and (Pb,Ca)TiO 3 ceramics caused by the mechanical clamping of spontaneous polarization switching [12] and so it lies outside our theoretical consideration. This may be related to the fact that possible evolution of the charge fluctuations caused by the relaxation/origin of internal stresses around defects was not taken into account in our model.
Let us briefly compare our theoretical results with classical ones concerning lightly donor doped BaTiO 3 . In many cases nonisovalent additives or unavoidable impurities containing even in very restricted quantity significantly attenuate piezoelectric properties of BaTiO 3 ceramics, increase electric conductivity, diffuse the peak of dielectric response and make hysteresis loop slim with low coercive field (see chapter 5 in [12]). That is why we consider the influence of charged defects on the ferroelectric properties of such "soft" materials, neglecting the contribution of inhomogeneous mechanical deformations arising due to intergranular stresses as well as local symmetry distortion appeared near the defects.
The (0.1-0.3)% La, Nb or Ce-doping changes wide band-gap intrinsic semiconductor BaTiO 3 to an extrinsic n-type semiconductor, achieved by replacing Ba 2+ or Ti 4+ with ions Le 3+ , Ce 3+ or Nb 5+ with higher valency [12], [39]. In the most of BaTiO 3 ceramics such unavoidable impurity as Fe 3+ is present. In such case, the sample has slightly brown hue. Such ceramics reveals "soft" hysteresis loops with small area and low coercive field, diffuse dielectric properties and disappearance of spontaneous displacement (see [12], Tab. 2), in contrast to BaTiO 3 single crystals which undergoes first order phase transition with abrupt displacement disappearance at T=110°C. We supposed, that coupled equations (18) evolved for ferroelectrics with second order phase transition, can be applied to the doped or "soft" BaTiO 3 ceramics. Obtained results are summarized in Fig. 9 and Tab.2.  It is follows from Tab. 2 and Fig.9 that our model can quantitatively describe typical ferroelectric hysteresis loops in "soft" BaTiO 3 ceramics, in particular it gives correct coercive field values.

Conclusion
We have proposed the phenomenological description of polarization switching peculiarities in ferroelectric semiconductors with charged defects and prevailing extrinsic conductivity. Exactly we have modified Landau-Ginsburg approach for the aforementioned inhomogeneous systems and firstly obtained the system of coupled equations (17). Solving the system (17) one can get the information about system ordering as a whole, without defining space distribution of the appeared displacement inhomogeneities, however the present model has the following advances over the conventional ones and our previously published works [15], [34]- [35]: • The coupled system (17) reveals inhomogeneous displacement switching in contrast to the Landau-Khalatnikov equation, which describes homogeneous one. • We have shown that the increasing of charged defect concentration (as well as its fluctuations) leads to the drastic decreasing of the coercive field, appearance of constricted and double hysteresis loops. • Obtained results quantitatively describe typical ferroelectric hysteresis loops in thick donor doped PZT films and BaTiO 3 ceramics, in particular our model gives correct coercive field values and its dependence over external field amplitude, in contrast to the conventional Landau-Ginsburg and nucleation theories.
. If charged defects are clusterized distance d can be several times greater than 3 1 d n . We consider the case when defects concentration n d does not exceed several percents, i.e. d>3r 0 (see Fig. 1). For the sake of simplicity, we approximate the defects charge densities by isotropic Gaussian distributions. The charge density of carriers screening clouds originated around charge defects can be estimated in Gaussian approximation using the Debye potential [40] and Green function method, namely at D R r0 we obtain: Here f n is the charge density of free electrons in the conduction band of extrinsic semiconductor [28], [38]. In our case high defects concentration % 1 Allowing for , intrinsic conductivity could not be neglected in comparison with extrinsic one and therefore the obtained expressions (A.6)-(A.7) become incorrect.

Appendix B
The equations for 2 D δ and s D δρ δ obtained directly from (12) have the form: Under multiplying (12)  (B.8) In accordance with (13) the last two terms in the right-hand side of (B.6) can be estimated as As it follows from (17a) , and the modified equation for D can be obtained directly from (C-1), namely Let us remind, that (18)).