# Experimental investigation of electric field distributions in a chaotic 3D
microwave rough billiard

Oleg Tymoshchuk, Nazar Savytskyy, Oleh Hul, Szymon Bauch, and Leszek Sirko
Institute of Physics, Polish Academy of Sciences, Aleja Lotników 32/46, 02-668
Warszawa, Poland

(January 28, 2007)

###### Abstract

We present the first experimental study of the electric field distributions
$E_{N}$ of a three-dimensional (3D) microwave chaotic rough billiard with the
translational symmetry. The translational symmetry means that the cross-
section of the billiard is invariant under translation along $z$ direction.
The 3D electric field distributions were measured up to the level number
$N=489$. In this way the experimental spatial correlation functions
$C_{N,p}({\bf x,s})\propto\langle E_{N,p}({\bf x}+\frac{1}{2}{\bf
s})E_{N,p}^{\ast}({\bf x}-\frac{1}{2}{\bf s})\rangle$ were found and compared
with the theoretical ones. The experimental results for higher two-dimensional
level number $N_{\bot}$ appeared to be in good agreement with the theoretical
predictions.

###### pacs:

05.45.Mt,05.45.Jn

In this paper we present the first experimental investigation of electric
field distributions of the chaotic 3D microwave rough billiard with the
translational symmetry. Due to experimental difficulties there are very few
experimental studies devoted to 3D chaotic microwave cavities Sirko1995 ;
Alt1997 ; Dorr1998 ; Eckhardt1999 ; Dembowski2002 . In a pioneering experiment
Deus et al. Sirko1995  have been measured eigenfrequencies of the 3D chaotic
(irregular) microwave cavity in order to confirm that their distribution
displays behavior characteristic for classically chaotic quantum systems,
viz., the Wigner distribution. Three-dimensional chaotic cavities as well as
properties of random electromagnetic vector field have been also scarcely
studied theoretically Primack2000 ; Prosen1997 ; Arnaut2006 .

In general, there is no analogy between quantum billiards and electromagnetic
cavities in three dimensions. However, for 3D cavities with the translational
symmetry the classification of the modes into transverse electric (TE) and
transverse magnetic (TM) is possible. The TM modes are especially important
because they allow for the simulation of 2D quantum billiards on cross-
sectional planes of 3D cavities. Furthermore, we show in this paper that the
distributions of the electric field of TM modes of the 3D chaotic rough cavity
can be experimentally measured.

Figure 1: Upper panel: Sketch of the chaotic half-circular 3D microwave rough
billiard in the $xy$ plane. Dimensions are given in cm. The cavity sidewalls
are marked by 1 and 2 (see text). Squared wave functions
$|\psi_{N,p}(R_{c},\theta)|^{2}$ were evaluated on a half-circle of fixed
radius $R_{c}=9.25$ cm. Billiard’s rough boundary $\Gamma$ is marked with the
bold line. Lower panel: White circle marks the position of the hole drilled in
the upper wall of the cavity. This hole was used to introduce the perturber
inside the cavity in order to measure the $z$-component of the electric field
distributions $E_{N,p}({\bf x})$.

In the experiment we used 3D cavity with the translational symmetry in the
shape of a rough half-circle (Fig. 1) with the height $h=60$ mm. The cavity
was made of polished aluminium.

We suppose that the direction of the translational symmetry of the cavity is
along the $z$-axis. The boundary conditions at $z=0$ and $z=h$ demand that the
$z$ dependence of the $z$-component of the electric and magnetic fields
$E_{N,p}({\bf x})$ and $B_{N,p}({\bf x})$ of TM modes be in the form
$E_{N,p}({\bf x})\equiv E_{N,p}(x,y,z)=A_{N,p}\psi_{N,p}(x,y)f_{p}(z)$, where
$f_{p}(z)=\cos(p\pi z/h)$, $p=0,1,2\ldots$, $A_{N,p}$ is the normalization
constant and $B_{N,p}({\bf x})=0$. The functional dependence of $E_{N,p}({\bf
x})$ on the plane cross section coordinates is denoted by the amplitude
$\psi_{N,p}(x,y)\equiv E_{N,p}(x,y)$. The amplitude $\psi_{N,p}(x,y)$
satisfies the Helmholtz equation

$(\bigtriangleup_{\bot}+k_{N,p}^{2})\psi_{N,p}(x,y)=0,$ $None$

where $\bigtriangleup_{\bot}$ is two-dimensional Laplacian operator and
$k_{N,p}=(k_{N}^{2}-(p\pi/h)^{2})^{1/2}$ is the effective wave vector. The
wave vector $k_{N}=2\pi\nu_{N}/c$, where $\nu_{N}$ is the resonance frequency
of the level $N$ and $c$ is the speed of light in the vacuum. The equation (1)
is equivalent to the Schrödinger equation (in units $\hbar=1$) describing a
particle of mass $m=1/2$ with the kinetic energy $k_{N}^{2}$ in an external
potential $V=(p\pi/h)^{2}$ Kim2005 . Therefore, microwave 3D cavities with the
translational symmetry simulate on the cross-sectional planes quantum
billiards with the external potential $(p\pi/h)^{2}$. In this way microwave
cavities can be effectively used beyond the standard 2D frequency limit (the
case $p=0$) Hans  in simulation of quantum systems. The amplitude
$\psi_{N,p}(x,y)$ fulfills Dirichlet boundary conditions on the sidewalls of
the billiard. Therefore, throughout the text the amplitudes $\psi_{N}(x,y)$
are also often called the wave functions $\psi_{N}(x,y)$. It is important to
note that the full electric field $E_{N,p}({\bf x})$ satisfies additionally
Neumann boundary conditions at the top and the bottom of the cavity.

Because of the relatively low quality factor of the cavity ($Q\simeq 4000$)
the value of the level number $N$ was evaluated from the Balian–Bloch formula
Balian

$N(k)=\frac{1}{3\pi^{2}}Vk^{3}-\frac{2}{3\pi^{2}}\int_{S}\frac{d\sigma_{\omega}}{R_{\omega}}k,$
$None$

where k is the wave vector, $V=(9.43\pm 0.01)\cdot 10^{-4}$ m3 is the volume
of the cavity and $\int_{S}\frac{d\sigma_{\omega}}{R_{\omega}}=0.932$ m $\pm
0.005$ m is the surface curvature averaged over the surface of the cavity.

The measurements allowed us for the first time to evaluate the spatial
correlation function Kaufman88

$C_{N,p}({\bf x},{\bf s})=\frac{1}{\langle|E_{N,p}({\bf
x})|^{2}\rangle}\langle E_{N,p}({\bf x}+\frac{1}{2}{\bf s})E_{N,p}^{\ast}({\bf
x}-\frac{1}{2}{\bf s})\rangle,$ $None$

where the local average $\langle\cdots\rangle$ is defined as follows

$\langle|E_{N,p}({\bf
x})|^{2}\rangle=\frac{1}{\Delta^{n}}\int_{-\Delta/2}^{\Delta/2}|E_{N,p}({\bf
x}+{\bf s})|^{2}d^{n}s.$ $None$

The 3D cavity sidewalls are made of 2 segments (see Fig. 1). The rough segment
1 is described on the cross-sectional planes by the radius function
$R(\theta)=R_{0}+\sum_{m=2}^{M}{a_{m}\sin(m\theta+\phi_{m})}$, where the mean
radius $R_{0}$=10.0 cm, $M=20$, $a_{m}$ and $\phi_{m}$ are uniformly
distributed on [0.084,0.091] cm and [0,2$\pi$], respectively, and
$0\leq\theta<{\pi}$. (Here, for the convenience, the polar coordinates $r$ and
$\theta$ are used instead of the Cartesian ones $x$ and $y$.)

The surface roughness of a billiard on the cross-sectional planes is
characterized by the function $k(\theta)=(dR/d\theta)/R_{0}$. For our billiard
we have the angle average
$\tilde{k}=(\left<k^{2}(\theta)\right>_{\theta})^{1/2}\simeq 0.400$. In such a
billiard the classical dynamics is diffusive in orbital momentum due to
collisions with the rough boundary because $\tilde{k}$ is much above the chaos
border $k_{c}=M^{-5/2}=0.00056$ Frahm97 . The roughness parameter $\tilde{k}$
determines also other properties of the billiard Frahm  on the cross-sectional
planes. The amplitudes $\psi_{N,p}(r,\theta)$ are localized for the two-
dimensional level number $N_{\bot}<N_{e}=1/128\tilde{k}^{4}$, where
$N_{\bot}=\frac{A}{4\pi}k_{N,p}^{2}-\frac{P}{4\pi}k_{N,p}$. $A=(1.572\pm
0.002)\cdot 10^{-2}$ m2 and $P=0.537$ m $\pm 0.001$ m are the cross-sectional
plane area and its perimeter, respectively. Because of a large value of the
roughness parameter $\tilde{k}$ the localization border lies very low,
$N_{e}\simeq 1$. The border of Breit-Wigner regime is
$N_{W}=M^{2}/48\tilde{k}^{2}\simeq 52$. It means that between
$N_{e}<N_{\bot}<N_{W}$ Wigner ergodicity Frahm  ought to be observed and for
$N_{\bot}>N_{W}$ Shnirelman ergodicity should emerge.

To measure the amplitudes $\psi_{N,p}(r,\theta)$ of the 3D electric field
distributions we used a very effective method described in Savytskyy2003 . It
is based on the perturbation technique Slater52  and preparation of the “trial
functions” Savytskyy2004 ; Hul2005 ; Hul2006 . In the perturbation method a
small perturber is introduced inside the cavity to alter its resonant
frequencies and in this way to evaluate the squared wave functions
$|\psi_{N,p}(R_{c},\theta)|^{2}$ (see Fig. 1). The perturber (4.0 mm in length
and 0.3 mm in diameter, oriented in $z$-direction) was moved by the stepper
motor via the Kevlar line hidden in the groove (0.4 mm wide, 1.0 mm deep) made
in the cavity’s bottom wall along the half-circle $R_{c}$. The measurements
were performed at 0.36 mm steps along a half-circle with fixed radius
$R_{c}=9.25$ cm.

In order to find the dependence of the electric field distributions
$E_{N,p}({\bf x})$ on the $z$ coordinate and to estimate the wave vector
$k_{3}=p\pi/h$ we measured the electric field inside the 3D cavity along the
$z$-axis. Also in this case the perturber (4.5 mm in length and 0.3 mm in
diameter) was attached to the Kevlar line and moved by the stepper motor. The
perturber entered and exited the cavity by small holes (0.4 mm) drilled in the
upper and the bottom walls of the cavity. The both holes were located at the
position: $r=9.11$ cm, $\theta=0.47$ radians.

Using the method of the “trial wave function” we were able to reconstruct 75
experimental wave functions $\psi_{N,p}(r,\theta)$, which belonged to TM modes
of the rough half-circular 3D billiard with the level number $N$ between 2 and
489. The range of corresponding eigenfrequencies was from $\nu_{2}\simeq 2.47$
GHz to $\nu_{489}\simeq 11.99$ GHz. The remaining wave functions belonging to
TM modes, from the range $N=2-489$, were not reconstructed because of near-
degeneration of the neighboring eigenfrequencies or due to the problems with
the measurements of $|\psi_{N,p}(R_{c},\theta)|^{2}$ along a half-circle
coinciding for its significant part with one of the nodal lines of
$\psi_{N,p}(r,\theta)$.

Figure 2:  The reconstructed wave function $\psi_{460,0}(r,\theta)$ of the
chaotic half-circular microwave rough billiard. The amplitudes have been
converted into a grey scale with white corresponding to large positive and
black corresponding to large negative values, respectively. Dimensions of the
billiard are given in cm. In the figure the $z$ dependence of the electric
field distribution
$E_{460,0}(r,\theta,z)|_{x=0}\propto\psi_{460,0}(r,\theta)|_{x=0}f_{0}(z)$ is
also shown.

In Fig. 2 and Fig. 3 we show two examples of reconstructed wave functions
$\psi_{460,0}(r,\theta)$ and $\psi_{463,4}(r,\theta)$, respectively. The
character of the wave functions predominantly depends on the effective wave
vector $k_{N,p}$. It is seen that the wave function $\psi_{463,4}(r,\theta)$
in Fig. 3 is more regular than the one presented in Fig. 2 in spite of having
the larger level number $N$.

Figure 3:  The reconstructed wave function $\psi_{463,4}(r,\theta)$ of the
chaotic half-circular microwave rough billiard. The $z$ dependence of the
electric field distribution
$E_{463,4}(r,\theta,z)|_{x=0}\propto\psi_{463,4}(r,\theta)|_{x=0}f_{4}(z)$ is
also shown.

In order to check ergodicity of the billiard’s wave functions
$\psi_{N,p}(R_{c},\theta)$, especially close to the ergodicity borders, one
should use some additional measures such as e.g., calculation of the
structures of their energy surfaces Frahm97 . For this reason we extracted
wave function amplitudes $C^{(N,p)}_{nl}=\left<n,l|N,p\right>$ in the basis
$n,l$ of a half-circular billiard with radius $r_{max}$, where $n=1,2,3\ldots$
enumerates the zeros of the Bessel functions and $l=1,2,3\ldots$ is the
angular quantum number. As expected, close to the border of the regimes of
Breit-Wigner and Shnirelman ergodicity the wave function
$\psi_{460,0}(r,\theta)$ ($N_{\bot}=65$) was found to be extended
homogeneously over the whole energy surface Hlushchuk01  (figure not shown
here). In contrary, the wave function $\psi_{463,4}(r,\theta)$, $N_{\bot}=16$,
which lies closer to the localization boarder, displays the tendency to
localization in $n,l$ basis (figure not shown here).

The measurement of 3D electric field distributions $E_{N,p}({\bf x})$ allowed
us for the first time to find the experimental spatial correlation function
$C_{N,p}({\bf x,s})$. It is easy to show Berry77  that for the 3D chaotic
cavity with the translational symmetry the spatial correlation function should
have the following form:

$C_{N,p}({\bf x},|{\bf s}|)\equiv C_{N,p}(|{\bf
s}|)=J_{0}(k_{N,p}s_{xy})\cos(p\pi s_{z}/h),$ $None$

where $|{\bf s}|=(s_{xy}^{2}+s_{z}^{2})^{1/2}$. For the cross-sectional planes
$z=const$ the correlation function $C_{N,p}(|{\bf s}|)\sim
J_{0}(k_{N,p}s_{xy})$ is reduced to the well known result of Berry Berry77
for chaotic 2D wave functions described by a random superposition of plane
waves.

Figure 4: Panels (a)-(c) show the experimental correlation function
$C_{460,0}({\bf x},|{\bf s}|)$ calculated at ${\bf x}=$ (-2.75 cm, 4.35 cm, 0
cm) for the three projection angles $\phi=0$, $\pi/4$, and $\pi/2$,
respectively. Experimental correlation function $C_{460,0}({\bf x},|{\bf s}|)$
(full line) is compared with the theoretical one (dashed line).

In Fig. 4(a)-(c) we show a representative example of the experimental
correlation function $C_{460,0}({\bf x},|{\bf s}|)$ ($N_{\bot}=65$) calculated
at ${\bf x}=$ (-2.75 cm, 4.35 cm, 0 cm) for the three different projection
angles $\phi=0$, $\pi/4$, and $\pi/2$, respectively, where
$\phi=\arcsin(s_{z}/|{\bf s}|)$. The local average $\langle\cdots\rangle$
required for the evaluation of $C_{N,p}({\bf x},|{\bf s}|)$ (see the formulas
(3) and (4)) was calculated on the cross-sectional plane $xy$ in the range
$\Delta/2=2\pi/k_{N,p}$. The experimental correlation functions
$C_{460,0}({\bf x},|{\bf s}|)$ are compared in Fig. 4 with the theoretical
ones. In all cases we find good agreement with the theoretical predictions
given by the formula (5). Small discrepancies observed in Fig. 4(a) for $|{\bf
s}|>1$ can be connected with the finiteness of the system and were
theoretically studied in Baecker02 .

Figure 5: Panels (a)-(c) show the experimental correlation function
$C_{463,4}({\bf x},|{\bf s}|)$ calculated at ${\bf x}=$ (-2.75 cm, 4.35 cm, 0
cm) for the three projection angles $\phi=0$, $\pi/4$, and $\pi/2$,
respectively. Experimental correlation function $C_{463,4}({\bf x},|{\bf s}|)$
(full line) is compared with the theoretical one (dashed line).

Fig. 5(a)-(c) shows the experimental correlation function $C_{463,4}({\bf
x},|{\bf s}|)$ ($N_{\bot}=16$) calculated at ${\bf x}=$ (-2.75 cm, 4.35 cm, 0
cm) for the three different projection angles $\phi=0$, $\pi/4$, and $\pi/2$,
respectively. The experimental correlation functions $C_{463,4}(|{\bf s}|)$
are compared in Fig. 5 with the theoretical ones. In Fig. 5 (a), even for
small $|{\bf s}|$, we find a significant departure of the experimental
correlation function from the theoretical prediction, which clearly suggests
that the wave function $\psi_{463,4}(r,\theta)$ is not chaotic. Also in Fig. 5
(b)-(c) the experimental correlation functions $C_{463,4}(|{\bf s}|)$ show for
larger $|{\bf s}|$ significant deviations from the theoretical ones. The
discrepancies between the correlation function $C_{463,4}({\bf x},|{\bf
s}|=z)$ for the $z$-component of the electric field distribution and the
theoretical prediction in Fig. 5(c) arise mainly due to the procedure of
averaging of the correlation function $C_{N,p}({\bf x},|{\bf s}|)$ in
$z$-direction, which was taken over the period of the cosine function.

In summary, we measured the wave functions of the chaotic 3D rough microwave
billiard with the translational symmetry up to the level number $N=489$. For
the first time the experimental correlation function $C_{N,p}({\bf x},{\bf
s})$ was estimated and compared with the theoretical prediction. For the
states with higher $N_{\bot}$ we find, especially for small values of the
parameter $|{\bf s}|$, good agreement with the theoretical predictions, which
show that the wave functions are chaotic. For the states with lower $N_{\bot}$
significant discrepancies between experimental and theoretical results are
observed.

Acknowledgments. This work was partially supported by the Ministry of Science
and Higher Education grants No. N202 099 31/0746 and 2 P03B 047 24.

## References

  * (1) S. Deus, P.M. Koch, and L. Sirko, Phys. Rev. E 52, 1146 (1995).
  * (2) H. Alt, C. Dembowski, H.-D. Gräf, R. Hofferbert, H. Rehfeld, A. Richter, R. Schuhmann, and T. Weinland, Phys. Rev. Lett. 79, 1026 (1997).
  * (3) U. Dörr, H.-J. Stöckmann, M. Barth, and U. Kuhl, Phys. Rev. Lett. 80, 1030 (1998).
  * (4) B. Eckhardt, U. Dörr, U. Kuhl, and H.-J. Stöckmann, Europhys. Lett. 46, 134 (1999).
  * (5) C. Dembowski, B. Dietz, H.-D. Gräf, A. Heine, T. Papenbrock, A. Richter, and C. Richter, Phys. Rev. Lett. 89, 064101 (2002).
  * (6) H. Primack and U. Smilansky, Phys. Rev. Lett. 74, 4831 (1995).
  * (7) T. Prosen, Phys. Lett. A 233, 323 (1997).
  * (8) L. R. Arnaut, Phys. Rev. E 73, 036604 (2006).
  * (9) Y.-H. Kim, U. Kuhl, H.-J. Stöckmann, and J. P. Bird, J. Phys.: Condens. Matter 17, L191 (2005).
  * (10) H.-J. Stöckmann, J. Stein, Phys. Rev. Lett. 64, 2215 (1990).
  * (11) R. Balian and C. Bloch, Ann. Phys. (N.Y.) 84, 559 (1974); Ann. Phys. (N.Y.) 64, 271(E) (1971).
  * (12) S.W. McDonald and A.N. Kaufman, Phys. Rev A 37, 3067 (1988).
  * (13) Y. Hlushchuk, A. Błȩdowski, N. Savytskyy, and L. Sirko, Physica Scripta 64, 192 (2001).
  * (14) K.M. Frahm and D.L. Shepelyansky, Phys. Rev. Lett. 78, 1440 (1997).
  * (15) K.M. Frahm and D.L. Shepelyansky, Phys. Rev. Lett. 79, 1833 (1997).
  * (16) N. Savytskyy and L. Sirko, Phys. Rev. E 65, 066202-1 (2002).
  * (17) L.C. Maier and J.C. Slater, J. Appl. Phys. 23, 68 (1952).
  * (18) N. Savytskyy, O. Hul, and L. Sirko, Phys. Rev. E 70, 056209 (2004).
  * (19) O. Hul, N. Savytskyy, O. Tymoshchuk, S. Bauch, and L. Sirko, Phys. Rev. E 72, 066212 (2005).
  * (20) O. Hul, N. Savytskyy, O. Tymoshchuk, S. Bauch, and L. Sirko, Acta Phys. Pol. A 109, 73 (2006).
  * (21) M.V. Berry, J. Phys. A 10, 2083 (1977).
  * (22) A. Bäcker and R. Schubert, J. Phys. A 35, 539 (2002).


remarkRemark hypothesisHypothesis A SINGLE PDE FOR THE MAXWELL’S EQUATIONSV.
G. Romanov and M. V. Klibanov

# Can a Single PDE Govern Well the Propagation of the Electric Wave Field in a
Heterogeneous Medium in 3D? ††thanks: Submitted to the editors DATE. The work
of V.G. Romanov was supported by Mathematical Center in Akademgorodok at
Novosibirsk State University (the agreement with Ministry of Science and High
Education of the Russian Federation number 075-15-2019-1613). The work of M.V.
Klibanov was supported by US Army Research Laboratory and US Army Research
Office grant W911NF-19-1-0044.

Vladimir G. Romanov Sobolev Institute of Mathematics, Novosibirsk 630090,
Russian Federation, and Mathematical Center in Akademgorodok, Novosibirsk
State University, Russian Federation<EMAIL_ADDRESS>Michael V.
Klibanov Corresponding author. Department of Mathematics and Statistics,
University of North Carolina at Charlotte, Charlotte, NC, 28223 USA ().
<EMAIL_ADDRESS>
###### Abstract

It is demonstrated in this paper that the propagation of the electric wave
field in a heterogeneous medium in 3D can sometimes be governed well by a
single PDE, which is derived from the Maxwell’s equations. The corresponding
component of the electric field dominates two other components. This justifies
some past results of the second author with coauthors about numerical
solutions of coefficient inverse problems with experimental electromagnetic
data. In addition, since it is simpler to work in applications with a single
PDE rather than with the complete Maxwell’s system, then the result of this
paper might be useful to researchers working on applied issues of the
propagation of electromagnetic waves in inhomogeneous media.

###### keywords:

Maxwell’s equations, geodesic lines, domination of one component, experimental
data for inverse problems

35Q61, 35R30

##  1 Introduction

In some previous works of the second author with coauthors coefficient inverse
problems were solved for frequency dependent microwave experimental
electromagnetic data using only the single Helmholtz equation, see, e.g. [6,
7, 8]. Reconstruction results were quite accurate ones. A similar observation
took place in [1, 15], although for time dependent experimental data. Thus, a
natural question to pose is:_Given that the propagation of the electromagnetic
wave field is governed by the Maxwell’s equations, why the use of only a
single equation has provided accurate reconstruction results?_ A positive
heuristic answer to this question can be found in the classical textbook of M.
Born and E. Wolf [3, pages 695,696] for the frequency domain case. In
addition, this question was positively addressed numerically in [2] for the
time domain case and in [10] for the frequency domain case. It was
demonstrated computationally in [2, 10] that if the incident electric wave
field has only a single non zero component, then this component dominates two
other components while propagating through the medium, and its propagation is
well governed by a wave-like PDE. That PDE is either the Helmholtz equation in
the frequency domain or the corresponding hyperbolic equation in the time
domain.

The goal of this paper is to investigate the above question rigorously. We
believe that the results of this paper might be useful not only for an
analytical explanation of the accuracy of imaging results of [6, 7, 8] but
also for applied mathematicians, physicists and engineers working on various
topics of electromagnetic waves propagation. Indeed, it is clear that it is
easier to work in applications with a single PDE rather than with the whole
Maxwell’s system.

In section 2 we work in time domain. These results are used then in section 3,
where we derive our desired conclusion for the case of the frequency domain.
In section 4, we link our main Theorem 2 of section 3 with the above cited
results of [6, 7, 8]. In section 5 (Appendix) we prove a certain energy
estimate.

##  2 Time Domain

Consider the Maxwell’s equations in a non magnetic medium

(1)
$\mathrm{curl}\,\mathbf{H}=\varepsilon(\mathbf{x})\mathbf{E}_{t},~{}\mathrm{curl}\,\mathbf{E}=-\mathbf{H}_{t},~{}\mathrm{div}\,\mathbf{H}=0,~{}\mathbf{x}\in\mathbb{R}^{3},t>0,$

where $\varepsilon(\mathbf{x})$ is the spatially distributed dielectric
constant. We work in this paper with dimensionless variables, since variables
were made dimensionless in the above cited works about the experimental data.
Thus, in vacuum $\varepsilon(\mathbf{x})=1,$ and we also assume that the
magnetic permeability $\mu\equiv 1$. Let
$\mathbf{\nu}\in\mathbb{S}^{2}=\left\\{\left|\mathbf{\nu}\right|=1\right\\}$
be a unit vector of the direction of propagation of the incident electric wave
field. If the space $\mathbb{R}^{3}$ is vacuum, then equations (1) admit the
following solution [13]:

(2)
$\mathbf{E}^{0}(\mathbf{x},t)=\mathbf{j}\delta(t+t_{0}-\mathbf{x}\cdot\mathbf{\nu}),~{}\mathbf{H}^{0}(\mathbf{x},t)=\left(\mathbf{\nu}\times\mathbf{j}\right)\delta(t+t_{0}-\mathbf{x}\cdot\mathbf{\nu}),$

where $\delta(t)$ is the Dirac delta function, $t_{0}$ is an arbitrary number,
$\,\mathbf{x}\cdot\mathbf{\nu}$ denotes the scalar product of these two
vectors and $\mathbf{j}\cdot\mathbf{\nu}=0$. The vector $\mathbf{j}$ defines
the polarization of this wave,

(3) $\mathbf{j}=\left(j_{1},j_{2},j_{3}\right),$

where $j_{1},j_{2},j_{3}$ are some constants. In the sequel we assume that
$\mathbf{j}\in\mathbb{S}^{2}$. The orthogonality of vectors $\mathbf{j}$ and
$\mathbf{\nu}$ is necessary to satisfy the equation
$\mathrm{div}\mathbf{H}=0$.

Below $\mathbf{j}$ and $\mathbf{\nu}$ are assumed to be arbitrary but fixed
vectors. Therefore we do not indicate dependence of the solution and some
functions on these parameters for brevity, unless this is really necessary.
Still, we use the parameter $\mathbf{\nu}$ to indicate some domains and a
plane wave for a more clear understanding. Below vectors $\mathbf{j}$,
$\mathbf{\nu}$, $\mathbf{E}$, $\mathbf{H}$, etc. are row vectors, see, e.g.
(3).

Let $R>0$ be an arbitrary number. Consider the ball $B$ with the center at
$\left\\{0\right\\}$ and the radius $R$. Let the sphere $S=\partial B.$ Then

$B=\\{\mathbf{x}\in\mathbb{R}^{3}:\,|\mathbf{x}|<R\\},~{}S=\\{\mathbf{x}\in\mathbb{R}^{3}:\,|\mathbf{x}|=R\\}.$

We assume that

(4) $1\leq\varepsilon(\mathbf{x})\leq\varepsilon_{0}\>\text{ in }B\text{,
}\varepsilon(\mathbf{x})=1\text{ in }\mathbb{R}^{3}\diagdown B,\text{ }$

where $\varepsilon_{0}\geq 1$ is a constant, and

(5) $t_{0}=\min_{\mathbf{x}\in S}(\mathbf{x}\cdot\mathbf{\nu})=-R.$

Let the incident plane wave propagates in the vacuum for $t<0$ and meets the
heterogeneous ball $B$ at a moment of time $t=0$. Then the propagation of the
electromagnetic wave field is governed by the following Cauchy problem:

(6)
$\mathop{\mathrm{c}url}\mathbf{H}=\varepsilon\left(\mathbf{x}\right)\mathbf{E}_{t},~{}\mathop{\mathrm{c}url}\mathbf{E}=-\mathbf{H}_{t},~{}\left(\mathbf{x},t\right)\in\mathbb{R}^{4},$
(7)
$\mathbf{E}|_{t<0}=\mathbf{E^{0}}(\mathbf{x},t),~{}\mathbf{H}|_{t<0}=\mathbf{H^{0}}(\mathbf{x},t).$

For the sake of convenience, we reduce now problem (6), (7) to the case when
only the vector function $\mathbf{E}(\mathbf{x},t)$ is unknown. Since
$\mathop{\mathrm{d}iv}\left(\mathop{\mathrm{c}url}\mathbf{U}\right)=0$ for any
appropriate vector function $\mathbf{U,}$ then applying the operator
$\mathop{\mathrm{d}iv}$ to both sides of the first equation (6), we obtain

(8)
$\mathop{\mathrm{d}iv}\left(\varepsilon\left(\mathbf{x}\right)\mathbf{E}_{t}\right)=0.$

Integrating (8) with respect to $t$ and using (7), we obtain

(9)
$\mathop{\mathrm{d}iv}\left(\varepsilon\left(\mathbf{x}\right)\mathbf{E}(\mathbf{x},t)\right)-\mathop{\mathrm{d}iv}\left(\varepsilon\left(\mathbf{x}\right)\mathbf{E^{0}}\left(\mathbf{x},\mathbf{\nu,}0\right)\right)=0.$

Note that $\mathrm{div}\,\mathbf{E}^{0}(\mathbf{x},t)=0$, since
$\mathbf{j}\cdot\mathbf{\nu}=0$. Indeed,

$\mathrm{div}\,\mathbf{E}^{0}(\mathbf{x},t)=-(\mathbf{j}\cdot\mathbf{\nu})\delta^{\prime}(t+t_{0}-\mathbf{x}\cdot\mathbf{\nu})=0.$

Hence,

(10)
$\mathrm{div}\,(\varepsilon(\mathbf{x})\mathbf{E}^{0}(\mathbf{x},0))=\nabla\varepsilon(\mathbf{x})\cdot\mathbf{E}^{0}(\mathbf{x},0)=(\nabla\varepsilon(\mathbf{x})\cdot\mathbf{j})\delta(t_{0}-\mathbf{x}\cdot\mathbf{\nu})=0.$

This is because by (5)
$\text{supp}\>\\{\delta(t_{0}-\mathbf{x}\cdot\mathbf{\nu})\\}$ is the tangent
plane to $S$, namely $(\mathbf{x}\cdot\mathbf{\nu})=-R$, along which
$\nabla\varepsilon(\mathbf{x})=0$. Hence
$\mathrm{div}\,(\varepsilon(\mathbf{x})\mathbf{E}(\mathbf{x},t))=0$ for any
$t$. Hence, problem (6), (7) is reduced to the following problem with
$n(\mathbf{x})=\sqrt{\varepsilon(\mathbf{x})}:$

(11)
$n^{2}(\mathbf{x})\mathbf{E}_{tt}-\Delta\mathbf{E}-\nabla(\mathbf{E}\cdot\nabla\ln
n^{2}(\mathbf{x}))=0,~{}(\mathbf{x},t)\in\mathbb{R}^{4},~{}$ (12)
$\mathbf{E}|_{t<0}=\mathbf{E}^{0}(\mathbf{x},t).$

Define two domains $D_{-}\left(\mathbf{\nu}\right)$ and
$D_{+}\left(\mathbf{\nu}\right)$ as

$D_{-}\left(\mathbf{\nu}\right)=\\{\mathbf{x}\in\mathbb{R}^{3}:\,\mathbf{x}\cdot\mathbf{\nu}+R<0\\},$
$D_{+}\left(\mathbf{\nu}\right)=\\{\mathbf{x}\in\mathbb{R}^{3}:\mathbf{x}\cdot\mathbf{\nu}+R\geq
0\\}.$

Note that the domain $D_{-}\left(\mathbf{\nu}\right)$ is situated outside of
$B$, while $B\subset D_{+}\left(\mathbf{\nu}\right)$.

To define geodesic lines, we partially follow our paper [9]. The function
$n(\mathbf{x})$ generates the Riemannian metric

$d\tau=n(\mathbf{x})\left|d\mathbf{x}\right|,\text{
}|d\mathbf{x}|=\sqrt{(dx_{1})^{2}+(dx_{2})^{2}+(dx_{3})^{2}}.$

For each vector $\mathbf{\nu}\in\mathbb{S}^{2}$ define the plane
$\Sigma(\mathbf{\nu})$ as

(13)
$\Sigma(\mathbf{\nu})=\\{\mathbf{\xi}\in\mathbb{R}^{3}:\mathbf{\>\xi\cdot\nu}=-R\\}.$

Observe that the plane $\Sigma(\mathbf{\nu})$ is tangent to $S$ at the point
$\mathbf{\xi}_{\tan}=-R\mathbf{\nu}$. Hence, $\Sigma(\mathbf{\nu})\cap
B=\varnothing$ and the vector $\mathbf{\nu}$ is a normal vector to the plane
$\Sigma(\mathbf{\nu}).$ Consider an arbitrary point
$\mathbf{y}\in\Sigma(\mathbf{\nu}).$ This point can be represented as

(14)
$\mathbf{y}=\mathbf{y}(a_{2},a_{3})=-R\mathbf{\nu}+a_{2}\mathbf{e}_{2}+a_{3}\mathbf{e}_{3},\quad(a_{2},a_{3})\in\mathbb{R}^{2},$

where unit vectors $\mathbf{\nu}$, $\mathbf{e}_{2}=\mathbf{j}$,
$\mathbf{e}_{3}=\mathbf{\nu}\times\mathbf{j}$ form an orthogonal triple. Note
that vectors $\mathbf{e}_{2}$, $\mathbf{e}_{3}$ are parallel to the plane
$\Sigma(\mathbf{\nu}).$

Let the function $\varphi(\mathbf{x},\mathbf{\nu})$ be the solution of the
Cauchy problem for the eikonal equation,

(15)
$|\nabla_{\mathbf{x}}\varphi(\mathbf{x},\mathbf{\nu})|^{2}=n^{2}(\mathbf{x}),\quad\varphi(\mathbf{x},\mathbf{\nu})|_{\mathbf{x}\in\Sigma(\mathbf{\nu})}=0$

satisfying the following conditions:

(16) $\varphi(\mathbf{x},\mathbf{\nu})\left\\{\begin{array}[]{c}<0\text{ if
}\mathbf{x}\in D_{-}\left(\mathbf{\nu}\right),\\\ >0\text{ if }\mathbf{x}\in
D_{+}\left(\mathbf{\nu}\right).\end{array}\right.$

The number $|\varphi(\mathbf{x},\mathbf{\nu})|$ is the Riemannian distance
between the point $\mathbf{x}$ and the plane $\Sigma(\mathbf{\nu})$. From the
Physics standpoint, $|\varphi(\mathbf{x},\mathbf{\nu})|$ is the travel time
between the point $\mathbf{x}$ and the plane $\Sigma(\mathbf{\nu})$. For
$\mathbf{\xi}\cdot\mathbf{\nu}<-R,$ i.e. in the domain
$D_{-}\left(\mathbf{\nu}\right)$, the function
$\varphi(\mathbf{x},\mathbf{\nu})$ has the form
$\varphi(\mathbf{x},\mathbf{\nu})=\mathbf{x}\cdot\mathbf{\nu}+R$. To find the
function $\varphi(\mathbf{x},\mathbf{\nu})$ in the domain
$D_{+}\left(\mathbf{\nu}\right),$ we need to solve problem (15), (16) in this
domain. It is known that to do this, we need to solve the following Cauchy
problem for a system of ordinary differential equations [14]:

(17)
$\frac{d\mathbf{x}}{ds}=\frac{p(\mathbf{x,\nu})}{n^{2}(\mathbf{x})},\quad\frac{d\mathbf{p}(\mathbf{x,\nu})}{ds}=\nabla\ln
n(\mathbf{x}),\quad\frac{d\varphi(\mathbf{x,\nu})}{ds}=1,s>0,$ (18)
$\mathbf{x}|_{s=0}=\mathbf{y},~{}\mathbf{p}|_{s=0}=\mathbf{\nu},~{}\varphi|_{s=0}=0,$

where $\mathbf{y}\in\Sigma(\mathbf{\nu})$ is an arbitrary point of the plane
$\Sigma(\mathbf{\nu}),$ see (14), $s$ is the Riemannian arc length and

(19) $\mathbf{p}(\mathbf{x,\nu})=\nabla\varphi(\mathbf{x,\nu}).$

Note that equations (17), (18) imply that $\varphi(\mathbf{\xi,\nu})=s.$ In
particular, this means that the second condition (15) is satisfied. Equations
(17), (18) define a geodesic line $\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$ in
$D_{+}\left(\mathbf{\nu}\right)$ which connects the point $\mathbf{x}\in
D_{+}\left(\mathbf{\nu}\right)$ with the point
$\mathbf{y}\in\Sigma(\mathbf{\nu})$ and orthogonal to $\Sigma(\mathbf{\nu})$
at $\mathbf{y}$.

Cauchy problem (17), (18) has the unique solution

(20)
$\mathbf{x}=\mathbf{f}(s,a_{2},a_{3}),~{}\mathbf{p}=\mathbf{g}(s,a_{2},a_{3}).$

These equations define the bundle of geodesic lines, which go out from
different points $\mathbf{y}\in\Sigma(\mathbf{\nu})$ in direction
$\mathbf{\nu}$. To find the geodesic line
$\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$, we need to invert the first
equation (20) and calculate $a_{2}$ and $a_{3}$ and then to find
$\mathbf{y}\in\Sigma(\mathbf{\nu})$, using formula (14). Set

$D_{+}(\mathbf{\nu},R)=\\{\mathbf{x}:\,-R\leq\mathbf{x}\cdot\mathbf{\nu}\leq
R\\},\text{ }T=\max_{\mathbf{x}\in
D_{+}(\mathbf{\nu},R)}\varphi(\mathbf{x,\nu}).$

Hence, $D_{+}(\mathbf{\nu},R)\subset D_{+}(\mathbf{\nu}).$ Note that both sets
$D_{+}(\mathbf{\nu})$ and $D_{+}(\mathbf{\nu},R)$ are closed ones, i.e.
$D_{+}(\mathbf{\nu})=\overline{D_{+}(\mathbf{\nu})}$ and
$D_{+}(\mathbf{\nu},R)=\overline{D_{+}(\mathbf{\nu},R)}.$ Here $T$ is a finite
number. Indeed, let $C(\mathbf{\nu},R)$ be the circular cylinder with the
circle of the radius $R$ and the axis
$\mathbf{x}=(-R+s^{\prime})\mathbf{\nu}$, $s^{\prime}\geq 0$, with generating
lines orthogonal to the plane $\Sigma(\mathbf{\nu}).$ Consider the
intersection $D_{0}(\mathbf{\nu},R)$ of $C(\mathbf{\nu},R)$ with
$D_{+}(\mathbf{\nu},R),$

(21) $D_{0}(\mathbf{\nu},R)=C(\mathbf{\nu},R)\cap D_{+}(\mathbf{\nu},R).$

Then the function $\varphi(\mathbf{x,\nu})=\mathbf{x}\cdot\mathbf{\nu}+R\leq
2R$ for $\mathbf{x}\in D_{+}(\mathbf{\nu},R)\diagdown D_{0}(\mathbf{\nu},R)$,
since $\mathbf{x}\cdot\mathbf{\nu}\in[-R,R]$ in $D_{+}(\mathbf{\nu},R)$. On
the other hand the domain $D_{0}(\mathbf{\nu},R)$ is finite. Therefore,

$T=\max\left(2R,\max_{\mathbf{x}\in
D_{0}(\mathbf{\nu},R)}\varphi(\mathbf{x,\nu})\right).$

Consider the Jacobian

(22)
$J(\mathbf{x})=\frac{\partial(x_{1},x_{2},x_{3})}{\partial(s,a_{2},a_{3})}.$

From relations (17) and (18) follows that $\partial\mathbf{x}/\partial
s=\mathbf{p}=\mathbf{\nu}$ and $\partial\mathbf{x}/\partial
a_{k}=\mathbf{e}_{k}$, $k=2,3$, at $s=0$. Then the Jacobian is the determinant
which rows are formed by components of three unite orthogonal vectors of the
positive orientation. Hence,

(23) $J(\mathbf{x})=1\text{ for }\mathbf{x}\in\Sigma(\mathbf{\nu}).$

Note that $J(\mathbf{x})=1$ for $\mathbf{x}\in D_{0}(\mathbf{\nu},R)$ as well
since $\mathbf{x}=\mathbf{y}+s\mathbf{\nu}$ in $D_{0}(\mathbf{\nu},R)$.

Denote by $c(\mathbf{x})=1/\sqrt{\varepsilon(\mathbf{x})}$ the speed of
propagation of electromagnetic waves. By (4)

$c_{0}\leq c(\mathbf{x})\leq 1,~{}\mathbf{x}\in\mathbb{R}^{3},$

where $c_{0}=1/\sqrt{\varepsilon_{0}}$.

Below we use the following assumptions:

Assumptions:

1. _The function_ $\varepsilon(\mathbf{x})\in C^{\infty}(\mathbb{R}^{3}),$ _satisfies conditions ( 4)._

2._There exists a positive constant_ $J_{0}$_such that_ $J(\mathbf{x})\geq
J_{0}$ for $D_{+}(\mathbf{\nu},R)$.

3\. _Any point_ $\mathbf{x}\in D_{+}(\mathbf{\nu},R)$_can be connected with
the plane_ $\Sigma\left(\mathbf{\nu}\right)$_by a single geodesic line_
$\Gamma(\mathbf{x},\Sigma(\nu))$_such that_
$\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$_is orthogonal to_
$\Sigma\left(\mathbf{\nu}\right)$_at a point_
$\mathbf{y}\in\Sigma\left(\mathbf{\nu}\right)$.

4\. _Any two points_ $\mathbf{x}$ _and_ $\mathbf{y}$ _in_ $\mathbb{R}^{3}$
_can be connected by a single geodesic line._

Under these Assumptions, the equality $\mathbf{x}=\mathbf{f}(s,a_{2},a_{3})$
is invertible in $D_{+}(\mathbf{\nu},R)$ and defines
$s=s(\mathbf{x},\mathbf{\nu})=\varphi(\mathbf{x,\nu})$ and parameters
$a_{k}=a_{k}(\mathbf{x,\nu})$, $k=2,3$, i.e. the point
$\mathbf{y}\in\Sigma\left(\mathbf{\nu}\right)$, see (14). By (17) and (18) if
$\mathbf{x}=\mathbf{y},$ then $d\mathbf{x}/ds=\mathbf{\nu}$. The latter vector
is directed along the geodesic line $\Gamma(\mathbf{x},\Sigma(\mathbf{\nu})).$
Hence, $\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$ is orthogonal to
$\Sigma\left(\mathbf{\nu}\right)$ at the point $\mathbf{y}$.

Define domains $G_{T}\left(\mathbf{\nu}\right)$ and
$G\left(T,\mathbf{\nu}\right)$ in $\mathbb{R}^{4}$ as

(24) $G_{T}\left(\mathbf{\nu}\right)=\\{(\mathbf{x},t):\,0\leq
t<\min(|\varphi(\mathbf{x,\nu})|,T)\\},$ (25)
$G\left(T,\mathbf{\nu}\right)=\\{(\mathbf{x},t)\,,|\varphi(\mathbf{x,\nu})|\leq
t\leq T\\}.$

To differentiate between notations of the Heaviside function $H(t)$ and the
magnetic wave field $\mathbf{H}(\mathbf{x},t),$ it is convenient to denote
$\theta_{0}(t):=$ $H(t),$

$\theta_{0}(t)=\left\\{\begin{array}[]{c}1,t\geq 0,\\\
0,t<0.\end{array}\right.$

Theorem 1. _Assume that the Assumption holds. Then for every vector_
$\mathbf{\nu\in}\mathbb{S}^{2}$ _the solution of problem ( 11), (12) can be
represented in _$\mathbb{R}_{T}^{4}=\\{(\mathbf{x},t)|\,0\leq t\leq T\\}$ _in
the form_

(26)
$\mathbf{E}(\mathbf{x},t)=\mathbf{\alpha}^{-1}(\mathbf{x})\delta(t-\varphi(\mathbf{x}))+\mathbf{\widehat{E}}(\mathbf{x},t)\theta_{0}(t-|\varphi(\mathbf{x})|),$

_where_ $\mathbf{\alpha}^{-1}(\mathbf{x})\in
C^{\infty}(D_{+}(\mathbf{\nu},R))$, $\mathbf{\alpha}^{-1}(\mathbf{x})=0$ _for_
$\mathbf{x}\in D_{-}$ _and_ $\mathbf{\widehat{E}}(\mathbf{x},t)\in
C^{2}\left(G(T,\mathbf{\nu})\right)\,$ _and_
$\mathbf{\widehat{E}}(\mathbf{x},t)=0$ _for_
$t=|\varphi(\mathbf{x},\mathbf{\nu})|$.

Proof. Introduce functions $\theta_{k}(t)$ as

(27) $\theta_{-3}(t)=\delta^{\prime\prime}(t),\text{
}\theta_{-2}(t)=\delta^{\prime}(t),\text{ }\theta_{-1}(t)=\delta(t),\text{
}\theta_{k}(t)=\frac{t^{k}}{k!}\theta_{0}(t),\quad k=1,2,\ldots.$

Observe that $\theta_{k}^{\prime}(t)=\theta_{k-1}(t)$ for all $k\geq-2$. We
seek the solution of problem (11), (12) in the form

(28)
$\mathbf{E}(\mathbf{x},t)=\sum\limits_{k=-1}^{r}\mathbf{\alpha}^{k}(\mathbf{x})\theta_{k}(t-\varphi(\mathbf{x}))+\mathbf{E}^{r}(\mathbf{x},t),$

where the natural number $r$ will be chosen later. Substituting representation
(28) in (11), using the eikonal equation (15) and equating coefficients at
$\theta_{k}(t)$ for $k=-2,-1,0,1,\ldots,r-1$, we obtain the following
reqursive formulas for finding coefficients $\mathbf{\alpha}^{k}(\mathbf{x})$:

(29) $\displaystyle
2(\nabla\varphi(\mathbf{x})\cdot\nabla)\mathbf{\alpha}^{k}(\mathbf{x})+\mathbf{\alpha}^{k}(\mathbf{x})\Delta\varphi(\mathbf{x})+(\mathbf{\alpha}^{k}(\mathbf{x})\cdot\nabla\ln
n^{2}(\mathbf{x}))\nabla\varphi(\mathbf{x})$
$\displaystyle=\Delta\mathbf{\alpha}^{k-1}(\mathbf{x})+\nabla(\mathbf{\alpha}^{k-1}(\mathbf{x})\cdot\nabla\ln
n^{2}(\mathbf{x})),~{}k=-1,0,1,...,r.$

Here we need to formally set

(30) $\mathbf{\alpha}^{-2}(\mathbf{x})=0.$

Since by (2) and (12)
$\mathbf{E}=\mathbf{j}\,\delta(t-(\mathbf{x}\cdot\mathbf{\nu}+R))$ for $t<0$,
then we obtain

(31) $\begin{array}[]{lll}\mathbf{\alpha}^{-1}(\mathbf{x})=0,\>\text{in
}\>D_{-};&\mathbf{\alpha}^{-1}|_{\Sigma({\mathbf{\nu)}}}=\mathbf{j},&\\\
\mathbf{\alpha}^{k}(\mathbf{x})=0,\>\text{in
}\>D_{-};&\mathbf{\alpha}^{k}|_{\Sigma({\mathbf{\nu)}}}=0,&k=0,1,\ldots,r.\end{array}$

Moreover, using (11), we obtain the following Cauchy problem for the residual
$\mathbf{E}^{r}(\mathbf{x},t)$ of expansion (28)

(32)
$n^{2}(\mathbf{x})\partial_{t}^{2}\mathbf{E}^{r}-\Delta\mathbf{E}^{r}-\nabla(\mathbf{E}^{r}\cdot\nabla\ln
n^{2}(\mathbf{x}))=\mathbf{F}^{r}(\mathbf{x},t),~{}\mathbf{E}^{r}|_{{t<0}}=0,$

where

(33)
$\mathbf{F}^{r}(\mathbf{x},t)=\left(\Delta\mathbf{\alpha}^{r}(\mathbf{x})+\nabla\left(\mathbf{\alpha}^{r}(\mathbf{x})\cdot\nabla\ln
n^{2}(\mathbf{x})\right)\right)\theta_{r}(t-\varphi(\mathbf{x})).$

We now construct solutions of equation (29) with the Cauchy data (31). We have
along the geodesic line $\Gamma(\mathbf{x},\Sigma(\mathbf{\nu})):$

(34)
$2(\nabla\varphi(\mathbf{\xi})\cdot\nabla)\mathbf{\alpha}^{k}(\mathbf{\xi})=2(\mathbf{p}(\mathbf{\xi})\cdot\nabla)\mathbf{\alpha}^{k}(\mathbf{\xi})=2n^{2}(\mathbf{\xi})\left(\frac{d\mathbf{\xi}}{ds}\cdot\nabla\right)\mathbf{\alpha}^{k}(\mathbf{\xi})$
$=2n^{2}(\mathbf{\xi})\frac{d\mathbf{\alpha}^{k}(\mathbf{\xi})}{ds},$

where $\mathbf{\xi}\in\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$ is an arbitrary
point. Moreover, it is stated in the paper [12] that the following formula
valid along $\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$ (see Lemma 1, Equation
(4.2)):

(35) $\frac{d\ln
J(\mathbf{\xi})}{ds}=\mathrm{div}\,\left(n^{-2}(\mathbf{\xi})\nabla\varphi(\mathbf{\xi})\right),$

where $J(\mathbf{\xi})$ is the Jacobian defined in (22), provided that
$\mathbf{x}=\left(x_{1},x_{2},x_{3}\right)$ is replaced with
$\mathbf{\xi}=\left(\xi_{1},\xi_{2},\xi_{3}\right)$. Calculating the right-
hand-side of (35) and using (35), we find

$\mathrm{div}\,\left(n^{-2}(\mathbf{\xi})\nabla\varphi(\mathbf{\xi})\right)=n^{-2}(\mathbf{\xi})\Delta\varphi(\mathbf{\xi})+\left(\nabla
n^{-2}(\mathbf{\xi})\cdot\mathbf{p}(\mathbf{\xi})\right).$

Hence, (35) implies

(36) $n^{-2}(\mathbf{\xi})\Delta\varphi(\mathbf{\xi})=\frac{d\ln
J(\mathbf{\xi})}{ds}-\left(\nabla
n^{-2}(\mathbf{\xi})\cdot\mathbf{p}(\mathbf{\xi})\right)=\frac{d\ln(J(\mathbf{\xi})n^{2}(\mathbf{\xi}))}{ds}.$

We have used here that, similarly (34),

$-\left(\nabla n^{-2}(\mathbf{\xi})\cdot
p(\mathbf{\xi})\right)=-n^{2}(\mathbf{\xi})\left(\nabla
n^{-2}(\mathbf{\xi})\cdot\frac{d\mathbf{\xi}}{ds}\right)=-n^{2}(\mathbf{\xi})\frac{d}{ds}n^{-2}(\mathbf{\xi})=\frac{d}{ds}\ln
n^{2}(\mathbf{\xi}).$

Using formulae (34) and (36) and replacing in (29) $\mathbf{x}$ with
$\mathbf{\xi}$, we obtain

$\displaystyle
2n^{2}({\mathbf{\xi}})\left[\frac{d\mathbf{\alpha}^{k}(\mathbf{\xi})}{ds}+\mathbf{\alpha}^{k}(\mathbf{\xi})\frac{d\ln(\sqrt{J(\mathbf{\xi})}n(\mathbf{\xi}))}{ds}\right]+(\mathbf{\alpha}^{k}(\mathbf{\xi})\cdot\nabla\ln
n^{2}(\mathbf{\xi}))\nabla\varphi(\mathbf{\xi})$ (37)
$\displaystyle=\Delta\mathbf{\alpha}^{k-1}(\mathbf{\xi})+\nabla(\mathbf{\alpha}^{k-1}(\mathbf{\xi})\cdot\nabla\ln
n^{2}(\mathbf{\xi})),~{}k=-1,0,1,...,r.$

Multiplying equation (37) by
$\sqrt{J(\mathbf{\xi})}/\left(2n(\mathbf{\xi})\right)$, we transform equation
(29) along $\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$ to the recursive form

$\displaystyle\frac{d}{ds}\left(\mathbf{\alpha}^{k}(\mathbf{\xi})n(\mathbf{\xi})\sqrt{J(\mathbf{\xi})}\right)-\frac{n(\mathbf{\xi})\sqrt{J(\mathbf{\xi})}}{2}\left(\mathbf{\alpha}^{k}(\mathbf{\xi})\cdot\nabla
n^{-2}(\mathbf{\xi})\right)\mathbf{p}(\mathbf{\xi})$ (38)
$\displaystyle=\mathbf{Q}^{k}(\mathbf{\xi}),k=-1,0,1,\ldots,r,$

were

(39)
$\mathbf{Q}^{k}(\mathbf{\xi})=\frac{\sqrt{J(\mathbf{\xi})}}{2n(\mathbf{\xi})}\big{[}\Delta\mathbf{\alpha}^{k-1}(\mathbf{\xi})+\nabla\left(\mathbf{\alpha}^{k-1}(\mathbf{\xi})\cdot\nabla\ln
n^{2}(\mathbf{\xi})\right)\big{]}.$

When we solve equations (38) going from $k$ to $k+1$, the function
$\mathbf{Q}^{k}(\mathbf{\xi})$ is always known from the previous step. It
follows from (30) and (39) that $\mathbf{Q}^{-1}(\mathbf{\xi})=0$. Functions
$\mathbf{\alpha}^{k}(\mathbf{\xi})$ satisfy on $\Sigma(\mathbf{\nu})$
conditions (31). We also recall that by (4) and (21)
$n\left(\mathbf{\xi}\right)|_{\Sigma(\mathbf{\nu})}=1$,
$J|_{\Sigma(\mathbf{\nu})}=1.$

Hence, integrating (38) with respect to $s$, we obtain a recursive integral
equation along the geodesic line $\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$

(40) $\displaystyle\mathbf{\alpha}^{k}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\frac{1}{n(\mathbf{x})\sqrt{J(\mathbf{x})}}\bigg{(}\mathbf{A}^{k}+\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}\left[\mathbf{Q}^{k}(\mathbf{\xi})\right.$
$\displaystyle\left.+\frac{n(\mathbf{\xi})\sqrt{J(\mathbf{\xi})}}{2}\left(\mathbf{\alpha}^{k}(\mathbf{\xi})\cdot\nabla
n^{-2}(\mathbf{\xi})\right)p(\mathbf{\xi})\right]\,ds\bigg{)},~{}$
$\displaystyle k$ $\displaystyle=$ $\displaystyle-1,0,1,\ldots,r,$

where $\mathbf{\xi}=\mathbf{f}(s,a_{2},a_{3})$, the vector function
$\mathbf{f}(s,a_{2},a_{3})$ is defined in (20) and

$\mathbf{A}^{-1}=\mathbf{j},~{}\mathbf{A}^{k}=0,~{}k=0,1,\ldots,r.$

Recall that $s$ is the Riemannian arc length of
$\Gamma(\mathbf{\xi},\Sigma(\mathbf{\nu}))$. Equation (40) is a Volterra-type
integral equation of the second kind along the curve
$\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$. Therefore, this equation can be
solved by the method of successive approximations which is rapidly converging.
We can solve equation (40) for different $k=-1,...,r$ step-by-step, starting
from $k=-1$.

As it will be clear in the sequel, the most important role plays the function
$\mathbf{\alpha}^{-1}(\mathbf{x})$. We now represent this function through the
resolvent $R(\mathbf{x},\mathbf{\xi})$ of equation (40) for $k=-1$. First, we
introduce the vector function $\mathbf{\beta(x)}$

$\mathbf{\beta(x)}={n(\mathbf{x})\sqrt{J(\mathbf{x})}}\mathbf{\alpha}^{-1}(\mathbf{x}).$

Then equation for this function has the form

$\mathbf{\beta(x)}=\mathbf{j}+\frac{1}{2}\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}\left(\mathbf{\beta(\xi)}\cdot\nabla
n^{-2}(\mathbf{\xi})\right)\mathbf{p}(\mathbf{\xi})\,ds.$

The more convenient form of this equation is:

$\mathbf{\beta}(\mathbf{x})=\mathbf{j}+\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}\mathbf{\beta}(\mathbf{\xi})K_{0}(\mathbf{\xi})\,ds,$

where $K_{0}({\mathbf{\xi}})$ is the $3\times 3$ matrix

$K_{0}({\mathbf{\xi}})=\frac{1}{2}(\nabla
n^{-2}(\mathbf{\xi}))^{\ast}\mathbf{p}(\mathbf{\xi})$

and $(\nabla n^{-2}(\mathbf{\xi}))^{\ast}$ is the transposed vector $(\nabla
n^{-2}(\mathbf{\xi}))$, i.e. column vector, while $\mathbf{p}(\mathbf{\xi})$
is row vector.

Represent $\mathbf{\beta(x)}$ as

$\mathbf{\beta(x)}=\sum\limits_{n=0}^{\infty}\mathbf{\beta}^{n}(\mathbf{x}),$

where

$\mathbf{\beta}^{0}(\mathbf{x})=\mathbf{j},~{}\mathbf{\beta}^{n}(\mathbf{x})=\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}\mathbf{\beta}^{n-1}(\mathbf{\xi})K_{0}(\mathbf{\xi})\,ds,~{}n=1,2,\ldots.$

Then

$\mathbf{\beta}^{1}(\mathbf{x})=\mathbf{j}\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}K_{0}({\mathbf{\xi}})\,ds,$

Next,

$\mathbf{\beta}^{2}(\mathbf{x})=\mathbf{j}\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}\int\limits_{\Gamma(\mathbf{\xi},\Sigma(\mathbf{\nu}))}K_{0}({\mathbf{\xi}}^{\prime})\,ds^{\prime}K_{0}({\mathbf{\xi}})\,ds.$

Changing here the repeated integration by place and then replacing
$s^{\prime}$ with $s$ and vice versa, we obtain

$\mathbf{\beta}^{2}(\mathbf{x})=\mathbf{j}\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}K_{1}({\mathbf{x}},{\mathbf{\xi}})\,ds,$

where

$K_{1}(\mathbf{x},\mathbf{\xi})=\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))\setminus\Gamma(\mathbf{\xi},\Sigma(\mathbf{\nu}))}K_{0}(\mathbf{\xi})K_{0}(\mathbf{\xi}^{\prime})ds^{\prime}$

Similarly,

$\displaystyle\mathbf{\beta}^{n}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\mathbf{j}\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}\int\limits_{\Gamma(\mathbf{\xi},\Sigma(\mathbf{\nu}))}K_{n-2}(\mathbf{\
\xi},{\mathbf{\xi}}^{\prime})\,ds^{\prime}K_{0}({\mathbf{\xi}})\,ds$
$\displaystyle=$
$\displaystyle\mathbf{j}\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}K_{n-1}({\mathbf{x}},{\mathbf{\xi}})\,ds,~{}n\geq
2,$

where

$K_{n-1}(\mathbf{x},\mathbf{\xi})=\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))\setminus\Gamma(\mathbf{\xi},\Sigma(\mathbf{\nu}))}K_{n-2}(\mathbf{\xi}^{\prime},\mathbf{\xi})K_{0}(\mathbf{\xi}^{\prime})ds^{\prime}$

and $K_{0}(\mathbf{x},\mathbf{\xi})=K_{0}(\mathbf{\xi})$. Thus, we obtain

$\mathbf{\beta(x)}=\mathbf{j}\left(I+\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}R(\mathbf{x},\mathbf{\xi})ds\right).$

Here $I$ is the identity matrix and $R(\mathbf{x},\mathbf{\xi})$ is defined as

(41)
$R(\mathbf{x},\mathbf{\xi})=\sum_{n=0}^{\infty}K_{n}(\mathbf{x},\mathbf{\xi}).$

Finally we obtain for $\mathbf{\alpha}^{-1}(\mathbf{x})$ the following formula

$\mathbf{\alpha}^{-1}(\mathbf{x})=\frac{{\mathbf{j}}}{n(\mathbf{x})\sqrt{J(\mathbf{x})}}\left(I+\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}R(\mathbf{x},\mathbf{\xi})ds\right).$

Applying the above technique, similar formulae can be easily derived for
$\mathbf{\alpha}^{k}(\mathbf{x})$, $k=0,1,\ldots,r$.

Denote

$P^{k}({\mathbf{x}})=\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}\mathbf{Q}^{k}(\mathbf{\xi})\,ds.$

Then

(42)
$\mathbf{\alpha}^{k}(\mathbf{x})=\frac{1}{n(\mathbf{x})\sqrt{J(\mathbf{x})}}\left(P^{k}({\mathbf{x}})+\int\limits_{\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))}P^{k}({\mathbf{\xi}})R(\mathbf{x},\mathbf{\xi})ds\right),~{}k=0,1,\ldots,r.$

We estimate $R(\mathbf{x},\mathbf{\xi})$ later in the proof of Theorem 2. The
uniform convergence of series (41) in $D_{+}(\mathbf{x},R)$ follows from that
estimate.

Since the function $n(\mathbf{x})\in C^{\infty}(\mathbb{R}^{3}),$ then
functions $\mathbf{f}(s,a_{2},a_{3})$ and $J(\mathbf{x})$ belong to
$C^{\infty}({D_{+}(\mathbf{\nu},R)})$. Hence, all functions
$\mathbf{\alpha}^{k}(\mathbf{x})\in C^{\infty}({D_{+}(\mathbf{\nu},R)})$.
Moreover, $\mathbf{\alpha}^{-1}(\mathbf{x})=\mathbf{j}$ and
$\mathbf{\alpha}^{k}(\mathbf{x})=0$ for $k=0,1,\ldots,r$, if $\mathbf{x}$ lies
outside $D_{0}(\mathbf{\nu},R)$ since
$\text{supp}\nabla\varepsilon(\mathbf{x})\subset B$. Therefore, functions
$\mathbf{\alpha}^{k}(\mathbf{x})$ for $k=0,1,\ldots,r$ are compactly supported
in $D_{+}(\mathbf{\nu},R)$. Hence, $\mathbf{F}^{r}(\mathbf{x},t)=0$ for
$\\{(\mathbf{x},t)|\,\mathbf{x}\in D_{0}(\mathbf{\nu},R),t\geq 0\\}$ (see
notation (21)). Moreover, $\mathbf{F}^{r}(\mathbf{x},t)=0$ for
$(\mathbf{x},t)\in G_{T}(\mathbf{\nu})$. These two facts imply that function
$\mathbf{F}^{r}(\mathbf{x},t)$ is compactly supported in domain
$\mathbb{R}_{T}^{4}\supset G(T,\mathbf{\nu})$. Hence, it follows from (32)
that the vector function $\mathbf{E}^{r}(\mathbf{x},t)$ vanishes in
$G_{T}(\mathbf{\nu})$ and it is compactly supported in $\mathbb{R}_{T}^{4}$
because the speed of electromagnetic waves is finite.

We now apply the method of energy estimates to the problem (32) in the domain
$\mathbb{R}_{T}^{4}\supset G(T,\mathbf{\nu})$ to estimate function
$\mathbf{E}^{r}(\mathbf{x},t)$. . This method is a powerful tool for
investigations of various problems of mathematical physics. Applications of
this methods for boundary value problems are given in [4], [11] and many other
books related to partial differential equations. So, main ideas of this method
are well known. At the same time, we can not give a reference of the exact
result that we need for our goal. Therefore, for the completeness of the proof
and for the reader’s convenience we formulate below a lemma related to the
estimate of solution to problem (32). Let
$Y(t,T)=\mathbb{R}_{T}^{4}\cap\\{t=const\\}$.

Lemma. _Let_ $\varepsilon(\mathbf{x})$ _satisfy the conditions_

(43)
$\|\varepsilon\|_{C^{r}(\mathbb{R}^{3})}\leq\mu,~{}\|\ln\varepsilon\|_{C^{r+1}(\mathbb{R}^{3})}\leq\mu,$

_with a positive constant_ $\mu$ _and function_ $\mathbf{F}^{r}(\mathbf{x},t)$
_belongs to_ $H^{r}(\mathbb{R}_{T}^{4})$ _and satisfies the inequality_

(44) $\|\mathbf{F}^{r}\|_{H^{r}(\mathbb{R}_{T}^{4})}\leq M.$

_Then the solution of problem_ (32) $\mathbf{E}^{r}\in H^{r+1}(Y(t,T))$ _for
all_ $t\in\left(0,T\right)$ _and_ _there exists a positive constant_
$C_{1}=C_{1}(\mu,T)$ _such that the following estimates hold:_

(45) $\|\mathbf{E}^{r}\|_{H^{r+1}(Y(t,T))}\leq
C_{1}M,~{}\|\partial_{t}\mathbf{E}^{r}\|_{H^{r}(Y(t,T))}\leq C_{1}M.$

The proof of this Lemma is given in the Appendix.

In our case $\mathbf{F}^{r}\in{H^{r}(\mathbb{R}_{T}^{4})}$ and $\varepsilon\in
C^{\infty}(\mathbb{R}^{3})$. Hence, conditions (43) and (44) valid with some
positive $\mu$ and $M$. Applying Lemma, we obtain that the solution of problem
(32) is such that $\mathbf{E}^{r}\in H^{r+1}(Y(t,T))$ and
$\partial_{t}\mathbf{E}^{r}\in H^{r}(Y(t,T))$ for all $t\in[0,T]$. Choosing
$r=4$ and applying Lemma, we obtain $\mathbf{E}^{4}\in H^{5}(Y(t,T))$ for all
$t\in[0,T]$ and, hence, $\mathbf{E}^{4}\in H^{5}(\mathbb{R}_{T}^{4})$.
Therefore the embedding theorem implies that $\mathbf{E}^{4}\in
C^{2}(\mathbb{R}_{T}^{4})$. Hence, the vector function $\mathbf{E}^{4}\in
C^{2}\left(\overline{G(T,\mathbf{\nu})}\right)$ and is continuous together
with space derivatives up to the second order across the characteristic wedge
$t=|\varphi(\mathbf{x})|$. In particular, $\mathbf{E}^{4}(\mathbf{x},t)=0$ for
$t=|\varphi(\mathbf{x})|$.

Setting

(46)
$\widehat{\mathbf{E}}(\mathbf{x},t)=\sum\limits_{k=0}^{4}\mathbf{\alpha}^{k}(\mathbf{x})\frac{(t-\varphi(\mathbf{x}))^{k}}{k!}+\mathbf{E}^{4}(\mathbf{x},t),$

we obtain (26) as well as the required smoothness $\widehat{\mathbf{E}}\in
C^{2}\left(\overline{G(\mathbf{\nu},T)}\right)$. $\square$

Remark 1. The equality (46) implies the following formula, which we use below:

(47)
$\lim_{t\rightarrow\varphi(\mathbf{x})^{+}}\widehat{\mathbf{E}}(\mathbf{x},t)=\mathbf{\alpha}^{0}(\mathbf{x}).$

##  3 Frequency domain

Consider the Fourier transform $\widetilde{\mathbf{E}}(\mathbf{x},k)$ of the
function $\mathbf{E}(\mathbf{x},t)$,

(48) $\displaystyle\widetilde{\mathbf{E}}(\mathbf{x},k)$ $\displaystyle=$
$\displaystyle\int\limits_{-\infty}^{\infty}\mathbf{E}(\mathbf{x},t)\exp\left(-ikt\right)dt$
$\displaystyle=$
$\displaystyle\widetilde{\mathbf{E}}^{0}(\mathbf{x},k)+\int\limits_{0}^{\infty}\mathbf{E}(\mathbf{x},t)\exp\left(-ikt\right)dt,$

where $k=2\pi/\lambda$ is the wave number, $\lambda$ is the dimensionless
wavelength and

$\widetilde{\mathbf{E}}^{0}(\mathbf{x},k)=\mathbf{j}\exp(i(\mathbf{x}\cdot\mathbf{\nu}+R))\theta_{0}(-\mathbf{x}\cdot\mathbf{\nu}-R)$

is the Fourier transform of $\mathbf{E}^{0}(\mathbf{x},t)$. The existence of
the integral in (48) follows from results of Vainberg [17] which claim that
the vector function $\mathbf{E}(\mathbf{x},t)$ decays exponentially together
with its appropriate derivatives as $t\rightarrow\infty$ while $\mathbf{x}$
runs over any bounded domain $\Omega\subset\mathbb{R}^{3}$. Next, theorem 3.3
of [16] and theorem 6 of Chapter 9 of [17] guarantee that
$\widetilde{\mathbf{E}}(\mathbf{x},k)$ is the solution to the equation

(49)
$(\Delta+k^{2}n^{2}(\mathbf{x}))\widetilde{\mathbf{E}}+\nabla(\widetilde{\mathbf{E}}\cdot\nabla\ln
n^{2}(\mathbf{x}))=0,~{}\mathbf{x}\in\mathbb{R}^{3},$

where the scattering field

$\widetilde{\mathbf{E}}^{sc}(\mathbf{x},k)=\widetilde{\mathbf{E}}(\mathbf{x},k)-\widetilde{\mathbf{E}}^{0}(\mathbf{x},k)$

satisfies the radiation condition as $|\mathbf{x}|\rightarrow\infty$.

We now consider the vector function $\widetilde{\mathbf{E}}(\mathbf{x},k)$ in
(48) for $\mathbf{x}\in D_{+}(\mathbf{\nu},R)$. Using representation (26), we
obtain

(50)
$\widetilde{\mathbf{E}}(\mathbf{x},k)=\mathbf{\alpha}^{-1}(\mathbf{x})\exp(-ik\varphi(\mathbf{x}))+\int\limits_{\varphi(\mathbf{x)}}^{\infty}\widehat{\mathbf{E}}\left(\mathbf{x},t\right)\exp\left(-ikt\right)\,dt.$

Integrating by parts in (50) and using formula, we obtain

$\displaystyle\widetilde{\mathbf{E}}(\mathbf{x},k)$ $\displaystyle=$
$\displaystyle\mathbf{\alpha}^{-1}(\mathbf{x})\exp(-ik\varphi(\mathbf{x}))+\frac{\exp(-ik\varphi(\mathbf{x}))}{ik}\mathbf{\alpha}^{0}(\mathbf{x})$
$\displaystyle+\frac{1}{ik}\int\limits_{\varphi(\mathbf{x})}^{\infty}\widehat{\mathbf{E}}_{t}\left(\mathbf{x},t\right)\exp\left(-ikt\right)dt,\quad$

Thus,

(51)
$\widetilde{\mathbf{E}}(\mathbf{x},k)=\mathbf{\alpha}^{-1}(\mathbf{x})\exp(-ik\varphi(\mathbf{x}))+O\left(\frac{1}{k}\right),\quad
k\rightarrow\infty,~{}\forall\,\mathbf{x}\in D_{+}(\mathbf{\nu},R).$

Consider now the equation

(52)
$(\Delta+k^{2}n^{2}(\mathbf{x})\overline{\mathbf{E}}=0,~{}\mathbf{x}\in\mathbb{R}^{3},$

with the incident plane wave
$\overline{\mathbf{E}}^{0}(\mathbf{x},k)=\mathbf{j}\exp(ik(\mathbf{x}\cdot\mathbf{\nu}+R))\theta_{0}(-\mathbf{x}\cdot\mathbf{\nu}-R)$
and the radiation condition for
$\left(\overline{\mathbf{E}}-\overline{\mathbf{E}}^{0}\right)$. Note that
$\overline{\mathbf{E}}\cdot\mathbf{\nu}=0$ and
$\overline{\mathbf{E}}\cdot(\mathbf{\nu}\times\mathbf{j})=0$ since
$\overline{\mathbf{E}}^{0}\cdot\mathbf{\nu}=0$ and
$\overline{\mathbf{E}}^{0}\cdot(\mathbf{\nu}\times\mathbf{j})=0$ because
$\overline{\mathbf{E}}^{0}$ is parallel to $\mathbf{j}$ and also
$\mathbf{j}\cdot\mathbf{\nu}=0$.

Consider the function
$\overline{u}(\mathbf{x},t)=\overline{\mathbf{E}}\cdot\mathbf{j.}$ Then

(53)
$(\Delta+k^{2}n^{2}(\mathbf{x})\overline{u}(\mathbf{x},k)=0,~{}\mathbf{x}\in\mathbb{R}^{3},$

with the incident plane wave
$\overline{u}^{0}(\mathbf{x},k)=\exp(ik(\mathbf{x}\cdot\mathbf{\nu}+R))\theta_{0}(-\mathbf{x}\cdot\mathbf{\nu}-R)$
and the radiation condition for $\left(\overline{u}-\overline{u}^{0}\right)$.

We impose conditions below, which guarantee that the electric wave field
$\overline{\mathbf{E}}_{\mathbf{j}}=\mathbf{j}\overline{u}(\mathbf{x},k)$ is
close to $\widetilde{\mathbf{E}}$ at the high values of the wave number $k$,
which is equivalent to small wavelengths $\lambda$. First, suppose that these
two electric wave fields are indeed close to each other in the norm of the
space $C\left(\overline{{D_{+}(\mathbf{\nu},R)}}\right)$. This means that

(54)
$\widetilde{\mathbf{E}}(\mathbf{x},k)=\overline{\mathbf{E}}(\mathbf{x},k)+\mathbf{V}(\mathbf{x},k),\text{
}\mathbf{x}\in{D_{+}(\mathbf{\nu},R)},$ (55)
$\left\|\mathbf{V}(\mathbf{x},k)\right\|_{C\left({D_{+}(\mathbf{\nu},R)}\right)}\leq\sigma,$

where $\sigma>0$ is a small number. Consider the component
$\widetilde{u}(\mathbf{x},k)=\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{j}$.
Then

(56)
$\widetilde{u}(\mathbf{x},k)=\overline{u}(\mathbf{x},k)+\mathbf{V}(\mathbf{x},k)\cdot\mathbf{j,}$
(57)
$\left\|\mathbf{V}(\mathbf{x},k)\cdot\mathbf{j}\right\|_{C\left({D_{+}(\mathbf{\nu},R)}\right)}\leq\sigma,$

which means that the function $\widetilde{u}(\mathbf{x},k)$ approximates well
the solution of the Helmholtz equation (53) with the above incident plane wave
and corresponding radiation conditions.

Next, since
$\mathbf{j\cdot\nu}=\mathbf{j}\cdot(\mathbf{\nu}\times\mathbf{j})=0,$ then
(64) implies that

$\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{\nu=V}(\mathbf{x},k)\cdot\mathbf{\nu,}$
$\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot(\mathbf{\nu}\times\mathbf{j})=\mathbf{V}(\mathbf{x},k)\cdot(\mathbf{\nu}\times\mathbf{j}).$

Hence, components $\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{\nu}$ and
$\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot(\mathbf{\nu}\times\mathbf{j})$ are
sufficiently small,

(58)
$\left\|\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{\nu}\right\|_{C\left({D_{+}(\mathbf{\nu},R)}\right)}\leq\sigma,\left\|\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{\nu}\right\|_{C\left({D_{+}(\mathbf{\nu},R)}\right)}\leq\sigma.$

Thus, it follows from (54)-(58) that the component
$\widetilde{u}(\mathbf{x},k)=\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{j}$
of the electric wave field dominates two other components and it is close to
the solution $\overline{u}(\mathbf{x},k)$ of the Helmholtz equation supplied
by the above incident plane wave and radiation conditions.

Remark 2. Since experimental data have noise, then it is sufficient to obtain
a good approximation $\overline{u}(\mathbf{x},k)$ for the component
$\widetilde{u}(\mathbf{x},k)=\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{j}$
of the electric wave field in the
$C\left(\overline{{D_{+}(\mathbf{\nu},R)}}\right)-$norm.

What is left to do is to prove (54), (55). And this is what is done in the
rest of this section.

It is easy to derive a complete analog of formula (51) for the function
$\overline{\mathbf{E}}$. To do this, one should consider the Cauchy problem
for the time dependent analog of (52) and repeat arguments of Theorem 1 for
the case when the term $\nabla(\mathbf{E}\cdot\nabla\ln n^{2}(\mathbf{x}))$ is
neglected in (11). Hence,

(59)
$\overline{\mathbf{E}}(\mathbf{x},k)=\widehat{\mathbf{\alpha}}^{-1}(\mathbf{x})\exp(-ik\varphi(\mathbf{x}))+O\left(\frac{1}{k}\right),\quad
k\rightarrow\infty,~{}\forall\,\mathbf{x}\in D_{+}(\mathbf{\nu},R),$

where

(60)
$\widehat{\mathbf{\alpha}}^{-1}(\mathbf{x})=\frac{\mathbf{j}}{n(\mathbf{x})\sqrt{J(\mathbf{x})}}.$

Theorem 2. _Suppose that the Assumptions hold. Let_ $\eta>0$_be such a
constant that_

(61) $\left\|\nabla
n(\mathbf{x})\right\|_{C\left(\overline{B}\right)}\leq\eta.$

_Then_

$\widetilde{\mathbf{E}}(\mathbf{x},k)-\overline{\mathbf{E}}(\mathbf{x},k)=\mathbf{V}(\mathbf{x},k)+O\left(\frac{1}{k}\right),~{}k\rightarrow\infty,~{}\mathbf{x}\in
D_{+}(\mathbf{\nu},R)$

_and_

(62)
$\begin{array}[]{ll}\displaystyle\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{j}=\frac{\exp(-ik\varphi(\mathbf{x}))}{n(\mathbf{x})\sqrt{J(\mathbf{x})}}+\mathbf{V}(\mathbf{x},k)\cdot\mathbf{j}+O\left(\frac{1}{k}\right),&\\\
\displaystyle\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot\mathbf{\nu}=\mathbf{V}(\mathbf{x},k)\cdot\mathbf{\nu}+O\left(\frac{1}{k}\right),&\\\
\displaystyle\widetilde{\mathbf{E}}(\mathbf{x},k)\cdot(\mathbf{\nu}\times\mathbf{j})=\mathbf{V}(\mathbf{x},k)\cdot(\mathbf{\nu}\times\mathbf{j})+O\left(\frac{1}{k}\right),&k\rightarrow\infty,~{}\mathbf{x}\in
D_{+}(\mathbf{\nu},R),\end{array}$

_where_

(63)
$\begin{array}[]{ll}\mathbf{V}(\mathbf{x},k)=(\mathbf{\alpha}^{-1}(\mathbf{x})-\widehat{\mathbf{\alpha}}^{-1}(\mathbf{x}))\exp(-ik(\mathbf{x}\cdot\mathbf{\nu})),&\\\
\displaystyle|\mathbf{V}(\mathbf{x},k)|\leq\frac{1}{\sqrt{J_{0}}}[\exp(\eta
T)-1],&\mathbf{x}\in D_{+}(\mathbf{\nu},R).\end{array}$

_Thus, if the number_ $\eta$_in ( 61) is sufficiently small, then (54) and
(55) hold._

Proof. Using formulae (26) and (51) we obtain the first relation (63).
Formulae (62) follow from (63) as well as from (51). Then, using (51) and
(60), we obtain

$|\mathbf{V}(\mathbf{x},k)|=|\mathbf{\alpha}^{-1}(\mathbf{x})-\widehat{\mathbf{\alpha}}^{-1}(\mathbf{x})|=\frac{1}{n(\mathbf{x})\sqrt{J(\mathbf{x})}}\left|\,\int\limits_{\Gamma(\mathbf{x},\Sigma)}R(\mathbf{x},\mathbf{\xi})ds\right|.$

We now estimate $R(\mathbf{x},\mathbf{\xi})$ for $\mathbf{x}\in
D_{+}(\mathbf{\nu},R)$ and
$\mathbf{\xi}\in\Gamma(\mathbf{x},\Sigma(\mathbf{\nu}))$. Note that
$R(\mathbf{x},\mathbf{\xi})=0$ for $\mathbf{x}\in
D_{+}(\mathbf{\nu},R)\setminus D_{0}(\mathbf{\nu},R)$ since $\nabla
n(\mathbf{x})=0$ and, hence, all $K_{n}(\mathbf{x},\mathbf{\xi})=0$ in this
domain. So, we need estimate $R(\mathbf{x},\mathbf{\xi})$ for $\mathbf{x}\in
D_{0}(\mathbf{\nu},R)$ only.

Introduce the matrix norm for a matrix
$K(\mathbf{x},\mathbf{\xi})=(k_{ij}(\mathbf{x},\mathbf{\xi})_{i,j=1}^{3}$ as

$\|K(\mathbf{x},\mathbf{\xi}\|=\max_{i,j=1,2,3}|k_{ij}(\mathbf{x},\mathbf{\xi})|.$

Let $\mathbf{x}\in D_{0}(\mathbf{\nu},R)$ and $\varphi(\mathbf{x})=s_{0}$ and
$\varphi(\mathbf{\xi})=s$. Then using (41) and formulae for
$K_{n}(\mathbf{x},\mathbf{\xi})$, $n=0,1,2,\ldots$, we obtain

$\begin{array}[]{llll}\|K_{0}(\mathbf{x},\mathbf{\xi})\|\leq\eta,&\mathbf{x}\in
D_{0}(\mathbf{\nu},R),&\mathbf{\xi}\in\Gamma(\mathbf{x},\Sigma(\mathbf{\nu})),&\\\
\|K_{1}(\mathbf{x},\mathbf{\xi})\|\leq\eta^{2}(s_{0}-s),&\mathbf{x}\in
D_{0}(\mathbf{\nu},R),&\mathbf{\xi}\in\Gamma(\mathbf{x},\Sigma(\mathbf{\nu})),&\\\
\displaystyle\|K_{n}(\mathbf{x},\mathbf{\xi})\|\leq\eta^{n+1}\frac{(s_{0}-s)^{n}}{n!},&\mathbf{x}\in
D_{0}(\mathbf{\nu},R),&\mathbf{\xi}\in\Gamma(\mathbf{x},\Sigma(\mathbf{\nu})),&n=2,3,\ldots\end{array}$

Hence

$\|R(\mathbf{x},\mathbf{\xi})\|\leq\eta\exp(\eta(s_{0}-s)),~{}\mathbf{x}\in
D_{0}(\mathbf{\nu},R),$

Since $n^{-1}(\mathbf{x})\leq 1$ and $J(\mathbf{x})\geq J_{0}$, then we arrive
at the estimate:

$|\mathbf{V}(\mathbf{x},k)|\leq\frac{1}{\sqrt{J_{0}}}[\exp(\eta
s_{0})-1]\leq\frac{1}{\sqrt{J_{0}}}[\exp(\eta T)-1],~{}\mathbf{x}\in
D_{0}(\mathbf{\nu},R).$

Thus, we obtain the estimate in the second line of (63)_._ This estimate
concludes the proof. $\square$

##  4 Relevance to Experimental Results of [6, 7, 8]

We now explain why Theorem 2 at least partially justifies the validity of
modeling of the propagation of electromagnetic waves in the frequency domain
by the single Helmholtz equation (53) in the works of the second author with
coauthors on experimental data [6, 7, 8]. We say “at least partially” because
a completely precise explanation is unlikely possible since we deal here with
a sort of a “mathematics-to-physics bridge”.

We recall that accurate reconstruction results were obtained in [6, 7, 8] when
solving coefficient inverse problems. Experimental data in these references
were collected for the cases when rather small inclusions mimicking land mines
and improvised explosive devices were embedded in an otherwise uniform
background (dry sand). The dielectric constant was not changing within such an
inclusion, although this was not an assumption in reconstruction algorithms.
Therefore,

$\nabla n(\mathbf{x})=\left\\{\begin{array}[]{c}0\text{ within such an
inclusion,}\\\ 0\text{ in the background.}\end{array}\right.$

The question remains now about the discontinuity of the function
$n(\mathbf{x})$ at the inclusion/background interface. Since any solution of
an elliptic equation, such as, e.g. (53), is sufficiently smooth outside of
discontinuities of its coefficients [5], then we conjecture that the medium
“percepts” the functions $n(\mathbf{x})$ in those inclusions as a smooth
function with rather non-small values of $\left|\nabla n(\mathbf{x})\right|$
near those interfaces. In fact, this has been observed in computed images of
[6, 7, 8]. Now, since values of $\left|\nabla n(\mathbf{x})\right|$ were not
small only in close proximities of those interfaces and volumes of those
proximities were small, then this means that norms $\left\|\left|\nabla
n(\mathbf{x})\right|\right\|_{L_{2}\left(B\right)}$ were actually small. On
the other hand, since finite differences with relatively small numbers of grid
points were used in [6, 7, 8] to solve inverse problems and since all norms in
a finite dimensional space are equivalent, then the smallness of the discrete
norm $\left\|\left|\nabla n(\mathbf{x})\right|\right\|_{L_{2}\left(B\right)}$
is equivalent to the smallness of the discrete norm $\left\|\left|\nabla
n(\mathbf{x})\right|\right\|_{C\left(\overline{B}\right)},$ which is close to
the smallness assumption imposed in Theorem 2 on the number $\eta$ in (61).
Note that smallness assumptions were not used in algorithms of [6, 7, 8].

Thus, Theorem 2 explains, at least partially, accurate reconstructions in [6,
7, 8].

##  5  Appendix

Proof of the Lemma.

Denote ${E}_{j}^{r}$, $j=1,2,3$, components of vector function
$\mathbf{E}^{r}$. Calculating scalar product of both sides of equation (32),
$2\partial_{t}\mathbf{E}^{r}$, using
$n^{2}(\mathbf{x})=\varepsilon(\mathbf{x})$ and applying the identity

$2\partial_{t}\mathbf{E}^{r}\cdot\Delta\mathbf{E}^{r}=2\mathrm{div}\left(\sum\limits_{j=1}^{3}(\partial_{t}{E}_{j}^{r})\nabla{E}_{j}^{r}\right)-\partial_{t}\left(\sum\limits_{j=1}^{3}|\nabla{E}_{j}^{r}|^{2}\right),$

we obtain:

$\displaystyle\partial_{t}\left(\varepsilon(\mathbf{x})|\partial_{t}\mathbf{E}^{r}|^{2}+\sum\limits_{j=1}^{3}|\nabla{E}_{j}^{r}|^{2}\right)-2\mathrm{div}\left(\sum\limits_{j=1}^{3}(\partial_{t}{E}_{j}^{r})\nabla{E}_{j}^{r}\right)$
(64)
$\displaystyle-2\sum\limits_{i,j=1}^{3}\partial_{t}{E}_{i}^{r}\left[(\partial_{x_{i}}{E}_{j}^{r})\partial_{x_{j}}\varepsilon(\mathbf{x})+{E}_{j}^{r}\partial_{x_{j}x_{i}}\ln\varepsilon(\mathbf{x})\right]$
$\displaystyle=$ $\displaystyle
2(\partial_{t}\mathbf{E}^{r})\cdot\mathbf{F}^{r}.$

Integrating identity (64) over the domain $\mathbb{R}_{t}^{4}$, $t\in(0,T]$
and taking into account that $\mathbf{E}^{r}$ and $\mathbf{F}^{r}$ are
compactly supported in $\mathbb{R}_{T}^{4}$ and the initial zero data, we
arrive to the equality

(65)
$\displaystyle\int\limits_{Y(t,T)}\left(\varepsilon(\mathbf{x})|\partial_{t}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla{E}_{j}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}$
$\displaystyle=$ $\displaystyle
2\int\limits_{\mathbb{R}_{t}^{4}}\sum\limits_{i,j=1}^{3}\left(\partial_{\tau}{E}_{i}^{r}(\mathbf{x},\tau)\right)\left(\partial_{\xi_{i}}{E}_{j}^{r}(\mathbf{x},\tau)\right)\partial_{\xi_{j}}\ln\varepsilon(\mathbf{x})\,d\mathbf{x}d\tau$
$\displaystyle+2\int\limits_{\mathbb{R}_{t}^{4}}\sum\limits_{i,j=1}^{3}\left(\partial_{\tau}{E}_{i}^{r}(\mathbf{x},\tau)\right)\left({E}_{j}^{r}(\mathbf{x},\tau)\right)\partial_{\xi_{i}\xi_{j}}\ln\varepsilon(\mathbf{x})\,d\mathbf{x}d\tau$
$\displaystyle+2\int\limits_{\mathbb{R}_{t}^{4}}\partial_{\tau}\mathbf{E}^{r}(\mathbf{x},\tau)\cdot\mathbf{F}^{r}(\mathbf{x},\tau)\,d\mathbf{x}d\tau.$

Transform this equality using assumption (4), (43), the algebraic inequalities
$2\mathbf{a\cdot b}\leq|\mathbf{a}|^{2}+|\mathbf{b}|^{2}$, we obtain

(66)
$\displaystyle\varepsilon_{0}\int\limits_{Y(t,T)}\left(|\partial_{t}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla{E}_{j}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}$
$\displaystyle\leq$ $\displaystyle
C_{1}\int\limits_{\mathbb{R}_{t}^{4}}\left[|\partial_{\tau}\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}+\sum\limits_{j=1}^{3}|\nabla{E}_{j}^{r}(\mathbf{x},\tau)|^{2}+|\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}\right]\,d\mathbf{x}d\tau$
$\displaystyle+\int\limits_{\mathbb{R}_{t}^{4}}|\mathbf{F}^{r}(\mathbf{x},\tau)\,d\mathbf{x}|^{2}\,d\tau.$

Using the inequality

$|\mathbf{F}^{r}(\mathbf{x},t)|^{2}=\left(\int\limits_{0}^{t}\partial_{\tau}\mathbf{E}^{r}(\mathbf{x},\tau)\,d\tau\right)^{2}\leq
T\int\limits_{0}^{t}|\partial_{\tau}\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}\,d\tau,$

we obtain the more general inequality

(67)
$\displaystyle\int\limits_{Y(t,T)}\left(|\partial_{t}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla{E}_{j}^{r}(\mathbf{x},t)|^{2}+|\mathbf{E}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}$
$\displaystyle\leq$ $\displaystyle
C_{1}\int\limits_{\mathbb{R}_{t}^{4}}\left[|\partial_{\tau}\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}+\sum\limits_{j=1}^{3}|\nabla{E}_{j}^{r}(\mathbf{x},\tau)|^{2}+|\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}\right]\,d\mathbf{x}d\tau$
$\displaystyle+\int\limits_{\mathbb{R}_{t}^{4}}|\mathbf{F}^{r}(\mathbf{x},\tau)|^{2}\,d\mathbf{x}d\tau.$

Applying Gronwall-Bellman to inequality (67), we find

$\int\limits_{Y(t,T)}\left(|\partial_{t}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla{E}_{j}^{r}(\mathbf{x},t)|^{2}+|\mathbf{E}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}\leq\|\mathbf{F}^{r}\|_{\mathbb{R}_{T}^{4}}^{2}\exp(C_{1}T).$

Thus, we have obtained the inequalities:

(68) $\|\mathbf{E}^{r}\|_{H^{1}(Y(t,T))}\leq
C_{1}M,~{}\|\partial_{t}\mathbf{E}^{r}\|_{L^{2}(Y(t,T))}\leq C_{1}M.$

where $M$ is defined in (44).

Differentiating equation (32) $k\leq r$ times with respect to $t$ and then
calculating scalar product of both sides of the resulting equation with the
vector function $2\partial_{t}^{k}\mathbf{E}^{r}$, we obtain relations
(64)-(68) with $\partial_{t}^{k}\mathbf{E}^{r}$ instead $\mathbf{E}^{r}$.
Therefore, the following estimates hold

(69) $\|\partial_{t}^{k}\mathbf{E}^{r}\|_{H^{1}(Y(t,T))}\leq
C_{1}M,~{}\|\partial_{t}^{k+1}\mathbf{E}^{r}\|_{L^{2}(Y(t,T))}\leq
C_{1}M,~{}k\leq r.$

Apply now the mathematical induction method to prove estimate (45). Suppose
that for some $n$, $1<n-1<r-1$, the estimates similar (68), (69) hold:

(70) $\begin{array}[]{ll}\|\mathbf{E}^{r}\|_{H^{n-1}(Y(t,T))}\leq
C_{1}M,~{}\|\partial_{t}\mathbf{E}^{r}\|_{H^{n-2}(Y(t,T))}\leq C_{1}M,&\\\
\|\partial_{t}^{k}\mathbf{E}^{r}\|_{H^{n-1}(Y(t,T))}\leq
C_{1}M,~{}\|\partial_{t}^{k+1}\mathbf{E}^{r}\|_{H^{n-2}(Y(t,T))}\leq
C_{1}M,&k\leq r-(n-2),\end{array}$

and prove that the similar estimates are valid when $n-1$ replaced with $n$.
Denote

$D^{\alpha}=\frac{\partial^{|\alpha|}}{\partial_{x_{1}}^{\alpha_{1}}\partial_{x_{2}}^{\alpha_{2}}\partial_{x_{3}}^{\alpha_{3}}},$

where $\alpha=(\alpha_{1},\alpha_{2},\alpha_{3})$ is a multi index,
$\alpha_{1},\alpha_{2},\alpha_{3}$ are integer nonnegative numbers and
$|\alpha|=\alpha_{1}+\alpha_{2}+\alpha_{3}$. We shall use the Leibnitz formula
for a product of two functions

$D^{\alpha}(uv)=\sum\limits_{\beta\leq\alpha}C_{\alpha}^{\beta}(D^{\beta}u)(D^{\alpha-\beta}v),$

where $\beta=(\beta_{1},\beta_{2},\beta_{3})$,
$C_{\alpha}^{\beta}=C_{\alpha_{1}}^{\beta_{1}}C_{\alpha_{2}}^{\beta_{2}}C_{\alpha_{3}}^{\beta_{3}}$
is product of the binomial coefficients and $\beta\leq\alpha$ means that
$\beta_{1}\leq\alpha_{1}$, $\beta_{2}\leq\alpha_{2}$,
$\beta_{3}\leq\alpha_{3}$. Applying the differential operator $D^{\alpha}$
with $|\alpha|=n$ to equation (32) and using the given above formula, we
obtain

(71)
$\displaystyle\varepsilon(\mathbf{x})\partial_{t}^{2}D^{\alpha}\mathbf{E}^{r}-\Delta
D^{\alpha}\mathbf{E}^{r}+\sum\limits_{\beta\leq\alpha,\,\beta\neq\alpha}C_{\alpha}^{\beta}\left(\partial_{t}^{2}D^{\beta}\mathbf{E}^{r}\right)\left(D^{\alpha-\beta}\varepsilon(\mathbf{x})\right)$
$\displaystyle+\sum\limits_{j=1}^{3}\sum\limits_{\beta\leq\alpha}C_{\alpha}^{\beta}\left[(\nabla
D^{\beta}{E}_{j}^{r}))(\partial_{x_{j}}D^{\alpha-\beta}\varepsilon(\mathbf{x}))+(D^{\beta}{E}_{j}^{r})\nabla(\partial_{x_{j}}D^{\alpha-\beta}\ln\varepsilon(\mathbf{x}))\right]$
$\displaystyle=$ $\displaystyle D^{\alpha}\mathbf{F}^{r}(\mathbf{x},t).$

Calculating scalar product of both sides of equation (71) and
$2\partial_{t}D^{\alpha}\mathbf{E}^{r}$, we obtain the relation similar in the
main part to (64), namely:

(72)
$\displaystyle\partial_{t}\left(\varepsilon|\partial_{t}D^{\alpha}\mathbf{E}^{r}|^{2}+\sum\limits_{j=1}^{3}|\nabla
D^{\alpha}{E}_{j}^{r}|^{2}\right)-2\mathrm{div}\left(\sum\limits_{j=1}^{3}(\partial_{t}D^{\alpha}{E}_{j}^{r})\nabla{E}_{j}^{r}\right)$
$\displaystyle+2\sum\limits_{\beta\leq\alpha,\,\beta\neq\alpha}C_{\alpha}^{\beta}\left[(\partial_{t}D^{\alpha}\mathbf{E}^{r})\cdot\partial_{t}^{2}D^{\beta}\mathbf{E}^{r}\right]\left(D^{\alpha-\beta}\varepsilon(\mathbf{x})\right)$
$\displaystyle-2\sum\limits_{i,j=1}^{3}\sum\limits_{\beta\leq\alpha}C_{\alpha}^{\beta}\left(\partial_{t}D^{\alpha}{E}_{i}^{r}\right)\left[\left(\partial_{x_{i}}D^{\beta}{E}_{j}^{r}\right)\partial_{x_{j}}\ln\varepsilon(\mathbf{x})+\left(D^{\beta}{E}_{j}^{r}\right)\partial_{x_{i}x_{j}}D^{\alpha-\beta}\ln\varepsilon(\mathbf{x})\right]$
$\displaystyle=$ $\displaystyle
2(\partial_{t}D^{\alpha}\mathbf{E}^{r})\cdot\mathbf{F}^{r}.$

Integrating this identity over domain $\mathbb{R}_{t}^{4}$, $t\in(0,T]$, we
arrive to the equality

(73)
$\displaystyle\int\limits_{Y(t,T)}\left(\varepsilon(\mathbf{x})|\partial_{t}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla
D^{\alpha}{E}_{j}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}$ $\displaystyle=$
$\displaystyle-2\sum\limits_{\beta\leq\alpha,\,\beta\neq\alpha}C_{\alpha}^{\beta}\int\limits_{\mathbb{R}_{t}^{4}}\left[(\partial_{t}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},\tau))\cdot\partial_{t}^{2}D^{\beta}\mathbf{E}^{r}(\mathbf{x},\tau)\right]\left(D^{\alpha-\beta}\varepsilon(\mathbf{x})\right)\,d\mathbf{x}d\tau$
$\displaystyle+2\sum\limits_{i,j=1}^{3}\sum\limits_{\beta\leq\alpha}C_{\alpha}^{\beta}\int\limits_{\mathbb{R}_{t}^{4}}\left(\partial_{t}D^{\alpha}{E}_{i}^{r}(\mathbf{x},\tau)\right)\left(\partial_{x_{i}}D^{\beta}{E}_{j}^{r}(\mathbf{x},\tau))\right)\partial_{x_{j}}D^{\alpha-\beta}\ln\varepsilon(\mathbf{x})\,d\mathbf{x}d\tau$
$\displaystyle+2\sum\limits_{i,j=1}^{3}\sum\limits_{\beta\leq\alpha}C_{\alpha}^{\beta}\int\limits_{\mathbb{R}_{t}^{4}}\left(\partial_{t}D^{\alpha}{E}_{i}^{r}(\mathbf{x},\tau)\right)\left(D^{\beta}{E}_{j}^{r}(\mathbf{x},\tau))\right)\partial_{x_{i}x_{j}}D^{\alpha-\beta}\ln\varepsilon(\mathbf{x})\,d\mathbf{x}d\tau$
$\displaystyle+2\int\limits_{\mathbb{R}_{t}^{4}}\partial_{\tau}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},\tau)\cdot
D^{\alpha}\mathbf{F}^{r}(\mathbf{x},\tau)\,d\mathbf{x}d\tau.$

Use now assumption (4), (43) and the inequality $2\mathbf{a\cdot
b}\leq|\mathbf{a}|^{2}+|\mathbf{b}|^{2}$. Then, taking into account that
$\sum_{\beta\leq\alpha}C_{\alpha}^{\beta}\leq 2^{3n}$ and
$C_{\alpha}^{\beta}\leq 2^{n}$, we obtain

(74)
$\displaystyle\varepsilon_{0}\int\limits_{Y(t,T)}\left(|\partial_{t}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla
D^{\alpha}{E}_{j}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}$
$\displaystyle\leq$ $\displaystyle
C_{2}\int\limits_{\mathbb{R}_{t}^{4}}\left[|\partial_{\tau}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}+\sum\limits_{j=1}^{3}|\nabla
D^{\alpha}{E}_{j}^{r}(\mathbf{x},\tau)|^{2}\right]\,d\mathbf{x}d\tau$
$\displaystyle+2^{n}\mu\sum\limits_{\beta\leq\alpha,\,\beta\neq\alpha}\,\int\limits_{\mathbb{R}_{t}^{4}}\left(|\partial_{t}^{2}D^{\beta}\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}\,d\mathbf{x}d\tau+\sum\limits_{j=1}^{3}|\nabla
D^{\beta}{E}_{j}^{r}(\mathbf{x},\tau)|^{2}\right)$
$\displaystyle+\int\limits_{\mathbb{R}_{t}^{4}}|D^{\alpha}\mathbf{F}^{r}(\mathbf{x},\tau)\,d\mathbf{x}|^{2}d\tau,$

where $C_{2}=2^{3n+1}\mu$.

Since the relations $\beta\leq\alpha,\,\beta\neq\alpha$ mean that
$|\beta|\leq|\alpha|-1=n-1$, then by the induction assumption (70), there
exists a positive constant $C_{1}$ such that

$2^{n}\mu\sum\limits_{\beta\leq\alpha,\,\beta\neq\alpha}\,\int\limits_{\mathbb{R}_{t}^{4}}\left(|\partial_{t}^{2}D^{\beta}\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}\,d\mathbf{x}d\tau+\sum\limits_{j=1}^{3}|\nabla
D^{\beta}{E}_{j}^{r}(\mathbf{x},\tau)|^{2}\right)\leq C_{1}M^{2}.$

Then we derive from (74) that

(75)
$\displaystyle\int\limits_{Y(t,T)}\left(|\partial_{t}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla
D^{\alpha}{E}_{j}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}$
$\displaystyle\leq$ $\displaystyle
C_{1}\int\limits_{\mathbb{R}_{t}^{4}}\left[|\partial_{\tau}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},\tau)|^{2}+\sum\limits_{j=1}^{3}|\nabla
D^{\alpha}{E}_{j}^{r}(\mathbf{x},\tau)|^{2}\right]\,d\mathbf{x}d\tau$
$\displaystyle+C_{1}M^{2}.$

Applying the Gronwall’s inequality, we obtain

(76)
$\int\limits_{Y(t,T)}\left(|\partial_{t}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla
D^{\alpha}{E}_{j}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}\leq
C_{1}M^{2},~{}|\alpha|=n.$

Differentiating equation (71) $k\leq r-(n-1)$ times with respect to $t$ and
then calculating scalar product of both sides of the obtained equation and
$2\partial_{t}^{k}\mathbf{E}^{r}$, we obtain relations (71)-(76) with
$\partial_{t}^{k}D^{\alpha}\mathbf{E}^{r}$ instead $D^{\alpha}\mathbf{E}^{r}$.
Therefore, the following estimates hold

$\int\limits_{Y(t,T)}\left(\partial_{t}^{k+1}D^{\alpha}\mathbf{E}^{r}(\mathbf{x},t)|^{2}+\sum\limits_{j=1}^{3}|\nabla\partial_{t}^{k}D^{\alpha}{E}_{j}^{r}(\mathbf{x},t)|^{2}\right)d\mathbf{x}\leq
C_{1}M^{2},~{}k\leq r-(n-1).$

Thus, inequalities (70) hold with $n-1$ replaced by $n$. This justifies the
mathematical induction method and we can set $n=r+2$ in (70). The latter
proves the required inequalities (45). $\Box$

## References

  * [1] L. Beilina and M.V. Klibanov, _Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems_ , Springer, New York, 2012.
  * [2] L. Beilina, _Energy estimates and numerical verification of the stabilized domain decomposition finite element/finite difference approach for the Maxwell’s system in time domain_ , Central European Journal of Mathematics, 11 (2013), pp. 702–733.
  * [3] M. Born and E. Wolf, _Principles of Optics_ , Cambridge University Press, 7th edition, 1999.
  * [4] R. Courant and D. Hilbert, _Methods of Mathematical Physics, Vol. 2: Partial Differential Equations_ , Wiley, 1966.
  * [5] D. Gilbarg and N.S. Trudinger, _Elliptic Partial Differential Equations of Second Order_ , Springer Verlag, Berlin, 1983\.
  * [6] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, A. J. Sullivan, Lam Nguyen and V. N. Astratov, _Convexification and experimental data for a 3D inverse scattering problem with the moving point source_ , Inverse Problems, 36 (2020), 085007.
  * [7] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, A. J. Sullivan, Lam Nguyen and V. N. Astratov, _An inverse problem of a simultaneous reconstruction of the dielectric constant and conductivity from experimental backscattering data_ , Inverse Problems in Science and Engineering, published online, https://doi.org/10.1080/17415977.2020.1802447, 2020.
  * [8] M. V. Klibanov, A. E. Kolesov and D.-L. Nguyen, _Convexification method for a coefficient inverse problem and its performance for experimental backscatter data for buried targets,_ SIAM J. Imaging Sciences _,_ 12 (2019), pp. 576-603.
  * [9] M.V. Klibanov and V.G. Romanov, _Two reconstruction procedures for a 3D phaseless inverse scattering problem for the generalized Helmholtz equation_ , Inverse Problems, 32 (2016), 015005.
  * [10] M.V. Klibanov, D.L. Nguyen and L.H. Nguyen,  _A coefficient inverse problem with a single measurement of phaseless data_ , SIAM Journal on Applied Mathematics _,_ 79 (2019), pp. 1-27.
  * [11] O.A. Ladyzhenskaya, _Boundary Value Problems of Mathematical Physics_ , Springer, New York, 1985.
  * [12] V.G. Romanov, _Phaseless inverse problems that use wave interference_ , Siberian Mathematical Journal, 59 (2018), pp. 494-504.
  * [13] V.G. Romanov, _Plane wave solutions to the equations of electrodynamics in an anisotropic medium_ , Siberian Mathematical Journal, 60 (2019), pp. 661-672.
  * [14] V.G. Romanov, _Investigation Methods for Inverse Problems_ , VSP, Utrecht, 2002.
  * [15] N.T. Thanh, L. Beilina, M.V. Klibanov and M.A. Fiddy, _Imaging of buried objects from experimental backscattering time-dependent measurements using a globally convergent inverse algorithm_ , SIAM J. Imaging Sciences, 8 (2015), pp. 757-786.
  * [16] B.R. Vainberg, _Principles of radiation, limiting absorption and limiting amplitude in the general theory of partial differential equations_ , Russian Math. Surveys, 21 (1966), pp. 115-193.
  * [17] B.R. Vainberg, _Asymptotic Methods in Equations of Mathematical Physics_ , Gordon and Breach Science Publishers, New York, 1989.