Physics-Based Modeling of Ferroelectric Hysteresis for Ceramic Capacitors in Inductively Coupled Microstimulators

In this article, we present a physics-based model for nonlinear and hysteretic ferroelectric capacitors in inductively coupled microstimulator circuits. The purpose of this model is to predict the system's dynamic response as a function of the dielectric material properties, with the aim of optimizing the performance in passive power control applications. We describe the workflow starting from the extraction of the dielectric material properties of commercial ceramic capacitors to the implementation into the physics-based model of the overall circuit. Selected capacitors were experimentally characterized at frequencies above 100 kHz by means of both small signals (5 mVrms, ±25 Vdc, and ±40 Vdc) and large signals (±25 Vac and ±40 Vac). Our results show that the developed model is well suited for accurately predicting the circuit's stimulation current for highly nonlinear capacitors and exhibits higher precision compared to commonly used models based on differential capacitance. Preliminary in vitro measurement results are finally described to provide proof of concept of the envisioned model-based passive power control in implantable microstimulators.


Physics-Based Modeling of Ferroelectric Hysteresis for Ceramic Capacitors in Inductively Coupled Microstimulators
Yves Olsommer , Frank R. Ihmig , and Gianluca Rizzello , Senior Member, IEEE Abstract-In this article, we present a physics-based model for nonlinear and hysteretic ferroelectric capacitors in inductively coupled microstimulator circuits.The purpose of this model is to predict the system's dynamic response as a function of the dielectric material properties, with the aim of optimizing the performance in passive power control applications.We describe the workflow starting from the extraction of the dielectric material properties of commercial ceramic capacitors to the implementation into the physicsbased model of the overall circuit.Selected capacitors were experimentally characterized at frequencies above 100 kHz by means of both small signals (5 mV rms , ±25 V dc , and ±40 V dc ) and large signals (±25 V ac and ±40 V ac ).Our results show that the developed model is well suited for accurately predicting the circuit's stimulation current for highly nonlinear capacitors and exhibits higher precision compared to commonly used models based on differential capacitance.Preliminary in vitro measurement results are finally described to provide proof of concept of the envisioned model-based passive power control in implantable microstimulators.

I. INTRODUCTION
I N the last few years, a new generation of leadless and battery- free implantable microstimulators has been approved for the market.An example of such systems is the RENOVA (BlueWind Medical LTD, Herzliya, Israel) tibial nerve stimulator, which received the European Conformity mark in 2016 for the treatment of overactive bladder [1].Its stimulation pulse amplitude, pulse width, and frequency can be up to 9 mA, 800 µs, and 40 Hz, respectively [2].A well-known way to improve reliability and enable further miniaturization of implantable microstimulators Yves Olsommer is with the Department of Molecular and Cellular Biotechnology/Nanotechnology, Saarland University, 66123 Saarbrücken, Germany (e-mail: yves.olsommer@ibmt.fraunhofer.de).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TPEL.2024.3361075.
Digital Object Identifier 10.1109/TPEL.2024.3361075consists of reducing their complexity and number of components [3].As a result, such microstimulators are simply designed to serve only one purpose, namely to convert an induced voltage into a stimulation pulse.In other words, the stimulation current is regulated indirectly by inductive power transfer (IPT) control.
Based on this principle, the leadless and battery-free bilateral hypoglossal nerve stimulator Genio (Nyxoah SA, Mont-Saint-Guibert, Belgium) was developed [4].This system was approved for the European market in 2019 for the treatment of therapy-resistant obstructive sleep apnea [5].In the same vein, a bioresorbable implantable pacemaker was developed for the treatment of postoperative temporary arrhythmias [3].Here, the microstimulator consists only of a coil, a p-i-n diode and the stimulation electrodes [3].The high parasitic capacitance of the p-i-n diode forms a resonant circuit with the coil, and the long reverse-recovery time is used to rectify the induced voltage in a stimulation pulse [3].The advantage of these leadless and battery-free microstimulators is twofold: 1) to prevent adverse events (AEs) related to lead breakage and migration [6] as well as battery depletion [7], and 2) to reduce the invasiveness of the stimulator to the patient, the complexity of the surgical procedure, and the risk of infection [8], [9].Additional surgeries due to AEs are costly and put patients who are already part of a frail population at additional risk for complications [10].Due to their increased reliability, leadless, and battery-free microstimulators are gaining popularity.A drawback of such stimulators, however, is that they operate as voltage (rather than current) sources.Consequently, the amplitude of the stimulation current is directly related to the inductive coupling factor and transmitted power, as well as the electrode impedance.In addition, there is no space to integrate the additional components required for conventional closed-loop control, such as sensors, microcontrollers, and communication channels, to set and maintain the stimulation current amplitude at a safe and stable value range.Therefore, a proper inductive power transfer control method is required, as an excessively high stimulation current amplitude could irreversibly damage the stimulation electrodes and the surrounding tissue [11].
In our research work, we aim to realize a novel approach for passive control of the stimulation current.As in [3] and [4], the microstimulator should consist of only a low number of components such as a coil, capacitors, resistors, and diodes.We intend to exploit the nonlinearity of ferroelectric dielectrics in class 2 multilayer ceramic chip capacitors (MLCCs) in the microstimulator's resonant circuits, particularly the plateaus occurring in their voltage-charge hysteresis curve, to selfregulate the stimulation current to a safe and stable value range.To identify the ideal shape of the ferroelectric hysteresis curve required to achieve this goal, a model-based approach is pursued.In the following Section I-A, we first provide an overview of active and passive IPT control strategies that are relevant for implantable systems.Common models of ferroelectric hysteresis and their limitations are presented in Section I-B.

A. IPT Control Strategies
Over the years, numerous control strategies have been developed to overcome the major limitations of IPT, related to the sensitivity of the inductive link to application-related fluctuations in load and coupling factor.Control methods are essential to adjust the voltage or current for proper operation of secondary-side inductively powered devices and to optimize the efficiency of the IPT, regardless of changes in the inductive coupling factor and load.Typically, conventional control methods require the input and output quantities to be actively controlled and measured, respectively.Common control parameters for IPT include the supply voltage of the primary-side output stage as well as the frequency [12], pulse width [13], and phase shift [14] of the output stage control signals used to drive the primary-side resonant circuit.Alternatively, the primaryand secondary-side resonant circuits can be actively tuned by trimming capacitors adjusted by a small stepper motor [15] or by switching capacitors using specific switching patterns, as shown in [16] and [17].Unfortunately, as described in [13], [14], [15], and [17], active control methods require additional components such as microcontrollers, communication channels, and sensors, increasing the complexity of the secondary-side circuit, space requirements, and cost [13], [18], [19].In addition, robust communication channels are an essential part of IPT closed-loop control as seen in [13], [14], [16], and [19].Any delay in data transmission or interruption of communication channels will result in increased response time or even failure of the closed-loop control [20], [21].The microstimulators in [3] and [4] described above are free of all these components.These consist of two main parts: 1) a resonant circuit for inductive power harvesting and 2) a rectifier to convert the induced voltage into a stimulation pulse.The amplitude of the stimulation current is controlled by IPT using open-loop control methods.Unlike active control approaches, passive control methods employ compensation circuits to optimize the transfer function of the IPT [22].To reduce the number of components needed to regulate the voltage on the receiver side and, thus, reduce the volume of the implantable electronics, Lin et al. [22] proposed an LCC/PS compensation circuit to design an IPT system with a load-independent output voltage at a frequency of 6.78 MHz.At a constant inductive coupling factor, the voltage could be kept within a range of roughly 4.6 to 4.9 V among a load range of 500 to 2500 Ω [22].However, if the distance between the two coils changes from 7 cm to 3 cm, the output voltage increases from 3 to 7.5 V [22].Similarly, a series-series compensation circuit has been implemented in [18] to achieve a load-independent current at a frequency of around 100 kHz.However, to maintain the condition of a load-independent current despite the application-related variation of the coupling factor, an additional active closed-loop control is required [18].
In [23] and [24], we introduced an alternative approach to develop a passive power control strategy by using ferroelectric materials in MLCCs as smart materials.We aim to exploit the MLCCs' nonlinear behavior to self-regulate the stimulation current within a safe and stable amplitude range by compensating for application-related fluctuations in the inductive coupling factor.An advantage of this method is that, ideally, only one MLCC with a defined nonlinearity is required in the receiver resonant circuit.

B. MLCC Modeling Approaches
Nowadays, ferroelectric MLCCs prevail in the market due to their high volumetric efficiency and reliability [25].These MLCCs are used in entertainment electronics, such as televisions and smartphones, as well as in applications with high reliability requirements, e.g., automotive powertrain and safety equipment as well as implantable medical devices [25], [26].In general, the nonlinear characteristics of the ferroelectric materials in MLCCs are rather considered as a disadvantage.However, there exist approaches that attempt to turn this nonlinearity into a benefit [27], [28], [29], [30].The most important application for ferroelectric materials is certainly in memory devices, such as ferroelectric random-access memories [31].The potential use of ferroelectric materials in various applications has led to an increased need for simulation models.The models of nonlinear MLCCs must provide reliable results, be numerically efficient, and should be easy to combine with other models.
The Preisach theory, originally used to model the hysteresis of ferromagnetic materials, has been used in [32], [34], [35], and [36] to develop a model for ferroelectric materials.Such models can be easily implemented into other models and are computationally efficient [33].However, due to the phenomenological nature of those models, the dynamic domain switching is either neglected [32], [36] or reproduced by time constants [34], [35].Thus, these models lack the physics behind dynamic domain switching.Further models that can be considered computationally efficient are those based on the phenomenological Landau-Ginzburg-Devonshire theory [37].According to this theory, the voltage across a ferroelectric material can be described as a function of electric charge by a fifth-order polynomial [38].Experimentally determined polynomial coefficients are used to model the nonlinear properties of ferroelectric materials without considering the underlying physics [38], [39].
SPICE models of voltage-dependent ferroelectric MLCCs are also reported in the literature [27], [28], [29], [40].Such models are popular among developers in the field of power electronics as they are readily available, user-friendly, and quick to set up and run [40].Since these models mimic the nonlinear behavior of ferroelectric MLCCs using current and voltage sources, the effect of the dielectric material properties on the nonlinear behavior of the MLCCs remains unclear.As an example, in [28] and [29], the nonlinear behavior of ferroelectric MLCCs was modeled with an equivalent experimentally measured differential capacitance, neglecting ferroelectric hysteresis.Recently, models were developed to describe the nonlinear properties of single-crystal piezoelectric transducers on a physical level [41].Conventional models for such devices are based on phenomenological approaches, whereby the insight into the physics behind the nonlinear properties is missing.In the physics-based model presented in [41], instead, ferroelastic and ferroelectric hysteresis are modeled by reproducing the switching processes of electric dipoles through a thermodynamic transition probability approach.The adopted free-energy framework makes it possible to easily extend the constitutive equations, and adapt them to further modeling purposes.Moreover, the energy formulation allows the device model to be easily coupled to other dynamical (i.e., electro-mechanical) systems in a thermodynamically consistent fashion.Inspired by the research in [41], we present a physics-based model of ferroelectric hysteresis for MLCCs.In order to quantitatively describe their highly nonlinear response, the single-crystal ferroelectric modeling approach from [41] is extended for the first time by incorporating polycrystalline mechanisms originally developed in [42] for shape memory alloy wire actuators.For a better understanding of the free energy framework and the thermodynamic transition probability approach, the reader is referred to the following literature [43], [44], [45].
The rest of this article is organized as follows.In Section II, we describe the selected MLCCs and the measurement setups used for the characterization experiments.Section III presents measured differential capacitances first, introduces the physicsbased model of ferroelectric MLCCs subsequently, and concludes with a presentation of the model predictions and with preliminary in vitro measurement results related to the envisioned application.The results based on the proposed model are then discussed and compared with existing models in Section IV.Finally, Section V concludes this article.

A. Description of the Selected MLCCs
Out of 40 MLCCs, seven nonlinear class 2 MLCCs from six different manufacturers (KEMET, AVX, TDK, Walsin, Multicomp Pro, Murata) with four case sizes and three temperature coefficients were investigated.A class 1 (C0G) MLCC, labeled #1, whose capacitance is voltage independent, was used as reference measurements.A detailed description of the selected MLCCs is given in Table I.

B. Measurement of Differential Capacitance
The differential capacitance of the selected MLCCs was measured using the precision impedance analyzer 4294A (Agilent Technologies Inc., Santa Clara, PA, USA, 4294A R1.11 Mar.25,2013) and the test fixture 16047D (Hewlett-Packard, Palo Alto, CA, USA).A small signal of 5 mV rms , superimposed to a bias voltage of ±40 V dc or ±25 V dc , was generated with the precision impedance analyzer [46].Prior to measurements, the impedance analyzer was calibrated for an open-and short-circuited test fixture.Each MLCC was measured twice by varying the bias voltage from the lowest to the highest value and vice versa.The relationship between the differential capacitance and the bias voltage was obtained from the average of both measurements.

C. Measurement of Ferroelectric Hysteresis
A measurement setup was realized to acquire the ferroelectric hysteresis of MLCCs at voltages and frequencies above 40 V ac and 100 kHz.The measurement setup consists of an output stage, which is inductively coupled with a Sawyer-Tower (S-T) circuit [47], [48].The voltage across a nonlinear and linear capacitor was measured with the oscilloscope MDO4104-6 (Tektronix Inc., Beaverton, OR, USA) and two probes TT-MF312-2-6 11020-2-6 (TESTEC Elektronik GmbH, Frankfurt, Germany).The electrical charge flowing through the series-connected nonlinear and linear capacitors is proportional to the capacitance of the linear capacitor.Measurement data were transferred from the oscilloscope to a computer using OpenChoice Desktop software from Tektronix.

D. Measurement Setup to Characterize Inductively Coupled Microstimulator Circuits
The inductively coupled microstimulator circuit was characterized using the measurement setup shown in Fig. 1.The microstimulator (right) is inductively powered through a transmitter (left).The transmitter consists of a half-bridge resonant converter GS-EVB-HB-61008P-ON (GaN Systems Inc., Ottawa, Canada) and a series resonant circuit.The half-bridge resonant converter is driven by a DG5102 (Rigol Technologies Inc., Suzhou, China) waveform generator and powered by a U8031A (Agilent Technologies Inc., Santa Clara, CA, USA) power supply.L 1 and R 1 represent the inductance and loss resistance of the transmitter coil (part number 760308100110).In a frequency range of 360 kHz and 463 kHz, L 1 is 6.23 µH and R 1 ranges from 0.22 to 0.26 Ω.To tune the series resonant circuit to resonate at different frequencies, the capacitance C 1 can be adjusted by a parallel connection of several high-voltage foil capacitors of type FKP 1 (WIMA GmbH and Co. KG, Mannheim, Germany).
For inductive power harvesting, the microstimulator circuit consists of a parallel resonant circuit (L 2 , R 2 , C 2 ), where C 2 corresponds to the nonlinear MLCC to be modeled.L 2 is 3.9 µH and R 2 ranges from 0.25 to 0.28 Ω in a frequency range of 360 kHz and 463 kHz (part number 760308101104).A half-wave rectifier with diode D (LL4148, onsemi LLC, Scottsdale, AZ, USA) and capacitor C 4 generate the stimulation pulse [see Fig. 2(a)].Thus, the width and repetition rate of the stimulation pulse are directly related to the duration and repetition rate of the induced voltage, i.e., these stimulation parameters are set by the transmitter.The amplitude of the stimulation current, however, depends on a large number of parameters.These include the inductive coupling factor k, the electrode impedance R L , the power fed into the transmitter, the frequency of IPT, and the nonlinear  voltage-dependent capacitance of C 2 .The electrode impedance, normally present in the application and resulting from the electrical properties of the stimulation electrodes themselves and the surrounding tissue, has been replicated by an ohmic load R L [11], [49].In the target application, the microstimulator circuit will be inductively powered by a battery-operated wearable.
Since the supply voltage is limited, a series resonant circuit is advantageous because it requires smaller voltage swings at its input as the phase of the inductor and capacitor voltages cancel at resonance [50].A parallel resonant circuit on the microstimulator side amplifies the induced voltage and helps to overcome the turn-ON voltage of the rectifier diode, which is particularly advantageous at low inductive coupling factors [50].
To characterize the inductively coupled microstimulator circuit, the stimulation current I Stim was measured as a function of the power P 1 fed into the transmitter.The power P 1 was calculated from the root mean square of the product of the current flowing through the series resonant circuit and the voltage at the output of the half-bridge converter in a steady state [see Fig. 2(c)].Current and voltage were measured using a TCP312A current clamp with its amplifier TCPA300 from Tektronix and a TT-MF312-2-6 11020-2-6 probe from TESTEC.To determine I Stim , the voltage at the load R L was measured [see Fig. 2(b)].The transmitter coil (L 1 and R 1 ) and the microstimulator circuit coil (L 2 and R 2 ) with a diameter of about 5 cm and 2 cm, respectively, were fixed on a spacer [see Fig. 1(a)].A coupling factor k of 5.6 % was measured using the impedance analyzer Agilent 4294A.For proof-of-concept measurements, the pork flesh in Fig. 1(b) was used as a spacer.A maximum k of 13 % was reached.

A. Differential Capacitance
The measured differential capacitances are illustrated in Fig. 3.An increasing degree of nonlinearity can be observed for MLCCs #1 through #8.To quantify this degree of nonlinearity, the differential capacitance at 0 V was divided by that at 25 V.For example, since the differential capacitance of linear MLCC #1 is independent of the voltage, the ratio is 1.00.The more the differential capacitance decreases as a result of increasing voltage, the larger the ratio.For MLCCs #2 through #8, the following ratios are obtained: 1.03, 1.24, 1.70, 2.61, 3.94, 5.05, and 7.58.

B. Physics-Based Model of Ferroelectric Hysteresis for MLCCs
The model analytically describes the voltage-charge hysteresis of the nonlinear MLCC and is formally represented as a set of ordinary differential equations (ODEs) expressed in causal impedance form (i.e., current input-voltage output).The model is grounded on the statistical thermodynamic approach presented in [41].We assume that the ferroelectric material is composed of a combination of dipoles, whose orientation with respect to the direction of the electric field is assumed to be either positive (+) or negative (-).The component of the polarization vectors of such dipoles along the direction of the electric field E, denoted as P + and P -for positively and negatively oriented dipoles, respectively, is given as follows: where ϵ 0 denotes the vacuum permittivity, χ e is the electric susceptibility of the material, while P +0 and P -0 represent the remnant polarizations of the two types of dipoles.For convenience, we introduce quantities x + and x -describing the volume (i.e., phase) fractions associated with dipoles having positive and negative orientations, respectively, with Then, the total average polarization P can be expressed as a function of P + and P -according to the following equation: By exploiting (1) and (2), we can conveniently rewrite (3) as Based on (4), we can then express the component of the average electric displacement field D along the direction of E as where ϵ r = χ e + 1 denotes the materials' relative permittivity.
According to (5), the knowledge of the phase fraction x + is needed in order to relate D and E. The evolution of x + over time  can be computed by integrating the following ODE [41]: where transition probabilities p +− and p −+ are derived based on thermodynamic considerations and are expressed as follows: In (7), τ x represents a time constant associated with thermal activation, V l is the volume of a mesoscopic material crystal, k B is the Boltzmann constant, T is the material temperature, while Δg +-and Δg -+ are the Gibbs free-energy density barriers associated to the x + → x -and x -→ x + transformations, respectively.By manipulating the driving force equations adopted in [41], it can be shown that the last two quantities can be analytically expressed as a function of E according to where E +0 and E -0 represent the remnant fields of the hysteresis.Such remnant fields are responsible for triggering the phase transformation between positively and negatively oriented dipoles, i.e., transformation x + → x -is favored if E ≤ E -0 while transformation x -→ x + is favored whenever E ≥ E +0 .Inspired by [42], we let such remnant fields be dependent on phase fraction x + .In here, we propose the following expressions: where e +0 , e +1 , e +n , e -0 , e -1 , e -n are constitutive material constants.Using ( 10)-( 11), we can capture inhomogeneities within the ferroelectric material by means of a representative single-crystal element, which "sees" a phase fraction-dependent effective remnant field, which changes progressively as more dipoles change orientation.As proven in [42], the introduction of this mechanism allows to implement a polycrystalline (i.e., smooth) hysteresis outer loop in a highly numerically efficient way.It is remarked that the considered model only provides an approximated description of the minor loops, and a further modification of the approach is needed to include smooth minor loops, as discussed in detail in [42].Since this modification is involved from the numerical standpoint, it is omitted from this work.To convert the average material quantities E and D into corresponding measurable quantities, namely the voltage u C2 and the charge q C2 , we use the following equations: where A el and l el represent the electrode surface area and the dielectric thickness, respectively.The applied current i C2 can then be related to the charge q C2 via the simple relationship By collecting ( 5)-( 7) and ( 12)-( 13), we can finally write the overall model of the MLCC as a set of two nonlinear ODEs in impedance form Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where i C2 is the input, u C2 is the output, while D and x + define the state variables of the model.Quantities Δg +-and Δg -+ appearing in ( 14) can be explicitly computed via ( 8)- (11).Before moving forward, we state the following property of ( 14).Lemma 1: If the values of V l and τ x -1 are chosen sufficiently large (i.e., τ x sufficiently small), then the electric field E tightly satisfies the following equation over the solutions of ( 14): Proof of Lemma 1: The proof can be obtained straightforwardly by recalling the result presented in [51], Section 4.1, i.e., by noting that ( 7)-( 9) behave approximately as threshold functions, replacing the stress with the electric field and the strain with the electrical displacement, and using the similarity among the phase fraction differential equations.
The conditions of Lemma 1 are usually satisfied for physically meaningful values of V l and τ x .Lemma 1 is instrumental in proving the following structural property of model (14).
Proposition 1: If the assumptions of Lemma 1 hold true, and under additional assumptions P +0 ≥ P -0 , E +0 (x + ) ≥ E -0 (x + ) x + ࢠ [0, 1], E +0 (y) ≥ E +0 (x) y ≥ x, and E -0 (y) ≥ E -0 (x) y ≥ x, then model ( 14) defines a passive operator at the port i C2 → u C2 , i.e., there exist a storage function Ψ(D, x + ): R2 → R+ such that the following inequality holds on the trajectories of ( 14) [52]: Proof of proposition 1: To prove the result, we propose the following candidate storage function: which satisfies Ψ(D, x + ) ≥ 0 for every D and x + since P +0 − P -0 ≥ 0 while E +0 (x + ) and E -0 (x + ) are assumed to be monotonically increasing functions of x + .By computing the time derivative of (17) and exploiting (14), we obtain Since the assumptions imply that A el l el (P +0 − P −0 ) ≥ 0, then ( 16) is satisfied if and only if the following inequality holds true: As Lemma 1 holds true, we can replace (15) in (19), and obtain which permits us to conclude that ( 16) holds true for any arbitrary trajectory of system (14).This concludes the proof.The importance of Proposition 1 is clarified in the sequel.The formal equivalence of model ( 14) to an electrical impedance (i.e., input current-output voltage) makes it possible to simply interconnect the physics-based MLCC model ( 14) to any external electric network formulated as an electrical admittance (i.e., voltage input-current output).Mathematically, coupling (14) to a generic electric circuit can be expressed as a negative feedback interconnection among two causal dynamic systems, i.e., the nonlinear impedance model ( 14) and the model of the external network, the latter generally consisting of a passive circuit (as in Fig. 1).It is well known from nonlinear systems theory that the negative feedback interconnection among two subsystems involving nonlinear and/or hysteretic components may lead to an unstable behavior [53], thus potentially hindering the simulation framework proposed in this article.However, in case both the MLCC model and the external network are passive at their corresponding coupling port, then their negative feedback interconnection will always result in a passive (and, thus, in practice, stable) system [52].The proposed MLCC model can then be effectively employed in simulations of nonlinear passive circuits without any risk of instabilities, thus making it a robust numerical tool for reliable predictions.

C. Material Parameter Extraction
For convenience, we report all the free parameters of model ( 14): V l , τ x , A el , l el , ϵ r , P +0 , P -0 , e +0 , e +1 , e +n , e -0 , e -1 , e -n .A systematic procedure for the calibration of such parameters is outlined in this section.First, since V l and τ x are practically hard to measure, they are usually fixed a priori based on the expected range for those parameters, as known from the literature.Since manufacturers do not disclose the exact structure of the MLCCs geometry and material properties, nominal parameters A el and l el are determined considering the package size of the MLCCs specified in the datasheets resulting in an "equivalent" material model.Additional analyses, such as scanning electron spectroscopy and energy dispersive X-ray spectroscopy, would be necessary to address this lack of information, although they are beyond the scope of this article.All the remaining parameters can then be identified via an experimentally measured voltage-charge hysteresis loop.To begin with, we can convert the experimental u C2 -q C2 hysteresis curve into the corresponding E-D loop based on (12).Once this plot is available, we can approximate the asymptotic behavior of the upper and lower branches of the hysteresis with two line segments having the same slope representing ϵ 0 ϵ r and different offsets, denoting P +0 and P -0 , respectively.Note that P +0 = −P -0 ≥ 0 commonly holds true.Next, we can use (5) to express x + as a function of D, E, and the identified material constants.If this calculation leads to an x + , which violates condition x + ࢠ [0,  1], then one must repeat the previous step and make a more suitable choice of the asymptotic lines used to determine ϵ r , P +0 , and P -0 .Once x + is known, we can assume that (15) holds true and, in turn, identify the parameters describing E +0 (E -0 ) by performing a least-square fitting of e +0 , e +1 , and e +n based on (10) and the increasing branch (e -0 , e -1 , and e -n based on (11) and the decreasing branch) of the obtained x + -E curve.In practice, it is noted that choosing e +0 = −e -0 ≥ 0, e +1 = e -1 ≥ 0, and e +n = e -n ≥ 0 generally leads to a satisfactory result, and automatically ensures that the conditions of Proposition 1 are met.After this step, all parameters have been successfully identified.The measured and modeled hysteresis were plotted in Fig. 4. The parameters of the "equivalent" material model are given in Table II.

D. Physics-Based Model of Ferroelectric MLCCs in Inductively Coupled Microstimulators
Inductively coupled microstimulators are modeled in Mathcad Prime 3.1 (PTC Inc., Boston, MA, USA) by a set of second-order ODEs.The ODEs can be derived from the circuit diagram shown in Fig. 1.
The square wave output signal of the transmitter half-bridge was approximated with its fundamental frequency and corresponding DC components.The current flowing through the rectifier diode D is modeled according to the nonlinear voltagecurrent relationship specified in the manufacturer's datasheet.
The inductive coupling factor between the coil of the transmitter and the microstimulator is expressed by k.
The hysteretic relationship between charge q C2 and voltage u C2 of the ferroelectric MLCC C 2 is obtained according to model (14).Note that, although (14) considers the current i C2 as the natural input signal, in the actual implementation we have instead used the charge q c2 as input of the C 2 sub-block.This is possible due to the particular structure of the resulting ODEs, which allows to couple the ferroelectric hysteresis model with the inductively coupled microstimulator circuit with no violation of implementation causality and without the need to include q C2 as an additional state variable.The backward differentiation formula method was used in Mathcad to solve the system ODEs.The convergence tolerance was set to 10 -7 .

E. Model Validation for Linear MLCC
Before modeling nonlinear ferroelectric MLCCs, the model was validated with a linear MLCC.For this purpose, a class 1 MLCC (C0G) of 47 nF with a rated voltage of 200 V, denoted as #1, was used.In Fig. 5, the measured stimulation current I Stim is plotted in black as a function of the power P 1 fed into the transmitter.P 1 was kept in a safe range to avoid excessive voltages across the MLCCs.For this reason, different ranges of P 1 can be observed for each case in Fig. 5.The calculations using the model based on the differential capacitance and the physicsbased model presented in this article are plotted in green and red, respectively.To quantify the consistency between measurements and calculations with both models, the percentage R-squared [54] value was calculated as follows: where I Stim , I Stim , and I Stim are the measured, calculated, and mean value of the measured stimulation current, respectively.The computed values of R 2 are summarized in Table III.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE III R 2 VALUES BETWEEN THE MEASUREMENTS AND BOTH THE DIFFERENTIAL CAPACITANCE AND PHYSICS-BASED MODEL
As expected, in the case of a linear MLCC, a high consistency between the measurements and both models is achieved, resulting in a R 2 value of 99.9 %.Dielectrics of class 1 MLCCs consist of paraelectric materials [55], whose dipoles spontaneously align with the applied electric field.Accordingly, the definition of the differential capacitance (C = dq / du) measured with the impedance analyzer and the absolute capacitance (C = q / u) measured with the S-T circuit are equivalent [40].In addition, for the linear MLCC #1, the power transfer efficiency (PTE) of the inductive link is approximately constant at 14.0 % over a range of P 1 from 2.3 W rms to 13.0 W rms and at an inductive coupling factor of 5.6 % and a frequency of 366 kHz.The PTE is defined as the ratio of the power supplied to the load R L to the power fed into the primary resonant circuit (L 1 , R 1 , and C 1 ).

F. Modeling of Ferroelectric MLCCs
The results for seven nonlinear ferroelectric MLCCs are shown in Fig. 5.The calculations for MLCC #7 using the differential capacitance model differ significantly from the measurements.For P 1 > 2.5 W rms , lower values of I Stim are calculated than measured.In contrast, the physics-based model well predicts the behavior of the system, even the steep increase in I Stim for P 1 < 2.5 W rms .A similar prediction of the system behavior is achieved for MLCC #8.For both MLCCs #7 and #8, the differential capacitance model yielded an R 2 value of 72.0 % and 51.3 %, respectively, while the physics-based model yielded an R 2 value greater than 90 %.With increasing P 1 , an increasing PTE was observed with MLCCs #7 and #8.The maximum PTE achieved with MLCCs #7 and #8 is 6.7 % and 8.2 % at a P 1 of 1.9 W rms and 0.3 W rms , respectively.If P 1 is further increased to 18.8 W rms and 45.7 W rms for MLCCs #7 and #8, respectively, the PTE drops from 6.7 % and 8.2 % to 2.5 % and 0.8 %.Interestingly, the calculations obtained using the physics-based model and the differential capacitance model of MLCC #6 show similar results, with respective R 2 values of 98.9 % and 94.3 %.In addition, it can be observed that the PTE reaches a maximum value of 9.2 % at a P 1 value of 1.2 W rms .The efficiency drops to 0.9 % when P 1 is increased to 159.7 W rms .The physics-based model provides a much better prediction of the system behavior for MLCC #5 compared to the differential capacitance-based model.Nevertheless, the calculated amplitude of I Stim is lower than the measured one for P 1 > 5 W rms and P 1 < 2.5 W rms .Using the differential capacitance-based model, the R 2 value was only 7.6 % compared to 87.4 % using the physics-based model.The sharp increase in I Stim in Fig. 5 can also be observed in the PTE.The PTE lies between 2.1 % and 2.7 % in a P 1 range between 0.3 W rms and 4.3 W rms .At a P 1 of 5.5 W rms , the PTE increases from 2.6 % to 6.1 % and rises to 6.9 % at a P 1 of 7.0 W rms .By increasing P 1 to 34.2 W rms , the PTE drops to 2.6 %.
For MLCC #3, the modeled amplitude of I Stim using the physics-based model is lower than the measured one for P 1 < 25 W rms .The steep increase in I Stim already occurs at 75 W rms in the physics-based model compared to 125 W rms in the measurements.Unlike the previous MLCCs, the PTE with MLCC #3 initially drops from 3.3 % to 1.0 % over a range of P 1 between 0.3 W rms and 31.7 W rms .For a P 1 range of 41.7 W rms to 67.6 W rms , the PTE remains nearly constant at 0.9 %.As The same observations hold for MLCC #4.The modeled amplitude of I Stim using the physics-based model is lower than the measured one for P 1 < 10 W rms .The steep increase in I Stim occurs for a slightly lower value of P 1 in the model than in the measurements.Here, the R 2 value was 66.9 % using the differential capacitance-based model, compared to 95.1 % using the physics-based model.At first glance, the PTE drops from 4.1 % to 2.1 % in the range of P 1 between 0.3 W rms and 8.9 W rms .In the P 1 range between 13.4 W rms and 16.6 W rms , the PTE increases to a value of 8.4 % and then decreases to a value of 6.8 % as P 1 increases from 16.0 W rms to 24.4 W rms .
In the case of MLCC #2, the above discrepancies between measurements and calculations, based on both differential capacitance and physics-based models, are illustrated in an even clearer manner.At first glance, the capacitance of MLCC #2 exhibits a linear behavior (see Figs. 3 and 4).However, as shown in Fig. 5, the calculations using the differential capacitance and physics-based models differ considerably from each other as well as from the measurements.The negative R 2 value for MLCCs #2 and #3 shows that the differential capacitance model does not reflect the behavior of the system in any way.Using the physics-based model, the R 2 value for MLCCs #2 and #3 was 43.5 % and 67.1 %, respectively.With MLCC #2, the PTE decreases from 5.1 % to 0.8 % as P 1 increases.

G. In Vitro Proof-of-Concept Measurements
As a preliminary investigation of the feasibility of the envisioned passive control of the stimulation current, the physicsbased model was used to further investigate the impact of the ferroelectric hysteresis on the IPT.The capacitor C 2 and the coil (L 2 and R 2 ) determine the resonant frequency of the microstimulator.First, the resonant frequency was measured with a small signal (5 mV rms ) by using the impedance analyzer.Then, the physics-based model was used to estimate the resulting stimulation current as a function of the IPT frequency and the coupling factor for a given power level in the primary-side resonant circuit.The sensitivity of the stimulation current I Stim as a function of the coupling factor was investigated with the highly nonlinear MLCCs #6, #7, and #8.These MLCCs exhibit the highest voltage-dependent capacitance (see Table IV, with C(u 90% ) / C(0 V) > 4).The lowest sensitivity, i.e., the strongest stabilization of I Stim , was achieved with the MLCC #6.By operating the system 10 kHz above its resonant frequency and inducing a voltage across MLCC #6 between 9.3 V rms and 9.7 V rms , I Stim was kept almost constant (see Fig. 6).For coupling factors ranging from 5.9 % to 11.9 %, we measured Fig. 6.Measured and modeled stimulation current I Stim as a function of inductive coupling factor k for linear MLCC #1 (black) and nonlinear MLCC #6 (red).The power level of the primary-side transmitter remains the same for both measurements.The system is operated 10 kHz above its resonant frequency.
an I Stim value between 11.0 mA and 11.7 mA, in contrast to a range of 11.8 mA to 15.7 mA using the linear MLCC #1.
Fig. 6 shows that the physics-based model adequately predicts the overall system behavior by using MLCC #6 to implement the passive control approach.With a k of 9.4 %, an I Stim of 11.37 mA was measured compared to a modeled I Stim of 11.03 mA.It should be noted that the outer loop hysteresis of MLCC #6 was measured at a voltage of 15.1 V rms and that the predicted stabilization of I Stim at a k above 6 % occurs at voltages between 9.3 V rms and 9.7 V rms .Since our physics-based model provides a rough approximation of the inner loops, system behavior is less accurately modeled at voltages below 15.1 V rms .This is particularly evident in the sharp increase in I Stim at k = 4.5 % in the measurements and k = 3 % in the model.Nevertheless, the overall system behavior in the region of interest (i.e., the current plateau) is well predicted by the model.

IV. DISCUSSION
A good consistency between the physics-based model and measurements was achieved for the nonlinear ferroelectric MLCCs #6, #7, and #8, with R 2 values of 98.9 %, 90.7 %, and 95.6 %, respectively.Nevertheless, from the above results, it appears that the system behavior could not be modeled properly for the nonlinear MLCCs #2, #3, #4, and #5.These discrepancies will be discussed in this section.
The hysteresis of the MLCCs was measured with the S-T circuit in different voltage ranges.For all the nonlinear MLCCs investigated, it can be observed that the minor loops of the hysteresis deviate from the outer loop.This is attributable to the dielectric manufacturing process.Domain switching is hindered by structural defects such as vacancies, interstitials, and dopants [56].The size of the domains is limited by the grain size used to manufacture the dielectrics [57].The probability of domain pinning increases with domain size [57].At lower ac voltages, the pinned domains will not align with the electric field and therefore disturb the kinetics of domain switching [57].At higher ac voltages, however, the previously fixed domains will also align with the electric field, causing the dielectric constant to Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
rise sharply [57].Further increases in the ac voltage cause the dielectric constant to decrease [57].
To give an idea of how the hysteresis minor loops change within a given voltage range and at different voltage ranges, the slope of the hysteresis, namely the differential capacitance (dq / du), is shown at 0 V and 90% of the applied voltage u in Table IV.
The ratio of capacitance at 0 V in a voltage range of ±5 V to that in a voltage range of ±25 V corresponds to 1.20, 1.23, and 1.10 for MLCCs #6, #7, and #8, respectively.
Therefore, it is clear that the shape of the hysteresis is predominantly determined by a change in the applied voltage rather than by a change in the voltage range itself.Note that a high consistency between the physics-based model calculations and the measurements was achieved for these MLCCs (R 2 > 90 %).
In contrast, the ratio of capacitance at 0 V in a voltage range of ±5 V to that in a voltage range of ±25 V for the respective MLCCs #2, #3, #4, and #5 is 0.86, 0.83, 0.92, and 0.91.This ratio is of the same order of magnitude as the ratio of the capacitance at 0 V to that at 90 % of u in a voltage range of ±25 V. Here, the shape of the hysteresis is determined both by a change in the applied voltage and by a change in the voltage range itself.
In our physics-based model, the material parameters of the ferroelectric MLCCs are extracted at the maximum applied voltage, corresponding to the outer loop of the hysteresis.Both the nonlinear hysteresis shape within a given voltage range and the dependence of the hysteresis minor loops on the voltage range impact the parallel resonant circuit in the microstimulator.When the detuning of the resonant circuit results rather from the nonlinear shape of the hysteresis outer loop (see MLCC #6, #7, and #8), i.e., when domain pinning can be neglected, the system behavior can be modeled accurately.However, for MLCCs with a lower degree of nonlinearity (see MLCC #2, #3, #4, and #5), the inner loops can no longer be neglected.Since the current version of our physics-based model only provides a rough approximation of the inner loops, the accuracy decreases in the case of MLCCs #2, #3, #4, and #5.The accurate description of minor hysteresis loops is possible by extending the physics-based model with the scaling approach from [42].However, since our passive control approach is based on highly nonlinear MLCCs at higher ac voltages (see the current plateau in Fig. 6), where the influence of the pinned domains on the domain switching kinetics is negligible, the minor loops can be neglected.
In almost all investigated nonlinear MLCCs, a significant deviation of the calculations using the differential capacitance model from the measurements could be observed.The variation of the bias voltage across the MLCC over the entire voltage range when measuring the differential capacitance occurs in a time window that is much longer than the period of the measurement signal itself.Note that, for a given bias voltage, the measurement time is at least about 1000 times longer than the period of the measurement signal [58].Hence, it can be assumed that the domains in the ferroelectric dielectric have a static orientation during the moment when the differential capacitance is measured.In contrast, the S-T circuit employs a large signal to measure the absolute capacitance.Here, the effect of domain switching kinetics on capacitance is included in the measurements.For this reason, differential capacitance is not sufficient to properly model ferroelectric MLCCs.
One exception, however, is MLCC #6.From Fig. 5, it can be seen that, although the physics-based model gives an overall better prediction of the system behavior, the calculations using the differential capacitance model are also close to the measured values of I Stim .These results suggest that the domain switching kinetics do not seem to play a major role in the resulting capacitance for MLCC #6 at a frequency of 360 kHz.
In Fig. 6, a maximum PTE of 10.6 % was measured for the linear MLCC #1 at a coupling factor k of 10.4 %.A similar PTE of 10.3 % was measured with the nonlinear MLCC #6 at a k of 10.7 %.A maximum PTE of 11.4 % was achieved with the MLCC #6 at a k of 11.9 %.For comparison, in [59], an efficiency of 18.9 % was achieved at a coupling factor of about 6 % and a frequency of 6.78 MHz.This has to be taken into account when applying this concept to other IPT systems.A low PTE can be problematic for continuous IPT in terms of both tissue heating and overall system energy balance.However, it is important to remember that power is delivered to the implantable microstimulator in bursts.The pulse width and repetition rate correspond to those of the stimulation pulses delivered by the implantable microstimulator.Based on the stimulation parameters for hypoglossal nerve stimulation in [4], the maximum duty cycle is approximately 1 %.This is beneficial when considering the power consumption of a battery-operated wearable to power the microstimulator and the energy absorption of human tissue due to changing electromagnetic fields, resulting in tissue heating [60].
The preliminary proof of concept investigations with the physics-based model show that it is in principle feasible to use ferroelectric materials in the dielectrics of nonlinear MLCCs as smart materials to enable passive control of the stimulation current.This physics-based model is the foundation for designing the passive control and can be further used to determine the nonlinear hysteresis shape to keep the stimulation current more constant and over a wider range of coupling factors.Furthermore, it offers the possibility to determine the material parameters of the dielectric corresponding to the exact structure of the MLCCs on the one hand and to establish the relationship between the nonlinear shape of the hysteresis and the material composition of the dielectric on the other hand.This opens the door to deeper investigations, for example in materials science, to design MLCCs with ideal nonlinear hysteresis.

V. CONCLUSION
Our results show that the physics-based model is well suited to predict the system behavior of inductively coupled microstimulators using highly nonlinear MLCCs.To model MLCCs with lower nonlinearities, it is necessary to extend the current model by improving the accuracy of the hysteresis inner-loop dynamics.In addition, our physics-based model can be coupled to other system models with reasonable effort and is generally applicable to any system powered wirelessly via resonant inductive coupling.However, particular attention should be paid to PTE.In addition, we demonstrate in this article that it is in Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
principle possible to use ferroelectric materials in MLCCs as smart materials to enable passive control of the stimulation current.We aim to determine a nonlinear hysteresis to compensate for application-related variations in the inductive coupling factor over a wider range, taking into account a fluctuation in the power level of the portable transmitter.
Beyond the scope of microstimulators, the proposed passive power control strategy should be a further step towards frugal engineering [61] leading to sustainable, affordable, and resource-efficient IPT systems using ferroelectric materials in MLCCs as smart materials.

Manuscript received 5
September 2023; revised 13 November 2023 and 20 December 2023; accepted 22 January 2024.Date of publication 1 February 2024; date of current version 19 April 2024.This work was supported by the Fraunhofer Center for Sensor Intelligence (ZSI).Recommended for publication by Associate Editor K.-H.Chen.(Corresponding author: Yves Olsommer.)

Fig. 1 .
Fig.1.Microstimulator circuit schematic (upper part) and measurement setup to characterize inductively coupled microstimulator circuits (lower part).The primary-side transmitter coil (L 1 and R 1 ) and the secondary-side microstimulator circuit coil (L 2 and R 2 ) were fixed on a spacer (a) with a constant coupling factor of 5.6 %.A detailed photo of the microstimulator circuit and MLCC #1 is shown on the right.For in vitro measurements, the microstimulator circuit coil and the transmitter coil were placed on the top and bottom of the pork flesh, respectively, in (b).The thickness of the flesh in (b) is about 5 mm.When both coils are positioned exactly on top of each other, the maximum coupling factor is 13 % and drops to around 5 % when the secondary-side coil is shifted about 2 cm.

Fig. 2 .
Fig. 2. Oscilloscope capture of: (a) induced voltage at the parallel resonant circuit in the microstimulator circuit (magenta) and with the half-wave rectifier generated stimulation pulse (green); (b) zoom in the induced voltage and stimulation pulse; (c) output voltage of the primary-side half-bridge (yellow) and the resulting current flowing through the primary-side series resonant circuit (blue).

Fig. 5 .
Fig. 5. Measured stimulation current I Stim as a function of the power P 1 fed into the transmitter (black) and calculations performed using the differential capacitance model (green) and the physics-based model (red) for MLCCs #1 through #8.Measurements for MLCCs #1 through #8 were performed at a frequency of 366, 382, 380, 360, 380, 373, 452, and 463 kHz, respectively.

TABLE I SELECTED
MLCCS WITH THEIR DESIGNATION, PART NUMBER, CAPACITANCE, RATED VOLTAGE, SIZE CODE, TOLERANCE, AND TEMPERATURE COEFFICIENT

TABLE II PARAMETERS
OF THE "EQUIVALENT" MATERIAL MODEL.THE PARAMETERS τ x AND V l WERE SET TO 10 NS AND 5•10 -23 M 3 , RESPECTIVELY

TABLE IV DIFFERENTIAL
CAPACITANCE C AT 0 V AND 90 % OF THE MAXIMUM APPLIED VOLTAGE U 90% IN A GIVEN VOLTAGE RANGE U P 1 increases from 67.6 W rms to 124.8 W rms , in addition to the stimulation current, the PTE increases from 0.9 % to 1.8 %.