Feasibility Study of Cooling a Bulk Acoustic Wave Resonator by Nanoparticle Enhanced Phase Change Material

: In the current study, the coupling of a cooling problem with the electromagnetic resonance of a bulk acoustic wave (BAW) material is investigated. As well, a new cooling method by the addition of nanoparticles to a phase change material surrounding the BAW resonator is presented. To solve the governing equations of piezoelectric charge and momentum balance, thermal balance, and ﬂuid ﬂow a code with the method of ﬁnite element is introduced. After validation of various features of the code with melting proﬁle, heat generation, charge curve, and dispersion curve with benchmarks, the eigenfrequency analysis of the system is done. The thermal behavior of the system at ﬁrst mode and various boundary conditions are studied. As well, the effect of nanoparticles in fastening the cooling of the BAW resonator is demonstrated.


Introduction
Nanotechnology could be a technique that searches about the behavior of materials on a nanometric scale to improve their abilities [1][2][3][4][5]. Nanomaterials have many industrial and engineering applications. One of the current challenges related to nanotechnology lies in how nanomaterials are produced [2]. In quickly creating a universe of IoT, the RFFE of shrewd gadgets should deal with higher information rates and access the full data transmission of 4G/5G remote innovation. The higher bandwidth and data rate of 5G requires smarter technology. The explanation behind this is the developing requests of pervasive low dormancy information (data should transfer fast) at higher working frequencies needed to oblige improved information transmission abilities and quickly develop quantities of clients [3]. Usual electro-acoustic resonators utilized for RF channel applications depend on two strategies for the excitation and proficient energy catching of an acoustic wave (1200 to 4000 m/s) proliferating in a piezoelectric material. The first strategy and technique depends on SAW, where the surface acoustic wave is created on a piezoelectric substrate by metal IDTs by all accounts. The second strategy and technique depends on BAW, where the mass acoustic wave in the bulk of the material is energized by the use of an electric field through terminals above and under a dainty piezoelectric plate [4]. Key execution markers and quality index for resonator configuration are frequency of resonance, quality factor, coupling coefficient, and impedance for RF channel applications [5].
A bulk acoustic wave (BAW) device is usually composed of a piezoelectric material (i.e., zinc oxide or aluminum nitride) between two electrodes [6][7][8][9][10]. The BAW devices are currently used in Wi-fi systems, wireless positioning systems, cell phones, and satellite navigation. The range of working frequency (10 8 Hz to 2 × 10 10 Hz) benefits the 6G technologies in small sizes (10 −6 m) [7]. For acoustic isolation of the BAW device from the surrounding or substrate, two methods usually are used: air cavity (free-standing) and (solidly mounted). The free-standing method is cheaper than using Bragg reflector [8]. Such acoustic mirrors usually used a series of coupled low and high acoustic impedance [9]. The thickness of the Bragg reflector is equal to the quarter wavelength for maximum acoustic reflectivity [10].
Many numerical methods are used in the layout design of the BAW devices to reduce the production costs and increase the accuracy of the device for requested resonance modes [11]. The network of BAW resonators makes a frequency filter that protects the electronic device from unwanted frequencies. When electrostatic discharge can occur in the system, in case of high power requests, lower insertion loss requests, high-frequency requests, ( f > 1. 5 2.5 × 10 9 Hz), or in case of steep stopband attenuation requests, the SAW device is replaced by a BAW device [12]. The design of a BAW device requires a scattering profile, impedance, stage, and Q-factor [13]. The scattering profile distinguishes all the potential modes, energy-related with every mode, the get over between the primary method of interest and other plate modes when the BAW resonator is energized. The impedance and stage contour measure the non-idealities in the BAW resonator plan and the general effect on the Q-values [14]. BAW resonator-based circuits require energy control in the piezoelectric layer dependent on the ideal choice of Bragg reflector layer, piezoelectric film, and terminal thickness. Figure 1 shows the BAW device geometry which characterizes the ideal resilience on the piezoelectric film and electrodes. The thickness of cathode and anode could be different [15]. That thicknesses will affect the frequency resonances of the system. Five significant decisions in the plan and advancement of high-quality acoustic wave resonator are numbered in the list as follows:
Pad between electrode and transmission line; 5.
Temperature expansion coefficient. The solution of governing equations is not straightforward. This is because of the intricacy of the piezoelectric conditions that depict the gadget and its stacking conditions, delivering scientific arrangements of the conditions for non-trivial 2D and 3D calculations [16]. The first challenge in the plan of thin-film resonators is the concealment of side resonances that can be energized around the recurrence of the ideal mode shapes. The BC applied by IDTs and electrodes can add many extra modes. Some tough layers are added at the end of fingers to control such effects. In any case, the anode measurements directed by electrical coordinating regularly are with the end goal that precise outcomes do not warrant 1D or 2D display. By 1D or 2D analysis, the effect of fingers and spurious modes or fake resonances cannot be detected [17]. Along these lines, one should examine the fake resonances through full 3D FEM calculations. Be that as it may, the solid advantage of FEM is in its natural ability to oblige convoluted calculations, various materials, piezoelectricity, and full precious stone anisotropy. Subsequently, just a few disentangling approximations are required, bringing about an exact and adaptable reproduction strategy [18].
Various nanoparticles such as Ag, TiO 2 , CuO, ZnO, Fe 3 O 4 , Indium, SiO 2 , SWCNTs, and MWCNTs, can be mixed by a base fluid throughout sonication to make some improvements in thermo-physical properties [9]. Usual PCMs in the industry are paraffin, wax, and poly-a-olefin [19]. In the process of sonication which takes from 30 min up to 24 h, some surfactants (such as ethanol and sodium dodecyl sulfate) are added for better mixing. In Figure 2, the geometry of the BAW resonator assisted by PCM box is presented. Figure 2. Geometry of the BAW resonator assisted by PCM box in 2D simulation.
As the thermal management and temperature compensation are important, the thermal design of BAW resonators is considered in its design [20]. Since the purpose of the current work is the use of NEPCM for cooling of BAW resonators in its design has significance in BAW development [21]. The main aim of the work is the investigation of nano-added PCM in thermal management of BAW which has the principal conclusions of durability and stability in electronic packaging application [22]. A novel design for PCM is a topic in BAW cooling [23]. Phase change of NEPCM could be used in BAW cooling [24].
Rational design calculations of surface acoustic wave devices were performed in many types of research [25]. However, the use of latent heat of phase change in a thermal energy storage systems especially PCM including nanoparticles for cooling on a smallscale is a new idea. In some batteries for temperature control purposes, the PCM and heat pipes are used throughout the discharge-charge cycle [26]. Property (temperaturedependent) variations in thermal energy storage and Brownian motion can affect freezing and melting problems in PCM. The big issues of thermal stability and compatibility within the context of materials are critical in latent heat storage systems. Before using for practical cases life-cycle assessment of phase change material is required. In such systems which are usually assisted by finned heat pipes (with highly conductive metal foams), hightemperature is a common problem [27]. Nowadays phase change materials are used as thermal energy storage for buildings. Thermal performance enhancement of NePCMs for solar energy heat exchangers and accumulators were tested. Since they can work with multiple phase change materials [28]. In this paper, a finite element method (FEM) code is developed for the mathematical arrangement of the electro-elastic conditions that oversee the directly constrained piezoelectric vibrations. Both methodologies, that of taking care of the field issue (consonant investigation) and that of taking care of the relating eigenvalue issue (modular examination), are portrayed. A FEM programming bundle has been made without any preparation. Significant angles fundamental to the proficient execution of FEM are clarified, for example, memory of the executives and tackling the summed up piezoelectric eigenvalue issue. Calculations for diminishing the necessary PC memory through enhancement of the framework profile, just as matrix calculation for the arrangement of the eigenfrequency issue are connected into the product from outer mathematical estimations. Current FEM programming is applied to to solve the mathematical governing equations of a thin-film mass BAW. Specifically, 3D reproductions are utilized to explore the impact of the top terminal shape. The legitimacy of the displayed strategy is shown by contrasting the recreated and estimated removal profiles at a few frequencies. The outcomes show that valuable data on the exhibition of the flimsy BAW acquired by generally large cross-sections and, subsequently, some assets. The novelty of the current study is that the NEPCM has not been used previously for cooling thin-film bulk acoustic wave resonance. The performance of that method of cooling and its side effects on electrical and mechanical parts are investigated in this paper.

Materials and Methods
In this section governing equations and boundary conditions of the system are presented.

Cooling Part
The governing equations of fluid, the non-slip boundary conditions of fluid on walls, and the boundary condition for thermal energy are presented in Tables 1-3 respectively. Enhancement of solidification rate of latent heat thermal energy storage using corrugated fins (shell-and-tube) was proved but here just a simple PCM cube is used. In some systems the PCM slurries' density changes at phase change temperature, but here a Newtonian model is used. The nanoparticles dispersed in a phase change material can also improve melting characteristics which is neglected here. Table 1. Governing equations of fluid [2].

Wall Equation
right ∂T(x=L) ∂x bottom Figure 1 shows the resonator configuration. The geometry of the BAW cell is equal to one of that circuit lines. Figure 2 illustrate the resonator and cooling part. The effective characteristics are shown in Table 4 while the material thermophysical properties are presented in Table 5. The effect of the addition of nanoparticle on the thermal behavior of PCM is obtainable from Table 4. Many fluids and solids such as organic materials, carb-oxy-methyl cellulose carbon foams, paraffin, and expanded graphite are used as phase change material previously for a thermal storage system. However, here TH29 was selected as [2,29].
Graphs of heat flow (in joule per second) as a function of temperature which is called in literature the differential scanning calorimeter (DSC), usually provided by PCM producers. TH29 (Calcium-chloride hexahydrate) performance in DSC form is presented in [29] for the pure material and its mixing with cellulose. As well, the TH29 shows good performance in aged cases which exposes to ambient air for two weeks. From the mass measurements, the hydrated salt-based TH29 absorbed moisture at about 50% of their weight and lost their ability in eight days.
As well, various nanoparticles were used such as Ag, carbon nanofiber, Al, aluminum, and carbon nanotubes. However, here copper nanoparticle was used. The governing equations and thermal boundary conditions for heat transfer in solid part are presented in Table 6. Table 4. Nano-fluid property formulas [4]. Table 5. Material thermo-physical properties [2]. Usually, the thermal delay method is used to determine the effective characteristics of micro-encapsulated phase change material. Table 5 mentions TH29, a commercially available PCM. Table 6. Governing equations and thermal boundary conditions for heat transfer in solid part [6].

Piezoelectric Induced Wave
The piezoelectric constitutive conditions that couple the mechanical and electrical amounts in the piezoelectric material are communicated utilizing the coupling equations. The conservation of momentum for the bulk acoustic wave in an electrical part is presented in Table 7. Table 7. Conservation of momentum for bulk acoustic wave in electrical part [6].

Zone Equation
periodic left and right walls u destination = u source e (−ik (r destination − r source)) , (26) top and bottom free surfaces T.n = 0, (27) The electrical boundary condition for the piezoelectric wave is depicted in Figure 3.
) which comes from conservation of momentum and presented in in Table 7. The coupled constitutive equations in strain-charge form [6] are and Consider next the case where the slab is electrically excited by applying a harmonic voltage, v = V 0 e jt , across the electrodes. In this case, it is convenient to consider a slab of length L with ends at x = L/2. The stresses at the ends are zero. The fields of interest are expressed by eigenfrequency analysis (−ρω 2 u = ∂T 11 ∂x ). The dielectric permittivity, elastic and piezoelectric constants of ZnO are considered in this study.
As the velocity of the electromagnetic wave is faster than the acoustic wave, the electrostatic equation of Laplace is applied for electric potential (∇ 2 φ = 0). As well, the electrical displacement is the function of strain and the electric field

Results
Here, predominantly, the instance of voltage stacking is considered, i.e., one terminal is electrically grounded, and a sinusoidal voltage with endorsed consistent plentifulness is associated between the grounded and "hot" anode. Arrangement of the eigenvalue issue at that point brings about arrangement reverberation frequencies (short-circuited terminals) at which the electrical impedance expects low qualities. The equal reverberation frequencies with enormous estimations of impedance can be acquired from the open-circuited case. For the open-circuit limit conditions, the potential is set consistent on the "hot" anode, yet not fixed (drifting potential).
In this section which describes the results, first, the confrontation of the obtained results with other research are presented through validations. Such a suitable comparison significantly raises the meaning of the presented paper. Then the effect of adding PCM to BAW is discussed.

Validation of Used Code
In this section, various parts of the code are validated against well-known benchmarks and commercial software.

Validation of Fluid Part
For the purpose of the validation of the fluid part, the previous work of Jamalabadi and Park [2] is used. The vertical axis in Figure 4 is solid content that has no unit. As revealed in Figure 4 the solid percent of PCM decreased by the increase of time (melting process). As seen in Figure 4 both methods are in good agreement.

Validation of Heat Generation
In this section, the validation of heat generation caused by resonance with the analytical solution is performed. The loss mechanisms in the thin film BAW resonators include leakage to the substrate, laterally escaping waves, ohmic losses, viscous losses, dielectric losses, and wave scattering. Additionally mentioned in the literature are eddy current losses and losses due to the formation of conducting channels in the high resistivity substrate interface. The analytical solution of the piezoelectric governing equation for heat generation is given in Appendix A. To calculate the heat generation of the electric source, both mechanical part (elastic energy) and electric part (electric displacement energy) should be considered. The Poynting vector (W/m 2 ) which is the directional energy flux (the energy transfer per unit area per unit time) usually used as a symbol of conservation of electromagnetic energy (S = E × H). The basic idea is that Poynting energy losses some parts as passed through the piezoelectric material and electrodes. Here the power loss distribution when the slab is electrically excited by an applied voltage is considered. The electrical excitation with the amplitude of 10 V is applied across the terminals. Material and geometrical parameters for the case of heat generation are presented in Figure 5 and Table 8.  The comparison of the current code with the analytical solution which is given in Appendix B is plotted in Figure 6. As seen in Figure 6 both numerical and analytical results are in good agreement. The validation results of Figure 6 has the error of simulation near zero as it designed to respond in 1D and the effect of transverse direction is not too much as the electrode extends all over the top and bottom surface. In engineering applications, the electrodes cover the surface partially and some materials are added to help the system to have these ideal mechanical support connections.

Validation of Charge Curve
In this section, the 3D problem is considered while in another part of the paper you can consider the 2D model as the third direction has not affected the results. Simulation of a conventional BAW resonator is performed in the appendix code presented at the end of the paper. After a variable definition on substrate, electro, piezo, and resonator parameters. Element types of solid226 and solid95 are used to model the piezoelectric parts. Materials parameters of aluminum nitride for various directions as anisotropic material in the form of tabular data are given. Some damping parameters are also considered for mechanical and electrical losses. Silicon and molybdenum data are entered as well. For building geometry, the bottom substrate, bottom electrode, piezo layers, and top electrode simple block are used. 3D regular mesh is created to apply the mechanical boundary conditions (two symmetry faces, right surface from top view symmetry, and up to surface from top view symmetry, down surface from top view hard boundary as a fixed plate) and electrical boundary conditions (electrode bottom, electrode top). The results of harmonic analysis in a range of 1100 and 1250 Mhz are used to validate the charge collected on the electrode surface. Material parameters definition and geometry of the compared case are presented in Tables 9 and 10, respectively. As seen the piezoelectricity and full anisotropy of the materials are taken into account. The comparison of the current code with the ANSYS-APDL code results which the code is provided in Appendix B is plotted in Figure 7. As seen in the figure, both results are in good agreement. The validation results of Figure 7 have the error of simulation near zero. As both used the same method (FEM) these results were anticipated.

Validation of Dispersion Curve
The behavior of the side resonances in a multilayer structure is a function of Lamb modes curves traveling in a multilayer structure is investigated in this section. The comparison of a dispersion curves presented by Makkonen et al. [3] is given in Figure 8. As depicted in the Figure 8, both methods are in good agreement.

Eigenfrequency Analysis
The first six modes of the BAW problem are presented in this part. The real part of the first six modes of BAW device is 221.4200 222. 02, 222.97, 224.16, 225.51, and 226.96 MHz. The imaginary parts are three orders of magnitude less than the real parts with the values of 83. 593, 84.982, 84.029, 90.701, 90.162, and 90.142 kHz. Detailed parameters are given in Table 11. To have a better understanding of given parameters they have been plotted in Figure 9. As shown, the coefficients have a fluctuation from one frequency to another. In participation factors of x and z, even frequencies (2, 4, 6, . . .) have higher values and y participation factors have maximum in odd frequencies (1, 3, 5, . . .).  Each of propagating modes and evanescent modes can also be visualized by its displacement field. As seen in Figure 10, the horizontal displacement of mode 1 has two peaks at the middle of terminals. Figure 11 presents four peaks in the vertical displacement of mode 1 where one of them is in opposite direction. Horizontal displacement of mode 2 as expected has four (2 × N) peaks which are symmetric around the middle of piezoelectric which is visible in Figure 12. As illustrated in Figure 13, vertical displacement of mode 2 has four maximum and four minimum which are aligned in three parallel lines. The presence of propagating and evanescent modes is clear. The exhibition of six (2 × N) simple maximum of mode 3 of horizontal displacement of mode 3 is demonstrated in Figure 14. As shown the eigenfunctions are aligned in three parallel lines and are symmetric conditions. Finally, the eigenfunctions of mode 3 of vertical displacement with six visible maximum are depicted in Figure 15. A detailed BAW design, general discussion on BAW resonators, and electric potential modes are discussed in [1,6]. In surface plotted figures including Figures 10-15 the color is automatically proportional to surface height.      The damping ration, quality factor, and S 11 are plotted in Figures 16-18 respectively. Damping ratio (ratio of imaginary of the frequency to its absolute value) and quality factor (ratio of absolute value to two times of the imaginary of the frequency) as a function of frequency are plotted in Figure 16. The quality factor for frequency as a function of frequency is depicted in Figure 17. The real part and the imaginary part of the scattering parameter as a function of frequency are illustrated in Figure 18 for the layer structure of a BAW resonator. Here the quality plot is against all selected intervals of frequency. The dispersion curve of imaginary (evanescent mode) and real (propagating mode) parts of the wavenumbers are figured in Figure 19 for a thin-film BAW composite resonator. The dispersion schematic includes both imaginary and real parts of the wavenumber. The plots actually are the contours of admittance |Y 11 | for real part of k, and contours of admittance |Y 22 | for imaginary the part of k. Figure 19 shows the calculated dispersion curve. The results for the the real and the imaginary parts of the wavenumber are combined by taking the positive k-axis with an optimum near (6 × 10 4 , 2 × 10 8 ) for the real part and the negative k-axis with a maximum near (−2 × 10 4 , 2 × 10 8 ) for the imaginary part. The horizontal axis of dispersion curve is k, wave number, with the unit of 1 m while the vertical axis is f , frequency, with aunit of 1 s . They correspond to the lower branch for the TS 2 mode (thickness-shear), the upper branch for the TE 1 mode (thickness-extension), and the left branch for the evanescent modes.    Total power dissipation as a function of frequency is plotted in Figure 20 for the specific device structure thin-film BAW composite resonator considered here. As shown, the maximum value of heat generation happened in the first mode while is in accordance with maximum admittance (Y 11 ). This graph will be used in the next part as the input for heat generation.

Fluid and Thermal Analysis
Profile development and contours of velocity magnitude distribution in the liquid part of the PCM cell are plotted in Figure 21. As shown, the liquid profile is developed fast near the initial condition as the absorbed heat in the PCM is converted to phase change in the liquid-solid boundary. Figure 21 does not have a color legend and is focused on showing the melting profile development. In the red zone, the velocity magnitude is maximum while in blue zones (near walls) the velocity magnitude is near zero. As time goes on, the melting process developed in the PCM box. As shown, the velocity is minimum near walls and is maximum at the middle of two flow channels created horizontally at the PCM cell. Detailed velocity distribution in PCM cell with cooling at the end of melting process is plotted in Figure 22. As shown in velocity vectors and magnitude in the PCM field, the maximum velocity happens at the middle of two flow channels with the horizontal direction. At the end of the flow stream, it turns vertically to turn inside by natural convection.  The power distribution inside BAW slab for first and second mode are plotted in Figure 23. The distribution is axially and matches the analytic solution given in Appendix A.
The temperature distribution inside BAW slab for first and second mode are plotted in Figure 24. The distribution is axially and for the case of constant temperature BC at both sides.  Temperature difference (∆T = T − T ∞ ) contours are plotted in Figure 25. Figure 25 does not have a color legend and is focused on showing the melting profile development. In the red zone, the velocity magnitude is maximum while in blue zones (near walls) the velocity magnitude is near zero. The figure shows the temperature profiles through the time for various BC. The initial condition of PCM cell is temperature equals the environmental temperature (T = T ∞ ). As time goes on, the cold zone is shortened through the volume, and the heated zone expands. If both sides are kept at environmental temperature the temperature profile inside the slab is parabolic. The case of constant temperature BC at both sides happens when the periodic structure of the BAW + PCM cell forced the melting temperature at both side until the end of phase change. The case of adiabatic BC at the left happens when each BAW + PCM cell is considered as an isolated domain from the other. Such conditions force the melting temperature on the right side until the end of phase change while adiabatic condition at the left of the cell. As shown, the periodic structure causes less temperature increase.
Temperature distribution inside BAW + PCM cell is illustrated in Figure 26. Time dependency of the temperature in the BAW + PCM cell shows that, by considering PCM, the thermal inertia happens in the system.  The effect of nanoparticle concentration on time dependency of maximum temperature is illustrated in Figure 27. By increasing the nanoparticle the rate of interface convection increased and faster phase change happens.

Conclusions
In this paper, a FEM code is introduced for the mathematical arrangement of the electro-elastic conditions that oversee the directly constrained vibrations of piezoelectric media. A consonant time reliance is accepted. Both of the methodologies, that of taking care of the field issue (consonant investigation) and that of taking care of the relating eigenvalue issue (modular examination), are portrayed. A FEM programming bundle has been made without any preparation. Significant angles fundamental to the proficient execution of FEM are clarified, for example, the memory of the executives and tackling the summed up piezoelectric eigenvalue issue. Calculations for diminishing the necessary PC memory through enhancement of the framework profile, just as Lanczos calculation for the arrangement of the eigenvalue issue are connected into the product from outer mathematical libraries. Current FEM programming is applied to nitty-gritty mathematical demonstrating of thin-film bulk acoustic wave (BAW) composite resonators. Correlation of results from 2D and full 3D recreations of a resonator is introduced. Specifically, 3D reproductions are utilized to explore the impact of the top terminal shape on the resonator electrical reaction. The legitimacy of the displaying strategy is shown by contrasting the recreated and estimated removal profiles at a few frequencies. The outcomes show that valuable data on the exhibition of the flimsy film resonators can be acquired even with generally coarse cross-sections and, subsequently, moderate computational assets. To highlight the advantages and disadvantages of using NEPCM for cooling application of BAW resonator it is good to notice that when comparing the current devised solution with other solutions from the scientific literature, the efficiency of the current design is more clear. While in common circuits the heat is removed from the semiconductor by conduction to support or with convection or radiative heat transfer modes to the environment, here the phase change plays the role of the heat sink with more thermal inertia. As well, the use of nanoparticles facilitates the heat removal, and also decreases the delay time and thermal inertia of the cooling system. As it would be useful to add information on further research of the authors related to the continuation of this research topic, future research could be conducted to consider the 3D design of a thin-film bulk acoustic wave composite resonators in conjunction with another electrical part to see the effectiveness of current study in an engineering case.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

2D
Two dimension 3D Three

Appendix A. Analytical Solution in One Dimension
The current model as presented in Figure 5 simulates a uniform layer structure extending horizontally to infinity. The computed dispersion diagram of such an infinite plate can be used to optimize device designs. For this tutorial, one can take a thin slice from the center of the complete device and use the periodic boundary condition to extend it laterally to infinity. If we apply φ(z = h) = φ m e jωt at the top electrode while bottom electrode is ground (φ(z = 0) = 0) then from Laplace equation in electrostatic approximation thus the electric field in z-direction is if we differentiate from the first equation of piezoelectric Considering the definition of strain from displacement field (S 11 = ∂u ∂x ) the following equation for the differential of the stress respect to x direction is obtained substituting in right hand side of the dynamic equation The above equation can solve by the free stress boundaries at the ends (T 11 (x = ±L/2) = 0) gives the displacement field as Considering the definition of strain from the displacement field Considering the definition of stress from the first piezoelectric equation Considering the definition of electric displacement from the second piezoelectric equation The analytical solution of electric displacement is obtained wheres 11 ,d 31 ,ε T 33 denote real and s 11 , d 31 , ε T 33 denote imaginary parts, respectively. An expression for the power dissipation density in the slab is obtained by considering hysteresis loops. Integrating to determine the areas under the loops and averaging over time yields: where Q E is