Large-Signal Equivalent-Circuit Model of Asymmetric Electrostatic Transducers

This article presents a circuit model that is able to capture the full nonlinear behavior of an asymmetric electrostatic transducer whose dynamics are governed by a single degree of freedom. Effects such as stress-stiffening and pull-in are accounted for. The simulation of a displacement-dependent capacitor and a nonlinear spring is accomplished with arbitrary behavioral sources, which are a standard component of circuit simulators. As an application example, the parameters of the model were fitted to emulate the behavior of an electrostatic MEMS loudspeaker whose finite-element (FEM) simulations and acoustic characterisation where already reported in the literature. The obtained waveforms show good agreement with the amplitude and distortion that was reported both in the transient FEM simulations and in the experimental measurements. This model is also used to predict the performance of this device as a microphone, coupling it to a two-stage charge amplifier. Additional complex behaviors can be introduced to this network model if it is required.

devices. Kaiser et al. [1], for instance, report the fabrication of an all-silicon, CMOS-compatible electrostatic microloudspeaker based on the nano e-drive (NED) actuation principle [2]. Electrostatic MEMS microphones have already found widespread use [26]. Grixti et al. [12] show the design of a diaphragmbased electrostatic MEMS microphone whose sensitivity has negligible temperature dependency (unlike traditional electret condenser microphones). Applications of this principle also extend to the ultrasound range, where capacitive micromachined transducers (CMUTs) have emerged as an interesting alternative to piezoelectric transducers [13], [14]. These devices intrinsically couple three physical domains-electrical, mechanical, and acoustical-of which the electrical one exhibits an inherent non-linearity due to the inverse-square law present in Coulomb forces. More sources of nonlinearity arise if one considers stress-stiffening effects [3], [5], [19] or, if relevant, squeeze-film air damping [20]. Given the expense of simulating such multidomain, nonlinear systems with finite-element (FEM) solvers, lumped-parameter models capable of predicting their main performance variables become a valuable analysis tool. Moreover, if such a lumped-parameter model can be represented as an electrical circuit, not only the capabilities of a circuit simulator can be exploited, but also the interaction between the MEMS component and the circuits required for its operation can be seamlessly predicted. Such "multidomain" electrical circuits are a standard analysis tool in electroacoustics [16], [17], but are often limited to a small-signal analysis.
In this article, we present a method to expand the classical small-signal representation of an electromechano-acoustical network towards a large-signal, nonlinear network capable of reproducing effects like signal distortion and instabilities (particularly the catastrophic collapse of electrostatic transducers, commonly known as pull-in [21]). This circuit representation is done in such a way that no code-specific "derivative" or "integral" elements are used, enabling thus its implementation in a standard circuit simulator-it only requires using controlled sources that compute algebraic equations. We also apply this network model to simulate the behavior of the MEMS loudspeaker designed by Kaiser et al. [1], which has been shown to exhibit nonlinear stiffening effects [3]- [5], and also make a prediction of its behavior as a microphone. Although laid out for the case of electrostatic transducers, the procedure that we follow here can also be applied to build nonlinear network models of other systems, such as piezoelectric or electrodynamic transducers.
This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ Fig. 1. Scheme of a linear electro-mechano-acoustical network. The loops for current and force are coupled by an electrostatic transducer characterized by T EM . Likewise, the loops for force and volume velocity interact through a piston whose area is S eff .

II. SMALL-SIGNAL MODELLING
Representing a loudspeaker or a microphone as a linear, electromechano-acoustical network is a classical textbook exercise [14]- [16]. The underlying concept by which a mechanical system can be represented by a circuit is that the differential equations stemming from Newton's laws are strikingly similar to those of RLC circuits, so viscous damping elements can be understood as resistors, lumped masses as capacitors, and ideal springs as inductances (this holds for the mobility analogy, where forces are represented as "currents" and velocities as "voltages"). This circuit analogy can also be applied to acoustic systems, representing volume velocities as "currents" and pressure changes as "voltage" differences. Consequently, slits can be modeled by resistors, closed cavities as capacitors, and air masses as inductances. These three domains (mechanical, electrical and acoustical) interact with one another via transducers, which are represented as four-port elements. As Lenk et al. [15] explain, electrostatic transducers can be described to a first approximation by two cross-coupled current sources: the motion of a charged electrode generates a current, and the application of a voltage across electrodes provokes a force. This same argument can be used to model an ideal piston as a pair of coupled sources-the motion of the piston translates into a displaced volume and application of pressure is countered by a punctual force. A typical linear network for a capacitive electroacoustic system is depicted in Fig. 1.
This model is very well suited for a quick frequency analysis, particularly if its operation is limited to small oscillations around the working point. Even in this linearized case, a circuit solver becomes very practical given the amount of impedances required to reproduce the acoustic behavior-listing all the differential equations and solving them numerically is not only impractical but may easily lead to errors. However, the real input of an electrostatic transducer is not merely an ac signal: it necessarily includes a dc voltage. Furthermore, if the dc voltage deviates from the one for which the network parameters were calculated, the model becomes invalid (unless a large-signal equation is used to adjust them correspondingly). This model also limits the simulation of any electronic components coupled to the transducer (e.g., amplifier or driving unit) to their small-signal simplification; hence, neither power consumption nor distortions can be predicted with it.
The main difficulty with standard circuit simulators is that the nonlinearity in electrostatic systems (at first sight) is not an arbitrary I-V curve-as is the case for semiconductor devices [30]-, but it entails a non-proportional relationship between charge and voltage and between force ("current") and displacement ("magnetic flux"). In other words, it requires the simulation of a nonlinear capacitor and a nonlinear inductor. Some researchers have resorted to specific simulation tools that offer the possibility to write a code for these complex behaviors, such as the equation-defined devices in Qucs [8] or the numerical integration (and derivation) incorporated in Verilog-A [7]. The work of Tsuchiya et al. [9] includes a circuit model fully implemented in SPICE. Their approach was to approximate the nonlinear capacitance with polynomial current sources (thereby using Taylor expansions) and introduce numerically negligible "dummy springs" from which the displacement could be calculated. A more general solution to simulate this inherent nonlinear capacitance was introduced by Pham and Nathan [10], who describe how to simulate a capacitor (or a resistor) that depends on any of the circuit variables, "x," according to ). We present here an application of their methodology and expand it to include the stress-stiffening effects (i.e., the nonlinear inductor). Nonetheless, we also introduce an alternative approach that does not require to express nonlinearities in this form. The calculation of displacement and charge is done as an integral part of the circuit, so dummy elements are not required.

III. LARGE-SIGNAL MODELLING
From the mechanical perspective, an electrostatic transducer is a membrane or a beam governed by partial differential equations and driven into oscillation by a Coulomb force [27]. From the electrical perspective, the device is merely a two-port element that exhibits a strong capacitive impedance [11]. From the acoustical perspective, said membrane or beam behaves as a piston [15]- [17]. The goal of a large-signal equivalent-circuit is to describe a device with two electrical ports that, upon excitation with an ac/dc signal, predicts the deflection of the mechanical component and interacts with the acoustic circuit.
A key assumption for modelling the mechanical behavior is that the distributed motion of the elastic element can be accurately described by a finite number of discrete variables [27], [28]. Furthermore, if (within the operation range) the elastic motion is governed by a single degree of freedom, a very convenient ansatz can be applied, namely w( r, t) ≈ a(t)φ( r), which reduces the task to solving an ordinary differential equation for the amplitude, a(t). This assumption is evident for the imaginary case of a rigid plate attached to a punctual spring (as illustrated in Fig. 2). Interestingly, Melnikov et al. show experimental and numerical evidence that this is also the case for a clamped-clamped, prismatic beam [29]. We limit the theoretical development in this article to mechanical systems governed by Diagram of an asymmetric (single-sided) electrostatic transducer coupled to a piston, forming thus a capacitive loudspeaker (or microphone). In between may lie a mechanical amplification (such as the NED effect explained in detail in [2]). a single degree of freedom, although our formalism naturally extends to multiple degrees of freedom, resulting e.g., from a Galerkin-type procedure.

A. Description of a Spring-Capacitor System
Asymmetric ("single-sided") electrostatic transducers are often modelled as a spring-capacitor system [15], [16], whose electrical and mechanical constitutive equations are presented in (1) and (2). A parallel-plate capacitor (C 0 ) has an initial gap (g 0 ) that is reduced when its moving electrode is displaced (x) towards the fixed one. Such a displacement is the result of applying a voltage (U in ) across the electrodes, but is opposed by the force of a polynomial spring (k i )-we introduce the higher-order terms with the aim to model stress-stiffening effects [3], [19]. Inertial (m) and damping (r m ) effects are also accounted for. The F out term captures any additional loads attached to the system. The system described by (2) is a Duffing oscillator 1 driven by a Coulomb force The key to implementing the variable capacitor and the Duffing spring in a standard circuit simulator is to use a suitable change of variables. This is done by incorporating two additional loops: one in which charge in the capacitor is represented as a current, and one in which the electrode's displacement is represented as a voltage (see Fig. 3) With this change of variables, a mere resistor would model an ideal capacitor or a linear spring, respectively. Following this logic, (1) and (2) can be directly implemented by introducing arbitrary behavioral current sources. The "current" source in the charge loop describes the displacement-dependent capacitance, whereas the "current" sources in the displacement loop define the force balance. In order to couple these two additional loops with  . Alternative approach to simulate the spring-capacitor system using the methodology proposed by Pham and Nathan [10]. The nonlinear capacitor is constructed with a series voltage source and the nonlinear inductor with a parallel current source.
the other domains, (where actual current is to be interpreted as current, and velocity as voltage) a numerical derivative is necessary. This is achieved by introducing a unitary 2 capacitor and a unitary inductor in two separate loops. The reader may recognise by inspecting the loops in Fig. 3 that this large-signal model is merely a graphical representation of (1) and (2).
It was mentioned that Pham and Nathan [10] had already proposed a methodology to simulate a nonlinear capacitor whose expression is of the form C (x) = C 0 /(1 + g(x)); we extended this principle to simulate a spring ("inductor") with Duffing coefficients and therefore also simulate the Duffing-Coulomb oscillator-the result is the circuit of Fig. 4. The displacementdependent capacitor is constructed with a current-dependent voltage source in series with a linear capacitor. The reader may verify that U in = u 0 (1 − x/g 0 ), and therefore the charge entering the capacitor follows (1). Notice that the displacement is read out from the force ("current") of the linear component of the spring (k 1 ). Likewise, the displacement-dependent spring is constructed with a current-dependent current source in parallel with a linear inductor. The total force ("current") on the spring , and the electrostatic force is expressed in terms of the charge F C = q 2 /(2C 0 g 0 ), which completes the force balance (2). In comparison to the circuit from Fig. 3, the approach of Pham and Nathan [10] is more compact in its representation, but it requires nonlinearities to have a specific form. In contrast, in the circuit introduced in Fig. 3, the equations for force and charge are directly written into the arbitrary behavioral sources. The advantage of having separate loops for displacement and current is that the option for calculating further effects remains open (e.g., saturation of the displacement after pull-in), whereas the advantage of having a single current loop is that spurious current impulses coming from the numerical derivative of the charge are avoided. Fortunately, these two solution approaches are not mutually exclusive: the electrical side can be modelled in the way of Fig. 4, while the mechanical domain as in Fig. 3, achieving both a smooth current signal and a robust displacement calculation. Both circuit models are suitable for dc, ac, and distortion analysis.

B. Generalisation of the Capacitance and Electrostatic Force
Even though the spring-capacitor system represents a good approximation for electrostatic transducers, a large-signal network model need not be based on this model, provided that an appropriate expression relating capacitance and deflection is found. Standard MEMS microphones [12], as well as CMUTs [13], are based on a membrane that bends upon application of an electric field, so the equation of a parallel-plate capacitor imposes certain limitations. A more recent development, the pull-in-free microphone from Ozdogan et al. [25], is based electrostatic levitation and holds no resemblance to the springcapacitor system. This transducer was therefore modelled with a displacement-varying capacitance following a polynomial expression. A more general implementation of (1) and (2) would thus be whereby f (x) portrays the variation of capacitance against a reference displacement. For instance, Schenk et al. [6] have recently demonstrated that the capacitance of a Coulomb-actuated, clamped-clamped, prismatic microbeam follows (3) and (4) with the σ-constants being σ 1 = ψ 2 77 , σ 2 = − ψ 2 38 , σ 3 = 15ψ 2 28 with ψ ≈ 1.5881. This general case can also be implemented as a large-signal equivalent circuit following the previous approaches. We show in Fig. 5    u 0 = U in + u C , so the charge in the capacitor follows (3). The displacement loop implements the force balance (4) directly.

C. Overall Electromechano-Acoustical Network
Having derived an equivalent-circuit model for the asymmetric electrostatic transducer, the complete electromechanoacoustical network would have the structure shown in Fig. 6. Notice that, unlike the linear network of Fig. 1, this model is much more robust: it allows simulating the complete input, which has both a dc and ac component, and the parameters of the nonlinear electromechanical network are independent of the dc level. Moreover, the system can be simulated in both transduction directions, that is, either as a loudspeaker or as a microphone.

IV. APPLICATION TO A MEMS TRANSDUCER
We have configured the circuit model of Fig. 5 to describe the behavior of the microspeaker designed by Kaiser et al. [1], a picture of which is shown in Fig. 7. The building block of this microspeaker is a clamped-clamped microbeam composed of  [1], showing two pairs of actuators (the beams are only shown until the midpoint). These beams deflect by virtue of the NED effect [2]. Given its asymmetric (single-sided) nature, it is termed an ANED transducer.
three arc-shaped electrodes-the outer electrodes are connected to the reference voltage, and the middle one is connected to the input signal. The attraction between electrodes is translated into an aggregate deflection of the beam by virtue of the nano e-drive effect demonstrated in [2]. One can surmise that this attraction between electrodes behaves similarly to a spring-capacitor system; however, the variable of interest in this case is the aggregate deflection of the beam, not merely the displacement of the electrodes. Spitz et al. [3], [4] proposed to couple these variables by introducing a lever factor, i.e., w = θx (forces are correspondingly scaled to preserve energy). Such a lever effect is schematically presented in Fig. 2. Even though the spring-capacitor system may provide a good approximation, the expressions from Schenk et al. [6] describe a system that bears a closer resemblance to the asymmetric NED ("ANED") beam: a prismatic, clamped-clamped microbeam. For this reason, (3) and (4), with the functions in (5) and (6), were used as a starting point to build the lumped-parameter model of this actuator. Given that the expressions from [6] stem from a solution of the Euler-Bernoulli equation for a Coulomb force, we term this model a "Coulomb-Bernoulli fit." With a model that describes the motion of the beam, it becomes possible to calculate the volume of air that it displaces as it oscillates, and with it the amplitude of the pressure waves it generates. This displaced volume is directly proportional to the deflection at the midpoint of the beam, the ratio between them being the effective piston area S eff . Twenty-eight beams connected in parallel (electrically, mechanically, and acoustically speaking) constitute the complete loudspeaker.

A. Parametrisation With FEM Simulations
In order to find the values of the parameters in (4) that best describe the ANED microbeam, a least-squares fit was applied to the FEM simulations of Melnikov et al. [5]. The values of the zero-volt capacitance and electrostatic gap are obtained directly from the simulation. The values of the elastic constants (k 1 , k 2 , k 3 ) and the NED lever factor (θ) are then adjusted to minimise the error with respect to the bifurcation diagram reported in [5] (see Fig. 8.)-more specifically, with respect to its stable branch. Notice that the implementation of the arc-length solver enabled the FEM simulation to accurately calculate the pull-in voltage (61 V) as the transition of the stable branch towards  (3) and (4), using the expression of capacitance suggested in [25]. An optical measurement of a fabricated microbeam is superimposed, as well as the experimental pull-in voltage. Fig. 9. Comparison of the capacitance of the ANED microbeam reported in [5] with the two model that was presented in Fig. 7. An indirect, experimental calculation of the beam capacitance is superimposed to the curves. the unstable one. It is worth mentioning that nonlinear effects arising from large deflections were already accounted for in the FEM model ("NLGEOM,ON") [5]. By fitting the three elastic constants to the stable branch of the curve, the deflection is predicted with high accuracy, although the critical voltage is slightly overestimated (62 V). The Coulomb-Bernoulli fit also accounts very well for the variation of the capacitance with the DC voltage (see Fig. 9), although the values close to the critical point are slightly underestimated. A high-fidelity parameter fit is essential to reproduce the simulation values-especially if the harmonic distortion is to be predicted-, and these values should in turn be distributed along the full stable branch. This requires not only performing a fine-grained voltage sweep (in this case with 242 points ranging from zero deflection until pull-in), but also choosing a suitable parametrised equation, which is not necessarily unique. Along the region of operation (30 -50 V), the prediction of the Coulomb-Bernoulli fit has a deviation within ±0.2%, which is a good indication that it can reproduce the results of the simulation.
In addition to the simulated values, Figs. 8 and 9 contain also the experimental measurements of the quasi-static deflection of the beam and its capacitance. A detailed description of the set-up is reported in [3], [4]. A digital-holographic microscope (DHMR2203, LyncéeTec.) was used to track the position of the midpoint, applying an AC signal of low frequency (50 Hz). Regarding the capacitance, an impedance analyser (Keysight E4990A) was used to determine the effective capacitance with a variable bias and an AC signal of 1V at 1 kHz. (An additional optical measurement was performed to verify that the actuator was still behaving in the quasi-static regime.) Given that the capacitance measurement also includes parasitic fields, the difference between the total capacitance and the simulated zero-volt capacitance of the actuators was subtracted from the measurement, and this value was divided by the number of connected beams. It was assumed, thus, that the only changes of capacitance with voltage came from the microbeams. The experimental measurement of the deflection fits very well with the simulated curve; nonetheless, the measured value of the critical voltage (55 V) is lower than the simulated one. This might be attributed to differences in the actual dimensions of the beam as compared to the simulated ones, whose impact on the stiffness is more notorious at high voltages. The measurement of capacitance does manifest an increase with the bias voltage, although the projection of the increase of capacitance within an individual beam is higher than in the simulation. Given that this is not a direct measurement, it should only be taken as a reference-the estimated parasitic capacitance per actuator (∼4 pF) is roughly four times higher than the active one.
The previous procedure allowed estimating the static parameters, but the dynamic ones (effective mass and damping coefficient) require further information. Melnikov [5] report that the frequency of the first mode of oscillation of this beam is at 9.28 kHz, followed by a second-mode oscillation at 31.9 kHz. The wide separation between them is a good indication that a single-mode approximation is justifiable-the second mode should not be excited in the application range [29]. The effective mass is thus computed as m = (2πf 1 ) 2 /k 1 . Due to the difference between the simulated (9.28 kHz [5]) and measured (8.5 kHz [4]) resonance frequency, the value of the effective mass was adjusted correspondingly when comparing the model with the experiments reported in [1]. The damping coefficient remains unfortunately as a free parameter: it represents several energy loss mechanisms that can affect the structure and can be very challenging to predict [18]. A common practice is to make an "order-of-magnitude" estimation of the quality factor Q m = √ k 1 m /r m . For the sake of comparison, a value of Q m = 5 (the same reported in [5]) was used when comparing the model to the transient FEM simulations, whilst a value of Q m = 2.6 (fitted by trial and error) was used when comparing to the experiments.
An interesting feature of the ANED microbeam is that its resonance frequency increases with the applied bias-an effect that has been both measured and simulated [3]- [5]. A classical spring-capacitor system cannot reproduce this effect, since the   [1] Coulomb force will always have a softening effect. As it can be seen in Fig. 10, the proposed lumped-parameter model is able to partially predict this increase in the resonance frequency, which is a result of the addition of Duffing terms to the stiffness. However, this increase could not be sustained after 30 V (half of the voltage range), which hints that further corrections should be done to the (3) and (4) in order to better capture the dynamic behavior at voltages close to the critical point.
The obtained network parameters of this asymmetric transducer are given in Table I. As explained before, the effective mass and mechanical quality factor were adjusted when comparing against experimental data.
The abstract calculation of an effective lever factor (θ) representing the NED amplification can be verified by further results of the quasi-static FEM simulations. The "electrode displacement" (x)-in the picture of a parallel-plate capacitor-can be understood as the average of the reduction in the electrostatic gaps upon application of a voltage (notice that the beam has three electrodes and thus two gaps). When the midpoint deflection is plotted against this average gap reduction (see Fig. 11) a direct dependency is observed, although not entirely linear. There appears to be a slight decrease in the ratio between midpoint displacement and gap reduction resembling a quadratic Fig. 11. Variation of the midpoint deflection with respect to the average gap reduction. The ratio between them corresponds to the mechanical amplification achieved by the NED mechanism [2]. behavior. It appears that the estimation of θ = 11.9 lies above the actual amplification ratio (∼6.6), so a further refinement of the governing equations should be explored in later studies.
A first test of the obtained parameter set can be performed against the transient FEM simulations from Melnikov et al. [5]. In this simulation, the ANED beam was subject to a sinusoidal excitation of 10 V in amplitude at a frequency of 1 kHz, superimposed to a bias of 30 V. Applying the same excitation to the network model resulted in a deflection that matches very well with the FEM simulation (see Fig. 12). A closer analysis of the harmonic components, as given in Table II, reveals that the constant level, as well as the first and second harmonics are almost exactly reproduced. The difference between simulations starts becoming apparent only from the third harmonic. Calculation of the total harmonic distortion (THD)-using the definition of the IEC-yields a very similar result for both cases.
The last element in the parametrization of the network model is the acoustic impedances and the effective piston area. The impedance of a rectangular slit includes a resistive and an   [1]. The input signal is set to 40 V dc plus 5 V ac. The acoustic load is the standardized network model of the IEC 60318-4 ear simulator [23].
inductive component that can be directly determined from its dimensions [15], [17]. These expressions were used to calculate the impedance of the front-and back-side slits. The clearance between the actuator and the bottom and handle wafers introduces a leakage that was also modelled as a rectangular slit. Finally, the effective piston area was directly calculated from the FEM simulations as the ratio between displaced volume and midpoint deflection. These values are given in Table III.

B. Simulation as a Loudspeaker
The complete network model of the microspeaker, as illustrated in Fig. 13, was implemented in SPICE and compared against the measurements of Kaiser et al. [1]. In this experiment, the dc voltage was set to 40 V, and the ac amplitude at 5 V. The chip was coupled to an IEC 60318-4 ear simulator [23], which closely reproduces the impedance of the ear canal. Figs. 14 and 15 show a comparison of the predicted and measured sound pressure level (SPL)-relative to 20 μPa rms-and the THD. Both the predicted and measured sound pressure levels are  [1] and simulated sound pressure level of the ANED loudspeaker reported in [1]. Save for the underestimated damping around the resonance peak of the ear simulator (around 11 kHz), a good agreement is observed between the network model and the experiment. Fig. 15. Measured [1] and simulated THD of the ANED loudspeaker reported in [1]. The critical peak around 5 kHz is overestimated by the simulation (an effect related to the resonance of the ear simulator), but the overall behavior is well represented by the network model. in good agreement. The most notorious difference is that the resonance peak of the ear simulator (at 11.2 kHz) is much higher in the simulation than in the measurement, which indicates that the damping in the circuit model of the IEC 60318-4 coupler differs from the damping obtained in the measurement. This circuit model is given by the fabricant [23], so an adjustment with the experimental measurements cannot be undertaken.
The prediction of the harmonic distortion also matches well with the experiments, although the critical peak was overestimated by the simulation. Both the low frequency and high frequency distortion are accurately predicted by the model, as well as the bulge around 1 kHz. The constant increase of the distortion after 1 kHz and its sudden drop after 5.6 kHz are also very well reproduced, but the peak right at 5.6 kHz seems to be exaggerated by the simulation. This is not a surprise if one considers that 5.5 kHz is exactly the half of the resonance brought about by the ear simulator: The quadratic (Coulomb) force inevitably generates a harmonic at twice the frequency, which is overly amplified when it coincides with the resonance of the ear simulator. This region between 4-6 kHz is expected to have the highest distortion because the mentioned second harmonic excites both the acoustical and mechanical resonances, which are near each other. A further effect, observed in the simulation but not in the measurement, is the very narrow peak in distortion at 3.7 kHz. This frequency is, again, an integer division of the resonance of the ear simulator (one third), which also indicates that the third harmonic was amplified by the ear simulator at this specific frequency. The absence of this peak in the measurement is also consistent with the mentioned difference in the damping of the actual and simulated acoustic load.

C. Simulation as a Microphone
The same network model of the transducer used for the simulation as a loudspeaker can be used to predict its behavior as a microphone, changing only the circuitry used for its operation and introducing pressure signal in the acoustical domain (see Fig. 16). A common implementation of a MEMS microphone includes the stiffening effect of a back volume [26], so an additional capacitance was introduced (N back ), emulating a cavity of 2 cm³. In the electrical domain, a two-stage amplifier was implemented to measure the charge in the transducer. The first stage makes an upstream measurement of the generated current with an instrumentation amplifier (AD8428), and the second stage integrates the signal with an operational amplifier (AD795). One further advantage of having a large-signal equivalent circuit model is evidenced in this simulation: the tested models of commercial electronic components can be simulated together with the MEMS transducer. Fig. 17 shows a linear frequency analysis of the microphone circuit, reporting the ratio between output voltage and input pressure (electrical sensitivity) and the ratio between midpoint deflection and input pressure (mechanical sensitivity). A peak is observed in the same frequency region (∼ 9 kHz) for both curves, corresponding to the resonance frequency. The sensitivity of this loudspeaker acting as a microphone is very low: roughly 5 mV are obtained from an input pressure of 1 Pa, even after considerable amplification. Further simulations in the time domain (not linearized), as given in Table IV, are consistent with the Fig. 17. Simulation of the behavior of the ANED microspeaker acting as a microphone. The mechanical sensitivity stands for the ratio between the amplitudes of the deflection and the input pressure, whereas the electrical sensitivity relates the (amplified) output voltage to the input pressure. frequency sweep. With these transient simulations, the harmonic distortion both before and after amplification can be calculated. One can observe how the amplifiers introduced a certain degree of distortion, since the THD of the measured current was always lower than that of the output voltage.

V. EXTENSIONS TO THE MODEL
The introduced network model allows adjusting the expressions in (1) and (2) to better account for more nonlinear effects, if required. The collapse of the electrodes after pull-in can be incorporated by placing a saturation block in the displacement loop: if the "voltage" (displacement) exceeds a certain value (e.g., g 0 ), it should not increase indefinitely but stay at this value until the forces are correspondingly reduced. The general expression of the electrostatic force in (3) and (4) can be applied to other cases like that of electrostatic levitation [25], which can be modelled as a sevent-degree polynomial. Although not explored in this article, nonlinear damping effects can be also taken into account in the network model, since they only require the simulation of non-proportional I-V curves (damping was modelled here as a mere resistive behavior).
It should also be said that asymmetric electrostatic transducers are not the only circuit topology used in commercial devices: their balanced variant is also of interest because of its higher linearity [16]. The analytic procedure presented in this article can be extended to describe such transducers by introducing an additional variable capacitance and a corresponding force pulling the spring in the opposite direction.
If the single-mode approximation is insufficient for describing the mechanical behavior, the presented network model could be expanded to include multiple oscillators. Each oscillation mode would also have an equivalent piston area, and all these "pistons" would be connected in parallel in the acoustic circuit, allowing the possibility for certain modes to subtract some of the volume velocity pumped by others.
This approach is not limited to electrostatic transducers: electrodynamic loudspeakers-which rely on an inductive transduction instead of a capacitive one-may also be modelled with such an equivalent circuit. It has been reported that these loudspeakers also exhibit a Duffing-like suspension stiffness, as well as a quadratic, displacement-dependent transduction coefficient [24]. The transduction analogy in this case, nonetheless, would be that of a "voltage"-dependent (velocity-dependent) voltage source in the electrical domain, coupled with a currentdependent "current" source (force source) in the mechanical domain.

VI. CONCLUSION
The constitutive lumped-parameter differential equations of an asymmetric electrostatic transducer can be equivalently described as an electric circuit with multiple, cross-coupled loops representing each physical domain. Nonlinearities in this circuit description require only the implementation of arbitrary behavioral sources-the same strategy by which semiconductor devices are simulated [30]. The hurdle of simulating a displacement-dependent capacitance and a force-dependent spring was overcome with two alternative approaches: either by introducing separate loops for the electric charge and mechanical displacement (our solution method), or by calculating corrections for the voltage in the capacitor and the force in the spring (a method proposed by Pham and Nathan [9]). The proposed large-signal equivalent circuit model was used to reproduce the behavior of an electrostatic MEMS loudspeaker reported in the literature [1], using the outcome of the FEM simulations reported in [5]. The introduction of nonlinearities in the elasticity as well as in the electrostatic force helped to accurately describe the variation of deflection and capacitance with the input voltage. As a result, both the sound pressure level and the distortion of the loudspeaker could be calculated in a standard circuit simulator, obtaining a good agreement with the measurements in [1]. In addition, the circuit model of this transducer was used to predict its behavior as a microphone, coupling it to an amplifying stage. It was thus illustrated that equivalent-circuit models need not be limited to a small-signal analysis and offer the opportunity to simulate an electrostatic transducer together with complex networks (such as an ear simulator or an amplifier stage).