CC BY
26Physics of Complex Systems, 2022, vol. 3, no. 3 _www.physcomsys.ru
R Check for updates
UDC 538.9 EDN QAMJZS
https://www.doi.org/10.33910/2687-153X-2022-3-3-122-136
Numerical simulations of nonlinear and chaotic order parameter responses in bulk antiferroelectrics using ammonium dihydrogen
phosphate parameter
S.-Ch. LimH1
1 Universiti Sains Malaysia, 11800 USM Penang, Malaysia
Authors
Lim Siew-Choo, ORCID: 0000-0001-8397-0886, e-mail: sclim@usm.my
For citation: Lim, S.-Ch. (2022) Numerical simulations of nonlinear and chaotic order parameter responses in bulk antiferroelectrics using ammonium dihydrogen phosphate parameter. Physics of Complex Systems, 3 (3), 122-136. https://www.doi.org/10.33910/2687-153X-2022-3-3-122-136 EDN QAMJZS " Received 22 April 2022; reviewed 9 June 2022; accepted 9 June 2022. Funding: The study did not receive any external funding.
Copyright: © S.-Ch. Lim (2022). Published by Herzen State Pedagogical University of Russia. Open access under CC BY-NC License 4.0.
Abstract. In this paper, the nonlinear and chaotic responses of bulk antiferroelectrics are elaborated phenomenologically and numerically. The first ordered phase of bulk antiferroelectrics is formulated by applying calculus of variations to Landau free energy density expansions of bulk antiferroelectrics. With applied time-dependent electric field, the antiferroelectrics dynamic responses are obtained by Landau-Khalatnikov equation of motion. The resulting dynamical equations are two nonlinearly-coupled second order differential equations corresponding to two inter-penetrating sub-lattices of antiferroelectrics, and these are solved numerically using forth-order Runge-Kutta methods and ammonium dihydrogen phosphate parameters in its first ordered phase. These calculated results are presented graphically for various frequencies and amplitudes in the applied electric fields.
Keywords: antiferroelectrics, ammonium dihydrogen phosphate, chaos, Landau free energy density, nonlinear, Poincare Sections
Introduction
Some switching processes and nonlinear chaotic dynamics in second order ferro- and anti-ferromagnetic systems have been studied theoretically using Landau-Lifshitz equations of motion (Chan 2010; Toh 2009). For ferroelectric system, i. e., triglycine sulfate, the similar studies have been carried out theoretically and numerically using Landau theory, and experimentally by measuring responses of ferroelectric capacitor in series resonance circuit (Diestelhorst 2003). The switching processes and some conventional nonlinear responses in ferro- and anti-ferroelectrics have been discussed and studied (Lines, Glass 1977; Tan 2001). However, the nonlinear chaotic behaviors of antiferroelectrics are still left behind. It is of interest to study nonlinear chaotic dynamics in antiferroelectric systems. The context of this paper is focused on theoretical and numerical studies of chaotic dynamics in bulk antiferroelectric system in its first ordered phase. In the numerical simulations, the fourth-order Runge-Kutta method is adopted (Press et al. 1996, 704; Strogatz 2018, 33), and is based on the parameters of ammonium dihydrogen phosphate (ADP) at 80 K (Ledzion et al. 2004). These numerical results are presented graphically.
Theoretical physics. Theory of condensed matter
Formalism
The formalism in this section is adopted and generalized from the Landau theory of first order ferroelectrics (FE) and second order antiferroelectrics (AFE) explained in Lines and Glass (1977). In general, the bulk AFE thermodynamic potential, G, is
G = U-TS-XjXJ-EiDi (1)
where subscript i represents all the existing field components in the AFE system. U is the internal energy of the system, Tis the temperature, S is the entropy, X. is the thermal stress, x. is the strain, E. is the electric field induced by the system's spontaneous polarizations and the applied electric field, and D. is the electric displacement.
For the study of AFE, it is convenient to describe the system as two interpenetrating sublattices, namely sublattice-A and sublattice-B. Let DA and DB represent electric displacements, and EA and EB represent electric fields for sublattice-A and -B, respectively. With the assumptions that the system has zero stress and strain, the volume and temperature of the system remains constant in the applied electric field, the thermodynamic potential is isobaric and isothermal, Equation (1) can be rewritten as a polynomial of order parameters of sublattices, i.e., displacements of the two sublattices, DA and DB:
G = a1(D2A+D2B)-2an(D4A+D4B) + am(D6A+D6B)+ VDADB -EADA -EBDB (2)
where a1, a11, and a111 are constants, referred to as Landau coefficients, and n is the interaction constant of the two sublattices. For AFE, n has positive value, which facilitates an antipolar situation. The negative signs of the fourth order terms represent a first order transition of the AFE system, and the constant 2 in front of a11 is for the convenience of the following derivations. In order to enable the fitting of experimental measurements, additional pairs of variables, i. e., staggered displacement R = DA - DB, staggered field ES = EA - EB, normal displacement Q = DA + DB, and the Maxwell field EM = EA + EB are defined (Lines, Glass 1977, 88). In terms of Q, R, EM and ES, Equation (2) becomes
v2 4.
R1-?1L(q*+6Q2R2 + Ra) +
n 11 (3)
I + ^(Q6 +15 Q*R2 +15 Q2R*+R6) - V-EMQ-V-ESR
For antiferroelectrics, the staggered field is necessarily zero, i. e., ES = 0. a1 is temperature dependent, and is defined as ax = -^- + /?(77 — rc). Then, Equation (3) becomes
G = ±\r, + p(T-Tc)]Q2+\p(T-Tc)R2-^(Q* + 6Q2R2+R*) +
+ ^{Q6+\5Q*R2+\5Q2R*+R6)-±EmQ
(4)
where ji is a constant, and TC is the Curie temperature of the AFE system.
In order to ease the numerical simulations, reduced variables, or dimensionless quantities, we introducedgA, t, e, q, r, and ty, corresponding to thermodynamic potential of AFE system, G, temperature, T, applied Maxwell field, EM, normal displacement, Q, and staggered displacement, R, and interaction constant of sublattices, n, as defined in Equations (5):
3
a,7
^ 16o[ll fc \
G = —y±gA (5a)
"■in
4 a2
P{T-Tc) = ^t (5b)
«111
2 «in
(5c)
Q = J^ (5d)
V «111
R= ^-r (5e)
V «ill
n = (5f)
«111
The substitution of Equations (5) in Equation (4) results in dimensionless thermodynamic potential:
gA = (w +1) q2 + tr2 - (q4 + V + / ) + + 15q4r2 + 15q2r4 + r6) - eq. (6)
Equation (6) is used in the formalism of nonlinear dynamical equations in the following sections. ty and t values are fitted by using material parameters of ammonium dihydrogen phosphate (ADP),
an order-disorder AFE, which always shatter at the transition temperature. The Curie-Weiss law, i. e.,
C
£ -s^--, is used to fit the value of ty, where s^, C1 and T1 are high frequency limit dielectric
T — Tx
constant, the Curie constant and the Curie-Weiss temperature, respectively (Ledzion et al. 2004). The assumptions here are that the value of ty is constant with respect to varying temperature and applied electric field.
For ADP, the Curie constant, the Curie-Weiss temperature and the Curie temperature are C1 = 10160 K, T1 = 22.7 K, and TC = 148 K, respectively (Ledzion et al. 2004; Milek, Neuberger 1972, 44). Below the Curie temperature, the fitting of ADP material constants with the Curie-Weiss law gives spontaneous sublattice polarization, P0« 0.081 C x m-2 and ty ~ 0.01233 (the corresponding n ~ 1.393 x 109). At T = 80K (t « -3.346 x 10-3), one unit of e is approximately equivalent to 1.1206 x 1010 V• m-1, and one unit of q is approximately equivalent to 0.1984 C x m-2.
Derivations of nonlinear dynamical equations
With the presence of time dependent electric field, the Lorentz force per unit volume exerting
8G
on the AFE system is given by , and the resulting Landau-Khalatnikov equation of motion from Newton's Second law is '
_d2Di _ dDi 8G /x
m-r^ + r—'- = -(p^= (7)
dz2 dz 8Dt
The first and second terms in Equation (7) represent acceleration and damping of the charge motion.
8G —
—= is the variation of AFE free energy density, G in Equation (4), with respect to quantity Dv
5Di _
in which Di represents the Maxwell and staggered displacements, i. e., Q and R. $ is the sublattice charge
density, m is the mass per unit charge, and Y is the mass per unit charge per second of the AFE system. For small damping, the dynamics of the Maxwell and staggered displacements are oscillatory.
I malu
In terms of dimensionless time variable s, and dimensionless damping constant g, defined as t = -5— s,
^ 11
Y=aI-S where m — ml(p and J — J '<P, and the previous dimensionless variables t, e, q, r, and ty,
V «in
the dimensionless oscillatory equations of motion are:
^l+g^. = -2(¥ + t)q + 4(q3 + 3qr2)-6(q5 + 10qV + 5qr4) + e (8a)
ds ds
^¡- + g— = -2tr + 4(3q2r + r3)-6(5q4r + 10q2r3+r5) (8b)
ds ds \ / \ /
Compare with the Duffing oscillator equations of motion (Goldstone, Garmire 1984), Equations (8) are more complicated with extra fifth order terms on the right side of the equations. Because of the existence of third and fifth order nonlinear terms, Equations (8) can be solved only numerically. These coupled equations are used to numerically simulate nonlinear and chaotic dynamics in the AFE system. In the numerical simulations, the fourth-order Runge-Kutta method is adopted (Press et al. 1996, 704; Strogatz 2018, 33), and is based on AFE material parameters and constants. The chosen antiferroelectric material is Ammonium Dihydrogen Phosphate (ADP) at temperature T = 80 K. For ADP, the dimensionless interaction constant between AFE sublattices and temperature are ty ~ 0.0123 and t « -3.346 x10-3, and dimensionless critical field is eC ~ 0.4166 (corresponding to critical applied electric to switch ADP from antiferroelectric to paraelectric states, i. e., EC « 4.668 x 109
V x m-1). The scaled natural frequency, a0, of ADP is obtained from Equation (10a) by dropping all the dissipating and nonlinear terms. For ADP, a>0 «¡/+t) «0.1341 per unit s(dimensionless time). This shows that the natural frequency is temperature dependent (Ledzion, Bondarczuk, Kucharczyk 2004).
Nonlinear and chaotic states of AFE from numerical simulations
The numerical simulations are carried out by applying the dimensionless time varying sinusoidal Maxwell field, e (corresponding to EM), to Equations (8). It is represented by e = e0 sin(as). As an example, this field is shown in Figure 1 for a = 1.0 a0, and e0 = 0.6 eC.
Fig. 1. e versus s with a = 1.0 a0, and e0 = 0.6 eC
The calculations are started with extremely low frequency and extremely small amplitude of e, i. e., a = 4.0 x 10-3 a0 and e0 = 1.0 x 10-13 eC. These values are close to the lowest limits computable by the numerical programming. The values of a and e0 are increases in the subsequent calculations to exhibit their effects on the AFE system. The value of damping is the same for all simulations, i. e., the value for reduced damping is g = 0.01. The numerical results for each applied field are plotted in four figures, (a) shows the reduced/dimensionless normal displacement, q, versus reduced time, s, (b) shows the reduced normal displacement, q, versus the reduced applied Maxwell field, e, (c) shows the phase diagram of the system, i. e.,
time derivative of reduced normal displacement, dq/ds, versus reduced normal displacement, q, and, (d) shows the Poincare Sections of the phase diagram started on 1/8 cycle for 200 cycles. The duration for numerical results in (a) to (c) are two cycles, as in Figure 1.
The numerical results for e0 = 1.0 x 1013 eC and a = 4.0 x 10-3 a0, are shown in Figures 2(a) to 2(d). At this small applied field amplitude and low frequency, the AFE system exhibits linear response to e, with the q versus s curve in Figure 2(a) is exactly the same pattern as the curve e versus s in Figure 1. This linearity is further affirmed in Figure 2(b), the curve of AFE state is a straight line through the origin of the q versus e axes. The phase diagram of the system is plotted in Figure 2(c). The elliptical closed curve shows that the attractor for the AFE system is at the origin of the amplitude of applied field and value of frequency. The regularity and periodicity in dynamical response of this AFE system under driven applied field are shown by the Poincare Sections plotted in Figure 2(d), in which 200 points are taken with 2n phase increases for each subsequent point, started from 1/8 cycle. Figure 2(d) shows almost an overlapping of these points.
tr
"U
-8.0E-18
c :> (C)
-l.E-14 0.E(00 l.E-14
q
Fig. 2. AFE's phase diagrams when a = 4.0 x 10 3 a0, and e0 = 1.0 x 10 13 eC 126 https://www.doi.org/10.33910/2687-153X-2022-3-3-122-136
For applied field amplitude and frequency increased to e0 = 1.0 x 10-5 eC and a = 0.1 a0, the numerical results are shown in Figures 3(a) to 3(d). At these values of applied field amplitude and frequency, the AFE system exhibits period multiplication centered on sinusoidal wave as shown in Figures 3(a). This period multiplication is further exhibit in Figure 3(b), with the wavy features centered on the line of AFE state. The corresponding phase diagram of the system is plotted in Figure 3(c). In the phase diagram, the period multiplication exhibits as multiple ellipses skewed to the right. The corresponding Poincare Sections are plotted in Figure 3(d). The Poincare Sections show an overlapping to only four points, this means the responses of the AFE system are quasi-periodic in these small values of applied field amplitude and frequency.
1.2E-14
q O.OEiOO
-1.2E-14
200
400
600
800 1,000
Fig. 3. AFE's phase diagrams when a = 0.1 a0, and e0 = 1.0 x 10-5 eC
With applied frequency maintains at a = 0.1 a0, the subsequent increases of applied field amplitude to (i) e0 = 0.3 eC, (ii) e0 = 0.95 eC, (iii) e0 = 1.0 eC, (iv) e0 = 100.0 eC, and, (v) e0 = 2.0 x 106 eC, are plotted in Figures 4 to 8. From the curves and Poincare Sections in Figures 4 to 7, the features of quasi-periodicity decrease as applied field amplitude increases.
Fig. 4. AFE phase diagrams when a = 1.0 a , and e0 = 0.3 e
Fig. 5. AFE phase diagrams when a = 1.0 a0, and e0 = 0.95 e
C
Fig. 6. AFE phase diagrams when a = 1.0 a0, and e0 = 1.0 eC
Fig. 7. AFE phase diagrams when a = 1.0 a0, and e0 = 100.0 eC
Fig. 8. AFE phase diagrams when a ■■
1.0 a0, and e0 ■■
: 2.0 X 106 e„
When the applied field amplitude is extremely high as in Figures 8, with e0 = 2.0 x 106 eC close to the upper limit computable by the numerical programming, the q versus s curve in Figure 8(a) is aperiodic, chaotic, and is highly distorted from sinusoidal shape, with little wavelike features in the first and third quadrants of each cycle. These features are more obvious in Figure 8(b). The curve of the AFE state, i. e., q versus e, shows wavelike features in between a double hysteresis loop. The corresponding phase diagram is plotted in Figure 8(c), which is also identified as butterfly curves of the AFE system. In the phase diagram, the curves are winding around two values of q, one is positive, the other one is negative, which are identified as a set of strange attractors (Strogatz 2018). The Poincare Sections in Figure 8(d) show one segment of line, this means the responses of the AFE system in subsequent cycles are very close but not overlapping.
With the applied field amplitude maintained at e0 = 0.6 eC, the subsequent increases of applied frequency to (i) a = 2.0 a0, and, (ii) a = 3.0 a0, are plotted in Figures 9 and 10, respectively. From the curves and Poincare Sections in Figures 9 and 10, the features of quasi-periodicity decrease as applied field frequency increases. This is obvious in the Poincare Sections, as Figure 10(d) shows more branches of points than Figure 9(d).
Fig. 9. AFE phase diagrams when a = 2.0 a0, and e0 = 0.6 eC
Fig. 10. AFE phase diagrams when a = 3.0 a0, and e0 = 0.6 eC Conclusion
From the results, i. e., the coupled oscillatory Equations of motion (8), and the numerically generated curves in Figures 2 to 10, the AFE system, i. e., ADP, in its first ordered phase exhibits high sensitivity to changes of initial conditions in frequencies and amplitudes of the applied electric field. This is due to the existence of nonlinear third and fifth order terms in Equations (8). As the amplitude of field slightly increases, the system exhibits quasi-periodicity with period multiplications on top of the applied frequencies. This feature is obvious in Figures 3 to 5, with the wavy structures on the slightly distorted patterns of q versus s, and q versus e curves. As the amplitude of applied field increases further, the system loses its periodicity, or the responses of the system become aperiodic, as shown in Figures 5 to 7. When the amplitude of applied field is extremely high, the system exhibits apriodic chaotic responses as shown in Figures 8. The main mechanism for the changes of responses of the AFE system from
periodic to quasi-periodic, to aperiodic, and to aperiodic chaotic, corresponding to the increases in amplitudes of the applied electric fields in Figures 2 to 8, is mainly due to repeated stretching and folding (Strogatz 2018) of the responses of the AFE system to the applied electric field.
The theoretical and numerical approaches here manage to unravel parts of the nonlinear and chaotic responses of the AFE systems, i. e., ADP, with respect to changes in frequency and amplitude of the applied electric field. The model described in this paper can be applied to other first order AFE which can be formulated by the same approaches. The results of nonlinear and chaotic responses using this approach can be further used to determine the transmission and reflection of optical waves through films and multilayers of nonlinear medium (Goldstone, Garmire 1984). This is an alternative to the conventional approach in derivations of nonlinear susceptibility tensors using series expansion of polarization in terms of electric fields, e. g. P = e0x(1) E + e0x(2) EE + e0x(3) EEE, and use the results in the derivations of linear and nonlinear dielectric functions (Murgan et al. 2002).
Conflict of Interest
The author declares that there is no conflict of interest, either existing or potential.
References
Chan, T. Y. (2010) Study of chaotic dynamics and hysteresis in bulk antiferromagnet and Antiferromagnetic Film.
MSc Thesis (Antiferromagnetism). George Town, Universiti Sains Malaysia, 125 p. (In English) Diestelhorst, M. (2003) What can we learn about ferroelectrics using methods of nonlinear dynamics? Condensed
Matter Physics, 6 (2), 189-196. https://doi.org/10.5488/CMP.6.2.189 (In English) Goldstone, J. A., Garmire, E. (1984) Intrinsic optical bistability in nonlinear media. Physical Review Letter, 53 (9),
910-913. https://doi.org/10.1103/PhysRevLett.53.910 (In English) Ledzion, R., Bondarczuk, K., Kucharczyk, W. (2004) Temperature dependence of the quadratic electrooptic effect andestimation of antipolarization of ADP. Crystal Research and Technology, 39 (2), 161-164. https://doi.org/10.1002/crat.200310165 (In English) Lines, M. E., Glass, A. M. (1977) Principles and applications offerroelectrics and related materials. Oxford: Clarendon Press, 664 p. (In English)
Milek, J. T., Neuberger, M. (1972) Handbook of Electronic Materials. Vol.8. Linear electrooptic modular materials.
New York: IFI/Plenum Publ., 264 p. (In English) Murgan, R., Tilley, D. R., Ishibashi, Y. et al. (2002) Calculation of nonlinear-susceptibility tensor components in ferroelectrics: Cubic, tetragonal, and rhombohedral symmetries. Journal of the Optical Society of America B, 19 (9), 2007-2021. https://doi.org/10.1364/J0SAB.19.002007 (In English) " Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. (1996) Numerical recipes in C: The art of scientific
computing. 2nd ed. Cambridge: Cambridge University Press, 537 p. (In English) Strogatz, S. H. (2018) Nonlinear dynamics and chaos with applications to Physics, Biology, Chemistry, and Engineering.
New York: CRC Press, 532 p. https://doi.org/10.1201/9780429492563 (In English) Tan, E. K. (2001) Static and dynamic properties of ferroelectric materials. MSc Thesis (Ferroelectricity). George
Town, Universiti Sains Malaysia, 164 p. (In English) Toh, P. L. (2009) Study of chaotic dynamics and hysteresis in bulk ferromagnet and ferromagnetic film based on yttrium iron garnet. MSc Thesis (Ferromagnetic Materials). George Town, Universiti Sains Malaysia, 76 p. (In English)
