Open Access
Add to BibSonomy
Abstract
Terahertz imaging is becoming a biological imaging modality in its own right, alongside the more mature infrared and X-ray techniques. Nevertheless, extraction of hyperspectral, biometric information of samples is limited by experimental challenges. Terahertz time domain spectroscopy reflection measurements demand highly precise alignment and suffer from limitations of the sample thickness. In this work, a novel hybrid Kramers-Kronig and Fabry-Pérot based algorithm has been developed to overcome these challenges. While its application is demonstrated through dielectric retrieval of glass-backed human bone slices for prospective characterisation of metastatic defects or osteoporosis, the generality of the algorithm offers itself to wider application towards biological materials.
Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
1. Introduction
Terahertz radiation grants us access to a unique view of the world. Superior resolution over microwave imaging, alongside non-ionising effects compared to X-rays, has resulted in the exploitation of terahertz radiation for applications [1] such as security screening [2,3] and quality control [4]. Particular interest has been focused on biological imaging [5], where terahertz energies match low frequency vibrations of biomolecules [6] and rotational modes of liquid water molecules, which can be utilised for sensing and imaging [7–9]. Such applications include characterisation of human bone metastatic defects or osteoporosis and can be extended to prostheses design [10–13]. The development of this powerful technique, however, remains limited by challenges in extraction of meaningful biometric information from measurements.
The typically high water content of biological materials, alongside its highly absorbing properties at terahertz frequencies, make reflection spectroscopic measurements far more practical than transmission. Accurate retrieval of both the amplitude and phase of the reflected radiation relies on a number of factors. The two critical factors, for time domain spectroscopy (TDS) in particular, are accurate alignment of the sample and reference planes [14,15] and the ability to resolve or model any Fabry-Pérot reflections that occur within the sample [16–19]. Here, we develop a hybrid retrieval algorithm to account for possible experimental misalignment and unresolvable reflections in thin samples through application of the Kramers-Kronig relations and modelling Fabry-Pérot effects. Our methodology improves upon existing algorithms incorporating Fabry-Pérot modelling for reflection measurements [16–19] by developing a self-reliant method of estimating sample thickness, without the need of an initial measurement and reducing computational complexity through two-dimensional scanning of dielectric properties in the frequency domain. Additionally, the generality of the approach means that its application is not limited by material thickness, absorption and reflectivity. Application of the algorithm is demonstrated on glass-backed heterotopic ossification (HO) human bone slices for prospective characterisation of disease. Needless to say, the algorithm can be applied for wider range of applications, but this is not within the scope of this paper.
2. Methods
2.1 Kramers-Kronig phase retrieval
TDS in reflection geometry is a common technique used to measure the interaction of a terahertz electric field with a sample upon reflection. Depending on the dielectric environment of the sample, the amplitude and phase of the field is modified. Practically, retrieval of the phase is challenging. Phase is highly sensitive to the sample position on the micron scale [20,21]. Any offset between a reference and sample plane $\Delta L$ induces an accumulated phase offset of $\Delta \theta (\omega )=\omega n (2\Delta L/c)$, where $n$ is the refractive index of air and $\omega$ and $c$ are the angular frequency and vacuum speed of the wave, respectively. Difficulties in precise positioning, despite experimental efforts [20–22], mean researchers have turned to post processing numerical methods to overcome this challenge, with particular focus on the principle of causality and its expression in the Kramers-Kronig relations [14,23,24].
The Kramers-Kronig relations are mathematical bidirectional equations that describe a fundamental relation between the real and imaginary parts of a linear complex function, for which the principle of causality stands. That is, the response of a system is zero until an impulse is applied to it. The refractive index of a material is a complex function which can be expressed in this way, where the real and imaginary parts describe the energy storage due to material polarisation and absorption, respectively. Hence, from measurements of the dispersion of a field propagating through a material, one can obtain the absorption, and vice versa. This fundamental relationship has found wide application in linear [25–27], and subsequently non-linear, optics [28–30]. A limitation of the Kramers-Kronig relations is that the knowledge of the spectrum must be across a semi-infinite frequency range. This limitation can be relaxed through application of Singly Subtractive Kramers-Kronig (SSKK) relations, for improved convergence [29–31].
This work implements an algorithm developed by V. Lucarini and colleagues [30], through application of SSKK relations to the real and imaginary parts of the reflected field defined as $r(\omega )=|r(\omega )|\exp {(i\theta (\omega ))}$, or its logarithmic amplitude $\ln {|r(\omega )|}$ and phase $\theta (\omega )$. The dispersion integrals are given by
The optimisation approach, presented in Fig. 1, involves the iterative application of the SSKK relations to the amplitude $\ln {|r(\omega )|}$ and a corrected phase $\theta (\omega )=\theta (\omega )_{measured}+\Delta \theta (\omega )$ for a range of guessed phase offsets $\Delta \theta (\omega )$. For each iteration, the degree of consistency between the two relations is calculated using the L2 norm $\|r(\omega )_{guess} - r(\omega )_{corrected}\|_2$. The application is repeated until the consistency has an incremental improvement less than $\epsilon$ (=0.01). Subsequently, the self-consistent reflectivity is calculated, from which the phase offset is extracted. Convergence is typically achieved after 10 iterations.
Figure 2 demonstrates the application of the algorithm to TDS measurements of paraffin wax with an offset of $\sim$ 200 $\mu$m between the sample and reference planes. Prior to the phase offset retrieval, one can observe a large change in phase with frequency in the left panel, attributed to the additional phase acquired over distance $\Delta L$. This results in highly distorted extracted dielectric properties, presented in the right panel. The corrected phase acquired from the algorithm is shown to have a much lower frequency dependence, enabling reliable extraction of the refractive index and extinction coefficient, consistent with values reported in the literature [19], illustrated by the transparent lines in Fig. 2(right).
Fig. 2. Phase response (left) and extracted dielectric properties (right) of paraffin wax obtained with a sample-reference plane offset of $\sim$200 $\mu$m, before and after the Kramers-Kronig retrieval algorithm is applied. The transparent lines in the dielectric results indicate the literature values [19]. The inset shows an illustration of the offset of the reference plane to the sample plane. $E_0$, $E_{s}$ and $E_{r}$ are the incident, sample and reference electric field amplitudes, respectively.
The extracted phase offset was found to be robust, independent of the anchor frequency selection, given that they are within the measured frequency range. The uncertainty of the output phase is largely dependent on the estimated phase offset increment, while the run time of the algorithm is limited by both the number of input estimated offsets and anchor frequencies. Calculations of errors of the Kramers-Kronig algorithm can be found in Ref. [23]. It should also be noted that while the algorithm is applied to cases of normal incidence, its application can be extended to oblique incidence, accounting for the additional path length travelled by the field due to the additional angle $\theta$, where the effective path length is given by $L_{\textit {eff}}=L/\cos {\theta }$.
For the case of optically thick samples, the fundamental causal basis of the Kramers-Kronig retrieval algorithm means it can be used for reliable extraction of dielectric properties for a vast range of samples, without prior knowledge of the dielectric response. For optically thin samples, however, where pulses overlap, Fabry-Pérot reflections come into play.
2.2 Fabry-Pérot dielectric property retrieval
Fabry-Pérot reflections of electromagnetic fields occur in cavities made from two parallel reflecting surfaces formed from refractive index boundaries. They are often exploited in spectroscopy to increase field-sample interactions. For TDS measurements of thin samples, however, they can become problematic. For pulsed excitation, individual reflected pulses will exit the cavity at different times depending on the cavity thickness. For thicker cavities, the reflections will be spread out further in time. The ability to resolve the reflected pulses can be crucial in TDS for reliable dielectric property extraction. It is worth noting that difficulties in pulse resolution extend to highly dispersive materials, where part of the pulse may overlap itself after a reflection.
Figure 3 demonstrates the effect of Fabry-Pérot reflections on the retrieved spectral response of the material under normal incidence conditions. The reflected time and frequency responses were simulated using the transient solver of CST Microwave Studio, which uses the finite integration technique to solve Maxwell’s equations, for structures (illustrated in the inset) of different thicknesses. The structures were designed to have dielectric properties comparable to the glass-backed HO bone structures the algorithm is subsequently applied to. For the optically thick sample, no Fabry-Pérot reflections are measured, since the first reflection occurs after the measurement window, and hence the true spectral response is extracted. When reflections are recorded in the waveform by reducing the thickness of the sample, oscillations emerge in the spectral response, with decreasing frequency for decreasing thicknesses. This, in turn, induces artifacts in the extracted dielectric properties.
Fig. 3. Time domain (left) and frequency domain (right) response of optically thick, thin and intermediate structures under normal incidence, simulated using the transient solver of CST Microwave Studio. Layer 1 was modelled as air, while layers 2 and 3 were modelled as bone and glass, with refractive index values of 1.8 and 2.6 and extinction coefficient values of 0.08 and 0.15, respectively. The temporal waveforms are offset along the $y$-axis with respect to each other by 0.2 a.u. for clarity. The inset illustrates the sample configuration, for which the sample (layer 2) sits on a reference layer (layer 3). $E_0$, $E_{s}$, $E_{r}$ and $E_{\textit {FP}}$ are the incident, sample, reference and Fabry-Pérot electric fields, respectively.
For the intermediate case, the Fabry-Pérot reflection at $\sim$13 ps can easily be resolved from the initial reflection, and hence the resultant artifacts removed by standard time-windowing. For the optically thin sample, however, the initial and subsequent reflections overlap completely. Separation of the terahertz pulses is a well known challenge in the TDS terahertz community for transmission measurement, where researchers have focused efforts on overcoming the problem through methods such as modelling the reflections [16,18,32–35]. These methods can be translated to overcome the same challenges presented in reflection measurements.
Here, an algorithm based on frequency domain modelling of the reflected field from TDS measurements alternative to those reported in the literature [16,18,32–35] is presented. The algorithm offers reduced computational complexity and relaxation of limits in sample thickness, absorption and reflectivity. The field is assumed to be a plane wave with linear polarisation. As mentioned, to extract the dielectric properties of the sample under investigation, TDS measurements typically involve a measurement of the reflected field for both the reference sample with known dielectric response and the investigated sample. Here, as illustrated in the inset of Fig. 3, the field is modelled for a thin sample (layer 2) fixed to a reference sample (layer 3). The versatility of this approach means that it can be applied to a range of configurations through appropriate field modelling. From the Fresnel coefficients, the sample and reference fields are given by
Providing the dielectric response of the reference material is known, the properties of the sample can be retrieved through a comparative approach. Here, medium 3 is thick enough such that Fabry-Pérot reflections are outside the measurement window and hence the terms within the layer do not need to be considered. If the layer was thin enough, however, the reflections could be incorporated.
The algorithm, presented in Fig. 4, utilises a frequency dependent two-dimensional scanning approach of a parameter space, involving the calculation of the transfer function $H(\omega )_{guess}$ for a range of estimates of the refractive index $n$ and extinction coefficient $\kappa$ of the sample. The sample thickness $L$ is given by the phase correction from the Kramers-Kronig algorithm. The transfer function is then discretised with respect to the frequency and the error between the measured and estimated transfer functions, $H(\omega )_{meas}$ and $H(\omega )_{guess}$, is then calculated using the combined L2 norm of the amplitude and phase. The output values of the refractive index and extinction coefficient are given by the minimum error, as a function of frequency.
While the simplicity of the approach enables straightforward implementation for applications, there exists limitations of the algorithm. The minimum detectable thickness, for example, is naturally system dependent, since it depends on the pulse duration. CST Microwave Studio was used to study the thickness limit through application of the Fabry-Pérot algorithm to the reflected field, for decreasing structure thicknesses. A 0.3 – 1 THz Gaussian pulse was used. The algorithm was able to successfully model the reflections until thicknesses reached as low as 0.1 $\mu$m. Uncertainty of the extracted values depends on the increments of the parametric sweep of dielectric properties. Additional uncertainty stems from the discretisation of the transfer function with respect to frequency, and in turn, the frequency resolution. Specifically for Gaussian beam propagation as in THz TDS systems and oblique incidence, an uncertainty in path length will be caused by the finite width of the Gaussian beam, that is not taken into account in the above derivation based on geometrical optics. All these factors affect computational run time, and hence compromises must be made depending on the application.
Figure 5 presents the reflected transfer function for normal incidence (i.e. $\theta _1 = \theta _2 = 0$, and $A_{wo,m}=1$) for a 100 $\mu$m structure, pre- and post-application of the Fabry-Pérot modelling algorithm. The temporal waveforms were simulated using the transient solver of CST Microwave Studio and the structure was designed to have dielectric properties comparable to the glass-backed HO bone (see sample configuration in the inset of Fig. 3). One can observe low frequency oscillations attributed to the Fabry-Pérot reflections within the structure. The algorithm successfully models the reflections, providing good agreement between the measured and modelled transfer function. This enables the removal of these oscillating artefacts, which propagate through to the refractive index and extinction coefficient results. The right panel of Fig. 5 presents the corrected results, demonstrating the successful extraction of the correct dielectric properties.
Fig. 5. (Left) Amplitude and phase of the transfer function $H(\omega )$ for the structure presented in Fig. 3, simulated using the transient solver of CST Microwave Studio. Layer 1 was modelled as air, while layers 2 and 3 were modelled as bone and glass, respectively. The dashed lines illustrate the measured transfer function, while the transparent lines illustrate the transfer function modelled using the Fabry-Pérot algorithm. (Right) Refractive index and extinction coefficient of the structure pre- and post-application of the Fabry-Pérot modelling algorithm, indicated by the dotted and dashed lines, respectively. The transparent lines illustrate the input dielectric properties in CST Microwave Studio.
It is important to note that the phase sensitivity of the dielectric response means that the performance of the Kramers-Kronig algorithm prior to Fabry-Pérot modelling is crucial for reliable dielectric properties retrieval. For the given sample setup, by applying the Kramers-Kronig algorithm, the plane of the reference reflection is essentially brought forward to the sample plane. The propagation $P_1$ is therefore equal to unity.
3. Results and discussion
To demonstrate the performance of the hybrid algorithm for reliable dielectric imaging, it is applied to TDS images of optically thin biological samples which exhibit an offset between the sample and reference planes. The structures are $\sim$100 $\mu$m thick HO bone slices, embedded in resin and fixed to a glass slide, as photographed in Fig. 6(a) (see Appendix A for details). The normal incidence measurements were taken in reflection geometry through single pixel raster scan imaging. The peak-to-peak value of the temporal E-field is presented in Fig. 6. The reference measurement was taken from a glass slide pixel, which sits further back from the bone (see the inset of Fig. 3), inducing a phase offset, while the thickness of the slices mean that the initial and subsequent reflections are unresolvable.
Fig. 6. (a) Photograph of a heterotopic ossification (HO) bone slice embedded in resin and fixed to a glass slide. (b) XRF image, indicating the elemental distributions of the bone slice, where calcium, phosphorus and sulphur (representing collagen) are represented by red, green and blue, respectively. The yellow regions indicate the colocalisation of calcium and phosphorus (representing CaP). Images of the (c) peak-to-peak value of the temporal reflected electric field, (d) extracted sample thickness. Normal incidence images of the extracted (e) refractive index and (f) extinction coefficient values at 0.3 THz, 0.4 THz and 0.5 THz, post application of the retrieval algorithm.
It is worth noting that for these samples in particular, no spectral fingerprints are expected in the probing frequency range. Spectral measurements of the constituent compounds of bone structures, including collagen and HA, are presented in Fig. 7 in Appendix B. The responses reveal reasonably featureless dielectric properties in this regime, for all compounds. Identification of compounds therefore relies on broadband differences in the refractive index and extinction coefficient.
Images of the dielectric properties of the structure calculated prior to the application of the hybrid algorithm are presented in Fig. 8 in Appendix C, at 0.3 THz, 0.4 THz and 0.5 THz. The images display unphysical dielectric properties, with refractive indices of less than one and negative extinction coefficients. These are attributed to both the phase offset between the bone and glass planes and the Fabry-Pérot reflections supported by the bone cavity. To demonstrate the performance of the hybrid algorithm for a single pixel, Fig. 9 (Appendix C) presents the hyperspectral dielectric properties of the structure pre- and post-application. One can observe the removal of the oscillating artifacts, attributed to phase accumulation and Fabry-Pérot reflections simultaneously, to reveal dielectric properties with low frequency dependence.
Figure 6 presents the imaged dielectric properties of the bone slice subsequent to the application of the algorithm. The composition of the bone sample is revealed in the XRF image in Fig. 6(b). The image identifies the bulk of the sample as collagen (blue), with small islands of hydroxyapatite (HA), a form of calcium phosphate (yellow). From measurements of the dielectric properties of the constituent compounds of bone structures presented in Fig. 7 in Appendix B, collagen is found to have a refractive index of $\sim$1.5, with low frequency dependence, and an extinction coefficient of $\sim$0.04 – 0.09, from 0.3 – 1 THz. The small deviation from the retrieved values of $\sim$1.8 and $\sim$0.1 for the refractive index and extinction coefficient respectively, is thought to be attributed to an increased hydration level in the bone structures, to which terahertz is highly sensitive to. Water represents up to 25$\%$ of bone [37]. Although room-temperature dehydrated bone samples lose most of the bulk water and surface water layer, strongly bound and structural water remains [38]. Indeed, even after a 225 $^\circ$C drying process, $\sim$40$\%$ of such tightly bound and structural water remains [38]. Hence, structural water is expected in the bone samples, while the collagen samples are dehydrated [39]. These findings demonstrate the success of the algorithm for dielectric extraction, with vast improvement upon the dielectric properties extracted prior to application.
It should be noted that the regions of HA are undetectable in the terahertz images due to large scattering of the field at the edges of the structure and low resolution. The resolution could be improved through excitation with higher frequencies. The dielectric images demonstrate marginal improvement upon resolution from 0.3 THz to 0.5 THz, however, pushing to higher frequencies would be beneficial. Unfortunately, the TDS system in use is limited at high frequencies due to low signal-to-noise ratio.
The dielectric images demonstrate the success of the hybrid algorithm for biometric imaging and its ability to deal with complex sample configurations. This is essential for practical biological application where controllability is limited. It should be noted that the free-standing sample set-up (i.e. not sandwiched between two windows, as commonly configured [19]), while making the system far more versatile, means that results are affected by surface roughness and resultant scattering [40]. This should be considered when dealing with structures of high roughness. The assumptions of plane wave excitation and a linear electromagnetic response of the sample should also be considered for application of the algorithm. The former assumption is important for TDS excitation where the divergence of the beam is large, particularly at low frequencies, and its profile deviates from a typical Gaussian distribution [41].
4. Conclusion
Application of this novel hybrid Kramers-Kronig and Fabry-Pérot based algorithm becomes highly valuable in the retrieval of biometric information of biological materials, where samples exhibit high terhertz absorption and sample placement and thickness are difficult to control. Its application has been demonstrated through retrieval of dielectric information of glass-backed HO human bone slices, for which bulk compound composition is identified. The causal basis of the Kramers-Kronig approach and lifted limitation of sample thickness attributed to the Fabry-Pérot approach grants great flexibility of measurements, and hence wide application to biological materials.
Appendix A: Preparation of bone samples
Heterotopic ossification (HO) samples were obtained from consented patients undergoing routine excision surgery at University Hospital Birmingham in United Kingdom. The HO samples were frozen immediately after excision surgery and stored at −80$^{\circ }$. Samples were thawed at ambient temperature for 1 h before use, and then embedded in optimal cutting temperature compound (VWR Chemicals, Leuven, Belgium) under vacuum (20 kPa) using a Struers CitoVac (Struers Ltd, Rotheram, United Kingdom). The embedded samples were equilibrated to −18$^{\circ }$C after refreezing at −80$^{\circ }$C overnight, and sectioned into 100 $\mu$m sample slices using a Leica microtome with a tungsten carbide blade. The obtained sample slices were mounted on a glass microscope slide and dried in ambient conditions for 72 h. Tissue transfer and handling was conducted under approval of the National Research Ethics Service (15/NW/0079) and in accordance with the Human Tissue Act 2004.
Appendix B: Dielectric properties of constituent compounds of bone
The dielectric properties of constituent compounds which make up bone structures were measured using transmission TDS. The compounds include amorphous calcium phosphate (ACP), calcium carbonate (CaCO$_3$), collagen and hydroxyapatite (HA). The compounds were in powdered form, and pressed into pellets.
Synthesis of hydroxyapatite (HA): A solution of (NH$_4)_2$HPO$_4$ (0.58 M) was prepared in deionised water (DIW, 125 ml). $\textrm {Ca}(\textrm {NO}_3)_2 \cdot 4\textrm {H}_2\textrm {O}$ solution (0.48 M) prepared in DIW (175 ml) was added dropwise into the prior solution under vigorous magnetic agitation at ambient temperature. Both of the solutions were adjusted to pH 11 prior to mixing, and the whole procedure proceeded in a sealed vessel. The white precipitate was harvested by centrifugation (4000 rpm, 1 min) and washed with DIW for 3 times. The resulting HA powder was obtained by freeze drying.
Synthesis of amorphous calcium phosphate (ACP): $\textrm {Ca}(\textrm {NO}_3)_2 \cdot 4\textrm {H}_2\textrm {O}$ (9.26 g) was dissolved in DIW (110 ml), containing ammonia solution (28 wt.%, 8 ml). Ammonium phosphate solution was prepared by dissolving (NH$_4)_2$HPO$_4$ (5.44 g) in DIW (260 ml) containing ammonia solution (28 wt.%, 8 ml). The prepared ammonium phosphate solution was poured into the prior calcium nitrate solution under vigorous magnetic agitation at ambient temperature. The white precipitate was immediately harvested by centrifugation (4000 rpm, 1 min) and washed with DIW for 3 times. The resulting ACP powder was obtained by freeze drying.
Synthesis of calcite: A Na$_2$CO$_3$ solution (0.1 M, 100 ml) was rapidly poured into CaCl$_2$ solutions (0.1 M, 100 ml) under intense magnetic agitation and kept stirring at 400 rpm for 1 day at ambient temperature. The white precipitate was obtained by centrifugation (4000 rpm, 1 min) and washed with DIW for 3 times. The resulting calcite powder was dried at 65 $^{\circ }$C for 5 days.
Preparation of sample pellets: The synthesised HA, ACP, calcite and commercially-purchased collagen powder (0.2 g) were compressed uniaxially in a stainless steel die with an internal diameter of 13.0 mm (Specac, UK) at a force of 7 kN by a universal material testing machine (Z030, Zwick Roell, Germany), respectively.
TDS settings: Focused illumination and detection was achieved using TPX50 lenses with 50 mm effective focal length. The dielectric properties were extracted using the commercial software Teralyzer, as in [43].
Fig. 7. The refractive index and extinction coefficient of constituent compounds of bone structures: ACP, CaCO$_3$, collagen and HA, measured using terahertz transmission TDS.
Appendix C: Terahertz dielectric images pre-algorithm application
Fig. 8. Images of the refractive index and extinction coefficient calculated from TDS reflections measurements of glass-backed HO bone slices, without application of the hybrid algorithm, at frequencies of 0.3 THz, 0.4 THz and 0.5 THz. Note the limits of the colour bar.
Fig. 9. Refractive index and extinction coefficient of a $\sim$100 $\mu$m thick glass-backed HO bone sample pre- and post-application of the Fabry-Pérot modelling algorithm.
Appendix D: TDS bone measurement method
The samples were characterized using an all fiber-coupled THz time-domain spectrometer TERA K15 from Menlo Systems with lock-in detection. The time-constant was set to 300 $\mu$s. The temporal length of the waveforms was 52 ps, providing a spectral resolution of 15 GHz. A standard reflection configuration was setup with a Si beam-splitter. The field was normally incident with accuracy of 0.3 deg and focused on the sample using a TPX35 lens with 35 mm effective focal length. To generate the images, a single pixel imaging technique was used, through translation of the sample fixed to an $xy$-translation stage, with steps of 0.5 mm.
Funding
University of Birmingham (Birmingham Fellowship); UK Research and Innovation (2137478); Royal Society (IEC/NSFC/191104, IES/R3/183131); Engineering and Physical Sciences Research Council (EP/S018395/1, EP/V001655/1).
Acknowledgments
The authors would like to thank Dr Pavel Penchev for measurements of the 3D profiles of the bone samples (not reported here) and V. Lucarini et al. for the use of the Kramers-Kronig MATLAB scripts developed in [42].
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.
References
1. G. Carpintero, L. E. García-Muñoz, H. L. Hartnagel, S. Preu, and A. V. Räisänen, Semiconductor Terahertz Technology: Devices and Systems at Room Temperature Operation (Wiley, 2015).
2. G. Tzydynzhapov, P. Gusikhin, V. Muravev, A. Dremin, Y. Nefyodov, and I. Kukushkin, “New real-time sub-terahertz security body scanner,” J. Infrared, Millimeter, Terahertz Waves 41(6), 632–641 (2020). [CrossRef]
3. D. Zimdars, J. S. White, G. Stuk, A. Chernovsky, G. Fichter, and S. Williamson, “Large area terahertz imaging and non-destructive evaluation applications,” Philos. Trans. R. Soc. Lond. A, Math. Phys. Eng. Sci. 48(9), 537–539 (2006). [CrossRef]
4. M. Karaliunas, K. Nasser, A. Urbanowicz, I. Kasalynas, D. Brazinskiene, S. Asadauskas, and G. Valusis, “Non-destructive inspection of food and technical oils by terahertz spectroscopy,” Sci. Rep. 8(1), 18025 (2018). [CrossRef]
5. J.-H. Son, Terahertz Biomedical Science and Technology (CRC Press, 2014).
6. P. Jepsen, D. G. Cooke, and M. Koch, “Terahertz spectroscopy and imaging – modern techniques and applications,” Laser Photonics Rev. 5(1), 124–166 (2011). [CrossRef]
7. C. Yu, F. S. Y. Sun, and E. Pickwell-Macpherson, “The potential of terahertz imaging for cancer diagnosis: A review of investigations to date,” Quant. Imaging Medicine Surg. 2(1), 33–45 (2012). [CrossRef]
8. S. Sy, S. Huang, Y.-X. Wáng, J. Yu, A. T. Ahuja, Y.-T. Zhang, and E. Pickwell-MacPherson, “Terahertz spectroscopy of liver cirrhosis: investigating the origin of contrast,” Phys. Medicine Biol. 55(24), 7587–7596 (2010). [CrossRef]
9. P. Knobloch, C. Schildknecht, T. Kleine-Ostmann, M. Koch, S. Hoffmann, M. Hofmann, E. Rehberg, K. Sperling, M. Donhuijsen, G. Hein, and K. Pierz, “Medical thz imaging: an investigation of histo-pathological samples,” Phys. Medicine Biol. 47(21), 3875–3884 (2002). [CrossRef]
10. M. R. Stringer, D. N. Lund, A. P. Foulds, A. Uddin, E. Berry, R. E. Miles, and A. G. Davies, “The analysis of human cortical bone by terahertz time-domain spectroscopy,” Phys. Medicine & Biol. 50(14), 3211–3219 (2005). [CrossRef]
11. M. Bessou, B. Chassagne, J. P. Caumes, C. Pradère, P. Maire, M. Tondusson, and E. Abraham, “Three-dimensional terahertz computed tomography of human bones,” Appl. Opt. 51(28), 6738–6744 (2012). [CrossRef]
12. B. A. Knyazev, V. V. Gerasimov, A. M. Gonchar, and N. G. Kolosova, “Using of terahertz radiation for monitoring of senile osteoporosis development,” IRMMW-THz2007 - Conference Digest of the Joint 32nd International Conference on Infrared and Millimetre Waves, and 15th International Conference on Terahertz Electronics pp. 563–564 (2007).
13. E. Y. Chao, N. Inoue, T. K. Koo, and Y. H. Kim, “Biomechanical considerations of fracture treatment and bone quality maintenance in elderly patients and patients with osteoporosis,” Clin. Orthop. Relat. Res. 425, 12–25 (2004). [CrossRef]
14. V. Lucarini, Y. Ino, K. E. Peiponen, and M. Kuwata-Gonokami, “Detection and correction of the misplacement error in terahertz spectroscopy by application of singly subtractive kramers-kronig relations,” Phys. Rev. B: Condens. Matter Mater. Phys. 72(12), 125107 (2005). [CrossRef]
15. E. M. Vartiainen, Y. Ino, R. Shimano, M. Kuwata-Gonokami, Y. P. Svirko, and K.-E. Peiponen, “Numerical phase correction method for terahertz time-domain reflection spectroscopy,” J. Appl. Phys. 96(8), 4171–4175 (2004). [CrossRef]
16. L. Duvillaret, F. Garet, and J. L. Coutaz, “A reliable method for extraction of material parameters in terahertz time-domain spectroscopy,” IEEE J. Select. Topics Quantum Electron. 2(3), 739–746 (1996). [CrossRef]
17. D. Liu, T. Lu, and F. Qi, “A reliable method for removing fabry-perot effect in material characterization with terahertz time-domain spectroscopy,” IEEE Trans. Terahertz Sci. Technol. 10(5), 443–452 (2020). [CrossRef]
18. I. Pupeza, R. Wilk, and M. Koch, “Highly accurate optical material parameter determination with thz time-domain spectroscopy,” Opt. Express 15(7), 4335–4350 (2007). [CrossRef]
19. M. Naftaly, Terahertz Metrology (Artech House, 2015).
20. S. Nashima, O. Morikawa, K. Takata, and M. Hangyo, “Measurement of optical properties of highly doped silicon by terahertz time domain reflection spectroscopy,” Appl. Phys. Lett. 79(24), 3923–3925 (2001). [CrossRef]
21. M. Khazan, R. Meissner, and I. Wilke, “Convertible transmission-reflection time-domain terahertz spectrometer,” Rev. Sci. Instrum. 72(8), 3427–3430 (2001). [CrossRef]
22. A. Pashkin, M. Kempa, H. Němec, F. Kadlec, and P. Kužel, “Phase-sensitive time-domain terahertz reflection spectroscopy,” Rev. Sci. Instrum. 74(11), 4711–4717 (2003). [CrossRef]
23. M. Bernier, F. Garet, J.-L. Coutaz, H. Minamide, and A. Sato, “Accurate characterization of resonant samples in the terahertz regime through a technique combining time-domain spectroscopy and kramers–kronig analysis,” IEEE Trans. Terahertz Sci. Technol. 6(3), 442–450 (2016). [CrossRef]
24. H. Son, D.-H. Choi, and G.-S. Park, “Improved thickness estimation of liquid water using kramers–kronig relations for determination of precise optical parameters in terahertz transmission spectroscopy,” Opt. Express 25(4), 4509–4518 (2017). [CrossRef]
25. K. E. Peiponen, E. M. Vartiainen, and T. Asakura, Dispersion, Complex Analysis and Optical Spectroscopy (Springer, Berlin, 1999).
26. J. S. Toll, “Causality and the dispersion relation: Logical foundations,” Phys. Rev. 104(6), 1760–1770 (1956). [CrossRef]
27. H. M. Nussenzveig, Causality and Dispersion Relations (Academic Press, New York, 1972).
28. F. Bassani and V. Lucarini, “Asymptotic behaviour and general properties of harmonic generation susceptibilities,” Eur. Phys. J. B 17(4), 567–573 (2000). [CrossRef]
29. K. E. Peiponen, V. Lucarini, J. J. Saarinen, and E. Vartiainen, “Kramers-kronig relations and sum rules in nonlinear optical spectroscopy,” Appl. Spectrosc. 58(5), 499–509 (2004). [CrossRef]
30. V. Lucarini, F. Bassani, K. E. Peiponen, J. J. Saarinen, C. Socie, and I. D. Fisica, “Dispersion theory and sum rules in linear and nonlinear optics,” La Rivista del Nuovo Cimento 26, 1 (2003). [CrossRef]
31. V. Lucarini, J. J. Saarinen, and K. E. Peiponen, “Multiply subtractive kramers–krönig relations for arbitrary-order harmonic generation susceptibilities,” Opt. Commun. 218(4-6), 409–414 (2003). [CrossRef]
32. W. Withayachumnankul, B. Ferguson, T. Rainsford, S. Mickan, and D. Abbott, “Direct fabry-pérot effect removal,” Fluct. Noise Lett. 06(02), L227–L239 (2006). [CrossRef]
33. M. Scheller, C. Jansen, and M. Koch, “Analyzing sub-100-μm samples with transmission terahertz time domain spectroscopy,” Opt. Commun. 282(7), 1304–1306 (2009). [CrossRef]
34. R. Fastampa, L. Pilozzi, and M. Missori, “Cancellation of Fabry-Perot interference effects in terahertz time-domain spectroscopy of optically thin samples,” Phys. Rev. A 95(6), 063831 (2017). [CrossRef]
35. R. Peretti, S. Mitryukovskiy, K. Froberger, M. A. Mebarki, S. Eliet, M. Vanwolleghem, and J. F. Lampin, “Thz-tds time-trace analysis for the extraction of material and metamaterial parameters,” IEEE Trans. THz Sci. Technol. 9(2), 136–149 (2019). [CrossRef]
36. P. U. Jepsen, U. Møller, and H. Merbold, “Investigation of aqueous alcohol and sugar solutions with reflection terahertz time-domain spectroscopy,” Opt. Express 15(22), 14717–14737 (2007). [CrossRef]
37. M. Granke, M. D. Does, and J. S. Nyman, “The role of water compartments in the material properties of cortical bone,” Calcif. Tissue Int. 97(3), 292–307 (2015). [CrossRef]
38. E. E. Wilson, A. Awonusi, M. D. Morris, D. H. Kohn, M. M. J. Tecklenburg, and L. W. Beck, “Three structural roles for water in bone observed by solid-state nmr,” Biophys. J. 90(10), 3722–3731 (2006). [CrossRef]
39. J. Einbinder and M. Schubert, “Binding of mucopolysaccharides and dyes by collagen,” J. Biol. Chem. 188(1), 335–341 (1951). [CrossRef]
40. S. M. Sabery, A. Bystrov, M. Navarro-Cía, P. Gardner, and M. Gashinova, “Study of low terahertz radar signal backscattering for surface identification,” Sensors 21(9), 2954 (2021). [CrossRef]
41. S. Freer, A. Gorodetsky, and M. Navarro-Cia, “Beam profiling of a commercial lens-assisted terahertz time domain spectrometer,” IEEE Trans. THz Sci. Technol. 11(1), 90–100 (2021). [CrossRef]
42. V. Lucarini, J. J. Saarinen, K.-E. Peiponen, and E. M. Vartiainen, Kramers-Kronig Relations in Optical Materials Research (Springer-Verlag Berlin Heidelberg, 2005).
43. M. Ma, Y. Wang, M. Navarro-Cía, F. Liu, F. Zhang, Z. Liu, Y. Li, S. M. Hanham, and Z. Hao, “The dielectric properties of some ceramic substrate materials at terahertz frequencies,” J. Eur. Ceram. Soc. 39(14), 4424–4428 (2019). [CrossRef]
