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

PLOS ONE, Aug 2018

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0201572&type=printable

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

August 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 Vyacheslav A. Trofimov 0 1 Svetlana A. Varentsova 0 1 Irina G. Zakharova 1 Dmitry Yu. Zagursky 1 0 Faculty of Computational Mathematics and Cybernetics, Lomonosov Moscow State University , Moscow , Russia , 2 Faculty of Physics, Lomonosov Moscow State University , Moscow , Russia 1 Editor: Luis Carretero, Universidad Miguel Hernandez de Elche , SPAIN 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. - Data Availability Statement: All relevant data are within the paper. Competing interests: The authors have declared that no competing interests exist. 1. 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?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?15], nondestructive inspection [16?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?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?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. 2 / 30 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 MaxwellBloch 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. 2. 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 3 / 30 Fig 1. The scheme of the THz pulse interaction with covered substance. The left coordinate axis shows the electric field strength in dimensionless units. The dielectric permittivity of layers is shown on the right axis. The crossed rectangle depicts the active medium. Lines marked as Erefl, Etrans denote coordinates at which the reflected and transmitted THz pulses are measured. 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 Erefl and Etrans; coordinates of the active substance faces zL and zR, and coordinates of the cover faces zc1 and zc2. 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 Erefl and Etrans. As a rule, these pulses may consist of a sequence of the THz sub-pulses. 2.1 Maxwell-Bloch equations We use the Maxwell's set of equations describing the THz pulse propagation in 1D case: D ? E ? 4pP?CGS?; D ? ?0E ? P?SI?; B ? mH?CGS?; B ? m0mH?SI?; ?1? ?2? ?3? Lt; Ll < z Lr; ?4? 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 108 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, Lt is a time interval during which the electromagnetic field propagation is analyzed. Parameters Ll and Lr (see Fig 1) denote spatial coordinates of the domain beginning and its ending. The total length of the domain is equal to Lz. = Lr?Ll. Coordinates zR and zL mark the positions of the medium boundaries, Ls. = zR?zL 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 Nl layers with the layer thickness hl (see Fig 1). Dielectric permittivities at each of the layers we denote as ?j, ?j0, j = [1, Nl], for left and right covering, correspondingly. Obviously, in vacuum (z 2 [Ll, zc1] [ [zc2, LR]) the electric field induction is equal to the electric field strength: D = E. In the cover layers (z 2 [zc1, zL] [ [zR, zc2]) the polarization is assumed to be proportional to the electric field strength and the Eq ( 3 ) is replaced by the following formula: D?z; t? ? ??z?E?z; t?; 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 [zL, zR]) we use the matrix density formalism. In this case the medium is described by the following set of equations: ?dmqrqn rmqdqn?; i XN E ? q?1 i XN E ? q?1 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 ?5? ?6? ?7? ?8? rmm ? 1 and obviously they are non-negative. The off-diagoHere E0 is the incident electric field strength amplitude, H0 is the incident magnetic field strength amplitude, zp the incident pulse center position, ?p and ?p = 2??p are the carrier and angular pulse frequency, az is the pulse length in space, ?p = az/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: r11 ? 1; rmn ? 0; m 6? 1; n 6? 1: 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. 2.2 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: ?9? ?10? ?11? ?12? 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: Lldim ? Ll=ct0; Lrdim ? Lr=ct0; Ltdim ? Lt=ct0: k ? EM02t?0 ? ??sc1m??c m3??c0m:5g20g:51ss 11??2 ?CGS?: Note that in the SI base units, this coefficient has the form [44]: k ? M? ? ?m 3??m2kg1s 1? ?0E02t0 ?m 3kg 1s4A2??m1kg1s 3A 1?2?s1? we write the following relations: : For the dimensionless values of E, H and P Edim ? E E0; Hdim ? H H0; Pdim ? P E0: Using Eqs ( 11 )?( 15 ), Eqs ( 1 )?( 10 ) can be rewritten in the dimensionless form: Lr; Lr; D ? E; z 2 ?Ll; zc1? [ ?zc2; LR?; D?z; t? ? ??z; t? E?z; t?; z 2 ?zc1; zL? [ ?zR; zc2?; D ? E ? 4pP; ?dmqrqn rmqdqn?; ?13? ?14? ?15? ?16? ?17? ?18? ?19? ?20? ?21? ?22? ?23? ?24? ?25? ?26? ?27? 2.3 Invariants 2.3.1 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: ?r L Ll ?r L Ll ?t 0 ?t 0 D?z; t?dz ? I1 ?H?Lr; t0? H?Ll; t0??dt0; H?z; t?dz ? I2 ?E?Ll; t0? E?Lr; t0??dt0: Here, I1 and I2 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: ?r L Ll ?r L Ll D?z; t?dz ? I1 ? const; H?z; t?dz ? I2 ? const: 2.3.2 Bloch invariant. Theorem 2.1. Provided relaxation rates Wmm and ?mn are negligible, the density matrix evolution described by Eqs ( 21 ) and ( 22 ) obeys the following invariant: XN Xm1 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: dmn ? dnm; rmn ? rnm; Omn ? Onm; m; n ? 1; N : Note here, that Omm = 0, Omn + Onm = 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: omn ? onm; gmn ? gnm; gmm ? 0; m; n ? 1; N : ?28? ?29? ?30? ?31? ?32? ?33? ?34? ?35? X N q?1 dnn?rmn dmndmn?; m; n ? 1; N : For brevity, we introduced above additional notations: X N ?39? Thus, multiplying Eq ( 37 ) by rmn and summing the result with its complex conjugate, and taking into account Eq ( 39 ), we obtain the following expression: ?36? ?37? ?38? ?40? ?41? ? ?Omn ? Onm?rmnrnm ? cmnrnm ? cnmrmn ? iE??dmm dnn?? ??dnn ? dmm??rmnrnm iE?dmndmnrnm ? dnmdnmrmn? ? 2gmnjrmnj2 ? cmnrnm ? cnmrmn ? ?mn; m; n ? 1; N ; ?mn ? iE?dmndmnrnm ? dnmdnmrmn?; m; n ? 1; N : It is now seen from Eq ( 40 ) that the non-zero relaxation (?mn6?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: X N 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: ? XN Xm1 XN Xm1 B@2dmn XN ?iE m?2 n?1 XN Xm1 XN m?2 n?1 q6?m;n q?1 ? XN Xm1 ?mn? dmqrqnrnm dqnrmqrnm ? dnqrqmrmn dqmrnqrmn ? m; n; q 6? b; a; c; m; n; q 6? c; a; b ; m; n; q 6? c; b; a ?iE?dbcrcarab ?iE?dcbrbarac ?iE?dcarabrbc dcarbcrab ? dacrcbrba dbarcbrac ? dabrbcrca dabrcarbc ? dbaracrcb dcbracrba?? dbcrabrca?? dacrbarcb?; 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: m?2 n?1 ? XN Xm1 Now, let us perform the same analysis for the population differences: ?dmqrqmdmn dqmrmqdmn ? dnqrqndmn dqnrnqdmn? 2??ba ? ?ca ? ?cb?: ?49? 4 XN Xm1 ?mn? ?2iE XN Xm1 XN Here it is important that definition of ?mn. Eq ( 38 ) provides that: dab ? dbc ? dac: Taking into account Eq ( 48 ) along with the definition ( 41 ), the Eq ( 47 ) can be transformed into: ?47? ?48? ?50? 4 XN Xm1 ?mn? ?2iE XN Xm1 XN 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 N3 binations 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: of various triples that the second term contains, in these com4 XN Xm1 which is simplified further: ?51? ?52? ?53? ?54? ?55? ?56? ?57? 3.1 Approximation of the equations The model (16)?( 27 ) is approximated on two offset grids O = ?z ? ?t and O~ ? o~z domain G = [0, Lt] ? [Ll, Lr]. o~t in the ot ? tl ? l t; l ? 0; Nt ; t ? Lt ; Nt o~ t ? tl?0:5 ? ?l ? 0:5? t; l ? 0; Nt 1; t ? Lt ; Nt oz ? zj ? j h Ll; j ? 0; Nz ; h ? Lz ; Nz o~ z ? zj?0:5 ? ?j ? 0:5? h Ll; j ? 0; Nz 1; h ? Lz : Nz 2N XN Xm1 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. Here Nz and Nt are numbers of nodes along spatial coordinate and time, correspondingly. Let us also introduce the following notations for nodes: jp ? ??zp Ll?=h?; jc1 ? ??zc1 Ll?=h?; jc2 ? ??zc2 Ll?=h?: jL ? ??zL Ll?=h?; jR ? ??zR Ll?=h?; nl ? ?hl=hz?; 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.) The mesh function ?(zj) for the dielectric permittivity is defined on the grid ?z: ??zj? ? > > > > > > > > > > > > > > > : 8 > > > > > > > >>> ?J ; j ? ?jc1 ? ?J > > > > <> ?~J ; j ? ?jR ? ?J ?1 ? ?Jc1 ?=2; j ? jc1; 1?nl ? 1?; ?jc1 ? J nl 1?; J ? 1; Nl ; 1?nl ? 1?; ?jR ? J nl 1?; J ? 1; Nl ; ??J ? ?J?1?=2; j ? jc1 ? J nl; J ? 1; Nl ??~J ? ?~J?1??=2; j ? jR ? J nl; J ? 1; Nl 1; 1; ??~Jc2 ? 1?=2; j ? jc2: 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 O~ : H?zj?0:5; tl?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(zj, tl), E(zj, tl), ?mn(zj, ti), P(zj, tl), 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 t H?zj?0:5; tl?0:5? H?zj?0:5; tl 0:5? ? E?zj?1; tl? h E?zj; tl? ; ?58? ?59? ?61? ?62? In vacuum, the electric field strength is equal to its induction: and in the cover layers their relation is defined as: E?zj; tl? ? D?zj; tl?=??zj?; l ? 0; Nt ; j ? jc1; jL 1 [ jR ? 1; jc2 : 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 approximation of Eqs ( 20 )?( 23 ), based on the Crank-Nicolson scheme, is made in the following form: rmn?zj; tl? ? ?gmn ? iomn?r~mn ? iE~ X N q?1 E~ ? 0:5?E?zj; tl?1? ? E?zj; tl??; r~mn ? 0:5?rmn?zj; tl?1? ? rmn?zj; tl??; P?zj; tl?1? ? wE?zj; tl?1? ? k dqrrrq?zj; tl?1?; j ? jL; jR ; XN XN 3.2 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 = Nz we use the Cabaret approximation [42] for the translation operator: H?zNz 0:5; tl?0:5? H?zNz 0:5; tl 0:5? ? H?zNz 1:5; tl 0:5? H?zNz 1:5; tl 1:5? ? t H?zNz 0:5; tl 0:5? H?zNz 1:5; tl 0:5? h t ? 0; l ? 2; Nt 1; h D?z0; tl? ? 0; l ? 2; Nt 1: ?73? t t h h 3.3 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: E?zj; t0? ? E0 exp? ?zj hjp?2=tp2? cos?op ?zj hjp??; j ? 2; Nz 2; H?zj?0:5; t0:5? ? E0 exp? ?zj?05 hjp?2=tp2? cos?op ?zj?0:5 hjp??; j ? 2; Nz 3; ?75? E?zj; t0? ? 0; j ? 0; 1; Nz H?zj?0:5; t0:5? ? 0; j ? 0; 1; Nz 1; Nz; 2; Nz 1: Here jp is the mesh node index corresponding to the position of the initial pulse center: jp ? ?zp Ll?=h: Note that the initial values of H are defined at t0.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): E?z1; t1? ? E?zNz ; t1? ? H?zNz 0:5; t1:5? ? H?z0:5; t1:5? ? 0: The initial state of the medium is approximated in the following way: NXz1 j?0 ?2 2 3.4 Proof of finite-difference scheme conservativeness Theorem 3.1. The finite-difference scheme (58)?(80) is conservative one and preserves difference analogues of the invariants ( 28 )?( 31 ). Proof. The invariants can be obtained by performing the summation of Eqs (59) and (60) with taking into account the boundary conditions (70)?(73). After that, we get the conservation law in the form: ?H?zj?0:5; tNt 0:5? H?zj?0:5; t1:5?? ? ?E?zNz 1; tl? E?z1; tl??? h l?2 Applying a similar procedure to the Eq (60), we obtain: ?D?zj; tNt ? D?zj; t0?? ? ?H?zNz 0:5; tl?0:5? H?z0:5; tl?0:5??? If there are the zero-value boundary conditions for the electromagnetic field then Eqs (81) and (82) transform into the equalities: D?zj; tl?1? ? D?zj; t0? ? const; l ? 0; Nt 1; H?zj?0:5; tl?0:5? ? H?zj?0:5; t0:5? ? const; l ? 0; Nt 1: Theorem 3.2. A difference analogue of the density matrix evolution invariant ( 32 ) holds in each node of the grid O: h l?2 h l?2 NXz1 j?0 XNz j?0 XN Xm1 h l?2 ?D?z1; tl 1? D?z0; tl 1?? ?D?z1; tNt 2? D?z1; t0?? ?D?zNz 1; tl 1? D?zNz 2; tl 1?? ?D?zNz 1; tNt 2? D?zNz 2; t0??: ?82? ?83? ?84? ?85? ?86? ?87? ?88? Definitions of E~ and r~ are the same as in Eq (65). With these notations, the Eqs (63) and (64) are written as: Next, we write the squared derivatives of the non-diagonal elements in the following way: t rmn?zj; tl? ? rmn?zj; tl?1? rmn?zj; tl? r~mn: ? t Using this definition, we repeat steps ( 34 ) and ( 40 ) of the analytical derivation and obtain an analogue of Eq ( 43 ): ? 2gmnjr~mnj2 ? c~ mnr~nm ? c~ nmr~mn ? ?mn: ?89? Similarly, we obtain the analog of the Eq ( 42 ) for the d~mn: jdmn?zj; tl?1?j2 jdmn?zj; tl?j2 t ? 2d~mn X q6?m;n 3.5 Approximation order investigation Theorem 3.3. The finite-difference scheme (59)?(73) approximates the set of differential Eqs (16)?( 23 ) with the second order in the inner nodes of the mesh. Namely, the Eq (59) approximates the Eq (16) with respect to the point (zj+0.5, tl). The Eq (60) approximates the Eq ( 17 ) with respect to the point (zj, tl+0.5). The density matrix Eqs ( 21 )?( 23 ) are approximated by (63)?(67) with the second order with respect to the point (zj, tl+0.5). The theorem 3.3 is proved in a standard way with the help of the Taylor series representation. The boundary conditions are approximated by the Cabaret finite-difference scheme (70)? (73) with the second order in the points ?z1; tl 0:5?; ?z0:5; tl?; ?zNz 0:5; tl?; ?zNz 1; tl 0:5?; l ? 2; Nt 1:. 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). 3.6 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 17 / 30 used: rmn?zj; tl? ? The functions on zero-iteration are chosen equal to their values on the previous time layer: E0?zj; tl?1? ? E?zj; tl?; rq0m?zj; tl?1? ? rqm?zj; tl?; j ? jL; jR ; l ? 0; Nt 1: The iteration process is stopped if both of the following conditions are satisfied: jEs?1?zj; tl?1? Es?zj; tl?1?j < e1jEs?zj; tl?1?j ? e2; jrsm?n1?zj; tl?1? rsmn?zj; tl?1?j < e1jrsmn?zj; tl?1?j ? e2; where e1 and e2 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(zj+0.5, tl+0.5) on the next time layer, then we have D (tl+1, zj+1) from Eq (60) and Es+1(zj, tl+1) from Eq (94). We substitute the last value in the Eqs (92) and (93) at the next iteration step, etc. Theorem 3.4. Provided that h<const and zero-iterations are chosen as (96), the process of iterations (92)?(95) converges to the unique solution of the finite-difference scheme (59)?(69), which exists in a vicinity of the mesh functions defined by the zero iteration (96). Convergence rate is as a geometric progression with the denominator q h. 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 kh: where CE and C? are constants, and that the following inequalities: jjEsjjh CE; jjrsjjh Cr; jjEs?1 Esjjh qjjEs jjrs?1 rsjjh qjjrs Es 1jjh; rs 1jjh; q < 1 ?94? ?95? ?96? ?97? ?98? ?99? ?100? ?101? are valid. We use the difference norm C: jjf jjC ? max?jf ?zj?j?: zj ?102? 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 kEs+1 ? EskC and k?s+1 ? ?skC, and using estimations similar to those applied for the iteration modules we get the inequalities (102) in the norm (102). 3.7 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 p?? dielectric permittivity is greater than unity and u ? 1= ? 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. 4. 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. 4.1 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 E0 was chosen to be 1.05 106 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 105 W/cm2. Another wellknown formula [46] gives the value of the magnetic field induction and the magnetic field strength: I ? c?0E02: jB0j ? n c jE0j; r???? m ?0ojH0j ? njE0j; ?103? ?104? where n ? p????m?? is the refractive index. Taking into account that in a vacuum ? = ? = 1, we obtain: B0 = 3.5 10?3 T, H0 = 2785 A/m (or 27.85 A/cm). Evidently, for practice the intensity less than 106 W/cm2 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. The dipole moment unit d0 is equal [47] to 1.47 10?18 cm2.5g0.5s-1 or 4.58 10?30 C m. The density of the molecules in the medium is equal to M = 5.67 1021 cm-3 = 5.67 1027 m-3. In this case the coefficient ? = 12.5 in the Eq ( 15 ). We chose the mesh steps as h = ? = 0.01. Iteration process parameters are chosen as e1 = 10?8 and e1 = 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. 4.2 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. 4.3 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: XN0 k?1 Ereflntran?on? ? E?zreflntran; tk?exp? 2pLitk n? ; on ? 2Np0 n; tk ? NL0 k: ?105? Here L is the time interval (in this case, L = Lt, N0 = Nt = 104) or the spatial computational domain length (L = Lz, N0 = Nz = 3200). This transform is applied to the measured signals Erefl, Etran 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. 20 / 30 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 Fig 2. Pulses Erefl and Etran reflected from (solid line) and transmitted through (dashed line) the medium without cover (a), and with cover (b) under averaging over 32 random realizations of disordered structure (dimensionless units). 21 / 30 Fig 3. Comparison of the incident pulse spectra (Einc) 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). 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 (E0 = 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: XN0 i?1 ?ai corr?a; b? ? amean??bi bmean? ?ai amean? ,0vu???????????????????????????????v???????????????????????????????1 u XN0 uu XN0 @t 2t ?bi bmean?2A; ?106? i?1 i?1 22 / 30 Fig 4. Time evolution of the energy level populations (a) in a section corresponding to the beginning of the medium z = 0. The detailed view of the initial stage for the second and third energy level populations normalized to their maximum values (b) (dimensionless units). N0 i?1 amean ? ai; bmean ? N0 i?1 bi; The correlation of the energy level populations with the instantaneous electric field strength in this section of the medium gives us the following values: where, amean, bmean 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: corr?r11; r22? ? 1; corr?r22; r33? ? 0:781; corr?r33; r11? ? 0:71: ?108? corr?r11; Ejz?0? ? 0:001; corr?r33; Ejz?0? ? 0:14: ?109? 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: X3 Xm1 ?107? ?110? 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. ?111? ?112? 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(Zrefl\tran, tk) 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] with window function G?t? ? exp? ??t tI;II;III...?=yI;II;III...?KG?; which is a near-rectangular window. Therefore, we obtain: EI;II;III;...?on? ? j reflntran XN0 k?1 G?tk? E?zreflntran; tk? exp? 2pLitk n?j; on ? 2Np0 n; tk ? NL0 k; where EI,II,III. . .(?n) is the spectrum of the sub-pulse located in the corresponding time interval marked in Fig 2. Time tI,II,III. . . denotes a centre the corresponding time interval and ?I,II,III. . . denotes the window half-length. As before, L is the time interval, where the signal is analyzed, and N' is the number of time nodes in this time interval (L = Lt, N' = 104). The parameter KG in Eq (111) is even and characterizes the steepness of the window function at the ends of a Fig 5. The first sub-pulse (area I in 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). 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 Erefl. 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 25 / 30 Fig 6. Spectrum of sub-pulses reflected from and transmitted through the covered medium and belonging to the time interval I (a), II (b), III (c), IV (d) in Fig 2B. The reflected sub-pulse spectra in the vicinity of the frequency 1.75 (emission frequency) corresponding to the energy level transition 3?1 (e). 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 26 / 30 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. 5. 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 27 / 30 spectral intensities of these frequencies. Therefore, they can be used with success for the substance detection and identification. Author Contributions Conceptualization: Vyacheslav A. Trofimov. Funding acquisition: Vyacheslav A. Trofimov. Investigation: Svetlana A. Varentsova, Irina G. Zakharova, Dmitry Yu. Zagursky. Methodology: Vyacheslav A. Trofimov. Software: Dmitry Yu. Zagursky. Supervision: Vyacheslav A. Trofimov. Visualization: Svetlana A. Varentsova, Dmitry Yu. Zagursky. Writing ? original draft: Irina G. Zakharova, Dmitry Yu. Zagursky. Writing ? review & editing: Svetlana A. Varentsova, Irina G. Zakharova. 28 / 30 16. Duling I, Zimdars D. Terahertz imaging: Revealing hidden defects. Nature Photonics. 2009; 3( 11 ): 630?632. 29 / 30 1. Federici JF , Schulkin B , Huang F , Gary D , Barat R , Oliveira F , et al. THz imaging and sensing for security applications?explosives, weapons and drugs . Semiconductor Science and Technology . 2005 ; 20 : S266? S280 . 2. Lu M , Shen J , Li N , Zhang Y , Zhang C , Liang L , et al. Detection and identification of illicit drugs using terahertz imaging . Journal of Applied Physics . 2006 ; 100 ( 10 ): 10310 . 3. Chen J , Chen Yu , Zhao H , Bastiaans GJ , Zhang X -C. Absorption coefficients of selected explosives and related compounds in the range of 0.1?2.8 THz . Optics. Express. 2007 ; 15 ( 19 ): 12060 . PMID: 19547570 4. Leahy-Hoppa MR , Fitch MJ , Zheng X , Hayden LM , Osiander R . Wideband terahertz spectroscopy of explosives . Chemical Physics Letters . 2007 ; 434 : 227 ? 230 . 5. Choi K , Hong T , Sim KI , Ha T , Park BC , Chung JH , et al. Reflection terahertz time-domain spectroscopy of RDX and HMX explosives . Journal of Applied Physics . 2014 ; 115 ( 2 ): 023105 . 6. Nickel DV , Ruggiero MT , Korter TM , Mittleman DM . Terahertz disorder-localized rotational modes and lattice vibrational modes in the orientationally-disordered and ordered phases of camphor . Physical Chemistry Chemical Physics . 2015 ; 17 ( 10 ): 6671 ? 7078 . 7. Van Rheenen AD , Haakestad MW . Detection and identification of explosives hidden under barrier materials?what are the THz-technology challenges ? Proc. SPIE . 2011 ; 8017 : 801719 . 8. Ergu?n S , So?nmez S. Terahertz technology for military applications . Journal of Military and Information Science . 2015 ; 3 ( 1 ): 13 ? 16 . 9. Davies AG , Burnett A D , Fan W , Linfield EH , Cunningham JE . Terahertz spectroscopy of explosives and drugs . Materials Today . 2008 ; 11 ( 3 ): 18 ? 26 . 10. Katz G , Zybin S , Goddard WA , Zeiri Y , Kosloff R. Direct MD simulations of terahertz absorption and 2D spectroscopy applied to explosive crystals . J. Phys. Chem . Lett. 2014 ; 5 ( 5 ): 772 ? 776 . https://doi.org/ 10.1021/jz402801m PMID: 26274066 11. Tanno T , Umeno K , Ide E , Katsumata I , Fujiwara K Ogawa N. Terahertz spectra of 1-cyanoadamantane in the orientationally ordered and disordered phases . Philosophical Magazine Letters . 2014 ; 94 ( 1 ): 25 ? 29 . 12. McIntosh AI , Yang B , Goldup SM , Watkinson M , Donnan RS . Terahertz spectroscopy: a powerful new tool for the chemical sciences? Chem. Soc. Rev . 2012 ; 41 : 2072 ? 2082 . https://doi.org/10.1039/ c1cs15277g PMID: 22143259 13. Parrott E , Zeitler JA , Friscic T , Pepper M , Jones W , Day GM , Gladden LF . Testing the sensitivity of terahertz spectroscopy to changes in molecular and supramolecular structure: a study of structurally similar cocrystals . Crystal Growth and Design . 2009 ; 9 ( 3 ): 1452 ? 1460 . 14. Li R , Zeitler JA , Tomerini D , Parrott E , Gladden L , Day G. A study into the effect of subtle structural details and disorder on the terahertz spectrum f crystalline benzoic acid . Phys. Chem. Chem. Phys . 2010 ; 12 : 5329 ? 5340 . PMID: 21491656 15. Kojima S , Shibata T , Igawa H , Mori T. Broadband terahertz time-domain spectroscopy: crystalline and glassy drug materials . IOP Conf. Series: Materials Science and Engineering . 2014 ; 54 : 012001 . 17. Kawase K , Shibuya T , Hayashi SI , Suizu K. THz imaging techniques for nondestructive inspections . Comptes Rendus Physique . 2010 ; 11 ( 7 ): 510 ? 518 . 18. Kraus E , Kremling S , Baudrit B , Heidemeyer P , Bastian M , Stoyanov OV , Starostina IA . Reflection time domain terahertz system for testing polymeric connections . Polymers Research Journal . 2015 ; 9 ( 3 ): 345 . 19. Ouchi T , Kajiki K , Koizumi T , Itsuji T , Koyama Y , Sekiguchi R . et al. Terahertz imaging system for medical applications and related high efficiency terahertz devices . Journal of Infrared , Millimeter, and Terahertz Waves . 2014 ; 35 ( 1 ): 118 ? 130 . 20. Ok G , Park K , Kim HJ , Chun HS , Choi SW . High-speed terahertz imaging toward food quality inspection . Applied optics . 2014 ; 53 ( 7 ): 1406 ? 1412 . https://doi.org/10.1364/AO.53.001406 PMID: 24663370 21. Shigekawa H , Yoshida S , Takeuchi O . Spectroscopy: Nanoscale terahertz spectroscopy . Nature Photon . 2014 ; 8 : 815 ? 817 . 22. Kemp MC . Explosives detection by terahertz spectroscopy?a bridge too far? Terahertz Science and Technology , IEEE Transactions on. 2011 ; 1 ( 1 ): 282 ? 292 . 23. Palka N. Identification of concealed materials, including explosives, by terahertz reflection spectroscopy . Opt. Eng . 2013 ; 53 ( 3 ): 031202 . 24. Ortolani M , Lee JS , Schade U , Hu?bers HW . Surface roughness effects on the terahertz reflectance of pure explosive materials . Appl. Phys. Lett . 2008 ; 93 ( 8 ): 08190 . 25. Schecklman S , Zurk LM , Henry S , Kniffin GP . Terahertz material detection from diffuse surface scattering . J. Appl. Phys . 2011 ; 109 ( 9 ): 094902 . 26. Choi J , Ryu SY , Kwon WS , Kim KS , Kim S. Compound explosives detection and component analysis via terahertz time-domain spectroscopy . Journal of the Optical Society of Korea . 2013 ; 17 ( 5 ): 454 ? 460 . 27. Withayachumnankul W , Fischer BM , Abbott D. Numerical removal of water vapor effects from terahertz time-domain spectroscopy measurements . Proceedings of the Royal Society of London A: Mathematical , Physical and Engineering Sciences . 2008 ; 464 : 2435 ? 2456 . 28. Puc U , Abina A , Rutar M , Zidans?ek A , Jegli? A , Valus?is G . Terahertz spectroscopic identification of explosive and drug simulants concealed by various hiding techniques . Applied optics . 2015 ; 54 ( 14 ): 4495 ? 4502 . https://doi.org/10.1364/AO.54.004495 PMID: 25967507 29. Trofimov VA , Varentsova SA . An effective method for substance detection using the broad spectrum THz signal: a ?terahertz nose? . Sensors . 2015 ; 15 ( 6 ): 12103 ? 12132 . https://doi.org/10.3390/ s150612103 PMID: 26020281 30. Trofimov VA , Varentsova SA . Essential limitations of the standard THz TDS method for substance detection and identification and a way of overcoming them . Sensors . 2016 ; 16 ( 4 ): 502 . 31. Trofimov VA , Varentsova SA . False detection of dangerous and neutral substances in commonly used materials by means of the standard THz time domain spectroscopy . Journal of the European Optical Society . 2016 ; 11 : 16016 . 32. Trofimov V , Zagursky D , Zakharova I. Propagation of laser pulse with a few cycle duration in multi-level media . Days on Diffraction (DD) , IEEE, 2015 ; 1 ? 5 . 33. Imura K , Nagahara T , Okamoto H . Near-field two-photon-induced photoluminescence from single gold nanorods and imaging of plasmon modes . J. Phys. Chem. B . 2005 ; 109 : 13214 ? 13220 . https://doi.org/ 10.1021/jp051631o PMID: 16852648 34. Marskar R , Osterberg U. Multilevel Maxwell-Bloch simulations in inhomogeneously broadened media . Optics Express . 2011 ; 19 : 16784 ? 16796 . https://doi.org/10.1364/OE.19.016784 PMID: 21935040 35. Castin Y , Molmer K . Maxwell-Bloch equations: A unified view of nonlinear optics and nonlinear atom optics . Phys. Review A . 1995 ; 51 ( 5 ): 3426 ? 3428 . 36. Taflove A. , Hagness CS . Computational electrodynamics: the finite-difference time-domain method: Artech house; 2005 . 37. Liu Y. Fourier analysis of numerical algorithms for the Maxwell equations . Journal of Computational Physics . 1996 ; 124 : 396 ? 416 . 38. Ziolkowski RW . The incorporation of microscopic material models into the FDTD approach for ultrafast optical pulse simulations . IEEE Transactions on Antennas and Propagation . 1997 ; 45 : 375 ? 391 . 39. Bidegaray B , Bourgeade A , Reignier D . Introducing Physical Relaxation Terms in Bloch Equations . Journal of Computational Physics . 2001 ; 170 : 603 ? 613 . 40. Bidegaray B . Time discretizations for Maxwell-Bloch equations . Numerical Methods for Partial Differential Equations . 2003 ; 19 : 284 ? 300 . 41. Born M. , Wolf E . Principles of optics: electromagnetic theory of propagation, interference and diffraction of light: Elsevier; 2013 . 42. Goloviznin VM . ; Samarskii AA . Finite approximation of convective transport with a space splitting time derivative . Matematicheskoe Modelirovanie (in Russian) . 1998 ; 10 ( 1 ): 86 ? 100 . 43. Boyd RW . Nonlinear optics: Elsevier; 2003 . 44. Thompson A , Taylor BN . ( 2008 ). Guide for the Use of the International System of Units (SI). Special Publication (NIST SP)-811: Ambler Thompson and Barry N. Taylor ; 2008 . 45. Yee K. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media . IEEE Transactions on Antennas and Propagation . 1966 ; 14 : 302 ? 307 . 46. Bakhvalov NS . Courant?Friedrichs?Lewy condition , in Hazewinkel M. Encyclopedia of Mathematics: Springer Science+Business Media B.V. Kluwer Academic Publishers; 2001 . 47. Frej H , Waclawek W , Kyziol JB . The dipole moments of carbazole derivatives . Chem. Papers . 1988 ; 42 ( 2 ): 213 ? 221 . 48. Chui K . Introduction to Wavelets: Academic; New York, 1992 . 49. Cohen L. Time-frequency distributions?a review . Proceedings of the IEEE , 1989 ; 77 ( 7 ): 941 ? 981 .


This is a preview of a remote PDF: https://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0201572&type=printable

Vyacheslav A. Trofimov, Svetlana A. Varentsova, Irina G. Zakharova, Dmitry Yu. Zagursky. 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, PLOS ONE, 2018, DOI: 10.1371/journal.pone.0201572