Gas discrimination by simultaneous sound velocity and attenuation measurements using uncoated capacitive micromachined ultrasonic transducers

Chemically functionalized or coated sensors are by far the most employed solution in gas sensing. However, their poor long term stability represents a concern in applications dealing with hazardous gases. Uncoated sensors are durable but their selectivity is poor or non-existent. In this study, multi-parametric discrimination is used as an alternative to selectivity for uncoated capacitive micromachined ultrasonic transducers (CMUTs). This paper shows how measuring simultaneously the attenuation coefficient and the time of flight under different nitrogen mixtures allows to identify hydrogen, carbon dioxide and methane from each other and determine their concentration along with identification of temperature and humidity drifts. Theoretical comparison and specific signal processing to deal with the issue of multiple reflections are also presented. Some potential applications are monitoring of refueling stations, vehicles and nuclear waste storage facilities.

Nowadays, most gas sensors present a chemical component such as a functionalized coating specifically engineered to target a particular analyte 1,2 . However, it is well known that such chemical coating presents long term stability issues 3 . This results in the need for recalibrating the sensor, typically every couple of months. For applications dealing with hazardous gases such as hydrogen, H 2 , or methane, CH 4 , this becomes dangerous for the operators. Alternatively, it is possible to compensate this drift digitally, for instance, with artificial intelligence 4 or sophisticated models 5 . Nevertheless, this naturally increases both the manufacturing and the development cost. For this reason and despite their typically lower selectivity and worse limits of detection (LOD) 2 , the development of uncoated gas sensors has increased in popularity over the past decade. Indeed, some parameters such as the resonant frequency of a cantilever 6 , the gas sound velocity 7 and the acoustic attenuation coefficient 8,9 depend on the gas physical properties. This allowed the development of gas density sensors 10 and in some cases these principles can be exploited to measure a binary gas mixture concentration 11,12 . Some of these sensors consist of capacitive micromachined ultrasonic transducers (CMUTs) 13 . Typical LODs associated to uncoated sensors go from 3% down to 100 ppm and they are usually considered as non selective. However, in some specific cases, uncoated sensors can show performances close to coated sensors such as a LOD of a few ppm 14 and even selectivities higher than 10 15 . In order to overcome the problem of selectivity, an alternative consists in discrimination by measuring multiple physical properties of a gas mixture such as mass density and viscosity 16 . In this paper, we propose to measure simultaneously two other gas properties; the acoustic attenuation coefficient, α , for frequencies ranging from 1 to 4.5 MHz, and the time of flight, τ , which is the time taken by an acoustic ultrasonic wave to travel from an emitter to a receiver passing through the gas to be characterized. In such a setup, α can be defined as: where d is the distance between the emitter and the receiver and P e and P r are the amplitudes of the acoustic wave at the emission and the reception ends, respectively. This work is the continuation of previous work showing the Temperature and humidity characterization. This section seeks to estimate the sensor cross-sensitivity to the temperature, T, between 20 and 50 • C of both the measurement of τ , S T τ , and of �α , S T α as well as the cross-sensitivities to the relative humidity, RH, noted S RH τ and S RH α . In order to achieve this, both τ and �α where measured at different temperatures under pure dry N 2 in that temperature range. The resulting characteristics are shown in Fig. 1a,b. Similarly, both quantities, τ and �α , were measured at 20 • C with RH going from 0% to about 50% (Fig. 1c,d).
In the case of the measurement of τ , the theoretical value can be estimated for an ideal gas by exploiting Newton-Laplace formula for the sound velocity, v, of an ideal gas 7 : where P is the atmospheric pressure, ρ the gas mass density, γ the heat capacities ratio (considered constant in this range of temperature) and M the molar mass of the gas. Differentiating Eq. (2) results in: where v and T are small variations of v and T, respectively. Since v = d/τ, where �τ represents small variations of τ . Finally, S T τ is given by: For pure nitrogen at 20 • C , τ = 21.35 µs (distance d ≈ 7.5 mm ), the theoretical value is S T τ = −36 ns/K against the experimental value of −31 ns/K . The difference might come from the relatively large temperature variation range or the fact that the sensor and the thermocouple used to measure the temperature are a few centimeters apart inside the oven. However, as it will be shown later, this difference is sufficiently small for the purpose of this article.
In the case of the humidity, one can express the variations of τ with respect to the molar fraction of vapor water x w in N 2 17 , where ρ w = 0.739 kg/m 320 , c w p = 1.97 kJ/K/kg 20 , c w v = 1.50 kJ/K/kg 20 are the mass density, the isobaric heat capacity and the isochoric heat capacity of water vapor, respectively and ρ N2 = 1.16 kg/m 320 , c N2 p = 1.04 kJ/K/ kg 20 , c N2 v = 0.743 kJ/K/kg 20 the same properties for N 2 . Since RH is the mass of water in the gas with respect to the maximum amount of water that can be diluted in N 2 , it is given by 20 : www.nature.com/scientificreports/ where P s = 2310 Pa 20 is the saturated water vapor pressure at 20 • C and M w = 18.02 g/mol 20 the molar mass of water. Therefore the sensitivity of τ to RH is: Once again the measured value of − 0.70 ns/%RH is consistent with the theoretical one of − 0.69 ns/%RH. In the case of �α , due to the complexity of its analytical expression and its dependence on an important number of gas properties, the temperature and humidity influence on the attenuation in this study is purely experimental. From Fig. 1, it can be concluded that, in the range between 2 and 4 MHz, the sensor's temperature cross-sensitivity increases with the frequency whereas the one of humidity follows the opposite behaviour. Their values will be exploited in the following section.
Gas characterization. The temperature and humidity influence on the sensor being known, this section aims to characterize the sensor response to different gas mixtures. For that, the sensor was exposed to sequences where pure N 2 is alternated with binary mixtures of N 2 with either H 2 , CO 2 or CH 4 at different molar concentrations, x of the targeted gas. An example of the measurements at 4 MHz is shown in Fig. 2 for both, (a) the shift in time of flight with respect to N 2 , �τ , and (b) �α . From these measurements one can extract the sensitivities through calibration curves, both of �τ (c) and �α (d) to the targeted gas ( H 2 , CO 2 or CH 4 ), S τ and S α , respectively. Although a complete study on the sensor drift (particularly degradation of the mechanical properties of www.nature.com/scientificreports/ the membrane) is out of the scope of this work, it should be noted that, at constant temperature with dry gases, no detectable drift was observed for several hours of experiment. This emphasizes one of the advantages of uncoated sensors with respect to functionalized sensors. Since S α depends on the excitation frequency, it was measured for the whole frequency range of 1 MHz up to 4.5 MHz (e) which is limited by the charge amplifier at the lower end and by the strong attenuation above 4.5 MHz. Such frequency range is convenient to prevent typical vibration sources to interfere with the measurements 21 although it might not be adapted for applications where other ultrasonic sources on the same frequency range as the measurements are present. Similarly to the derivation of the expression of the humidity cross-sensitivity, the expression of S τ can be obtained 17 : where c G p and c G v are the isobaric and isochoric heat capacities per unit mass of the targeted pure gas, respectively, and c N2 p and c N2 v the ones of pure N 2 . To derive an exact expression of S α would be unnecessarily complex. Therefore, in this study, a simple model is used from the expression of the attenuation of a pure gas α p 22 : where η is the viscosity of the gas, η v its volume viscosity, K its thermal conductivity, c p its isobaric heat capacity per unit mass and α vr the vibrational relaxation attenuation contribution modeled as 23 : where A is a constant depending on the strength of the relaxation process and f vr the frequency corresponding to the maximum of the normalized attenuation coefficient αv/f . When dealing with a binary gas mixture, an additional contribution of the diffusion process α δ is needed, especially for molecules with very different molar masses M such as H 2 and N 2 24 : where a is the thermal diffusion factor of the mixture, δ its diffusion coefficient and M is the difference in molar masses between both constituents of the mixture. The shift in attenuation with respect to N 2 is then approximated to: where α N2 and α G are the attenuation coefficients in the pure state of N 2 and the targeted gas respectively, as given by Eq. (11). Finally, S α is simply �α(x = 1%).
All the necessary constants to compute both S τ and S α for the three mixtures studied in this article are given in Table 1  www.nature.com/scientificreports/ there seems to be some additional losses than the ones predicted by the linear model from Eq. (14). This is due to the beginning of the relaxation attenuation peak of H 2 , which is already an important and non-linear component of the attenuation 8 . Nevertheless, in a practical application, this can be solved with a calibration mapping �α to the molar fraction of H 2 (Fig. 2d).
Discrimination. So far, the measurements of �τ and �α have been displayed separately. However they are always measured simultaneously. Therefore, one can exploit this in order to identify the different mixtures. Indeed, as it is shown in Fig. 3a, in theory (Eqs. 10 and 14), when plot �α against �τ , the three different types of mixtures are placed in different parts of the plane. Additionally, since in absolute value, S α increases with the frequency, it is easier to distinguish between them at higher frequencies (4 MHz rather than at 2 MHz). Figure 3b shows the corresponding measurements resulted from the previous characterization as well as the influence of the temperature and humidity, taking into account the cross-sensitivities measured in "Temperature and humidity characterization" section. In the case of the humidity the values above 50% RH were extrapolated to estimate the influence up to 100% RH assuming that the cross-sensitivities of both τ and �α remain constant. Overall the humidity influence remains weak compared to the one of the temperature. In order to illustrate what would happen with a contaminant which is not a pure gas but a mixture of two gases with a fixed constitution (1:1 ratio in this case), measurements of a mixture of H 2 and CO 2 at 4 MHz at different concentrations (1% of the mixture contains 0.5% of each gas) are shown (Fig. 3b). As one would expect, the points fall in between the curves at 4 MHz of each of the constituents. In Fig. 3c, preliminary results at 4 MHz of how the measurements using N 2 as reference gas are modified by using a different reference such as a mixture of N 2 and CO 2 are shown. In fact, if the reference gas is changed, the sensitivities are modified and consequently the slopes of the curves in the plane ( �α , �τ ) change as well. However, if the reference gas is known, this can be taken into account during the cali- www.nature.com/scientificreports/ bration step (Fig. 2c, d). Finally, additional measurements detecting 4% of each of the three mixtures ( N 2 with either H 2 , CO 2 or CH 4 ) at 30 °C and at 1, 2 and 3 MHz were performed and displayed with the corresponding measurements at 20 °C in Fig. 3d. From them it can be observed that the main influence of the temperature is an horizontal displacement (shift in τ ) with a relatively small displacement along the axis of �α.

Discussion
The CMUT ultrasonic gas sensor presented in this paper allows to measure simultaneously both the time of flight and the attenuation coefficient in a binary gas mixture. This allows to discriminate, at a fixed temperature and fixed humidity, between the proposed binary mixtures of N 2 with either H 2 , CO 2 or CH 4 . Remarkably, a change in the sensor output due to the temperature or humidity can be distinguished by a small change in the attenuation coefficient and a large shift in the time of flight especially for the temperature. This means that, to some extent the sensor presents also temperature or humidity sensor capabilities that can help for numerical compensation. This can be exploited even further by the fact that a slight increase in T from 20 to 30 °C or change in relative humidity results mainly in a substantial decrease in τ with a small shift in �α . This means that if a slight loss in accuracy (due to the minor distortion in attenuation caused by either the temperature or the humidity) can be afforded, whenever a change in T or RH is detected (delimited by the dark and light red zones in Fig. 3b), the sensor can perform a self recalibration compensating such effect. Additionally, after the mixture is identified, the concentration can be determined with calibration curves. Concretely, the simultaneous measurement of time of flight and attenuation coefficient with CMUTs results in a robust three binary gas mixtures sensor with built in temperature or humidity sensor capabilities and long term stability. Some potential applications are monitoring of refueling stations, vehicles and nuclear waste storage facilities. Preliminary results in the case of more complex mixture have been presented as well which showed that the calibration step has to be adapted to the application.

Methods
CMUT details. An optical microscope picture of the CMUT array used in this study is shown in Fig. 4 along with 3D schematics from different perspectives of single CMUT. The microfabrication protocol used is similar to the one presented in 25 . Their geometrical characteristics along with polarization and excitation details are given in Table 2.
Setup. The setup to measure both the shift in attenuation, �α , and the sound velocity, v, for different gases at different concentrations is shown in Fig. 5. The gases from industrial grade bottles (1) are mixed by a system of flowmeters (2) controlled by a computer (3) at the desired concentrations. The mixture passes through an oven (4) where the desired temperature is set and verified by a thermocouple (5). The gas enters the measurement cell (6) where a continuous ultrasonic sine wave is sent from an emitter CMUT array to a receiver CMUT array after travelling a distance d. Thus, the wave arrives attenuated exponentially with d. Additionally its phase is shifted with respect to the emission by an amount proportional to the time of flight, τ = d/v , required to go through the measuring cell. The receiver is connected to a charge amplifier (7) in order to improve the signal to noise ratio www.nature.com/scientificreports/ before being fed along with the emission signal to a gain/phase network analyzer (8). Finally, for the humidity measurements, a vapor generator PUL110 was used to humidify N 2 . It should be noted that although the measuring setup presented is rather bulky, more compact components such as the network analyzer can be built 34 .
No reflections model. Far from the resonance of the CMUTs, on the side of lower frequencies, one can assume that their electromechanical transfer function is a constant. Thus, using the complex formalism, the pressure wave p can be considered of the form: where P 00 is the amplitude of the wave at t = 0 on the emission end, ω the excitation angular frequency, j = √ −1 , t the time, z the coordinate orthogonal to the CMUT array (Fig. 5), α the attenuation coefficient of the gas and k the wave number which, in the case of a plane wave, is linked to the gas sound velocity v as follows: (15) p(t, z) = P 00 e −αz e j(ωt−kz) ,  www.nature.com/scientificreports/ Since the wave is continuous, at any given time t 0 ≥ τ , the value of p at the receiving end ( z = d ) is: Thus, its phase is shifted by a value of −ωτ ( τ = d/v ) with respect to the value of p at the emission end at the same time t 0 : The transfer function due only to the setup, H 0 , independent of the gas can be defined as: where φ 0 is a constant and accounts for the shifts due to the setup and |H 0 | is the modulus of H 0 , which depends only on the angular frequency. Therefore the total phase shift φ is: As a consequence τ is simply the slope of the curve φ(ω) . Additionally, the change in amplitude is simply e −αd . Hence the total gain |H| is given by: In order to overcome the need of knowing |H 0 | , a calibration step was performed to measure |H| under N 2 , |H| N2 . Indeed, the shift in attenuation with respect to N 2 , �α , is given by: Dealing with multiple reflections. In order to measure a high signal to noise ratio, especially for higher frequencies, it is possible to increase the excitation amplitude. However, this will cause reflections that will interfere with the main emission, particularly for lower frequencies, where the attenuation is smaller. Therefore, they will modify both φ and |H|. This section seeks to understand the extent of this effect and, particularly, how to retrieve information about both α and τ . Assuming normal incidence and a reflection coefficient of r identical for both CMUT arrays, at a given position z and time t, the acoustic wave in the cell is the superposition of a wave travelling forward p f (towards increasing z) and one travelling backwards p b (towards decreasing z). The former is the sum of all the reflections that happen an even number of times (0 included) between the emission time and the measurement time i.e.: where n is the echo order (initial 0). Similarly, p b is the sum of the waves that have been reflected an odd number of times.
Notice that the term of order 0 is not to be accounted since the emission happens only towards the increasing z. Both p f and p b are geometric sums of same quotient q given by: Since |q| < 1 , p f and p b are simply given by: (17) p(t 0 , d) = P 00 e −αd e jωt 0 e −jωd/v . (18) p(t 0 , 0) = P 00 e jωt 0 .   www.nature.com/scientificreports/ The total acoustic signal is finally One can verify that, the case n = 0 corresponds exactly to the model presented in "No reflections model" section. In permanent regime, one can consider the number of echoes to be infinite, resulting in the final acoustic signal p ∞ : The total transfer function H is then given by: which equals: It should be noticed once again that without any reflections ( r = 0 ) H(jω) is equal to the transfer function with a single peak: In order to extract information about α some additional calculations are needed. Indeed, by using the fact that for |y| < 1: and some basic series manipulations Eq. (33) can be rewritten as: where It is important to point out that all the P n are independent of τ . From this form it can be seen that H is the sum of complex exponentials modulated by the functions P n . Applying the inverse Fourier transform to H (noted Ĥ ) defined as: will result in: where P n is the inverse Fourier transform of P n . Therefore, Ĥ is the sum of peaks, where the envelope is the shifted inverse Fourier transform of P n . In particular, the amplitude of the second peak is given by: Therefore, by isolating |P 1 | it is possible to apply the same formula from the single reflection case (Eq. 22) to retrieve �α . Figure 6a shows the measured result of Eq. (38) under pure N 2 as a function of the delay with respect to the emission, t, which is consistent with Eq. (39). It should be noted that, apart from the main peak P 1 (at t = τ ) and the reflected echo P 2 (at t = 3τ ) there is an additional peak at t = 0 . This is due to an instantaneous electrical coupling between the emitter and the receiver resulting from the high voltages applied 17 . However, since the objective is to isolate P 1 , the coupling won't influence the final result. Figure 6b shows the wavelet transform of the measurement from Fig. 6a and shows that the main peak contains more spectral information, since the highest frequencies are no longer present in the reflected echo. This supports the fact that in order to have information about the attenuation at the highest frequencies it is useful, and sometimes necessary, to accept the presence of reflections of the lower frequency components constituting the acoustic wave. However, this perturbs strongly H. This can be appreciated in Fig. 6c,d, which show the comparison between the measured transfer function (29) u 0 = 1 r P 00 e j(ωt+kz)+αz .
(31) p ∞ (z) = p f ,0 (z) + qu 0 (z) 1 − q .  www.nature.com/scientificreports/ (both with multiple reflections and electrical coupling) and the result of the data processing after retrieving only the contribution of P 1 by filtering out components of Ĥ , which are outside of the range t ∈ [15,30] µs . This numerical criterion should be modified according to the application (v and d). In our case, τ is not expected to vary more than 1 µs . It should be noted that, after filtering, |H| is smooth as one would expect from the case without reflections. Moreover, the unwrapped (2π-modulo removed) phase after filtering is consistent with the model of a single peak as shown in "Results" section its slope corresponds to −2πτ.