Dataset of the Optimization of a Low Power Chemoresistive Gas Sensor: Predictive Thermal Modelling and Mechanical Failure Analysis

Abstract

Over the last few years, employment of the standard silicon microfabrication techniques for the gas sensor technology has allowed for the development of ever-small, low-cost, and low-power consumption devices. Specifically, the development of silicon microheaters (MHs) has become well established to produce MOS gas sensors. Therefore, the development of predictive models that help to define a priori the optimal design and layout of the device have become crucial, in order to achieve both low power consumption and high mechanical stability. In this research dataset, we present the experimental data collected to develop a specific and useful predictive thermal-mechanical model for high performing silicon MHs. To this aim, three MH layouts over three different membrane sizes were developed by using the standard silicon microfabrication process. Thermal and mechanical performances of the produced devices were experimentally evaluated, by using probe stations and mechanical failure analysis, respectively. The measured thermal curves were used to develop the predictive thermal model towards low power consumption. Moreover, a statistical analysis was finally introduced to cross-correlate the mechanical failure results and the thermal predictive model, aiming at MH design optimization for gas sensing applications. All the data collected in this investigation are shown.

Dataset:http://doi.org/10.5281/zenodo.4590243

Dataset License: CC BY-NC-ND

1. Summary

The market for gaseous sensors is continuously growing. There is a significant demand for small, low-cost and high-performance gas sensors in several sectors, including chemical process control, precision agriculture, environmental monitoring, and medical applications [1,2,3,4,5]. Among the various new classes of smart gas detectors, one of the most investigated and used are the chemoresistive gas sensors, which show high sensitivity, allowing to detect low concentrations of gases (ppm-ppb range) with high accuracy [6,7,8,9]. These devices are composed of an active sensing material, typically nanostructured metal oxide semiconductors (MOS), and a heated substrate [6,10,11,12]. The substrate plays a key role in the gas sensors devices because it acts as mechanical support of the sensing material and it hosts the heater and electrodes, allowing to heat the sensing material at its best working temperature and to read the electrical resistance of the sensing material, respectively [13,14]. Recently, silicon microheaters (MHs) have replaced the alumina ones, becoming popular in both research prototype development and commercial devices [13,15]. The great advantage of silicon is the use of established microfabrication techniques, which allow very small, reproducible and energy-efficient devices to be produced on an industrial scale [13,15,16]. The low power consumption is obtained by developing very thin suspended membranes of insulating material, obtained by selective etching of silicon in the last microfabrication step [13,15]. These membranes allow to almost eliminate heat dissipation completely by conduction along the substrate [13,15]. However, they are extremely fragile, and their breakage is one of the major causes of failure for silicon-based gas sensors. The use of predictive models to optimize device design is crucial to achieve low power consumption and high mechanical stability [17,18]. This work presents the experimental dataset used in our previous work [17] to develop a specific combined thermal-mechanical model useful to tune MH design for gas sensing application to the required performance. To this aim, three different MH layouts over three different membrane areas were microfabricated and analyzed. The thermal characterization was carried out using a manual probe, while a shear/pull tester was used for the mechanical failure analysis. The dataset was divided into two excel files, both in tidy format (with one row corresponding to one observation), one for the thermal analysis and the other for the mechanical analysis. Both files were supplied with a second sheet reporting the raw thermal and mechanical data in series from which the Yields were obtained. The first sheet of each dataset can be directly imported in R for further analysis. Finally, a statistical method was used to cross-correlate thermal and mechanical analysis. The scheme reported in Figure 1 summarizes the approach used in this work for the development of the predictive thermal-mechanical model for the optimization of MHs, useful for gas sensing applications.

The thermal characterization (Table S1) allows for developing a thermal model for the evaluation of the microheater power consumption, by taking into consideration the heat transfer contributions and a constant empirically calculated (K_EAD, reported in Table S2) [17]. The power consumption values calculated through the developed thermal model were in good agreement with the experimental data [17]. The failure analysis (Table S3) showed that MH layout and membrane size greatly impact on the failure forces [17]. The statistical analysis, used to cross-correlate the thermal and mechanical results, allowed to define the best combination of membrane area and MH layout for the optimization of the microheater performance [17]. To this aim, the desirability factor was introduced for the predictive evaluation of the MH efficiency [17]. A case study is also here reported, in which the maximization of the mechanical stability and the minimization of the MH power consumption were considered for the desirability factor.

This methodology could be further consolidated to readers’ liking by increasing the number of factors defining MH geometry. On this note, our future experimental analysis will expand by varying area layout and area membrane outside the factor domain considered in this work.

2. Data Description

Thermal and mechanical data are organized in separate excel files (Tables S1 and S3), each file is easily accessible to conduct further data analysis. The first sheet of each file is reported in the tidy format: every single row represents an observation; the use of this standard allows the rapid importation in R as well as easy analysis of the data. In both files, the second sheet contains in a progressive manner the collected curve from which the data of the first sheet was extrapolated by analysis. The progressive number is also present in the first sheet (reported as sample number) to easily relate the tidy results with the correspondent curves. In the thermal database, the first sheet reports the Power in mW needed to reach the temperature of 400 °C for each combination of membrane and heater. The second sheet includes the Voltage (V) applied as an independent variable and as dependent variables the measured resistance (Ohm), the measured temperature (°C), and the measured Power (mW). Each sample has then five columns. In the first sheet of the mechanical database, the Force in N at which the rupture of the device occurred is reported for each combination of membrane and heater. The second sheet presents all the Force-Deflection curves in a progressive fashion within their respective sample number. These sheets contain two columns per sample, the first being the deflection (µm) and the second being the measured force (N). The first data points are registered after the indenter, approaching the MH membranes, meets the membrane and ending in the brittle failure after the maximum force is registered. Moreover, the values of the corrective thermal constant KEAD are organized in a separate file (Table S2) and are easily accessible to conduct further data analysis. The first sheet of the file is in tidy format, every single row represents an observation, the use of this standard allows the rapid importation in R as well as easy analysis of the data. The progressive number is also present in the first sheet (reported as sample number) to easily relate the tidy results with the correspondent curves. The experimental values of the corrective thermal constant KEAD are reported as a function of temperature (°C) and of the percentage area occupied by the heater (A%) for sample L1 in sheets 2 and 3, respectively. The calculation and the physical meaning of KEAD is discussed and explained by Gaiardo et al. in [17]. In addition, the calculated values of the corrective thermal constant KEAD are reported as a function of temperature (°C) and of the percentage area occupied by the heater (A%) for samples L2 and L3 in sheets 4 and 5.

3. Methods

3.1. Device Fabrication

The MHs were developed in the clean rooms of the Bruno Kessler Foundation [19]. A Centrotherm E 1200 HT 260-4 4 diffusion and LPCVD furnace was used for the thermal growth of SiO₂ and deposition of Si₃N₄ and SiO₂ layers. An ulvac EBX-16C with Ferrotec EV S-6 e-gun was used for the deposition of Ti and Pt layers. The thermal annealing of the silicon wafers was carried out using an Expertech CTR 200. The SiO₂ intermetal passivation layer was deposited through an STS-MPS PECVD. SiO₂ was etched over pad areas of the heaters by means of an OEM Tegal 903 ACS reactive ion etcher. Layouts of the deposited layers were defined by photolithography technique, using an SVG 8600 Photoresist Coat Track and a Karl Suess automatic mask aligner. Overall, 9 different combinations of layout and membrane were tested (listed in

Table 1. Tested combinations of layout and membrane within their area.

Group	Factor A-Area Layout (mm2)	Factor B-Area Membrane (mm2)
1	0.15	0.81
2	0.15	1.74
3	0.15	3.04
4	0.56	0.81
5	0.56	1.74
6	0.56	3.04
7	0.77	0.81
8	0.77	1.74
9	0.77	3.04

).

3.2. Thermal Analysis

An automatic prober Accretech UF200R, equipped with ATT LOW TEMP System L200T, was used to perform temperature coefficient of resistance (TCR) measurements on all the produced devices. This system provides a nominal Temperature Stability of ±0.1 °C, an Accuracy of ±0.5 °C and a uniformity (along the chuck surface) less than 0.5%. To extract TCR parameters on heather and bulk resistors, resistance measurements were performed at different temperatures: 20 °C, 60 °C, 100 °C, and 140 °C. Electrical measurements were performed using Agilent/Keysight equipment including a low leakage Switching Matrix Mainframe B2201A, with four modules B2211A, and a SMU Mainframe E5270B, with four medium-power high-resolution SMU E5281B, two high-power high-resolution SMU 5280B, and a Ground Unit. Both SMU types have a nominal resolution of 10 fA and 20 uV. Manual measurements were performed using a Karl Suss Manual probing station PM8, equipped with an Agilent 4156C Precision Semiconductor Parameter Analyzer with a nominal resolution of 1 fA and 2 uV,

Table 2. Evaluation of the reported measurements.

Reported Measurement	Derivate From
Temperature	Resistance measurement
Power consumption	Current-voltage measurement
Failure force	Peak of force-deflection curve
Work of fracture	Integral of force-deflection curve

and

Table 3. Measurement sensitivity of thermal and mechanical curves.

	Measurement Sensitivity
Manual IV probe	Current: 1 fA Voltage: 2 uV
Indentation test	Deflection: 30 nm Force: 0.1 mN

. Each measurement was repeated three times. Resistive temperature detection method can display only spatial averages of the temperature distribution while being very reliable and convenient to implement. Furthermore, several resistive elements can be integrated into the device for temperature control in different locations [20]. Conversely, other methods, such as optical thermal imaging methods [21,22] that measure thermal radiation spectrum, can precisely capture the thermal profile of the specific MH geometry [23]. However, instrumentation for precise thermal measurements requiring micrometric pixel resolutions and high confocality is expensive and a proper calibration for emissivity could be complex [24].

3.3. Mechanical Analysis

Tests are intended to evaluate the resistance of different membrane sizes and layouts to screen-printing pressure. The testing apparatus for destructive indentation tests was a shear/pull tester (Condor Sigma, XYZTEC, Panningen, The Netherlands) equipped with a stainless steel flat conical tip with a diameter of 290 µm and an axial load sensor of 20 N and high displacement sensitivity. The test procedure consisted of a manual approach along the z-axis ending within 50 µm of distance from the membrane’s surface. Successively, the tip was optically aligned in the center of the membrane using heater geometry as xy reference. The mechanical characterization was in a quasi-static regime with a vertical speed of 10 µm/s, and data were recorded at approximately 15 Hz for better data visualization. Test distance was set in the range of 150 ÷ 400 µm to reach fracture in all samples depending on the membrane’s deformability, Table 2 and Table 3. Stainless steel tip was cleaned from debris in between tests. Membrane testing was performed on silicon pieces containing arrays of several membranes. Each silicon piece was secured in place on a thick aluminum plate by means of a thin bi-adhesive Kapton tape. The sample size for each membrane set was the following: 20 MHs of 1A, 11 MHs of 1B, 30 MHs of 1C, 44 MHs of 2A, 11 MHs of 2B, 19 MHs of 2C, 8 MHs of 3A, 16 MHs of 3B, and 30 MHs of 3C. The indentation tests are used to register force-deflection curves until fracture of the clamped square membranes. Each curve was analyzed to extract the MH failure force, defined as the maximum of the force deflection curve, and the work of fracture defined as the integral of the curve [17].

3.4. Statistical Methods

A statistical method was used to cross-correlate the results obtained from the thermal predictive model and the mechanical failure analysis, with the aim to define a comprehensive model to tune future sensor layouts to the required performance. The entire statistical analysis has been done by the use of the programming language R [2] following the statistical strategy described in previous works [3,4,5,6]. An initial comparison by verifying the presence of significant difference among the different groups by using the analysis of variance (ANOVA) followed by a Turkey multi-comparison test [7]. The levels of significance were assigned as follows p ≤ 0.1 [19], p ≤ 0.05 (*), p ≤ 0.01 (**), p ≤ 0.001 (***). A Response Surface Methodology (RMS) has been adopted to model the empirical equations relating the considered factors to the yields. In this case, we considered two continuous factors, the area covered by the heater (Factor X) and the membrane surface (Factor Y) and two yields, the produced power (mW), and the mechanical strength (N). Table 1 shows the considered group. The complete model is reported in Equation (1). An ANOVA test followed by a Turkey multi-comparison was conducted to verify the significance of each term of the reported equation. Only the terms with a significant effect (p ≤ 0.01) were included in the model. The function F has been chosen to both normalize the model residues and to make them pattern less. The model was considered significant with a p-value ≤ 0.05. To the goodness of fit of the model, the coefficient of determination (r²) was calculated. Models with a perfect fitting have an r² = 1. A desirability approach has been used to perform the numerical optimization, this method is quite general and has been applied to solve different typology of problems, from the inkjet printing optimization [25,26,27] to biopolymer processing [28,29,30]. This method is based on the desirability functions which value can vary in the [1, 0] range, where 1 represents the optimum solution. A single desirability function was assigned to each of the yield (Y_i). For the maximization of Y_i, the function is reported in Equation (2), for the minimization in Equation (3). In the equations Y_i is the specific yield, d_i the corresponding desirability function, Ui and Li the maximum and the minimum value of the yield, respectively. The overall desirability (D) is the geometric mean of all these functions reported in Equation (4) with k equal to the total number of yields (in our case 2). Then D is plotted against the process factors to find its minimum value and thus the best solution.

$$\mathrm{F}\left(\mathrm{Yield}\right)={ \mathrm{c}}_0 +{ \mathrm{c}}_1\mathrm{X}+{ \mathrm{c}}_2\mathrm{Y} +{\mathrm{c}}_3\mathrm{XY}+{ \mathrm{c}}_4{X}^2 +{ \mathrm{c}}_5{Y}^2 +{ \mathrm{c}}_6{X}^2\mathrm{Y} +{ \mathrm{c}}_7\mathrm{X}{Y}^2 +{ \mathrm{c}}_8{X}^2{Y}^2$$

(1)

$$d_i=\{\begin{array}{cc}\hfill 1& if Y_i\ge U_i\hfill \\ \hfill \frac{Y_i-L_i}{U_i-L_i}& if L_i\ge Y_i\ge U_i\hfill \\ \hfill 0& if Y_i\le L_i\hfill \end{array}$$

(2)

$$d_i=\{\begin{array}{cc}\hfill 0& if Y_i\ge U_i\hfill \\ \hfill \frac{Y_i-L_i}{U_i-L_i}& if L_i\ge Y_i\ge U_i\hfill \\ \hfill 1& if Y_i\le L_i\hfill \end{array}$$

(3)

$$D={\left(d_1{d}_2{d}_3\dots d_k\right)}^{\frac{1}{k}}$$

(4)

4. Data Analysis

4.1. Preliminary Analysis

From a preliminary screening on the results, it is clear that in case of the power (Figure 2A) all the measurements performed on the single group were close to each other, thus between almost all the groups the differences resulted statistically significant. Instead, in case of the mechanical response (Figure 2B) due to the presence of a higher variability on the collected, only some of the differences were significant. Interestingly, a trend could be recognized in both datasets, a decrease in each triplet passing from the first to the last group (1 to 3, 4 to 6, 7 to 8), that could be referred to the effect of the membrane area, and a general increase between the different triplets (1 to 4 to 7, 2 to 5 to 8, 3 to 6 to 9) that could be referred to the effect of the layout area. All these differences were statistically significant. The levels of significance were assigned as follows: p ≤ 0.1 [19], p ≤ 0.05 (*), p ≤ 0.01 (**), p ≤ 0.001 (***).

4.2. Thermal Data Analysis

The ANOVA analysis (

Table 4. ANOVA table of the Thermal Model factors. The levels of significance were assigned as follows: p ≤ 0.1, p ≤ 0.05 (*), p ≤ 0.001 (***).

Source	Sum of Squares	df	Mean Square	F-Value	p-Value	Significance
Model	16,514.52	8	2064.31	1741.77	<0.0001	***
A-Layout Surf.	3947.71	1	3947.71	3330.88	<0.0001	***
B-Membrane Surf.	633.00	1	633.00	534.09	<0.0001	***
AB	225.33	1	225.33	190.12	<0.0001	***
A2	92.46	1	92.46	78.01	<0.0001	***
B2	0.4885	1	0.4885	0.4122	0.5290
A2B	28.66	1	28.66	24.18	0.0001	***
AB2	0.0172	1	0.0172	0.0146	0.9053
A2B2	5.89	1	5.89	4.97	0.0387	*
Pure Error	21.33	18	1.19
Cor Total	16,535.85	26

) indicates that both the considered factors were significant. In addition, the terms AB, A², A²B, A²B² resulted significantly. In order to maintain the model hierarchy, the term B² and AB² has been added. The significance of the two mixed terms indicates the interaction between the effects of A and B on the power outcome (this could be further analyzed by a physical interpretation). The model is shown in Figure 3A. The power increases diagonally with the increasing of Al and the decreasing of L_m. The predicted versus actual plot is shown in Figure 3B, the model well fits the data (the calculated r² = 0.998). Coefficients of the model and its 95% confidence intervals are reported in

Table 5. Estimation of the coefficients of the model and of the 95% high and low confidence interval. The coefficients were estimated in the coded forms in which each factor is normalized in the [−1,1] range.

Factor	Coefficient Estimate	95% CI Low	95% CI High
Intercept	77.49	75.9	79.08
A-Layout Surf.	26.6	25.63	27.57
B-Membrane Surf.	−11.90	−12.98	−10.82
AB	−4.33	−4.99	−3.67
A2	−8.07	−9.99	−6.15
B2	0.5921	−1.35	2.53
A2B	3.06	1.76	4.37
AB2	0.0679	−1.11	1.25
A2B2	2.49	0.1438	4.83

.

4.3. Mechanical Data Analysis

The ANOVA analysis (

Table 6. ANOVA table of the Mechanical model factors. The levels of significance were assigned as follows: p ≤ 0.1, p ≤ 0.05 (*), p ≤ 0.01 (**), p ≤ 0.001 (***).

Source	Sum of Squares	df	Mean Square	F-Value	p-Value	Significance
Model	7.03	8	0.8785	83.09	<0.0001	***
A-Layout Surf.	0.1156	1	0.1156	10.94	0.0039	**
B-Membrane Surf.	3.07	1	3.07	290.17	<0.0001	***
AB	0.0874	1	0.0874	8.27	0.0101	*
A2	0.0992	1	0.0992	9.38	0.0067	**
B2	0.3511	1	0.3511	33.21	<0.0001	***
A2B	0.3574	1	0.3574	33.80	<0.0001	***
AB2	6.436 × 10−6	1	6.436 × 10−6	0.0006	0.9806
A2B2	0.2120	1	0.2120	20.05	0.0003	***
Pure Error	0.1903	18	0.0106
Cor Total	7.22	26

) indicates that only the factor B was significant. In addition, the terms AB², A²B² were significant. In order to maintain the model hierarchy, all the remaining terms have been added. The significance of the two mixed terms indicates the interaction between the effects of A and B on the mechanical outcome (as with the previous case, this could be interpreted with a physical explanation). The model is shown in Figure 4A. The mechanical response increases with the increasing of the layout surface and the decreasing of the membrane surface. The predicted versus actual plot is shown in Figure 4B: as shown, the model is more complex than the previous. Furthermore, in this case, the model well fits the data (the calculated r² = 0.924). This high value of r² is partially due to the higher value of the Yields, in fact by magnifying the predicted versus actual plot on the lower range some scattering of the data points around the diagonal is present. However, in this case a trend is also well recognizable. Coefficients of the model and its 95% confidence intervals are reported in

Table 7. Estimation of the coefficients of the model and of the 95% high and low confidence interval. The coefficients were estimated in the coded forms in which each factor is normalized in the [−1,1] range.

Factor	Coefficient Estimate	95% CI Low	95% CI High
Intercept	−2.88	−3.03	−2.73
A-Layout Surf.	0.144	0.0525	0.2354
B-Membrane Surf.	−0.8283	−0.9305	−0.7261
AB	0.0854	0.023	0.1477
A2	0.2644	0.0831	0.4458
B2	0.502	0.319	0.685
A2B	0.3422	0.2185	0.4659
AB2	−0.0013	−0.113	0.1104
A2B2	−0.4722	−0.6938	−0.2507

.

4.4. Optimization—Minimization P, Maximization s

In the specific case of maximization of the mechanical response and minimization of the power, the problem results as non-trivial and particularly complex since these two constraints are in countertendency. The best point if we treat both the yields in the same manner (they are considered both with the same importance) is the one depicted in Figure 5 and reported in

Table 8. Results of the optimization. Among the different optimizations with the same results the first was selected as solution.

Number	Layout Surf.	Membrane Surf.	Power	Strength	Desirability
1	0.148	0.810	50.333	0.148	0.746
2	0.148	0.857	49.892	0.146	0.745
3	0.148	1.000	48.630	0.140	0.741
4	0.148	1.077	47.985	0.137	0.738

with a desirability of 0.746. One of the tested points gave the highest desirability.

5. User Note

Datasets are reported in excel format easily convertible in csv for a further importation in third software to perform statistical analysis. The data reported in the first sheet of both the Thermal and the Mechanical datasets are in tidy format and importable in R to perform statistical analysis. In the tidy format each row corresponds to a single observation.