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
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.
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 , food
quality inspection . Currently, nanoscale terahertz spectroscopy is being actively developed
as well .
Despite the obvious advantages of THz TDS using for the substance detection and
identification, there is some shortcoming of its application . 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)
. 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
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 . 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 .
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 , .
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
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  with the step-split method , . 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 . However, the Crank-Nicholson scheme does not possess the
positiveness property of the density matrix diagonal elements for three and more energy levels.
In  it was shown that negative values of the energy levels populations may even occur in
numerical simulation. In the same paper  and later in  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 . Using
the Cabaret scheme , 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
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
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?;
Lt; Ll < z
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 (
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
. Below we consider the case of a non-magnetic medium, thus, in Eq (
) 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
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,
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 (
) is replaced by the
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:
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
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 (
) to a dimensionless form.
We use the following normalization, where the index ?dim? means dimensionless variables:
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 :
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 (
), Eqs (
) can be rewritten in the dimensionless form:
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;
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 (
along z and t coordinates.
After the integration, the following invariants can be written:
D?z; t?dz ? I1
H?z; t?dz ? I2
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
) and (
) are transformed:
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 (
) and (
) obeys the following invariant:
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:
onm; gmn ? gnm; gmm ? 0; m; n ? 1; N :
dmndmn?; m; n ? 1; N :
For brevity, we introduced above additional notations:
Thus, multiplying Eq (
) by rmn and summing the result with its complex conjugate, and
taking into account Eq (
), we obtain the following expression:
? ?Omn ? Onm?rmnrnm ? cmnrnm ? cnmrmn ? iE??dmm
iE?dmndmnrnm ? dnmdnmrmn? ?
2gmnjrmnj2 ? cmnrnm ? cnmrmn ? ?mn; m; n ? 1; N ;
iE?dmndmnrnm ? dnmdnmrmn?; m; n ? 1; N :
It is now seen from Eq (
) that the non-zero relaxation (?mn6?0) causes exponential decay
of each of terms in the second sum in the formula (
). Therefore, to derive the invariant, we
assume below that ?mn = 0. In general case, one can write a law of density matrix elements
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 (
) in the
Now, let us perform the summation of Eqs (
) at zero relaxation (?mn = 0) using Eq (
and also the summation of Eq (
) for all the pairs of indices m and n such that 1 n<m N:
XN Xm1 B@2dmn XN
XN Xm1 XN
m?2 n?1 q6?m;n
dqnrmqrnm ? dnqrqmrmn
m; n; q 6? b; a; c; m; n; q 6? c; a; b ; m; n; q 6? c; b; a
dcarbcrab ? dacrcbrba
dbarcbrac ? dabrbcrca
dabrcarbc ? dbaracrcb
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 (
) contains all such combinations. Thus, the second term in
) is equal to zero as well. Hence, the all-pairs sum of the squared off-diagonal elements is
Now, let us perform the same analysis for the population differences:
dqmrmqdmn ? dnqrqndmn
2??ba ? ?ca ? ?cb?: ?49?
XN Xm1 XN
Here it is important that definition of ?mn. Eq (
) provides that:
dab ? dbc ? dac:
Taking into account Eq (
) along with the definition (
), the Eq (
) can be transformed
XN Xm1 XN
We see that our analysis demonstrates the appearance of an additional term ?2(?ba + ?ca +
?cb) in Eq (
), 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 (
) can be transformed into:
of various triples that the second term contains, in these
which is simplified further:
3.1 Approximation of the equations
The model (16)?(
) is approximated on two offset grids O = ?z ? ?t and O~ ? o~z
domain G = [0, Lt] ? [Ll, Lr].
o~t in the
tl ? l t; l ? 0; Nt ; t ? Lt ;
o~ t ?
tl?0:5 ? ?l ? 0:5? t; l ? 0; Nt
1; t ? Lt ;
zj ? j h
Ll; j ? 0; Nz ; h ? Lz ;
o~ z ?
zj?0:5 ? ?j ? 0:5? h
Ll; j ? 0; Nz
1; h ? Lz :
To solve Eqs (16)?(
), we use a combination of the well-known Yee's scheme , 
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
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
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? ? >
>>> ?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
??~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
Eqs (16) and (
) are approximated using the Yee's method. Therefore, in our notations the
corresponding finite-difference scheme is written as
H?zj?0:5; tl 0:5? ?
E?zj; tl? ;
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 (
), based on the Crank-Nicolson scheme, is made in the
rmn?zj; tl? ? ?gmn ? iomn?r~mn ? iE~
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 ;
3.2 Boundary conditions approximation
At the edges of the computational domain, Eqs (
) and (
) 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  for the translation
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? ?
H?zNz 0:5; tl 0:5?
H?zNz 1:5; tl 0:5?
? 0; l ? 2; Nt
D?z0; tl? ? 0; l ? 2; Nt
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
H?zj?0:5; t0:5? ? E0 exp? ?zj?05
hjp?2=tp2? cos?op ?zj?0:5
hjp??; j ? 2; Nz
E?zj; t0? ? 0; j ? 0; 1; Nz
H?zj?0:5; t0:5? ? 0; j ? 0; 1; Nz
Here jp is the mesh node index corresponding to the position of the initial pulse center:
jp ? ?zp
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:
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 (
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?
Applying a similar procedure to the Eq (60), we obtain:
?D?zj; tNt ?
D?zj; t0?? ?
?H?zNz 0: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
H?zj?0:5; tl?0:5? ?
H?zj?0:5; t0:5? ? const; l ? 0; Nt
Theorem 3.2. A difference analogue of the density matrix evolution invariant (
) holds in
each node of the grid O:
?D?z1; tl 1?
D?z0; tl 1??
?D?z1; tNt 2?
?D?zNz 1; tl 1?
D?zNz 2; tl 1??
?D?zNz 1; tNt 2?
D?zNz 2; t0??:
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:
rmn?zj; tl? ? rmn?zj; tl?1?
rmn?zj; tl? r~mn:
Using this definition, we repeat steps (
) and (
) of the analytical derivation and obtain
an analogue of Eq (
? 2gmnjr~mnj2 ? c~ mnr~nm ? c~ nmr~mn ? ?mn:
Similarly, we obtain the analog of the Eq (
) for the d~mn:
3.5 Approximation order investigation
Theorem 3.3. The finite-difference scheme (59)?(73) approximates the set of differential Eqs
) 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 (
respect to the point (zj, tl+0.5). The density matrix Eqs (
) 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
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
The iteration process is stopped if both of the following conditions are satisfied:
Es?zj; tl?1?j < e1jEs?zj; tl?1?j ? e2;
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:
rs 1jjh; q < 1
are valid. We use the difference norm C:
jjf jjC ? max?jf ?zj?j?:
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 , 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 . 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= ? 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  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  gives the value of the magnetic field induction and the magnetic field
I ? c?0E02:
?0ojH0j ? njE0j;
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  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 (
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:
2pLitk n? ; on ? 2Np0 n; tk ? NL0 k:
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
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 . To eliminate these mistakes, the sub-pulses should be
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)
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  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:
corr?a; b? ?
u XN0 uu XN0
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).
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
corr?r11; r22? ?
1; corr?r22; r33? ? 0:781; corr?r33; r11? ?
corr?r11; Ejz?0? ? 0:001; corr?r33; Ejz?0? ?
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 (
), 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:
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.
In our previous papers ,  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  is applied to the signal E(Zrefl\tran, tk) in the corresponding time
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 , 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
 with window function
G?t? ? exp? ??t
which is a near-rectangular window. Therefore, we obtain:
EI;II;III;...?on? ? j
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
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
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.
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
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.
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(
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 .