Spurious Absorption Frequency Appearance Due to Frequency Conversion Processes in Pulsed THz TDS Problems.

The appearance of the spurious absorption frequencies caused by the frequency conversion process at the broadband THz pulse propagation in a medium is theoretically and experimentally discussed. The spurious absorption frequencies appear due to both the frequency doubling and generation of waves with sum or difference frequency. Such generation might occur because of the nonlinear response of a medium or its non-instantaneous response. This phenomenon is confirmed by the results of a few physical experiments provided with the THz CW signals and broadband THz pulses that are transmitted through the ordinary or dangerous substances. A high correlation between the time-dependent spectral intensities for the basic frequency and generated frequencies is demonstrated while using the computer simulation results. This feature of the frequency conversion might be used for the detection and identification of a substance.


Literature review and introduction
The THz time-domain spectroscopy (THz TDS) is a non-destructive tool that is very popular nowadays for both basic research and the security screening and for industrial applications also. Is well known that the THz pulsed TDS is traditionally used for the detection and identification of explosive and illicit drug [1][2][3][4][5][6][7][8][9][10]. Moreover, it is currently considered to be a new powerful tool for application in the chemical sciences and molecular structure investigations [11][12][13][14], nondestructive inspections [15][16][17][18], thickness measurements [19,20], art investigations [21,22], pharmaceutical and biomedical applications [23][24][25][26], and food quality control [27][28][29]. In [30], a comprehensive review of the history and state of the art of the GaN-based devices for THz applications is provided. One important limitation of the most THz sources and detectors should also be noted. Namely, they are restricted by their low conversion efficiency. Great efforts have therefore been made to improve the efficiency of THz sources and detectors [31,32].
As a rule, substance identification is carried out based on a comparison of the absorption frequencies of a substance under investigation with a set of known absorption frequencies from a database. We call this identification technique the standard THz TDS method. Obviously, this method has well-known disadvantages. For example, if the object is masked by some non-opaque covering or hidden under the packaging material, then the spectrum of the measured THz signal might be significantly distorted [33,34]. Moreover, such factors as inhomogeneity of the substance surface and atmospheric humidity (and many other factors) also significantly distort the spectra of the measured THz pulse [35][36][37][38]. Moreover, the spurious artifacts associated with the experimental setups and with TDS THz system as well as with signal processing and sample geometries, also bring difficulties at measuring of the optical properties of the substance under investigation [39][40][41][42]. In all of these cases, the substance identification by means of the standard THz TDS is ineffective.
In addition, some non-opaque packaging materials (paper, plastic material, clothes, etc.) usually cover an examined substance [43]. Often, these covers possess both density fluctuations and sample thickness variations in the sub-millimeter scale of length and, therefore, they can be considered to be a disordered photonic structure in the THz frequency range. Consequently, one can appear the false absorption frequencies, which are caused by inhomogeneities, in the spectrum of the broadband THz pulse transmitted through or reflected from these ordinary materials and these frequencies may coincide with those of the dangerous substances [44][45][46]. Thus, the appearance of the false absorption frequencies under the THz TDS is an ordinary phenomenon and recognizing of the true absorption frequencies inherent to a substance under consideration is an urgent problem. We have also proposed an efficient tool for overcoming of influence of most of these factors on the detection and identification of a substance. Our approach is based on an analysis of the time-dependent spectral intensities (SDA-method) on both the absorption frequencies and the emission ones. Earlier, the SDA-method was successfully applied for the detection and identification of neutral and dangerous substances under real conditions [47][48][49][50][51][52][53][54], including the cases of very noisy THz signal and a sample with inhomogeneous surface. Therefore, we also follow this approach below.
There are also other physical mechanisms, which are not associated with the factors mentioned above and they lead to observing of the spurious absorption frequencies in the spectrum of the pulse transmitted through or reflected from a substance. They may appear due to the frequency conversion processes at resonant or non-resonant interaction of the THz pulse with a substance or due to non-instantaneous response of a medium. One of them is the cascade mechanism of the high energy level excitation of molecules [55,56]. It leads to the generation of the frequencies not belonging to the incident pulse spectrum. Another one, which is considered below, is a frequency conversion, because of the nonlinear response of a medium. For example, the second harmonic generation (SHG), or generation of waves with the frequencies, which are a linear combination of the frequencies that belong to the incident pulse spectrum, may appear due to the quadratic nonlinear response of a medium. It is important to stress that the frequency conversion or interaction of the Fourier harmonics of a wave packet can simultaneously occur because the incident THz pulse possesses broadband spectrum and the pulse interaction with a medium is essentially non-instantaneous. Moreover, a part of them might be observed after the action of a THz pulse transmitted through/reflected from a medium because of the possible delay of the emission of a medium. This feature distinguishes the frequency conversion of the broadband THz pulse from processes occurring in nonlinear optics, where it occurs for monochromatic optical wave, or from frequency changing when the monochromatic wave interacts with two-level medium, for example.
It should be stressed that the spurious absorption frequencies that appear under the frequencies up-conversion process (for example, at SHG) may be used as a new effective tool for the detection and identification of the substance. We illustrated such a possibility in [57,58], where a description of the nonlinear response of a medium was provided while using density matrix formalism [59] and the THz pulse propagation is described by one-dimensional (1D) Maxwell's equations. A single substance layer was placed between two coverings, which consisted of linear layers with random dielectric permittivity. The appearance of the doubled absorption frequency under non-stationary interaction of the THz broadband pulse with multi-levels molecules was demonstrated in both transmission and reflection mode. Analyzing time-dependent spectral intensity at the doubled frequency, we showed a possibility for the substance detection and identification despite covering with random parameters.
In this paper, we demonstrate the appearance of the spurious absorption frequencies under the frequency conversion process that is caused by nonlinear (quadratic) response of a medium. In contrast to the previous papers, in the theoretical analysis, we consider the instantaneous response of a medium under the action of different types of signals consisting of several sub-pulses with narrow and broad spectra. Such a signal is a typical one for the pulsed THz TDS of the substance made similar to a tablet. A high correlation between the time-dependent spectral intensities for the basic frequency and generated frequencies (doubled frequency, or sum frequency, or difference frequency) is demonstrated for various types of the signals. This is the basis for using spurious absorption frequencies or the emission frequencies as the frequencies for the detection and identification of a substance.
We discuss the results of a few physical experiments carried out with broadband THz pulses and narrowband CW radiation, transmitted through the ordinary materials, to verify our theoretical conclusions about the appearance of the spurious absorption frequencies caused by the frequency conversion processes due to the nonlinear or non-stationary process and to demonstrate their appearance (soap, paper, chocolate, and plastic). It should be stressed that the experiments with ordinary materials were carried out in the room atmosphere with relative humidity of about 30%. It means that the measured THz signals are noisy, and this brings our investigation close to reality. Using the standard THz TDS method under such a condition can lead to a large number of false positive alarms during security screening, which renders its effectiveness low. At the same time, a possibility of the substance detection and identification by means of SDA-method together with several integral correlation criteria (ICC) using noisy THz signals with SNR about 0.5 was demonstrated in [47,48,51]. Additionally, we use this technique in the current paper.
The measurements for narrowband THz CW radiation transmitted through the samples with aspirin, soap, paper, and chocolate were provided in the University of Dayton (Dayton, OH, USA). The measurements for broadband THz pulse transmitted through the polyethylene and paper were made in the Military University of Technology (Warsaw, Poland) and in the South China Normal University (Guangzhou, China), respectively. The measurements for the broadband THz pulse interaction with 2,4-DNT were made at the Capital Normal University (Beijing, China). The physical experiments with transmission of the broadband THz pulses through the soap and chocolate were measured in the Institute of Spectroscopy RAS (Moscow, Troitsk, Russian).
We use the results of several physical experiments performed at different times and in different labs to illustrate that the conclusions do not depend on particularities of the THz setup. We stress that for the measurements of the broadband THz pulse the similar setups in various labs were used. We specify their types and characteristics in Section 3.

Theoretical Analysis
As we mentioned in the Introduction, the process of a broadband THz pulse interaction with a substance because of both non-instantaneous response of a medium, and its nonlinear response, as well as its resonance response, can lead, in particular, to the appearance of both the spurious absorption frequencies and the substance emission frequencies in the spectrum of the pulse transmitted through a medium. Let us remind that frequency conversion efficiency depends on a spectral brightness of the corresponding Fourier harmonics ν 0 . Therefore, if we consider a process of the frequency doubling under the propagation of the broadband pulse and the frequency is the absorption one of a substance, then the spectrum of the pulse, which is transmitted through a medium, contains a minimum near the frequency 2ν 0 (see Figure 1a, an accordance between the intensities is schematically depicted by arrows). Thus, we will observe the induced absorption frequency, if the spectral resolution of the measurement device is high enough.
In Figure 1 the scheme of verification of the substance absorption frequencies ν 0 (the basic frequency) and 2ν 0 (the doubled frequency) (whether it is true or false), as observed in the spectrum of the incident broadband pulse, transmitted through the substance, is shown. With this aim, the narrowband pulse should be used (black line in Figure 1a). While shifting its current frequency along the frequency axis, one can observe the presence or absence of the pulse energy absorption at the frequency 2ν 0 . If the pulse energy absorption is present at the frequency (Figure 1a, ν 0 , 2ν 0 ), then it Sensors 2020, 20, 1859 4 of 36 means that we observe the true absorption frequency of a substance. If there is no minimum in the spectrum at the doubled frequency 2ν 0 at using the probing narrowband THz CW radiation (Figure 1b), then the absorption frequency 2ν 0 occurring in Figure 1a at the broadband pulse transmission through a medium is induced by the frequency doubling process. Therefore, this frequency is the spurious one.
Sensors 2020, 20, x FOR PEER REVIEW 4 of 37 frequency 0 2ν . If the pulse energy absorption is present at the frequency (Figure 1a, , 0 2ν ), then it means that we observe the true absorption frequency of a substance. If there is no minimum in the spectrum at the doubled frequency 0 2ν at using the probing narrowband THz CW radiation ( Figure   1b), then the absorption frequency 0 2ν occurring in Figure 1a at the broadband pulse transmission through a medium is induced by the frequency doubling process. Therefore, this frequency is the spurious one. Before the discussion of results of the theoretical analysis of an appearance of the spurious absorption (or, may be, emission) frequencies we make two remarks. Firstly, for definitely, below we consider a transmission of the incident pulse through a medium with the quadratic nonlinear response in the approximation of the so-called optical thin layer. In this case, at the supposing of the instantaneous response of a medium, we can write the following expression: for the electric field strength of the pulse transmitted through a medium. A dimensionless parameter (2) χ characterizes the quadratic susceptibility of a medium. To illustrate the appearing of the basic frequency conversion, we choose this parameter equal to unity. However, in practice, this parameter can be significantly less than unity; its value depends on both the properties of a medium and the incident pulse maximal intensity. Secondly, we consider two types of the artificial signal while taking into account the structure of the real THz pulse transmitted through the layer of a medium. Depending on its thickness, this pulse can contain separate or overlapping sub-pulses.

Frequency Conversion Near the Spectral Maximum of Artificial Signals
Let us consider the artificial incident pulse ( ) inс E t that is the sum of four terms and has the following structure: The absorption frequencies ν 0 , 2ν 0 occurring in the spectrum of the incident broadband pulse, transmitted through a substance, are seen in the spectrum of the incident narrowband signal (E_tran), transmitted through this substance, at changing its carrier frequency (a). Thus, the frequency 2ν 0 is a true absorption frequency. (b) demonstrates that the frequency 2ν 0 is induced by the up-conversion process (the spurious frequency), because it is not observed in the spectrum of the incident narrowband signal transmitted through a substance. Black line denotes the spectrum of the incident narrowband signal (E_inc).
Before the discussion of results of the theoretical analysis of an appearance of the spurious absorption (or, may be, emission) frequencies we make two remarks. Firstly, for definitely, below we consider a transmission of the incident pulse through a medium with the quadratic nonlinear response in the approximation of the so-called optical thin layer. In this case, at the supposing of the instantaneous response of a medium, we can write the following expression: for the electric field strength of the pulse transmitted through a medium. A dimensionless parameter χ (2) characterizes the quadratic susceptibility of a medium. To illustrate the appearing of the basic frequency conversion, we choose this parameter equal to unity. However, in practice, this parameter can be significantly less than unity; its value depends on both the properties of a medium and the incident pulse maximal intensity. Secondly, we consider two types of the artificial signal while taking into account the structure of the real THz pulse transmitted through the layer of a medium. Depending on its thickness, this pulse can contain separate or overlapping sub-pulses.

Frequency Conversion Near the Spectral Maximum of Artificial Signals
Let us consider the artificial incident pulse E inc (t) that is the sum of four terms and has the following structure: Here, A ki is the amplitude of the corresponding sub-pulse, t ki is its center, and τ ki is the half-duration of the sub-pulse, where i = 1, 4. The carrier frequencies ν k are fixed for all artificial signals and have the values: ν 1 = 0.8, ν 2 = 1.2, ν 3 = 1.4, and ν 4 = 1.8. We will call them as the basic frequencies. In Table 1, the parameters of two artificial signals E inc (t) are shown. Let us note that a time t and the frequency ν and A ki are dimensionless variables.  The pulse E inc (t), its spectrum, and the spectra of the signals E inc 2 (t) and E inc (t) + E inc 2 (t) are shown in Figure 2 for different values of parameters A ki , t ki , τ ki , as depicted in Table 1.
A is the amplitude of the corresponding sub-pulse, ki t is its center, and ki τ is the halfduration of the sub-pulse, where 1, 4 i = . The carrier frequencies k ν are fixed for all artificial signals and have the values: 1 0.8 ν = , 2 1.2 ν = , 3 1.4 ν = , and 4 1.8 ν = . We will call them as the basic frequencies. In Table 1, the parameters of two artificial signals are shown. Let us note that a time t and the frequency ν and ki A are dimensionless variables.  The pulse ( ) inс E t , its spectrum, and the spectra of the signals Figure 2 for different values of parameters ki A , ki t , ki τ , as depicted in Table 1.
Artificial pulses E inc (t) (a,e), their spectra (b,f), and the signal E inc 2 (t) spectra (c,g) and the It is easy to see that the artificial signal E inc (t) (a) consists of four intersecting sub-pulses with three narrow spectra and one broad spectrum with the pulse duration τ 12 = 2. The sub-pulses spectra also partially overlap. Nevertheless, the signal E inc (t) spectrum (b) possesses four pronounced maxima at the basic frequencies and small maximum near the frequency 1.5. In the signal E inc there are maxima at the doubled basic frequencies ν = 1.6, 2.4, 2.8, 3.6, due to the quadratic nonlinear response of a medium. Moreover, one can also see the spectral maxima at the frequencies ν = 0.4, 0.6, 2.2, 3.2, which are the linear combinations of the basic frequencies (b). An appearance of these frequencies in the spectrum is discussed below. In (c), we can also observe the broadening of the incident pulse spectrum in the high (ν > 2.0) and low (ν < 0.6) frequency ranges. The second artificial signal E inc (t) (e) also consists of intersecting sub-pulses, but their spectra essentially overlap that results in appearance of the additional spectral maximum at the frequency ν = 0.96 (b1). The signal E inc 2 (t) spectrum (g) obviously contains maxima at the doubled frequencies ν = 1.6, 2.4, 2.8. At the same time, in (g), the maximum at the doubled frequency ν = 3.6 is absent, and one can see additional spectral maximum at the frequency ν = 3.28. This is a result of a competition Sensors 2020, 20, 1859 6 of 36 between processes of SHG (ν = 3.6 = 1.8·2) and sum frequency wave generation (SFG) ν sum = 3.2 = 1.8 + 1.4. The additional spectral maxima can be also observed as a result of the difference frequency wave generation (DFG): ν di f = 0.4 = 1.2 − 0.8, 0.6 = 1.8 − 1.2, 2.2 (b1). However, some of them can be obtained due to the generation of waves with half basic frequencies. For example, the frequency ν = 0.4 might occur as a generation of a wave with the frequency ν = 0.8/2. (Figure 2(g)). We will denote these frequencies as ν sum , ν di f , respectively, and ν prod -as a linear combination of the basic frequencies. Thus, various processes of wave generation give contribution to the appearance of a wave at certain frequency.
Obviously, the spectral maxima of the signal E inc (t) + E inc 2 (t) (Figure 2(d,h) coincide with the maxima of the corresponding signals E inc (t) (b,b1) and E inc 2 (t) (c,g). Therefore, in the following Sections, we will analyze the evolution of the time-dependent spectral amplitudes (spectral line dynamics) at the basic and doubled frequencies as well as at the additional emission frequencies for the signals E inc (t) and E inc 2 (t). The correlation between these spectral intensities is of most interest for us. In this Section, we investigate the evolution of the spectral amplitudes-square root from spectral intensities-at the basic, doubled, and other emission frequencies for the signals E inc (t) and E inc 2 (t) (Example 1). The main aim of this investigation is a demonstration of high similarity between the spectral dynamics at these frequencies.
As an example, the modulus of such spectral amplitude |P ν (t)|, |P 2ν (t)| evolution (their computation is described in Appendix A and in [44][45][46]) at the frequencies Figure 3. A high degree of similarity between time-dependent spectral amplitudes is observed in all cases (a-d), and we may expect their high correlation. Indeed, the corresponding correlation coefficient between them are close to unity: Here we use the well-known expression: where P ν 1 , P ν 2 are the non-zero mean values of the spectral amplitudes |P ν 1 | and |P ν 2 |. Subtracting of these values from the spectral dynamics in expression (2) is necessary to eliminate their influence on the correlation coefficient c ν 1 ,ν 2 . Here M denotes the number of time moments in the spectral dynamics |P ν 1 |, |P ν 2 |, which depends on the dynamics construction parameters-the window length T and its shift ∆ along the signal under investigation. With accordance to our previous articles ( [44][45][46]), these parameters are chosen as T = 2.8 and ∆ = 0.2 dimensionless units. The DFG may be obtained by two pairs of Fourier harmonics interaction: between the harmonics with frequencies ν = 1.2 and ν = 0.8; as well as ν = 1.8 and ν = 1.4, respectively: dif Each of these pairs might give a contribution to the frequency dif ν = 0.4 generation. To clarify The DFG may be obtained by two pairs of Fourier harmonics interaction: between the harmonics with frequencies ν = 1.2 and ν = 0.8; as well as ν = 1.8 and ν = 1.4, respectively: ν di f = 0.4 = 1.2 − 0.8 = 1.8 − 1.4. Each of these pairs might give a contribution to the frequency ν di f = 0.4 generation. To clarify which contribution occurs and in which time interval, we depict the spectral dynamics at the emission frequency ν = 0.4 together with the spectral dynamics at the basic frequencies ν = 0.8, 1.2, and ν = 1.4, 1.8 in Figure 4a,b, respectively. Table 2 provides the corresponding correlation coefficients c ν base ,ν di f (Example 1). The maximal correlation occurs between the DFG with ν di f = 0.4 and the spectral line dynamics computed at the basic frequencies: The DFG may be obtained by two pairs of Fourier harmonics interaction: between the harmonics with frequencies ν = 1.2 and ν = 0.8; as well as ν = 1.8 and ν = 1.4, respectively: dif ν = 0.4 = 1.2 − 0.8 = 1.8 − 1.4. Each of these pairs might give a contribution to the frequency dif ν = 0.4 generation. To clarify which contribution occurs and in which time interval, we depict the spectral dynamics at the emission frequency ν = 0.4 together with the spectral dynamics at the basic frequencies ν = 0.8, 1.2, and ν = 1.4, 1.8 in Figure 4a,b, respectively. Table 2 provides the corresponding correlation coefficients  When comparing the spectral line dynamics, one can see that the contribution of various basic frequencies to generation process changes in different time intervals. Accordingly, one can see in When comparing the spectral line dynamics, one can see that the contribution of various basic frequencies to generation process changes in different time intervals. Accordingly, one can see in Figure 4 that the generation of the DFG with ν di f = 0.4 due to the interaction of the harmonics at the basic frequencies ν = 0.8, 1.  Table 2. Correlation coefficients c ν base ,ν di f between the basis frequency ν = v base and the difference frequency ν = ν di f or the sum frequency ν = ν sum . Empty squares mean that the corresponding values of the correlation coefficient were not computed. For further confirmation, we consider the DFG with frequency ν di f = 0.6, which is the linear combinations of three generation processes: (i) frequency generation via interaction of two basic (carrier) Sensors 2020, 20, 1859 8 of 36 frequencies ν base = 1.4, 0.8; (ii) similar process with other pair of carrier frequencies ν base = 1.8, 1.2; and, (iii) generation of wave with half frequency from the basic one ν base = 0.6 = 1.2/2. Each of these processes may contribute to the frequency ν di f = 0.6 generation process. Additionally, we consider the SFG ν sum = 3.2 = 1.8 + 1.4. Figure 5 depicts the evolution of modulus of the spectral amplitudes of all frequencies ν = 0.6. We see the same features of the generation processes as in the previous case depicted in Figure 4. Therefore, different frequencies play determining role in various time intervals. The maximal correlation coefficients are for the spectral line dynamics at the basic frequency ν base = 0.8 (a), 1.2 (b) and the difference frequency ν di f = 0.6. They are equal to c 0.8,0.6 = 0.839, c 1.2,0.6 = 0.688 ( Table 2, Example 1), respectively. Thus, the half frequency generation from the basic one ν base = 1.2 can also contribute to the frequency ν di f = 0.6 generation process. This Section is similar to the previous one. However, a principal distinction of the signal under consideration from previous one lies in the structure of the pulse: it consists of three broadband subpulses and a sub-pulse with narrow spectrum (Example 2, Table 2). Therefore, the pulse spectrum differs markedly from the spectrum of the pulse investigated in Section 2.1.1 (see Figure 2 (e,g)). One can see that the spectral intensity of the signal 2 ( ) inc E t spectrum insignificantly changes in the frequency range ν= [3.0, 3.8] units. This is a consequence of the spectrum overlapping of two subpulses. In this case, we do not see a pronounced maximum at the doubled frequency ν= 3.6 in Figure  2 (g). Instead, the maximal spectral intensity is seen at the frequency ν= 3.28.

Example
In Figure 6, the evolution of the spectral amplitude at the basic and doubled frequencies, and at the shifted frequency ν = 3.28 (e) is shown. One can see that an appearance of the wave with doubled frequency is delayed with respect to basic wave in any case and its trailing edge vanishes early than the spectral intensity at the basic frequency. This fact is in good agreement with the theory of SHG.
Computed correlation coefficient between the corresponding spectral line dynamics demonstrates very high value:  In Figure 5c, we see a high correlation between the spectral line dynamic at the basic frequency ν base = 1.8 and the sum frequency ν sum = 3.2 if a time is less than 40 units. However, if a time is greater than 50 units, then an evolution of the spectral line dynamics at the basic frequency ν base = 1.4 defines a process of the sum frequency generation. However, in the time interval t = [35,40] units, the sum frequency generation is defined by the spectral intensity evolution at the basic frequency ν base = 1.8. Nevertheless, the correlation coefficient between the spectral line dynamics at the basic frequency ν base = 1.4 and the sum frequency ν sum = 3.2 is equal to c 1.4,3.2 = 0.66.

Correlation between the Spectral Line Dynamics at the Basic and Doubled Frequencies for the Signal Consisting of Sub-Pulses with Three Broad Spectra and One Narrow Spectrum (Example 2)
This Section is similar to the previous one. However, a principal distinction of the signal under consideration from previous one lies in the structure of the pulse: it consists of three broadband sub-pulses and a sub-pulse with narrow spectrum (Example 2, Table 2). Therefore, the pulse spectrum differs markedly from the spectrum of the pulse investigated in Section 2.1.1 (see Figure 2 (e,g)). One can see that the spectral intensity of the signal E 2 inc (t) spectrum insignificantly changes in the frequency range ν = [3.0, 3.8] units. This is a consequence of the spectrum overlapping of two sub-pulses. In this case, we do not see a pronounced maximum at the doubled frequency ν = 3.6 in Figure 2 (g). Instead, the maximal spectral intensity is seen at the frequency ν = 3.28.
In Figure 6, the evolution of the spectral amplitude at the basic and doubled frequencies, and at the shifted frequency ν = 3.28 (e) is shown. One can see that an appearance of the wave with doubled frequency is delayed with respect to basic wave in any case and its trailing edge vanishes early than the spectral intensity at the basic frequency. This fact is in good agreement with the theory of SHG.
Computed correlation coefficient between the corresponding spectral line dynamics demonstrates very high value: c ν,2ν = 0.978 (a), 0.917 (b), 0.946 (c), 0.884 (d). It should be stressed that, for the case (e), we obtain the following correlation coefficient value c 1.8,3.28 = 0.881 and this value is very close Sensors 2020, 20, 1859 9 of 36 to c 1.8,3.6 = 0.884. In our opinion, it means that two sub-pulse with overlapping spectra give their contributions to the process of SHG at the basic frequency ν = 1.8 and there is a generation of frequency band simultaneously.
Note that the spectral amplitudes modulus P ν (t), P 2ν (t) are computed according Equations (A1) and (A2) from Appendix A and are presented in Figure 6 in the non-normalized form. Therefore, it is not evident why the correlation coefficient for the spectral line dynamics depicted in (c-e) is so high. To clarify this fact, in Figure 6b,c, the corresponding spectral line dynamics are depicted in the additional panels that are normalized in L2-norm for convenience. In these additional panels, the high similarity in the shapes of the normalized spectral dynamics is evident and, therefore, the high values of the corresponding correlation coefficients are true. The same result takes place for the spectral line dynamics in (d), (e).
We see a high correlation between time-dependent spectral intensities at chosen frequencies:   Table 2 (Example 2).
We see a high correlation between time-dependent spectral intensities at chosen frequencies: Figure 7, we clearly see that these carrier frequencies give the main contribution to the generation process during time interval t ≈ [10,20] units. Outside of this interval, the frequency generation occurs due other frequency mentioned above. Similar results were obtained for the artificial signal E inc (t), which consists of four narrow sub-pulses with non-overlapping spectra. Due to the simplicity of this case, no figures are given. inc Figure 7 shows the spectral line dynamics at the difference frequency dif ν = 0.4 = 1.2−0.8 = 1.8−1.4 (a), dif ν = 0.6 = 1.4−0.8 = 1.8−1.2 (b), and sum frequency sum ν = 2.2 = 0.8 + 1.4 (c), and basic frequencies ν = 0.8, 1.4, because these frequencies demonstrate the maximal correlation coefficients. All of their values are shown in Table 2 (Example 2).
We see a high correlation between time-dependent spectral intensities at chosen frequencies:  (black line, signal E inc (t)) and difference frequencies ν di f = 0.4 (a), 0.6 (b) and sum frequency ν sum = 2.2 (c) (red line, signal E inc 2 (t)), respectively.

Frequency Conversion Near Minimum of the Pulse Spectrum
Below, we investigate the frequency conversion if the frequency lies near the spectrum minimum (Examples 1 and 2 in Table 2).

The Artificial Signals E inc (t) Consisting from the Sub-Pulses (Examples 1 and 2)
We see in the artificial signal E inc (t) spectra, as depicted in Figure 2 (b,f), several minima at the basic frequencies ν base = 0.96, 1.345, 1.46, 1.675 (b), 1.06, 1.31, and 1.53 (f). These minima are a consequence of the pulse structure: the pulse under consideration consists from the four sub-pulses. At the same time, in the signals E inc 2 (t) and E inc (t) + E inc 2 (t) spectra (see Figure 2) one can observe the minima at the frequencies close to the corresponding doubled frequencies: ν = 1.8, 2.74, 2.86, 3.34 (c), (d), and 2.1, 2.67, 3.05 (g), (h). We see that, in some cases, the difference between values for doubled frequencies and observed frequencies might be significant. This is a result of the pulse structure: it consists of sub-pulses. The corresponding spectral line dynamics are shown in Figure 8 for the signals E inc (t) and E inc 2 (t) computed for the Example 1 (a)-(c) and Example 2(d-f). We see high similarity between these spectral line dynamics. This fact is confirmed by the correlation coefficients, which are shown in Table 3.
In most cases, they are very high: c ν,2ν ≥ 0.9 (a), (b), (d)-(f). Maximal correlation coefficient occurs for time-dependent spectral intensities (e), (f): c ν,2ν = 0.94 (e), 0.97 (f). As above (see Figure 6), the inset in (f) shows the same spectral line dynamics normalized in L2-norm, to clarify the high values of the corresponding correlation coefficients in (d-f). The additional panel in (f) confirms the high similarity in the shapes of the normalized spectral dynamics and, therefore, the high values of the corresponding correlation coefficients. The same result is valid for the spectral line dynamics in (d), (e). The dip in the spectral line dynamics that was computed at the doubled frequency ν = 2.86 and depicted in Figure 8b was caused by the pulse structure: the presence of the sub-pulses. The same results occur for the signal E inc (t) + E inc 2 (t) The spectral amplitude evolution computed for the doubled frequency is very similar to those computed for the basic ones. Thus, they (as well as the emission frequencies) can be used for the detection and identification of the substances, because the THz radiation at the frequencies belonging to the high-frequency range can propagate through the disordered structure with minimal distortion. base consequence of the pulse structure: the pulse under consideration consists from the four sub-pulses. At the same time, in the signals Figure 2) one can observe the minima at the frequencies close to the corresponding doubled frequencies: ν = 1.8, 2.74, 2.86, 3.34 (c), (d), and 2.1, 2.67, 3.05 (g), (h). We see that, in some cases, the difference between values for doubled frequencies and observed frequencies might be significant. This is a result of the pulse structure: it consists of sub-pulses.   Table 3. In most cases, they are very high:  In this Section, we consider the artificial signal E inc,min (t) (Example 3), as depicted in Figure 9a, which does not possess the sub-pulse structure and contains two minima in its spectrum at the frequencies ν base = 0.5, 0.8 ( Figure 9b). Let us notice that such a structure of signal is present as a rule if we only consider the main THz pulse reflected from or transmitted through a medium. Below, we briefly describe a way of obtaining such a signal. For this aim, we consider the Gaussian pulse: (3) on the time interval t = [0, 100] with the time step h t , make the discrete Fourier transform and subtract from the sineand cosine-Fourier transforms of the discrete function E base (t i ), i = 1, . . . , N base , N base = 100/h t + 1 the discrete functions G k (ν l ), where: Here, the parameter max|P b | denotes the maximal modulus of the spectral amplitude of the pulse E base (t) spectrum at the frequency ν k . The parameter B k describes the depth of the minimum at the frequency ν k and it is equal to the following values: B 1 = 0.95, B 2 = 0.3. The parameter ∆ν k characterizes the absorption line spectral width and is equal to ∆ν 1 = 0.1, ∆ν 2 = 0.05. After the inverse Fourier transform, we obtain the pulse, being denoted as E inc,min (t), with the spectral minima at the given frequencies ν 1 = 0.5 and ν 2 = 0.8. Additionally, in the signal E inc,min (t) spectrum that is depicted in Figure 9b, one can see three maxima at the frequencies ν = 0.41, 0.63, 0.95.
Here, the parameter max | | b P denotes the maximal modulus of the spectral amplitude of the In each of the spectra, an emission at the doubled frequency ν = 0.5·2 is present. The spectral line dynamics computed at the basic frequencies that correspond to minimal spectrum intensities and at the frequencies, which are produced by interaction of these spectral harmonics, are depicted in Figure 9e,f. We see that the main contribution gives the spectral harmonic at the frequency ν = 0.8. The corresponding correlation coefficients are shown in Table 4. All of them are high and greater than 0.9.  In each of the spectra, an emission at the doubled frequency ν = 0.5·2 is present. The spectral line dynamics computed at the basic frequencies that correspond to minimal spectrum intensities and at the frequencies, which are produced by interaction of these spectral harmonics, are depicted in Figure 9e,f. We see that the main contribution gives the spectral harmonic at the frequency ν = 0.8. The corresponding correlation coefficients are shown in Table 4. All of them are high and greater than 0.9. Thus, we see that the spurious absorption frequencies that belong to the spectrum of the transmitted broadband pulse, which initially does not contain a sub-pulse structure, may appear in complicated way. They are more pronounced if the generated frequencies lie out of the incident pulse spectrum. We also see induced spurious absorption frequencies, which are the linear combinations of various frequencies. They have high similarity with the true absorption frequencies of a substance. This fact allows for us to use both the emission frequencies and the induced spurious absorption frequencies for the detection and identification of substances instead of the true absorption frequencies of a medium. Obviously, it is necessary to take the spectral intensity of generated frequencies into account, which depends on the incident pulse power density, in particular.

Experimental Setups
Obviously, from provided consideration the following question arises: how can we verify the effect of the spurious absorption frequency generation in physical experiment? We mentioned above that it could be done by using two types of the physical experiments providing with narrowband and broadband THz signals transmitted through the substance.
In the first experiment, the absorption frequencies of 2,4-DNT sample are investigated. The measurements of a broadband THz pulse transmitted through the sample with 2,4-DNT were made using a Ti:sapphire regenerative amplifier (Spectra-Physics Mai Tai) delivering optical pulses with a duration of 100 fs and a central wavelength of 800 nm at a repetition rate of 1 kHz, with noise less than 0.15% [48,[60][61][62]. The transmitted THz signal was detected by free-space electro-optic sampling in a 1-mm-thick <110> ZnTe crystal with the sampling pulse. The DNT sample was 0.5 mm thickness.
In the second experiment, we investigate the ordinary materials-aspirin, soap, paper, chocolate, and polyethylene bag. The measurements were made using a narrowband CW THz radiation spectrometer Spectra 400 (TeraView Ltd., Cambridge, UK) in the room atmosphere. The spectral range, in which the measurements were provided, is 0.15-1.6 THz, frequency resolution is 100 MHz, signal-to-noise in 100-1000 GHz range is less than 60 db [63]. The frequency step at the measurement is equal to ∆ν = 0.001 THz. Obviously, the spectral resolution is the same. The distance between the sample and the receiver does not exceed 20 cm. Below we will call the spectra measured while using the narrowband CW THz radiation as narrowband spectra for brevity.

Results of the Experiment Using a Broadband THz Pulse
There are a number of papers in which the substance absorption frequencies were measured in transmission mode using the broadband THz pulse. For example, in [6], the THz spectra of 15 ERCs were measured in transmission mode with a broadband THz pulse and FTIR spectrometer (Bruker IFS 66 V/S) allowing for providing measurements in the 0.1-21 THz frequency range made in a vacuum or a nitrogen atmosphere. Their absorption frequencies of the RDX, HMX, 2,4-DNT, and 2,6-DNT are the following: In the spectra of each of these substances, there are the absorption frequencies, which are equal or close to the linear combinations of other ones. Some of these linear combinations contain the doubled absorption frequency. For example, for RDX, one can see the following combinations: ν = 13.86 − 12.33 = 1.53 THz (close to 1.50 THz), ν = 14.52 − 12.33 = 2.19 THz (close to 2.20 THz), and ν = 6.73·2 − 3.08 = 10.38 THz (close to 10.35 THz). For the substance HMX, we see: ν = 1.78·2 + 2.51 = 6.07 THz (close to 6.06 THz). Additionally, for 2,4-DNT, there are two exact linear combinations containing the doubled absorption frequency: ν = 4.98·2 − 1.08 = 8.88 THz, ν = 10.56·2 − 1.08 = 20.04 THz. Note that both linear combinations are formed with the participation of the same absorption frequency ν = 1.08 THz. For 2,6-DNT, there are also two linear combinations with the doubled absorption frequency: ν = 1.56·2 + 2.5 = 5.62 THz (close to 5.61), ν = 11.43·2 − 17.25 = 5.61 THz. The last one is exact.
In [48] an investigation of the 2,4-DNT detection we used the measurements of the absorption frequencies of 2,4-DNT made in transmission mode at relative humidity 4.2% and in the time interval t = [0, 20] ps with time step h t = 0.05 ps. Because we use this substance as example in current paper then, in Figure 10, we show its spectrum in the frequency ranges ν = [0, 3.0] THz (a), [4.6, 9.0] THz (b) made with the spectrum resolution ∆ν = 0.05 THz. The spectral minima at the frequencies ν = 0.4, 0.65, 1.1, 1.4, 2.55, 4.95, and 8.8 THz are in a good agreement with the absorption peaks mentioned above [6]. It should be stressed that, in the 2,4-DNT spectrum (a), (b) one can also see the linear combination of the frequencies containing the doubled absorption frequency, which is close to the linear combination of frequencies mentioned above: ν = 4.95·2 − 1.1 = 8.8 THz. The minimum at the doubled frequency ν = 9.9 THz = 4.95·2 THz is also present in (b). Therefore, these measurements are correct.
In (c-f), an evolution of the spectral amplitude modulus at the frequencies ν = 1.1, 4.95, 8.8, and 4.95·2 = 9.9 THz are depicted, which is computed with the window length T = 2.8 ps and the window shift ∆ = 0.05 ps.
Appendix B demonstrates the computation of the integral spectral intensity |P dyn (ν)| (i.e., the brightness of the spectral line) using the time-dependent spectral brightness. We show that the value of the integral spectral intensity |P dyn (ν)| computed at the chosen frequency ν does not depend on the window length T and its shift ∆, and it is equal to the spectral amplitude |P(ν, t b ) L t | of the Fourier transform of the function F(t) at this frequency in the time interval [t b , t e ] length L t = t e − t b . The second remark refers to choosing of the window shift ∆ = 0.05 ps. This value follows from a requirement of computation of the spectral dynamics at high frequencies ν = 4.95, 8.8, 9.9 THz ≤ ν max = 10 THz. According to the Nyquist-Shannon-Kotel'nikov sampling theorem [64][65][66], the maximal frequency, which can be resolved for the time step h t , is ν max ≤ 1/(2·h t ). At construction of the spectral line dynamics, the window shift ∆ plays the role of the time step h t . Thus, the window shift ∆ must satisfy the condition ∆ ≤ 1/(2·ν max ) = 0.05 ps. For convenience, the influence of window shift ∆ on the accuracy of constructing the spectral line dynamics is investigated in Appendix C.
We see in Figure 10 that the shapes of the time-dependent spectral amplitudes for the frequencies ν = 4.95, 8.8, 9.9 THz are very similar, and the corresponding correlation coefficient are close to unity:  (Figure 10a,b). Taking into account that the additional absorption frequency ν = 8.8 THz is the linear combination of the frequencies ν = 4.95, 1.1 THz (ν = 8.8 = 4.95·2 − 1.1 THz), it might mean that, during the process of the frequency ν = 8.8 THz generation, the part of the pulse energy corresponding to the base frequencies ν = 4.95, 1.1 THz can be transferred to the harmonics at the frequency ν = 8.8 THz. To clarify this, it is necessary to analyze the beginning of the spectral intensity evolution at these frequencies.
| (9.9) | P =9.69·10 −4 (Figure 10a,b). Taking into account that the additional absorption frequency ν = 8.8 THz is the linear combination of the frequencies ν = 4.95, 1.1 THz (ν = 8.8 = 4.95·2 − 1.1 THz), it might mean that, during the process of the frequency ν = 8.8 THz generation, the part of the pulse energy corresponding to the base frequencies ν = 4.95, 1.1 THz can be transferred to the harmonics at the frequency ν = 8.8 THz. To clarify this, it is necessary to analyze the beginning of the spectral intensity evolution at these frequencies.  Therefore, one can see high correlation between the spectral line dynamics at the frequencies ν = 4.95 THz and ν = 4.95·2 = 9.9 THz. Thus, the time delay at the beginning of the generation at the frequency ν = 4.95·2 = 9.9 THz with respect to spectral intensity at the basic frequency ν = 4.95 THz, as well as the growth rate of the spectral amplitude evolution at ν = 4.95·2 = 9.9 THz, which is less  Therefore, one can see high correlation between the spectral line dynamics at the frequencies ν = 4.95 THz and ν = 4.95·2 = 9.9 THz. Thus, the time delay at the beginning of the generation at the frequency ν = 4.95·2 = 9.9 THz with respect to spectral intensity at the basic frequency ν = 4.95 THz, as well as the growth rate of the spectral amplitude evolution at ν = 4.95·2 = 9.9 THz, which is less Therefore, one can see high correlation between the spectral line dynamics at the frequencies ν = 4.95 THz and ν = 4.95·2 = 9.9 THz. Thus, the time delay at the beginning of the generation at the frequency ν = 4.95·2 = 9.9 THz with respect to spectral intensity at the basic frequency ν = 4.95 THz, as well as the growth rate of the spectral amplitude evolution at ν = 4.95·2 = 9.9 THz, which is less than at ν = 4.95 THz, confirm the appearance of a doubled absorption frequency in the 2,4-DNT spectrum due to the frequency conversion process.
Because, in our data of the physical experiments, an influence of humidity is present, we need to discuss this influence on the 2,4 DNT absorption frequencies. For this aim, we compare the absorption frequencies ν = 1.1, 4.95, 8.8, and 9.9 THz of the 2,4-DNT signal measured at 4.2% relative humidity with the corresponding values in the spectra of reference signal measured in the ambient air with increasing (from 4.2% up to 12 %) relative humidity. In Figure 12a-c, the spectrum of the 2,4-DNT (4.2% relative humidity) is depicted together with the reference spectra measured at 4.2% (a-c) and 12% (b,c) relative humidity in the frequency ranges containing the frequencies ν = 1.1 (a), 4.95 (b), 8.8, and 9.9 THz (c). In (a), at relative humidity 4.2%, the water vapor practically does not influence on the 2,4-DNT spectral minimum at the frequency ν = 1.1 THz. In (b,c), one can see the minima at the frequencies ν = 4.95, 8.8, and 9.9 THz in the 4.2% humidity reference spectrum (red line). However, these minima are not preserved in the 12% humidity reference spectrum (blue line). In (d-f), the same spectra are computed at decreased spectral resolution Δν = 0.01 ps. One can see that, for both values of humidity 4.2% (red line) and 12% (blue line), the minima at the 2,4-DNT absorption frequencies ν = 4.95, 8.8 and 9.9 THz are absent in the reference spectra (e,f). In (d) one can see a minimum at the frequency ν = 1.09 THz in the reference spectrum (blue line), but it is weakly pronounced. Thus, an increase of spectral resolution confirms that the minima at the frequencies ν = 1.1, 4.95, 8.8, and 9.9 THz in the 2,4-DNT spectrum do not appear because of the THz radiation absorption by water vapor and they are true 2,4-DNT absorption frequencies.

Results of the Experiment Using a Narrowband CW THz Radiation
In this Section, we investigate the absorption frequencies that appear at the narrowband THz signal transmission through the neutral substances and compare them with the absorption frequencies in the spectra appearing at the broadband THz signal transmission through the similar substances. In (d-f), the same spectra are computed at decreased spectral resolution ∆ν = 0.01 ps. One can see that, for both values of humidity 4.2% (red line) and 12% (blue line), the minima at the 2,4-DNT absorption frequencies ν = 4.95, 8.8 and 9.9 THz are absent in the reference spectra (e,f). In (d) one can see a minimum at the frequency ν = 1.09 THz in the reference spectrum (blue line), but it is weakly pronounced. Thus, an increase of spectral resolution confirms that the minima at the frequencies ν = 1.1, 4.95, 8.8, and 9.9 THz in the 2,4-DNT spectrum do not appear because of the THz radiation absorption by water vapor and they are true 2,4-DNT absorption frequencies.

Results of the Experiment Using a Narrowband CW THz Radiation
In this Section, we investigate the absorption frequencies that appear at the narrowband THz signal transmission through the neutral substances and compare them with the absorption frequencies in the spectra appearing at the broadband THz signal transmission through the similar substances.
The aim of this consideration is the finding of the absorption frequencies that are absent at measurement with CW THz radiation. Figure 1 shows the scheme of providing of the physical experiment. If a minimum at the doubled absorption frequency in the spectrum of the narrowband THz signal transmitted through the sample is present, then this is a true absorption frequency (Figure 1a). However, the up-conversion process induces this frequency if the spectrum minimum at the doubled frequency is absent at the action of the narrowband THz pulse (Figure 1b).
Let us notice that the first sample is the standard aspirin pellet, the second one is a piece of soap, and the third one is a standard printer paper. The fourth example is a chocolate bar and powdered chocolate in a polyethylene bag. The empty polyethylene bag is also investigated.

Aspirin Pellet
In Figure 13a-c, a standard aspirin pellet is shown at different angles to the receiver. In Figure 13d, the pellet is placed in the thin polyethylene bag. The measured narrowband spectra were averaged over a number of measurements (i.e., four) and the resulted spectrum (Aspirin_AVG, for brevity) is depicted in the frequency ranges ν = For finding of influence of the water vapor on the measured spectrum, in Figure 13g,h, the reference spectrum is depicted in the same frequency ranges. It contains the spectral minima at the frequencies close to those of Aspirin_AVG Figure 13e,f, but the most part of them lies at the spectral distances exceeding the spectral resolution. Nevertheless, the Aspirin_AVG spectral minima at the frequencies ν = 0.532 THz and 0.959 THz in Figure 13e can be caused by environment influence, because, in the reference spectrum Figure 13g, there are spectral minima at the frequencies ν = 0.533 THz and 0.959 THz, which are within the spectral resolution.
Let us notice that one of the absorption frequencies of Aspirin, being measured at using the broadband THz pulse (RIKEN), is a linear combination of the basic absorption frequencies, and we name this frequency as ν minus = 0.44 THz = (1.4-0.96) THz. For verification, if this absorption frequency is true or spurious one, it is necessary to show the presence (or absence) of such a minimum in the Aspirin_AVG spectrum, measured using the narrowband signal. We see in Figure 13e that there is a minimum at the Aspirin_AVG frequency ν = 0.441 THz. However, its appearance might be caused by the environment influence, because there is the minimum at the close frequency ν = 0.439 THz in the reference spectrum (Figure 13g). For clearing, the corresponding spectra are shown in Figure 13e,g in magnified view. Therefore, in this case it is difficult to draw an unambiguous conclusion, since there is a reference spectral minimum at the close frequency ν = 0.44-0.001 THz.
reference spectrum is depicted in the same frequency ranges. It contains the spectral minima at the frequencies close to those of Aspirin_AVG Figure 13e,f, but the most part of them lies at the spectral distances exceeding the spectral resolution. Nevertheless, the Aspirin_AVG spectral minima at the frequencies ν = 0.532 THz and 0.959 THz in Figure 13e can be caused by environment influence, because, in the reference spectrum Figure 13g, there are spectral minima at the frequencies ν = 0.533 THz and 0.959 THz, which are within the spectral resolution.   Let us notice that one of the absorption frequencies of Aspirin, being measured at using the broadband THz pulse (RIKEN), is a linear combination of the basic absorption frequencies, and we name this frequency as νminus = 0.44 THz = (1.4-0.96) THz. For verification, if this absorption frequency is true or spurious one, it is necessary to show the presence (or absence) of such a minimum in the Aspirin_AVG spectrum, measured using the narrowband signal. We see in Figure 13e that there is a minimum at the Aspirin_AVG frequency ν = 0.441 THz. However, its appearance might be caused by the environment influence, because there is the minimum at the close frequency ν = 0.439 THz in the reference spectrum (Figure 13g). For clearing, the corresponding spectra are shown in Figure  13e,g in magnified view. Therefore, in this case it is difficult to draw an unambiguous conclusion, since there is a reference spectral minimum at the close frequency ν = 0.44-0.001 THz.

Piece of Soap
The second example is a piece of soap placed at different distances and angles to the receiver, see Figure 14a-e. In (e), the sample has an inhomogeneous surface. We note that, in [52], we investigated the absorption frequencies of soap (Soap_BB for brevity, BB means "broadband") in the transmission mode using the broadband THz pulse. Soap_BB spectrum possesses the following absorption frequencies in the frequency range ν < 1.6 THz: ν = 0.25, 0.35, 0.6, 0.7, 0.8, 0.9, 1.15 THz, as seen in Figure 14f. There are also several spectral minima in the frequency range ν > 1.6 THz: ν = 1.7, 3.05, 3.35, 3.65, 3.95 THz (f), (g).

Piece of Soap
The second example is a piece of soap placed at different distances and angles to the receiver, see Figure 14a-e. In (e), the sample has an inhomogeneous surface. We note that, in [52], we investigated the absorption frequencies of soap (Soap_BB for brevity, BB means "broadband") in the transmission mode using the broadband THz pulse. Soap_BB spectrum possesses the following absorption frequencies in the frequency range ν < 1.6 THz: ν = 0.25, 0.35, 0.6, 0.7, 0.8, 0.9, 1.15 THz, as seen in Figure 14f. There are also several spectral minima in the frequency range ν > 1.6 THz: ν = 1.7, 3.05, 3.35, 3.65, 3.95 THz (f), (g).
The second example is a piece of soap placed at different distances and angles to the receiver, see Figure 14a-e. In (e), the sample has an inhomogeneous surface. We note that, in [52], we investigated the absorption frequencies of soap (Soap_BB for brevity, BB means "broadband") in the transmission mode using the broadband THz pulse. Soap_BB spectrum possesses the following absorption frequencies in the frequency range ν < 1.   (Figure 15a,b). It is important that, in the reference spectrum, the minima at these frequencies are also absent (Figure 15c,d), and the closest spectral minima are: ν = 0.6 ± 0.003 THz (c), 0.9 − 0.002 THz (d). It is far enough, keeping in mind the spectral resolution. Therefore, the appearance of the linear combinations of frequencies νminus = 0.6, 0.9 THz in the spectrum of the broadband pulse Soap_BB is induced by the frequency conversion process. Thus, these frequencies are not the true Soap_BB absorption frequencies.  The measured spectra using the narrowband THz signals transmitted through the soap samples were averaged over a number of measurements (i.e., five), and the resulting spectrum (Soap_AVG, for brevity) is depicted in the frequency ranges ν =  (Figure 15a,b). It is important that, in the reference spectrum, the minima at these frequencies are also absent (Figure 15c,d), and the closest spectral minima are: ν = 0.6 ± 0.003 THz (c), 0.9 − 0.002 THz (d). It is far enough, keeping in mind the spectral resolution. Therefore, the appearance of the linear combinations of frequencies ν minus = 0.6, 0.9 THz in the spectrum of the broadband pulse Soap_BB is induced by the frequency conversion process. Thus, these frequencies are not the true Soap_BB absorption frequencies.
Sensors 2020, 20, 1859 20 of 36 narrowband spectrum, the minimum at the frequencies ν = 0.6 ± 0.001 THz and ν = 0.9 ± 0.001 THz are absent (Figure 15a,b). It is important that, in the reference spectrum, the minima at these frequencies are also absent (Figure 15c,d), and the closest spectral minima are: ν = 0.6 ± 0.003 THz (c), 0.9 − 0.002 THz (d). It is far enough, keeping in mind the spectral resolution. Therefore, the appearance of the linear combinations of frequencies νminus = 0.6, 0.9 THz in the spectrum of the  It should be noted that one Soap_BB minimum at the frequency νminus = 0.6 THz corresponds to three different linear combinations of the basic absorption frequencies. Therefore, in Figure 16  It should be noted that one Soap_BB minimum at the frequency ν minus = 0.6 THz corresponds to three different linear combinations of the basic absorption frequencies. Therefore, in Figure 16 In our opinion, the appearance of the Soap_BB spectral minima at the frequencies ν = 3.35, 3.65, and 3.95 THz (Figure 14g) can be due to the energy levels transitions that are caused by the photons with the frequencies ν = 3.05 THz and ν = 0.3, 0.6 THz (see the spectral dynamics in (d)), because there is the high correlation between the corresponding spectral line dynamics. The physical mechanism of these energy level transitions is the cascading mechanism that is described in [57,58], which occurs under the action of the broadband THz pulse. We analyzed such a mechanism of the energy level excitation in the framework of one-dimensional (1D) Maxwell's equations describing the pulse propagation and matrix density formalism for the description of the response of a medium. However, the direct experimental verification of our supposing might be made only using the narrowband CW radiation.

Sheet of Paper
The third example is a standard sheet of paper. One, two, or four sheets are placed at different distances and angles with respect to the receiver (see Figure 17a-e). We stress that, in [53], the absorption frequencies of the paper sheet (Paper_BB, for brevity) were found out in transmission mode using the broadband THz pulse. They are depicted in the Figure 18a  In our opinion, the appearance of the Soap_BB spectral minima at the frequencies ν = 3.35, 3.65, and 3.95 THz (Figure 14g) can be due to the energy levels transitions that are caused by the photons with the frequencies ν = 3.05 THz and ν = 0.3, 0.6 THz (see the spectral dynamics in (d)), because there is the high correlation between the corresponding spectral line dynamics. The physical mechanism of these energy level transitions is the cascading mechanism that is described in [57,58], which occurs under the action of the broadband THz pulse. We analyzed such a mechanism of the energy level excitation in the framework of one-dimensional (1D) Maxwell's equations describing the pulse propagation and matrix density formalism for the description of the response of a medium. However, the direct experimental verification of our supposing might be made only using the narrowband CW radiation.

Sheet of Paper
The third example is a standard sheet of paper. One, two, or four sheets are placed at different distances and angles with respect to the receiver (see Figure 17a-e). We stress that, in [53], the absorption frequencies of the paper sheet (Paper_BB, for brevity) were found out in transmission mode using the broadband THz pulse. They are depicted in the Figure 18a,b for the frequency range ν < 1.6 THz ν = 0. 36 us for comparison of those absorption frequencies with the absorption frequencies measured while using the narrowband THz CW radiation.

Sheet of Paper
The third example is a standard sheet of paper. One, two, or four sheets are placed at different distances and angles with respect to the receiver (see Figure 17a-e). We stress that, in [53], the absorption frequencies of the paper sheet (Paper_BB, for brevity) were found out in transmission mode using the broadband THz pulse. They are depicted in the Figure 18a,b for the frequency range ν < 1.6 THz ν = 0. 36 Figure 17. Location of the paper sheets in the THz spectrometer (a-e) . One (a,b), two (c,d) and four sheets (e) placed at different distances and angles with respect to the receiver.
The spectra of the narrowband THz signals that were transmitted through the paper samples were averaged over the number of measurements (i.e., six) and the resulted spectrum (Paper_AVG for brevity) is depicted in Figure 18a for the frequency range ν = [0.35, 1.5] THz. It possesses the spectral minima at the frequencies ν = 0.362, 0.56, 0.76, 0.959, 1.199, and 1.477 THz. One can see that the Paper_AVG spectral minima and the Paper_BB ones are in a good agreement. At the same time, Figure 17. Location of the paper sheets in the THz spectrometer (a-e). One (a,b), two (c,d) and four sheets (e) placed at different distances and angles with respect to the receiver.
The spectra of the narrowband THz signals that were transmitted through the paper samples were averaged over the number of measurements (i.e., six) and the resulted spectrum (Paper_AVG for brevity) is depicted in Figure 18a Figure 19). The minima, which are closest to these values, are located at the frequencies ν = 0.36 + 0.002 THz (a) and ν = 1.48 ± 0.003 THz (c). It is very important that the minimum at the frequencies ν = 0.36 ± 0.001 THz and ν = 1.48 ± 0.001 THz are also absent in the reference spectrum (Figure 19b,d) Figure 19). The minima, which are closest to these values, are located at the frequencies ν = 0.36 + 0.002 THz (a) and ν = 1.48 ± 0.003 THz (c). It is very important that the minimum at the frequencies ν = 0.36 ± 0.001 THz and ν = 1.48 ± 0.001 THz are also absent in the reference spectrum (Figure 19b,d).

Polyethylene Bag
Below, we investigate the spectral properties of the empty polyethylene bag (Figure 20a) using a narrowband THz pulse (Bag signal, for brevity). The spectral properties of polyethylene are wellknown. Accordingly, with accordance to RIKEN THz database [67], the polyethylene (PE) does not demonstrate any pronounced absorption frequencies in the frequency range ν < 3.0 THz, but, in the frequency range ν > 3.0 THz, it possesses two frequencies: ν = 6.72, 7.88 THz. In [59], we investigated the broadband THz pulse at non-zero humidity. In Figure 20, the corresponding spectrum of the first sub-pulse of the THz signal reflected from a pure polyethylene (PE_Thick, for brevity) is shown in the frequency ranges ν = [0, 2.0] THz (b) and [6.4, 8.0] THz (c). In (d), (e), the reference spectrum is depicted in the same frequency ranges. These measurements were made with the spectral resolution Δν = 0.04 THz. Therefore, the frequency conversion process causes the appearance of the frequencies ν minus = 0.36 THz and 1.48 THz, which are linear combinations of absorption frequencies, in the spectrum of the broadband Paper_BB signal.

Polyethylene Bag
Below, we investigate the spectral properties of the empty polyethylene bag (Figure 20a) using a narrowband THz pulse (Bag signal, for brevity). The spectral properties of polyethylene are well-known. Accordingly, with accordance to RIKEN THz database [67], the polyethylene (PE) does not demonstrate any pronounced absorption frequencies in the frequency range ν < 3.0 THz, but, in the frequency range ν > 3.0 THz, it possesses two frequencies: ν = 6.72, 7.88 THz. In [59], we investigated the broadband THz pulse at non-zero humidity. In Figure 20, the corresponding spectrum of the first sub-pulse of the THz signal reflected from a pure polyethylene (PE_Thick, for brevity) is shown in the frequency ranges ν = [0, 2.0] THz (b) and [6.4, 8.0] THz (c). In (d), (e), the reference spectrum is depicted in the same frequency ranges. These measurements were made with the spectral resolution ∆ν = 0.04 THz.
When comparing PE_Thick and reference spectra in the range ν = [0, 2.0] THz (b), (d), one can see that the part of spectral minima is caused by the environment or noise influence. Only two spectral minima at the frequencies ν = 0.52, 1.16 THz (b) are free from this influence, because the corresponding spectral minima in (d) are absent. In the frequency range ν = [6.4, 8.0] THz (c) one can see two minima ν = 6.72, 7.88 THz, which are equal to the PE absorption frequencies (RIKEN). Let us note that the minima at these frequencies in the reference spectrum (e) are also absent. Once more minimum in the PE_Thick spectrum, which is not caused by the environment or noise, is located at the frequency ν = 7.24 THz (c), (e). Additionally, the spectral minima in (b) at the frequencies ν = 0.52, 1.16 THz are the caused by the frequency conversion process for appearance of the broadband PE_Thick absorption frequencies: ν minus = (7.24 − 6.72) =0.52 THz; and, ν minus = (7.88 − 6.72) = 1.16 THz.
For clarity, in (f), (h), the spectrum of the narrowband Bag signal is depicted in the frequency ranges ν = [0.5, 0.54] THz (f), [1.156, 1.164] THz (h). The spectral minima at the frequencies ν = 0.52 THz and 1.16 THz are absent in (f), (h), which refers to the narrowband Bag signal spectrum. Taking the spectral resolution ∆ν = 0.001 THz and the absence of the corresponding minima in the reference spectrum (g), (i) into account, we conclude that an appearance of these frequencies ν minus = 0.52 THz and 1.16 THz in the spectrum of the broadband PE_Thick pulse (Figure 20b,c) is caused by the frequency conversion process. Therefore, these frequencies are spurious ones in the PE_Thick pulse spectrum.
corresponding spectral minima in (d) are absent. In the frequency range ν = [6.4, 8.0] THz (c) one can see two minima ν = 6.72, 7.88 THz, which are equal to the PE absorption frequencies (RIKEN). Let us note that the minima at these frequencies in the reference spectrum (e) are also absent. Once more minimum in the PE_Thick spectrum, which is not caused by the environment or noise, is located at the frequency ν = 7.24 THz (c), (e). Additionally, the spectral minima in (b) at the frequencies ν = 0.52, 1.16 THz are the caused by the frequency conversion process for appearance of the broadband PE_Thick absorption frequencies: νminus = (7.24 − 6.72) =0.52 THz; and, νminus = (7.88 − 6.72) = 1. 16   For clarity, in (f), (h), the spectrum of the narrowband Bag signal is depicted in the frequency ranges ν = [0.5, 0.54] THz (f), [1.156, 1.164] THz (h). The spectral minima at the frequencies ν = 0.52 THz and 1.16 THz are absent in (f), (h), which refers to the narrowband Bag signal spectrum. Taking the spectral resolution Δν = 0.001 THz and the absence of the corresponding minima in the reference spectrum (g), (i) into account, we conclude that an appearance of these frequencies νminus = 0.52 THz and 1.16 THz in the spectrum of the broadband PE_Thick pulse (Figure 20b,c) is caused by the frequency conversion process. Therefore, these frequencies are spurious ones in the PE_Thick pulse spectrum.

Chocolate Bar
In the last example, we consider the samples with chocolate bar and powdered chocolate in the polyethylene bag (see Figure 21a-d). While using the broadband THz pulse in [46], we investigated the absorption frequencies of chocolate (Choc_BB, for brevity) measured in transmission mode. The Choc_BB pulse spectrum possesses the following pronounced spectral minima in the frequency range ν < 1.  (RIKEN [67]). Therefore, these spectral minima are due to the presence of sugar in the sample. ν < 1.6 THz: ν = 0.44, 0.56, 0.76, 1.0, 1.16, 1.32, and 1.44 THz (Figure 22a) and several spectral minima in the frequency range ν>1.6 THz: ν = 1.8, 1.96, 2.4, 2.6, 2.84, 2.96, 3.12, and 3.48 THz (b). The spectral resolution in (a), (b) is equal to Δν = 0.04 THz. In those experiments, the obtained spectral minima at the frequencies ν = 0.56, 0.76, 1.0 THz appear due to water vapor absorption. A part of the spectral minima ν = 1.44, 1.8, 2.6, 2.84 THz (a), (b) are close to the absorption frequencies of pure Sucrose: ν = 1.43, 1.82, 2.55, 2.8 THz (RIKEN [67]). Therefore, these spectral minima are due to the presence of sugar in the sample. The spectra of the narrowband THz signals were averaged over the measurements (two for both the chocolate bar and the powdered chocolate). We will call these signals Choc_Bar_AVG and Choc_Bag_AVG, for brevity. In Figure 22  The spectra of the narrowband THz signals were averaged over the measurements (two for both the chocolate bar and the powdered chocolate). We will call these signals Choc_Bar_AVG and Choc_Bag_AVG, for brevity. In Figure 22  The Choc_Bar_AVG signal and Choc_Bag_AVG signal spectra demonstrate the spectral minima at the frequencies ν = 1.321, 1.439 THz (c), and ν = 1.322, 1.437 THz (e), which are in a good agreement with the spectral minima of the broadband Choc_BB signal (a). However, in the reference spectrum (g), one can see the spectral minima at the close or equal frequencies ν = 1.322, 1.437 THz. Taking the spectral resolution ∆ν = 0.001 THz into account means that, in both spectra (c), (e) the spectral minima at the frequencies ν = 1.321 THz (c) (1.322 THz (e)) is due the atmospheric air influence. The same situation occurs at the spectral minimum ν = 1.437 THz in the Choc_Bag_AVG signal spectrum.

Chocolate Bar
In the last example, we consider the samples with chocolate bar and powdered chocolate in the polyethylene bag (see Figure 21a-d). While using the broadband THz pulse in [46], we investigated the absorption frequencies of chocolate (Choc_BB, for brevity) measured in transmission mode. The Choc_BB pulse spectrum possesses the following pronounced spectral minima in the frequency range ν < 1.6 THz: ν = 0.44, 0.56, 0.76, 1.0, 1.16, 1.32, and 1.44 THz (Figure 22a) and several spectral minima in the frequency range ν>1.6 THz: ν = 1.8, 1.96, 2.4, 2.6, 2.84, 2.96, 3.12, and 3.48 THz (b). The spectral resolution in (a), (b) is equal to Δν = 0.04 THz. In those experiments, the obtained spectral minima at the frequencies ν = 0.56, 0.76, 1.0 THz appear due to water vapor absorption. A part of the spectral minima ν = 1.44, 1.8, 2.6, 2.84 THz (a), (b) are close to the absorption frequencies of pure Sucrose: ν = 1.43, 1.82, 2.55, 2.8 THz (RIKEN [67]). Therefore, these spectral minima are due to the presence of sugar in the sample. The spectra of the narrowband THz signals were averaged over the measurements (two for both the chocolate bar and the powdered chocolate). We will call these signals Choc_Bar_AVG and Choc_Bag_AVG, for brevity. In Figure 22    The Choc_Bar_AVG signal and Choc_Bag_AVG signal spectra demonstrate the spectral minima at the frequencies ν = 1.321, 1.439 THz (c), and ν = 1.322, 1.437 THz (e), which are in a good agreement with the spectral minima of the broadband Choc_BB signal (a). However, in the reference spectrum (g), one can see the spectral minima at the close or equal frequencies ν = 1.322, 1.437 THz. Taking the spectral resolution Δν = 0.001 THz into account means that, in both spectra (c), (e) the spectral minima at the frequencies ν = 1.321 THz (c) (1.322 THz (e)) is due the atmospheric air influence. The same situation occurs at the spectral minimum ν = 1.437 THz in the Choc_Bag_AVG signal spectrum.
In Figure 22i, the spectrum of the Bag narrowband THz signal is shown in the frequency range ν = [1.2, 1.5] THz. One can see two spectral minima at the frequencies ν = 1.321, 1.436 THz, which are also close (within the spectral resolution) to the Ref_Choc_Bar spectral minima at the frequencies ν= 1.322, 1.437 THz (g). Taking the absence of polyethylene bag absorption frequencies in the frequency range ν < 3.0 THz (RIKEN database [67]) into account, the appearance of these two minima in the signal Choc_Bag_AVG spectrum can also be caused by the influence of the atmospheric air. That is, the spectral minimum in the signal Choc_Bar_AVG spectrum ν= 1.439 THz (c) can be regarded as an absorption frequency of Chocolate.
In the broadband Choc_BB spectrum Figure 22 (a), (b), one can see the absorption frequencies, which is caused by the frequency conversion process: νminus = (2.96 − 1.8) = (3.12 − 1.96) = 1.16 THz. At the same time, in the Choc_Bar_AVG and Choc_Bag_AVG narrowband spectra this minimum is absent (Figure 22d,f). Moreover, in the reference spectrum (Figure 22h), the minimum at this frequency is also absent. Therefore, this absorption frequency νminus = 1.16 THz is a spurious one observing in the spectrum of the broadband Choc_BB signal that was measured in the transmission mode.
Other possible spurious absorption frequency, including a doubled absorption frequency, is the following: νminus = (1.96·2 − 3.48) = (2.4 − 1.96) = 0.44 THz (Figure 22a,b). The corresponding spectral minimum at this frequency is absent in the spectrum of the narrowband Choc_Bar_AVG signal and in the reference spectrum (Figure 23a,c). Note that water vapor absorption is absent at the frequency ν = 0.44 THz [6,67]. Thus, the frequency νminus = 0.44 THz is not a true absorption frequency that belongs to the broadband Choc_BB spectrum. At the same time, in the Choc_Bag_AVG and Bag spectra (b), (d) there is a spectral minimum at the frequency ν = 0.441 THz = 0.44 + 0.001 THz. Taking the spectral resolution Δν = 0.001 THz into account, this minimum in the Choc_Bag_AVG signal spectrum can be caused by the noise influence in this experiment. Therefore, we believe that a conclusion regarding the spurious absorption frequency is right. In Figure 22i, the spectrum of the Bag narrowband THz signal is shown in the frequency range ν = [1.2, 1.5] THz. One can see two spectral minima at the frequencies ν = 1.321, 1.436 THz, which are also close (within the spectral resolution) to the Ref_Choc_Bar spectral minima at the frequencies ν= 1.322, 1.437 THz (g). Taking the absence of polyethylene bag absorption frequencies in the frequency range ν < 3.0 THz (RIKEN database [67]) into account, the appearance of these two minima in the signal Choc_Bag_AVG spectrum can also be caused by the influence of the atmospheric air. That is, the spectral minimum in the signal Choc_Bar_AVG spectrum ν = 1.439 THz (c) can be regarded as an absorption frequency of Chocolate.
In the broadband Choc_BB spectrum Figure 22 (a), (b), one can see the absorption frequencies, which is caused by the frequency conversion process: ν minus = (2.96 − 1.8) = (3.12 − 1.96) = 1.16 THz. At the same time, in the Choc_Bar_AVG and Choc_Bag_AVG narrowband spectra this minimum is absent (Figure 22d,f). Moreover, in the reference spectrum (Figure 22h), the minimum at this frequency is also absent. Therefore, this absorption frequency ν minus = 1.16 THz is a spurious one observing in the spectrum of the broadband Choc_BB signal that was measured in the transmission mode.
Other possible spurious absorption frequency, including a doubled absorption frequency, is the following: ν minus = (1.96·2 − 3.48) = (2.4 − 1.96) = 0.44 THz (Figure 22a,b). The corresponding spectral minimum at this frequency is absent in the spectrum of the narrowband Choc_Bar_AVG signal and in the reference spectrum (Figure 23a,c). Note that water vapor absorption is absent at the frequency ν = 0.44 THz [6,67]. Thus, the frequency ν minus = 0.44 THz is not a true absorption frequency that belongs to the broadband Choc_BB spectrum. At the same time, in the Choc_Bag_AVG and Bag spectra (b), (d) there is a spectral minimum at the frequency ν = 0.441 THz = 0.44 + 0.001 THz. Taking the spectral resolution ∆ν = 0.001 THz into account, this minimum in the Choc_Bag_AVG signal spectrum can be caused by the noise influence in this experiment. Therefore, we believe that a conclusion regarding the spurious absorption frequency is right.

Conclusions
One of the physical mechanisms of the spurious absorption frequencies appearance caused by the frequency conversion process (the doubled frequency, sum frequency generation, and difference frequency generation) because of the nonlinear response of a medium at the broadband THz pulse interaction with a medium is predicted while using the computer simulation results. The results of a few physical experiments provided with the narrowband THz CW radiation and broadband THz pulses transmitted through the ordinary and dangerous materials confirm this phenomenon. They are more pronounced in the frequency range not belonging to the incident pulse spectrum.
Based on computer simulation results, we show that, if the pulse consists of the intersecting subpulses with substantially overlapping spectra, then it is possible to observe the shift of the spurious absorption frequency. At the same time, the dynamics of the spectral lines at the basic frequency and the shifted spurious absorption frequency are shown to be very similar. The correlation between them ranges from 0.88 to 0.98. Let us notice that such a structure of the THz pulse corresponds to its transmission through a layer of a medium (tablet, for example). If the pulse possesses the sub-pulse structure, which corresponds to consideration of the main THz pulse only, then the frequency shift at the conversion process is not observed. The same processes also occur for the emission frequency conversion.
The spurious absorption (emission) frequency appearance might be used as an additional effective tool for the detection and identification of the substance based on the integral correlation criteria due to the high degree of similarity between the spectral line dynamics at the basic and converted frequencies. This method can be especially effective if noise or environment influence mask the basic absorption frequency.
It is desirable to conduct a physical experiment under laboratory conditions, where there is no influence of water vapor and other factors, in order to illustrate the presence of spurious frequencies in measured THz signals and as possible future direction of the research. In this article, we used the results of the physical experiments conducted under real conditions to emphasize the presence of the discussed effect in this case, which are most important for a practice, and the possibility of its use for the substance detection and identification. We also note that the effect of the carrier frequency shift, which occurs if the THz signal consists of several pulses with overlapping, must be taken into account in both the diagnostics of the substances and in the construction of the expert systems based on artificial intelligence.
The approach that was developed by us (analysis of the spectral line dynamics) can be used under real conditions together with other THz sensors platforms (for example, [68][69][70][71]) used for the identification and detection of the substance. It might essentially increase a probability of the substance detection. Let us remind that, in contrast to the standard THz TDS method, the analysis of the spectral line dynamics allows for us to identify the substance using a noisy THz signal or in the case of inhomogeneous sample surface, or at an analysis of the mixture of substances with similar spectral properties.

Conclusions
One of the physical mechanisms of the spurious absorption frequencies appearance caused by the frequency conversion process (the doubled frequency, sum frequency generation, and difference frequency generation) because of the nonlinear response of a medium at the broadband THz pulse interaction with a medium is predicted while using the computer simulation results. The results of a few physical experiments provided with the narrowband THz CW radiation and broadband THz pulses transmitted through the ordinary and dangerous materials confirm this phenomenon. They are more pronounced in the frequency range not belonging to the incident pulse spectrum.
Based on computer simulation results, we show that, if the pulse consists of the intersecting sub-pulses with substantially overlapping spectra, then it is possible to observe the shift of the spurious absorption frequency. At the same time, the dynamics of the spectral lines at the basic frequency and the shifted spurious absorption frequency are shown to be very similar. The correlation between them ranges from 0.88 to 0.98. Let us notice that such a structure of the THz pulse corresponds to its transmission through a layer of a medium (tablet, for example). If the pulse possesses the sub-pulse structure, which corresponds to consideration of the main THz pulse only, then the frequency shift at the conversion process is not observed. The same processes also occur for the emission frequency conversion.
The spurious absorption (emission) frequency appearance might be used as an additional effective tool for the detection and identification of the substance based on the integral correlation criteria due to the high degree of similarity between the spectral line dynamics at the basic and converted frequencies. This method can be especially effective if noise or environment influence mask the basic absorption frequency.
It is desirable to conduct a physical experiment under laboratory conditions, where there is no influence of water vapor and other factors, in order to illustrate the presence of spurious frequencies in measured THz signals and as possible future direction of the research. In this article, we used the results of the physical experiments conducted under real conditions to emphasize the presence of the discussed effect in this case, which are most important for a practice, and the possibility of its use for the substance detection and identification. We also note that the effect of the carrier frequency shift, which occurs if the THz signal consists of several pulses with overlapping, must be taken into account in both the diagnostics of the substances and in the construction of the expert systems based on artificial intelligence.
The approach that was developed by us (analysis of the spectral line dynamics) can be used under real conditions together with other THz sensors platforms (for example, [68][69][70][71]) used for the identification and detection of the substance. It might essentially increase a probability of the substance detection. Let us remind that, in contrast to the standard THz TDS method, the analysis of the spectral line dynamics allows for us to identify the substance using a noisy THz signal or in the case of inhomogeneous sample surface, or at an analysis of the mixture of substances with similar spectral properties. Acknowledgments: We thank Joseph W. Haus from the University of Dayton, Dayton, USA, for the opportunity to conduct the experiment with a narrowband THz pulse and fruitful discussion.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Fourier-Gabor Window Transform
Below we briefly remind the main notations and formulas at the construction of the evolution of the spectral amplitude (dynamics of spectral line) of the signal F(t) at the chosen frequency ν [44][45][46].
Let F(t) be the transmitted or reflected THz signal, which is measured in the time interval [t b , t e ] length L t = t e − t b , where t b and t e denote its beginning and the end. In order to get information about the evolution of the spectral amplitude at the frequency ν in the time interval [t b , t e ], we use the time window length T ≤ L t , which moves along the signal. At each step, the time window is shifted on the chosen time interval ∆ (window shift). To avoid the spectrum "spreading", in each window length L t the signal F(t) is multiplied by window function g(t), which quickly tends to zero at the left and right ends of the window. Thus, the Fourier-Gabor transform is carried out in each time window: where t j is a time of window beginning, j is a serial number of window, ν is a frequency of interest. The units of t j , T, ∆ and ν are measured in ps and THz, respectively. In order to align the beginning of the physical THz pulse and its representation in the evolution of the spectral amplitude under construction, the spectral amplitude value |P ν (t j )| in the j-th time window is associated with the end of this window: To start the computations, the first window is placed in time interval [t b − T, t b ] and the function F(t) is continued by zeros to the left to point t b − T. The last window is placed in the time interval [t e − T, t e ]. The integrals (A1) are computed using trapezoid method with time step h t = 0.001 ps.

Appendix B Computation of the Integral Spectral Intensity Using the Time-Dependent Spectral Brightness
As it is well-known, the spectral intensity at the chosen frequency ν is equal to quadratic power from the Fourier harmonic of the function F(t) at this frequency in this time interval. However, we use the Fourier-Gabor transform at a computation of the spectral line dynamics with a sliding window. Obviously, parts of the signal can be involved many times at its computation. Therefore, there is a question about relation between the spectral intensity of the Fourier transform and the spectral intensity determined from the Fourier-Gabor transform and how does this spectral intensity get. Below we analyze this problem. We show that the spectral intensity does not depend on the window length T and its shift ∆, and it is equal to the spectral intensity of the Fourier harmonic computed on the time interval [t b , t e ] under analysis. The result, obtained in this Appendix, is used in Section 4 to estimate the contribution of the pulse energy corresponding to the doubled frequency harmonic in the process of the frequency conversion.
Obviously, if we choose in (A1) the following values of window parameters: t j = t b , T = L t and the window function , we get the Fourier transform of the signal F(t) in the time interval [t b , t e ] at the frequency ν: and, therefore, the spectral amplitude at the frequency ν is equal to |P(ν, Further, for simplicity, we assume that t b = 0 and the window function g(t) has the form Then from expression (A3) we have: Let define N W = T/∆ -the number of shifts ∆ in the time window length T, and N D is a number of the window shifting along the signal F(t) from the time point t b = 0 till the time moment t e = L t with window shift ∆ to cover the time interval [0, L t ]. Thus, one can write L t = ∆·(N D − 1). Then t j = −T + ( j − 1)·∆, j = 1, . . . , N D . Therefore, the integral value P(ν, t j ) defined in (A1) can be represented as a sum of integrals over the interval [t j , t j + T]: From expression (A5) we obtain the following set of linear equations with respect to variables X n : The schematic representation of the window moving along the signal F(t) is shown in Figure A1, where the window length and the window shift is equal to T = 10 ps and ∆ = 5 ps, respectively. A time interval is equal to L t = [t b , t e ] = [0,20] = 20 ps. A number of the window shifts ∆ covering the time window is equal to N W = T/∆ = 2. The parameter N D is equal to 5. The set of equations (A6) has the following form: ………………………………………………..

j ND =
The matrix of the system (A9) is lower-triangular one with units on the main diagonal. Accordingly, the set of equations has a unique solution: ……………………………………………….. The set of equations (A6) has the following form: The system (A8) has N D equations with N D + N W − 1 variables. As we mentioned above, the THz signal F(t) is continued by zeros to the left till the time moment t 1 = −T in order to start the computation of the spectral amplitude evolution. Thus, we can assume that the following equations X 1 = 0, X 2 = 0, . . . , X N W −1 = 0 are valid. Let us note if the value F(0) 0, then we put the value X N W = F(0)/2. Thus, a number of the variables in the set of equations (A8) becomes equal to the number of the equations, and they take the form: The matrix of the system (A9) is lower-triangular one with units on the main diagonal. Accordingly, the set of equations has a unique solution: (A10) Our next step consists in finding the sum of variables X 1 , X 2 , . . . ,X N W , X N W +1 , X N W +2 , . . . , X N D +N W −1 , which gives contribution to the integral spectral amplitude |P dyn (ν)| at the chosen frequency ν in the time interval [0, L t ]. It can be found directly using expressions (A10). Indeed, taking into account that X 1 = 0, X 2 = 0, . . . , X N W −1 = 0 and X N W does not contribute to the integral spectral amplitude |P dyn (ν)| because it is computed in the time interval [−∆, 0] [0, L t ] (in Figure A1 this is denoted as X 2 ), we get: Further, the expression (A11) can be transformed to the simpler one using the expressions (A7), (A3): Thus, the integral spectral amplitude can be found from (A12): It is easy to see that the same result occurs if t b 0: We should emphasize an important conclusion following from expressions (A13), (A14). Namely, the integral spectral amplitude |P dyn (ν)|, which is found from the set of equations (A9), does not depend on the window length T and its shift ∆ and it is equal to the spectral amplitude |P(ν, t b ) L t | of the Fourier transform of the function F(t) at the frequency ν in the time To illustrate the obtained result let us consider the artificial signal with one minimum in its spectrum obtained using the functions E base (t) (3) and G 1 (ν) (4) with the following parameters: k = 1, B 1 = 0.8, ν 1 =0.9, ∆ν = 1.0 in the time interval t = [40,60] which length is equal to L t = 20 units (Example 4). Its Fourier spectrum is depicted in Table A1. The corresponding spectral amplitude modulus |P(ν)| L t =20 at the frequencies ν = 0.7, 0.9, 1.05 are shown in Table A1 in the second line. The integral spectral amplitude |P dyn (ν)| at the chosen frequency ν computed using Equation (A11) for window length T = 2.8 and shift ∆ = 0.2 and h t = 0.05; 0.01 is shown in Table A1. In both cases, the computational error in calculating the values of |P dyn (ν)| is about O(h 2 t ). We see good agreement between both values of the spectral amplitudes.  Figure B1 this is denoted as 2 X ), we get: Further, the expression (A11) can be transformed to the simpler one using the expressions (A7), (A3): Thus, the integral spectral amplitude can be found from (A12): It is easy to see that the same result occurs if We should emphasize an important conclusion following from expressions (A13), (A14).
Namely, the integral spectral amplitude | ( )| dyn P ν , which is found from the set of equations (A9), does not depend on the window length T and its shift Δ and it is equal to the spectral amplitude | ( , ) |  Table A1 in the second line.

Appendix C Influence of Window Shift on the Accuracy of Constructing Spectral Line Dynamics
Using as example the 2,4-DNT substance, we illustrate below an influence of the window shift ∆ on the accuracy of constructing spectral line dynamics at frequencies ν = 1.1, 4.95, 8.8 and 9.9 THz because it is very important question for the problem under consideration. In Figure A2  Using as example the 2,4-DNT substance, we illustrate below an influence of the window shift Δ on the accuracy of constructing spectral line dynamics at frequencies ν = 1.1, 4.95, 8.8 and 9.9 THz because it is very important question for the problem under consideration. In Figure A2  The shapes of all the dynamics (a-e) are very similar. However, in (a-c) they contain the oscillations at the frequency ν = 5 THz (close to ν = 4.95 THz) due to the small shift of window Δ ≤ 0.2 ps. If the window shift increases up to value Δ = 0.2, 0.4 ps, then these oscillations disappear. This fact is illustrated by the Figure A3(a), where the magnified view of the spectral line dynamics computed at the frequency ν = 1.1 THz at the increasing window shift Δ = 0.1, 0.2, 0.4 ps is depicted. Let us note that the spectral amplitude corresponding to the frequency ν = 5 THz decreases by 10 times at increasing window shift from Δ = 0.05 ps to 0.4 ps. For illustration of this, the corresponding part of the spectra is shown in Figure A3 Figure A3 (b). In order to estimate the frequency values of these oscillations, we apply the Fourier transform to the evolution of spectral amplitudes modulus | ( )| j P t ν (they are computed using Equations (A1) and (A2)). The spectrum of the corresponding spectral line dynamics computed at Δ=0.005 ps is shown in Figure A3 (e) in the frequency range ν= [8.7, 10.1] THz and it contains a number of spectral maxima corresponding to small oscillations of the time-dependent spectral intensity. Moreover, one can see the spectral maxima at the frequencies ν = 8.82 THz, 9.87 THz, which are close to 2,4-DNT absorption frequencies: ν = 8.8 THz and 9.9 THz. It should be stressed that the spectral resolution in (c-e) is equal to Δν = 0.01 THz.  The shapes of all the dynamics (a-e) are very similar. However, in (a-c) they contain the oscillations at the frequency ν = 5 THz (close to ν = 4.95 THz) due to the small shift of window ∆ ≤ 0.2 ps. If the window shift increases up to value ∆ = 0.2, 0.4 ps, then these oscillations disappear. This fact is illustrated by the Figure A3(a), where the magnified view of the spectral line dynamics computed at the frequency ν = 1.1 THz at the increasing window shift ∆ = 0.1, 0.2, 0.4 ps is depicted. Let us note that the spectral amplitude corresponding to the frequency ν = 5 THz decreases by 10 times at increasing window shift from ∆ = 0.05 ps to 0.4 ps. For illustration of this, the corresponding part of the spectra is shown in Figure A3 At the same time, if the window shift is decreased up to ∆ = 0.05, 0.005 ps, then one can observe the small oscillations occurring at the higher frequencies in the spectral line dynamics. The magnified view of them is depicted in Figure A3 (b). In order to estimate the frequency values of these oscillations, we apply the Fourier transform to the evolution of spectral amplitudes modulus |P ν (t j )| (they are computed using Equations (A1) and (A2)). The spectrum of the corresponding spectral line dynamics computed at ∆ = 0.005 ps is shown in Figure A3 (e) in the frequency range ν = [8.7, 10.1] THz and it contains a number of spectral maxima corresponding to small oscillations of the time-dependent spectral intensity. Moreover, one can see the spectral maxima at the frequencies ν = 8.82 THz, 9.87 THz, which are close to 2,4-DNT absorption frequencies: ν = 8.8 THz and 9.9 THz. It should be stressed that the spectral resolution in (c-e) is equal to ∆ν = 0.01 THz.
The similar result occurs for the frequency ν = 4.95 THz, which spectral line dynamics are depicted in Figure A4 computed at the window shift ∆ = 0.005 ps (a) and 0.4 ps (b). The magnified view of them is also depicted in Figure A4  Using as example the 2,4-DNT substance, we illustrate below an influence of the window shift Δ on the accuracy of constructing spectral line dynamics at frequencies ν = 1.1, 4.95, 8.8 and 9.9 THz because it is very important question for the problem under consideration. In Figure A2 the 2,4-DNTspectral line dynamics at the frequency ν = 1.1 THz computed at the window length T = 2.8 ps and with increasing window shift Δ = 0.005, 0.05, 0.1, 0.2 and 0.4 ps in the shortened time interval t = [5, 17.5] ps are shown. The shapes of all the dynamics (a-e) are very similar. However, in (a-c) they contain the oscillations at the frequency ν = 5 THz (close to ν = 4.95 THz) due to the small shift of window Δ ≤ 0.2 ps. If the window shift increases up to value Δ = 0.2, 0.4 ps, then these oscillations disappear. This fact is illustrated by the Figure A3(a), where the magnified view of the spectral line dynamics computed at the frequency ν = 1.1 THz at the increasing window shift Δ = 0.1, 0.2, 0.4 ps is depicted. Let us note that the spectral amplitude corresponding to the frequency ν = 5 THz decreases by 10 times at increasing window shift from Δ = 0.05 ps to 0.4 ps. For illustration of this, the corresponding part of the spectra is shown in Figure A3 Figure A3 (b). In order to estimate the frequency values of these oscillations, we apply the Fourier transform to the evolution of spectral amplitudes modulus | ( )| j P t ν (they are computed using Equations (A1) and (A2)). The spectrum of the corresponding spectral line dynamics computed at Δ=0.005 ps is shown in Figure A3 (e) in the frequency range ν= [8.7, 10.1] THz and it contains a number of spectral maxima corresponding to small oscillations of the time-dependent spectral intensity. Moreover, one can see the spectral maxima at the frequencies ν = 8.82 THz, 9.87 THz, which are close to 2,4-DNT absorption frequencies: ν = 8.8 THz and 9.9 THz. It should be stressed that the spectral resolution in (c-e) is equal to Δν = 0.01 THz. The similar result occurs for the frequency ν = 4.95 THz, which spectral line dynamics are depicted in Figure A4 computed at the window shift Δ = 0.005 ps (a) and 0.4 ps (b). The magnified view of them is also depicted in Figure A4     The similar result occurs for the frequency ν = 4.95 THz, which spectral line dynamics are depicted in Figure A4 computed at the window shift Δ = 0.005 ps (a) and 0.4 ps (b). The magnified view of them is also depicted in Figure A4    Further, the spectrum of the spectral line dynamics computed at the window shift ∆ = 0.4 ps (j-l) also contains the spectral maxima at the frequencies ν = 1.08, 4.92, 8.82, 9.92 THz, which are close to 2.4-DNT absorption frequencies ν = 1.1, 4.95, 8.8, 9.9 THz. However, their amplitudes are several orders smaller than in (d-f) and (g-i). Therefore, using of the window shift ∆ = 0.4 ps > 0.05 ps does not allows us to see the response of a medium at the frequencies appeared due to the frequency conversion process, i.e., the spurious absorption frequency.
Thus, to see the spectral line dynamic at the high frequencies (for, example ν = 4.95, 8.8, 9.9 THz) it is necessary to decrease the window shift if the window length is constant.