Super-resolution time-domain spectroscopy method and apparatus for sample characterization

Abstract

A method for determining a set of physical parameters of a sample, comprising the following steps: - A acquisition of a measured sample time trace Es(t), - B acquisition of a measured reference time trace Eref(t), - C determination of a broadened reference time trace, called Eref0(t), and determination of the discrete Fourier transform of the broadened reference time trace - D determination of a modeling of the impulse response of the sample in the frequency domain, called sample frequency model, from the Fourier transform of the broadened reference time trace and a physical behavior model of the sample, according to a set of physical parameters (pi), - E application of an optimization algorithm to the set of physical parameters (pi), comprising the following sub-steps: - El initialization of the physical parameters (pi), - iteratively implementing the following sub-steps: - E2 computation of the inverse discrete Fourier transform of the sample frequency model, called estimated sample time trace E est {p i}(t), - E3 computation of an error function (ε er {pi}) until obtaining a set of values of the physical parameters (pi opt ) that minimize the error function.

Description

TECHNICAL FIELD

[0001] The present invention relates generally to a method for determining physical parameters of a sample by time-domain spectroscopy (TDS). In particular, the present invention relates to a method for fitting a measured time trace of a sample and extracting physical parameters therefrom. BACKGROUND

[0002] Time-domain spectroscopy methods are today a well-established technique in the terahertz time-domain (THz-TDS) regime, i.e. by using an excitation beam with frequencies between 0.1 THz and 30 THz, have shown the ability to investigate different materials such as semiconductors, ferroelectrics, superconductors, liquids, gases, biomolecules, molecular crystals such as carbohydrates, bandpass filters, metasurfaces embedded in microfluidic circuits, etc.

[0003] THz-TDS is based on a direct measurement of the electric field in the THz frequency range, in contrast to Fourier-transform infrared spectrometers where the measured time-domain data are autocorrelated.

[0004] The working principle of a common THz-TDS setup is depicted in Figure 1 The THz excitation beam EB is emitted by means of the optical rectification effect of near-infrared pulses generated by a femtosecond laser L by means of a photoconductive THz antenna or a nonlinear crystal NL as shown in Figure 1 The excitation beam EB is shown in Figure 2 It consists of periodically emitted electromagnetic / optical pulses 20 (left side) with a period T having a frequency spectrum 21 (right side) comprising a frequency comb. Both periodicities are important for the TDS method to work.

[0005] Today electromagnetic (THz) pulses can be obtained by optical rectification of femtosecond laser pulses or by illuminating a photoconductive THz antenna with femtosecond laser pulses.

[0006] In order to be in the THz range, the pulses should last from a few hundred femtoseconds to a few picoseconds. The repetition rate frep of the laser, equal to 1 / T, can vary from GHz to sub-Hz repetition rates.

[0007] Next, a lens or a parabolic mirror is used to collimate the pulse beam and direct it to the sample S under investigation. The transmitted (or reflected) pulses originating from the sample are then collected by an optical system and aligned onto a detector D. The detector measures the electric field of the electromagnetic wave as a function of time on a scale ranging from femtoseconds to a few hundred picoseconds or even a few nanoseconds. This can be done by means of photoconductive or electro-optic sampling. The time sampling is usually performed by means of a delay line DL made of mirrors mounted on motorized translation stages, by coherent detection as shown in Figure 3As shown, a delay At with a maximum time offset tmax is introduced. The time sampling can also be performed via the beat of two slightly different repetition rate frequency combs, a technique known as HASSOPS. The time interval between two measurements is called the sampling period ts, with a sampling frequency fs = 1 / t. In the following, we call tmax the time range over which the measurement of the electric field is achieved, called the time trace.

[0008] A typical sampling period ts is between 10 fs and 50 fs. The ability to directly measure the electric field of a THz pulse instead of the average energy yields both the phase and the amplitude of the waveform, thus providing information about the absorption coefficient and the refractive index of the sample. We call the measured sample time trace Es(t) and the measured reference time trace Eref(t) without the sample.

[0009] For material analysis, the usual method to retrieve the material parameters is to perform a Fourier transform of the recorded (with and without the sample) pulse time traces Es(t) and Eref(t). The ratio between these two spectra in a linear approximation is called the complex transmission coefficient and can be written as:

[0010]

[0011] and are the Fourier transforms of the time domain signals E s (t) and E ref (t), respectively, is the complex refractive index, where the real part corresponds to the delay and the imaginary part corresponds to the absorption in the material, d is the thickness of the sample that has to be measured, and co is the angular frequency, which is related to the frequency by the formula co = 2pf. The term is the product of the Fresnel coefficients for normal incidence of the two air / material interfaces, and is a term that takes into account the Fabry-Pérot multiple reflections in the sample.

[0012] Equation (1) sets the so-called “forward problem”: knowing and one can obtain Because the experiment gives and the actual interest is the “inverse problem”, i.e. using the information of and one can determine

[0013] In general, a first method 10 for determining the physical parameters of a sample of interest has the following principle steps:

[0014] - performing experiments with and without sample, that is, measuring Es(t) and Eref(t) (in Figure 4 An example of lactose is given in the middle, which has a region 5) corresponding to an oscillation following a damped sinusoid,

[0015] - calculating the discrete Fourier transform of the two time traces and over a time window equal to tmax, where there are measurements (see Figure 5 which is based on the signal of Figure 4 the image 6 corresponds to the principal component of the Fourier transform of the damped sinusoid 5), and determining the transfer function of the measurement

[0016] - measuring the thickness of the sample;

[0017] - extracting the real and / or imaginary part of the refractive index by minimizing an error function (using for example a gradient-free or quasi-Newton minimization algorithm) The error function is defined based on the modulus error and the phase error between the modeled transmission coefficient and the measured transmission coefficient Figure 6 (see For lactose, d = 900 pm ± 20 pm;

[0018] - fitting the refractive index to obtain material parameters, where the refractive index is assumed to follow a model such as the Drude-Lorentz model. The material parameters are typically the strength and resonance frequency and / or the line width of the spectral line (see below).

[0019] Regarding step ii), the discrete Fourier transform (DFT) is a mathematical tool for determining the Fourier transform of a periodic discrete function. Most commonly, the calculation of the DFT is performed by means of the FFT algorithm (fast Fourier transform).

[0020] The electromagnetic pulse that interacts with the sample undergoes a delay related to the real part of the refractive index of the sample and an absorption related to the imaginary part of the refractive index. Figure 4 The delay R of the oscillation 5 is a time delay related to The oscillation 5 shows the wave points in the view of the absorption of the sample and the peak 6 is the absorption peak in the frequency domain related to Figure 4 The oscillation 5 shows the wave points in the view of the absorption of the sample and the peak 6 is the absorption peak in the frequency domain related to The oscillation 5 shows the wave points in the view of the absorption of the sample and the peak 6 is the absorption peak in the frequency domain related to

[0021] Regarding step (v), the Drude Lorentz model is based on the assumption that, in the sample, light interacts with electrons in the form of damped harmonic oscillators. Its equation of motion can be written in the time domain by an exponential decayed sinusoid (see Figure 4 in 5) which corresponds to a Lorentz curve in the frequency domain. The absorption line thus follows a Lorentz curve as illustrated in Figure 7 . This model is used exclusively in gas spectroscopy, provided the pressure is greater than a few tens of microbar.

[0022] Each line following a Lorentz curve is characterized by 3 parameters: the maximum M, the line width g, the resonance frequency f0, corresponding respectively to the physical parameters pi, p2 and p3 (in this case i = 1 to 3).

[0023] This first method has shown promising results, but has several drawbacks:

[0024] First, it requires an accurate measurement of the sample thickness. Second, the result does not take into account the causality (in this case, it takes the form of the Kramers-Kroenig relation), which can lead to erroneous results.

[0025] To overcome some of the drawbacks of the above method, the publication “THz-TDS time-trace analysis for the extraction of material and metamaterial parameters”, by Peretti et al., IEEE transactions on Terahertz science and technology, vol. 9, no. 2, March 2019, referred to as the Peretti publication, describes a second method based on a comparison of the initial time-domain data of the measured THz pulse. The theoretical approach for implementing this method would be:

[0026] - (i) acquiring the measured sample time trace Es(t) and the measured reference time trace Eref(t) over a time period tmax,

[0027] - (ii) determining a model E i {p model}(t) from a set of parameters {p i} that depict how the sample transforms the reference pulse into a modeled pulse,

[0028] - (iii) determining an objective function Obj{p iwhich is equal to the L 2 norm of the difference between the modeled pulse and the measured (sampled) pulse

[0029]

[0030] (iv) determining the set of values of the parameters p, that minimizes the objective function.

[0031] This method does not require information on the thickness d of the sample and takes into account the causality, since the thickness appears as a parameter and because the error calculation is done in the time domain.

[0032] Unfortunately, the modeling that describes how the sample transforms the reference pulse into the modeled pulse is unknown in the time domain. All known models describe the sample behavior in the frequency domain:

[0033] But the real advantage of the representation of formula (2) is given by the Parseval theorem, which states that the norm of a function is the same as the norm of its Fourier transform, which represents the objective function Obj{p i} can be calculated by using the following formula:

[0034]

[0035] where:

[0036] ω max = π / dt, where dt corresponds to the sampling period ts

[0037] ω min = -π / dt

[0038] Thus, in formula (3), the objective function is represented as a function of the difference between the modeling of the sample behavior in the frequency domain and the discrete Fourier transform of the measured sample time trace Es(t).

[0039] This is very convenient, which allows to calculate the objective function in both the time and the frequency domain. This method can thus be implemented directly in the frequency domain, avoiding to calculate the Fourier transform at each iteration.

[0040] By using the DFT, the two above methods assume that their data is periodic, with a period that corresponds to the length of time for which the recording is made. Thus, they assume that the time distance T between the two pulses is exactly equal to the maximum shift tmax of the delay line, which is usually not true.

[0041] The second method 20 described in Peretti's publication is shown in Figure 8 , which is implemented in the frequency domain but corresponds to a fit in the time domain, and comprises the following steps:

[0042] - Step A0: acquisition of the measured sample time trace Es(t);

[0043] - Step B0: acquisition of the measured reference time trace Eref(t),

[0044] - Step C0: determination of the discrete Fourier transform of Es(t) on a time window equal to tmax

[0045] - Step C0': determination of the discrete Fourier transform of Eref(t) on a time window equal to tmax

[0046] These steps are identical to steps (i) and (ii) of the first method 10.

[0047] - Step D0: determination of a sample frequency model based on a physical behavior model of the sample according to a set of (pre-determined) physical parameters pi This model performs a modeling of the impulse response of the sample in the frequency domain.

[0048] Typically, the model is the Fourier transform of the reference time trace is multiplied by a transfer function related to the complex refractive index according to formula (1) .

[0049]

[0050] - Step E0: application of an optimization algorithm to the set of physical parameters pi, comprising the following steps:

[0051] - E1 initialization of the physical parameters

[0052] - then iteratively implementing a step E2 of determination of the difference between the discrete Fourier transform of the measured sample time trace and the sample frequency model determines an objective function Obj{pi} until the values of the parameters pi that minimize the objective function are obtained. For example, the iterations are performed until the objective function becomes less than a predetermined threshold.

[0053] For example, the objective function is determined based on the square of the difference according to the following (least squares algorithm)

[0054]

[0055] This method performs the calculation of several steps, implementing a fitting on the time domain, but using simple calculations on the frequency domain.

[0056] For bulk solid materials, the publication by Peretti implements a multi-Drude Lorentz model which defines the dielectric permittivity of the sample as a collection of electronic harmonic oscillators (matrix vibrations, oscillating charges, etc.) and yields the following permittivity function:

[0057]

[0058] where ε ∞ is the dielectric permittivity at high frequencies relative to the range of interest, ω p is the plasma frequency, γ p is the damping rate, k max is the number of oscillators considered, ω 0k , γ k and Δε k are respectively the resonance angular frequency, the damping rate (line width) and the strength (maximum value expressed in units of permittivity) of the kth oscillator.

[0059] The term is the Drude term due to the presence of free carriers.

[0060] Substituting equation (6) into equation (1) of the model is defined by equation (5).

[0061] Thus, in the case of a single homogeneous layer without free carriers (Drude term equal to 0), the propagation will be modeled by the Fresnel coefficients and a multi-Lorentz harmonic oscillator model using a set of 3k+1 physical parameters pi.

[0062] The optimization algorithm will thus provide a set of values for the 3 parameters p1k opt , p2k opt , p3k opt of each spectral line k plus the value of ε ∞ which minimizes the objective function.

[0063] Another model for super-surfaces and called Time-Domain Coupled Mode Theory (TDCMT) is also described in the publication by Peretti.

[0064] The above method proposes to directly model the time trace of the wave using the above oscillator model and to adjust the parameters of this oscillator (natural frequency, decay and duration of coupling with the electromagnetic wave) so that the model matches the sample time trace. To do this, an optimization algorithm is implemented in the frequency domain due to the Parceval theorem.

[0065] The method described in Peretti's publication has several advantages: (i) No need for an accurate measurement of the thickness of the sample. Indeed, obtaining an accurate sample thickness of a material such as a carbohydrate or a semiconductor wafer at sub-micron precision is challenging. Therefore, avoiding this step is a real improvement. (ii) The problem of modeling the refractive index is analyzed globally, therefore, the fit requires a small number of parameters compared to the usual two values per frequency, the method is much less sensitive to noise. This enables a very high precision on the refractive index to be achieved. (iii) Because the residual fit error is in amplitude units, one can clearly interpret this error, so a better understanding of the achieved model and possible oversights. (iv) Because the fit is performed in the time domain (the optimization is performed in the frequency domain), the phase is not lost in the presence of strong absorption and no additional step is needed (v) Finally, it allows to obtain material parameters, also including metamaterials' parameters, accurately, reliably and consistently using the Drude-Lorentz model.

[0066] Like all classical methods, the method described in Peretti's publication also has some drawbacks. The main drawback is the resolution, which is limited by the value of the maximum time delay tmax.

[0067] Indeed, the optimization is performed in the frequency domain thanks to the Parceval theorem. (By the FFT algorithm) a discrete Fourier transform is performed on the time window tmax and The calculation of the has a direct impact on the spectral resolution, which is the ability to represent the width of a measured line or to distinguish two close lines (see below).

[0068] The spectral resolution δf, that is to say the ability to measure the width of a line or to distinguish two close lines, corresponding to the frequency of the inverse of the total delay time of the delay line, is limited by the so-called Fourier uncertainty:

[0069] δf = 1 / tmax

[0070] With tmax of the order of 1 ns (best value achievable), the spectral resolution is in the range of 1 GHz, which is optimal.

[0071] This spectral limitation in the form of 1 / tmax directly results from the DFT calculation.

[0072] The standard Fourier transform of a function f is:

[0073]

[0074] For a function f constant on the interval [n, n+1[, the discrete Fourier transform is

[0075]

[0076] For a periodic function f, N discrete values y1 to yN on the period are chosen and the DFT is given by the following equation:

[0077]

[0078] Because time and frequency are two conjugated magnitudes, the spectral resolution of the DFT is limited by the Heisenberg Fourier criterion. The DFT has a resolution δf limited by the Fourier transform, which is:

[0079]

[0080] In fact, because the DFT computation has a periodicity exactly equal to tmax, it is considered that the information corresponding to the instant just after tmax is equal to the information corresponding to the instant just after 0. This is the well-known spectral aliasing. For example, the publication Xu, J., Yuan, T., Mickan, S. and Zhang, X. C. (2003), "Limit of spectral resolution in terahertz time-domain spectroscopy", Chinese Physics Letters, 20(8), 1266, explains that the frequency limitation of the resolution of equation (7) is due to the fact that the detection of the electromagnetic pulse is limited in time by tmax, assuming that the SNR (signal-to-noise ratio) is greater than 1.

[0081] There is therefore a need for a new method for determining the physical parameters of a sample based on the results of the measurements of the electromagnetic field originating from the sample on which a TDS technique is performed, which overcomes the spectral limitations explained above and which therefore has a higher spectral resolution. SUMMARY

[0082] According to a first aspect, there is provided a method for determining a set of physical parameters of a sample, comprising the steps of:

[0083] - acquiring a measured sample time trace Es(t),

[0084] The measured sample time trace Es(t) has been obtained by time-domain spectroscopy by irradiating the sample with an excitation beam periodically emitting electromagnetic pulses having a period T and exhibiting a comb-like frequency, and by detecting the electromagnetic field from the sample as a function of time by coherent detection, the time period during which the sample time trace is measured being tmax, with tmax < T,

[0085] - B acquiring a measured reference temporal trace Eref(t) obtained by irradiation and detection under the same conditions as in step A but in the absence of the sample,

[0086] - C determining a broadened reference temporal trace, called Eref0(t), which extends over the period T and is obtained by assuming a value of 0 for the instants for which no measurement is performed, and determining the discrete Fourier transform of the broadened reference temporal trace calculated over a time window equal to T

[0087] - D determining, from the Fourier transform of the broadened reference temporal trace and a physical behavior model of the sample, a modeling of the impulse response of the sample in the frequency domain, called sample frequency model ,

[0088] - E applying an optimization algorithm to a set of physical parameters, comprising the following sub-steps:

[0089] - El initializing the physical parameters,

[0090] - iteratively implementing the following sub-steps:

[0091] - E2 computing the inverse discrete Fourier transform of the sample frequency model called estimated sample temporal trace E est {p i}(t),

[0092] - E3 computing an error function from the difference between the measured sample temporal trace Es(t) and the estimated temporal trace Eest(t) until a set of values of the physical parameters is obtained that minimizes the error function.

[0093] According to a development of the first aspect, the excitation beam is in the THz domain having frequencies comprised between 100 GHz and 30 THz.

[0094] According to a further development of the first aspect, the maximum time delay tmax is chosen to include more than 95% of the energy of the measured reference temporal trace.

[0095] According to a further development of the first aspect, the sample frequency model consists in multiplying the Fourier transform Eref0(ω) by a transfer function T(ω) characterizing the behavior of the sample.

[0096] Preferably, the transfer function T(ω) depends on a complex refractive index n(ω).

[0097] According to a further development of the first aspect, for each spectral line, the square of the complex refractive index, called permittivity ε(ω), follows a Drude-Lorentz model, the spectral line being characterized by a set of three parameters called amplitude, width called damping rate, and central frequency.

[0098] According to a further development of the first aspect, the error function is defined as:

[0099]

[0100] According to a second aspect, there is provided a characterization device for characterizing a sample, the device comprising:

[0101] - a memory storing a measured sample time trace Es(t) and a measured reference time trace Eref(t),

[0102] The measured sample time trace Es(t) has been obtained by time-domain spectroscopy by illuminating a sample by means of an excitation beam periodically emitting electromagnetic pulses having a period T and presenting a comb of frequencies, and by detecting the electromagnetic field coming from the sample as a function of time by coherent detection, the time period during which the sample time trace is measured being tmax, with tmax < T,

[0103] The measured reference time trace Eref(t) has been obtained by performing the illumination and the detection in the same conditions as for the measured sample time trace Es(t) but in the absence of the sample,

[0104] - a processing unit configured to:

[0105] - determine a broadened reference time trace, called Eref0(t), extending over the period T and obtained by assuming as 0 values the instants for which no measurement has been performed, and determine the discrete Fourier transform of the broadened reference time trace computed over a time window equal to T

[0106] - determine, from the Fourier transform of the broadened reference time trace Eref0(ω) and a physical behavior model of the sample, a modeling of the impulse response of the sample in the frequency domain, called sample frequency model - apply an optimization algorithm to the physical parameters, comprising the following steps:

[0107] - initialize the physical parameters

[0108] - iteratively implement the following sub-steps:

[0109] - compute the error function

[0110] * computing said sample frequency model of the estimated sample time trace Eest(t) est {p i} the inverse discrete Fourier transform of

[0111] * computing an error function from the difference between said measured sample time trace Es(t) and the estimated time trace Eest(t) until a set of values of the physical parameters is obtained that minimizes said error function.

[0112] According to a third aspect, there is provided a spectrometer comprising:

[0113] - a characterization device according to the second aspect of the application,

[0114] - a measurement device comprising:

[0115] - a source configured to illuminate a sample by said excitation beam,

[0116] - a detector configured to detect said measured sample time trace Es(t) and said measured reference time trace Eref(t). BRIEF DESCRIPTION OF DRAWINGS

[0117] Embodiments of the application and further advantageous objectives thereof will be described in detail below with reference to the drawings, in which:

[0118] [ Figure 1 ] Figure 1 shows the principle of THz time-domain spectroscopy.

[0119] [ Figure 2 ] Figure 2 shows the time (left) and frequency (right) behavior of the excitation beam illuminating the sample.

[0120] [ Figure 3 ] Figure 3 shows the time sampling of the electromagnetic pulse.

[0121] [ Figure 4 ] Figure 4 shows an example of the measured sample time trace Es(t) and the reference time trace Eref(t) as a function of time.

[0122] [ Figure 5 ] Figure 5 shows the discrete Fourier transform of both the time trace Es(t) and the time trace Eref(t) over a time window equal to tmax and

[0123] [ Figure 6 ] Figure 6 shows the real part of the refractive index that minimizes the error function and imaginary part

[0124] [ Figure 7 ] Figure 7 Lorentz curve is shown.

[0125] [ Figure 8 ] Figure 8 Method described in the prior art is shown.

[0126] [ Figure 9 ] Figure 9 Reference signal and sample signal in time domain (left) and in spectral domain (right) obtained by FFT of the time signal are shown. Figure 9 a (bottom) shows the time signal with scale tmax and the frequency signal obtained with a time window equal to tmax, and Figure 9 b (top) shows the time signal with scale different T and the frequency signal obtained with a time window equal to T.

[0127] [ Figure 10 ] Figure 10 Method according to the application is shown.

[0128] [ Figure 11A ] Figure 11A Transmission spectrum of the sample at a pressure of 1 kPa is shown.

[0129] [ Figure 11B ] Figure 11B Zoom of the absorption spectrum around f0= 0.573 THz is shown.

[0130] [ Figure 12 ] Figure 12 Comparison of the central frequency values as a function of the pressure P obtained by different methods (method according to the application, prior art, database) is shown.

[0131] [ Figure 13 ] Figure 13 Damping rate γ values (also called FWHM) of as a function of the pressure P obtained by different methods (method according to the application, prior art, database) are shown.

[0132] [ Figure 14 ] Figure 14 Characterization device 50 for characterizing a sample according to the application is shown.

[0133] [ Figure 15 ] Figure 15 Spectrometer according to the application is shown.

[0134] [ Figure 16 ] Figure 16An embodiment of the application is shown in which the measuring device is connected to a computer I / O interface and the characterizing device is located in the computer.

[0135] [ Figure 17 ] Figure 17 Another embodiment is shown in which the measuring device is connected to a computer via an I / O interface and the characterizing device is located in a remote server connected to the computer Comp via an Internet communication subsystem. DETAILED DESCRIPTION

[0136] The inventors have carried out a detailed analysis of TDS experiments, initially finding that in the prior art method the DFT is always performed on a time window tmax. As mentioned above, the consequence is that the period of the light pulses of the excitation beam EB is assumed to be equal to tmax.

[0137] This assumption is false in real experimental conditions. The period T between two pulses is usually greater than tmax. For currently available lasers, T is in the range of s to ns (corresponding to a repetition rate of the order of GHz to Hz), that is, tmax is usually lower than 20% of T. The repetition rate of commercial lasers cannot be easily modified, and lasers with a very high repetition rate, that is, with T smaller than 10 ns, are very expensive and difficult to handle.

[0138] Therefore, in a certain sense, the calculation of the DFT performed in the prior art is counter-intuitive. The main consequence is that the spectral resolution is equal to 1 / tmax.

[0139] As Figure 9 shown, if the measurement is performed on the whole period T, an improved resolution can be obtained. The left side shows the signal in the time domain and the right side shows the curve in the spectral domain obtained by FFT of the time signal.

[0140] The signal Tr is an example of the simulated reference time trace Eref(t) and the signals T1 to T3 are examples of simulated sample time traces Es(t) obtained with samples presenting an absorption line with a frequency width v decreasing from T1 to T3 (that is, the trace in the time domain is spread from T1 to T3).

[0141] In the lower picture a), only the signal on a time tmax of 10 ps can be used, while in the upper drawing b) the signal on a time period (300 ps) greater than the period T (equal to 100 ps) can be used.

[0142] The spectrum of the lower a) is calculated by FFT using a time window equal to tmax. It can be seen that the resolution Res in the frequency domain is limited and the decrease of v cannot be observed: the spectral resolution is limited to 1 / tmax = 100 GHz.

[0143] The spectrum of the upper part b) is calculated by using an FFT with a time window equal to T. It can be seen that the resolution in the frequency domain allows to observe a decrease of v: the spectral resolution is 1 / T = 10 GHz.

[0144] The lower part of the figure shows the limitation of a real experiment performing a measurement on a time tmax smaller than T. In particular, in the experiment, the resolution is decreased according to a coefficient equal to the ratio of T / tmax (here 10). This spectral limitation is a huge drawback to characterize narrow lines of a gas.

[0145] In Figure 10 A method 30 for determining a set of physical parameters of a sample according to the application is shown in Fig. 1. The method 30 comprises a step A of acquisition of a measured sample time trace Es(t) and a step B of acquisition of a measured reference time trace Eref(t). Preferably, the two measured time traces, the sample time trace and the reference time trace (which is obtained by illuminating and detecting in the same conditions as the step A but in the absence of the sample) are obtained by a classical TDS device as previously described. The measured sample time trace Es(t) has been obtained by illuminating the sample S to be characterized by means of an excitation beam EB of electromagnetic pulses with a period T and presenting a frequency comb (which can be called optical pulses when the wavelength is less than 10 pm), and by detecting the electromagnetic field from the sample (transmitted or reflected) as a function of time by means of a coherent detection. The coherent detection generally comprises a delay line or is based on the Hassops experiment (double frequency comb).

[0146] Preferably, the excitation beam EB is in the THz domain with frequencies comprised between 0.1 THz and 30 THz, but the claimed method 30 can be extended to higher frequencies (IR, visible) when in the future illumination and detection devices will be available.

[0147] The time period during which the two time traces are measured is called tmax, with tmax smaller than T.

[0148] The steps A and B of the claimed method 30 are identical to the steps AO and BO of the method 20 of the Peretti’s publication previously described.

[0149] In step C, a broadened (also called padded) reference time trace, called Erefo(t), is determined. Erefo(t) is prolonged over the period T and is obtained by assuming a value of 0 for the instants where no measurement is performed, that is to say the instants after t = tmax and until t = T. Any arbitrary value could be assumed at these instants, but zero is the optimal value for indicating that almost no or no pulse energy can be found in the interval ]tmax, T]. For this reason, the maximum time delay tmax is preferably chosen to be large enough to include more than 95% (preferably 99%) of the energy of the measured reference time trace (the precision of the method increases according to the percentage value). This method of perfecting the unknown values by 0 before calculating the data is called "0-padding".

[0150] In step C, the discrete Fourier transform of the broadened reference time trace Erefo(t) is determined The calculation is performed on a time window equal to T, which is possible because Erefo(t) spans the entire period T.

[0151]

[0152] Thus, the DFT of is calculated, which is different from the DFT of of step C0' of the method 20, which is performed on the time window tmax.

[0153] Then, in step D, a modeling of the impulse response of the sample in the frequency domain is determined (called sample frequency model ). The sample frequency model depends on a set of physical parameters pi (i being the index of the parameters) and is determined based on the Fourier transform of the broadened reference time trace and on a model of the physical behavior of the sample. The modeling of step D is the same as the modeling of step D0 of the method 20, which differs in that the modeling of the claimed method 30 is defined based on which is different from of step D0 of the method 20.

[0154] Preferably, the sample frequency model is present in the Fourier transform multiplied by a transfer function T(co) characterizing the behavior of the sample:

[0155]

[0156] Preferably, the transfer function T(co) depends on the complex refractive index n(co), as previously described (see equation (1)).

[0157] The choice of the model depends on the sample, but is constrained because the line shape results from physical laws.

[0158] According to one embodiment, for each spectral line, the square of the complex refractive index, called permittivity ε(ω), follows (but is not limited to) the Drude-Lorentz model, that is to say, each spectral line follows a Lorentzian distribution shape. In this case, as mentioned in equation (6), the spectral line k is characterized by a set of three parameters: the amplitude p1= ε κ (the permittivity unit), the width p2= γ κ , called damping rate, and the center (resonance) frequency p3= ω 0k . This model is adapted to model the spectral lines of a gas.

[0159] In step E, an optimization algorithm is applied to the set of physical parameters pi. First, in sub-step El, the physical parameters pi are initialized, then the following sub-steps E2 and E3 are implemented iteratively:

[0160] In E2, the inverse discrete Fourier transform of the sample frequency model is computed. This inverse discrete Fourier transform is called the estimated sample time trace E est {p i}(t) and is determined according to:

[0161]

[0162] ts is the sampling time.

[0163] The real period of the experiment, that is to say the repetition rate T of the laser source, is used here for the computation. Thus, the model and the experiment follow the same periodicity, preventing any DFT aliasing artifact and reproducing the real. Therefore, because of the additional information introduced by the modeling, a narrower line than the Fourier-Herzberg limit will produce a time signal that reproduces over the right edge and at the beginning of the left edge of the time window. This signal can still be fitted, thus producing a better resolution than the Fourier-Herzberg criterion.

[0164] Then, in E3, the error function ε er {pi} is computed from the difference between the measured sample time trace Es(t) and the estimated time trace Eest(t), that is to say, ε er {pi} depends on [Es(t)-Eest(t)] according to a function f, where f is a function defining a topological distance:

[0165]

[0166] The iteration continues until a set of values of the parameters pi is obtained that minimizes the error function. The iteration stops, for example, when the error function becomes smaller than a predetermined threshold or stops to evolve.

[0167] The result of the final optimization is a set of values of the parameters pi.

[0168] For a gas spectrum following the Lorentz model, a set of 3 parameters values (Δεk γk, ω0k) is provided for each line, p3 = ω0k allows the gas identification, pi = Δεk is related to the gas pressure and p2 = γk is related to the gas temperature.

[0169] It is important to note that the optimization is performed for values of t between 0 and tmax. The recorded data Es(t) can only be obtained on tmax (as explained previously, a measurement over the whole time frame T is not possible). As it is not possible to add information that is not present on Es(t), only a part of the time frame of the model can be compared to the recorded data Es(t). Typically, about 10% of the available calculated data of Eest(t) (ratio tmax / T) is used for the error calculation.

[0170] Therefore, as the hypothesis is not fulfilled, the Parseval theorem cannot be used to perform the comparison of the curves and the error calculation in the frequency domain as it was done previously in the method 20 of the Peretti publication.

[0171] As the Parseval theorem cannot be used, a Fast Fourier Transform must be performed at each iteration in step E2 to get the modeled data in the time domain, then the error is calculated (step E3).

[0172] The error function selection depends on the type of optimization. An example of error function is the root mean square difference:

[0173]

[0174] As an example, the Augmented Lagrangian Particle Swarm Optimizer is used to perform the constrained optimization.

[0175] If applicable, the sample behavior, that is to say The formula followed by T(ω) is known. This "a priori" information is important as it allows to implement this constrained reconstruction algorithm without losing any information, so resulting in a much higher resolution (super-resolution) than the method 20 of the Peretti publication.

[0176] For example, the knowledge of the lines of a gas given by theoretical physics allows to indicate that the recorded line follows a Lorentz distribution shape in the frequency domain, corresponding to an exponentially decaying sinusoid in the time domain. Super-resolution relies on the fact that each decaying sinusoid is described by only three parameters: the amplitude, the central frequency and the damping rate. One only needs a few time domain points, combined with a constraint on the number of oscillators which is less than a fraction of the total number of points in the time trace (generalized Fourier-Heisenberg uncertainty), to get the parameters. In other words, a TDS system records the whole time trace as information, and it would like to get 3k+1 parameters from it. Performing the calculation in the time domain means that instead of searching for information of each δf (δf is the Fourier transform in frequency step) in the frequency domain, one searches for much less information, also called sparse information. Using this constraint, the resolution will be limited only by the signal to noise ratio and the calculation method used to optimize the fit, similarly to super-resolution microscopy.

[0177] The claimed method 30 of determining sample parameters can use the amount of information present in the TDS spectrum completely by performing high resolution optical spectroscopy over a very wide spectrum. For example, this allows to monitor events by looking at several hundreds of gas lines (of the same gas or mixture) and thus by measuring the relative concentrations of the different components.

[0178] To improve the speed of the claimed method which consumes many computer times (one DFT per iteration), in one embodiment of the claimed method, a harmonic inversion method is implemented instead of a constrained algorithm. The harmonic inversion can be used if all the following requirements are fulfilled: it is particularly important to ensure that the excitation pulse is very close to a Dirac pulse, so that the system is in full relaxation during a significant fraction of the time trace.

[0179] To test the performance of the method, a commercial THz-TDS TERASPEC®TERASPEC®TERASMART The ammonia (NH3) time trace was recorded. It uses a femtosecond laser with a repetition frequency of 100 MHZ (90 fs pulses). The THz pulse time is about 400 fs for a spectrum extending from 200 GHz to 5 THz. A Brewster angle silicon window gas cell with an optical path length of 8 cm was set in the optical path to make the measurements. The experiment of 13 minutes long was repeated for different pressures from 3 mb to 100 mb. Figure 11A The transmission spectrum of the sample at 1 kPa (10 mb) is given. The main absorption peak is the peak at the frequency f0= 0.573 THz.

[0180] Figure 11BAn enlarged view of the FFT showing the peak around f0= 0.573 THz. It can be seen that, using a sampling frequency fs equal to 1.2 GHz, the peak is seen only at one point.

[0181] Figure 12 The central frequency value f0as a function of the pressure P (curve 121, triangles) obtained by the claimed method 30 is shown, and Figure 13 The damping rate γ value (also referred to as FWHM) as a function of the pressure P (curve 131, triangles) obtained by the claimed method 30 is shown.

[0182] Also plotted on Figure 12 and Figure 13 are the data for NH3 extracted from the database Hitran 2016 (curves 122 and 132 respectively, squares) for comparison. The values of the physical parameters f0and γ determined by the method 30 fit the data from Hitran for both γ and f0very well, thus proving the very good accuracy of the claimed method.

[0183] Also plotted for comparison are the values of f0and damping rate obtained using the Peretti method 20 (curves 123 and 133 respectively, solid circles). It can be seen on Figure 13 that the method 20 gives good results for high FWHM values corresponding to high pressures, but the FWHM values become very inaccurate for pressures lower than 10 mb due to the insufficiently high frequency resolution of the method 20.

[0184] In addition to focusing on the line of NH3 around 530 GHz, analysis was also performed on other spectral lines, and super-resolution was achieved along the entire spectrum from 530 GHz up to more than 3 THz. More specifically, the line around 1250 GHz was also focused on. This line is in fact a doublet line spaced by 350 MHz (one third of the delay line "Fourier transform" resolution limit). Thanks to the claimed super-resolution method, the frequency of each line of this doublet was acquired with better than 30 MHz accuracy down to pressures of a few mBar. Thus, the doublet separation at the resolution limit has been achieved, which is the strongest evidence of super-resolution.

[0185] Thanks to the high resolution of the claimed method, the TDS super-resolution method 30 applied to THz pulses extends the use of THz-TDS for gas spectroscopy used for atmospheric or health purposes, such as in breath analyzers, as well as for controlling industrial environments.

[0186] Another aspect according to the application relates to Figure 14A characterization device 50 is shown for characterizing a sample S. The device 50 comprises: a memory MEM storing the measured sample time trace Es(t) and the measured reference time trace Eref(t) as previously defined; and a processing unit PU configured to:

[0187] - determining a broadened reference time trace (called Eref0(t)) extending over the period T and obtained by assuming as 0 values the instants for which no measurement was performed, and determining a discrete Fourier transform of the broadened reference time trace computed over a time window equal to T

[0188] - determining from the Fourier transform of the broadened reference time trace and from a set of physical parameters pi(ω) a modeling of the impulse response of the sample in the frequency domain (called sample frequency model

[0189] applying to the set of physical parameters (pi) an optimization algorithm comprising the following sub-steps:

[0190] - initializing the physical parameters (pi),

[0191] - iteratively implementing the following sub-steps:

[0192] - computing the inverse discrete Fourier transform of the sample frequency model called estimated sample time trace Eest(t), est {p i}(t) from the difference between the measured sample time trace Es(t) and the estimated time trace Eest(t),

[0193] - computing an error function (ε er {pi}) from the difference between the measured sample time trace Es(t) and the estimated time trace Eest(t),

[0194] until a set of values of the physical parameters (pi opt ) is obtained that minimizes the error function.

[0195] Another aspect according to the application relates to a spectrometer Spectro as shown in Figure 15 comprising i) a measurement device MeD comprising a light source LS configured to illuminate a sample S by an excitation beam EB, and a detector D configured to detect a measured sample time trace Es(t) and a measured reference time trace Eref(t); and ii) a characterization device 50.

[0196] It should be understood that the previous embodiments are only non-limiting examples. In particular, the measurement device MeD and the characterization device 50 can be located in different pieces and used together in any combination.

[0197] In one embodiment, such as Figure 16 As shown, the MeD can be connected to a computer Comp via I / O interface 703, indicating that the device 50DBP is located in that computer. In another embodiment, as... Figure 17 As shown, the MeD can be connected to the computer Comp via I / O interface 703, and the characterization device 50 can be located in a remote server 76 of a communication subsystem 720 connected to the computer Comp via the Internet 75.

[0198] On the other hand, the present invention relates to a computer program suitable for implementing the steps of the claimed method. On the other hand, the present invention relates to a computer-readable medium including said computer program.

Claims

1. Method for determining a set of physical parameters of a sample, comprising the following steps: - A. acquiring a measured sample temporal trace Es(t), said measured sample temporal trace Es(t) having been obtained by time-domain spectroscopy by illuminating a sample (S) by means of an excitation beam (EB) periodically emitting electromagnetic pulses having a period T and presenting a comb-like frequency, and by detecting the electromagnetic field coming from said sample as a function of time by means of coherent detection, the time period during which said sample temporal trace is measured being tmax, with tmax < T, - B. acquiring a measured reference temporal trace Eref(t), said measured reference temporal trace Eref(t) having been obtained by illumination and detection in the same conditions as step A but in the absence of said sample, - determining a broadened reference time trace called Eref0(t) extending over said period T and obtained by assuming a 0 value for the instants where no measurement is performed, and determining a discrete Fourier transform of said broadened reference time trace computed over a time window equal to T , - D according to a set of physical parameters from a Fourier transform of said widened reference time trace and a physical behavior model of said sample, determining a modeling of the impulse response of said sample in the frequency domain, called sample frequency model ​ - E applies an optimization algorithm to the set of physical parameters including the following sub-steps: - E1 initializes physical parameters , - iteratively implementing the following sub-steps: - E2 calculates the sample frequency model of the estimated sample time trace the inverse discrete Fourier transform of the estimated sample time trace,​ - E3 calculates an error function from the difference between the measured sample time trace Es(t) and the estimated time trace Eest(t) until a set of values of the physical parameters is obtained that minimizes the error function.

2. The method of claim 1, wherein, said excitation beam (EB) being in the THz domain having a frequency comprised between 100 GHz and 30 THz.

3. The method of claim 1 or 2, wherein, The maximum time delay tmax is chosen so as to comprise more than 95% of the energy of said measured reference temporal trace.

4. The method of claim 1 or 2, wherein, The sample frequency model In the Fourier transform With the transfer function characterizing the behavior of the sample Is multiplied.

5. The method of claim 4, wherein, The transfer function Depending on the complex refractive index .

6. The method of claim 5, wherein, For each spectral line, the square of the complex refractive index, called the dielectric constant follows a Drude-Lorentz model, the spectral line being characterized by a set of three parameters: an amplitude (M), a width, called the damping rate and a central frequency .

7. The method of claim 1 or 2, wherein, said error function is defined as: 。 8. Characterization device (50) for characterizing a sample (S), said device comprising: - a memory (MEM) storing a measured sample temporal trace Es(t) and a measured reference temporal trace Eref(t), said measured sample temporal trace Es(t) having been obtained by time-domain spectroscopy by illuminating said sample (S) by means of an excitation beam (EB) periodically emitting electromagnetic pulses having a period T and presenting a comb-like frequency, and by detecting the electromagnetic field coming from said sample as a function of time by means of coherent detection, the time period during which said sample temporal trace is measured being tmax, with tmax < T, said measured reference temporal trace Eref(t) having been obtained by illumination and detection in the same conditions as for the measured sample temporal trace Es(t) but in the absence of said sample, - a processing unit (PU) configured to: - determining a widened reference time trace, called Eref0(t), extending over said period T and obtained by assuming as 0 value the instants where no measurement is performed, and determining the discrete Fourier transform of said widened reference time trace computed on a time window equal to T , -- a set of physical parameters from a Fourier transform of the widened reference time trace and a physical behavior model of the sample, a modeling of the impulse response of the sample in the frequency domain, called sample frequency model is determined, - on physical parameters An optimization algorithm is applied, comprising the following steps: --- initialize physical parameters, --- iteratively implement the following sub-steps: computing the sample frequency model the inverse discrete Fourier transform of the sample time trace referred to as estimation calculating an error function from the difference between the measured sample time trace Es(t) and the estimated time trace Eest(t) until a set of values of the physical parameters is obtained that minimizes the error function.

9. Spectrometer (Spectro) comprising: - a characterization device (50) according to claim 8, - a measurement device (MeD) comprising: - a source (LS) configured to illuminate a sample (S) by means of said excitation beam (EB), - a detector (D) configured to detect said measured sample temporal trace Es(t) and said measured reference temporal trace Eref(t).

10. Computer program product comprising a computer program adapted to implement the steps according to any one of claims 1 to 7.

11. Computer-readable medium comprising the computer program product according to claim 10.

Images