Patrick Fagan

Benjamin Ducharne

Stanislaw Zurek

Mathieu Domenjoud

Anastassios Skarlatos

Laurent Daniel

Christophe Reboud

Iterative Methods for Waveform Control in Magnetic Measurement Systems

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.



Abstract-Magnetic losses in a ferromagnetic lamination can be separated into three contributions. Bertotti theoretically assessed this distribution at the end of the 20th century in the Statistical Theory of Losses (STL), triggering significant progress in understanding the dissipation mechanisms. Recent studies have shown the possibility of reconstructing a hysteresis cycle from the high-frequency Barkhausen noise signal. Applying STL to the Barkhausen noise cycles has never been done before. Still, it could help establish a parallel with the measurement of the magnetization cycle versus frequency and the energy loss. However, STL analysis in its ultimate description requires sinusoidal flux density, while Barkhausen noise measurements are usually done with a constant excitation slope. Multiple magnetic flux density control methods were described in the literature and are reviewed in this manuscript. However, the Barkhausen noise context, requiring high-frequency sampling during the magnetization cycle, is more constraining. Therefore, specific performance criteria were considered, followed by numerical tests to determine the most adapted method to a Barkhausen STL description. Eventually, the Proportional Iterative Learning Control (P-ILC) gave the highest satisfaction rate and was chosen for experimental tests. Some of these experimental results are provided in the manuscript discussion together with suggestions for convergence speed improvement. It is, for instance, recommended to increase the gain near saturation, where the system response is poor.

Index Terms-Magnetic Barkhausen noise energy, waveform control, feedback, magnetic losses

I. INTRODUCTION

Magnetic cores are omnipresent in electrical energy conversion and transport. Losses inevitably happen while magnetic cores operate and are a significant cause of inefficiency. During one magnetization cycle, these losses are equivalent to the area of the Ba(Hsurf) hysteresis cycle obtained Patrick Fagan, Anastasios Skarlatos and Christophe Reboud are with LIST, CEA, CNRS, Université Paris-Saclay, CEA Saclay, 91191 Gif-sur-Yvette, France.

Benjamin Ducharne is with ELyTMaX IRL3757, CNRS, Univ Lyon, INSA Lyon, Centrale Lyon, Université Claude Bernard Lyon 1, Tohoku University, Sendai, Japan and with LGEF, Laboratoire de Génie Electrique et Ferroélectricité, INSA Lyon, France. (e-mail: benjamin.ducharne@insalyon.fr). by plotting the flux density Ba averaged through the tested specimen cross-section as a function of the tangent surface excitation field Hsurf [START_REF] Fiorillo | Measurement and characterization of magnetic materials[END_REF]:

P = ∮ B (H ) • dH (1)
The first attempts for the hysteresis losses prediction trace back to Steinmetz's model (Eq. 2) [2][3]. This empirical relation was limited to sinusoidal Ba conditions and stated that iron losses followed a fractional power "law" of both frequency f and B

(peak value of sinusoidal Ba):

P = k • 𝑓 • B (2)
where k, a, and b were three constant parameters depending on the nature and geometry of the magnetic circuit. Loss separation was initially proposed in 1924 with Jordan's work [START_REF] Jordan | Die ferromagnetischen Konstanten für schwache Wechselfelder[END_REF], who assumed that magnetic core losses could be divided (Eq. 3) in a static contribution Pst and a dynamic contribution related to classical eddy current loss Pcl [START_REF] Jordan | Die ferromagnetischen Konstanten für schwache Wechselfelder[END_REF]- [START_REF] Krall | Analysis and implementation of algorithms for calculation of iron losses for fractional horsepower electric motors[END_REF]:

P = P + P = α • 𝑓 + β • 𝑓 (3)
where α and β were fit parameters. This simple approach was later improved by adding an excess loss term to match the experimental data, wrongly approximated by Eq. 3 (grainoriented electrical steels FeSi GO, etc. [START_REF] Goss | [END_REF]). Initially, the physical justification for the additional contribution was unknown, and it was even referred to as "anomalous" loss. It was obtained by either modifying β or simply adding a third contribution or "excess" loss (Pexc) [START_REF] Shilling | Magnetic properties and domain structure in grain-oriented 3% Si-Fe[END_REF]: P = P + P + P (4) This empirical method found theoretical foundations from Bertotti Statistical Theory of Losses (STL) [START_REF] Bertotti | Hysteresis in magnetism: for physicists, materials scientists, and engineers[END_REF], expressed in terms of power in Eq. 5 or energy in Eq. 6 (ast, acl and aexc are fitting parameters):

P = a • 𝑓 + a • 𝑓 + a • 𝑓 (5)
W(𝑓) = = a + a • 𝑓 + a • 𝑓 (6)
STL is an advanced theoretical method that brought significant progress in understanding the magnetization mechanisms [START_REF] Bertotti | The prediction of power losses in soft magnetic materials[END_REF]- [START_REF] Zhao | Modified loss separation in FeSi laminations under arbitrary distorted flux[END_REF]. It is, however, worth mentioning that STL supposes a full flux penetration which restricts its domain Recent works [START_REF] Ducharne | Characterization and modeling of magnetic domain wall dynamics using reconstituted hysteresis loops from Barkhausen noise[END_REF]- [START_REF] Fagan | Multiscale modelling of the magnetic Barkhausen noise energy cycles[END_REF] have demonstrated the feasibility of plotting hysteresis cycles from Barkhausen noise measurement. For this, the so-called Magnetic Barkhausen Noise energy MBNenergy (Eq. 7) is plotted as a function of Hsurf.

MBN (t) = ν • ∫ sign • V (s) • ds (7)
ν is a normalization coefficient and VMBN the magnetic Barkhausen noise voltage drop across the sensor coil.

Comparably to the classic Ba(Hsurf), MBNenergy(Hsurf) hysteresis cycles depend on the excitation frequency and thus reflect physical properties. It is, therefore, tempting to apply STL to these alternative cycles and get additional insights regarding the physics of the magnetization mechanisms. MBNenergy is linked predominantly to the magnetic domain wall motions; thus, in the MBNenergy STL, the classical loss contribution, related to the macroscopic eddy currents and first-order-frequencydependent, should be negligible. Hence, the frequency dependency of the MBNenergy hysteresis cycle energy is expected to be written as in Eq. 8:

W (𝑓) = b + b • 𝑓 (8)
with bst and bexc, two constants. Eq. 8 is purely hypothetical and has never been validated by comparison to experimental results. By increasing the magnetization frequency, the available frequency band of the MBN spectral density is expected to shift upward. Still, this effect remains unclear and the comparison with STL is expected to bring clarification. Even if STL was originally developed, for the sake of simplicity, by assuming that the magnetization process occurs under controlled macroscopic constant induction derivative (triangular induction), the ensuing formulation was successively modified to comply with sinusoidal and generic induction waveform. Eq. 6 supposes this condition to be respected [START_REF] Ragusa | Loss decomposition in non-oriented steel sheets: the role of the classical losses[END_REF]. In the case of Barkhausen noise measurements, either the excitation current or the magnetic excitation Hsurf is usually imposed triangular [START_REF] Stupakov | A system for controllable magnetic measurements of hysteresis and Barkhausen noise[END_REF][22] and measurements rarely done under sinusoidal flux density. Ba(Hsurf) and MBNenergy(Hsurf) can be obtained from the same experimental setup. Therefore, a suitable method for the flux density control during Barkhausen noise measurement can be inspired by published work related to standard hysteresis characterizations [START_REF] Grote | Measurement of the magnetic properties of soft magnetic material using digital real-time current control[END_REF]- [START_REF] White | A multichannel magnetic flux controller for periodic magnetizing conditions[END_REF]. Still, additional constraints owed to Barkhausen noise measurement have to be considered before setting the most adapted method. A magnetic characterization setup contains two non-linear elements (the inductor yoke used to drive the magnetic field Hsurf and the tested sample). Analytical solutions are sometimes proposed for setting the induction control system. Still, the effect of hysteresis and the to-be-measured properties of the specimen makes them approximative hence iterative methods prevail. The focus of this paper is to examine a wide range of digital feedback methods found in the literature to assure a sinusoidal flux density during hysteresis measurements and select the most appropriate one for the specific context of the Barkhausen noise STL application. Comparisons exist in the literature [START_REF] Baranowski | Comparison of digital methods of the control of flux density shape[END_REF], but they are limited to two or three methods and never deal with the specific Barkhausen noise perspective. The study is restricted to digital feedback methods [START_REF] Zurek | Practical implementation of universal digital feedback for characterisation of soft magnetic materials under controlled AC waveforms[END_REF]- [START_REF] White | A multichannel magnetic flux controller for periodic magnetizing conditions[END_REF]. Even if widespread before the proliferation of computers, the analog feedback methods [START_REF] White | A multichannel magnetic flux controller for periodic magnetizing conditions[END_REF][32] [START_REF] Pólik | Measuring and control the hysteresis loop by using analog and digital integrators[END_REF] are poorly tunable and less robust since they rely on discrete components whose values can be challenging to set, especially if heated. Analog systems work in real-time. Perturbations cannot be anticipated, and high precisions cannot be reached on a wide range of frequencies and materials. Other problems can arise if the system is strongly non-linear and has unstable feedback [20]. It should also be noted that high-gain high-bandwidth analog (real-time) feedback can suppress large-amplitude Barkhausen noise activity [START_REF] Fiorillo | Measurement and characterization of magnetic materials[END_REF][20]. For the same reasons, this paper will not consider hybrid methods (obtained by combining digital and analog feedback methods [34][35]).

II. REVIEW OF ITERATIVE FEEDBACK METHODS

A. Performance criteria

This study aims to define and test the most efficient magnetic flux density control method in the context of the MBNenergy hysteresis cycle characterization and STL application [START_REF] Stupakov | A system for controllable magnetic measurements of hysteresis and Barkhausen noise[END_REF][22] [START_REF] Patel | A new adaptive automated feedback system for Barkhausen signal measurement[END_REF]. To obtain adapted comparisons and reach our objective, a specific series of criteria has been defined:  Number of iterations: a reduced number of iterations is important, especially in the low-frequency range where a measure can take several minutes and generate large data files complex to process. It is also critical in the high-frequency range, where thermal transfers due to the magnetic losses can affect the experimental conditions and the magnetic response of the material.  Accuracy: Convergence should be reached with a minimum error. Error estimations can take different forms, including relative Euclidean difference, form factor difference, Pearson dissimilarity, total harmonic distortion, etc.  Number of parameters: Feedback control parameters must be tuned for each new experimental situation. Optimizing a large number of parameters requires a lot of experimental data.  Robustness: the feedback method should remain undisturbed by external stimuli, including white noises, drifts, and offsets.  Memory allocation and computation time efficiency: The ideal feedback method computes an iteration with reduced time and limited memory capacity. This criterion is especially detrimental for techniques based on square matrix inversion requiring high computation capacity and large memory allocation.  Universality: This criterion is related to the capability of providing satisfactory results in different experimental conditions without extensive calibrations processes. A feedback method providing a positive answer to all the criteria listed above does not exist. Many approaches have been described in the literature [START_REF] Stupakov | A system for controllable magnetic measurements of hysteresis and Barkhausen noise[END_REF]- [START_REF] Matsubara | Acceleration technique of waveform control for single sheet tester[END_REF], and each technique can perform well in specific conditions. This study aims to compare these methods to find the most suitable feedback technique in the particular context of the MBNenergy(Hsurf) hysteresis cycles characterization.

B. Detailed description of iterative feedback methods

Let us introduce the feedback notation and a general feedback scheme (Fig. 1)  YG(t) is the reference, i.e. the ideal desired goal output (at time t).  yM(t,j) is the measured output of iteration j (at time t).  ϵ(t,j)= yG(t) -yM(t,j) is the error of j-th iteration (at time t).  x(t,j) is the system input of j-th iteration (at time t).

1) Iterative Learning Control (ILC) A straightforward iterative method for the control of a nonlinear system can be derived from the classical real-time Proportional Integral Derivative (PID) technique:

x(t) = K • ϵ(t) + K • ∫ ϵ(s) ds + K • (9)
with Kp, KI, and KD, the proportional, integral, and derivative gains. The iterative PID method has been described by several authors, including Gruebler et al. in [START_REF] Gruebler | Losssurface-based iron loss prediction for fractional horsepower electric motor design[END_REF], and consists in:

x(t, j + 1) = x(t, j) + Δx(t, j) (10)
Δx(t, j) = K • ϵ(t, j) + K • ∫ ϵ(s, j) ds + K • (t, j) (11)
In their simplest form (proportional correction only), the above equations can be simplified, which leads to the Proportional-Iterative Learning Control formulation P-ILC [START_REF] Gruebler | Losssurface-based iron loss prediction for fractional horsepower electric motor design[END_REF]:

x(t, j + 1) = x(t, j) + K • ϵ(t, j) (12)
The phase-lead ILC method is similar to P-ILC but involves the addition of a constant delay τ in the error term [START_REF] Freeman | Experimental evaluation of iterative learning control algorithms for non-minimum phase plants[END_REF][39]:

x(t, j + 1) = x(t, j) + K • ϵ(t + τ, j) (13) P-ILC is simple; the inputs are reduced to ϵ(t,j), and parameters to Kp. Its implementation is very straightforward, and, like classic PID, it can be very robust with the right choice of Kp. However, the choice between high gain/fast convergence and small gain/no divergence makes the optimization tricky, typically ending with a slower convergence speed at the benefit of better robustness.

2) Fourier Series Proportional-Iterative Learning

Control (FSP-ILC) Switching from the time domain to the frequency domain can be highly beneficial by simplifying mathematical operations. The Fourier transform being linear, Eq. 14 becomes 15:

x(t, j + 1) = x(t, j) + K • [y (t) -y (t, j)] (14) X(𝑓, j + 1) = X(𝑓, j) + K • [Y (f) -Y (𝑓, j)]
(15) Where X(f,j) is the Fourier transform of x(t,j). Fourier transforms lead to complex numbers, and FSP-ILC works with complex number formalism. While Fourier Distribution Iterative Learning Control FDP-ILC would apply Eq. 15 to the whole frequency spectrum, FSP-ILC is limited only to the excitation frequency's multiples (fexc) and can even be reduced to those of substantial contribution. Eq. 16 gives FSP-ILC's equation when k

∈ [1 -M]: X(k ⋅ 𝑓 , j + 1) = X(k ⋅ 𝑓 , j) + K • [Y (k ⋅ 𝑓 ) -Y (k ⋅ 𝑓 , j)]
(16) FSP-ILC performance is relatively close to those of P-ILC. The main advantage of the former is its ability to ignore all the high-frequency components (> M•fexc), including white noise and power source oscillations, especially when fexc is low. Its main drawback is the two Fourier transforms and the sum of complex numbers required per iteration, slowing down the control speed and inducing limitations in the low-frequency range. The computation time can be reduced significantly if Fast Fourier Transform (FFT) is employed, but at the additional restriction that the number of samples in the waveform becomes a power of 2.

3) Phase Correction by determination of Measure-Goal Delay (PhC-MGD) P-ILC performances depend on the power supply dynamic response and capability to generate x without undesired phase lag. A phase correction is sometimes needed to increase the feedback performance. For this, several methods have been proposed, including PhC-MGD, a technique described by Stupakov et al. in [START_REF] Stupakov | A system for controllable magnetic measurements of hysteresis and Barkhausen noise[END_REF]. In this method, a preliminary step increases x until yM reaches the required amplitude. This operation is achieved by a P-ILC correction of gain modulated by yM's amplitude identified at the previous stage (Eq. 17):

x(t, j + 1) = x(t, j) • 1 + •[ ( ) ( ( ))] ( ) (17)
Once yM reaches the required threshold, x is recalculated through a sum of two weighted contributions (Gph is the weight), the corrected phase contribution xph and the corrected amplitude contribution xampl:

x(t, j + 1) = G • x (t, j + 1) + 1 -G • x (t, j + 1) (18)
where xph is calculated as follows:  the measure is normalized according to the targeted waveform:

y = y • ( ) ( ) (19)
 the resulting signal is divided into sections where the targeted waveform is monotonic,  the delay ø(t) (see Fig. 2) between the targeted waveform and the measurement is estimated.  ø(t) is applied simultaneously to yM and x and leads to the delayed versions yreq and xph:

x (t, j + 1) = x(t + ø(t), j) y (t, j + 1) = y (t + ø(t), j)

Fig. 2 -PhC-MGD delay illustration.

Xampl is obtained from a P-ILC method of yreq targeted waveform instead of the usual yG:

x (t, j + 1) = x(t, j) + K • y (t, j) -y (t, j) (

In [START_REF] Stupakov | A system for controllable magnetic measurements of hysteresis and Barkhausen noise[END_REF], Gph has a non-zero value exclusively in the highfrequency range when phase shift exists. The main drawback of PhC-MGD is the normalization step, which flattens the measurement signal in the presence of undesired noise peaks and makes the phase delay estimation complicated. This noise issue has been solved partially in [START_REF] Stupakov | A system for controllable magnetic measurements of hysteresis and Barkhausen noise[END_REF] by taking one thousand points per cycle, interpolating Hsurf with a cubic spline, and smoothing Ba with a numerical filter.

4) Phase Correction by P-ILC on angles (PhC-P-ILC)

In 2005, Zurek et al. in [START_REF] Zurek | Use of novel adaptive digital feedback for magnetic measurements under controlled magnetizing conditions[END_REF] proposed a correction method described as the combination of two steps:  An amplitude correction leading to x * (an intermediary variable) and assimilated to a P-ILC of effective gain inversely proportional to yM's amplitude:

x * (t, j + 1) = x(t, j) + K ( , ) ( ) [START_REF] Stupakov | Analog of the induction law for the magnetic Barkhausen noise[END_REF] x * (t, j + 1) = x(t, j) +

( ( )) [y (t) -y (t, j)] (23)
 A phase correction applied through a phase delay function F:

x(t, j + 1) = F[x * (t, j + 1)] (24)
During this phase correction, every involved quantity is expressed by its Fourier series (Eq. 25 to 28):

y (t, j) = a + ∑ a (j) sin[2π k𝑓 t + ϕ (j)] (25)
y (t) = b + ∑ b sin[2π k𝑓 t + ψ ] (26)
x * (t, j) = c + ∑ c (j) sin[2π k𝑓 t + θ (j)]

I(t, j) = d + ∑ d (j) sin [2π k𝑓 t + β (j)] ( (27)
) 28
where M is the number of considered harmonics, 𝑓 is the targeted waveform frequency, ai, bi, ci, and di are the Fourier amplitudes, and I the output current. I and x are linked through Eq. 29. They are identical if the source is an ideal unity gain amplifier (Gsource = 1):

I(t, j) = G (s) • x(t, j) (29)
øi, ψi, θi and βi, are the Fourier phase lags associated to ai, bi, ci, and di. With this formalism:

x(t, j + 1) = F[x * (t, j + 1)] = ∑ c (j + 1) • sin[2π kf t + θ (j + 1) + α (j + 1)] (30) and, α (j + 1) = α (j) + K • (ψ -β (j)) (31
) αk is the k th harmonic applied phase lag, and Kα a proportional gain. PhC-P-ILC reduces the influence of the power source but works in the Fourier domain, which means time-consuming direct and inverse transformations (especially when the number of considered data points is large). The implementation is complex compared to P-ILC. 5) Non-linear correction with a Quasi-Newtonian Method (QNM) In 2008, Yamamoto et al. in [START_REF] Yamamoto | Waveform control for magnetic testers using a quasi-Newton method[END_REF] described QNM, a nonlinear iterative control method derived from the BFGS-like technique published by Li and Fukushima [START_REF] Li | A derivative-free line search and global convergence of Broyden-like method for non-linear equations[END_REF] a few years before. The BFGS algorithm (Broyden-Fletcher-Goldfarb-Shanno) can be described as a line search optimization method, and Li's method is the derivative-free version of this algorithm.

At each iteration, several variables need to be evaluated: an N•N Hessian matrix Bj (where N is the size of vector x), yM, yG, and ϵM. ϵ (t, j) = f[x(t, j)] (32) B0 (initialization at j = 0) is defined as an identity matrix, and f a function to minimize. For each iteration j, the following algorithm is run:  the search direction pj is found by solving:

B • p = -ϵ (33)
 the step-size αj is set to 1, a new measure is made with:

x = x + α • p (34)
which leads to obtaining yMtemp and ϵ .  if the error is not small enough ϵ

> ρ • ϵ - σ • p the previous step is repeated, but αj is contracted by a factor β ∈ [0 -1]: α = β • α ( 35
)
This step is repeated until ϵ

> (1 + η ) • ϵ -σ • α • p
 Bj is updated according to Eq. 36 below:

B = B + θ • • • x - x ( 36
)
where θj is set to ensure Bj is not singular.  xj, and yj are set from the temporary variables: x = x , y = y QNM differs from [START_REF] Li | A derivative-free line search and global convergence of Broyden-like method for non-linear equations[END_REF] by a different initialization of αj:

α = (37)
QNM converges faster than P-ILC, especially near saturation, where P-ILC requires hundreds of iterations. But, QNM suffers from several weaknesses, including the excessive number of parameters: the fixed parameters ρ, β, θ , σ1, σ2 and the variable ones such as η and θ modified for every iteration:

∑ η < η < ∞ (38) θ -1 < θ ( 39
)
Another drawback is the computationally expensive inversion of square matrix. Convergence speed is also impacted by the multiple measurements needed per iteration. Hence, QNM is inadequate in the very-low-frequency range. 6) Least-Square Method for a Hsurf(Ba) Polynomial Identification (LSM-PI) In 2011, Anderson [START_REF] Anderson | Constant dB/dt DC characterisation through digital control of magnetic field[END_REF] proposed LSM-PI, an alternative iterative method which can be summarized as follows:  Like in PhC-MGD, the data set y is divided into monotonic sections, and an offset is applied to obtain x(yM=0)=0 for each section.

 X(yM) is defined for every section and approximated by a high-order polynomial (up to order of 30 in [START_REF] Zurek | Use of novel adaptive digital feedback for magnetic measurements under controlled magnetizing conditions[END_REF]):

x(y ) ≈ ∑ a • y (t) (40)  ai coefficients are determined with a least-squares method, and a phase term is considered for the hysteretic behavior.  x is calculated from Eq. 40 by replacing yM by yG (Eq. 41):

x(t) = ∑ a • y (t) (41
) This method converges with a minimal number of iterations (just three as claimed in [START_REF] Anderson | Constant dB/dt DC characterisation through digital control of magnetic field[END_REF]), but its performance relies heavily on the order of the polynomial function. Many oscillation issues are noticed (especially near saturation, where the magnetic permeability µ is low, but the Hsurf(Ba) slope is high). The leastsquares optimization is computationally expensive, especially if many sampling points are considered. This problem can be lessened by expressing Hsurf(Ba) in a different orthonormal system. Moreover, this method assumes that Hsurf(Ba) is bijective, which is not the case if the maxima of Hsurf and Ba are not simultaneous (as it is in the high-frequency range). In that case, a phase delay must be considered to avoid wrong results from the least-square optimization.

7) Other methods

It is not possible to provide an exhaustive list of all feedback methods and their modifications described in the scientific and technical literature. The main techniques introduced in the sections above have been numerically implemented and tested in this study. They have been chosen for their singularities and originalities, but more methods exist, and even if not detailed nor tested, they are worth mentioning in this manuscript:  In 2016, Zhang et al. in [START_REF] Zhang | Measurement of threedimensional magnetic properties with feedback control and harmonic compensation[END_REF] described a proportional corrector, working in the frequency domain and in which both magnitudes and phases are corrected. For every harmonic, the correction can be written as:

X = K • Y -Y + K • ∫ Y - Y dt (42)
X = K • Y -Y + K • ∫(Y -Y )dt (43)
Then, x(t,j+1) is written as a Fourier series thanks to the Xmag and Xph coefficients. This method shares the same strengths and weaknesses as FSP-ILC but is also very sensitive to the nonlinear behavior of the ferromagnetic sample. In [START_REF] Zhang | Measurement of threedimensional magnetic properties with feedback control and harmonic compensation[END_REF], this issue is solved by correcting the calculated phases based on a lookup table. Unfortunately, no details are provided about the method for constructing such a look-up table.  In [START_REF] White | Control of flux in magnetic circuits for Barkhausen noise measurements[END_REF], White et al. use a proportional derivative PD-ILC method to control the excitation current I (assuming that Hsurf is proportional to I and the resistances and inductances values perfectly known). Good results are obtained, but compared to P-ILC the implementation is complex and requires detailed knowledge of the experimental conditions and their evolution during the test, which is not trivial considering that the inductance varies significantly with the level of excitation.  In [START_REF] Bosack | Closed loop control of hysteretic magnetization[END_REF], Bosack et al. start from Jiles-Atherton's model and assumes the magnetization M can be written as:

= g(H , M, t) + f(H , M, u) (44)
where f and g are two known functions, and u is a control variable, defined by:

= + = + u (45)
H0 is the ambient field excitation, and Hc is the corrected contribution. The resolution of the system gives:

u = - -K • f(M) • (M -M ) (46)
Eq. 46 looks like a Kp proportional correction. No details about the practical implementation are given in [START_REF] Bosack | Closed loop control of hysteretic magnetization[END_REF], except the use of a real-time PID corrector. The estimations of f and g rely on identification steps, and calibration must be made each time the whole system changes, which can be time-consuming, especially if recalibration is to be applied at each new measurement frequency. Finally, non-linear algorithms have also been described in contexts unrelated to magnetic waveform control (electrohydraulic molding machine in [START_REF] Havlicsek | Non-linear control of an electrohydraulic injection molding machine via iterative adaptive learning[END_REF] or lithographic apparatus in [START_REF] Heertjes | Non-linear iterative learning control with applications to lithographic machinery[END_REF]). Like QNM, these methods require matrices inversion, limiting the experimental sampling rate and leading to feedback control incompatible with MBNenergy characterizations.

C. Required precision criterion

All the methods described in this section have been developed to comply with international magnetic characterization standards, and different criteria have been proposed for their validation. It is worth noting that some of these criteria apply to the time derivative z of the targeted waveform x.

These criteria include:  The relative Euclidean Difference :

d (x, y) = ∫[ ( ) ( )] ∫ ( ) (47)
 The Form Factor (applied only to z. It is worth noting that z criteria are particularly difficult to meet, as any minor distortion in x gets amplified due to the derivative):

FFD(z, y′) = |𝐹𝐹(𝑧) -𝐹𝐹(𝑦′)| = ( ) (| |) - ( ) (| |) (48)  The Pearson coefficient: d y) = ∫[ ( ) ( )][ ( ) ( )] ∫[ ( ) ( )] ∫[ ( ) ( )] (49)
 The Total Harmonic Distortion (applied only to z):

THD(z) = ∑ (50)
Tab. 1 -Accuracy criteria of the methods described in this section.

Table 1 provides the target values as applied in the literature. Even if different waveforms can have the same form factor [20] [START_REF] Yamamoto | Waveform control for magnetic testers using a quasi-Newton method[END_REF], the IEEE standards recommend the use of this criterion for the magnetic hysteresis and losses characterization [START_REF]magnetic properties of electrical steel strip and sheet by means of an Epstein frame[END_REF][50]. It is unsurprising to find it applied in many studies. For the MBNenergy characterization, we found it relevant to apply the following criteria:

 Relative Euclidean Difference < 0.5%  z Form Factor = Fgoal ± 0.5%  Pearson coefficient > 1 -3 • 10 -5  z THD < 0.5% III. NUMERICAL IMPLEMENTATION
In our quest toward the "best" iterative method for MBNenergy hysteresis cycles, characterization, and STL application, all the techniques described in section 2 have been numerically implemented using MATLAB®. A sigmoid-type anhysteretic behavior (Eq. 52) has been used to simulate the material's answer:

y (t, j) = f[x(t, j)] = arctan [x(t, j)] (52)
Eq. ( 52) is convenient as saturation is taken into account, and xG can be expressed analytically:

x (t) = tan y (t) ( 53
)
The objective is to find xG leading to a sinusoidal yG. A preliminary test consists in plotting the spectral content of xG as a function of yG's amplitude (Fig. 3). Fig. 3 -xG harmonic content vs. yG amplitude. When yG's amplitude is large, high amplitude harmonics are generated, triggering issues if the power source dynamic performance is limited. Waveform control is easier at low amplitude (no saturation and quasi-linear material behavior). THD of xG can reach 0% at very low amplitudes, it however increases up to 18% at 0.75•max(yG) and even 51% at 0.95•max(yG). THD values exceeding 100% are possible if no control is applied [20], and this is expected to occur for even deeper saturation. No noise has been considered in all the following tests. The power amplifier is supposed to be ideal (infinite bandwidth, etc.), with perfect impedance matching.

The sampling frequency has been reduced to 500 Hz to limit the memory allocation and reach convergence even with QNM.

A. Shallow saturation

The tested methods are firstly compared in a shallow saturation case: 𝑦 (𝑡) = 0.75 sin(2𝜋𝑡) (54) The simulations are stopped when the Relative Euclidean Difference or the Form Factor Difference falls below 10 -10 . The maximum iteration number is set to 600. On the one hand, such low error is only achievable in simulation; experimental conditions are affected by white noise, drifts, etc. On the other hand, such high accuracy allows testing the methods with yM extremely close to yG. Table 2 gives the simulation parameters, and Fig. 4 Convergence is obtained for all the methods tested. LSM-PI and PhC-MGD are the fastest, with approximately 10 iterations. Still, for both these methods, the amplitude correction step requires a lot of intermediary measurements, which can be problematic in the low-frequency range. QNM converges after Table 3 "calculation time" only considers the waveform identification computation time, i. e., it does not include additional times associated with virtual measurement simulation (equivalent to the measurement time in the experimental setup).

B. Deep saturation

In the next test, the iteration methods are tested closer to a fully saturated configuration, where larger non-linearity is present:

𝑦 (𝑡) = 0.95 sin(2𝜋𝑡) (55)
Fig. 5 -Relative Euclidean and form factor differences for the deep saturation test.

The maximum iteration number is raised to 1500 since overall convergence is slower in this case. Fig. 5 gives the simulation results, and Table . 4 the simulation parameters.

Tab. 4 -Deep saturation simulation parameters

Again, QNM and PhC-MGD show the fastest convergence speed. PhC-P-ILC is also very fast, outclassing QNM, exhibiting temporary convergence errors. P-ILC converges slowly for this test, but unlike PhC-P-ILC, it never gets stuck on a precision plateau. Table 5 compares the performances. The precision criterion is set to < 10 -7 .

Tab. 5 -Performance comparison for the deep saturation test. A high number of iterations limits LSM-PI and PhC-MGD performances. QNM's iterations number is lower, but each requires a significant calculation time. The PhC-P-ILC method converges with a reduced number of iterations, eight times lower than P-ILC, but the calculation time for the latter is extremely short.

C. Overall simulation results

Table 6 compiles the comparisons based on the performance criteria defined in sub-section 2.1.

Tab. 6 -Performances comparison based on the performance criteria defined in sub-section 2.1. MBNenergy measurements require high sampling frequency, up to several hundreds of kHz, and long-time measurements leading to huge memory size for numerical feedback variables. A pre-treatment on Ba is required to reach convergence.

Oppositely, P-ILC does not require a perfect drift compensation to reach convergence. After all the numerical tests performed in this study, and based on Table 6 analysis, P-ILC appears to be the most adapted method in the context of the MBNenergy(Hsurf) hysteresis cycles characterization. In the next section, improvements are proposed for even better and faster convergence.

IV. P-ILC EXTENSIONS P-ILC is an excellent method for magnetization control digital feedback. P-ILC is simple to implement and tune. It is robust and fast. Its only limitation comes from the convergence speed, especially near saturation when the permeability falls and where a weak variation of Hsurf generates an even lower variation of Ba. This problem can be partially solved by increasing the proportional gain, but the response will diverge in the high permeability zones. A better solution consists in modulating Kp according to the system answer:

x(t, j + 1) = x(t, j) + K (t, j)ϵ(t, j) (56) This method requires additional parameters and should be considered with special attention. It can be implemented from an error array based on the two previous iterations (P-ILC-2 [START_REF] Zheng | A novel iterative learning control method and control system design for active magnetic bearing rotor imbalance of primary helium circulator in high-temperature gas-cooled reactor[END_REF]) and give Eq. 57:

x(t, j + 2) = x(t, j + 1) + K ϵ(t, j + 1) + K ϵ(t, j) (

The resulting error becomes a weighted sum of j+1 and j errors iterations. A generalized version (P-ILC-N) considering all the previous state N can even be written by extending Eq. 57:

x(t, j + N) = x(t, j + N -1) + ∑ K ϵ(t, j + N -q) (58) P-ILC-N convergence is faster. It is also more robust than standard P-ILC [START_REF] Zheng | A novel iterative learning control method and control system design for active magnetic bearing rotor imbalance of primary helium circulator in high-temperature gas-cooled reactor[END_REF]. But a minimum of N measurements are necessary for the corrector to be fully working. All K coefficients need to be optimized individually, which can be complex and demanding in experimental data. Hence, N should be kept as small as possible unless a reliable model is available for simulation. Another possibility consists of structuring the P-ILC iterative law as a Taylor approximation (assuming x is a yM smooth function):

x(y ) = x(y ) + A(y )(y -y )

If x(t,j+1)= xG(t), higher orders Taylor approximation gives P-ILC-TA, Eq. 60:

x(t, j + 1) = x(t, j) + ∑ K ϵ(t, j)

It is also possible to replace A in Eq. 59 with its optimal value, as obtained by the Taylor's approximation (P-ILC-TD):

𝑥(𝑦 ) = 𝑥(𝑦 ) + (𝑦 ) (𝑦 -𝑦 ) ⇒ = ∝ (61)
⇒ 𝐾 (𝑡, 𝑗) = (𝑡, 𝑗)

P-ILC-TD can reach high-speed convergence rates. Kp being inversely proportional to the system reactivity, the correction will be significant when dyM/dx is small. However, relying on derivatives, P-ILC-TD requires exact measurement, no noise, delays, or bandwidth limitations. Otherwise, this method diverges very quickly. Finally, the deep saturation test (subsection 3.2) was repeated, and all P-ILC new variants were tested. Fig. 6 shows the simulation results, and Table 7 shows the corresponding iteration numbers. All alternative methods converge faster than P-ILC. P-ILC-TD outclasses all the proposed methods. An iterative process close to P-ILC-TD robust enough to handle white noise would be by far the most indicated method.

V. EXPERIMENTAL SETUP AND IMPLEMENTATION

The experimental setup used for the Barkhausen noise characterization has been precisely described in [START_REF] Fagan | Multiscale modelling of the magnetic Barkhausen noise energy cycles[END_REF]. The excitation of the sample is based on single C-yoke with the magnetizing winding, and with the sensors attached to the sample under test. An overall 2D view of this experimental setup is depicted in Fig. 7. Fig. 7 -Overall 2D view of the Barkhausen noise experimental setup. The power amplifier was a Kepco BOP100-10MG. The excitation coil was made out of 10 turns. The studied specimens were all grain-oriented electrical steels (FeSi GO 3wt%, M140-27). Their dimensions were 280 x 30 x 0.3 mm, with the length in the easy magnetization direction. Two 120 turns coils were wound around the specimen and plugged in opposite directions as recommended in [START_REF] Moses | AC Barkhausen noise in electrical steels: Influence of sensing technique on interpretation of measurements[END_REF]. The distance separating the sensor coils was set arbitrarily to 10 mm, as used in the previous work [19][53][54] by the authors. The influence of this distance has not been investigated in detail. Two Krohn-Hite 3362 amplifiers-filters were used for the signal conditioning and a National Instruments DAQ USB-6346 acquisition card, controlled through a GUI in Python and of 500 kHz sampling frequency for their acquisitions. Python® and Matlab® were used for the numerical treatment. Eq. ( 63 The 200 mHz current peak is unexpectedly high as compared to other frequencies. However, this difference could be caused simply by the non-linearity of the magnetic material because of the larger difference in amplitude as evident from Table 8. Table 8 summarizes the accuracy of the experimental results by comparing them using sub-section 2.3 criteria.

Tab. 8 -Experimental performances comparison based on the accuracy criteria defined in sub-section 2.

Criteria met

Criteria not met Some Fig. 8 experimental results do not reach the accuracy targeted by the international standards. On the first hand, these standards imposed tight specifications of the experimental conditions (geometry, measurements, etc.), far from the experimental setup depicted in Fig. 7. IEC 60404-3 [50] related to the single sheet tester imposed by instant to use a top and a bottom yoke of large dimensions as displayed in Fig. 9. On the other hand, Fig. 7 setup has been designed based on usual MBN observation methods, like in a non-destructive testing context where magnetization waveform control is not required and never done (see [START_REF] Stewart | Magnetic Barkhausen noise analysis of stress in steel[END_REF]- [START_REF] Ding | Analysis of domain wall dynamics based on skewness of magnetic Barkhausen noise for applied stress determination[END_REF] for examples). This setup delivers magnetic excitation significantly less homogeneous compared to those of [50], and the volume of the tested specimen is reduced. These limitations result in a considerable increase in the magnetization control complexity, justifying a lower accuracy in Table 8 results. A phase delay was also noticed between I(t) and Ba(t). This delay could be a source of divergence for some ILC algorithms (like QNM, whose core assumption relies on the I(Ba) bijectivity). Finally, Fig. 10 shows, as examples, the 20 Hz, 1.8 T sinus flux density, Ba(Hsurf), and MBNenergy(Hsurf) measured cycles. Differences can be observed between the experimental MBNenergy(Hsurf) and the Ba(Hsurf) hysteresis cycles. Those differences were expected and can be explained as follows:  All magnetization contributions are involved in the Ba(Hsurf) hysteresis cycles, mainly the domain wall motions and the magnetization rotation.  For the MBNenergy(Hsurf) cycle, the contribution is limited to the domain wall motions.

It is worth noting the difference at saturation once the cycle is closed. On the first hand, the MBNenergy(Hsurf) cycle reaches a flat saturation. No more variation of the MBNenergy is observed, reflecting the entire disappearance of the domain wall motions. On the other hand, the Ba(Hsurf) still varies. The magnetization rotation remains active and increases the magnetic flux density. The differences are expected to be even more pronounced in the higher frequency range. When Ba(Hsurf) cycles reflect all STL contributions, the MBNenergy(Hsurf) ones are limited to the domain wall motion contributions excluding, notably, the classical loss contribution.

VI. CONCLUSION

Studying the magnetization mechanisms in magnetic cores is a genuine problem that has generated substantial research efforts. A fine study of the MBNenergy(Hsurf) hysteresis cycle excitation frequency dependency and its prediction through an STL-like theory is expected to bring insights into the physical behavior of the magnetization mechanisms.

For this, the flux density has to be imposed sinusoidal from the quasi-static state up to approximately a few hundred Hz, depending on the nature of the tested specimen. Such control might be seen as a simple problem. Still, because of the strong linearities, the practical aspect happens to be especially complex, hence the proliferation of the feedback algorithms in the literature ( [START_REF] Zurek | Use of novel adaptive digital feedback for magnetic measurements under controlled magnetizing conditions[END_REF]- [START_REF] Anderson | Constant dB/dt DC characterisation through digital control of magnetic field[END_REF] are good examples).

In this study, the theoretical problem of Ba control on a classic setup (yoke and sample) has been established, and the ILC (the iterative version of the classic PID controller) has been explained. Different ILC settings exist, and a detailed review of these methods was done in the second section of this manuscript providing valuable insights generalizable to every waveform control environment. For a proper choice, the experimental conditions and the final objective have to be perfectly defined from the very beginning.

Then, six performance criteria have been proposed to identify the most adapted method in the specific context of the MBNenergy(Hsurf) hysteresis cycle characterization, and numerical tests were performed for a comparison purpose followed by conclusions.

P-ILC gave the best performance and the highest satisfaction rate. It was therefore chosen for experimental implementation. Experimental tests were realized on a wide range of amplitude and frequency. We noticed, as expected, a more significant error for higher frequency (limitations of the practical setup bandwidth) and amplitude (stronger non-linear behavior of the tested specimen).

Fig. 1 -

 1 Fig. 1 -Feedback structure.

45 Fig. 4 -

 454 Fig. 4 -Relative Euclidean and form factor differences for the shallow saturation test.

  on matrix inversions like QNM or parameter optimization like LSM-PI are unsuitable. This issue can be partially solved by downsampling the signals, applying correction, and upsampling the resulting waveforms by interpolation (PhC-MGD[START_REF] Stupakov | A system for controllable magnetic measurements of hysteresis and Barkhausen noise[END_REF]). But it means complexity and uncertainty in the measurement treatment. White noise's consequence on the working signals is another issue to consider. Methods like PhC-MGD require an intense averaging process to reach convergence which means extended time acquisition (several cycles) or a sliding window filter. These treatments bring complexity in the signals processing and potentially additional phase delays. Since Ba is obtained by integrating a noisy signal, a drift is always expected. Such a drift can be problematic on methods that normalize signals or expect yM to have a specific amplitude (PhC-MGD, LSM-PI).

Fig. 6 -

 6 Fig. 6 -Relative Euclidean differences, deep saturation test for the different P-ILC methods. Tab. 7 -Comparative results based on the iterations number before convergence. Number of iterations P-ILC 2135 P-ILC-TA (7th degree) 1835 P-ILC-3 979 P-ILC-9 613 P-ILC-TD 55

  Fig.7-Overall 2D view of the Barkhausen noise experimental setup. The power amplifier was a Kepco BOP100-10MG. The excitation coil was made out of 10 turns. The studied specimens were all grain-oriented electrical steels (FeSi GO 3wt%, M140-27). Their dimensions were 280 x 30 x 0.3 mm, with the length in the easy magnetization direction. Two 120 turns coils were wound around the specimen and plugged in opposite directions as recommended in[START_REF] Moses | AC Barkhausen noise in electrical steels: Influence of sensing technique on interpretation of measurements[END_REF]. The distance separating the sensor coils was set arbitrarily to 10 mm, as used in the previous work[START_REF] Fagan | Multiscale modelling of the magnetic Barkhausen noise energy cycles[END_REF][53][START_REF] Fagan | Effect of stress on the magnetic Barkhausen noise energy cycles: A route for stress evaluation in ferromagnetic materials[END_REF] by the authors. The influence of this distance has not been investigated in detail. Two Krohn-Hite 3362 amplifiers-filters were used for the signal conditioning and a National Instruments DAQ USB-6346 acquisition card, controlled through a GUI in Python and of 500 kHz sampling frequency for their acquisitions. Python® and Matlab® were used for the numerical treatment. Eq. (63) -(65) summarized all the tests carried out: B (goal) = B sin(2πf t) (63) B = 1.8 T (64) f = {0.2, 2, 20, 200} Hz (65) P-ILC was used to set the current waveforms. Fig. 8 depicts the experimental results obtained on three decades of frequency and Ba equal to 1.8 T, i.e., the worst-case analyzed scenario in terms of non-linear behavior.

Fig. 8 -

 8 Fig. 8 -a I(t) experimental waveforms, Fig. 8 -b Ba(t) experimental waveforms,

Fig. 9 -

 9 Fig. 9 -IEC 60404-3 single sheet tester illustration and yoke dimensions

Fig. 10 -

 10 Fig. 10 -a Ba(Hsurf) experimental hysteresis cycle, Fig. 10b MBNenergy(Hsurf) renormalized experimental hysteresis cycles.

  40 iterations but needs long calculation times. P-ILC follows with around 60 iterations and minimal calculation times. Finally, FSP-ILC converges after almost 180 iterations. It is worth noting the residual error on P-ILC and FSP-ILC inherent to those methods and impossible to remove. Table3concludes this first set of tests by comparing the methods based on the criteria described in sub-section 2.1:

	Tab. 3 -Performance comparison for the shallow saturation
				test.
		Num. of iter.	Num. of meas.	Fin. rel. euc. diff.	Final form factor	Calc. time per iter. (s)	Total calc. time (s)
	P-ILC	66	66	< 10 -10	8.9•10 -10	4.7•10 -4	0.031
	FDP-ILC	66	66	< 10 -10	8.9•10 -10	4.2•10 -4	0.023
	FSP-ILC	60	60	2.5•10 -8	1.6•10 -9	0.008	0.47
	PhC-P-ILC	174	174	2.2•10 -9	7.6•10 -10	0.01	1.7
	PhC-MGD	11	196	9.9•10 -5	< 10 -10	0.031	0.338
	QNM	40	90	1.4•10 -8	< 10 -10	0.52	20.92
	LSM-PI	7	192	9.0•10 -10	< 10 -10	0.39	2.743

  3.

		200 mHz	2 Hz	20 Hz	200 Hz	Goal (Int. stand.)
	Euc. Dist. (%)	0.678	0.35	0.57	0.34	< 0.5
	Pears. Coef.	0.9999897	0.99999394	0.9999841	0.99999438	> 0.99997
	Form
	Fact. Diff.	0.74	0.29	0.14	0.31	< 0.5
	(dB/dt)(%)
	Amp. Diff. (%)	0.96	0.27	0.05	0.3	< 1
	THD
	(dB/dt)	5.33	1.4	2.95	3.7	< 0.5
	(%)


Thanh Nguyen

Gaspard Tung

Hung Hoang Taxil

Benjamin Nguyen

Mickaël Null Ducharne

Elie Lallart

Hiroki Lefeuvre

Gaël Kuwano

Sebald

Tung

Gaspard Taxil

Hung Hoang Nguyen

Benjamin Ducharne

Mickaël Lallart

Elie Lefeuvre

Hiroki Kuwano

Gael Sebald

; Nguyen

Thanh Tung

Gael Tung

Mickaël Sebald

Benjamin Lefeuvre

Ducharne

Ultimate electromechanical energy conversion performance and energy storage capacity of ferroelectric materials under high excitation levels

Keywords: Ericsson cycle, energy density, energy storage, compressive stress, energy harvesting, ferroelectric

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Ferroelectric materials are considered potential materials for numerous energy harvesting [START_REF] Bowen | Piezoelectric and ferroelectric materials and structures for energy harvesting applications[END_REF][START_REF] Mane | Cyclic energy harvesting from pyroelectric materials[END_REF] and energy / information storage [START_REF] Wu | Lead-free BaTiO 3 -Bi(Zn 2/3 Nb 1/3 )O 3 weakly coupled relaxor ferroelectric materials for energy storage[END_REF][START_REF] Zhao | High-Performance Relaxor Ferroelectric Materials for Energy Storage Applications[END_REF][START_REF] Kumar | High Energy Storage Properties and Electrical Field Stability of Energy Efficiency of (Pb0.89La0.11)(Zr0.70Ti0.30)0.9725O3 Relaxor Ferroelectric Ceramics[END_REF] applications, such as vibrational microgenerators and non-volatile random-access memory. These types of materials are typically pyroelectric and piezoelectric when operating below their Curie temperature. Ferroelectric materials exhibit spontaneous polarization, which can be reoriented or reversed by applying an appropriate and significant electric field in the opposite direction. Based on this property, information could be stored as ferroelectric polarization for non-volatile memory applications [START_REF] Guo | Non-volatile memory based on the ferroelectric photovoltaic effect[END_REF][START_REF] Garcia | Ferroelectric tunnel junctions for information storage and processing[END_REF]. In addition to the paraelectric phase (cubic phase), numerous potential ferroelectric phases exist depending on the shape and polarization of the unit cell, with rhombohedral, orthorhombic, and tetragonal phases being the most common. When a ferroelectric material is in the paraelectric phase above its Curie temperature, it loses its piezoelectric properties.

Energy harvesting consists of converting ambient energy sources into electrical energy, accumulating and storing it for later use [START_REF] Papagiannakis | Energy Harvesting from Roadways[END_REF]. This process has emerged as a promising technique for supplying sensors and sensor networks, limiting the use of natural resources while providing a reliable long-term energy source, especially compared to batteries, which suffer from self-discharge issues. With improvements in nano/microelectronics technology, electronic devices continue to decrease in size, resulting in continuous reduction in power consumption for given operations. Thus, harvested energy may be sufficient to partially or entirely replace batteries in the frameworks of wireless sensors and sensor networks that operate discontinuously. Recently, flexible piezoelectric nanogenerators have shown great potential for the next generation of wearable electronic-based devices for human health and movement monitoring [START_REF] Wu | Self-powered skin electronics for energy harvesting and healthcare monitoring[END_REF][START_REF] Zhao | Flexible nanogenerators for wearable electronic applications based on piezoelectric materials[END_REF]. Energy harvesting can generate electricity from several sources, including kinetic, thermal, solar, or electromagnetic energy [START_REF] Prauzek | Energy Harvesting Sources, Storage Devices and System Topologies for Environmental Wireless Sensor Networks: A Review[END_REF][START_REF] Kim | Application of ferroelectric materials for improving output power of energy harvesters[END_REF][START_REF] Liu | A novel configuration for high power-output and highly efficient vibration energy harvesting[END_REF][START_REF] Na | Wind energy harvesting from a magnetically coupled piezoelectric bimorph cantilever array based on a dynamic magneto-piezo-elastic structure[END_REF][START_REF] Vocca | Kinetic energy harvesting with bistable oscillators[END_REF]. Human activities such as walking or opening a door are also considered an attractive energy source for mechanical energy harvesting using dielectric or electromagnetic transduction methods [START_REF] Romero | Energy scavenging sources for biomedical sensors[END_REF][START_REF] Defosseux | Highly efficient piezoelectric micro harvester for low level of acceleration fabricated with a CMOS compatible process[END_REF]. For piezoelectric energy harvesting, most of the proposed devices work as resonators harvesting energy from vibrations. The energy density per cycle can range from 10 µJ/cm3 to 25 mJ/cm 3 [START_REF] Elfrink | Vibration energy harvesting with aluminum nitride-based piezoelectric devices[END_REF] of the active material at most. However, such a value is solely based on empirical analysis of achieved experimental devices. Until now, no investigation has been conducted on the theoretical limit in terms of convertible and harvestable energy, although such a study would: 1) allow a comparison of the performance of energy harvesters with respect to the ultimate energy that can actually be harvested (i.e., in a similar way to the Carnot efficiency for thermal devices) and 2) provide insights into the design of optimized harvesters to work close to their theoretical limits in terms of output energy. In contrast, research efforts have been made to investigate the pyroelectric energy conversion capabilities of ferroelectric materials through Olsen or Ericsson cycles, with output energies far above those given by the linear pyroelectric coefficient [START_REF] Starner | Human-powered wearable computing[END_REF][START_REF] Vats | Enormous energy harvesting and storage potential in multiferroic epitaxial thin film hetrostructures: an unforeseen era[END_REF][START_REF] Sebald | Energy harvesting based on Ericsson pyroelectric cycles in a relaxor ferroelectric ceramic[END_REF][START_REF] Unruan | Estimation of energy density of PMN-PT ceramics utilizing mechanical stress[END_REF][START_REF] Patel | Enhanced energy harvesting in commercial ferroelectric materials[END_REF][START_REF] Bakshi | Effect of Stress on Energy Conversion and Storage Characteristics of (1-x-y)PIN-xPMN-yPT Single Crystals[END_REF][START_REF] Smith | Cascaded pyroelectric conversion: Optimizing the ferroelectric phase transition and electrical losses[END_REF]. However, the maximum amount of energy that can be harvested through electromechanical coupling using ferroelectric materials is still unknown; however, it is of primary importance [START_REF] Chen | A high density piezoelectric energy harvesting device from highway traffic -System design and road test[END_REF].

Ericsson cycles include two iso-fields and two iso-stresses loadings, which were used in this study to estimate the energy conversion performance of ferroelectric materials. Previous studies have demonstrated that coupling an electric field and mechanical excitation through Ericsson cycles yielded a significant increase in harvested energy [START_REF] Zhang | Energy harvesting based on piezoelectric ericsson cycles in a piezoceramic material[END_REF][START_REF] Zhang | Energy harvesting based on piezoelectric Ericsson cycles in a piezoceramic material[END_REF]. This study investigated the energy conversion performance of various ferroelectric ceramics and single crystals that underwent ferroelectric and paraelectric phases through Ericsson cycles. We also developed a simple and accurate semi-empirical estimation of energy density based on bipolar hysteresis curves (P-E loops) under compressive stresses of 1 and 100 MPa. The ferroelectric phase of lead zirconate titanate (PZT) ceramics (hard PZT C203, medium PZT C6, soft PZT C9), single crystal lead magnesium niobate-lead titanate (PMN-25PT), single crystal lead zinc niobate-lead titanate (PZN-8PT), and the paraelectric phase of ceramic PMN 15 were characterized. The objective of this study was to estimate the ultimate energy conversion capability of ferroelectric materials under intermediate and high levels of stresses and electric fields.

Finally, as a complementary investigation, these materials were also assessed for electrical energy storage applications.

Experimental methods

Materials and testbench

Disk-shaped electrostrictive PMN 15 ceramics (10 mm in diameter and 1 mm thick, with silver electrodes) were purchased from TRS Technologies, State College, PA, USA. PMN 15 exhibited a maximal dielectric constant at 15 o C and was in the paraelectric phase at room temperature. Hard, medium, and soft piezoelectric lead zirconate titanate (Pb(ZnTi)O3) (PZT) ceramics C203, C6, and C9 were purchased from Fuji Ceramics Corporation, Fujinomiya, Japan. In addition, single crystal [001]oriented lead-magnesium niobate-lead titanate PMN-25PT was purchased from Otsuka Corporation, Japan, and [001]-poled single crystal lead zinc niobate-lead titanate PZN-8PT was purchased from Microfine Materials Technologies, Singapore. All samples had the same thickness of 1 mm. Table 1 shows the properties of these samples. PZT C203, C6, C9, and PMN-25PT with a section of 6 x 6 mm² and PZN-8PT (2.5 x 2.5 mm²) were cut by a dicing machine (DISCO, DAD3240). Photographs of the samples are shown in Figure 1c. Then, the cut samples were poled in a silicone oil bath using the field cooling technique [START_REF] Kounga | High-temperature poling of ferroelectrics[END_REF]. The treated samples were heated to the desired temperature (200, 180, 130, 90, and 150 o C for PZT C203, C6, C9, PMN-25PT, and PZN-8PT, respectively), where an electric field of 3 kV/mm was applied for 30 min and maintained while cooling down. When the sample temperature reached room temperature, the electric field was powered off, and the sample was removed. The Ericsson cycles and bipolar cycle measurements could also act as a poling process, which would be sufficient for most of the tested materials. However, for hard PZT, it was necessary to pole it beforehand, as the material could not be poled at room temperature. For the other materials, if the electric field was not very high compared to the coercive field, we could have partially maintained a depolarized state during the cycle, which would be detrimental to the piezoelectric activity. As a result, we chose to prepare the fully poled samples as a reference status before applying the electric field and mechanical stress Figure 1a shows a schematic diagram of the measurement system, where a Tektronix® AFG1022 arbitrary function generator and Trek® 10/10B-HS high voltage amplifier were used to apply the electric field. A Shimadzu® AGS-X series traction machine was used to apply compressive stress to the sample, as shown in Figure 1b. The current through the sample was monitored using a low-noise current amplifier (Stanford research®, Sunnyvale, CA, United States of America), and a Digilent® data acquisition system (Pullman, WA United States of America) was used to monitor and record the applied voltage and electrical current. An overload protection circuit was used to protect the measuring system, and the sample under test was immersed in silicone oil to avoid electric arcing during measurement. In processes 1-2, the electric field was increased from 0 to 1.5 kV/mm under constant compressive stress of 1 MPa. In processes 2-3, the compressive stress was increased from 1 to 100 MPa under a constant 1.5 kV/mm electric field. In processes 3-4, the electric field was decreased from 1.5 to 0 kV/mm under a constant 100 MPa compressive stress. In processes 4-1, compressive stress was released from 100 MPa to 1 MPa under a constant 0 kV/mm electric field.

Materials

Ericsson cycle application and analysis methods

Ericsson cycles are used to convert thermal or mechanical energy into electrical energy with ferroelectric materials [START_REF] Unruan | Estimation of energy density of PMN-PT ceramics utilizing mechanical stress[END_REF][START_REF] Patel | Enhanced energy harvesting in commercial ferroelectric materials[END_REF][START_REF] Bakshi | Effect of Stress on Energy Conversion and Storage Characteristics of (1-x-y)PIN-xPMN-yPT Single Crystals[END_REF][START_REF] Smith | Cascaded pyroelectric conversion: Optimizing the ferroelectric phase transition and electrical losses[END_REF]. Figure 2a shows the electric field and compressive stress waveforms during the Ericsson cycle, where compressive stresses of 1 and 100 MPa were applied to the PMN-25PT sample. A low-stress value of 1 MPa was chosen to ensure good electrical contact, and the applied electric field was swept from 0 to 1.5 kV/mm for all samples. Polarization was calculated by integrating the current in time. Figure 2b shows a PMN-25PT Ericsson cycle as a representative example. Between steps 1 and 2, the electric field increased from 0 to a high electric field (1.5 kV/mm) for 1 second under a constant low compressive stress level of 1 MPa.

The polarization was increased following the electric field variations. Electrical energy was transferred from the voltage amplifier to the sample, whose polarization changed from P1 to P2 (equivalent to capacitor charging). From steps 2 to 3, the compressive stress increased from 1 to 100 MPa for 6 s under a constant electric field of 1.5 kV/mm. The applied compressive stress induced the reduction of polarization from P2 to P3, which could be used to harvest electrical energy [START_REF] Dw | Compressive stress-induced depolarization in ferroelectric ceramics at room and high temperatures: experiment and prediction[END_REF][START_REF] Dw | Nonlinear behavior of ferroelectric ceramics during mechanical depolarization at room and high temperatures: experiment and prediction by an experimental formula[END_REF][START_REF] Yimnirun | Stress-dependent ferroelectric properties of Pb(Zr1/2Ti1/2)O3-Pb(Zn1/3Nb2/3)O3 ceramic systems[END_REF][START_REF] Unruan | Effects of Parallel and Perpendicular Compressive Stresses on the Dielectric and Ferroelectric Properties of Soft PZT Ceramics[END_REF]. This polarization drop originated from the switching phenomenon of the ferroelectric domains during mechanical depolarization and the phase transitions induced by mechanical stress. From step 3 to step 4, the applied electric field decreased from 1.5 kV/mm to 0 for 1 second under constant compressive stress of 100 MPa. The polarization decreased following the reorientation of the ferroelectric domains, resulting in a small amount of discharging current. From steps 4 to 1, the compressive stress was released to the original value of 1 MPa to complete the Ericsson cycle. All Ericsson cycles were performed at room temperature (20 o C). The area of the Ericsson cycle (1-2-3-4) was equivalent to the converted mechanical energy density of the sample under test. Similar to the Olsen cycle [START_REF] Nguyen | Pyroelectric energy converter using co-polymer P(VDF-TrFE) and Olsen cycle for waste heat energy harvesting[END_REF][START_REF] Siao | Pyroelectric waste heat energy harvesting using the Olsen cycle on Pb(Zr, Ti)O 3 -Pb(Ni, Nb)O 3 ceramics[END_REF], this area could be calculated by Eq. 1 [START_REF] Unruan | Estimation of energy density of PMN-PT ceramics utilizing mechanical stress[END_REF][START_REF] Patel | Enhanced energy harvesting in commercial ferroelectric materials[END_REF][START_REF] Bakshi | Effect of Stress on Energy Conversion and Storage Characteristics of (1-x-y)PIN-xPMN-yPT Single Crystals[END_REF][START_REF] Smith | Cascaded pyroelectric conversion: Optimizing the ferroelectric phase transition and electrical losses[END_REF]:

𝑊 = ∮ 𝐸•𝑑𝑃 𝐴 = ∮ 𝑃•𝑑𝐸 𝐴 , (Eq. 1)
where P, E, and A denote the polarization, electric field, and sample area, respectively. All tests consisted of three consecutive cycles to avoid polarization irreversibility issues, and only data from the last cycle was monitored.

Figure 3: Bipolar P-E loops under lower and higher compressive stresses, with the estimated Ericsson cycle P1-Pr-C-D-P1 using descending-descending curves of bipolar P-E loops (Estimation 1). Our semi-empirical estimated Ericsson cycle P1-g(E)-C-D-P1

was based on the first-order reversal curves.

Bipolar P-E loops-based estimated Ericsson cycles

In Uruan et al. and Patel et al. [START_REF] Unruan | Estimation of energy density of PMN-PT ceramics utilizing mechanical stress[END_REF][START_REF] Patel | Enhanced energy harvesting in commercial ferroelectric materials[END_REF], the Ericsson cycle energy densities were estimated using the two descending curves of the bipolar P-E loops under low and high stresses, respectively (P1-Pr-C-D-P1 cycle in Figure 3). However, this technique was imprecise due to the irreversibility of the polarization decrease under such high stress (hysteresis issue). In this work, an alternative semiempirical estimation technique was considered, in which the ascending curve of the P-E loop under low-stress f(E) shifted upward to become g(E), which was equivalent to a first-order reversal curve [START_REF] Goh | First-order reversal curve diagrams for characterizing ferroelectricity of Hf 0.5 Zr 0.5 O 2 films grown at different rates[END_REF]. Hence, the estimated Ericsson cycle became the P1-g(E)-C-D-P1 cycle, where g(E) was calculated from f(E) according to the following relation:

𝑔(𝐸) = 𝑓(𝐸) • 𝑃 2 -𝑃 1 𝑃 2 -𝑃 3 + 𝑃 1 -𝑃 3 • 𝑃 2 -𝑃 1 𝑃 2 -𝑃 3 , (Eq. 2)
The bipolar P-E loops of all samples were obtained by applying a sinusoidal voltage of 1 Hz and an electric field of ±1.5 kV/mm.

Energy storage

In addition to the energy conversion assessment, the potential applications of these ferroelectric and paraelectric materials as electrostatic capacitors were considered. The energy storage capacity could be estimated by the surface integral of the bipolar P-E loops [START_REF] Bakshi | Effect of Stress on Energy Conversion and Storage Characteristics of (1-x-y)PIN-xPMN-yPT Single Crystals[END_REF][START_REF] Chauhan | Mechanical confinement for improved energy storage density in BNT-BT-KNN lead-free ceramic capacitors[END_REF][START_REF] Patel | Enhancing electrical energy storage density in antiferroelectric ceramics using ferroelastic domain switching[END_REF]:

𝑈 = ∮ 𝐸•𝑑𝑃 𝐴 𝑃 𝑚𝑎𝑥 𝑃 𝑟 , (Eq. 3)
where Pr and Pmax denote the remnant and maximum polarization values, respectively.

Results

Energy conversion performance was first investigated (from 1 to 100 MPa compressive stress, electric fields of 0 to 1.5 kV/mm) using Ericsson cycles for PZT C203, C6, C9, PMN-25PT, PZN-8PT, and PMN 15 under significant compressive stress variations (Figure 4). This was followed by intermediate compressive stress variations (from 1 to 30 MPa compressive stress, electric fields of 0 and 1.5 kV/mm) (Figure 5). Comparisons of the converted energy densities for the high and intermediate This polarization decrease was the smallest for the paraelectric phase PMN 15, smaller for hard PZT C203, larger for medium PZT C6, much larger for soft PZT C9, and the greatest for the PMN-25PT, PZN-8PT single crystal. The reduction in polarization at an electric field of 1.5 kV/mm, by increasing the compressive stress from 1 to 100 MPa, could be roughly estimated by an equivalent effective piezoelectric coefficient d33_eff (at 1.5 kV/mm):

∆𝑃 = 𝑑 33_𝑒𝑓𝑓 • ∆𝑇 3 , (Eq. 4)
where ΔT3 is the applied stress difference in the 3-direction (z-direction). This value was smaller for hard PZT, larger for soft PZT, and the largest for the single crystal ferroelectric materials, resulting in different influences of compressive stresses on polarization reduction. The ferroelectric materials in the paraelectric phase lost their piezoelectric properties; thus, the energy density of the PMN 15 was the smallest. The same tendency was observed at the intermediate level (Figure 6b). The energy densities of PMN-25PT and PZN-8PT single crystals were much higher than the other materials. For high-level variations, the energy density of the PZT ceramic C9 was comparable to the PMN-25PT single crystal. In the framework of a realistic energy harvester, considering the cost and complexity of fabrication process, PZT C9 (or similar soft PZT materials) appeared to be a plausible candidate for future applications. To finalize the analysis, much lower stresses and electric field levels were also considered. When limiting the electric field to 100 V/mm and the mechanical stress to 1 MPa, we could assume that the material would remain in its linear region, both in terms of the electric and mechanical responses.

where d33 (C/N) denotes the piezoelectric coefficient, T (Pa) the stress, and E (V/m) the electric field. The piezoelectric coefficients of PZT C203, C6, C9, PMN-25PT, PZN-8PT, and PMN 15 are shown in Table 1. The piezoelectric coefficient of PMN 15 was deduced to be around 50 (pC/N), and the energy densities in this case are given in Figure 6c. The energy densities of these samples at low levels manifested the same tendency as those observed at high and intermediate levels.

Predicting these values from the bipolar P-E loops may also be of significant interest, to quickly assess the energy conversion and harvesting performance. Hence, an Ericsson cycle energy density estimation was conducted for PZT C6 from the descending curves [START_REF] Unruan | Estimation of energy density of PMN-PT ceramics utilizing mechanical stress[END_REF][START_REF] Patel | Enhanced energy harvesting in commercial ferroelectric materials[END_REF] and compared to the experimental Ericsson cycles. were inaccurate and much larger than those obtained by the actual Ericsson cycles. This difference was attributed to the hysteretic nature of the sample response. An alternative semi-empirical estimation technique based on the first-order reversal curve was proposed to address this issue. Figure 8 shows the estimated Ericsson cycles for PZT C6, C9, PMN-25PT, PZN-8PT, and PMN 15, under high variations of 1 to 100 MPa compressive stresses, using this approach. Of note, we expect that testing the Ericsson cycles at a much higher electric field would make both estimations similar to the actual Ericsson cycle. More precisely, if the electric field amplitude was For PMN-25PT and PMN 15, the DC loss effect was stronger, resulting in estimated values that were slightly higher than the experimental values. At an intermediate level, the effect of DC losses was significant compared to the effect of 30 MPa of compressive stress, resulting in energy densities obtained from real Ericsson cycles that were smaller than those obtained from the estimated densities.

The energy conversion performance of hard PZT C203 was not estimated under current experimental conditions (1.5 kV/mm electric field), because the bipolar hysteresis curves of hard PZT C203 were unsymmetrical; therefore, this would require a larger electric field or higher stress, which was beyond the scope of this study. Following this analysis, we further studied the potential applications of these ferroelectric materials in the field of electrical energy storage, and the stored energy in the samples was estimated, as shown in Figure 10, while the combined data are shown in Figure 10f. The results indicated that the static mechanical compressive stress could enhance the stored energy for all samples, except PZN-8PT at 30 MPa. In the case of 100 MPa, the gain in energy storage ranged from +20% to +100%, depending on the materials, and the mechanical stress tended to decrease the remnant polarization. If the saturation polarization was still of a high value, the energy density increased as it was roughly proportional to the difference between the saturation polarization and the remnant one. The highest gain was observed for PZT C9 (+100%), with its behavior under 100 MPa favoring a significant difference between remnant and saturation polarization. Unlike electromechanical energy conversion, where the paraelectric phase PMN 15 exhibited the lowest performance, the paraelectric material showed the highest level of energy storage. Because the stored energy was estimated by the surface integral of the bipolar P-E loops with polarization values from remnant to saturation polarization, the paraelectric phase PMN 15 had the smallest remnant polarization compared to the other materials; therefore, the stored energy of PMN 15 was the largest.

Discussion

Energy conversion ability

For applications of piezoelectric materials featuring low excitation level, such as sensors, energy harvesting, or to some extent power transducers, performance can usually be assessed by considering the materials' linear properties. For example, in resonant energy harvesting systems, the electromechanical coupling factor k 2 of the system along with the mechanical quality factor Qm of the resonator will play important roles [START_REF] Lefeuvre | Materials, structures and power interfaces for efficient piezoelectric energy harvesting[END_REF][START_REF] Guyomar | In press. Energy Harvesting using Nonlinear Techniques[END_REF]. However, in other analyses, the products of piezoelectric constants d and g were more relevant, depending on the mechanical excitation of the material [START_REF] Deutz | Analysis and experimental validation of the figure of merit for piezoelectric energy harvesters[END_REF]. In Ericsson cycles used for electromechanical energy conversion, with a fully controlled electric field and mechanical stress, the converted energy density, i.e., the area enclosed by points 1 to 4 in Figure 2b, was only dependent on the effective d33_eff coefficient, as given in Eq. 5. Of note, the ranking of the materials followed their d33 values for low levels and even for intermediate levels. However, at the higher electric field and mechanical stress levels, the output energies converged due to polarization saturation, which was mainly dependent on the dipole moment and dipole density. For example, for low-level energy harvesting, the energy ratio of PZT C9 and PMN-25PT was given by the ratio of d33, which was approximately 2.2 times larger for PMN-25PT. For energy cycles of 100 MPa and 1.5

kV/mm, this ratio dropped to roughly 1.1.

PZN-8PT requires additional analysis, as it did not follow this trend. This material exhibited multiple phase transitions under the conditions of mechanical stress, electric field, and temperature [START_REF] Davis | Electric-field-, temperature-, and stress-induced phase transitions in relaxor ferroelectric single crystals[END_REF][START_REF] Janolin | Uniaxial-stress induced phase transitions in [001]C-poled 0.955Pb(Zn1∕3Nb2∕3)O3-0.045PbTiO3[END_REF].

From the experimental results obtained in this study, it appeared that the effective d33_eff value remained very high, even under 1.5 kV/mm. Hence, the polarization decrease was significant due to the mechanical stress. This likely occurred because of the phase transition induced by the stress, although further investigations will be necessary to verify this assumption.

The overall convergence of all materials was further analyzed by considering the weak scattering of the spontaneous polarization among the ferroelectric materials, ranging from 0.15 to 0.4 C/m 2 for typical materials [START_REF] Damjanovic | Ferroelectric, dielectric and piezoelectric properties of ferroelectric thin films and ceramics[END_REF][START_REF] Wongdamnern | Comparative studies of dynamic hysteresis responses in hard and soft PZT ceramics[END_REF]. If the electric field was high enough to entirely saturate the material, and if the mechanical stress was sufficient to completely depolarize the material, then the converted energy would be similar for all tested materials. In this work, both the electric field and mechanical stress were far from these ideal conditions, resulting in a ranking between the materials by a factor of four, in terms of the converted energy from the lowest to the highest, especially in the vicinity of the phase transitions induced by both the electric field and mechanical stress.

Considering an electric field value of 1.5 kV/mm and saturation polarization of 0.4 C/m², the ultimate energy density would be 600 mJ/cm 3 . The experimental results obtained in this study ranged from 50 mJ/cm 3 to 236 mJ/cm 3 , with the latter representing 40% of the ultimate energy density.

Furthermore, the output energy was dependent on the chosen electric field of the Ericsson cycle.

Applying an electric field that was two times larger would approximately double the output energy.

By contrast, depolarization induced by the mechanical stress would not be proportional to the stress level, and would quickly saturate when most of the polarization was already lost. Regarding the limit of the electric field, it was directly related to the dielectric strengths of the materials. In the case of bulk ferroelectric materials, applications working close to the limit can be related to power transducers and electrocaloric materials [START_REF] Valant | Electrocaloric temperature change constrained by the dielectric strength[END_REF]. The dielectric strength can be related to the probability of finding defects in ceramics and single crystals and will be lower for bulk materials than thick and thin films.

For bulk ceramics, the values range from 5 kV/mm to 15 kV/mm [START_REF] Takeda | Mechanism study of the residual stress evaluation of low-carbon steels using the eddy current magnetic signature method[END_REF][START_REF] Gerson | Dielectric Breakdown of Porous Ceramics[END_REF][START_REF] Yang | Dielectric Breakdown Properties of Zr-Rich Lead Zirconate Titanate Ceramics[END_REF], while higher values have been successfully tested for single crystals [START_REF] Park | Ultrahigh strain and piezoelectric behavior in relaxor based ferroelectric single crystals[END_REF]. In the case of thick films (typically those deposited to form multilayer capacitors), to maintain a safe level of the electric field, piezoelectric layers will usually be limited to 3 kV/mm [START_REF] Esteves | Electric-field-induced structural changes in multilayer piezoelectric actuators during electrical and mechanical loading[END_REF], although some attempts on multi-layered capacitors have allowed tests up to 15 kV/mm [START_REF] Nair | Large electrocaloric effects in oxide multilayer capacitors over a wide temperature range[END_REF]. Considering the previous limits, a limit value of 1.5 kV/mm could possibly be increased up to 3 to 10 times higher for some compounds, although the fatigue life might significantly decrease [START_REF] He | A Probabilistic Approach to Evaluate the Reliability of Piezoelectric Micro-Actuators[END_REF]. Regarding the mechanical stress limit, a combined application of an electric field and stress could possibly lead to premature fatigue with stress as low as 100 MPa [START_REF] Makino | Effects of dc Electric Field on Mechanical Properties of Piezoelectric Ceramics[END_REF].

Consequently, in this study, the electric field was limited to 1.5 kV/mm and the stress was 100 MPa.

However, further detailed analysis will be necessary to explore the upper limits of the different materials. Nonetheless, the orders of magnitude should not be significantly modified as a margin of improvement, for both limits might remain weak for bulk materials when considering high cycle capability. Furthermore, the selected limits are in line with targeted application of left-behind devices operating without any maintenance.

Interestingly, the order of magnitude in terms of energy density for all materials was in the same range as pyroelectric energy harvesting under significant temperature variations, reaching a few hundred mJ/cm 3 /cycle [START_REF] Bowen | Pyroelectric materials and devices for energy harvesting applications[END_REF]. Generally speaking, this convergence was attributed to the trade-off between the input energy (higher in thermal input and lower for mechanical excitation) and the conversion abilities (higher electromechanical coupling coefficient compared to electrothermal).

However, the use of Ericsson cycles on piezoelectric materials would not provide particular advantages in all applications of mechanical to electrical energy conversion. In the case of energy harvesting from vibration using resonant systems, provided that the electromechanical coupling k² and mechanical quality factor Qm were sufficiently high (k²Qm>>1), energy harvesting would result in a strong damping effect. In such a case, the standard energy harvesting processes could convert as much energy as possible with known strategies. By contrast, in the case of weakly coupled problems, or for a non-resonant mechanical source, the best strategy would be to maximize the energy conversion per cycle. A similar situation was already encountered in the case of pyroelectric energy harvesting. Due to low electrothermal coupling, the energy conversion strategy had little impact on the input thermal energy; thus, the Ericsson cyclesalso known as Olsen cycles in several worksbrought much larger energy conversion density compared to the use of the solely linear pyroelectric coefficient [60][61].

Not only could Ericsson cycles be employed, but several thermodynamic cycles could also be envisaged, such as Stirling or Brayton cycles, each providing specific advantages about its feasibility and energy density [START_REF] Pandya | New approach to waste-heat energy harvesting: pyroelectric energy conversion[END_REF]. Looking back to the electromechanical energy conversion, it is possible that a case study could consist of energy conversion from a single mechanical source displacement at a given force, like those encountered on roads [START_REF] Chen | A high density piezoelectric energy harvesting device from highway traffic -System design and road test[END_REF][63], pathways, or sea waves [START_REF] Zhang | Experimental sea wave energy extractor based on piezoelectric Ericsson cycles[END_REF].

Practical implementation of Ericsson cycles

Even if we exposed and discussed the energy conversion abilities of the ferroelectric or dielectric materials, the actual implementation of the associated electrical interface was outside the scope of this study. Some preliminary considerations could still be drawn regarding the practical realization of these cycles. More precisely, practical implementation of long-term energy harvesting strategies cannot include external battery unless the scheme encompasses recharging.

In this view, a potential solution would involve circuits with a bidirectional energy flow, including controlled switches for electrically charging and discharging the active element (transducer).

According to the previously exposed constraints regarding energy self-consistency, a pre-charged battery or capacitor, with a value much larger than the maximum transducer capacitance, could also harvest energy when discharging the transducer. For example, a possible implementation would consist of DC/DC switched converters (Figure 11). Specifically, the application of an electric field in the transducer could be performed using a buck converter that was activated as the stress was applied, while the energy recovery phase (electric field release) could be achieved using a boost converter.

Additionally, a diode could be placed in an anti-parallel configuration with the transducer, to prevent inverse conduction. As a first approximation, the transducer was modeled as a variable capacitor whose capacitance was related to the mechanical excitation, and Figure 12 depicts the obtained simulation implementation and results. Thus, the obtained cycles were consistent with those experimentally achieved in previous investigations. However, a storage capacitor much larger than the maximum transducer capacitance was required to maintain a constant voltage during the charging phase. Otherwise, the applied voltage would change during the process.

Furthermore, according to the energy cycle (Figure 12c), the clamped voltage increased with time.

This was explained by the harvesting process, where the additional energy that was transferred to the storage stage, made of the capacitor, actually increased the voltage of the capacitor. This type of voltage increase was eventually limited by the losses within the circuit or simply by the energy consumed by the load (e.g., sensor), and this type of a voltage-constrained approach could also be implemented without an inductor, but at the cost of more complex switching topology and control [START_REF] Torres | A 0.7-µm BiCMOS Electrostatic Energy-Harvesting System IC[END_REF].

However, the previously exposed scheme required active switches, for instance those composed of transistors, driven by signals that were synchronized to mechanical excitation. Although the generation of these control signals and their synchronization required minimal energy, their implementation could potentially yield some complexity that may compromise the energy balance, as the harvested energy may not be sufficient. Another approach that did not directly involve Ericsson cycles consisted of self-synchronous converters as the Bennett doubler, where an overlap between mechanical excitation and electrical changes appeared [START_REF] Ghaffarinejad | Bennet's doubler working as a power booster for triboelectric nano-generators[END_REF]. This type of approach showed a significant advantage, i.e., once pre-charged, it provided a virtuous energy cycle loop in a purely passive way, without requiring any switch control; therefore, it provided significant energy savings. A simulation example is shown in Figure 13, which was still based on a variable capacitor. Hence, while not precisely an Ericsson cycle due to non-zero voltage discharge and overlap between mechanical excitation and electric field, the obtained results still shared strong similarities with the previously presented principles. Hence, we inferred that the relative material performance would be similar. Finally, we noted that such a self-synchronous structure could be enhanced by including multi-stage feedback loops that can handle limited capacitance variations [START_REF] Lefeuvre | Self-Biased Inductor-less Interface Circuit for Electret-Free Electrostatic Energy Harvesters[END_REF].

Conclusions

In this work, the electromechanical energy conversion performance of various ferroelectric (hard PZT C203 ceramic, medium PZT C6 ceramic, soft PZT C9 ceramic, and PMN-25PT and PZN-8PT single crystal) and paraelectric (PMN 15 ceramic) materials was investigated using Ericsson cycles under intermediate and high compressive stress levels (30 MPa and 100 MPa) with an electric field of 1.5

kV/mm and at room temperature (20 o C). The ferroelectric PZN-8PT single crystal exhibited the highest potential for energy conversion among the assessed materials, with an energy density reaching 235.6 mJ/cm 3 . However, under high electric field and stress levels, the energy conversion performance of the PZT C9 ceramic was found to be comparable to the PMN-25PT single crystal, with an energy density of 118.7 mJ/cm 3 . PZT C9 was found to be an excellent candidate for practical applications due to its well-established commercial fabrication routes, yielding an affordable transducer. Measurement of the Ericsson cycle remained a time-consuming process, as it required complete control of the electric field and stress with time, along with their synchronization. Thus, a simple and accurate semiempirical estimation technique was proposed to evaluate the energy conversion performance of the ferroelectric materials via the compressive stress-induced bipolar hysteresis curves. This was found to be much more precise, compared to solely relying on the descending hysteresis P-E curves.

Estimations using this approach were accurate for PZT C6, C9, and PZN-8PT and slightly larger for PMN 15 and PMN-25PT compared to the experimental results. The acquisition of bipolar hysteresis curves was easier than the Ericsson cycles; thus, this method permitted the fast screening of ferroelectric materials. The results obtained in this study were dependent on the electric field amplitude used during the Ericsson cycle measurements. Consequently, a higher dielectric strength for a given material might allow for a larger electric field and output energy. Nevertheless, utilizing a safe operating condition (which is necessary for ensuring reliable operation of the device in the long term) in the kV/mm range for bulk ferroelectrics led to energy in the hundreds of mJ/cm 3 . Increasing the electric field amplitude by a factor of 10 (highest dielectric strength on bulk ferroelectrics) would lead to an energy density in the J/cm 3 range at the highest.

Lastly, we further studied the potentials of these materials for electrical energy storage applications.

The maximum obtained stored energy was found to be 104.7 mJ/cm 3 for the paraelectric PMN 15 ceramic at a high level of 100 MPa compressive stress. In contrast to electromechanical energy conversion, in which paraelectric PMN 15 exhibited the lowest performance compared to the other ferroelectric materials, in terms of energy storage, PMN 15 had the highest amount of stored energy.

This study paves the way for future applications of these materials for energy harvesting and energy storage.

Figure 1 :

 1 Figure 1: (a) Schematic diagram of the measurement system to obtain the Ericsson cycles and bipolar P-E loops, (b) photograph of the sample under compressive stress using a Shimadzu electromechanical test frame, and (c) photographs of the tested samples.

Figure 2 :

 2 Figure 2: (a) Waveforms of the applied electric field and compressive stress for the Ericsson cycles (1 MPa and 100 MPa, which were applied on PMN-25PT for illustration); (b) Ericsson cycle of PMN-25PT for illustration.In processes 1-2, the electric field was increased from 0 to 1.5 kV/mm under constant compressive stress of 1 MPa. In processes 2-3, the compressive stress was increased from 1 to 100 MPa under a constant 1.5 kV/mm electric field. In processes 3-4, the electric field was decreased from 1.5 to 0 kV/mm under a constant 100 MPa compressive stress. In processes 4-1, compressive stress was released from 100 MPa to 1 MPa under a constant 0 kV/mm electric field.

  mechanical variations are shown in Figures6a and 6b. For the high variations (Figure6a), and considering the ferroelectric phase ceramics, the energy density was smaller for hard PZT C203, larger for medium PZT C6, and much larger for soft PZT C9. The ferroelectric phase single crystals (PMN-25PT, PZN-8PT) exhibited excellent energy density, especially PZN-8PT, which had tremendous values. Regarding the paraelectric phase material, the energy density of PMN 15 was smaller than the ferroelectric phase materials. The application of mechanical compressive stresses under a constant electric field of 1.5kV/mm induced polarization decreases as low as ∆𝑃 = 𝑃 2 -𝑃 3 of 49.7, 66.0, 70.2, 83.5, 208, and 36.6 mC/m 2 for ferroelectric phase ceramics of hard PZT C203, medium PZT C6, soft PZT C9, single crystal PMN-25PT and PZN-8PT, and paraelectric phase PMN 15, respectively.

Figure 4 :

 4 Figure 4: The Ericsson cycles of PZT C203, C6, C9, PMN-25PT, PZN-8PT, and PMN 15 at high levels of compressive stress of 1 and 100 MPa, with 0 and 1.5 kV/mm electric fields.

Figure 5 :

 5 Figure 5: The Ericsson cycles of PZT C203, C6, C9, PMN-25PT, PZN-8PT, and PMN 15 at intermediate level of compressive stress of 1 and 30 MPa, and 0 and 1.5 kV/mm electric fields.

Figure 6 :

 6 Figure 6: Energy densities of samples (a) at high (1 and 100 MPa, 0 and 1.5 kV/mm), (b) intermediate (1 and 30 MPa, 0 and 1.5 kV/mm), and (c) low (1 MPa and 100 V/mm) levels.

Figure

  Figure 7a and 7b show the descending curves based on the estimated Ericsson cycles at intermediate level (under compressive stresses of 1 and 30 MPa) and high level (under compressive stresses of 1 and 100 MPa), respectively. The energy densities were calculated and compared to those obtained from the actual experimental Ericsson cycles, as shown in Figure 7c. For the tested samples, the energy densities obtained by the descending curves based on the estimated Ericsson cycles (Estimation 1)

Figure 7 :

 7 Figure 7: Descending-descending curves based on the estimated Ericsson cycles for PZT C6 under compressive stresses of (a) 1 and 30 MPa and (b) 1 and 100 MPa, and (c) comparison of energy density estimation using the actual Ericsson cycles and descending-descending curves based on the estimated Ericsson cycles (Estimation 1) of PZT C6. The energy densities were calculated and compared to those obtained by the experimental Ericsson cycles, as shown in Figure 8f. These results indicated that the energy densities of the samples obtained using this method were much closer to those obtained by the real measured cycles for PZT C6 and C9, and slightly larger for PMN-25PT, PZN-8PT, and PMN15. Figure 9 depicts the estimated Ericsson cycles for PZT C6, C9, PMN-25PT, PZN-8PT, and PMN 15 under intermediate variations of 1 to 30MPa of compressive stress. The energy densities were calculated and compared to those obtained from the actual experimental Ericsson cycles at an intermediate level, as shown in Figure9f. The results indicated that the energy densities obtained by the estimated Ericsson cycles were 2 to 3 times larger than those obtained experimentally, and the estimation appeared to be accurate for the 100 MPa case, but somewhat inaccurate for the 30 MPa case. The reason for the discrepancy at 30 MPa was due to two opposite effects on the polarization change in processes 2-3 (Figure2).

Figure 8 :

 8 Figure 8: Proposed semi-empirical estimation technique based on the first-order reversal curve, for estimating energy densities of (a) PZT C6, (b) PZT C9, (c) PMN-25PT, (d) PZN-8PT, (e) PMN 15, and (f) combined data under compressive stress values of 1 and 100 MPa.

Figure 9 :

 9 Figure 9: Estimation technique based on the first-order reversal curve for estimating the energy densities of (a) PZT C6, (b) PZT C9, (c) PMN-25PT, (d) PZN-8PT, and (e) PMN 15, and (f) combined data under compressive stress values of 1 and 30 MPa.

Figure 10 :

 10 Figure 10: Energy storage capacities of PZT C203, C6, C9, PMN-25PT, PZN-8PT, and PMN 15 under compressive stress values of 1, 30, and 100 MPa.

Figure 11 :

 11 Figure 11: Implementation of the Ericsson cycle using DC/DC switched converters.

Figure 12 :

 12 Figure 12: Simulation of the electrical interface based on the switched converters: (a) schematic, (b) time-domain waveform of the transducer capacitance C(t), (c) Vp and Vs, the transducer and energy storage terminals, (d) Vc1 and Vc2, the switch commands for charging and discharging circuits, and (e) energy cycle Q vs Vp, where Q denotes the transducer capacitance charge.

Figure 13 :

 13 Figure 13: Simulation of the electrical interface using self-synchronous converter: (a) schematic, (b) time-domain waveform of the transducer capacitance C(t), (c) currents in the diodes (steady-state current in diode D3 is the same as in diode D2), (d) voltages of the transducer Vp and of the storage capacitance Vs, (e) energy cycle (transient), (f) energy cycle (steady-state), where Q denotes the charge.

Table 1 :

 1 Properties of PZT C203, C6, C9, PMN-25PT, PZN-8PT, and PMN 15, as obtained from the suppliers

		C203	C6	C9	PMN-25PT	PZN-8PT	PMN 15
					[30][31]		[32][33]
	Coupling factor	71	76	69	91-94.9	94	60
	𝑘 33 × 10 -2
	Dielectric constant	1450	2130 6640	5400-5493	7700	~ 24000
	𝜀 33 𝑇 /𝜀𝑜
	Piezoelectric coefficient	325	427	718	1540-2060	2890	700
	𝑑 33 × 10 -12 (C/N)						Peak value

Acknowledgments

This work was performed under the framework of the ANR-FIESTA project, funded by the French Agence Nationale pour la Recherche, grant #ANR-20-CE05-0026, and under the framework of the International Research Network ELyT Global.

Data availability

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

Conflict of interests

The authors declare that they have no competing interests.