A two-port electrothermal model for suspended MEMS device structures with multiple inputs

Advances in micromachining have led to the development of microelectromechanical systems (MEMS) devices with suspended structures used in a variety of sensors. Of note for this work are sensor types where two elements exist on the suspended membrane, including examples like air flow and differential pressure detectors, gas detection, and differential scanning calorimetry sensors. Intuitively one would argue that some thermal loss exists between the two elements. However, surprisingly little is documented about this electrothermal interaction. The work presented here addresses this shortcoming by defining a new parameter set, a matrix of thermal coupling coefficients. They are used within our novel two-port electrothermal model based on the heat flow equation adapted as a linear system of equations. However, the model is only effective with knowledge of these coefficients. We introduce an approach to extract the coefficients using finite-element method (FEM)-based multiphysics simulation tools and revisit and extend our previous method of non-ideal power coupling, this time to extract the coefficient matrix from measured data. Both specialist simulation tools and device manufacturing are very expensive. However, they are the only choices in the absence of an analytic model. A major contribution of this work is the derivation of a model to predict the coefficients by analytic means from the device dimensions and material properties. The research contribution and paper culminate in a comparison of analytic, simulated, and experimentally extracted values of two different devices to verify and demonstrate the effectiveness of the proposed models. The values compare well and show that the best results achieved are approximately 90 % and 70 % thermal linkage respectively for vacuum and atmospheric pressure conditions.

1 Introduction

Electrothermal circuits (ETC) are circuits where the electrothermal interaction between components and the substrate is considered. Their analysis has been present in the literature since the early sixties (Matzen et al., 1964; Gray and Hamilton, 1971), but researchers have pursued it with less vigour than the purely electrical counterparts. This is mostly because of the large power consumption constraint associated with their operation (Louw et al., 1977). It was the commercialisation of MEMS manufacturing technologies that gave rise to a vast range of applications for ETC, notably because of the advances in micromachining. This led to high-yield, low-cost manufacturing and massive improvements in electrothermal performance. These applications have penetrated the industrial, military, and consumer markets, showing their diversity and the extensive marketing opportunities they present. Examples include the likes of on-chip sensors for diagnostics (Sisto et al., 2010; Zhang et al., 2009), chemical detection (Corsi et al., 2012; Barritault et al., 2013), thermal conductivity sensors (heat flow sensors and micro-hotplates) (Senesac et al., 2009; Arndt, 2002; Elmi et al., 2008), air flow sensors (Baltes et al., 1998; Johnson and Higashi, 1987; Simon et al., 2001), and thermal imaging microbolometers (Akula et al., 2011; Dulski et al., 2013; Castaldo et al., 1996; Liddiard, 2013).

The last two application areas are especially of interest in this work. Researchers have invested significant effort in improvements in the analytic modelling of the electrothermal effect, particularly for suspended membrane structures. Two recent introductions offer substantial improvements and are worth mentioning. First, joining a narrow and wide region results in an additional radial thermal conduction component. The restriction resistance model proposed is very effective and accounts for the additional component under vacuum conditions (Topaloglu et al., 2010). Second, an elliptical model better accounts for the radial temperature distribution on a suspended membrane in atmospheric pressure conditions (Schoeman and du Plessis, 2016). Both methods improve conventional models by up to 40 %.

Despite the continued active contributions in the modelling of traditional microbolometers, the detailed modelling of more generic dual-element devices, especially the interaction between these elements labelled as C in Fig. 1a, is lacking. This is surprising considering that these devices have existed for many years (Arndt, 2002; Yoon and Wise, 1994; Klaassen et al., 1995; Baltes et al., 1998). The approach in analytic electrothermal modelling of the detailed thermal interaction between the two devices or elements is never very rigorous and the interactions of the suspended elements with one another are mostly unclear. Authors typically assume perfect thermal linkage without mentioning the possibility of imperfect linkage. Furthermore, the electrical feedback network will often compensate for the imperfection. One technique offers an analytic solution to the temperature distribution of a thin-film resistor sensor pasted onto a large mass. However, the solutions are complex and time-consuming to solve. The author concludes that the model is useful for creating CAD tools (von Arx et al., 2000). This limits the practicality of the solutions.

It is interesting to note that the absence of an accurate analytical model for the electrothermal interaction can likely be attributed to the strength of the FEM simulation methodologies. However, these tools are very expensive and mostly used to simulate device-level behaviour, although some tools allow the capability of extracting the device model to a circuit-level model that can be used in SPICE-based simulation, but at an additional cost. As such, simpler SPICE models have been developed to couple the electrothermal parameters to electrical ones used in readout circuits for infrared microemitters and microbolometers (Kiran and Karunasiri, 1999; Nazdrowicz et al., 2015; Kim and Ko, 2015). However, we should highlight that the lack of an existing model hampers (i) the understanding and design of optimal thermal coupling devices and (ii) the design of effective readout circuits for ETC while accounting for thermal losses. This leaves a research gap that this paper will address. In Sect. 2 of this work we present a novel method to transform complex suspended device structures with two resistive elements into electrothermal equivalent microbeams that can be used to determine a matrix of thermal coupling coefficients. This is derived for both vacuum and atmospheric pressure conditions. We also propose that a multidimensional heat flow problem can be simplified and described using this set of coefficients within a linear system of equations modelling an equivalent electrothermal two-port system. Section 3 introduces a new simulation methodology to extract the matrix of coefficients. We also extend our experimental method introduced previously (Schoeman and du Plessis, 2012) for Ti1 at atmospheric pressure to account for the complete set of parameters required by the latest two-port model. In Sect. 4 we present and compare analytic, simulated, and experimental results for two different devices at two different pressure levels. Finally, we present some conclusions based on the results.

Figure 1 The operating principles and the thermal model for a generic dual-resistive element device. (a) The device topology highlighting the required components of the device structure and the thermal linkage that exists between the resistive elements. (b) An equivalent heat flow diagram for a generic two-element device. The elements are labelled as the heater and detector, although each element can serve as both.

2 Theory and modelling of dual-element suspended plate structures

2.1 Operating principles of suspended devices

The device structure is similar to conventional microbolometers that are created by stacking multiple layers and removing a sacrificial layer by wet etching. The structure can be divided into two regions, (i) a membrane that is suspended over a cavity after etching and (ii) support beams to hold the membrane in position. This is illustrated in Fig. 1a. The suspended device has a thermal capacity that stores absorbed energy over time and acts as a heat reservoir. It can be heated by applying electrical power to a resistive material deposited on top of the suspended membrane. However, the resulting temperature gradient between the membrane and substrate will cause heat to flow from the warmer membrane to the cooler substrate through the support beams and the gas filling the cavity underneath the membrane. These thermal conductive components are referred to as the beam thermal conduction and the gaseous thermal conduction respectively. These interactions between the various components are illustrated in Fig. 1b and, with the exception of the thermal linkage components, are described by a heat flow model for conventional microbolometers blind to IR radiation stated mathematically as (Shie et al., 1996; Galeazzi and McCammon, 2003)

$$\begin{array}{}\text{(1)}& H{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{m}}}{\mathrm{d}t}}+\underset{G_{\mathrm{eff}}\mathrm{\Delta }T}{\underbrace{G_{\mathrm{atm}}\left(T_{\mathrm{m}}-T_{\mathrm{a}}\right)+G_{\mathrm{leg}}\left(T_{\mathrm{m}}-T_{\mathrm{sub}}\right)}}=P,\end{array}$$

where H is the thermal capacity of the detector with unit J K⁻¹, G_leg is the thermal conductance of the detector to the substrate via the supporting legs in W K⁻¹ analogous to the thermal resistance that results in losses in heat, G_atm is the thermal conductance of the detector through the atmosphere to the substrate with unit W K⁻¹, T_m is the temperature of the membrane in Kelvin, T_sub is the temperature of the substrate in Kelvin, T_a is the ambient temperature surrounding the membrane in Kelvin, and P is the electrical power applied to bias the resistive element in Watt. The equation is often simplified by assuming that (i) the substrate is at room temperature (T_a=T_sub) and/or (ii) the device operates at steady state (G_effΔT=P). An additional assumption is that G_atm=0 under vacuum conditions, simplifying G_eff even further.

G_eff consists of two components (Eriksson et al., 1997), one of which is a thermal conductance component attributed to heat loss between the suspended membrane and the substrate via the support beams referred to as the beam thermal conductance, G_leg. This is the dominant heat loss mechanism under vacuum conditions and is given by (Topaloglu et al., 2010; Niklaus et al., 2007)

$$\begin{array}{}\text{(2)}& G_{\mathrm{leg}}=N\sum _{i=\mathrm{1}}^n{k}_{\mathrm{b},i}{\displaystyle \frac{W_{\mathrm{B},i}{d}_{\mathrm{th},i}}{L_{\mathrm{B},i}}},\end{array}$$

with i an index to a specific layer, N the number of support beams (typically two or four), k_b the effective thermal conductivity of the beam measured in W m⁻¹ K⁻¹, W_B the width of the beam measured in metres, d_th the thickness of the beam measured in metres, and L_B the length of the beam measured in metres. Accounting for the radial thermal conductive component can significantly improve the estimation of the thermal conductivity. The transformed thermal conductivity is given by (Topaloglu et al., 2007)

$$\begin{array}{}\text{(3)}& k_{\mathrm{b}}{}^{\prime }={\displaystyle \frac{\mathit{\pi }{k}_{\mathrm{b}}{L}_{\mathrm{M}}}{\mathit{\pi }{L}_{\mathrm{M}}+\mathrm{2}{W}_{\mathrm{M}}\ln \left(\frac{\mathit{\pi }{W}_{\mathrm{M}}}{\mathrm{2}{W}_{\mathrm{B}}}\right)}}={\displaystyle \frac{k_{\mathrm{b}}}{\mathrm{1}+\frac{\mathrm{2}{W}_{\mathrm{M}}}{\mathit{\pi }{L}_{\mathrm{M}}}\ln \left(\frac{\mathit{\pi }{W}_{\mathrm{M}}}{\mathrm{2}{W}_{\mathrm{B}}}\right)}},\end{array}$$

where L_M and W_M are the membrane length and width respectively.

The second component contributing to G_eff is the atmospheric or gaseous thermal conductance, G_atm. This is the dominant heat loss mechanism at atmospheric pressure and is given as (Niklaus et al., 2007)

$$\begin{array}{}\text{(4)}& G_{\mathrm{atm}}=k_{\mathrm{air}}{A}_{\mathrm{eff}}/d_{\mathit{\lambda }},\end{array}$$

where $A_{\mathrm{eff}}=L_{\mathrm{M}}\times W_{\mathrm{M}}$ is the area of the membrane in metres squared used in conventional methods, k_air is the thermal conductivity of air at room temperature and pressure (0.00284 W m⁻¹ K⁻¹), and d_λ is the cavity depth measured in metres.

Up to this stage, we have only considered the modelling of a conventional microbolometer device. This is a device with two terminals and a single port where power can be applied. Therefore, thermal linkage does not exist. It would be natural to wonder what the impact on the heat flow equation would be when we introduce a second (or multiple) resistive element(s) to the suspended device. This would allow for an additional port to the device through which we can control the plate temperature. The steady-state equivalent of Eq. (1) can be changed to allow for the perspective from the kth input as

$$\begin{array}{}\text{(5)}& \underset{G_{\mathrm{eff}}\mathrm{\Delta }T}{\underbrace{\left(G_{\mathrm{atm}}+G_{\mathrm{leg}}\right)\left(T_{\mathrm{m},k}-T_{\mathrm{sub}}\right)}}=P_k,\end{array}$$

where P_k is the biasing electrical power of the k_th detector element in Watt and T_m,k is the average temperature increase of the membrane in Kelvin resulting from the applied power component P_k. Furthermore, we should be able to observe a portion of the power applied to one port on every other port of the device, since the resistive elements all share the same thermal reservoir. This suggests that the elements are thermally linked to each other, and we introduce a parameter, β_i,j, to describe this thermal linkage. Since each resistive element will contribute towards a portion of the total change in the average membrane temperature, we propose to model the effective or total thermal increase as a system of linear equations, best described in vector form as

$$\begin{array}{}\text{(6)}& G_{\mathrm{eff}}\cdot \mathbf{\Delta }\stackrel{\mathrm{‾}}{T}=\stackrel{\mathrm{‾}}{\mathit{\beta }}\cdot \mathit{P},\end{array}$$

where $\mathbf{\Delta }\stackrel{\mathrm{‾}}{T}$ is the average temperature increase of the membrane consisting of incremental temperature increases resulting from the biasing of the respective resistive elements, P is a power vector consisting of the power applied to each individual resistor, and $\stackrel{\mathrm{‾}}{\mathit{\beta }}$ is the matrix of weighted directional thermal coupling coefficients that exist between the various resistive elements. Here it is assumed that G_eff is a constant scalar value for the developed structure which holds true for symmetric designs.

Therefore, the heat balance equation for a dual-element device in the steady state can be given by

$$\begin{array}{}\text{(7)}& G_{\mathrm{eff}}\cdot \left[\begin{array}{c}\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}\\ \mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}\end{array}\right]=\left[\begin{array}{cc}\stackrel{\mathrm{‾}}{{\mathit{\beta }}_{\mathrm{1},\mathrm{1}}}& \stackrel{\mathrm{‾}}{{\mathit{\beta }}_{\mathrm{1},\mathrm{2}}}\\ \stackrel{\mathrm{‾}}{{\mathit{\beta }}_{\mathrm{2},\mathrm{1}}}& \stackrel{\mathrm{‾}}{{\mathit{\beta }}_{\mathrm{2},\mathrm{2}}}\end{array}\right]\cdot \left[\begin{array}{c}{P}_{\mathrm{1}}\\ P_{\mathrm{2}}\end{array}\right],\end{array}$$

from which we can easily solve the respective individual element average temperature increases if $\stackrel{\mathrm{‾}}{\mathit{\beta }}$ is known. The latter is relatively simple to determine from simulation and experimentally, as shown in Sect. 3. However, the analytic approach proves more challenging, as elaborated on next.

2.2 An analytic solution for the thermal coupling coefficients

Researchers have published comprehensive techniques for modelling the electrothermal behaviour of device structures with various degrees of similarity to our own. Examples include temperature profiling by solving the general heat flow partial differential equation (PDE) first developed and solved by Fourier, reducing dimensionality for specific applications by transformation, applying lumped thermal element analysis, and considering thermal coupling in the context of Seeback and Peltier effects. As expected, the intent is often to find one or more mechanical properties. We use similar methods that culminate in a robust approximate model to describe the electrothermal effect. The novel model consists of three parts: (i) finding the thermal coupling coefficients for two points on a microbeam, (ii) transforming between the microbeam and a more complex structure, and (iii) increasing the number of points by averaging over the entire resistive length. Each part is elaborated on next in a separate subsection.

2.2.1 Determine the coupling coefficient between two points on a microbeam

A number of examples exist where the steady-state heat transfer, well known to be governed by a second-order PDE, is solved analytically to find the temperature profile of a microbeam, cantilever, or similar rectangular MEMS type structures, which is then in turn further developed to solve some mechanical property, be it displacement, deflection, dampening, stresses, or elasticity (Chow and Lai, 2009; Kokkas, 1974; Huang and Lee, 1999; Yan et al., 2004; Mankame and Ananthasuresh, 2001). We use a similar approach to derive an expression for the newly proposed thermal coupling coefficients.

Figure 2 A microbeam with two resistive strips and the relevant dimensions, along with the expected temperature profiles over the length of the bar. (a) The temperature profile through a microbeam when heat flow occurs only through the microbeam is typical under vacuum conditions. (b) The temperature profile through a microbeam when heat flow occurs through both the solid components of the microbeam, as well as the medium (either a gas or a fluid) between the microbeam and a heat sink.

First, consider a microbeam with two resistive components as illustrated in Fig. 2a. R₁ is heated to a specific temperature T_m,1 that results in two heat flow vectors. One of these, q₂, flows through R₂ that experiences an increase in temperature in response. We know that under vacuum conditions the simplest form of a second-order linear homogeneous differential equation describes heat diffusion, given by $\frac{{\partial }^{\mathrm{2}}T}{\partial x^{\mathrm{2}}}=\mathrm{0}$. The differential equation is solved by double integration and a particular solution is found from the general solution. The former is given by

$$\begin{array}{}\text{(8)}& T\left(x\right)=T_{\mathrm{m},\mathrm{1}}-\left({\displaystyle \frac{T_{\mathrm{m},\mathrm{1}}-T_{\mathrm{sub}}}{L_{\mathrm{2}}}}\right)x.\end{array}$$

We define a new parameter, the thermal coupling coefficient β_2,1, as a ratio of temperature differences. Substitution of the parameters in Fig. 2a into Eq. (8) yields

$$\begin{array}{}\text{(9)}& {\mathit{\beta }}_{\mathrm{2},\mathrm{1}}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{2}}}{\mathrm{\Delta }{T}_{\mathrm{1}}}}={\displaystyle \frac{T_{\mathrm{m},\mathrm{2}}-T_{\mathrm{sub}}}{T_{\mathrm{m},\mathrm{1}}-T_{\mathrm{sub}}}}={\displaystyle \frac{L_{\mathrm{1}}}{d+L_{\mathrm{1}}}}.\end{array}$$

The problem is more complex when the device operates at atmospheric pressure since an additional gaseous thermal conduction component exists. This is illustrated in Fig. 2b. The best approach to derive a solution is to model the problem as a distributed thermal resistance network with r_s and g_air the unit components for the solid thermal resistance and the gaseous thermal conduction respectively. The temperature distribution is found by solving the standard form of a second-order linear homogeneous differential equation if the system is in steady state. This is given by

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{d^{\mathrm{2}}T\left(x\right)}{\mathrm{d}{x}^{\mathrm{2}}}}-K\cdot T\left(x\right)=\mathrm{0},\end{array}$$

with K=r_sg_air. The roots of the quadratic equation may be solved and, after defining boundary conditions, the particular solution can be shown to be given by

$$\begin{array}{}\text{(11)}& T\left(x\right)=\left(T_{\mathrm{m},\mathrm{1}}-T_{\mathrm{sub}}\right)e^{-\sqrt{K}\cdot x}.\end{array}$$

It is insightful to consider the parameter L_ch, the characteristic thermal length governing the exponential decay observed for the temperature over distance, given by

$$\begin{array}{}\text{(12)}& L_{\mathrm{ch}}={\displaystyle \frac{\mathrm{1}}{\sqrt{K}}}={\displaystyle \frac{\mathrm{1}}{\sqrt{r_s{g}_{\mathrm{air}}}}}=\sqrt{{\displaystyle \frac{d_{\mathit{\lambda }}}{k_{\mathrm{air}}}}\sum _i{k}_{\mathrm{b},i}{d}_{\mathrm{th},i}}\end{array}$$

for a multilayered device structure and approximately 10 µm for our devices. The thermal coupling coefficient follows from the definition earlier in Eq. (9) and is given by

$$\begin{array}{}\text{(13)}& {\mathit{\beta }}_{\mathrm{2},\mathrm{1}}={\displaystyle \frac{T\left(L_{\mathrm{1}}\right)}{T\left(L_{\mathrm{2}}\right)}}=e^{-\sqrt{K}\cdot d}=e^{-d/L_{\mathrm{ch}}}.\end{array}$$

It should be clear that the coupling coefficient depends on the ratio of the separation distance between the resistive elements and the characteristic length of the microbeam in an exponential manner. The thermal coupling can be maximised by reducing the separation distance of the resistive elements in a design.

2.2.2 Transform between microbeams and suspended plate structures

It is also desirable to calculate thermal linkage for more complex device geometries. This should be possible if the more complex geometry can be transformed into a simpler microbeam and the thermal coupling coefficients of the previous section is applied. In general, transformations between device geometries are not new and examples include multidimensional structure transformation into cylindrical coordinates (Khan and Falconi, 2013; Sberveglieri and Hellmich, 1997; Simon et al., 2001). Our approach is quite different and is based on an equivalent heat flow vector analysis.

Figure 3 A graphical illustration of the heat flow vector of interest that results in thermal linkage for a dual-element suspended plate with four support beams. A simulation model with a simplified elemental heat source supports the proposal. (a) A structural overview of the device with the relevant dimensions indicated. The proposed vectors are indicated by the orange arrows, while a single heat source element is indicated as a red square and the target element on the opposing resistor is indicated by a blue square. (b) A simulation result showing the heat flow vectors from a single source element as well as the proposed heat flow vector (indicated by the orange arrow) resulting in the thermal linkage between the heat source and target.

We propose that the equivalent heat flow of interest from a specific heat source point through a specific target element will occur only along a single vector as illustrated in Fig. 3b. This should be intuitive, but can also be observed from the supporting FEM simulation result of Fig. 3b. The simulation result also shows an infinite amount of additional heat flow vectors, but recall that the desired solution is the single vector flowing through the target element. Therefore, we only need to consider the vector V₂.

Closer inspection reveals that V₂ extends from the heat source towards the closest point on the second resistive element where thermal linkage occurs, after which the vector extends by some fraction C_m past R₂. It changes direction towards the connection between the support beam and the membrane and finally extends along the support beam towards the substrate. V₂ is stated mathematically as

$$\begin{array}{}\text{(14)}& {\mathit{V}}_{\mathrm{2}}=\mathit{d}+\mathit{C}+\mathit{D}+{\mathit{L}}_{\mathrm{B}}\end{array}$$

and the distance between the heat source point and the substrate is calculated by

$$\begin{array}{}\text{(15)}& \begin{array}{rl}∥{\mathit{V}}_{\mathrm{2}}∥& =d+∥〈C_x,C_y〉∥+∥〈D_x,D_y〉∥\\ & +∥〈L_{B,x},L_{B,y}〉∥\\ & =d+C_m\left({\displaystyle \frac{W_{\mathrm{M}}}{\mathrm{2}}}-{\displaystyle \frac{d}{\mathrm{2}}}\right)+\\ & \sqrt{{\left(\left(\mathrm{1}-C_m\right)\left({\displaystyle \frac{W_{\mathrm{M}}}{\mathrm{2}}}-{\displaystyle \frac{d}{\mathrm{2}}}\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{L_{\mathrm{M}}}{\mathrm{2}}}\right)}^{\mathrm{2}}}+L_{\mathrm{B}}\\ & =L_{\mathrm{2}},\end{array}\end{array}$$

where C_m is a proportionality constant (a typical value is 0.5) that describes the depth of penetration of C in the x direction before V₂ changes direction, W_M is the width of the membrane, L_B is the length of the support beam, and d is the separation distance between the source and target points. It should become immediately clear after inspection of Eq. (15) that this vector is in the form $L_{\mathrm{2}}=d+L_{\mathrm{1}}$. This can be used to relate a suspended plate structure, like a microbolometer, to the microbeam of the previous section and allows for the calculation of L₁, a microbeam parameter, as a function of the plate dimensions.

A complete transformation between the two structures is possible if the effective thermal resistance experienced by V₂ is the same as that of the equivalent microbeam. This occurs if the microbeam has an effective width given by

$$\begin{array}{}\text{(16)}& W_{\mathrm{eff}}={\displaystyle \frac{W_{\mathrm{M}}\left(\mathrm{\Delta }{x}_{\mathrm{1}}+\mathrm{\Delta }{x}_{\mathrm{2}}\right)}{\mathrm{\Delta }{x}_{\mathrm{1}}+\mathrm{\Delta }{x}_{\mathrm{2}}\left(\frac{W_{\mathrm{M}}}{W_{\mathrm{B}}}\right)}},\end{array}$$

where $\mathrm{\Delta }{x}_{\mathrm{1}}=\frac{L_{\mathrm{M}}-d}{\mathrm{2}}+d=\frac{L_{\mathrm{M}}+d}{\mathrm{2}}$ and Δx₂=L_B.

2.2.3 Determine an average coupling coefficient over the entire resistor length

The model components developed thus far for thermal linkage between two single points are not very useful until the concepts are extended to apply to heaters and detectors with more practical lengths. As such, we introduce the notion of dividing the heater and detector elements into unique regions and applying an averaging mechanism. To classify a region as unique, we identify the parts of the heater that have a constant separation distance to the detector. Every such length with a unique separation distance would make up a unique region. A number of k such regions can exist for a device, but the thermal linkage is optimised by reducing the number of regions, with k=1 being optimal. We illustrate the concept of regions in Fig. 4 and show examples in Sect. 2.3. We can use the entire length of a region as part of the weighting coefficient, $l_{\mathrm{reg},i}/L_{R\mathrm{1}}$.

Furthermore, a length of the heater element might be associated with more than one region. This increases computational complexity, but when considering Eqs. (9) and (13) we realise that the increase in separation distance dramatically reduces the thermal linkage. We can simplify the problem and find an approximate solution by only considering the region with the smallest separation distance for any length of the heater. Therefore, a region only needs to be coupled to its closest neighbouring region on the alternate resistor, with the remaining regions considered insignificant. These considerations result in an averaging algorithm given by

$$\begin{array}{}\text{(17)}& {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{2},\mathrm{1}}\approx \sum _{i=\mathrm{1}}^k{\displaystyle \frac{l_{\mathrm{reg},i}}{L_{R\mathrm{1}}}}{\mathit{\beta }}_{\mathrm{2},\mathrm{1}}\left(d_{i,j_{\mathrm{n}}}\right),\end{array}$$

with j_n being the single nearest neighbour only of the particular ith region. This results in a computational reduction, often to the extent that the thermal coupling can be computed quickly by hand without the need for simulation. An example is provided at the end of the next section.

2.3 Device designs

The material properties of the materials used to construct the suspended devices play a significant role in the device performance parameters. Notably, the thermal conductivity, k_i, of the materials used to construct the support beams plays a role in the thermal conductance, G_leg, especially that of silicon nitride in typical suspended plate designs. The various material densities, ρ_i, and specific heat, c_i, play a major role in the thermal time response, τ, and the thermal capacity, H, of the devices. The relevant properties are summarised in Table 1.

Table 1 A summary of the most important material properties that are used in the manufacturing of the test devices.

The manufactured Ti thin-film dual resistive element devices are shown in Fig. 4 with the dimensions relevant to Eq. (17) superimposed on the images. Earlier, the plate dimensions were defined and shown in Fig. 3a. Take note of the additional parameters: (i) the length of the first Ti thin-film resistor, L_R,1, (ii) the width of the first resistor, W_R,1, and (iii) the length of the coupling region between the resistive elements, d_i,j. The length and the width of the second Ti thin-film resistor are not indicated, but they are defined in the same way as for the first resistor. A summary of the device designs and their dimensions and surface characteristics are provided in Table 2, while the linkage parameters are compared in Table 3.

Figure 4 SEM micrographs of the various manufactured test structures with the relevant thermal linkage parameters indicated. Regions are identified and defined as parts of the design that have a constant separation distance d_i over the entire resistance length of that region. (a) A suspended plate structure (Ti1) with optimal thermal linkage. (b) A suspended plate structure (Ti2) with high thermal linkage.

Table 2 A summary of the Ti thin-film dual-element device dimensions.

Table 3 A summary of the device parameters required by the proposed analytic model to determine the thermal linkage.

A device (Ti1) optimised for a very high thermal coupling coefficient is illustrated in Fig. 4a. The device features two Ti thin-film resistors with a separation distance d_A=3 µm over the entire length of the resistors $L_{R,\mathrm{1}}=L_{\mathrm{A}}=\mathrm{460}$ µm. A second device, Ti2, offers a high thermal coupling coefficient and is illustrated in Fig. 4b with the device dimensions specified in Table 2. The device has a resistance length $L_{R,\mathrm{1}}=\mathrm{210}$ µm and is divided into three unique regions (A, B, and C) with region lengths L_A=155, L_B=38, and L_C=8.5 µm and separation distances d_A=3, d_B=9, and d_C=15 µm.

3 Methodology

3.1 Approach followed in the simulation study

Multiphysics simulation platforms are highly effective in simulating 3-D device structures. We opted to use Coventorware, an industry leading multiphysics MEMS simulator, for the study. The standard material properties database of the simulation software is used and the values correspond closely to those in Table 1. A thermomechanical solver is selected to simulate thermoelectric physics in the steady state by applying an electrical potential to the gold contacts of the device as a boundary condition and then to investigate the resulting thermal effects. The substrate is considered a thermal heat sink with constant temperature, which can be applied as a fixed boundary condition where the substrate temperature is selected as 300 K. During the steady-state analysis, the input potentials of both resistors are swept independently over their respective voltage ranges and the resulting average temperature changes of the resistive elements are noted.

The parameters of the heat balance equation of Eq. (6) are solved in two parts. First, recall that each value in the coupling coefficient matrix, $\stackrel{\mathrm{‾}}{\mathit{\beta }}$, is in the form of the proposed definition of Eq. (9). These coefficients are solved by applying a controlled average temperature stimulus to only one element and then extracting the average element temperature increase of the alternate element to find the ratios as

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& \mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}={\left.{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{1},\mathrm{2}}{T}_{\mathrm{m},\mathrm{2}}\right|}_{T_{\mathrm{m},\mathrm{1}}=\mathrm{0}}\Rightarrow {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{1},\mathrm{2}}={\displaystyle \frac{\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}}{T_{\mathrm{m},\mathrm{2}}}},\\ & \mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}={\left.{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{2},\mathrm{1}}{T}_{\mathrm{m},\mathrm{1}}\right|}_{T_{\mathrm{m},\mathrm{2}}=\mathrm{0}}\Rightarrow {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{2},\mathrm{1}}={\displaystyle \frac{\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}}{T_{\mathrm{m},\mathrm{1}}}}.\end{array}\end{array}$$

Secondly, once $\stackrel{\mathrm{‾}}{\mathit{\beta }}$ is known, the effective thermal conductance is solved from Eq. (7) when rewritten as

$$\begin{array}{}\text{(19)}& \left[\begin{array}{c}{G}_{\mathrm{eff},\mathrm{1}}\\ G_{\mathrm{eff},\mathrm{2}}\end{array}\right]=\left[\begin{array}{cc}{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{1},\mathrm{1}}& {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{1},\mathrm{2}}\\ {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{2},\mathrm{1}}& {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{2},\mathrm{2}}\end{array}\right]\cdot \left[\begin{array}{c}{P}_{\mathrm{1}}\\ P_{\mathrm{2}}\end{array}\right]\cdot {\left[\begin{array}{c}\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}\\ \mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}\end{array}\right]}^{-\mathrm{1}},\end{array}$$

where $G_{\mathrm{eff}}=G_{\mathrm{eff},\mathrm{1}}=G_{\mathrm{eff},\mathrm{2}}$ for symmetric device designs.

3.2 Approach followed for the experimental method

The experimental method requires three parts to solve the parameters of the heat balance equation of Eq. (6). First, we use a popular purely electrical method that is widely used in industry for microbolometers to solve the temperature coefficient of resistance (TCR), α, in a controlled oven with variable temperature (Shie et al., 1996). Second, the device voltage is measured when sweeping the input current using a parameter analyser. The device resistance can then be calculated from the measured results. In turn, G_eff may be extracted by finding the slope of the inverse resistance as a function of the square of the biasing current that is given by

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{1}}{R_{Bi}}}={\displaystyle \frac{\mathrm{1}}{R_{Bi,\mathrm{0}}}}-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}}{G_{\mathrm{eff}}}}{I}_B^{\mathrm{2}}\end{array}$$

as a linear equation with $i=\mathrm{1},\mathrm{2}$ the different ports of the device.

The third and final step is to notice from Eq. (7) that the average temperature increase of one element can be solved by applying a controlled power stimulus to the alternate resistive element while ensuring that no power is applied to the element under investigation, from which the set of equations reduces to provide the thermal coupling coefficients as

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& \mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}={\left.\mathrm{1}/G_{\mathrm{eff}}\left(\stackrel{\mathrm{‾}}{{\mathit{\beta }}_{\mathrm{1},\mathrm{2}}{P}_{\mathrm{2}}}\right)\right|}_{P_{\mathrm{1}}=\mathrm{0}}\Rightarrow {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{1},\mathrm{2}}=\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}{G}_{\mathrm{eff}}/P_{\mathrm{2}},\\ & \mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}={\left.\mathrm{1}/G_{\mathrm{eff}}\left({\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{2},\mathrm{1}}{P}_{\mathrm{1}}\right)\right|}_{P_{\mathrm{2}}=\mathrm{0}}\Rightarrow {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{2},\mathrm{1}}=\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}{G}_{\mathrm{eff}}/P_{\mathrm{1}},\end{array}\end{array}$$

which is equivalent to

$$\begin{array}{}\text{(22)}& \begin{array}{r}{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{1},\mathrm{2}}={\displaystyle \frac{\mathrm{\Delta }{R}_{B\mathrm{1}}{G}_{\mathrm{eff}}}{R_{B\mathrm{1},\mathrm{0}}\mathit{\alpha }{P}_{\mathrm{2}}}},\\ {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{2},\mathrm{1}}={\displaystyle \frac{\mathrm{\Delta }{R}_{B\mathrm{2}}{G}_{\mathrm{eff}}}{R_{B\mathrm{2},\mathrm{0}}\mathit{\alpha }{P}_{\mathrm{1}}}}.\end{array}\end{array}$$

This implies that the resistance of the first element can be increased by a value equivalent to $\mathit{\alpha }{R}_{B\mathrm{1},\mathrm{0}}{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{12}}{P}_{\mathrm{2}}/G_{\mathrm{eff}}$ compared to the reference resistance. The modulation of the resistance R_B1 is controlled by P₂ and the effect is illustrated in Fig. 5. The figure also serves to further clarify the differences in the contributions of T_m,1 (or P₁) and T_m,2 (or P₂). Although not indicated in the figure, the same modulation concept will hold true for the resistance of the second element, $R_{B\mathrm{2}}\left(\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}\right)$, by applying a control signal to the first element.

Figure 5 Effect on the resistance in terms of the current stimuli. The bottom graph represents the reference curve when biasing the first resistive element alone, while the top graph shows the additional increased resistance when the second resistor is also biased.

4 Results and discussion

The analytically predicted results of the devices are provided in Table 3. However, to verify the validity of the theoretical model and reach a meaningful conclusion, those results need to be compared to results of simulation and experimental methods that extract the coupling coefficients. As such, the simulation and experimental results obtained for the two thermally coupled MEMS devices are presented next for both vacuum and atmospheric pressure conditions.

Table 4 A comparison of the thermal beam conductance of the Ti thin-film dual-element devices under vacuum conditions.

^a Calculated using the standard approach with Eq. (2). ^b Calculated using the compensated method with Eq. (3).

Table 5 A comparison of the total thermal conductance of the devices at atmospheric pressure.

^a Calculated using the standard approach with Eq. (2) and Eq. (4). ^b Calculation includes the effect of Eq. (3).

Table 6 A comparison of the average thermal coupling coefficient, ${\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{12}}$, of the devices under vacuum conditions.

Table 7 A comparison of the average thermal coupling coefficient, ${\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{\mathrm{12}}$, of the dual-element devices at atmospheric pressure.

4.1 Simulation results of the thin-film metal devices

The first results presented are of a simulated parametric voltage sweep applied to the gold contacts of R_B1, ranging from 0 to 0.25 V, while R_B2 remains unbiased. This offers a benchmark result similar to what is measured for conventional microbolometers. Once the benchmark curve is obtained, the potential applied across R_B2 is increased in 0.05 V increments up to 0.25 V. These increments are indicated with labels A through E in Fig. 6. At each increment, the original potential sweep applied to R_B1 is repeated. The average temperatures of the thin-film Ti resistors are extracted at each point of the multivariable sweep, with the results for the first element plotted in Fig. 6c. These results are also used to extract the electrical power dissipated by the device. These are combined as the P–T curve illustrated in Fig. 6d that serves to extract the beam thermal conductances from the slopes. Note that these are intermediate results. Therefore, only the results of Ti1 under vacuum conditions are presented for brevity. The same method is followed for atmospheric pressure conditions and repeated for device Ti2 under both conditions.

Figure 6 Simulation results of the thermoelectric interaction that exists for the Ti devices under vacuum conditions. (a) XZ projection of the thermal profile of R_B1 for Ti1. (b) YZ projection of the thermal profile of R_B1 for Ti1. (c) 3-D thermal profile of R_B1 for Ti1. (d) Temperature curves at different power levels to extract the beam thermal conductance.

Figure 7 Experimental results of a parametric current sweep under vacuum conditions for the Ti thin-film devices indicating the resulting resistance of the first resistive element, R_B1, while R_B2 is biased at different values as indicated in the legend entries. These results are also used to extract the beam thermal conductance. Note that the legends of (b) also apply to (a) but that Ti1 is presented at the top of (a) and at the bottom of (b), as expected when comparing R to 1∕R. (a) Experimental R–I curves. (b) Experimental 1∕R–I² curves.

4.2 Measurement results of the thin-film metal devices

The V–I curves are useful to calculate the power dissipation of the devices, since G_eff is extracted from this. The R–I curves also provide additional graphical insight into the resistive modulation effect. The next results show the experimental measurements of both devices in Fig. 7a. The sensed element voltage is measured when applying an input current sweep under vacuum conditions with the Hewlett Packard 4155B parameter analyser. The input current I_B1 is swept up to 500 µA for the devices within the evacuated dewar as long as a maximum power consumption limit of 100 µW is met. The lower curves serve as the reference since $I_{B\mathrm{2},\mathrm{0}}=\mathrm{0}$ µA. Once the reference curves have been obtained, the experiment is repeated at increments of approximately $I_{B\mathrm{2},i}=\mathrm{50}$ µA and $I_{B\mathrm{2},i}=\mathrm{62.5}$ µA respectively for devices Ti1 and Ti2. A very high thermal increase is possible even at moderate power consumption levels for Ti1 because of the low beam thermal conductance. This explains the need to limit the power consumption of the experiment. The slope of the curves suggests a positive TCR as expected for metal devices. The 1∕R–I² curves are also indicated and used to extract the effective thermal conductance using Eq. (20). The results are plotted in Fig. 7b.

4.3 Comparison and discussion of the results

A comparison of the effective thermal conductances for the devices is provided in Table 4. The theoretically derived values using both Eqs. (2) and (3) are compared to the simulation results extracted from Fig. 6d and the experimental results extracted from Fig. 7b. The results compare very well overall, but the experimentally obtained value for Ti2 is slightly smaller than anticipated. This is due to mask alignment difficulties encountered during manufacturing, as seen in the SEM micrograph in Fig. 4b. This led to the left support beams playing a smaller role in the thermal isolation of the membrane than expected. We do notice a large discrepancy between the analytic and simulated results for device Ti2 at atmospheric pressure conditions. This warranted a deeper investigation into the particular design. It revealed that the entire plate of the simulated device does not reach the average temperature of the heater, which means A_eff is much smaller than the value used in the prediction. This leads to a smaller simulated thermal conductance. This is similar to an effect we first observed in conventional microbolometer designs (Schoeman and du Plessis, 2016). Also notice that the experimental values are smaller than those simulated. The simulation is set up with d_λ=2 µm. However, microscopic inspection revealed that the devices are bending upwards. This was confirmed by a profilometric analysis, and it was found that Ti1 is bending upwards with d_λ in excess of 3 µm and Ti2, with its shorter supports, measuring d_λ=2.8 µm. This reduces the thermal conductance considerably, as expected from Eq. (4).

We use the methods of Sect. 3.1 and Sect. 3.2 to determine the thermal linkage between the resistive elements. Three biasing points, [0, 20, 100] % of the maximum I_B1, are selected on each of the I_B2,2–4 curves plotted in Fig. 7a. This results in nine data points for each device. We can use each result to extract a thermal coupling coefficient with Eq. (22) at that biasing point. The nine results are then averaged. The analytic, simulation, and experimental results obtained under vacuum conditions are presented in Table 6, while the same results at atmospheric pressure are presented in Table 7. As can be seen from the tabulated average thermal coupling coefficients, the predicted results of the Ti thin-film devices compare very well with the simulation and experimental results. The largest observed errors relative to the developed model are 7.6 % and 16.1 % respectively for the simulation and experimental results.

5 Conclusions

In this work, an analytic method is presented to quantify the thermal linkage expected between two resistive elements that share the suspended membrane of a MEMS device. This is a crucial component currently lacking in the complete modelling of the electrothermal interaction of such devices. We have confidence in the intermediate results, since the simulated and experimentally obtained thermal conductance values achieved are reasonable and compare well with what has been reported in the literature for conventional microbolometers. These devices compare well with our suspended structures. However, our thermal conductance values are larger within reason since our beam widths are increased to allow for the additional electrical connections of the second resistor. We have demonstrated, with both simulation and experimental measurement using manufactured devices, that the proposed model adequately predicts the thermal coupling effect. The main reason for the success of the model lies in the definition of Eq. (17). The best results are achieved and performance is optimised for coupling over very long lengths of resistors while maintaining a very small separation distance between them. This is typically achieved with long and narrow resistors, or L_R≫W_R, and hence, thin-film resistive materials are very suitable. Ti-based devices, for example, are ideal and outperform alternative material choices like VO_x from our experience. We conclude that the developed model offers a reasonable output that will prove very useful in future MEMS designs of applications making use of thermal coupling. Other than the applications mentioned in the introduction, this effect holds potential in the low-frequency domain for solutions of MEMS-based on-chip electrothermal isolators and even signal processing. Designers using alternative design simulation methods, for example at a circuit level using SPICE to simulate ETC, also stand to benefit from the proposed analytic model, since it allows for a simple method to calculate the thermal linkage between elements without the need to use expensive device-level FEM software. This will be the focus of future work. Even our best thermal linkage results of approximately 90 % and 70 %, depending on the pressure conditions, are far from 100 %. As such, designers should be cautious to assume ideal thermal linkage, which is often assumed in previous work, and consider incorporating thermal coupling coefficients into their ETC designs.