Conservative finite-difference scheme for the problem of THz pulse interaction with multilevel layer covered by disordered structure based on the density matrix formalism and 1D Maxwell’s equation

On the basis of the Crank-Nicolson method, we develop a conservative finite-difference scheme for investigation of the THz pulse interaction with a multilevel medium, covered by a disordered layered structure, in the framework of the Maxwell-Bloch equations, describing the substance evolution and the electromagnetic field evolution. For this set of the partial differential equations, the conservation laws are derived and proved. We generalize the Bloch invariant with respect to the multilevel medium. The approximation order of the developed finite-difference scheme is investigated and its conservatism property is also proved. To solve the difference equations, which are nonlinear with respect to the electric field strength, we propose an iteration method and its convergence is proved. To increase the computer simulation efficiency, we use the well-known solution of Maxwell’s equations in 1D case as artificial boundary condition. It is approximated using Cabaret scheme with the second order of an accuracy. On the basis of developed finite-difference scheme, we investigate the broadband THz pulse interaction with a medium covered by a disordered structure. This problem is of interest for the substance detection and identification. We show that the disordered structure dramatically induces an appearance of the substance false absorption frequencies. We demonstrate also that the spectrum for the transmitted and reflected pulses becomes broader due to the cascade mechanism of the high energy levels excitation of molecules. It leads to the substance emission at the frequencies, which are far from the frequency range for the incident pulse spectrum. Time-dependent spectral intensities at these frequencies are weakly disturbed by the disordered cover and, hence, they can be used for the substance identification.


Introduction
The detection and identification of hazardous chemical, biological and other substances by using the remote sensing is in research focus. One of the ways to achieve this aim is THz TDS because the THz radiation is non-ionizing and most common substances are transparent to it. Moreover, many hazardous substances have spectral fingerprints in the THz range of frequencies [1][2][3][4][5][6][7][8][9][10]. As a rule, the substance identification is based on a comparison of the absorption frequencies for a substance under investigation with the corresponding frequencies from a database. This method is used not only for security applications but also for a molecular structure investigation [11][12][13][14][15], nondestructive inspection [16][17][18], medical application [19], food quality inspection [20]. Currently, nanoscale terahertz spectroscopy is being actively developed as well [21].
Despite the obvious advantages of THz TDS using for the substance detection and identification, there is some shortcoming of its application [22]. Apparently, under real conditions, the factors like opaque packaging, inhomogeneity of the substance surface, atmospheric humidity and others affect the THz spectroscopy results [23][24][25][26][27]. In addition, a substance is usually covered by ordinary materials (paper, clothes or some packaging) [28]. Being made from materials transparent for a THz radiation these covers often possess a microscopic mechanical structure in the sub-millimeter scale of lengths. This leads to cover density fluctuations and also to its thickness variations. Therefore, the dielectric permittivity of such covers is not homogeneous and can be considered as a disordered photonic structure in the THz frequency range. Thus, observing a THz pulse reflected from a neutral material, we can reveal the false absorption frequencies, which are not inherent to this material, and the material may be wrong considered as dangerous one [29][30][31][32]. The overcoming of this serious drawback of the THz TDS method mentioned above and its physical mechanism understanding are the challenging problems for developing the real detection systems.
Another significant feature of the pulsed THz TDS is the THz pulse broadening after its interaction with a substance. Usually, the spectrum of the THz pulse transmitted through or reflected from a medium contains frequencies, which are, for example, greater than the maximal frequency of the incident pulse spectrum. This effect cannot be explained by multi-photon absorption of the THz radiation, because the THz pulse intensity is usually insufficient for its appearance. Using the computer simulation we have demonstrated a new spectral harmonic generation mechanism due to the absorption of a photon sequence (cascade mechanism), which results in the molecule higher energy levels excitation [32]. This mechanism may be observed in a great number of materials, possessing of vibration or rotation energy levels. We note that a photon sequence absorption belonging to visible frequency range was observed at studying the two-photon luminescence in a semiconductor [33].
It is well-known that the excited molecules relax to the equilibrium states. Therefore, the emission frequencies in the reflected or transmitted signal spectrum appear. The substance emission at high frequencies can be used for the substance identification because a radiation with a shorter wavelength is less affected by the cover material. The computer simulation is based on the Maxwell-Bloch equations [34], [35].
We emphasize that in our knowledge a problem statement under consideration is novel one and has not been investigated earlier. There are, at least, three features, which allow to state this. The first of them is an analysis of the interaction of the THz pulse containing a few cycles with disordered (or ordered) structure. The well-known Transfer Matrix Method for the description of the electromagnetic radiation interaction with the periodic structure is inapplicable in the case under consideration.
Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure The second one is the following. An active medium (i.e. the medium capable of absorbing and emitting the THz waves), which is placed inside the disordered structure, can be considered as an impurity of the periodical structure. In optics, the sizes of such impurities are much less than the sizes of the photonic crystal elements. In our case, the active medium length is comparable with or greater than a length of the disordered structure layers.
The third of the mentioned features is the description of the active medium response to the THz pulse action. We describe its characteristics evolution in the framework of a matrix density formalism.
Usually, computer simulation of the electromagnetic field interaction with a substance describing in the framework of the Maxwell-Bloch equations is based on a combination of the FDTD method [36] with the step-split method [34], [37]. As we believe, for the first time, the Crank-Nicolson scheme for the time discretization of the Bloch equations was applied for a two-level medium in [38]. However, the Crank-Nicholson scheme does not possess the positiveness property of the density matrix diagonal elements for three and more energy levels. In [39] it was shown that negative values of the energy levels populations may even occur in numerical simulation. In the same paper [39] and later in [40] a step-split method was used for solution of the matrix density elements equations with respect to the multilevel medium. This approach guarantees the positiveness of the density matrix diagonal elements but it leads to the decoupling of the Maxwell-Bloch equations. The latter feature leads to non-preserving the integrals of motion (invariants).
Below we derive a conservation law, which is a generalization of the Bloch invariant with respect to the multilevel medium. For computer simulation, we use a finite-difference scheme based on the Crank-Nicolson scheme supplemented by artificial boundary conditions, which approximate the well-known solution of the 1D wave equation [41]. Using the Cabaret scheme [42], we achieve the second order of accuracy for its approximation. Essentially that the proposed finite-difference scheme also possesses the conservativeness property with respect to the derived Bloch invariant and the invariants of the Maxwell's equations. To realize the finite-difference scheme we apply an iteration process and prove its convergence. Computer simulations show that the positiveness property of the density matrix diagonal elements is satisfied by choosing a mesh time step, which is defined by the problem parameters.
This article is organized in the following way. In the Section 2, we state the problem and derive its invariants.
In the Section 3, we develop the conservative finite-difference scheme for the Maxwell-Bloch equations formulated in the previous Section. Since the finite-difference scheme is a nonlinear one, we solve it by using an iterative method. We prove the theorem of the iteration convergence.
Finally, in the Section 4 we discuss the computer simulation results for the THz pulses transmitted through (or reflected from) an uncovered medium or a medium covered by disordered structure. Then we discuss the differences in spectra of these pulses and a physical mechanism, which causes such differences.

Problem statement
Initially, a THz pulse with a few cycles is located in vacuum and falls on an optically active medium (i.e. a medium, capable of electromagnetic field energy absorbing or emitting at specific frequencies). This medium is placed inside a disordered layered structure. Each of layers is characterized by its inherent dielectric permittivity and all layers are transparent for the THz radiation. The latter means that the layers do not possess any absorption frequencies in the Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure THz frequency range. Below, the term "medium" will denote the active medium and the term "cover" or "covering" will denote a disordered layered structure.
A scheme of the laser pulse propagation is shown in Fig 1. Several important coordinates are marked in the figure, namely: positions of electric field detectors E refl and E trans ; coordinates of the active substance faces z L and z R , and coordinates of the cover faces z c1 and z c2 . The pulse propagates in vacuum, then transmits through the left cover, a medium, and finally, through the right cover, and exits into vacuum to the right. The pulse is partly reflected from various boundaries between layers and from the medium faces too. A certain part of the pulse energy is absorbed by the active medium. Generally speaking, a part of the absorbed energy can be emitted by the medium due to radiative transitions (emission) between excited energy levels. Reflected and transmitted pulses are detected near the left and right boundaries of the computational domain at the sections E refl and E trans . As a rule, these pulses may consist of a sequence of the THz sub-pulses.

Maxwell-Bloch equations
We use the Maxwell's set of equations describing the THz pulse propagation in 1D case: together with equations with respect to the matrix density elements (see below). For convenience, we write above the set of equations in both systems of physical units. Here H, E are the magnetic field strength and the electric field strength, correspondingly. B, D are the magnetic field and electric field induction strength, correspondingly, z and t are the spatial coordinate and time, c = 3Á10 8 mÁs -1 is the light velocity in a vacuum. Parameters ε 0 and μ 0 in Eqs (1)- (3) are the permittivity of free space and the magnetic permeability of free space and their values are equal to ε 0 = 8.85Á10 −12 FÁm -1 , and μ 0 = 4πÁ10 −7 HÁm -1 = 4πÁ10 −7 TÁm A -1 , correspondingly [43]. Below we consider the case of a non-magnetic medium, thus, in Eq (3) the magnetic permeability is equal to μ = 1. P describes the medium polarization, L t is a time interval during which the electromagnetic field propagation is analyzed. Parameters L l and L r (see Fig 1) denote spatial coordinates of the domain beginning and its ending. The total length of the domain is equal to L z . = L r -L l . Coordinates z R and z L mark the positions of the medium boundaries, L s . = z R -z L is used to denote its length. We consider a TE mode of the THz pulse propagation. The THz pulses reflected from and transmitted through the medium with covering are computed for a number of random realizations for the dielectric permittivities and then the result is averaged. This procedure is usually applied in physical experiments also to suppress a random fluctuation influence on the measured THz pulse characteristics. Let us notice that, each of two covering structures contains N l layers with the layer thickness h l (see Fig 1). Dielectric permittivities at each of the layers we denote as ε j , ε j 0 , j = [1, N l ], for left and right covering, correspondingly.
Obviously, in vacuum (z 2 [L l , z c1 ] [ [z c2 , L R ]) the electric field induction is equal to the electric field strength: D = E. In the cover layers (z 2 [z c1 , z L ] [ [z R , z c2 ]) the polarization is assumed to be proportional to the electric field strength and the Eq (3) is replaced by the following formula: where the dielectric permittivity ε(z) for each of the layers varies randomly in the range ε l = [ε min , ε max ] and it is uniformly distributed.
To describe the medium polarization (z 2 [z L , z R ]) we use the matrix density formalism. In this case the medium is described by the following set of equations: ðd mq r qm À r mq d qm Þ; ð7Þ Here, ρ mn = ρ mn (z,t) is the element of the density matrix describing the state of a molecule, which possesses N energy levels, in a time moment t and at the medium section z. Diagonal elements of the matrix ρ mn describe the number of molecules at the energy level m. These Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure elements obey the condition X N m¼1 r mm ¼ 1 and obviously they are non-negative. The off-diagonal elements ρ mn are related to the medium polarization corresponding to a transition between n and m energy levels. γ mn , ω mn and d mn are the relaxation rate of non-diagonal elements, the transition frequency, and the dipole moment. They are associated with a transition from energy level m to n, correspondingly, W mn is the population decay rate. Parameter M is the density of "active" molecules in the medium, χ corresponds to the non-resonant linear electric susceptibility of the substance, ћ is the reduced Plank constant. It should be noted that the Eqs (6)-(8) describe both linear and nonlinear interaction of the electromagnetic field with the multilevel medium in the dipole approximation. The term ðd mn r nm Þ in the expression for the medium polarization (8) corresponds to a resonant or near-resonant response of the medium to the THz pulse action and causes the electromagnetic field energy absorption and emission of the medium. The second term χE in the same equation describes a linear non-resonant polarization of the medium arising because of the pulse interaction with the medium which has transition frequencies distant from the pulse carrier frequency. This term is responsible for the Fresnel reflection from the medium boundaries and plays the key role in the formation of reflected pulses.
The initial condition for the electromagnetic field is as follows: Here E 0 is the incident electric field strength amplitude, H 0 is the incident magnetic field strength amplitude, z p the incident pulse center position, ν p and ω p = 2πν p are the carrier and angular pulse frequency, a z is the pulse length in space, τ p = a z /c is the pulse duration. The amplitudes for E and H are connected by well-known formulas (see below) at the electromagnetic wave propagation in a linear medium. We suppose that the medium is in the ground state until an action of the THz pulse: Obviously, this condition corresponds to the absolute zero value of temperature. Nevertheless, for the investigation of the cascade mechanism manifestation we consider "pure conditions". Under room temperature, the effect will be still pronounced but its interpretation is difficult because a number of possible energy level transitions take place at the electromagnetic pulse action on the medium.

Dimensionless equations
For computer simulation, it is useful to transform the Eqs (1)-(10) to a dimensionless form. We use the following normalization, where the index "dim" means dimensionless variables: Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure Here τ 0 is chosen to be equal to 10 −12 s as a unit of time measurement. This value corresponds approximately to 0.3 mm. We also introduce the following dimensionless coefficient: Note that in the SI base units, this coefficient has the form [44]: For the dimensionless values of E, H and P we write the following relations: Using Eqs (11)-(15), Eqs (1)-(10) can be rewritten in the dimensionless form: Dðz; tÞ ¼ εðz; tÞ Á Eðz; tÞ; z 2 ½z c1 ; z L [ ½z R ; z c2 ; ð19Þ ðd mq r qm À r mq d qm Þ; ð22Þ For brevity, an index "dim" is omitted below.

Maxwell's equations invariants.
To verify the computer simulation accuracy we follow the problem invariants, which can be easily obtained by integrating Eqs (16) and (17) along z and t coordinates.
After the integration, the following invariants can be written: Here, I 1 and I 2 are the constants determined by the incident pulse parameters. If the electromagnetic field does not reach the computational domain boundaries, the invariants (28) and (29) are transformed:

Bloch invariant. Theorem 2.1.
Provided relaxation rates W mm and γ mn are negligible, the density matrix evolution described by Eqs (21) and (22) obeys the following invariant: Proof. Firstly, let us introduce auxiliary definitions: the purpose of which is to reduce the complexity of the forthcoming operations. Although we derive the invariant supposing the influence of relaxation determined by γ mn negligible, it is still present in our consideration in order to demonstrate its role in invariant deterioration. Secondly, let us remind that the matrices d, ρ and and O are Hermitian ones: Note here, that O mm = 0, O mn + O nm = 2γ mn , because the matrices ω and γ are antisymmetric and symmetric real-valued matrices correspondingly, and the main diagonal elements of matrix γ is equal to zero: Using theses notations, Eqs (6) and (7) can be re-written in the following form: Let us further transform Eq (36) separating the terms in the sum, which contain diagonal and off-diagonal elements: For brevity, we introduced above additional notations: where δ mn = -δ nm is the difference of the energy levels populations. Now, we use well-known formula and Eq (35) to obtain the time derivative for the squared matrix elements, which are complex values: Thus, multiplying Eq (37) by r Ã mn and summing the result with its complex conjugate, and taking into account Eq (39), we obtain the following expression: Where It is now seen from Eq (40) that the non-zero relaxation (γ mn 6 ¼0) causes exponential decay of each of terms in the second sum in the formula (32). Therefore, to derive the invariant, we assume below that γ mn = 0. In general case, one can write a law of density matrix elements decaying.
Next, let us construct evolution equations for the squared population differences between two energy levels δ mn . After trivial transformations, it can be obtained from Eq (39) in the following form: Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure Now, let us perform the summation of Eqs (40) at zero relaxation (γ mn = 0) using Eq (38) and also the summation of Eq (42) for all the pairs of indices m and n such that 1 n<m N: ðd mq r qn r nm À d qn r mq r nm þ d nq r qm r mn À d qm r nq r mn Þ þ φ mn Þ; ð43Þ We chose such summation limits using the analogy with the derivation of the Bloch invariant in the case of two energy levels medium and keep the time derivatives as separate sum for clearness of further transformations.
Let us notice that the sum contains all combinations of three distinct integer numbers in the range 1; N . In combinatorics, the term "combination" means "a way to choose k distinct elements from a set of l elements without regard to the order of selection. And word "triple" below will be used to denote "a combination of three elements from the range 1; N ".
For each trio of distinct numbers, we can assign names a, b, c to them in such a way so that a<b<c. Let us now examine one such trios that is spanned according to the sum limits. There exists only three possible ways to assign numbers a, b, c to indices m, n, q. The possible assignments are listed in the Table 1. Now let us examine the terms in Eq (43) related to the specific trio a, b, c: ¼m;n q¼1 d mq r qn r nm À d qn r mq r nm þ d nq r qm r mn À d qm r nq r mn þ m; n; q 6 ¼ b; a; c; m; n; q 6 ¼ c; a; b ; m; n; q 6 ¼ c; b; a þiEðd bc r ca r ab À d ca r bc r ab þ d ac r cb r ba À d cb r ac r ba Þþ þiEðd cb r ba r ac À d ba r cb r ac þ d ab r bc r ca À d bc r ab r ca Þþ þiEðd ca r ab r bc À d ab r ca r bc þ d ba r ac r cb À d ac r ba r cb Þ; m; n ¼ 1; N : The last three terms form the only possible set described by a particular combination of three numbers. These three terms cancel out each other for any arbitrary combination of numbers and the second sum in Eq (45) contains all such combinations. Thus, the second term in Eq (45) is equal to zero as well. Hence, the all-pairs sum of the squared off-diagonal elements is equal to: Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure Now, let us perform the same analysis for the population differences: Here it is important that definition of δ mn . Eq (38) provides that: Taking into account Eq (48) along with the definition (41), the Eq (47) can be transformed into: We see that our analysis demonstrates the appearance of an additional term −2(φ ba + φ ca + φ cb ) in Eq (49), which contains φ mn indexed by all three possible two number combinations in the triple a, b, c. There are N 3 of various triples that the second term contains, in these combinations of three elements each possible combination of two elements is repeated N-2 times in total. So, φ mn associated with each particular pair of indices appears in the sum N-2 times as well and Eq (49) can be transformed into: Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure which is simplified further: Now it is evident, that we can obtain the required invariant in time (32) by multiplying Eq (46) by 2N and summing it with Eq (51).

Finite-difference scheme
To solve Eqs (16)-(27), we use a combination of the well-known Yee's scheme [36], [45] for Maxwell's equations and a Crank-Nicolson scheme for the equations with respect to the density matrix with an iteration process. Therefore, to realize the positiveness property of the density matrix diagonal elements we use a sufficiently small mesh steps in computer simulation.

Approximation of the equations
The model (16) Here N z and N t are numbers of nodes along spatial coordinate and time, correspondingly. Let us also introduce the following notations for nodes: which correspond to the incident pulse centre coordinate; the coordinates of the disordered structure beginning and ending and coordinates of the medium beginning and ending and the number of nodes inside single cover layer, correspondingly (see Fig 1). The square brackets in Eqs (56) and (57) mean the integer part of the fraction. (We always use parameters such that these indices are integer without rounding.) Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure The mesh function ε(z j ) for the dielectric permittivity is defined on the grid ω z : We note that the dielectric permittivity is defined in the same nodes as the electric field induction and electric field strength. Mesh function for the magnetic field strength is defined on the gridÕ : Hðz jþ0:5 ; t lþ0:5 Þ. The mesh functions for the electric field strength and its induction as well as for the mesh functions for density matrix elements and medium polarization are defined on the grid O: D(z j , t l ), E(z j , t l ), ρ mn (z j , t i ), P(z j , t l ), correspondingly. For brevity, we save the same notations for the mesh functions as for the differential functions introduced above.
Eqs (16) and (17) are approximated using the Yee's method. Therefore, in our notations the corresponding finite-difference scheme is written as Hðz jþ0:5 ; t lþ0:5 Þ À Hðz jþ0:5 ; t lÀ 0: Dðz j ; t lþ1 Þ À Dðz j ; t l Þ t ¼ À Hðz jþ0:5 ; t lþ0:5 Þ À Hðz jÀ 0: In vacuum, the electric field strength is equal to its induction: and in the cover layers their relation is defined as: A relation between the mesh functions for the electric field strength and its induction for the medium is written below after writing the finite-difference scheme for the density matrix.
The novel feature of this well-known finite-difference scheme consists of the multilevel Maxwell-Bloch equations approximation for a medium response description. So, the Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure approximation of Eqs (20)- (23), based on the Crank-Nicolson scheme, is made in the following form: ðd mqrqn Àr mq d qn Þ; ð63Þ ðd mqrqm Àr mq d qm Þ; ð64Þ

Boundary conditions approximation
At the edges of the computational domain, Eqs (1) and (2) are replaced with equations supplying a one-directional propagation, which makes boundaries transparent and reflectionless.
In the boundary node j = N z we use the Cabaret approximation [42] for the translation operator:

Initial condition approximation
We see that these Eqs (70)-(73) are not applicable on the first time step l = 1. Since we consider the finite incident electromagnetic pulse, therefore we can define the initial mesh functions for the electromagnetic field distribution in the following way: Hðz jþ0: Here j p is the mesh node index corresponding to the position of the initial pulse center: Note that the initial values of H are defined at t 0.5 node, which corresponds to the propagation of a wave packet in positive direction of z-coordinate. Therefore, taking into account the wave equation we can define the zero-value mesh functions in the boundary nodes on the second time layer (l = 1): The initial state of the medium is approximated in the following way:

Approximation order investigation
We stress that the difference analogues of the invariants (82), (83) and (86) preserve their values during the computation time with the accuracy not worse than 10 −6 for the mesh steps being equal to 0.01 (the full set of parameters are shown in the next section). Invariant (86) takes place only in the case of zero relaxation rate of off-diagonal elements. In our simulation, relaxation is non-negligible and therefore, the invariant value decreases over time as expected. In current paper, we have to introduce this invariant to explain peculiarities in population dynamics (section 2.3.2).

Iterations convergence proving
The difference Eqs (59)-(69) are nonlinear ones in respect to the electric field strength E and the density matrix elements ρ mn . To solve the equations, the following iteration method is E s ¼ 0:5ðE s ðz j ; t lþ1 Þ þ Eðz j ; t l ÞÞ;r s mn ¼ 0:5ðr s mn ðz j ; t lþ1 Þ þ r mn ðz j ; t l ÞÞ; s ¼ 1; 2; . . . ; ð95Þ The functions on zero-iteration are chosen equal to their values on the previous time layer: The iteration process is stopped if both of the following conditions are satisfied: jr sþ1 mn ðz j ; t lþ1 Þ À r s mn ðz j ; t lþ1 Þj < e 1 jr s mn ðz j ; t lþ1 Þj þ e 2 ; ð98Þ where e 1 and e 2 are positive constants, characterizing the maximal absolute and relative accuracy of the iteration process. In our computer simulations, their values were chosen as 10 −8 and 10 −10 correspondingly. The procedure for the solution of the system (59)-(62) and (92)-(95) is the following. First, from Eq (59) we obtain the value of H(z j+0.5 , t l+0.5 ) on the next time layer, then we have D (t l+1 , z j+1 ) from Eq (60) and E s+1 (z j , t l+1 ) from Eq (94). We substitute the last value in the Eqs (92) and (93) at the next iteration step, etc. Proof. To prove this theorem, we have to demonstrate that (92)-(95) is a contraction, i.e. that all iterations are uniformly limited in a certain difference norm kÁk h : where C E and C ρ are constants, and that the following inequalities: jjr sþ1 À r s jj h qjjr s À r sÀ 1 jj h ; q < 1 ð101Þ Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure are valid. We use the difference norm C: In Eq (102) the letter "f" denotes mesh function defined on the mesh. We estimate the function modules on different iteration steps and show that they are uniformly bounded in the norm (102). Then considering the iteration differences kE s+1 − E s k C and kρ s+1 − ρ s k C , and using estimations similar to those applied for the iteration modules we get the inequalities (102) in the norm (102).

Numerical stability of the finite-difference scheme
A stability condition for the finite-difference scheme (59)-(69) is the Courant-Friedrichs-Levy (CFL) condition [46], because the developed scheme involves the Yee's one. When the Maxwell's equations are solved, the CFL condition is well-known for the linear wave equation and has the form u h/τ, where u is the wave propagation velocity [46]. In Eqs (59) and (60) the wave velocity u for a linear medium with ε = 1 (vacuum) is equal to unity (u = 1), the mesh steps are chosen to be equal: h = τ, and the CFL condition reduces to the inequality u 1. Thus, for the wave propagation in a vacuum the CFL condition is valid: u = 1. In the active medium as well as in the disordered structure, the wave velocity is less than unity because the dielectric permittivity is greater than unity and u ¼ 1= ffiffi ffi ε p in the dimensionless units. Therefore, the CFL condition is also fulfilled.
It should be stressed that the preservation of the invariants (Theorem 3.1) is an additional confirmation of the solution correctness.

Computer simulation results
The aim of investigation provided below consists in an analysis of the pulses reflected from and transmitted through a three energy-level medium covered by disordered structure. We analyze the pulse structure (it consists from several sub-pulses) and its spectrum, and the spectra of sub-pulses. We pay attention to the substance emission at the transition 3-1 (high frequency), which appears due to the cascade mechanism of high energy level excitation and we demonstrate the evidence of its manifestation.

Computer simulation parameters
The computer simulation results for substance with covering are averaged over 32 random realizations for the cover dielectric permittivity random distribution. Usually, in physical experiments [29] we made also 16 or 32 realizations at the measurements of THz signals with duration about 100 ps. If we use more random realizations, then the measurement time increases up to 1 min, which is not acceptable in practice.
The electric field strength unit E 0 was chosen to be 1.05Á10 6 V/m. Using the well-known formula, we can estimate the peak intensity of the pulse in SI units as This formula yields the maximal intensity of the pulse equal to 3Á10 5 W/cm 2 . Another wellknown formula [46] gives the value of the magnetic field induction and the magnetic field strength: Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure where n ¼ ffiffiffiffiffi ffi εm p is the refractive index. Taking into account that in a vacuum ε = μ = 1, we obtain: B 0 = 3.5Á10 −3 T, H 0 = 2785 A/m (or 27.85 A/cm). Evidently, for practice the intensity less than 10 6 W/cm 2 is more reasonable but in computer simulation, the effect of second harmonic generation is better pronounced for greater than this value of the pulse intensity.
We chose the mesh steps as h = τ = 0.01. Iteration process parameters are chosen as e 1 = 10 −8 and e 1 = 10 −10 , their decreasing does not result in notable changes of the computer simulation results at chosen mesh steps. Other dimensionless coordinates are given in the Table 2. The parameters of the medium and the pulse are provided in the Table 3.

Transmitted and reflected pulse structure
The computer simulation is made for two cases: with the cover presence and without it to clarify its influence on the reflected and transmitted pulse spectra. The typical signals are depicted in Fig 2. One can see that in both cases the spatial pulse length is comparable with the length of a medium. The pulse splitting into a few sub-pulses is clearly seen for both covered and uncovered media. However, if the cover is absent then there are only three noticeable sub-pulses, which are caused by THz radiation multi-reflections from the left and right medium layer faces and interference pattern induced by these reflected signals. The first reflected sub-pulse in Fig  2A is usually named as the main sub-pulse and it occurs directly due to the THz pulse reflection from a medium layer face. Then, we observe much more sub-pulses in the reflected signal. They can usually be exploited with success for the substance detection and identification.

Fourier spectra
The Fig 2 demonstrates the comparison of the signal spectra, which are calculated as the absolute value of the Fourier transform using the well-known formula: Here L is the time interval (in this case, L = L t , N 0 = N t = 10 4 ) or the spatial computational domain length (L = L z , N 0 = N z = 3200). This transform is applied to the measured signals E refl , E tran and similar Fourier transform is also applied to the incident pulse electric field strength distribution E(z, t = 0). Note that in physical experiments the frequency ν = ω/2π is usually measured, that is why in the figures below we use this ordinary frequency instead of the frequency ω.
As follows from the Fig 3A and 3B, the spectra are significantly disturbed because they possesses a strong modulation and it is even more powerful if a cover is absent. This is a consequence of the interference of different sub-pulses depicted in Fig 2. Obviously, the spectrum minima can be considered as the absorption frequencies and, thus, they may be treated as the false absorption frequencies [31]. To eliminate these mistakes, the sub-pulses should be analyzed separately.  Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure At the same time, if covering of the medium is present, the THz pulse spectrum is not so much periodically modulated because the reflected signal is averaged over 32 random realizations. However, they exhibit irregular minima of various depth, which could be also interpreted as the absorption spectral lines of some substances while they are simply caused by the pulse multiple reflections from the faces of various layers and, consequently, an interference of appeared sub-pulses produces this spectrum modulation.
Other way for the substance identification may be achieved at using the substance emission frequency analysis. For the substance under consideration, a medium emission is present due to the energy level transition 3-1. We see that this frequency is absent in the incident pulse spectrum. Its appearance in the reflected pulse spectrum is due to the cascade mechanism of high energy level excitation of substance molecules. Because the spectrum shape near the frequency 1.75 (Fig 3C and 3D) is not distorted despite the presence of covering it is possible to use it for the substance identification. In Fig 4 the energy level population evolutions in the chosen section (z = 0) is shown. We stress that in our computer simulation, the parameters correspond to optically thin layer. Thus, the general behavior of the energy level populations in any medium sections is qualitatively similar to each of them. Hence, we depict a time evolution of the density matrix elements for a single section only. In these figures, it may be not evident that the third energy level is populated as a result of the 2-3 energy-level transition. Indeed, in Fig 4A, we see notable  Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure oscillations of the third energy level population and these oscillations are not synchronous to the second energy level population changing. Let us notice that in [32] we have presented results corresponding to strong electric field strength amplitude (E 0 = 2) at studying the transmitted signal only. In that case, the third energy level population increases after the second energy level population becomes high enough. At significantly lower electric field strength, it is necessary to use the correlation coefficient between a pair of the energy level populations as an estimation of the cascade mechanism presence for the high energy level excitation.
The correlation coefficient between two series of measurements containing N' numbers of data is given by the following formula: Comparison of the incident pulse spectra (E inc ) with the spectra reflected from (solid line) or transmitted through (dashed line) the medium without cover (a) and spectra averaged over 32 random realizations for the dielectric permittivity in layers of a disordered structure (b). The spectra near the frequency, corresponding to 3-1 energy level transition, for uncovered and covered medium correspondingly (c, d) (dimensionless units). https://doi.org/10.1371/journal.pone.0201572.g003 Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure where, a mean , b mean are the mean values of the measured variables. This formula application to the energy level populations in the time interval [10,35] with N' = 2500 yields the following correlation magnitudes: Thus, from Eq (106) it can be seen that the populations of the first and second energy levels are tightly connected and are uninfluenced by the instantaneous value of the electric field strength (107). On the other hand, the third level population is significantly correlated with the electromagnetic field due to non-resonant interaction. That means, when the electromagnetic pulse propagates in a medium, the electric field causes its polarization, related to nondiagonal elements of the density matrix ρ.
In accordance with invariants of Maxwell-Bloch equations (32), the appearance of the polarization induces temporary shifts in the energy level populations. In the case of three level medium, the invariant is written as follows:  Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure As soon as external electric field strength decreases, the substance polarization decreases also, and consequently the populations return to their initial values. This is non-resonant process unrelated with the medium emission.
Nevertheless, the medium polarization is not only the process influencing the third energy level population evolution. High correlation between the second and third level corr(ρ22, ρ33) is more notable. That means that the second energy level population is significant factor governing the third energy level population evolution and the third energy level excitation takes place due to the cascade mechanism.
Thus, there are three processes that influence the population of the third level. One of them is the medium non-resonant polarization. Other two of them are resonant transitions between the second and third energy levels and between the third and first energy levels. The first process can cause the first growth of the third energy-level population ρ 33 before of the second level population ρ 22 growth. This process is not accompanied by electromagnetic field energy absorption and does not result in non-zero value of the third energy level population changes after the pulse action ending. Therefore, as soon as the field strength decreases, the population returns to the initial state without any medium emission. Other two processes are responsible the population relaxation process, which is several orders of magnitude slower than the energy level excitation processes.

Spectra of sub-pulses
In our previous papers [29], [30] we have shown that the detection and identification of a substance using the THz radiation measured in the reflection mode can be effective at analyzing the sub-pulses following the main pulse reflected from the substance face. In this case, we get information about the substance absorption frequencies. Below we provide a similar analysis for a medium covered by disordered structure. However, another our purpose is to observe also a substance emission at the frequencies corresponding to transitions from the high energy levels excited due to the cascade mechanism. To obtain the sub-pulse spectra, the Fourier-Gabor transform [48] is applied to the signal E(Z refl\tran , t k ) in the corresponding time intervals.
It should be noted that if the signal amplitudes at the left and right ends of the time interval differ then the effect, which is known as the spectrum spreading [49], occurs. As a result, false frequencies appear in the signal spectrum. They can be suppressed by using the window weight function that tends to zero toward the ends of the time interval. For this purpose, at each time interval, which contains the corresponding sub-pulse, we use the Fourier-Gabor transform [48]  which is a near-rectangular window. Therefore, we obtain: E I;II;III;... reflntran ðo n Þ ¼ j where E I,II,III. . . Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure time interval. In our simulation, this parameter was taken to be equal to 20 for the window function G(t) to fall off sharply. Fig 5 shows the sub-pulse spectra for the signals reflected from and transmitted through an uncovered medium. The transform window edges and the signal partitioning can be seen in the Fig 2A. In Fig 5A the spectrum modulation is absent at all. Apparently, the first reflected sub-pulse depicted in the Fig 5A is formed due to the reflection by the left face of the medium and thus, does not propagate inside the medium and, therefore, does not interact with it. This reflected pulse spectrum contains no information about the medium absorption frequencies, and the spectrum shape is the same as that for the incident pulse (Fig 3A).
The second sub-pulse ( Fig 5B) appears due to multi-reflections from the right and left faces of the medium before being detected at the left section E refl . The absorption frequency corresponding to the transition 1-2 (the frequency 0.8) is noticeable in its spectrum as well as weak spectral intensity maximum at the frequency 1.75 corresponding to the cascade mechanism of high energy level excitation. The absorption frequency corresponding to the energy level transition 2-3 is not seen (Fig 5B) because of weak increasing of the second energy level population. As follows from Fig 4A, the latter remains quite low: thousand times less than the  Fig 2) spectrum (a, c) and the second sub-pulse (area II in Fig 2) spectrum for the signal reflected from and transmitted through a medium without covering. The spectra (c, d) are shown in the higher frequency range (dimensionless units). https://doi.org/10.1371/journal.pone.0201572.g005 Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure population of the first energy level. Therefore, the spectrum dip, corresponding to this energy level transition, is quite small.
In opposite, both transmitted sub-pulses exhibit notable spectrum peculiarities. The first one contains the absorption frequency 0.8. The second sub-pulse spectrum (Fig 5B and 5D) exhibits a substance emission at the frequencies 0.8 and 1.5. Hence, we observe the spectral maximum instead of spectral minimum at the frequency corresponding to the 1-2 energy level transition. The frequencies corresponding to the energy level transition 2-3 or 3-2 is not observable at all here. It should be stressed that in the spectrum of the second sub-pulse transmitted through a medium, some additional false absorption frequencies occur.
Let us discuss the corresponding spectra of the sub-pulses reflected from or transmitted through the medium covered by disordered structure. The temporal partitioning of the subpulses can be seen in the Fig 2B. In Fig 6 the corresponding sub-pulse spectra are depicted. We see that the main reflected pulse spectrum in Fig 6A is slightly different from the incident pulse spectrum.
The main transmitted pulse, located in the time region II (Fig 2B), exhibits the absorption frequency 0.8 presence in the spectrum (Fig 6B). Nevertheless, the pulse spectrum contains the false absorption frequency, which is equal to 0.7. The second reflected sub-pulse spectrum ( Fig  6B) contains only the false absorption frequencies, which do not belong to a medium, and they are a result of superposition of the pulses caused by electromagnetic field multi-reflections from the layered structure.
In the time area III (Fig 2B), the second transmitted and the third reflected pulses are located (Fig 6C). The electric field strengths of these pulses are 10 times less than the ones for previous sub-pulses. Nevertheless, the emission at the frequency 0.8 related to 2-1 energy Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure level transition can be observed in the transmitted sub-pulse spectrum. At the same time, the absorption frequency corresponding to 2-3 energy level transition is noticed together with several absorption frequencies. We stress that the third reflected sub-pulse is the first reflected one that contains the absorption frequencies related to the medium. However, it contains also the false frequencies about of 0.5 and 0.6 dimensionless units, and there is a spectrum minimum at the frequency 1.01, which is shifted away from the 2-3 transition resonance frequency.
The spectrum of the third transmitted sub-pulse, located in time interval IV, contains maximum at the frequency about of 0.8 (this means a substance emission presence (Fig 6D)) and it contains also a number of false absorption frequencies. The fourth reflected sub-pulse is distorted as well. In the frequency range ν<0.7 the spectrum shape is similar to that of the third reflected sub-pulse. Thus, the same mechanism of the spectrum distortions exists. However, the other part of the spectrum is quite different: the spectrum minimum at the frequency 0.8 is still present but it is shifted to the range of lower frequencies, and the spectrum minimum at the frequency 0.9 is observed, unlike the third sub-pulse in which the spectrum minimum is located at the frequency 1.01. Fig 6E demonstrates comparison of reflected sub-pulses spectra near the frequency 1.75 corresponding to 3-1 energy level transition. For the first and second sub-pulse, the spectral intensity maximum at this frequency is absent. The third and fourth sub-pulse spectra have the maximum near the frequency 1.75. It means that those sub-pulses can be used for the substance identification. However, the fourth sub-pulse is very weak and is significantly distorted even in this frequency range.
Thus, we have demonstrated the cascade mechanism of high energy level excitation and the emission frequency observation corresponding to the transition 3-1 in the spectra of both pulses reflected from or transmitted through the substance. Since this frequency is observed in the case of the covered medium, it may be used with success for the detection and identification of the substance.

Conclusions
We develop the conservative finite-difference scheme for the problem describing an interaction of the THz pulse, containing a few cycles, with the multilevel medium in the framework of the Maxwell-Bloch equations. This finite-difference scheme is based on the Crank-Nicolson scheme with using a new approximation of the artificial boundary conditions on the basis of the Cabaret scheme. For a solution of the corresponding nonlinear difference equations, the iterative process is proposed and its convergence condition is written.
We generalize the Bloch invariant for a multilevel medium and the conservativeness of the developed finite-difference scheme with respect to this invariant as well as to the invariants of the Maxwell's equations were proved. The proposed finite-difference scheme allows to control the positiveness property of the energy level populations at computer simulation by choosing the mesh step in appropriate way.
We have demonstrated that the disordered structure covering of a medium results in appearing of the false absorption frequencies. This effect arises from complicated interference of the multiple sub-pulses reflected from various faces between layers of the disordered structure.
We showed the appearing of the emission frequencies belonging to the range of high frequencies in the sub-pulses spectra due to the cascade mechanism of high energy level excitation for both transmitted and reflected THz signals. The disordered cover less distorts the Finite-difference scheme for the problem of THz pulse interaction with layer covered by disordered structure spectral intensities of these frequencies. Therefore, they can be used with success for the substance detection and identification.