A Normalized Model of a Microelectromechanical Relay Calibrated by Laser-Doppler Vibrometry

This work presents a behavioral model for a microelectromechanical (MEM) relay for use in circuit simulation. Models require calibration, and other published relay models require over a dozen parameters for calibration, many of which are difficult to extract or are only available after finite element analysis. This model improves on prior work by taking advantage of model normalization, which often results in models that require fewer parameters than un-normalized models. This model only needs three parameters extracted from experiment and one dimension known from device fabrication to represent its non-contact behavior, and two additional extracted parameters to represent its behavior when in contact. The extracted parameters–quality factor, resonant frequency, and the pull-in voltage–can be found using laser Doppler vibrometry. The device dimension is the actuation gap size, which comes from process data. To demonstrate this extraction process, a series of velocity step responses were excited in MEM relays, the measured velocity responses were used to calibrate the model, and then then simulations of the model (implemented in Verilog-A) were compared against the measured data. The error in the simulated oscillation frequency and peak velocity, two values selected as figures of merit, is less than 10% across many operating voltages.


Introduction
Vertical-gap microelectromechanical (MEM) relays have wide ranging applications in circuit design, especially in radio frequency applications [1]. For instance, MEM relays have been used as transmit/receive switches and antenna tuners in radios [2]. Recent work has also explored the application of MEM relays to digital circuit design, examining their potential as logic gates [3], power gates [4], and back-end-of-line integrated switches [5,6].
Circuit designers rely on behavioral models of MEM relays, often written in SPICE or Verilog-A, to design relays using electronic design automation software [3,[7][8][9][10], but these models often require the user to specify many material and device parameters, which can limit their use. For instance, the model in [3] requires two material properties, nine layout dimensions, two vertical device fabrication dimensions, and three values that are extracted from measurements. As another example, the model in [8] requires six material parameters, two layout dimensions three vertical device fabrication dimensions, and one parameter extracted from finite element analysis (FEA). Models that have many parameters or difficult-to-extract parameters, particularly if the parameters require FEA or device measurements, are more difficult to develop.
Model normalization simplifies relay models, resulting in fewer, easier-to-extract parameters. Normalizing models is a technique that was originally developed for circuit simulation to reduce rounding errors [11]. The technique requires dividing a behavioral model by an upper bound to change the scale of the numbers representing the model's state [12]. Properly selected normalization terms can also reduce the number of parameters in a model [13]. No normalized relay models have been reported to date.
Normalized models still have parameters, and this work suggests that laser Doppler vibrometry (LDV) is a good tool to measure the model parameters the normalized model in this work. This is because typical LDV characterization captures the dynamic behavior of a device from its step or harmonic response, making it well-suited to directly measure quality factor and natural frequency. This is not a new observation because LDV is a proven technique for extracting models from microdevices [14,15], and LDV has also been used for switch models in RF switches [16,17]. However, the RF switch models in [16,17] rely on FEA models in conjunction with LDV. Combining LDV and normalized analytical models offers the prospect of directly measuring relay model parameters and eliminating the need for FEA.
This work contributes a normalized dynamic model of a vertical-gap MEM relay that uses a small number of easily extracted parameters in the behavioral model. The normalized model relies on three parameters that can be extracted from measurementsquality factor, Q 0 , natural frequency, ω 0 , and pull-in voltage, V pi -and one vertical device fabrication dimension-the actuation gap, g-to describe the relay's motion when the drain does not contact the source. The model can be extended to describe contact by adding two more extracted parameters and a few tuning parameters to affect simulator behavior. All parameters are calibrated using only simple LDV measurements, and calibrating them does not require FEA. This model was validated by implementing it in Verilog-A, simulating it using a standard circuit simulator (hSpice) and comparing the simulated velocity against measurements.
Further description of this model can be found in the remainder of this work. Section 2, describes the device being modeled, the normalized model, the experiments used to extract model parameters, and the algorithm for extracting model parameters from measured data. Section 3 shows the measured data, the extracted parameters, and the comparison of the simulated model to measurements. Section 4 concludes the paper and suggests future work.

Device Description
The relay analyzed in this work was first demonstrated in [18], and a micrograph and cross-section of a device are shown in Figure 1. The device is a parallel-plate MEM switch in which a movable shuttle, called the gate, is supported by folded flexures above an electrode on the substrate, called the body. The gate and the body are separated by a distance g = 220 nm, called the gap. . The cross-section shows that when a voltage difference is applied between the Gate and Body, the device is closed and current (I ds ) can flow. The actuation gap, g, is also labeled.
During operation, a Voltage V g is applied to the gate and a different Voltage V b is applied to the body. The difference between these Voltages, V gb = V g − V b , creates an electrostatic force that causes the gate to move toward the body. After the device has moved a distance g d = 22 nm, called the dimple gap, an electrode affixed to the gate, called the drain, makes contact with an electrode affixed to the substrate, called the source. A drain-source current, I ds , can flow when the drain and source are in contact. No current flows between the gate and the drain because they are separated by a non-conductive oxide. The Voltage at which the source adn the drain make contact is called the contact Voltage, V con .
The displacement of this parallel plate actuator from its resting position, x, is determined by a balance of electrostatic forces from V gb and spring forces from the supporting flexures. This electrostatic force balance results in the well-known pull-in instability: if x > g/3, the electrostatic force increases with decreasing gap size faster than the flexures' restoring force. This instability occurs at the pull-in voltage, V pi . The devices in this paper have g d < g/3, so the drain electrode prevents the gate from displacing to g/3. As a result, the relay always operates in non-pull-in mode. However, it is still possible to calculate V pi for a non-pull-in device -V pi is the Voltage that would be required to displace the gate by g/3 if drain electrode were not present-and doing so is important in the following analysis.
Note that using g d = 22 nm is a deviation from [18], which lists the as-designed g d as 60 nm. All of the measured V con values in this work, which are approximately 14 V, are consistent with g d = 22 nm. Changing g d affects the value of V con without changing V pi or other parameters, all of which appear correct because they were derived both experimentally and theoretically. Changing g d made the data consistent with the least possible change to published process parameters.

Free Space Model
The relay can be modeled as a second order dynamical system when it is not in contact with the drain and source electrodes, where m is the effective mass of the relay, b is the damping coefficient of the relay, k is the spring constant of the folded flexures, and x is the displacement of the gate. The right hand side of the equation is an electrostatic force, denoted F el in future equations, where 0 is the relative vacuum permittivity, assumed to be 8.85 pF/m, V gb is the voltage between the gate and the body, A is the active electrode overlap area between the gate and body, and g is the size of the gap between the gate and body. The expression for the dynamic system can be normalized by kg, to givë x wherex = x/g is the normalized gap size, ω 0 = √ k/m is the system's natural frequency, is the system's quality factor, and V pi = 8 27 kg 3 0 A is the pull-in Voltage of a parallel plate actuator. This model is modified in Section 2.5 to describe relay behavior in contact with the surface, and additional practical details required to implement this model are discussed in Appendix A.

Experiment
Calibrating and verifying the model in Section 2.2 required measurements of the MEM relay's dynamic behavior. Those measurements were obtained by configuring the relay as a pseudo-inverter, a common test structure [3], and manipulating the voltages on the device to induce a series of mechanical step responses. The velocity of the gate during these step responses was measured with LDV.
The details of the experimental configuration are as follows: the drain was connected to a 2 V Voltage source through a 47 kΩ pull-up resistor, the source was grounded, the body was connected to a bias Voltage that varied from 6.5 to 13.0 V in different trials, and the gate was driven by a ±1 V square wave. This resulted in small displacements of the device, and the velocity corresponding to those displacements was measured using LDV-specifically, a Polytec OFV-534 sensor head with an OFV-5000 controller and a VD-09 velocity decoder. The velocity decoder was set to a resolution of 200 mm/s/V, so it had a bandwidth of 2.5 MHz. A 10x lens on the OFV 534 sensor head was used to focus the laser spot onto the gate electrode to measure micromechanical movement. The laser spot was adjusted to achieve maximum focus on the device gate. A velocity response of the gate to a step input in V g and the displacement curve calculated by integrating the velocity are shown in Figure 2. . This is higher than the deviation observed in the measured velocity curve above, and may represent a worst case measurement. The LDV displacement resolution, estimated by integrating the velocity noise deviation over the sampling window of 50 ns, is 40 pm.
The Voltage ranges used in this experiment were selected deliberately. The swing value of V g was set to minimize stress on the gate oxide and to create small displacements of the gate, which allowed for trials to be carried out at many different V b biases. The minimum V b bias was chosen to make the gate move quickly enough to overcome the LDV noise floor. The maximum V b was limited to prevent the drain from contacting the source. Though the drain and source make contact in normal operation, avoiding contact here was desirable for two reasons: contact introduces a risk of device failure into an otherwise non-destructive test, and contact introduces non-linear forces that made it harder to extract model parameters from the measured velocity.

Extracting Model Parameters
Equation (2) is described by four parameters: the pull-in voltage, the natural frequency, the quality factor, and the gap size. Though the gap size does not appear directly in (2), it is necessary to normalize the x variable. The first three parameters can be extracted experimentally, while the gap size requires knowledge of the fabrication process.
In order to find V pi , two states of static equilibrium can be compared before and after the Voltage step on the gate. Several values are known in these static equilibria: before the gate moves, the system has a gate-body voltage called V gb,open , after the gate moves to its next position the gate-body voltage is called V gb,closed , and velocity and acceleration are both zero in both static equilibria. These known values were substituted into Equation (2), which resulted in the following system of equations in V pi ,x, and a small shift in the normalized displacement, ∆x:x These equations can be solved forx and V pi using measured ∆x values.
Extracting ω 0 and Q 0 from measurements is complicated by the non-linear behavior of the electrostatic force. If the relay is moving only a short distance while affected by an electrostatic force, the electrostatic force can be represented as a constant force and a nonlinear, electrostatic contribution to the spring constant, k es . The electrostatic contribution to the spring constant can be found by taking a derivative of the electrostatic force with respect to the gate displacement: Combining the mechanical spring constant with the electrostatic spring contribution given in Equation (5) gives an effective spring constant for a small displacement around the actuator's equilibrium point, given by This effective spring constant produces an effective quality factor, Q e f f and an effective resonant frequency, ω e f f . Expressions for Q e f f and ω e f f can be found by substituting (6) into the expressions for Q 0 and ω 0 : ω The expression 8 27 V gb V pi 1 (1−x) 3 is called the voltage-displacement term. The effective natural frequency and quality factor are extracted from measurements by observing the period of oscillation of the step response, T osc , and the amplitude of successive velocity peaks,ẋ max,i , where i is an index. Figure 2 shows these parts of the velocity curve. Standard second order differential modeling can be used to derive equations that relate T osc andẋ max,i to ω e f f and Q e f f : 2π whereẋ max,0 is the height of the first peak in the velocity measurement.

Contact Model
Careful examination of LDV measurements when the drain and source make contact, which can be seen in Figure 3, indicate that the free-space model does not describe this regime well. The free space model fails to accurately predict ω e f f and Q e f f when the drain and the source make contact. The techniques in Section 2.4 can be used to find ω e f f and Q e f f for the closing transient in Figure 3, and these extracted quantities are called ω con and Q con . Though Equations (7) and (8) indicate that ω e f f and Q e f f should decrease with high V gb andx, the extraction process indicates that ω con and Q con are higher than free-space ω e f f and Q e f f . This behavior can be captured in simulation by updating the model presented in Equation (2): Several terms have been added to this updated model: exp(p con (x − g d /g)) is a force, called F con in future equations, that represents contact with the surface, p con is a tunable convergence parameter that is discussed below, and ω(x) are Q(x) are a natural frequency and a quality factor that change as a function ofx.
There are multiple ways to implement F con . The simplest implementation is a piecewise function that is large whenx > g d , but piece-wise functions can cause convergence issues for simulators. Another implementation is using high order polynomials to accurately represent a three dimensional integral of Lennard-Jones potentials, and that implementation appears in [8]. Such a high order polynomial can include a third order term to represent an attractive Van Der Waals force. Instead of these other implementations, F con is implemented as an exponential here for three reasons: it is smooth, monotonic and continuously-differentiable; Verilog-A has a limexp function that provides a safe representation of exponentials to simulators; and the proportionality constant, p con , is a simple parameter to tune if the model is not converging. p con was set to a large value (500) to make F con approximate a step in force.
The contact model changes ω 0 and Q 0 to new values, ω(x) and Q(x), when the drain and source electrodes are in contact. These new values reflect changes in the spring constant, k, and the damping constant, b. These changes can be caused by a variety of effects: squeeze film damping (which may explain the "bump" in the closing transients in Figure 3 based on preliminary analysis from [20]), energy lost to impact with the surface, the increased deflection of the gate relative to the drain and source electrodes, etc. This model does not attempt to capture each of those phenomena in detail. Instead, the model is designed to mimic behavior extracted from the LDV by using a transition function to modify ω 0 and Q 0 asx becomes large: These functions use a hyperbolic tangent to approximate a smooth change in the value of ω 0 and Q 0 . The hyperbolic tangent is described by new convergence parameters: p ω and p Q describe how quickly ω 0 and Q 0 change asx changes, and f ω and f Q , which are fractions between 0 and 1, describe how muchx must displace before the ω 0 and Q 0 begin to change. p ω and p Q were set to 100, and f ω and f Q were set to 0.9. Figure 3. Displacement of the relay during transients that make contact with the surface. V B is higher during these measurements to allow the relay to make contact. The bump in the contacting behavior during the closing transient reflects additional displacement of the gate structure of the relay after the drain and source are in contact, and that displacement may be influenced by squeeze film damping. (left) Displacement when the relay is openining. (right) Displacement when the relay is closing and making contact.

Velocity and Displacement Data
Several trials of the experiment in Section 2.3 were conducted using different body voltages, and LDV produced plots of the transient velocity of the MEM relay for each trial. Integrating the velocity plots produced plots of displacement vs. time, which are pictured Figure 4. LDV does not produce a signal that is perfectly centered at 0 mm/s, so it was necessary to correct this artefact during integration. The curves in Figure 4 have a correction factor of c(t − t0) added to each point of the displacement curve, where t 0 is the time at the start of the measurement, t is time, and c is a correction factor that was extracted from the slope of the uncorrected position data.

Parameter Extraction
The first step in parameter extraction was calculating ω e f f , Q e f f and ∆x for each of the transient responses. ω e f f and Q e f f varied slightly from cycle to cycle because the amplitude of the oscillation decayed, which affected the quality of the small-displacement approximation that was used to calculate the electrostatic spring. The mean of the ω e f f and Q e f f values from the first four cycles was used for the analysis. ω con and Q con were calculated from transient responses where the drain and source made contact, but ω con and Q con are not included in calculations of ω 0 and Q 0 .
The second step of parameter extraction was finding V pi using the measured value for ∆x, the known values for V gb,open and V gb,closed , and Equations (3) and (4). In addition to a V pi estimate, solving these equations produced ax − V gb curve, pictured in Figure 5. The extractedx − V gb points are compared against an analytical estimate of the position that was found by using Equation (2) using process data. The extracted curves match analytical expectations closely, deviating only at high V gb values. At those values,x is relatively large, so fringing capacitance, which is not modeled, may be increasing the electrostatic force. The final step of parameter extraction was finding ω 0 and Q 0 . The squares of the extracted ω e f f and Q e f f were plotted against the voltage-displacement term from Equations (7) and (8). Those plots are pictured in Figure 5. Equations (7) and (8) show that ω 2 e f f and Q 2 e f f both have a linear relationship with the voltage-displacement term, which is reflected in the figure. Linear fits to the data points in the plots for ω 2 e f f and Q 2 e f f are used to calculate ω 2 0 and Q 2 0 , which are the intercepts of the ω 2 e f f and Q 2 e f f fits.

Model Validation
The extracted parameters were used in a Verilog-A implementation of Equation (2), and the Verilog-A model was simulated using hSpice. In the simulation, the Verilog-A model was connected as pseudo-inverter using the same parameter values described in Section 2.3. The model's simulated velocity is compared with LDV results in Figure 6, which overlays simulated and measured velocity data. Though the measured and simulated velocity data in Figure 6 appear to agree well, figures of merit were necessary to compare the simulation to the measurement. Two such figures of merit were selected: frequency of oscillation and the peak velocity. Figure 6 shows the error in frequency of oscillation and peak velocity as V gb is varied. The error in these figures of merit is never larger than 10%.

Conclusions
A normalized model of a MEM relay was developed and the parameters of that model were extracted from LDV measurements. The model's simulated velocity matched measurements over a wide range of operating voltages even though the free-space model only had four parameters.
This improved model could immediately benefit simulations and analysis in publications that follow in the rich tradition of using relay simulations to explore circuit designs, as in [3,[8][9][10]. The low complexity and good convergence properties of the model make it particularly suitable for simulating large circuits, as in [10].
This model is tightly tied to parallel plate electrostatic actuators because it relies so heavily on the pull-in voltage as a normalizing term. A natural path for future work is extending this normalization procedure to other capacitance-vs-displacement profiles. Such work could also account more accurately for fringing capacitance, an outstanding inaccuracy in this model at high V gb .

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

Abbreviations
The following abbreviations are used in this manuscript:

MEM Microelectromechanical FEA
Finite element analysis LDV Laser Doppler vibrometry