Diffuse terahertz spectroscopy in turbid media using a wavelet-based bimodality spectral analysis

Current terahertz (THz) spectroscopy techniques only use the coherent light beam for spectral imaging. In the presence of electromagnetic scattering, however, the scattering-mitigated incoherent beams allow for flexible emitter-detector geometries, which enable applications such as seeing through turbid media. Despite this potential, THz spectroscopy using diffuse waves has not been demonstrated. The main obstacles are the very poor signal to noise ratios of the diffused fields and the resonance-like spectral artifacts due to multiple Mie scattering events that obscure the material absorption signatures. In this work, we demonstrate diffuse THz spectroscopy of a heterogeneous sample through turbid media using a novel technique based on the wavelet multiresolution analysis and the bimodality coefficient spectrum, which we define here for the first time using the skewness and kurtosis of the spectral images. The proposed method yields broadband and simultaneous material characterization at detection angles as high as 90° with respect to the incident beam. We determined the accuracy of the wavelet-based diffuse spectroscopy at oblique detection angles, by evaluating the area under the receiver operating characteristic curves, to be higher than 95%. This technique is agnostic to any a priori information on the spectral signatures of the sample materials or the characteristics of the scattering medium, and can be expanded for other broadband spectroscopic modalities.

2n + 1 n(n + 1) (a n π n + b n τ n ), and 2n + 1 n(n + 1) (a n τ n + b n π n ), representing the parallel and perpendicular polarizations of the scattered field, respectively. Here, a n and b n are the complex coefficients of Riccati-Bessel functions and depend on the complex refractive index ratio and the size parameter x = 2πr/λ . Additionally, π n = p 1 n /sinθ and τ n = dp 1 n /dθ , where p 1 n is associated with Legendre polynomials, represent the Mie angular functions. We calculated the angular scattered intensity using the Mie scattering amplitude functions for a LDPE particle with r = 180 µm at frequencies in the range 0.2 -1 THz. Figure 1(a) illustrates the angular scattering patterns at frequencies f = 0.2, 0.4, 0.6, 0.8, and 1 THz, corresponding to size parameters x = 0.6, 1.3, 1.9, 2.5, and 3.1, respectively, which are all located in the Mie regime. It can be noted that the lowest frequency, 0.2 THz, yields the highest amounts of backscattering and scattering at oblique angles θ = 0 • , while as frequency increases, most of the radiation is scattered in the forward direction. This corroborates our measurements, where at higher detection angles, we could only extract resonances located at the lower frequencies. Figure 1(b) illustrates the scattered intensity 1 2 (|S 1 (θ )| 2 + |S 2 (θ )| 2 ) between θ = −90 • to 90 • , demonstrating that a particle-size scattering scenario results in different amounts of scattered intensity at different frequency components present in a broadband THz wave.
At high concentrations of scatterer particles (> 10%), multiple scattering becomes the dominant process. In other words, the incident field on each particle represents a random superposition of the primary incident wave and the waves scattered by other particles in the medium. Effective medium theories incorporating the Mie formalism have been developed to predict the higher-order scattering contributions [3][4][5] . Here, we demonstrate the implementation of a self-consistent effective medium theory 4-7 for describing the frequency-dependent scattering loss by non-absorbent granular particles. Using the dynamic effective medium theory of Stroud and Pan 7 , an effective dielectric constant ε e f f of a two-component composite material can be derived using an iterative scheme given by, where  1. (a) A polar plot representing the angular scattering patterns generated by a LDPE particle with r = 180 µm at frequencies f = 0.2, 0.4, 0.6, 0.8, and 1 THz, corresponding to size parameters x = 0.6, 1.3, 1.9, 2.5,and 3.1, respectively, which are all located in the Mie regime. (b) The scattered intensity, 1 2 (|S 1 (θ )| 2 + |S 2 (θ )| 2 ), for the LDPE particle at different size parameter values, demonstrating that the same particle generates different amounts of scattered intensity at different frequency components present in a broadband THz wave. and with V representing the volume fraction of the scatterers in the medium, λ representing the incident beam's wavelength, ε 1 representing the scatterers' dielectric constant, ε 2 representing the medium's dielectric constant, n is the number of Mie partial waves, a n and b n are the complex coefficients of Riccati-Bessel functions, and ρ(r) represents the scatterer's size distribution. By assuming a fixed particle size, the size distribution can be written as, Because LDPE is non-absorbent at THz frequencies, this attenuation in the pulse can be attributed to scattering. It is also evident that the model accurately predicts the scattering-induced attenuation in the measurement.

2/9
where δ (r) indicates the Kronecker delta function. The first approximation for this iterative implementation is obtained from the ordinary Bruggeman mixing rule given by, As an example, Fig. 2 illustrates the attenuation coefficient measured from loosely-packed LDPE powders with r = 25 µm and volume density V = 0.47. Figure 2(a) illustrates the spectrum of a signal measured through ambient air, the reference measurement, and a signal measured through the LDPE powders. Figure 2(b) compares the measured attenuation coefficient, shown by the solid balck line with squared blue markers, with the attenuation coefficient obtained from the Stroud and Pan's dynamic effective medium model, shown by the solid red line. Because LDPE is non-absorbent at THz frequencies, the pulse attenuation can be attributed to scattering. It is also evident that the red line derived from the model can precisely predict the scattering-induced loss in the measurement.
For the more complicated scenario of having a pressed pellet with smaller particle size and larger scatterer density hidden beneath a turbid medium (0.3 volume fraction) with larger particle size, i.e., the test scenario studied in our manuscript, the effects of scattering by different media are convoluted. However, the discussed effective medium model can still predict the scattering-induced attenuation originated from the two media. Figure 3 shows the attenuation coefficients measured from four example HDPE pixels in the sample. Along each measured attenuation coefficient, shown by the blue lines, the red solid line demonstrates the attenuation coefficient derived from the effective medium model. It can be noted that the model accurately predicts the aggregate scattering attenuation caused by the two media.

Particle size
The particle size and scatterer density are not uniform across the sample. Hence, the mean value might not be representative of the correct particle size at different locations. However, Stroud and Pan's effective medium model can still be used to find (c) (d) Figure 3. (a-d) Blue lines illustrate the attenuation coefficients measured through the HDPE pellet and LDPE particles from four random pixels selected from HDPE region of the sample, along with the attenuation coefficients calculated using the Stroud and Pan's effective medium model, shown by the red lines. Because both HDPE and LDPE are non-absorbent at THz frequencies, this high attenuation can be attributed to the Mie scattering. The average particle radius that yielded the minimum mean squared error between the measurement-and model-derived attenuation in each pixel is given in (a-d).
a more accurate particle size in each pixel by minimizing the mean squared error between the model and the measurement. Figure 3 demonstrates four examples of such pixels with the particle size derived from the model. It can be noticed that the radius (i.e., particle size/2) at each pixel is slightly different. Additionally, Fig. 4 illustrates the size parameter x = 2πnr/λ , which governs the scattering regime (x between 0.2 and 20 indicates the Mie regime 2 ), for the size values obtained from the effective medium theory in Fig. 3. It confirms that across the measurements' useful bandwidth of 0.1 -2 THz, these particles scatter according to the Mie theory.

Impact of the spectral properties of the resonant modes and noise on diffuse THz spectroscopy
We have shown that using the thresholded local maxima in the bimodality spectrum of the THz spectral images, the signature resonant frequencies of all the sample constituent chemicals can be extracted. Moreover, by applying this technique to the MODWT MRA of the THz extinction spectra, we recovered the resonant frequencies of the α-lactose monohydrate and PABA from diffuse THz waves at very low signal to noise ratios (SNR). Here, we further investigate the robustness of this technique in identifying the characteristic resonant frequencies from two overlapping resonant signatures. We also numerically explore the effects of the dielectric resonance properties, including the absorption mode's height and full width at half maximum (FWHM), on the performance of the bimodality-based diffuse spectroscopy technique. To study these effects, we used the α-lactose's dielectric function in the range 0.2 -1.7 THz shown in Fig. 5(a). This dielectric function was acquired from a sample pellet purely made from α-lactose monohydrate. Additionally, using the Lorentz oscillator model of the classical dispersion theory, we created artificial resonant signatures in the proximity of the α-lactose's resonance at 0.53 THz. The Lorentz oscillator model describes a medium as a compound of atomic dipole oscillators 8,9 , and thereby the absorption lines in the dielectric function represent the uncoupled resonant modes of these oscillators. Accordingly, the dielectric function of a material with K resonant modes is given by where ε ∞ is the dielectric constant at the higher-frequency limit of the bandwidth, and f k is the kth resonant frequency, while ∆ε k and γ k represent the height and FWHM of the resonant signature at frequency f k . The solid lines in Fig. 5(a) illustrate the Lorentz model fitted to the α-lactose's dielectric function using a 10-element parameter vector given by where the subscripts 1, 2, and 3 correspond to the resonant modes at 0.53, 1.2, and 1.38 THz, respectively. Therefore, by excluding the parameters associated with the resonant modes at 1. resonant signatures with desired properties can be created. For example, Fig. 5(b) illustrates two extinction spectra resultant from including only the 0.53 THz resonant signature and shifting its center position to 0.53+2δ f THz in the α-lactose's Lorentz model, where δ f = 0.015 THz represents the sampling period of the ε( f ). Therefore, we have simulated the dielectric function of an artifical material having a vibrational modes at 0.56 THz with variable resonance height and FWHM. Next, we simulate the THz extinction spectra of a 2 mm-thick sample comprised of both lactose and the new artificially modeled material in the presence of various noise scenarios. Fig. 5(c) shows the simulated THz time-domain signals corresponding to the extinction spectra of the two materials depicted in the sample, illustrated in Fig. 5(d). The left half of the sample is composed of the material with a resonant frequency at 0.53 THz, while its right half is made from the material with a resonance at 0.56 THz, and the surroundings are made from a material similar to polyethylene with no resonant frequencies, corresponding to the black lines in Figs. 5(b-c). Furthermore, using additive colored-noise models, we recreated wavelength-dependent low SNR scenarios often encountered in diffuse spectroscopy measurements. The color of a random noise model depends on the characteristics of its power spectral density (PSD), which is proportional to 1/ f β per unit of the bandwidth, where β = 0 produces white noise and β = 1 yields a pink noise model 10 . Because the scattering effects depend on the wavelength of the illuminating light, using a colored noise model with a frequency-dependent PSD, such as pink noise, better represents the experimental conditions described in this paper. Figures 6(a) and 6(c) compare the time-domain THz signals contaminated with additive pink noise at SNR = 10 dB and SNR = 0 dB, respectively. Note that although the model of the noise remains the same, the added noise is random and different for each pixel. Figures 6(b) and 6(d) illustrate the extinction spectra corresponding to these time-domain signals. It can be seen that the identification of the resonant features at 0.53 and 0.56 THz among other spectral artifacts is challenging, especially at SNR = 0 dB. At SNR = 0 dB, the noise contribution is dominant, which results in obscured resonant signatures. For example, the resonant feature at 0.56 THz is obscured in the green trace in Fig. 6(d). This corresponds with the measurements taken at the higher detection angles discussed in the main manuscript. Figures 6(e-f) illustrate that the proposed broadband bimodality spectrum of the sample shown in Fig. 5(d), when it is composed of the noise-contaminated THz signals, can still reveal the location of the resonant frequencies, even at a very low SNR. For example, Fig. 6(e) compares the bimodality spectrum when the time-domain signals are corrupted with white noise at SNR = 0, 10, and 20 dB. In all three cases, the local maxima above the 5/9 threshold correspond to the resonant frequencies at 0.53 and 0.56 THz, while the bimodality at all other frequencies is below 5/9, even at SNR = 0 dB. Moreover, with decreasing SNR the bimodality at resonant frequencies also decreases, which is in accordance with the results obtained in the oblique measurement geometries in the main manuscript. Additionally, Fig. 6(f) illustrate the bimodality spectrum when signals are contaminated with additive pink noise at different SNRs. Although at SNR = 0 dB the local maxima at resonant frequencies are very close to the 5/9 threshold, they are still capable of revealing the resonant frequencies. Therefore, results in Fig. 6(e-f), considering the destructive noise effects on resonant signatures at lower SNRs shown in Figs. 6(b), confirm the robustness of the bimodality-based material detection.
Here, we address the problem of having more than two materials in the imaging target. We show the results obtained for the identification of three substances with different resonant frequencies. The resonances are simulated using the α-lactose monohydrate's Lorentz oscillator parameters at 0.53 THz. The first substance has a resonant frequency at 0.53 THz, similar to α-lactose, as shown in the dielectric function in Fig. 7(a). The dielectric function of the second and third substances are formed by shifting the resonant frequency to 0.6 and 0.75 THz, while the height and width of the resonant features (i.e., ∆ε k and γ k in Eq. (8)) are changed. Figure 7(b-c) show the refractive indices of the two new substances. Figure 7(d) shows an artificial sample made from these three substances, where the circle represents a pellet embedded in the background material.  8)) are changed, (d) an artificial sample pellet made from substances in (a-c), the black, cyan, and orange pixels designate the resonances at 0.53, 0.6, and 0.75 THz, whereas the background red pixels do not contain any resonances, (e) an example extinction coefficient from each region of the modeled sample, random pink noise at SNR = 20 dB is added to the extinction coefficients to replicate the frequency-dependent scattering artifacts, (f) the bimodality coefficient spectrum of the sample, where the frequency of the three local maxima above the 5/9 threshold, marked using black, cyan, and orange asterisks, correspond to the resonant frequencies at 0.53, 0.6, and 0.75 THz.
The black, cyan, and orange pixels designate the resonances at 0.53, 0.6, and 0.75 THz, respectively. The background red pixels do not exhibit any spectral resonances. We also added random pink noise to the extinction coefficients to artificially replicate the frequency-dependent scattering artifacts. Figure 7(e) shows an example extinction coefficient from each region. Figure 7(f) shows the bimodality coefficient spectrum of this sample. It can be noted that the frequency of the three local maxima above the 5/9 threshold, marked using black, cyan, and orange asterisks, correspond to the three simulated resonant frequencies.
Finally, we further investigate the effect of the height, ∆ε k , and FWHM, γ k , of a resonant mode on the bimodality spectrum. With this aim, we varied the ∆ε k and γ k of the resonant signature at 0.56 THz from half to two times the value of the same  Fig. 8 indicate the number of distinct vibrational modes identified in the overlapping bimodality spectra under each simulation scenario. These results indicate that the success of the proposed technique can be limited by the spectral proximity and overlap of the two very close resonant modes. However, at reasonable experimental noise levels (better than SNR = 0 dB) the new computational approach almost always identified two distinct but overlapping vibrational modes.