1. Introduction

Gas sensors are in high demand in different fields such as environmental monitoring, industrial safety and quality control [1]. Available sensors can be categorized based on the way they interact with the gas: optical, chemical, acoustical, conductive and calorimetric sensors to name a few [1–4]. An optical sensor, such as the one presented here, has the benefit of high selectivity due to the characteristic, specific gas absorption lines, that function as a fingerprint for the gas. Moreover the interaction with the gas is reversible (as opposed to e.g. some chemical sensors), allowing for a long sensor lifespan. Optical spectroscopy can be carried out for a wide range of electromagnetic wavelengths, ranging from the visible to the mid-infrared and the terahertz (THz) spectral regions. The terahertz frequency band is of particular interest due to the presence of many rotational energy levels in a wide range of polar gases [5]. Mid-IR is also often used, typically to excite vibrational energy levels [6]. An important drawback to the THz band is that typical THz photon generation and detection processes require instruments that need to be cooled down to a couple of Kelvins [7]. This is due to the small photon energy compared to the thermal energy $k_B T$. The need for cooled detectors can be mitigated by transducing the optical signal, for example using photoacoustics. Here we designed a photoacoustic sensor that is sensitive enough to operate at low power. The photoacoustic technique has been studied most successfully in the so-called quartz-enhanced-photoacoustic spectroscopy configuration (QEPAS [8]). Here an electromagnetic beam is focused between two prongs of a small quartz tuning fork. Optical absorption is transduced into sound which excites a vibration in the fork. This vibration in turn is read out as a voltage through the quartz’s piezoelectric properties. This was initially done at mid-infrared frequencies (achieving ppt detection limits [8]). Terahertz waves have also been employed [8], using quantum cascade laser sources that have to be cooled down to cryogenic temperatures and the setup is very bulky. The device proposed in this paper is an on-chip photoacoustic gas sensor, that can detect concentrations down to 234 ppb of methanol using 1 mW of power and 1 ms of integration time (room temperature, compact terahertz sources with a power of 1.6 mW at 640 GHz have been shown [9]). The device does so by employing resonators at each newly generated form of energy: an optical resonator, an acoustical, and a mechanical one. The geometry of the sensor is shown in Fig. 1 where the locations of the different resonators are pointed out: Firstly externally generated terahertz waves are coupled into the waveguide through an inverted taper. The photonic platform consists of an undoped high-resistivity air-cladded device layer. The optical resonator is made by appropriately positioning holes (radius $R=29$ $\mu$m) etched into the waveguide, forming a 1D photonic crystal cavity. The gas in each of these holes will expand upon absorption of the terahertz waves. When the terahertz source is modulated a pressure wave is generated. The cylindrical holes form individual acoustic cavities for these waves. The acoustic power will be the strongest in the central hole (where most of the optical field is, see Fig. 1(b)), the other holes will hence not be considered. The holes are covered on one side with a thin poly-crystalline silicon membrane. This membrane gets excited by the acoustic pressure and its eigenfrequency is matched to the acoustic resonance frequency, such that its displacement is resonantly enhanced. The flow of energy is also illustrated in Fig. 2. The signal is read out externally through the mechanical resonator’s displacement. By employing three different resonators the signal is greatly enhanced. This high sensitivity allows one to excite the sensor with low-power photomixing sources [10] operating at room temperature. Moreover we avoid the use of large, cooled photodetectors by transducing the signal. Also due to this transduction the sensor suffers less from signal to noise issues at low gas concentrations compared to classical absorption spectroscopy. The latter aims to detect a small drop in the signal intensity, whereas this sensor generates a sound signal on a silent background. Currently we propose using a commercially available laser Doppler vibrometer (LDV) [11]. In future we hope to be able to use an integrated version of the LDV. A proof of concept can be found in [12], though it is far from sensitive enough to achieve the sub-ppm detections claimed here. The chip can be fabricated at a very low cost by leveraging the, already widely used, CMOS and MEMS fabrication technologies, although an external source and read-out are also required. This paper explains the physics of the sensor and provides a model for the signal in function of the concentration of the trace gas (methanol). In Sec. 2 the transduction from light absorption into pressure, into a displacement of the membrane is explained and the three resonators are discussed individually. Sec. 3 concerns the results from finite element simulations. Firstly the optical absorption is modelled in Sec. 3.1 and then the acousto-mechanical repsonse is given in Sec. 3.2 (more details on the simulation can be found in Appendix C). These results are used to provide a single formula that gives the displacement of the membrane in function of the THz source power and methanol concentration in Sec. 4, where a sensitivity limit is also given when using a commericially available laser Doppler vibrometer.

 

Fig. 1. A drawing of the sensor. a) The geometry: a high-resistivity silicon waveguide (cross section $210\times 90\,\mu m^2$) is tethered to the surrounding structure using anchors. At the edge an inverted taper allows for more efficient coupling of externally generated terahertz waves into the waveguide. Holes (radius 29 $\mu$m) are etched into the waveguides and covered with thin (100 nm) poly-Si membranes (red). The holes are etched at specific positions such that a photonic crystal cavity is formed. b) A 2D cut of the optical field profile of the 636 GHz resonant mode of the photonic crystal cavity, generated using COMSOL. c) The resonant acousto-mechanical mode formed in the half-open central hole. The coupling of both resonators defines the eigenfrequencies of the system (Appendix C). The fundamental one is shown here ($f_1=$446.5 kHz). Only the central hole is considered here, as this is where most of the THz-gas interaction and hence heat production (through absorption) takes place. The transduction of the optical absorption into mechanical motion is schematized in Fig. 2.

Download Full Size | PDF

 

Fig. 2. The operation principle of the photoacoustic gas sensor. Terahertz waves are coupled into the chip and their intensity is locally enhanced in the photonic cavity. Upon absorption the energy is turned into heat and hence a pressure increase in the gas. This is modulated at $f_{mod}$ (through modulation of the source) to create pressure waves with the same frequency $f_{mod}$. The holes in the photonic cavity are covered with a movable membrane such that they form a coupled acousto-mechanical system, $f_{mod}$ is tuned to their fundamental frequency. This means that acoustic waves will be trapped in the cylindrical hole which will in turn excite a resonant motion in the membrane. The movement is then read out externally with a laser Doppler vibrometer (LDV). This transduction of terahertz into motion happens in all the holes of the photonic crystal cavity, but only the central one is considered as it has the highest concentration of THz radiation and will hence generate the strongest signal.

Download Full Size | PDF

2. Operation principle: a triply resonant structure

This section is dedicated to the description of the signal generation in the passive sensor, as shown in Figs. 1 and 2. A THz source is coupled into an optical waveguide on-chip (Fig. 1(a)). A one dimensional photonic crystal cavity is defined into the waveguide, through the etching of holes in specific locations, which traps waves of a specific terahertz frequency $f_{THz}=$636 GHz (corresponding to an absorption line of methanol). This trapping will result in an enhanced optical field and hence more THz-matter interaction with the trace gas that naturally fills the holes. The THz field will have its highest intensity in the central hole (Fig. 1(b)). The trace gas will upon absorption of the electromagnetic wave heat up and increase its pressure. This mechanism is exploited by modulating the incoming electromagnetic wave at a frequency $f_{mod}$, which generates a pressure wave with the same frequency. This pressure wave is generated mostly in the central cylindrical hole of height $H=90$ $\mu$m and radius $R=29$ $\mu$m (where the field intensity and hence heat production through absorption is stronger). The hole is covered by a thin poly-crystalline silicon (poly-Si) membrane with a thickness of $100\,nm$ and internal tensile stress of 1.5 MPa (Fig. 1(c)). Such a half-open tube will have acoustic resonant frequencies, the fundamental one $f_{ac}$ having a wavelength which is roughly four times the height: $\lambda _{ac}=4H$ (Sec. 2.2). The thickness and tensile stress of the membrane, covering the cylindrical hole, are tailored such that their fundamental resonant frequency $f_{mech}$ matches $f_{ac}$. The coupling of the mechanical and acoustic resonators leads to two coupled resonances at $f_1=$446.5 kHz and $f_2=$834.2 kHz. The modulation frequency is chosen to excite the fundamental coupled resonance $f_{mod}=f_1$ (Fig. 1(c)). In summary terahertz waves of frequency $f_{THz}$ are modulated at frequency $f_{mod}$. The terahertz waves are locally enhanced in the central hole of the photonic crystal cavity. The trace gas will absorb the terahertz waves to excite an acoustic resonance in the central hole at frequency $f_{mod}$, which in turn excites a resonant mechanical motion in the poly-Si membrane of frequency $f_{mod}$. The movement is read out externally through a laser Doppler vibrometer, which determines the membrane’s displacement. In the rest of this section the mechanisms of the three resonators are evaluated independently. The acoustic and mechanical eigenfrequencies and the losses are discussed. The calculation of the acoustic losses can be found in Appendix A. These acoustic losses are included in the acousto-mechanical simulation in Sec. 3.2. The mechanical losses are discussed in Sec. 2.3 (for details see Appendix B) and serve as a parameter for the simulation in Sec. 3.2. In Sec. 3.1 the amount of terahertz radiation absorbed in the central hole of the photonic cavity is determined. This absorbed power then serves as a thermal source for the acousto-mechanical model (Sec. 3.2), details on this transduction be found in Appendix D. The governing equations for the acousto-mechanics and details on the solver can be found in Appendix C.2. In Sec. 4 the simulation results are combined to give a function for the membrane’s displacement in relation to the concentration of methanol for a given optical source power.

2.1 Terahertz resonator: a one dimensional photonic crystal

The structure consists of an air-clad waveguide made in undoped high-resistivity (>10 k$\Omega$.m) [13] silicon. This is because impurities in silicon introduce free carriers and these can efficiently absorb terahertz frequencies. The index contrast, necessary for the guiding of the mode, is provided by the air cladding. The waveguide is suspended in air by being tethered to a surrounding structure (Fig. 1(a)). The waveguide is 210 $\mu$m wide and 90 $\mu$m tall, resulting in a confinement factor of 0.9 for the fundamental TE-mode at $636$ GHz (meaning that 90% of the optical power is confined to the waveguide). Externally generated terahertz waves are coupled into the waveguide by using an inverted taper with a length 3 mm and tip width 66 $\mu$m, as the mode at this tip is less confined and hence has more overlap with larger modes. The actual coupling efficiency depends on the source profile (e.g. free-space vs. guided). The coupling loss is not considered in this discussion. A photonic crystal effect can be achieved by etching holes periodically into the waveguide. This will modulate the refractive index and introduce Bloch modes. This means that only electromagnetic waves of specific frequencies can propagate along the photonic crystal, while other (forbidden) frequencies will be reflected. In Fig. 3 the first two allowed modes are shown for an illustrative period of 120 $\mu$m and hole radius $R=29$ $\mu$m. It is clear that the fundamental one (left) has most of its field concentrated in the high index material (silicon, dielectric mode), while the other is concentrated mostly in the holes (air mode). This air mode will be of special interest as it has a better overlap with the gas that penetrates these holes. In between these allowed modes we have a set of forbidden modes. The propagating frequencies can be tailored by designing the radius $R$ or period $P$ of the holes. In Fig. 4 the first two allowed modes are shown, as well as the center of the bandgap, in function of the period in between holes with radius $R=29\mu$m. To make a cavity we alter the distances between the holes such that the THz waves are reflected on both ends but propagate in between them. The algorithm used to define the positions of the holes is similar to the one described by Bazin et al.[14], namely a quadratic decrease of the inter-distances $P_i$ from the center outwards. This method results in a Gaussian field profile and less loss [14]

 

Fig. 3. The first two (propagating) Bloch modes for a photonic crystal with period of 120 $\mu$m and a hole radius of 29 $\mu$m. The fundamental one (522 GHz) has most of its field concentrated in the higher index material (dielectric mode), while the next (620 GHz) has most of its field concentrated at the holes (air mode). In between forbidden modes are found, which will be reflected.

Download Full Size | PDF

 

Fig. 4. The first two propagating TE modes and bandgap center for a photonic crystal made up of silicon waveguide with dimensions $210$ $\mu$m $\times$ $90$ $\mu$m and cylindrical holes with radius $R=29$ $\mu$m in function of the period $P$ of the holes. The dielectric mode is shown in blue, the air mode in green. The red dashed line is centrally located between them, those frequencies will be maximally reflected.

Download Full Size | PDF

2.2 Acoustic resonator

Following up on the discussion about the optical resonator, this section details the acoustic resonator. The eigenfrequency $f_{ac}$ is important to determine, as the mechanical eigenfrequency $f_{mech}$ has to be matched to it. The loss mechanisms are included in the acousto-mechanical simulation discussed in Appendix C. The acoustic quality factor $Q_{ac}$ is analytically discussed and referenced to literature in Appendix A. This value is not used in further simulations however, but allows one to evaluate how efficiently the cavity can trap acoustic waves. The acoustic resonator is formed by the half-open cylindrical central hole in the photonic crystal cavity (Fig. 1(c)). Neglecting the effect of the finite width of the resonator ($R\rightarrow 0$), its resonant frequencies are odd multiples of its fundamental one $f_0$

2.3 Mechanical resonator

In this section the resonance frequency of the poly-Si membrane and its loss mechanisms (in vacuum) are discussed. The knowledge of the resonance frequency as a function of the membrane’s thickness and internal stress is of importance as it has to match the resonance frequency of the acoustic wave $f_{ac}$. The discussed loss mechanisms result in a mechanical quality factor $Q_{mech}$. The losses in vacuum ($Q_{vac}$) are calculated separately from the membrane’s movement in air (discussed in Sec. 3.2). Since the damping mechanisms are not influenced by the type of excitation, we can keep computational costs manageable by including these damping mechanisms in further simulations through a redefining of the material’s Young modulus $E$ as

3. Finite element model

In the previous section the individual resonators were evaluated. Here the numerical model is discussed. First an optical simulation is performed to calculate how much THz radiation is absorbed by the trace gas (Sec. 3.1). Next the performance of the acousto-mechnical resonator is discussed (Sec. 3.2). The losses of the mechanical membrane calculated in Sec. 2.3 are included through Eq. (4). An equation relating the displacement of the membrane to the thermal power applied in the cavity is presented (Eq. (16)). In Sec. 4 these results are combined to discuss the performance of the sensor, resulting in an expression for the displacement of the membrane in function of the concentration of trace gas (methanol) $x_{trace}$ and the applied source power $P_{THz}$.

3.1 Optical absorption

In this section the fraction of THz radiation absorbed by the trace gas in the central hole of the photonic crystal cavity is calculated as a function of its concentration. This will operate as a thermal source in Sec. 3.2. The absorption of the trace gas of interest (methanol) is modeled by giving the material inside the central cylindrical hole the appropriate imaginary refractive index $\kappa _{gas}$. This approximation holds when operating at concentrations where the absorption line is broader than the photonic resonator’s linewidth, such that the frequency dependence of $\kappa _{gas}$ can be neglected. The absorption in the other holes of the photonic cavity is not taken into account as the effect is negligible. The quality factor $Q_{THz}$ is calculated numerically as a function of $\kappa _{gas}$. The quality factor due to the gas absorption $Q_{gas}$ can be separated from all other losses (scattering and absorption in silicon) $Q_0$ through

 

Fig. 5. The contribution of trace gas absorption $Q_{gas}$ to the total optical quality factor $Q_{THz}$ as extracted from simulations using Eq. (8) and the fit according to Eq. (9).

Download Full Size | PDF

3.2 Acousto-mechanical response

A homogeneous thermal source is applied to the central cavity shown in Fig. 1(c) and the displacement is extracted using COMSOL. The power of the thermal source will be approximated as equal to the absorbed power in the gas (Sec. 4), so we will denote it by $P_{abs}$ as well. The governing equations that are used as well as the implementation in COMSOL can be found in Appendix C. A sweep over the thermal source frequency is performed for a thermal power of $1\mu W$ (using Eq. (15) this corresponds to a THz source of $P_{THz}=1\,mW$ and a methanol concentration of 292 ppm). The displacement at the center of the membrane is extracted and plotted in Fig. 6. Two resonances $f_1=446.5$ kHz and $f_2=834.2$ kHz appear due to the coupling of the mechanical and acoustic oscillators. Figure 7 shows the displacement $\Delta x$ when the source is excited at $f_{mod}=446.5$ kHz and a Parametric Sweep over the thermal power $P_{abs}$ is performed. A linear behaviour can be seen for the considered range of thermal powers, to which the following proportionality is fitted:

 

Fig. 6. The displacement at the center of the poly-si membrane in function of the modulation frequency of the thermal source. Two resonances appear due to the coupling of the acoustic and mechanic resonators.

Download Full Size | PDF

 

Fig. 7. The central displacement of the membrane under excitation at $446.5\,kHz$ in function of the thermal power.

Download Full Size | PDF

4. Results: performance of the sensor

In Sec. 3.1 the optical absorption by the trace gas in the region of interest was determined in function of the concentration (Eq. (15)). In Sec. 3.2 the displacement of the poly-Si membrane in function of a homogeneously applied heat source in the central hole was simulated. There, the loss mechanisms discussed in Sec. 2.3 were included through Eq. (4). The thermal power generated in the acoustic cavity can be approximated as equal to the optical power absorbed $P_{abs}$. This is due to the fast relaxation of rotational excitation into translational energy for gases. This approximation is also used in infrared photoacoustics where the conversion is typically three orders of magnitude slower [8]. The profile of the heat generation reflects the profile of the optical power in Fig. 1(a), but is approximated as homogeneous and confined to the central hole (Fig. 1(c)). This approximation is motivated in Appendix D, where a comparison between the homogeneous and non-homogeneous sources is made and their difference found to be minimal. Combining all these results leads to a function that relates the displacement of the membrane (Fig. 1(c)) to the concentration of methanol $x_{trace}$ and the optical power coupled into the chip $P_{THz}$:

 

Fig. 8. The performance of the sensor according to Eq. (17) (blue) for an optical source of power $P_{THz}=1$ mW. The dashed black line represents the noise floor of a commercial laser Doppler vibrometer with sensitivity $50\,\textrm {fm}/\sqrt {\textrm {Hz}}$ and integration time 1 ms [11]. The detection limit for this source is $234$ ppb.

Download Full Size | PDF

5. Conclusion

A new sensor configuration for terahertz gas spectroscopy is proposed that avoids the need for high-power cooled sources and cooled detectors by employing a transduction and triple enhancement of the signal. The performance for methanol is shown as an example, showing a minimal detectable concentration of $234$ ppb for a source power of $1$ mW and integrating the signal measured by a laser Doppler vibrometer for $1$ ms. The cost of this passive sensor can be minimal by leveraging the already developed CMOS and MEMS fabrication technologies and due to it being made purely out of silicon. This sensitivity, combined with the small size and low price allows for the sensor to be employed ubiquitously for e.g. air quality monitors or quality control in transport of perishables. In the last case we are only interested in the quality at the time of arrival, this means that we can put the sensor inside a transparent package and be able to measure its gas content by scanning it under a THz source and an LDV. One possible application is detecting degradation of packaged food without destroying the controlled environment.

Appendix

A. Evaluation of the acoustic resonator

Here the loss mechanisms of the acoustic resonator and the resulting quality factor $Q_{ac}$ calculated analytically and compared to estimates based on comparable structures in literature. There are three dominant ways through which acoustic energy is lost [25]: firstly sound power will radiate away at the open end, secondly there is viscous drag at the sidewalls. Lastly there are thermal boundary effects, these originate from the fact that the surrounding solid has a much higher heat conductivity than the fluid [26,27].

• • Radiation loss $Q_1$

  The radiation loss for a half-open tube can be calculated through [25]

  (18)$$Q_1=\frac{Hc_s}{\pi R^2 f}\;,$$
  where $f$ is the resonant frequency, $c_s$ is the speed of sound and $R$ is the radius of the cylinder and $H$ its height. This formula shows that the radiation loss can be decreased by decreasing the radius of the cylinder. For $H=90$ $\mu$m and $R=29$ $\mu$m this results in $Q_1=16.5$.

• • Viscosity loss $Q_2$

  Since air particles near the sidewall will adhere to it, we have a zero-velocity condition for the flow of air at these boundaries. The drag this introduces results in the following quality factor function [25]

  (19)$$Q_2=\sqrt{\frac{\rho\pi R^2f}{\eta_d}},$$
  with $\rho =1.225$ kg/m$^3$ the air’s density and $\eta _d=18.3\,\mu$.Pa.s its dynamic viscosity [28]. This quality factor on the other hand is improved when the hole’s radius is increased. Plugging in the known values for this sensor gives $Q_2=11.2$.

• • Thermal effects $Q_3$

  The surrounding walls are made of (poly-)silicon, whose thermal conductivity is three orders of magnitude greater than that of air [26,27]. This allows one to model the walls as perfect heat-sinks at constant temperature. This in turn leads to the creation of a thermal boundary layer which will damp the (thermo-)acoustic eigenmode. Using this approximation the resulting quality factor is [25]

  (20)$$Q_3=\sqrt{\dfrac{\pi R^2f \rho C_p}{k}}\left(\dfrac{C_v}{C_p-C_v}\right),$$
  with $C_p=1005$ J/kg.K and $C_v=718$ J/kg.K the specific heat capacity of air at constant pressure and constant volume respectively and $k=2435$ mW/kg.K its thermal conductivity [17,27,29]. Which gives $Q_3=24.4$. Similar to the viscosity loss, this value can be increased by increasing the radius of the holes.

Increasing the radius hence increases the quality factor relating to the boundary effects ($Q_2$ and $Q_3$), at the cost of more radiation loss ($Q_1$). The result is that one has to settle for a single digit acoustical quality factor: $Q_{ac}=(Q_1^{-1}+Q_2^{-1}+Q_3^{-1})^{-1}=5.2$.

B. Evaluation of the mechanical resonator

Here the simulations which result in the mechanical quality factor in vacuum $Q_{mech}$ (Sec. 2.3) are discussed. All simulations are performed with COMSOL (version 5.3a).

• • Anchor losses

  The loss is numerically calculated using the Solid Mechanics module (wave equation see Eq. (35)) and employing the Eigenfrequency solver with Geometric nonlinearities enabled, as required for plates with non-zero stiffness [30]. The simulation domains are shown in Fig. 9. $D_1$ (red) represents the poly-Si membrane ($R=29$ $\mu$m), while the blue regions $D_2$ and $D_3$ are crystalline silicon half-spheres with radii $300$ $\mu$m and $900$ $\mu$m. $D_3$ is modelled as a perfectly matching layer (PML), providing a model for an open boundary by damping waves coming into it and minimizing reflections at its boundary with $D_2$ by impedance matching [30]. The axial symmetry is exploited by using a 2D Axisymmetric model, greatly reducing the number of mesh cells. For region $D_1$ and $D_2$ a Free Triangular mesh is used, while for the PML $D_3$ a Mapped mesh is defined. This is to reduce discretization errors in the PML [30]. The constraints for the mesh elements in the different domains are summarized in Table 1. The quality factor is extracted from half the ratio of the real part of eigenfrequency to its imaginary part and found to be $Q_1=7.2\times 10^6$. The eigenfrequency is $f_{mech}=706$ kHz, which is close to $f_{ac}$ due to the choice of internal stress of $1.5$ MPa.

• • Thermoelastic damping

  The Thermoelasticity module in COMSOL couples the local stretching (and compression) and heating to each other. It does so by using the local mechanical tension and compression as a heat source in the thermal equations and simultaneously using the temperature from these equations as a thermal load to cause expansion and compression of the solid. An Eigenfruequency analysis was performed (again with Geometric nonlinearities enabled). The simulation region was restricted to domain $D_1$ in Fig. 9 with zero-displacement conditions imposed on the boundary $S$ by defining it as a Fixed Constraint in COMSOL. The mesh for $D_1$ is the same mesh as described in Table 1.

 

Fig. 9. A schematic breakdown of the 2D axi-symmetric mechanical model into separate domains $D_i$ and a boundary $S$. For the simulation of anchor loss region $D_1$ is the poly-Si membrane (red), $D_2$ represents the silicon waveguide. $D_3$ is defined as a perfectly matching layer. In the simulation for thermoelastic damping $D_1$ is modelled with zero-displacement condition at $S$.

Download Full Size | PDF

Table 1. The mesh settings for simulating the anchor loss and eigenfrequency of the mechanical membrane.

View Table | View all tables in this article

C. Acousto-mechanical response

C.1. Acousto-mechanical dynamics

Here the physics behind the displacement of the poly-Si membrane in response to the application of a thermal source in the central cavity are discussed. The model is built by partitioning the region in and around the acousto-mechanical cavity (Fig. 1(c)) in 6 different domains, shown in Fig. 10. Regions $D_1$, $D_3$, $D_4$, $D_5$ and $D_6$ are filled with fluids: air at ambient pressure 1 atm and a concentration of methanol. The thermal power is considered to be homogeneously distributed across $D_1$. The effects of a non-homogeneous source are investigated in Appendix D. In region $D_1$ we use a thermoviscous acoustical model. The heat conductivity of silicon is three orders of magnitude higher than that of air [26,27]. Which is why it is modeled as a heat sink, meaning that boundaries $S_3$ and $S_4$ are defined as isothermal. Operating at room temperature this means

(21)$$T_{\textrm{air}}=20^\circ \textrm{C},$$

at $S_3$ and $S_4$.This leads to the creation of a thermal boundary layer $\delta _{th}$. Likewise we have a viscous boundary layer $\delta _v$ (again close to $S_3$ and $S_4$) due to the constraint of zero acoustic velocity at interfaces with solids (no-slip condition). The size of these layers is frequency dependent, operating at $f_1=446.5$ kHz gives [30]

(22)$$\delta_v=\sqrt{\frac{2\eta_d}{\omega\rho_{\textrm{air}}}}=3.2\,\mu m,\qquad \delta_{th}=\sqrt{\frac{2k}{\omega\rho_{\textrm{air}}C_p}}=3.8\,\mu m$$

with $\rho _{\textrm {air}}=1.225$ kg/m$^3$ the ambient density of air, $\omega$ the angular frequency, $\eta _d$ the dynamic viscosity of air, $k$ its thermal conductivity and $C_p$ its heat capacity at constant volume (the values can be found in Sec. 2.2). Compared to the radius of the cylindrical cavity ($29$ $\mu$m) these boundary layers are not to be neglected. These layers effectively change the dimensions of the acoustic cavity and hence also the eigenfrequencies. The conservation equations for mass, momentum and density for a viscous compressible fluid are (in order)

(23)$$ \frac{d\rho}{dt}+\rho\left(\nabla\cdot\vec{u}\right)=0, $$

(24)$$ \rho\frac{d\vec{u}}{dt}=\nabla\cdot\vec{\vec{\sigma}}, $$

(25)$$ \rho C_p\frac{dT}{dt}-\alpha_p T\frac{dp}{dt}={-}\nabla\cdot\vec{q}+\frac{P_{th}}{V}, $$

(26)$$\textrm{with}\quad \frac{dA}{dt}=\frac{\partial A}{\partial t}+\vec{u}\cdot\vec{\nabla}A, $$

with $\rho$ the local fluid density, $T$ its temperature, $\vec {u}$ the acoustic velocity, $\vec {\vec {\sigma }}$ the fluid’s stress tensor, $\vec {q}$ the heat flux into the domain, $p$ the pressure and $\alpha _p=3.43\times 10^{-3}/K$ the heat expansion coefficient of air [27]. $P_{th}$ is the thermal source which is generated by the absorption of the terahertz wave by the trace gas and finally $V$ is the volume of the cylindrical cavity. The stress tensor contains the material dependant viscosity parameters, namely the dynamic and bulk viscosity coefficients $\eta _d$ and $\eta _b=10\,\mu$Pa.s [31] through the definition

(27)$$\vec{\vec{\sigma}}={-}p\unicode{x1D7D9}+\eta_d\left[\vec\nabla\vec u+(\vec\nabla\vec u)^T\right]-\left(\frac23\eta_d-\eta_b\right)\nabla\cdot\vec u\unicode{x1D7D9},$$

where $\unicode{x1D7D9}$ is the unity matrix. The gas in region $D_3$ and $D_4$ is likewise treated as a thermoviscous fluid, but they don’t contain a thermal source. Further away from the moving membrane ($D_2$ and $D_5$) the gas can be approximated as an isobaric frictionless fluid as waves entering this domain will not reflect back (due to scattering boundary conditions on $S_5$ and $S_6$) and hence their loss is not relevant to the operation of the sensor. This means that the equation of state will be approximated as

(28)$$p=c_s^2\rho,$$

since $\alpha _p=0$ in this model. The fluid’s stress tensor $\vec {\sigma }$ reduces to $-p\unicode{x1D7D9}$ (Eq. (27)) as both the dynamic and bulk viscosity coefficients reduce to zero.Equations (23) and 24 become

(29)$$ \frac{dp}{dt}+p\left(\vec{\nabla}\vec u\right)=0, $$

(30)$$ p\frac{d\vec u}{dt}={-}c_s^2\nabla\cdot p\unicode{x1D7D9}. $$

Assuming harmonic waves of angular frequency $\omega$ and eliminating the acoustic velocity $\vec u$ results in the following wave equation for the pressure field

(31)$$\vec\nabla\left(\vec\nabla p\right)-\left(\frac\omega {c_s}\right)^2 p=0.$$

On boundaries $S_7$ and $S_8$ continuity of normal stress and normal acceleration is enforced through

(32)$$ \vec n\frac{\partial\vec u}{\partial t}=\frac{\vec n}{\rho}\nabla \cdot p\unicode{x1D7D9}, $$

(33)$$ \vec{\vec{\sigma}}\vec{n}={-}p\vec{n}. $$

Furthermore, since the fluid in $D_5$ and $D_6$ is assumed isothermal, the boundary condition for the temperature reads as

(34)$$\vec{n}\left(k\vec{\nabla}T\right)=0.$$

All the variables in Eqs. (32), (33) and (34) are evaluated in domains $D_3$ and $D_4$ and those on the right in $D_5$ and $D_6$. Region $D_2$ contains the poly-Si membrane and is modeled as a linear elastic material. The conservation of momentum equation for harmonic oscillation states

(35)$$\rho_{\textrm{Si}}\omega^2\vec x=\nabla \cdot\vec{\vec S}+ \nabla \vec F+\nabla\cdot\vec{\vec\sigma},$$

where $\vec {x}$ is the displacement vector, $\rho _{Si}=2320$ kg/m$^3$ is the mass density of poly-Si, $\vec {\vec {\sigma }}$ is again the stress tensor but here defined by two material constants: $E$ the Young modulus and $\nu =0.22$ the Poisson modulus. $F$ is the boundary force per area applied to the membrane’s surface, which will be provided by the acoustic pressure wave. $S$ is a $3\times 3$ matrix corresponding to the initial internal tension in each direction. Assuming isotropy, this becomes a diagonal matrix with 1.5 MPa in every entry. As mentioned in Sec. 2.3 the mechanical loss is included into the simulation by defining the Young’s modulus $E$ according to Eq. (4) ($E=(1+jQ_{\textrm {mech}}^{-1})166$ GPa). This has a damping effect since a term linearly proportional to the frequency appears in Eq. (35). The boundaries $S_1$ and $S_2$ allow for the transfer of momentum between the acoustic wave and the mechanical membrane through

(36)$$\vec{u}=\frac{\partial \vec x}{\partial t}.$$

The solid membrane is also modelled as a heat sink at constant temperature, meaning that the temperature variation of the gas at the boundary needs be zero (Eq. (34)). The solution to these wave equations and boundary conditions is found using the numerical solver software COMSOL. The details of this simulation will now be discussed.

 

Fig. 10. A 2D schematic breakdown of the model into separate domains ($D_i$) and boundaries ($S_i$). This region represents a 2D cut of the central hole in the photonic crystal cavity. $D_1$, $D_3$, $D_4$, $D_5$ and $D_6$ are filled with a mixture of air (ambient pressure 1 atm) and the trace gas. $D_2$ is the moving poly-Si membrane. Boundaries $S_5$ and $S_6$ are open boundaries (allowing pressure waves to leave the simulation), $S_3$, $S_4$ are walls with no-slip and isothermal conditions. $S_1$, $S_2$ are the interface between the fluid and solid waves. $S_7$ and $S_8$ are the boundaries between the barotropic frictionless regions ($D_5$, $D_6$) and the thermoviscous fluid. The region in blue represents the surrounding silicon waveguide, but is not part of the simulation. Instead it’s boundaries are set as no-slip and isothermal where it borders regions $D_3$ and $D_4$ and where it borders $D_5$ and $D_6$ the simpler symmetric boundary condition is used as thermal and viscous effects are not of interest here.

Download Full Size | PDF

C.2. Acousto-mechanical simulation

The simulation region is shown in Fig. 10 (not to scale). In fact only half of this 2D surface has to be simulated as a 2D Axisymmetric model is used. The region $D_1$ has a radius of $29$ $\mu$m and height of $90$ $\mu$m. Domains $D_4$ and $D_3$ are half-spheres with radius $200$ $\mu$m as are domains $D_5$ and $D_6$ but with radius $1200$ $\mu$m). Finally domain $D_2$ shares the radius of the cylinder ($29$ $\mu$m) and is $100$ nm thick. A Free Triangular mesh is employed for all domains, it’s feature sizes are presented in Table 2.

Table 2. Mesh size settings for the Free Triangular mesh in the different domains, used in Appendix C. The domains refer to the ones in Fig. 10

View Table | View all tables in this article

Domains $D_{1}$, $D_{4}$ and $D_{3}$ are governed by Eqs. (23), 24 and 25, which are contained in the Thermoviscous Acoustics, Frequency Domain physics module of COMSOL. The Heat Source is modeled as homogeneous across $D_{1}$ with power $P_{th}/V$ with $V$ the volume of $D_{1}$. The boundaries $S_{3}$ and $S_{4}$ are set as Isothermal and No-slip, as are the other edges of the simulation region, since it represents the solid silicon waveguide with much higher thermal conductivity. The fluid in regions $D_{5}$ and $D_{6}$ is modeled as frictionless and barotropic fluid which corresponds to modeling it as a linear elastic fluid in the Pressure Acoustics, Frequency Domain module, with wave equation Eq. (31). The boundaries $S_{5}$ and $S_{6}$ are set as Spherical Wave Radiation, allowing spherical waves to leave the simulation with minimal reflection. The other boundaries are set as Sound Hard Boundary (Wall). Region $D_{2}$ is evaluated using the Solid Mechanics physics module, using an Linear Elastic Material with Initial Stress $S_{M}=2$ MPa. Its boundaries with the edge of the simulation are modeled as Fixed Constraint, forcing the displacement $\vec {x}$ to be zero. Boundaries $S_{7}$ and $S_{8}$ are defined by the Acoustic-Thermoviscous Acoustic Boundary multiphysics module, which enforces the continuity equations Eqs. (32) and 33. Finally boundaries $S_{1}$ and $S_{2}$ use the Thermoviscous Acoustic-Structure Boundary multiphysics module, with Isothermal boundary condition, such that Eqs. (36) and 21 are respected. The Frequency Domain solver is used with all the described physics and multiphysics modules and Include Geometric Nonlinearities enabled.

D. Transduction of optical absorption into heat

The heat source in Sec. 3.2 is considered homogeneous and confined to the central hole. In fact it has the profile in Fig. 11. In this section the difference between the acoustic pressure generated by a heat source as modelled in Sec. 3.2 and one reflecting the profile in Fig. 11 is investigated. The reason for this approximation is that a full 3D simulation is computationally very intensive due to the very thin membrane ($100$ nm) compared to the acoustical wavelength ($480$ $\mu$m for $710$ kHz). The combination of the small mesh size necessary for the membrane and the large volume needed for the acoustic wave propagation makes for a huge number of mesh cells. Instead the simulations in Sec. 3.2 were done in a 2D Axisymmetric model (see Appendix C). This greatly reduces the number of elements and hence allows for a much finer mesh setting, while retaining a large simulation volume. What follows is a comparison in 3D between the acoustic pressure performed onto the membrane by a homogeneous source and one representing the optical cavity profile (Fig. 11). When the membrane is approximated as a fixed wall (infinite thickness) the mesh becomes manageable in 3D. This means that the boundary condition on $S_1$ in Fig. 10 changes to a No-Slip condition. The optical profile was imported from the the COMSOL simulation in Sec. 2.1. This power is renormalized such that the total heat power generated in the cavity is 1 $\mu$W, the same power is used for the homogeneous source. All meshes use a Tetrahedral Mesh and their dimension parameters are summarized in After solving the equations the pressure generated at $S_1$ is integrated over that surface to give the total force generated onto it , plotted in Fig. 12. The ratio of the force in both situations is also plotted, from which it is clear the the difference is minimal: 3% at most and the difference even changes sign around $800$ kHz. Hence we conclude that the optical field distribution may be approximated as homogeneous. The source is also confined to the region $D_1$ in Fig. 10. The motivation here is based on the small optical modal volume (Eq. (2)) and the high confinement factor of the waveguide mode (Sec. 2.1), showing that the field intensity in the surrounding fluid is negligible. All the domains are meshed using Free Tetrahedral elements, their settings can be found in Table 3.

 

Fig. 11. 3D image of the optical field profile in the central hole of the photonic cavity shown in Fig. 1(b).

Download Full Size | PDF

 

Fig. 12. Left: The force generated in the acoustic cavity onto the membrane, which is approximated as rigid, in function of the frequency of the source. In blue when the thermal source has the same profile as the optical field, in green for a homogeneous source. Right: The ratio of both in function of the frequency.

Download Full Size | PDF

Table 3. Mesh size settings for the Free Tetrahedral mesh in the different domains, used in Appendix D. The domains refer to the ones in Fig. 10

View Table | View all tables in this article

Funding

Interreg (1.1.11); Fonds Wetenschappelijk Onderzoek.

Disclosures

The authors declare no conflicts of interest.

References

1. X. Liu, S. Cheng, H. Liu, S. Hu, D. Zhang, and H. Ning, “A survey on gas sensing technology,” Sensors 12(7), 9635–9665 (2012). [CrossRef]  

2. H. Bai and G. Shi, “Gas Sensors Based on Conducting Polymers,” Sensors 7(3), 267–307 (2007). [CrossRef]  

3. K. Arshak, E. Moore, G. Lyons, J. Harris, and S. Clifford, “A review of gas sensors employed in electronic nose applications,” Sens. Rev. 24(2), 181–198 (2004). [CrossRef]  

4. S. Sharma and M. Madou, “A new approach to gas sensing with nanotechnology,” Philos. Trans. R. Soc., A 370(1967), 2448–2473 (2012). [CrossRef]  

5. N. Shimizu, H.-J. Song, Y. Kado, and Y. Muramoto, “Gas Detection Using Terahertz Waves,” NTT Tech. Rev. 7, 1–6 (2009).

6. F. K. Tittel, G. Wysocki, A. Kosterev, and Y. Bakhirkin, “Semiconductor laser based trace gas sensor technology: Recent advances and applications,” in Mid-Infrared Coherent Sources and Applications, M. Ebrahim-Zadeh and I. T. Sorokina, eds. (Springer, Netherlands, 2008), pp. 467–493.

7. P. H. Siegel, “Terahertz Technology,” IEEE Trans. Microwave Theory Tech. 50(3), 910–928 (2002). [CrossRef]  

8. P. Patimisco, G. Scamarcio, F. K. Tittel, and V. Spagnolo, “Quartz-enhanced photoacoustic spectroscopy: a review,” Sensors 14(4), 6165–6206 (2014). [CrossRef]  

9. V. Radisic, S. Member, K. M. K. H. Leong, X. Mei, S. Sarkozy, W. Yoshida, W. R. Deal, and S. Member, “Power Ampli fi cation at 0 . 65 THz Using InP HEMTs,” IEEE Trans. Microwave Theory Tech. 60(3), 724–729 (2012). [CrossRef]  

10. F. Hindle, C. Yang, G. Mouret, A. Cuisset, R. Bocquet, J. F. Lampin, K. Blary, E. Peytavit, T. Akalin, and G. Ducournau, “Recent developments of an opto-electronic THz spectrometer for high-resolution spectroscopy,” Sensors 9(11), 9039–9057 (2009). [CrossRef]  

11. POLYTEC MSA-600 Micro System Analyzer, Preliminary datasheet, Waldbronn, Germany.

12. Y. Li, S. Verstuyft, G. Yurtsever, S. Keyvaninia, G. Roelkens, D. V. Thourhout, and R. Baets, “Heterodyne laser doppler vibrometers integrated on silicon-on-insulator based on serrodyne thermo-optic frequency shifters,” Appl. Opt. 52(10), 2145–2152 (2013). [CrossRef]  

13. J. Dai, J. Zhang, W. Zhang, and D. Grischkowsky, “Terahertz time-domain spectroscopy characterization of the far-infrared absorption and index of refraction of high-resistivity, float-zone silicon,” J. Opt. Soc. Am. B 21(7), 1379–1386 (2004). [CrossRef]  

14. A. Bazin, R. Raj, and F. Raineri, “Design of silica encapsulated high-Q photonic crystal nanobeam cavity,” J. Lightwave Technol. 32(5), 952–958 (2014). [CrossRef]  

15. M. Verstuyft, E. Akiki, G. Ducournau, B. Walter, M. Faucher, E. Mairiaux, J.-F. Lampin, M. Vanwolleghem, and B. Kuyken, “Integrated Photonic Crystal Cavity for Terahertz Waves on a Free-Standing Silicon Platform,” in Optical Terahertz Science and Technology (OTST’19), 10-15 March, Santa Fe, NM, USA, (Santa Fe, United States, 2019).

16. E. Akiki, M. Verstuyft, G. Ducournau, B. Walter, E. Mairiaux, M. Faucher, J. Lampin, B. Kuyken, and M. Vanwolleghem, “A suspended silicon terahertz platform with low loss waveguide and high q photonic crystal cavities,” in 2019 44th International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz), (2019), pp. 1–2.

17. S. W. Rienstra and A. Hirschberg, An Introduction to Acoustics (Eindhoven University of Technology, 2015).

18. H. Hosaka, K. Itao, and S. Kuroda, “Damping characteristics of beam-shaped micro-oscillators,” Sens. Actuators, A 49(1-2), 87–95 (1995). [CrossRef]  

19. J. Rodriguez, S. Chandorkar, D. D. Gerrard, G. M. Glaze, L. Comenecia, Y. Chen, I. B. Flader, and T. W. Kenny, “Direct Measurements of Anchor Damping in MEMS Resonators,” IEEE Sensors pp. 2–4 (2017).

20. M. Gologanu, C. G. Bostan, V. Avramescu, and O. Buiu, “Damping effects in mems resonators,” Proceedings of the International Semiconductor Conference, CAS pp. 67–76 (2012).

21. B. Chouvion, “Vibration Transmission and Support Loss in MEMS Sensors,” Ph.D. thesis, University of Nottingham (2010).

22. A. L. Kimball and D. E. Lovell, “Internal friction in solids,” Phys. Rev. 30(6), 948–959 (1927). [CrossRef]  

23. S. Ghaffari, E. J. Ng, C. H. Ahn, Y. Yang, S. Wang, V. A. Hong, and T. W. Kenny, “Accurate modeling of quality factor behavior of complex silicon MEMS resonators,” J. Microelectromech. Syst. 24(2), 276–288 (2015). [CrossRef]  

24. Spectraplot, “Atomic/Molec. Absorption,” http://www.spectraplot.com/absorption.

25. M. J. Moloney, D. L. Hatten, M. J. Moloney, and D. L. Hatten, “Acoustic quality factor and energy losses in cylindrical pipes Acoustic quality factor and energy losses in cylindrical pipes,” Am. J. Phys. 69(3), 311–314 (2001). [CrossRef]  

26. H. R. Shanks, P. D. Maycock, P. H. Sidles, and G. C. Danielson, “Thermal conductivity of silicon from 300 to 1400k,” Phys. Rev. 130(5), 1743–1748 (1963). [CrossRef]  

27. The Engineering Toolbox, “Thermal Conductivity of selected Materials and Gases,” https://physics.nist.gov/PhysRefData/ASD/levels_form.html.

28. Engineers Edge, “Viscosity of Air, Dynamic and Kinematic,” https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm.

29. M. Daschewski, R. Boehm, J. Prager, M. Kreutzbruck, and A. Harrer, “Physics of thermo-acoustic sound generation,” J. Appl. Phys. 114(11), 114903 (2013). [CrossRef]  

30. COMSOL Multiphysics® User’s Guide v. 5.3a. COMSOL AB, Stockholm, Sweden. 2017.

31. J. Lin, C. Scalo, and L. Hesselink, “Bulk viscosity model for near-equilibrium acoustic wave attenuation,” (2017). http://arxiv.org/abs/1707.05876.