Calibration-free wavelength modulation spectroscopy based on even-order harmonics

This paper proposes a novel and rapid calibration-free wavelength modulation spectroscopy algorithm based on even-order harmonics. The proposed algorithm, analytically deduced from Voigt line-shape function, only involves simple algebraic operations to describe the actual gas absorption spectra, thus eliminating the time-consuming simulations and line-shape fitting procedures adopted in traditional algorithms. Instead of acquiring the entirely scanned absorption line-shape, the proposed technique only requires extraction of the peak values of the harmonics. This characteristic significantly benefits gas diagnosis at elevated pressure and/or temperature, in which the entirely scanned absorption is very difficult to be obtained due to the broadened line-shapes. The proposed algorithm is validated by both numerical simulation and condition-controlled experiment, indicating millisecond-level calculation of gas parameters with the relative error less than 4% in the experiments.


Introduction
Tunable diode laser absorption spectroscopy (TDLAS) has been widely employed to measure gas properties, such as temperature and species concentration, given its advantages of contactless, rapid, and high sensitivity [1][2][3][4][5]. As a modality of the TDLAS technique, wavelength modulation spectroscopy (WMS) can achieve very high sensitivity with strong robustness to background noise, and therefore, has been extensively applied in harsh environments with, for example, strong turbulence, high pressure and/or temperature. Although WMS technique offers many benefits, quantitative WMS measurement is challenging since harmonic signals that are used for calculating the gas parameters strongly depend on the absorption line-shape. Consequently, the measurement results generally need to be calibrated with a known gas mixture, which is impractical in case of (a) poorly-known gas conditions in industrial applications, (b) difficulty in deploying a well-controlled calibration equipment, in particular, with high-temperature and high-pressure gases, and (c) optical inconsistency, e.g. transmission of the optical windows, between the calibration equipment and the target reactors. In the past decade, the above-mentioned challenges have stimulated the development of multiple calibration-free WMS methods. Implementation of the calibrationfree WMS can be mainly categorized as model-dependent and model-independent methods.
The model-dependent method is to calculate the gas parameters using the analytical model of WMS, alongside the well-characterized laser parameters. Li et al. [6] pioneered the development of the model-dependent method by simulating WMS signals as a function of laser parameters and absorption spectra. To eliminate the influence of light intensity fluctuations, Rieker et al. [7] normalized the second-order harmonic signal of WMS by the first-order harmonic signal, i.e. WMS-2f/1f, and further included laser-specific tuning characteristics in the spectral-absorption model to infer gas properties. However, the retrieval accuracy of gas parameters using these methods highly depends on comprehensive characterization of both the line-shape parameters, e.g. collisional broadening coefficients, of the target transition, and the laser parameters, e.g. amplitudes and phases of laser intensity modulation and frequency modulation. In practice, these parameters measured in the laboratory can be very different from those measured on-site, due to the complex gas composition and wide ranges of temperature and pressure.
In contrast, the model-independent method, achieved by waveform fitting of the experimentally measured WMS-2f/1f (or WMS-nf/1f) signal with a simulated one, is more suitable for practical applications since less strict or even no requirement is placed on characterizing the line-shape parameters and laser parameters. This advantage significantly speeds up its application towards in situ diagnosis of reactive flows in harsh environments, for example, at elevated pressure and/or temperature [8][9][10], as the retrieval is independent of the prior information of line broadening. However, the waveform fitting requires time-consuming iterations until the residual is converged to an acceptable level. When the target absorption is interfered by nearby transitions, the result of the multi-parameter fitting cannot be always unique, as the iterations can be converged to the local optimum instead of the global one. To address these issues, Chen et al. [11] recently proposed to use the peak values of zero-, second-, and fourth-order harmonics, enabling the direct and fast measurement without complex simulations and curve fittings. However, Chen's algorithm is derived from the Lorentzian lineshape function, which is not applicable in low-pressure scenarios. Furthermore, the zero-order harmonic signal is ease of noise contamination, leading to reduced fidelity when the measurement is deployed at the presence of environmental noise owing to mechanical vibration, beam steering and window fouling.
Although higher order harmonics have lower signal strength [12,13], they can benefitted from larger signal-to-background ratios (SBRs) [14][15][16][17] and stronger robustness at the extreme conditions where absorption spectra are significantly broadened and large modulation index is required. In this work, we propose, for the first time, a calibration-free WMS algorithm based on even-order harmonics (order  2). The amplitude of each harmonic signal, used for retrieving the gas properties, is analytically derived from the Voigt line-shape function, enabling not only speedy measurement down to milliseconds, but also a general suitability for various degrees of line-shape broadening. The rest of the paper is organized as follows: Section 2 details the fundamentals of the proposed algorithm. Numerical simulation and proof-ofconcept experiment are carried out to validate the proposed algorithm in Section 3 and Section 4, respectively. The paper is finally summarized in Section 5.

Theoretical background
In this subsection, some basics are presented to facilitate the derivation of the proposed method. WMS is implemented by imposing a high-frequency sinusoidal modulation with angular frequency ω on a slowly varying diode laser injection current. As a result, the instantaneous laser frequency v and the laser intensity It can be expressed as: (1 cos( )), where v0 is the laser center frequency, a the modulation depth. Ῑ0 is the averaged laser intensity. ik is the k-th order Fourier coefficient of the laser intensity. ψk is the k-th order phase shift between intensity modulation and frequency modulation. For an isolated absorption line, the modulation index m is defined by 2a/λ, where λ represents the full width at half maximum (FWHM) of the absorption line. According to the Beer-Lambert law, the spectroscopic absorbance can be written as: where It and I0 are the transmitted and incident laser intensities, respectively. A is the integral absorption area. φ(v) is the line-shape function, described by Voigt line-shape function, is a convolution of Lorentzian and Gaussian line-shape functions, noted as φL(v) and φG(v), respectively. φ(v) can be approximated by [18] ( ) ( ) ( ), where the line-shape parameter d is defined as ( ) / ( ).
A dimensionless parameter Δ = 2(v -v0)/λ is introduced to describe the relative offset of the laser center frequency, and therefore, φ(v) can be reformulated as where L(Δ, m) and G(, m) are the time-dependent peak normalized Lorentzian and Gaussian linear functions respectively, and are expressed as [14,19], the analytical expressions of the n-th harmonic of L(, m) can be written as 22 22 [ (1 ) (1 )] ( ) . ., 2 (1 ) where c.c. denotes the conjugate complex number of the preceding term. 0 = 1 and n = 2 for n ≥ 1. Given no offset of the laser center frequency, i.e.,  = 0, the amplitude of each harmonic obtained from Eq. (14) can be analytically expressed by By using infinite series theory [20], the amplitude of each harmonic of G(, m) can be analytically expressed by 22 2 ln 2 / 2 ln 2 / 2 exp( ) ( ) for even , 0 for odd where In/2(z)=i n/2 Jn/2(iz). Jn/2 is the first kind of modified Bessel functions of order n/2. As indicated by Eq. (11), φ(v) is a linear function of L(Δ, m) and G(Δ, m). When the linearity of Fourier transform is employed, the amplitude of the n-th harmonic hn of the given the spectral absorbance α(v) can be expressed by the following linear combination: where pn(m) and qn(m) are defined as  (19) respectively. Dependence of the amplitudes of pn(m) and qn(m) on the modulation index m is shown in Fig. 1. For a given m, both pn(m) and qn(m) decrease as the even harmonic order n increases. Consequently, the first several nonzero even harmonics should be adopted to achieve higher signal-to-noise ratios (SNRs) of the WMS measurement.

Fundamentals of the proposed algorithm
The purpose of the proposed algorithm is to calculate the spectral parameters from the measured amplitudes of nonzero even harmonics. As indicated by Eq. (17), the amplitudes of the first N nonzero even harmonics adopted can be rewritten as:  (21) where E is the N-dimensional identity matrix. Theoretically, with H in hand, the modulation index m can be calculated from the nonlinear equation (21). However, the measurement noise deviates the results of H -H' in Eq. (21) from its desired value of 0. In this work, the modulation index m* is calculated by solving the following minimization problem: When N = 2, M ℝ 2 × 2 is non-singular matrix, and thus the projection matrix, M(M T M) -1 M T , is reduced to the identity matrix, E. In this case, H -H' always equals to 0, and arbitrary modulation index m* can be chose. To make the above minimization problem solvable, the dimension of H must be no less than 3.
With m* in hand, kp and kq in Eq. (20) can be calculated by using the least square method: where M * = (P(m*), Q(m*)). kp/kq is given by: 23 Obviously, kp/kq is only determined by the line-shape parameter d. Fig. 2 shows kp/kq is monotonically dependent on d when -0.85 ≤ d ≤ 0.85. That is to say, for any given d, there is a unique value of kp/kq. It is worth mentioning the line-shape is described as the Voigt function when -0.85 ≤ d ≤ 0.85. The line-shape functions can be characterized by the Lorentzian and the Gaussian functions for 0.85 < d ≤ 1 and -1 ≤ d < -0.85, respectively, resulting into simplified calculation of kp/kq to be detailed in subsection 2.3. Therefore, the FWHM of the Voigt line-shape can be calculated by * 2 / , am   (25) and thus the integrated absorbance area can be calculated by The integrated absorbance areas obtained from two individual transitions can be used to infer the gas temperature using two-line thermometry [21,22]. With the temperature in hand, the gas concentration can be calculated from one of the integrated absorbance areas. The flow chart for calculating gas properties using the proposed algorithm is shown in Fig. 3.

Simplified version for Lorentzian and Gaussian line-shape functions
As noted above, the line-shape function is Lorentzian when 0.85 < d ≤ 1, while that is Gaussian when -1 ≤ d < -0.85. In these two cases, Eq. (20) can be reformulated as where kp = kq = 4A/(πλ). Obviously, H is linear-dependent on P(m) or Q(m). Thus, the modulation index m * can be solved from the following minimization problems: With kp or kq in hand, the line-shape parameters, such as FWHM, integrated absorbance area, and thus the gas parameters can be obtained in the same way described in subsection 2.2.

Parameter optimization
Eq. (20) indicates H can be represented linearly by P(m) and Q(m). The coefficients kp and kq, therefore, are critical to solve the gas properties. According to the theory of linear algebra, the numerical stability for solving kp and kq can be characterized by the condition number of M T M. In case of an inverse problem with larger condition numbers, a small perturbation in H will cause large perturbations in kp and kq. Therefore, it is necessary to choose an appropriate modulation index m* to make kp and kq noise-insensitive. Fig. 4 shows the dependence of the condition number of M T M on modulation index m. It can be seen the condition number decreases as m starts to increase from 0, and reaches a local minimum value when m is around 1.5. Then, condition number increases and tops when m is around 3. After that, the condition number monotonically decreases as m increases. When m  3.5, the condition number is below the previous local minimum value, indicating the everstronger numerical stability for solving kp and kq. However, as indicated in Fig. 1, an excessively large m leads to the reduction of harmonic amplitude. Therefore, the optimal m is selected to be 3.5.

Simulation setup
The proposed algorithm is numerically validated for temperature measurement using water vapor (H2O) transitions. The absorption transitions at 7185.60 cm -1 and 7444.36 cm -1 , that are widely used in two-line thermometry [21,22], are employed in both the simulation and experiment afterwards. In the simulation, the atmospheric-pressure gas concentration is set as 0.05. The absorption length is 20 cm. The gas temperature ranges from 500 K to 2,500 K. The modulation depths for the transitions at 7185.60 and 7444.36 cm -1 are set to 0.131 cm -1 and 0.127 cm -1 , respectively. The scanning depths for the both transitions are set to 0.15 cm -1 . The frequencies of sinusoidal scan and modulation in the simulation are 100 Hz and 10 kHz, respectively. According to the HITRAN database [23], there are several neighbor transitions, with parameters detailed in Table 1, interfering the target transitions. Although each of the transitions appeared as a single absorption line-shape, the combined absorbance shown in Fig.  5 exhibits overlap of absorption line-shapes of nearby transitions with different lower energy levels. The overlap is particularly significant for the transition at 7444.36 cm -1 , covering two interference transitions spaced up to 0.0133 cm -1 . As shown in Fig. 5 (b), the two interference transitions result in significant asymmetry of the target absorbance with additional broadening of the line shape.

. For each transition, the table shows its wavenumbers, line strengths (S(T0)) at T0 = 296 K, pressure broadening-induced halfwidths for itself (ξself) and air (ξair), energies of the lower-state (E), and coefficients of temperature dependent broadening (nair) [23].
line index wavenumber (cm -1 )  With the above-mentioned settings of the WMS simulation, the modulated absorbance signal at the temperature of 1000 K for the two target transitions are obtained, and thus the corresponding second-to tenth-order even harmonics are acquired by the FFT spectral analysis [24]. As shown in Fig. 6, the higher-harmonics give lower amplitude, but have the same order of intensities with the lower-order ones. In the proposed algorithm, only the peak values of the even harmonics are needed, as illustrated in Eqs. (20) -(30), thus significantly simplifying and fastening of data processing in comparison with the traditional waveform fitting algorithm.

Simulation results with noise-free and noise-contaminated data
With the modulation index m* and the amplitudes of the even-order harmonics in hand, we can calculate the integrated absorbance areas for the two transitions according to Eq. (26), and their ratio to infer the gas temperature. We will evaluate retrieved temperature using the proposed algorithm with noise-free and noise-contaminated measurements, respectively.
For the above three cases, Fig. 8 shows the relative error of the retrieved temperature using noise-free data. As the target transitions are interfered by the nearby transitions, minor errors exist in the retrieved temperature even if no noise is imposed on the measurement. For all the three cases, the relative errors are less than 0.5% in the temperature range of 500 -1250 K. When the temperature is above 1250 K, the retrievals degenerate in all the three cases, with the errors in cases 1 and 2 less significant than that in case 3. The worst scenarios at 2500 K give the relative errors of -2.135%, -2.406%, and -3.175% for cases 1, 2, and 3, respectively. Then, white noise is imposed on both the transmitted and incident laser intensities, i.e. It and I0, according to the following formula: where noiset and noise0 are white noise generated by the randn function of MATLAB with a mean value of 0 and a standard deviation of 0.01% × average laser intensity, equivalent to a signal-to-noise ratio (SNR) of 60 dB. This noise level is dominated by the thermal noise from the detector, signal amplification, and digitization circuits [25]. To characterize the noise performance of the proposed algorithm, the average relative error and the relative standard deviation (RSD) of the temperature retrieved from 100 repetitive simulations are shown in Figs. 9 (a) and (b), respectively. The average relative error shown in Fig. 9 (a) indicates similar results as those using the noise-free data. As the temperature sensitivities of the chosen H2O transitions decrease as the temperature increases [22,26], the retrieved temperature suffers from more significant fluctuations, indicated by the larger values of RSDs shown in Fig. 9 (b), at higher temperatures. In addition, involvement of more harmonics contributes to lower RSDs, and thus stronger robustness of the proposed algorithm.

Experimental setup
In this section, experiments, with setup shown in Fig. 10, were carried out to validate the proposed algorithm. The quartz optical cell has three zones, all fitted with wedged windows to suppress etalon interference. The two outer zones were purged with dry air to prevent unwanted H2O absorption along the optical path, while the center zone, with the path length of 20 cm, was filled with a uniform gas mixture at a controlled temperature between 773 and 1273 K. The gas mixture is obtained from the dry air going through H2O liquid, and therefore, contains the target H2O molecules for the absorption measurement. Three K-type thermocouples were equally spaced along the optical cell to measure the temperature of the target gas. The optics, detector and the intermediate open paths were purged with N2 to remove interfering H2O absorption in the room air. Two distributed feedback (DFB) diode lasers operating near 7185.60 cm -1 (NTT, NLK1E5GAAA) and 7444.36 cm -1 (NTT, NLK1B5EAAA) were used to probe the H2O transitions. Each laser was controlled independently by a laser controller (Stanford Research Systems, LDC501), which adjusts the temperature and the injection current of the laser. To achieve wavelength modulation, the injection current was driven by a function generator (RIGOL, DG1062Z) which provided a 100 Hz sinusoidal scan signal superimposed with a 10 kHz sinusoidal modulation signal. The two lasers were operated in time division multiplexing (TDM) mode [27]. To monitor DC background signals such as the dark current of the detector and thermal radiation, there is a 10 ms non-laser interval between the neighboring wavelength scan of the two lasers. The two lasers were combined and then split by a 2×2 single-mode fiber-coupler, with one output being guided through an etalon with a free spectral range (FSR) of 0.01 cm −1 for wavelength monitoring and the other output penetrating the target gas cell. Each laser beam was detected by a photodetector (Thorlabs, PDA10CS-EC) and further sampled by a 14-bit data acquisition card (NI Corporation, PXIe-5170R). The digitization rates for sampling the etalon signal and the absorption signal are 25 Mega samples per second (Msps) and 1 Msps, respectively.
The optical cell was firstly heated to the desired temperature and filled with dry air to obtain the absorption-free transmission. Then, the center zone was evacuated and filled with the H2O/air mixture. When the adsorption and transient effects was stabilized, the absorption transmissions were recorded at the given furnace temperature. The temperature of the optical cell was then adjusted to repeat the entire measurement process. In the experiment, the flow rates for both the purging gas and the target gas mixture were small and stabilized to reduce the influence of convective heat transfer on the temperature set by the furnace. Therefore, two mass flow meters were used to control both gas flow rates at 100 standard cubic centimeter per minute (sccm). Since the optical cell is open to the air, the gas pressure inside the cell is same as the atmospheric pressure, that is 0.998 atm, measured by the pressure gauge (MKS, type660).

Experimental results and discussions
Figs. 11 (a) and (b) shows the sampled background transmission I0(t) and the absorption transmission It(t) at 973 K, respectively. Both I0(t) and It(t) are averaged for 100 repetitive measurements to reduce the interference of measurement noise. Although the background transmission is not an ideal absorption-free signal, caused by the small purging flow of the dry air, logarithmic transformation in Eq. (3) neutralizes any impact of the background absorption on the temperature retrieval. In addition, the rising edges of both the second periods of I0(t) and It(t) are adopted in the harmonic analysis in the experiment, since the first period, just after switching the lasers, is under stabilization. Subsequently, the even-order harmonics demodulated by FFT spectral analysis are shown in Figs. 12 (a) and (b) for the transitions at 7185.60 cm -1 and 7444.36 cm -1 , respectively. Finally, the amplitude of each harmonic is extracted by conic fitting of the partial waveform nearby the peak.  With the optimal modulation index determined in subsection 3.2 and the extracted amplitudes of the even-order harmonics, the temperature is retrieved using the proposed algorithm by following Eqs. (20)-(30). To examine the performance of the proposed algorithm, the retrieved temperature values are compared with those measured by the thermocouple in the range of 773 -1273 K with an interval of 100 K. Similar as the simulation, the three cases with different combinations of the even-order harmonics are verified in the experiment. As shown in Fig. 13 (a), the temperature retrieved in all the three cases agree well with that measured by the thermocouple. The maximum relative errors, characterized in Fig. 13 (b), are 2.59%, 3.60% and 2.82% for cases 1, 2 and 3, respectively. To evaluate the superiority of the proposed algorithm in terms of computational cost, the waveform fitting algorithm [28] was also employed to analyze the same acquired data. Fig. 14 shows the measured and fitted 1f-normalized WMS-2f (WMS-2f/1f) waveforms with the fitting residuals for the transitions of 7185.60 cm -1 and 7444.36 cm -1 at 973 K. The waveform fitting is implemented using the samples on the wavenumber free from the interference transitions. Therefore, 4000 samples on the left wing of absorbance are used for fitting for the transition 7185.60 cm -1 (Fig. 14 (a)), while 2800 samples of the right wing are used for fitting for the transition 7444.36 cm −1 .
It is worth noting the WMS-2f/1f waveform fitting algorithm is sensitive to the initial value of parameters to be fitted. Improper initial values may increase the computational time, and even make the solution fall into the local optimum. In this work, the furnace temperature given by thermocouples and other needed line-shape parameters retrieved by the proposed algorithm were taken as the initial values of the fitting algorithm. In contrast, the proposed algorithm is always converged with no restriction attached to the initial value. All the computational process in this section were conducted by a computer with an Intel Core i7-1165G7 @ 2.80GHz CPU.
The computational cost results for these two algorithms are summarized in Table 2. The proposed algorithm in the temperature range of 773 -1273 K gives the average computational time of 0.137 s, 0.142 s, and 0.144 s for cases 1, 2, and 3, respectively. Although the appropriate initial values were given, the average computational time of the WMS-2f/1f waveform fitting algorithm still reaches the average time of 22.902 s. That is to say, the computational efficiency of the proposed algorithm is improved by at least two orders of magnitude.

Conclusions
A novel calibration-free WMS algorithm, analytically deduced from the Voigt line-shape function, is proposed for accurate retrieval of gas properties using the amplitudes of the evenorder harmonics. As only the peak values of the harmonics are needed, the proposed algorithm substantially benefits WMS measurement with large line-shape broadening, for example, measurements in high-pressure and high-temperature shock-wave wind tunnels. For the two pre-selected H2O transitions, numerical simulation was firstly carried out to determine the optimal modulation index for different combinations of the even-order harmonics, and examine their accuracy and noise resistance. Simulation results indicate the relative errors are less than 0.5% in the temperature range of 500 -1250 K and stronger robustness when more harmonics are involved. Furthermore, experimental validation of proposed algorithm shows a maximum relative error of 3.6% in the temperature range of 773 -1273 K. Compared with the traditional WMS-2f/1f waveform fitting algorithm, the computational efficiency achieved by the proposed algorithm is improved by at least two orders of magnitude.

Funding. National Key Research and Development Program of China (2017YFB0603204); Global Challenges
Research Fund (LMIC Travel and Partnerships Fund).

Disclosures. The authors declare no conflicts of interest.
Data availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.