Thermoelastic damping in MEMS gyroscopes at high frequencies

Microelectromechanical systems (MEMS) gyroscopes are widely used, e.g., in modern automotive and consumer applications, and require signal stability and accuracy in rather harsh environmental conditions. In many use cases, device reliability must be guaranteed under large external loads at high frequencies. The sensitivity of the sensor to such external loads depends strongly on the damping, or rather quality factor, of the high-frequency mechanical modes of the structure. In this paper, we investigate the influence of thermoelastic damping on several high-frequency modes by comparing finite element simulations with measurements of the quality factor in an application-relevant temperature range. We measure the quality factors over different temperatures in vacuum, to extract the relevant thermoelastic material parameters of the polycrystalline MEMS device. Our simulation results show a good agreement with the measured quantities, therefore proving the applicability of our method for predictive purposes in the MEMS design process. Overall, we are able to uniquely identify the thermoelastic effects and show their significance for the damping of the high-frequency modes of an industrial MEMS gyroscope. Our approach is generic and therefore easily applicable to any mechanical structure with many possible applications in nano- and micromechanical systems.


I. INTRODUCTION
Microelectromechanical systems (MEMS) gyroscopes are well established and indispensable in modern consumer and automotive electronics [1,2].Especially in automotive applications, where gyroscopes operate in safety-critical systems, device reliability is of utmost importance.Functionality has to be ensured under various harsh environmental conditions [3] and the sensor signal stability has to be maintained despite many adverse linear and nonlinear effects [4,5].Most importantly, sensors have to withstand temperatures ranging from −40 • C to 120 • C and should be insensitive against external vibrations [2].Therefore, the ability to predict the sensitivity of the sensor to such external conditions is crucial during MEMS design.In the past, vibrational robustness was mainly concerned with frequencies up to a few tens of kHz [2,6].State of the art applications, however, are faced with ever-increasing requirements.Among these requirements is the robustness against large external loads, at frequencies much higher than the operational frequency of the oscillatory gyroscope.High eigenfrequency modes, far beyond the operational frequency, can be decisive for the response of the sensor.The response of the corresponding high frequency modes is, among other quantities, determined by their quality factors.At typical pressures of around a few millibar, the quality factors of low frequency modes are known to be limited by gas damping [7,8].However, to the authors' knowledge, there has been no exhaustive research on the damping of high frequency modes in MEMS gyroscopes.Known damping mechanisms that can contribute to the quality factors of MEMS resonators are gas damping, thermoelastic damping (TED), anchor losses, surface losses, material losses and Akhiezer damping [8][9][10][11][12][13][14][15][16].The first three are usually considered as the dominant mechanisms in polysilicon MEMS resonators.Material losses are considered negligible for silicon, as it exhibits very linear material behavior, and surface losses are mainly relevant in nanoresonators [10][11][12]15].Akhiezer damping is only expected to be relevant for frequencies above 10 MHz [16] and for very high quality factor and frequency (Q-f ) products [17].
In this paper, we compare measured and simulated quality factors of industrial MEMS gyroscopes over a wide range of eigenfrequencies.The aim of this work is to illuminate the significance of the thermoelastic damping contributions for high frequency modes.We show that thermolastic damping indeed limits the quality factor of high-frequency modes of the gyroscope and is thus crucial for the gyroscope's response to high frequency vibrations.
In Section II we introduce the MEMS devices and the measurement method.In Section III the governing equations of thermoelasticity are introduced and an efficient method to simulate thermoelastic damping, based on the finite element method (FEM), is derived.We then verify the validity of our method in Section IV, by comparing our simulation results to measured data.Finally, in Section V, we summarize our results and conclude that thermoelastic damping is highly relevant for high frequency modes in our devices.

II. EXPERIMENTAL SETUP
Two different industrial three-axis MEMS gyroscope designs (A and B), developed by Bosch, were investigated.The devices are made of polycrystalline silicon and are therefore assumed to exhibit isotropic material behavior.The designs were measured with two different scanning laser Doppler vibrometers (SLDV) from Polytec.The oscillation modes of the gyroscopes were excited in the linear regime by broadband signals (see details below).The measurements were performed on a dense grid of points over the structures (see details below) and the spectra of velocity and displacement were obtained from a fast Fourier transform at each point.The measured displacement maps obtained from the grid allowed to identify the vibrational modes by comparing with the simulated mode shapes.The quality factors of the modes were obtained from the linewidths of the resonance peaks.The frequency resolution of ∼1 Hz was sufficient for an accurate resolution of the peaks.
Design A is an unencapsulated single chip, that was held at 1 mbar and 25 • C inside a vacuum chamber.The excitation was realized dominantly in out-of-plane direction via a piezo-shaker.A chirp signal in the range from 10 kHz to 2 MHz was applied to the piezo-shaker.The measurement was performed with a 1D SLDV.Therefore, only out-of-plane modes were detected for design A. The measurement was performed on a grid of around 400 points over the structure.The measured out-ofplane modes were semi-automatically matched to simulated modes.Although this is prone to errors, it enables the investigation of quality factor trends over a wide range of eigenfrequencies.
Design B was measured on the wafer and not encapsulated.The excitation of design B was realized electrostatically.The design contains a capacitive comb-drive as well as three different electrode pairs for capacitive sensing.Applying an electric broadband signal to one of the drive or sense electrode pairs, while grounding the remaining electrodes, enabled the excitation of various in-plane or out-of-plane modes.The applied signal was a pseudo-random broadband signal in the range from 30 kHz to 1.25 MHz.The wafer containing design B was mounted on a thermal chuck inside a vacuum chamber.Thus, temperature and pressure could be varied.A 3D SLDV was used to measure in-plane and out-ofplane modes of design B. For the measurement of the out-of-plane modes a grid of around 100 points over the structure was used.The in-plane modes were measured inidividually with a grid of around 30 points each, which were only placed on the relevant oscillating parts of the structure.The measured modes of design B were manually identified with simulated modes, based on mode shapes and eigenfrequencies.
Fig. 1 shows the measured quality factors of design A for several out-of-plane modes up to an eigenfrequency of 1.8 MHz at 1 mbar and 25 • C. The pressure of 1 mbar is a typical operational value for MEMS gyroscopes [6].Ad- ditionally, simulated quality factors based on gas damping are also included in the figure.The gas damping simulations have been done using a Bosch internal gas damping simulation tool based on molecular flow simulations in COMSOL [18].The validity and precision of the gas damping simulation is highlighted in the inset of Fig. 1 showing a closeup of the frequency regime up to 200 kHz.Up to 200 kHz, the measured quality factors follow the trend of the simulations.However, for higher eigenfrequencies the measured quality factors clearly saturate.Gas damping quality factors, on the other hand, increase approximately linearly with eigenfrequency [7,12].This motivates the incorporation of additional damping mechanisms into the simulation, to identify and accurately predict the damping contributions and the total quality factor.In this work, we will investigate the influence of TED on the two MEMS gyroscope designs.

III. NUMERICAL ANALYSIS OF THERMOELASTIC DAMPING
Thermoelastic damping arises naturally from the coupling of the displacement and temperature fields.Therefore, any material with a non-zero thermal expansion coefficient exhibits TED.When a thermoelastic structure with a positive thermal expansion coefficient oscillates, regions under compression will heat up and regions under tension cool down.Thus, oscillatory temperature gradients arise across the structure.The resulting periodic heat flow along the temperature gradients is irreversible and leads to dissipation of energy.Zener pioneered the research on TED and derived an approximate analytic equation for the corresponding quality factor of a beam's fundamental bending mode [19,20].Lifshitz and Roukes later derived a refined solution for the same problem [21].In order to obtain quality factors for arbitrary geometries, the finite element method can be employed.It has been shown, that TED quality factors can be obtained from a complex eigenvalue problem of the thermoelastic equations [22] or from the calculation of dissipated and stored energy [23][24][25].
We show how TED of industrial scale problems can be modelled by deriving a modified equation of motion for the mechanics, which can then be transferred into a mechanical reduced order model (ROM).In Section III A the governing differential equations of continuum mechanics are introduced along with their FEM formulation.In Section III B an efficient simulation method for the evaluation of TED quality factors is derived.

A. Governing Equations
We start with the fundamental equations of thermoelasticity, which can be found e.g. in [26,27].The governing equation of the mechanical response, i.e. the equation of motion, is given by the linear momentum balance.For small deformations and in the absence of body forces it reads div(σ) = ρa, where σ is the stress tensor, ρ is the density and a is the acceleration vector.The coupling to the temperature field affects Eq. (1) via thermal expansion.As this work is concerned with structures made of polysilicon, linear and isotropic material behavior will be assumed for mechanical and thermal properties.The constitutive equation accompanying Eq. ( 1) is given as where C is the fourth-order elasticity tensor, ε is the total strain tensor, 1 is the second-order unit tensor, α is the thermal expansion coefficient and ∆T is the difference between the temperature field T within the body and the ambient temperature T 0 , i.e. ∆T = T − T 0 .The second term in the bracket of Eq. ( 2) signifies the strain due to thermal expansion.The temperature changes that result from the thermoelastic coupling are generally very small.Therefore, the heat equation, which determines ∆T , is linearized around T 0 as with heat flux vector q t , specific heat C V and time derivatives denoted by dots above the symbols.It is assumed that no additional heat sources are present within the body.The coupling to the stress field in Eq. (3) manifests itself in the heating of regions under compression and cooling of regions under tension, if the thermal expansion coefficient is positive.The constitutive equation for the heat flux vector is given by Fourier's law where κ is the thermal conductivity.The global FEM equations can be obtained in the usual way, by deriving and discretizing the weak forms of the local equations ( 1) and (3), leading to where M is the mass matrix, K u the stiffness matrix, K ut the thermoelastic coupling matrix, C t the specific heat matrix, K t the thermal conductivity matrix, u and ∆T are the nodal displacement and temperature change vectors and f is the external force vector.Only the oscillating structure is considered in the simulations.In Eqs. ( 5) and ( 6) we assume that the displacement and temperature change are zero at the connection of the oscillating structure to the substrate.Furthermore, in Eq. ( 6) we assumed insulating boundary conditions on the boundary that isn't fixed.See e.g.[28] for the definitions of the FEM matrices.

B. Solution Method
Several approaches exist to evaluate the thermoelastic damping of mechanical modes based on Eqs. ( 5) and (6).Common but computationally expensive methods solve the coupled Eqs. ( 5) and ( 6) simultaneously.However, these methods require the solution of non-symmetric equation systems with 4N degrees of freedom for a mesh with N nodes.In the modelling of MEMS gyroscopes one usually deals with models where N > 10 6 and quality factors have to be calculated for many modes over a wide frequency range.Therefore, solving the coupled problem of Eqs. ( 5) and ( 6) is time consuming and numerically expensive.Instead, we will take a different approach, where we eliminate the heat equation and arrive at an effective equation of motion, which can then be efficiently evaluated in a ROM.
We consider the case where Eq. ( 5) is harmonically driven at a frequency ω.Thus, the steady-state oscillations of displacement and temperature change are given as where u 0 and ∆T 0 are the complex steady-state amplitudes.Equations ( 7) and ( 8) are inserted into the heat equation ( 6), which can then be formally solved for ∆T 0 .Consequently, one can then express the temperature change ∆T , based on Eq. ( 8), in dependence of displacement u and velocity u as where A = (K t + iωC t ) −1 (K ut ) T .Substituting Eq. ( 9) into the equation of motion (5), we obtain the modified equation of motion with damping matrix and stiffness matrix Note that Eq. ( 10) is still exact in the sense that it fully incorporates the effect of the thermoelastic coupling on the mechanics for harmonic forcing.For oscillatory structures, such as MEMS gyroscopes, the equation of motion is usually solved in a modal ROM.The mechanical modes are obtained from the purely mechanical eigenvalue problem with eigenfrequency ω n and mode shape φ n of the nth mode.The mode shapes are mass-normalized, i.e. φ T n M φ n = 1.The displacement is then expressed as a superposition of the modes where q n is the modal coordinate of mode n and Φ = [φ 1 φ 2 ... φ m ] is a matrix, which contains the massnormalized eigenvector of mode n in the n-th column.
The index m indicates the mode at which the superposition is truncated, leading to an approximation of the actual u.Inserting the modal superposition given by Eq. ( 14) into Eq.( 10) and multiplying by Φ T from the left one obtains The effect of the thermoelastic coupling thus influences the mechanical modes by a damping contribution as well as a change in stiffness, i.e. a change of the eigenfrequencies.Furthermore, the modal damping and stiffness matrices Φ T CΦ and Φ T KΦ are not diagonal, i.e. they lead to a linear coupling between modes.This is simply a manifestation of the two-way coupling of Eqs. ( 5) and (6).The temperature field that results from the motion of a mechanical mode and is determined by Eq. ( 6) may also impose forces on other mechanical modes in Eq. ( 5), providing an intermodal coupling.We assume that the effect of this coupling is weak and thus only consider the diagonal entries in Eq. (15).Furthermore, the change of eigenfrequency due to thermoelastic coupling is very small and therefore only a very small error is made by neglecting it.
The damping matrix C depends on the oscillation frequency ω.In this work, we are interested in the damping of a mode at its eigenfrequency ω n .Hence, to obtain the quality factor of mode n, one can set ω = ω n .The reciprocal quality factor due to thermoelastic damping T ED,n is found by dividing the n-th diagonal entry of Φ T CΦ by ω n , leading to ) Remarkably, Eq. ( 16) allows us to determine the quality factors by only having to solve a symmetric linear equation system of size N , i.e. the size of the temperature degrees of freedom, per mode.Therefore, this approach is much more efficient than solving the coupled equations directly and is suitable for large models.We have implemented the assembly of the FEM matrices and the evaluation of Eq. ( 16) in a self-written Matlab code.
We note that Eq. ( 16) is equivalent to the result obtained with a perturbation method in [29].Furthermore, we note that the same expression can be obtained by calculating the quality factor as the ratio of stored to dissipated energy, if one calculates the dissipated energy due to the temperature field given by Eq. ( 9) and neglects the effect of temperature on the stored energy.

IV. RESULTS
The main damping mechanisms in MEMS resonators are gas damping, thermoelastic damping and anchor losses.Other possible damping mechanisms include material losses and surface losses.Material losses are known to be negligible for silicon and surface losses are mainly relevant for nanoresonators [10,11].Additionally, Akhiezer damping has been observed in silicon MEMS resonators, but is only expected to be relevant for frequencies above 10 MHz [16] and for very high Q-f products [17].
From here on, when we refer to temperature, we mean the temperature T 0 of the atmosphere surrounding the oscillating part of the structure, i.e. the temperature inside the vacuum chamber.
Gas damping depends on temperature T 0 and pressure p, thermoelastic damping only depends on temperature and anchor losses are assumed to be independent of pressure and temperature.The total reciprocal quality factor is obtained from the sum of the reciprocal quality factors of the individual damping mechanisms with total quality factor Q, gas damping quality factor Q gas , thermoelastic damping quality factor Q TED and anchor loss quality factor Q anchor .In principal, as already mentioned, there are also other damping mechanisms that contribute to Eq. (17).We assume that these other damping mechanisms are negligible compared to the gas damping, TED and anchor losses.We note, however, that other temperature-and pressure-independent damping mechanisms would not be distinguishable from anchor losses in our measurements.The dependence of Q gas on experimental conditions is particularly simple.At very low pressures, in the molecular regime, it scales as Q −1 gas ∝ p.At higher pressures, a transition into the viscous gas damping regime occurs, where the dissipation scales as Q −1 gas ∝ √ p [12].In the molecular regime, if the pressure isn't controlled, the dissipation scales with temperature as Q −1 gas ∝ √ T 0 [30].In order to verify that Q TED is determined by Eq. ( 16), we measured 7 different modes of design B. In contrast to design A, which was excited via a piezo-shaker, design B was excited electrostatically.Due to the placement of the electrodes, only certain mode shapes were excitable.Thus, it wasn't possible to excite as many modes for design B as for the out-of-plane measurements of design A. Out of the measured modes, we chose those that could be identified unambiguously with simulated mode shapes and exhibited a clear resonance peak in our measurements.This lead to the 7 modes, which are enumerated by letters a to g, from lowest to highest eigenfrequency.
The lowest measured mode is mode a with a simulated eigenfrequency of f 0 = 118.98kHz and the highest measured mode is mode g with a simulated eigenfrequency of f 0 = 733.46kHz.Out of the 7 measured modes, 2 are in-plane modes and the remaining 5 are out-of-plane modes.

A. Gas Damping
Since Q −1 gas ∝ p in the molecular regime, gas damping can be made negligible by reducing the pressure sufficiently.Figure 2 shows the measured quality factors of the 7 measured modes of design B over pressure at a temperature of 20 • C. For each pressure the quality factor of every mode was measured at 8 different spots on the MEMS structure.The spots were chosen individually for each mode according to the mode's anti-nodes.From the 8 measurements the mean value was calculated and the standard deviation was used for the vertical errorbars.Furthermore, a fit is shown, which was obtained for each mode from the linear relationship of the reciprocal quality factor and pressure, i.e Q −1 = m • p + b, where m is the slope and b is the pressure-independent offset.It can be seen that all modes follow this expected trend, which confirms that the measurements were performed in the molecular regime.The quality factors only show very little pressure dependence below 10 −2 mbar.Subsequent measurements were performed at 10 −3 mbar, to ensure that the gas damping contribution is negligible and the measured quality factors are in good approximation equal to the contributions from thermoelastic damping and anchor losses.

B. Material Parameters
In order to verify Q −1 TED according to Eq. ( 16), the correct temperature dependence has to be taken into account.At first sight Eq. ( 16) appears to be linear in T 0 .However, Q −1 TED also depends on thermal expansion coefficient α, thermal conductivity κ and specific heat C V , which exhibit significant temperature dependencies.On the other hand, Young's modulus E, Poisson's ratio ν and density ρ have much smaller temperature dependencies, which are negligible in this context.To make the dependence on temperature-dependent material properties more explicit, we rewrite Eq. ( 16) as where we defined K ut = α Kut , K t = κ Kt and C t = C V Ct , so that Kut , Kt and Ct are then independent of α, κ and C V .It is clear that Eq. ( 18) scales with α 2 .Therefore, the thermal expansion coefficient α affects every mode in the same way.Thermal conductivity κ and specific heat C V , on the other hand, affect every mode in a different way, due to their appearance within the inverse matrix in Eq. (18).To predict the quality factors accurately over temperature, the correct temperature dependencies of the material parameters have to be taken into account.For the purely mechanical properties, we assumed constant values of E = 161 GPa, ν = 0.22 and ρ = 2330 kg m −3 , which are the standard values used at Bosch for polycrystalline silicon.Due to a lack of reported data for polysilicon, the temperature-dependent specific heat C V was calculated from the Debye model with a Debye temperature of 645 K for silicon [31].The Debye model for silicon has also been used by others in the context of TED [30], albeit for monocrystalline silicon.We assume that the polycrystallinity has no significant impact on the specific heat.The value of κ depends strongly on doping concentration and film thickness.Reported room temperature values for polysilicon samples of various doping concentrations and film thicknesses lie between 15 W m −1 K −1 and 60 W m −1 K −1 [32].However, our samples have a film thickness of a few dozen micrometers, while the reported samples in [32] are significantly thinner.The thermal conductivity is expected to increase with film thickness and decrease with doping concentration [32].Therefore, a thermal conductivity above 60 W m −1 K −1 would be realistic for sufficiently low doping concentration.The thermal expansion coefficient α of monocrystalline silicon over temperature is well documented [33].However, there exists no conclusive data for polycrystalline silicon.It has been suggested that the thermal expansion coefficient of polycrystalline silicon thin films might be significantly higher than that of bulk monocrystalline silicon [34,35].Other researchers have performed measurements that found the thermal expansion coefficient of polycrystalline silicon to be constant over temperature and only slighty higher than that of monocrystalline silicon [36].Furthermore, it has been indicated in [37] that the thermal expansion coefficient of polycrystalline silicon differs from that of monocrystalline silicon depending on residual stresses.We conclude that there is ambiguous data on the temperature dependence of κ and α for polysilicon thin films.Therefore, we will treat them as fit coefficients and estimate them from measured quality factors.
Quality factors of design B were measured over temperature at a pressure of 10 −3 mbar.The measurements span the range from −21 • C to 114 • C. The bond pads of the device that was measured over pressure were already damaged from contacting them multiple times.Therefore, we performed the temperature measurements on a second device of design B on the same wafer.Due to the low pressure, the measured Q −1 is approximately equal to the contributions from TED and anchor losses.We added the mean difference between simulated Q −1 TED and measured Q −1 for each mode to Q −1 TED , to emulate the effect of temperature-independent anchor loss contributions Q −1 anchor .We then determined α and κ such that the simulated Q −1 TED over temperature matched the measured Q −1 after the mean difference was added.The measured Q −1 (blue circles), simulated Q −1 TED (red dots) according to Eq. ( 16) and Q −1 TED with added offset (black dots) are shown in Fig. 3 over temperature.It can be seen that the simulated Q −1 TED is lower than the measured Q −1 for all modes, as expected.The corresponding Q values are shown in Fig. 4, revealing that the Q changes by a factor of ∼2 over the application-relevant temperature range.In order to find suitable values for α and κ, we assumed simple temperature-dependencies.In our temperature range, the thermal expansion coefficient α of monocrys-FIG.3. Reciprocal quality factors over temperature for the 7 measured modes of design B, device 2. The measurements were performed at a pressure of 10 −3 mbar.The blue circles show the experimental values.The red dots were obtained from Eq. ( 16) with temperature-dependent material properties as shown by the solid lines in Fig. 5.The black dots were obtained by adding the mean difference between simulation and measurement to the simulated red dots, in order to emulate the effect of temperature-independent anchor losses.The plots in (c) and (e) show in-plane modes, while the remaining measurements show out-of-plane modes.The corresponding modes and their simulated eigenfrequencies are: talline silicon increases with a declining slope over temperature [33].We assume a qualitatively similar behavior for polycrystalline silicon.For simplicity, we choose a quadratic function for α(T 0 ) where α 0 , α 1 and α 2 are positive fit parameters and T RT = 25 • C. The thermal conductivity κ of silicon tends to decrease with 1/T 0 in our temperature range [32].To mimic this behavior, we choose the fit function for κ(T 0 ) in analogy to [38] as where κ 0 and κ 1 are positive fit parameters and T low = −21 • C is the lowest measured temperature.The material parameters were estimated from least-square fits.For that purpose, mode b was excluded, as it showed irregular behavior.The measured quality factor of mode b was three times larger for the device that was measured over temperature than for the device that was measured over pressure, as can be seen from Fig.Although we can not explain this peculiarity, we expect it to not be related to thermoelastic damping.
In order to estimate the five unknown parameters α 0 , α 1 , α 2 , κ 0 and κ 1 , we proceeded as follows: As a first estimate, we assumed the α(T 0 ) of monocrystalline silicon according to [33], which is shown as the blue dashed line in Fig. 5.We then chose κ(T 0 ) according to Eq. ( 20) and determined κ 0 and κ 1 from a least-square fit.For that purpose, we added the mean difference between measurement and simulation to the simulated Q −1 TED of each mode, to account for anchor losses.We then calculated the squared differences between the resulting values and the measured Q −1 values.The squared differences for all modes over all temperatures, except mode b, were then added and the material parameters were determined from the minimization of this sum.The κ(T 0 ) that was obtained from this first fit, based on the α(T 0 ) of monocrystalline silicon, is shown as a red dashed line in Fig. 5. Subsequently, this κ(T 0 ) was assumed and α(T 0 ) was chosen as in Eq. (19).The coefficients α 0 , α 1 and α 2 were then determined from a least-square fit, similar to the previous one, only differing in the unknown parameters.The resulting α(T 0 ) is shown as a solid blue line in Fig. 5 and can be seen to be lower than the α(T 0 ) of monocrystalline silicon.Finally, to obtain a corrected estimate for κ(T 0 ), we assumed this new α(T 0 ) and treated κ 0 and κ 1 as fit parameters again.From this final fit we then obtained the κ(T 0 ) according to the solid red line in Fig. 5.For verification, we performed a final iteration, where we assumed κ(T 0 ) according to the solid red line in Fig. 5 and again obtained the α(T 0 ) according to the solid blue line.Therefore, we conclude that the fitting process has converged to some minimum.The final material parameters, according to the solid red and blue lines in Fig. 5, will be assumed for the remainder of this paper.Figures 3 and 4 show the Q TED based on these final material parameters.

C. Damping Contributions
If α(T 0 ) were to be rather constant over temperature, as suggested in [36], then the simulated values would not be steep enough to match the measured Q −1 slope over T 0 .The α(T 0 ), that was determined via fitting, leads to a good agreement with the measured Q −1 over T 0 slope and is in a reasonable range compared to the monocrystalline values.The κ(T 0 ) values are higher than the reported values in [32], but seem reasonable due to the few dozen micrometer thickness of our samples.The change of around 10 W m −1 K −1 of κ(T 0 ) over our temperature range is similar to the data in [32].However, we note that there might be other possible combinations of α 0 , α 1 , α 2 , κ 0 and κ 1 , which also lead to good agreement with the measured data.In prinicipal, this is related to the unknown temperature-independent anchor loss contributions, which manifest themselves as an offset between Q −1 TED and Q −1 .There might be other curves for α(T 0 ) and κ(T 0 ) that yield a similar slope of Q −1 TED over T 0 and only differ in their offset to Q −1 .Additionally simulating the anchor losses might help to circumvent this issue and will be subject of a future publication.Nevertheless, we conclude that the temperature dependence of the measured Q −1 can be very well reproduced by the Q −1 TED according to Eq. ( 16) with reasonable material parameters.This clearly shows the validity of Eq. ( 16) and confirms the assumption that the temperature dependence of Eq. ( 17), at pressures where gas damping is negligible, is due to Q −1 TED for the investigated MEMS gyroscope.In Fig. 4 it can be seen that the Q at 10 −3 mbar increases by roughly a factor of two from the highest to FIG. 5. Thermal expansion coefficient α and thermal conductivity κ over temperature.The dashed blue line shows the α of monocrystalline silicon according to [33] and was used as an initial guess.Subsequent material parameters were obtained by least-square fitting the black dots to the measured blue circles in Fig. 3.The dashed red line shows an initial fit for κ.The red and blue solid lines show the converged result of our fitting method for the thermal conductivity κ and the thermal expansion coefficient α, respectively.the lowest temperature.The temperatures are all within the application-relevant temperature range.To ensure device functionality over the whole temperature range, knowledge of quality factors is crucial.The significant change over temperature highlights the relevance of TED for the development of MEMS gyroscopes.
In Fig. 6 the Q −1 contributions are shown for a pressure regime where gas damping is relevant.For that purpose, the contributions were obtained from the fit gas ∝ √ T 0 , if the pressure isn't held constant, the gas damping can be calculated for other temperatures based on the known value at 20 • C.This is relevant if the device would be encapsulated, where the pressure would also change over temperature.Based on Fig. 2, we extracted the gas damping at 1 mbar and 20 • C and then scaled it up to 114 • C, which also corresponds to a higher pressure.Q −1 TED is also shown for 20 • C and 114 • C, based on Eq. ( 16).Q −1 anchor is temperature-independent and therefore identical for both temperatures.Looking at Fig. 6, one can see that modes a to d are gas damping dominated.However, TED gains significance at the higher temperature of 114 • C.This is due to the fact that, at least in our temperature range and with our devices, Q −1 TED scales roughly linear with T 0 , as seen in Fig. 3, while Q −1 gas only scales with √ T 0 .On the other hand, for the higher modes e to g, one can see that even at 20 gas and for modes f and g even larger than Q −1 gas .At FIG. 6. Contributions of gas damping, TED and anchor losses to the total Q −1 for the 7 modes of design B. The Q −1 gas was obtained from the fit in Fig. 2 at 1 mbar and 20 • C (light green).Assuming non-constant pressure it was then scaled to 114 • C (dark green), according to Q −1 gas ∝ √ T0.Q −1 TED was obtained from Eq. ( 16) for 20 • C (light red) and 114 • C (dark red).The temperature-and pressure-independent Q −1 anchor (gray) was calculated by subtracting Q −1 TED at 20 • C from the yintercept of the linear fit for Q −1 , that was shown in Fig. 2. 114 • C this becomes even more pronounced.The Q anchor corresponding to Fig. 6 are between 70000 and 110000.The only exception is mode b with a Q anchor of around 8500.Again, we note that the Q-factor of mode b from the temperature and pressure measurements, which were performed on different but nominally identical devices, differ significantly.One can also obtain Q −1 anchor as the offset in Fig. 3.This leads to comparable values as the ones in Fig. 6, except for mode b, for which the extracted Q anchor would then be 48000.Mode b is a symmetric mode, whereas mode a is the corresponding antisymmetric mode.Therefore, the deviation for mode b can not be attributed to a simple measurement error, as the entire spectrum was measured simultaneously and one would then expect a discrepancy for at least mode a as well.We also exclude that a MEMS device was faulty, as we would then also expect to see significant deviations for the other modes.One possible explanation could be a substrate-related effect, which was only present in one of the two devices.

D. Damping Contributions over a Wide Frequency Spectrum
The same quality factors of design A as in Fig. 1  TED quality factors are included as red dots.Design A is made of the same material and features the same outof-plane thickness as design B. Therefore, the material parameters according to the solid blue and red lines in Fig. 5 at 25 • C were used.The quality factors resulting from gas damping plus TED are shown as black dots.Furthermore, trend lines are included as a guide to the eye.The red horizontal line signifies that there is no clear up-or downward trend for TED quality factors over the measured range.It was calculated as the median of the simulated TED quality factors, to avoid the bias by large outliers.The green line is a linear fit of the simulated gas damping quality factors over frequency.The black line was obtained from the previous two trend lines, i.e. based on the trend of TED and gas damping.Additionally, in order to estimate possible anchor losses, a constant Q anchor of 75000 was added to the contributions of the black line, to obtain the blue line.
One can see that up to 200 kHz the influence of TED is rather small compared to gas damping, as Q TED is much larger than Q gas .Therefore, Q is well approximated by Q gas alone.However, the incorporation of TED into the model significantly decreases the simulated Q for higher frequencies.For many modes this leads to a reduction of the simulated Q by a factor of around 2. In fact, the simulation based on gas damping and TED is quite close to the measurements up to around 800 kHz.Above 800 kHz almost all measured modes exhibit a lower quality factor in the measurement than in the simulation.For the highest modes the measurement is still around 25% below the simulation of gas damping and TED.For the modes above 800 kHz anchor losses play an increasing role.TED quality factors show no clear up-or downward trend over our measurement range.Gas damping quality factors increase approximately linearly with frequency.Therefore, for sufficiently high modes the gas damping quality factors approach the value of the TED quality factors and the effect of TED becomes relevant.For the modes above 800 kHz, the gas damping quality factors are large enough, such that further damping mechanisms, which tend to have larger quality factors than TED, become relevant as well.For illustration, a Q anchor of 75000 was therefore included in the blue trend line.The blue trend line then leads to a stronger saturation of Q over mode frequency compared to only including gas damping and TED and therefore gives a better approximation of the measured data at high frequencies.However, this is only a phenomelogical remedy and it could also be that Q anchor has some trend over frequency.Clarifying the influence of anchor losses and investigating whether the incorporation of an anchor loss simulation helps to explain the remaining deviation between simulation and measurement at high frequencies will be part of future research.
Note that the simulated Q gas is a simplified approximation.Furthermore, the measured data of design A in Fig. 7 is not as accurate as the measurements on design B. This is due to the fact that measuring and identifying as many modes as in Fig. 7 is very elaborate and the semiautomated evaluation of the quality factors of design A from the measured data is prone to errors.Therefore, deviations between measurement and simulation are to be expected for individual modes.Regardless, the measurement in Fig. 7 shows a clear trend.The simulations and the measured trend over the broad spectrum clearly show the significance of TED with increasing frequency.For higher temperatures or lower pressures TED would become even more important.The relevance of TED for higher modes is also in agreement with the observations made from design B, as was shown in Fig. 6.

V. CONCLUSION
We reported quality factor measurements of an industrial MEMS gyroscope (design A) for a multitude of outof-plane modes over a wide frequency range.Although gas damping matches the observations for the first few modes, up to 200 kHz, we found that for higher modes the measured quality factors are significantly lower than the pure gas damping model predicts.The measured quality factors saturate for high frequency modes above 800 kHz, while gas damping quality factors increase approximately linearly with eigenfrequency.The deviation starts to become notable above 200 kHz.In order to account for this deviation, we introduced thermoelastic damping into our model.
We demonstrated an efficient way to simulate thermoelastic damping, by eliminating the heat equation and deriving an effective equation of motion for the harmonically driven mechanics.The FEM equations were implemented in a self-written code.In order to verify our thermoelastic damping simulations, we measured the quality factors for 7 different modes of a second MEMS gyroscope design (design B) over temperature in a vacuum chamber.We found a good agreement between simulated and measured quality factors over temperature.We showed that the temperature dependence of thermoelastic material parameters, in particular the thermal expansion coefficient, has to be taken into account to reproduce the observed behavior over temperature.Fitting our simulation results to the measurements, we were able to es-timate the thermal conductivity and thermal expansion coefficient over temperature.However, we note that the presence of damping mechanisms which are independent of pressure and temperature, such as anchor losses, leads to some uncertainty of the fitted material parameters.
Having validated the thermoelastic damping simulations and having obtained the relevant material parameters, we then applied our method to the measurement over a wide frequency range of gyroscope design A. Taking thermoelastic damping and gas damping into account, we then found good agreement with measured quality factors up to 800 kHz.For even higher frequencies, we found that additional damping mechanisms might be relevant, as the simulated quality factors were still higher than the measured ones.Nevertheless, the simulated quality factors were significantly closer to the measured ones over the entire frequency range, as compared to pure gas damping simulations.
Our results clearly show the significance of thermoelastic damping for high frequency modes in MEMS gyroscopes.On the lower end of the mode spectrum, gas damping dominates and is the sole relevant damping mechanism.However, for increasing eigenfrequency, thermoelastic damping gains relevance and appears to be the limiting damping mechanism, i.e. it has the smallest quality factor.At the upper end of the measured spectrum, additional damping mechanisms might have to be taken into account, although thermoelastic damping still seems to be the primary contribution.Expanding our model with anchor loss simulations, in order to improve the accuracy and also investigate the origin of the additonal damping contributions at high eigenfrequencies, will be subject of future research.Finally, as our method is based on FEM, it can be readily applied to other MEMS and NEMS structures where quality factors and damping mechanisms are also subject of ongoing research [39][40][41].We would like to emphasize that, although our discussion was focused on frequency modes, our method is generally applicable also for low frequency modes.

FIG. 1 .
FIG. 1. Quality factors for all out-of-plane modes up to 1.8 MHz of design A plotted over the eigenfrequencies of the modes at 1 mbar and 25 • C. The green dots show the simulated quality factors based on gas damping.The solid green line is a linear fit of the simulated quality factors over frequency and indicates the trend of the gas damping quality factors.The blue dots show the measured quality factors.The inset shows a magnified plot up to 200 kHz.

FIG. 2 .
FIG. 2. Quality factors over pressure for the 7 measured modes of design B, device 1.The measurements were performed at a temperature of 20 • C. The blue circles show the experimental values.The red curves were obtained from a linear fit of the reciprocal quality factors, i.e.Q −1 = m • p + b.The plots in (c) and (e) show in-plane modes, while the remaining measurements show out-of-plane modes.The corresponding modes and their simulated eigenfrequencies are: (a) Mode a with f0 = 118.98kHz, (b) Mode b with f0 = 121.77kHz, (c) Mode c with f0 = 131.91kHz, (d) Mode d with f0 = 188.44kHz, (e) Mode e with f0 = 331.01kHz, (f) Mode f with f0 = 629.11kHz, (g) Mode g with f0 = 733.46kHz.
FIG.3.Reciprocal quality factors over temperature for the 7 measured modes of design B, device 2. The measurements were performed at a pressure of 10 −3 mbar.The blue circles show the experimental values.The red dots were obtained from Eq. (16) with temperature-dependent material properties as shown by the solid lines in Fig.5.The black dots were obtained by adding the mean difference between simulation and measurement to the simulated red dots, in order to emulate the effect of temperature-independent anchor losses.The plots in (c) and (e) show in-plane modes, while the remaining measurements show out-of-plane modes.The corresponding modes and their simulated eigenfrequencies are: (a) Mode a with f0 = 118.98kHz, (b) Mode b with f0 = 121.77kHz, (c) Mode c with f0 = 131.91kHz, (d) Mode d with f0 = 188.44kHz, (e) Mode e with f0 = 331.01kHz, (f) Mode f with f0 = 629.11kHz, (g) Mode g with f0 = 733.46kHz.
4 (b)  and Fig.2 (b).The devices are nominally identical and one would expect similar quality factors, as was the case for all other modes.Furthermore, the quality factor measurement of mode b showed a clear outlier over temperature, as seen in Fig.3(b) and Fig. 4 (b) at 114 • C. We performed a measurement of mode b at 10 −3 mbar and 25 • C on a third device and obtained a quality factor of 27000±1000.This is closer to the value of 23000 ± 1200 in Fig. 4 (b).

FIG. 4 .
FIG. 4. Quality factors over temperature for the 7 measured modes of design B, device 2. The measurements were performed at a pressure of 10 −3 mbar.The plots show the quality factors corresponding to the reciprocal values from Fig. 3. Blue circles show measured values, red dots show simulated TED and black dots show simulated TED with an additional temperature-independent damping contribution included.The plots in (c) and (e) show in-plane modes, while the remaining measurements show out-ofplane modes.The corresponding modes and their simulated eigenfrequencies are: (a) Mode a with f0 = 118.98kHz, (b) Mode b with f0 = 121.77kHz, (c) Mode c with f0 = 131.91kHz, (d) Mode d with f0 = 188.44kHz, (e) Mode e with f0 = 331.01kHz, (f) Mode f with f0 = 629.11kHz, (g) Mode g with f0 = 733.46kHz.
FIG. 7. Quality factors for several out-of-plane modes of design A plotted over the eigenfrequencies of the modes at 1 mbar and 25 • C. The green dots show the simulated quality factors based on gas damping.The red dots show the simulated quality factors based on TED, i.e.Eq. (16).The black dots show the resulting quality factors based on the gas damping and TED simulations.The blue dots show the measured quality factors.The solid lines indicate the trends of the data.The solid green line is a linear fit of the simulated gas damping quality factors over eigenfrequency.The solid red line was obtained as the median of the TED quality factors.The solid black line is the resulting quality factor due to the solid green and red lines.The solid blue line is equal to the solid black line with an additional constant Q anchor = 75000 included.The inset shows a magnified plot until 200 kHz without trend lines.