theory.tex

%% Things to be optionally fixed:
%% *  fix that f is used as frequency (replace with omega) and also as arbitrary function f(t) 
%% *  ensure more consistent usage of "dispersion curve", "dispersion relation" and "solution of Maxwell equations"
% TODO 2015-01-23: 
% 2) investigate how mode anticrossing relates to the significant spatial dispersion
% 3) find out why in effparam.py the data need to be conjugated, when MEEP seems to use the exp(-iwt) convention
% 4) possibly write into the numeric.tex about using the Debye model, with a scan for the damping coefficient (0.1 omega0 ... 100 omega0)

\section{Electrodynamics of local homogeneous media} %{{{
\subsection{Electromagnetic wave in vacuum}
\paragraph{Maxwell equations}  %{{{
In the realm of classical physics, the electromagnetic phenomena are governed by the \textit{Maxwell equations}\index{Maxwell equations} in the following form.
We assume here that no free charges and no sources of currents are present: 
\begin{equation} \nabla \cdot  \D = 0, \label{eq_me1}\end{equation}  
\begin{equation} \nabla \cdot  \B = 0, \label{eq_me2}\end{equation}  
\begin{equation} \nabla \times \E = -\frac{\partial \B} {\partial t}, \label{eq_me3}\end{equation}  
\begin{equation} \nabla \times \HH =  \frac{\partial \D} {\partial t}, \label{eq_me4}\end{equation}  
where $\E$ and $\HH$ are the electric and magnetic vector fields, and $\D$ and $\B$ are the electric and magnetic displacements,
 respectively. These two pairs of field and displacement are related in a similar way as a force is related to the deformation. % note that Simovski2007 regards the current J as force, and E as the deformation!
 The \textit{constitutive relations}\index{constitutive relations} depend on the properties of the medium the wave propagates in, and in vacuum they take the simplest possible form:
\begin{equation}		\D = \varepsilon_0	\E, \quad\quad\quad						\B = \mu_0			\HH,				 \label{eq_ce}\end{equation}
$\varepsilon_0 \approx 8.85\cdot10^{-12}$ F/m being the \textit{vacuum permittivity}\index{permittivity!of vacuum} and $\mu_0 = 4\pi \cdot 10^{-7} \approx 1.25\cdot10^{-6}$ H/m the \textit{vacuum permeability}\index{permeability!of vacuum}. 
Pages \pageref{starttext}--\pageref{endtext} of this thesis will be concerned with computation, interpretation, and experimental verification of the solutions of Eq. (\ref{eq_me1}--\ref{eq_me4})  for specific choices of constitutive relations.
%}}}
\paragraph{Wave equation in vacuum} %{{{
The pair of first-order differential equations (\ref{eq_me3}, \ref{eq_me4}) can be converted to a single second-order differential equation. To this end, we apply an extra curl operator $\nabla\times$ from the left, and substituting one equation into the another, we obtain two curl operators on the left hand side and two derivatives on the right hand side: 
\label{starttext}
\begin{equation} \nabla\times (\nabla\times \E) = \nabla\times \left(- \frac{\partial \B} {\partial t}\right) = -\mu_0 \frac{\partial}{\partial t} \left(\nabla\times \HH\right) 
= -\mu_0 \frac{\partial^{2} \D}{\partial t^{2}} = -\mu_0 \varepsilon_0 \frac{\partial^{2} \E}{\partial t^{2}}.  \label{eq_elim}\end{equation}
Using the vector calculus identity
\begin{equation} \nabla\times (\nabla\times \E) \equiv \nabla (\nabla \cdot \E) - \nabla^2 \E, \label{eq_rotrot}\end{equation}
we obtain the \textit{wave equation}\index{wave equation!in vacuum} for the electric field in vacuum: 
\begin{equation}  \nabla (\nabla \cdot \E) - \nabla^2 \E = -\mu_0 \varepsilon_0 \frac{\partial^{2} \E}{\partial t^{2}}.  \label{eq_wave}\end{equation}
%Note this holds for each component of electric field independently. -- NOT: the divergence is not independent of other components
Starting with Eq. (\ref{eq_me4}) instead of (\ref{eq_me3}), an analogous result can also be easily obtained for the magnetic field $\HH$.

%}}}
\paragraph{Plane wave} %{{{
The solutions of the linear wave equation (\ref{eq_wave}) can be decomposed into a sum of \textit{harmonic plane waves}\index{wave!harmonic}, where \textit{harmonic} means that the amplitude depends on the time $t$ as a harmonic function (e.g. $\sin(\omega t)$ or $e^{\ii \omega t}$). As a \textit{plane wave}\index{wave!plane wave} we denote a spatial shape of the fields that is a function of a single scalar parameter $\kk\cdot\rr$, where $\kk$ is an arbitrary vector and $\rr$ is the position vector. Assuming the wave propagates with a nonzero constant velocity, it follows that also the spatial dependence of the fields is harmonic. Any other complicated shape of the fields can be decomposed into a linear superposition of more such waves and treated separately \cite{jackson1962book}. 

The electromagnetic field can be described as a complex exponential, i.e. as a superposition of two waves differing by a quarter-period phase shift, one defining the real part, one the imaginary part of the field. In comparison with the straightforward description of a plane wave in terms of a cosine (or sine) function, the complex notation formally simplifies some mathematical operations, namely, it allows one to identify easily the \textit{phase of a wave}\index{phase of a wave} (divided by the imaginary unit) with the exponent. 

We define the electric field as a function of time $t$ and position in space $\rr$, corresponding to a plane wave in the complex notation:
\begin{equation} \E(t, \rr) := \E_0\, e^{\ii\omega t - \ii\kk\cdot\rr} \label{eq_pw}\end{equation}
The plane wave is fully characterised by its \textit{amplitude vector}\index{amplitude vector} $\E_0$, \textit{frequency}\index{frequency} $f$ or \textit{angular frequency}\index{frequency!angular} $\omega = 2\pi f$ and \textit{wave vector}\index{wave vector} $\kk$. Note that no restrictions were put to the amplitude vector $\E_0$ so far, thus Eq. (\ref{eq_pw}) can describe both \textit{transverse}\index{wave!transverse} wave of any polarization, with $\E_0 \perp \kk$, and \textit{longitudinal wave}\index{wave!longitudinal} with $\E_0 || \kk$.  As discussed in the introduction, the time dependence of complex fields with positive evolution in time ($e^{+\ii\omega t}$)  was chosen without any impact on the physical conclusions.

%}}}
\paragraph{Dispersion relations in vacuum} %{{{
Only some combinations of ($\E_0, \omega, \kk$) provide a physical solution of the wave equation (\ref{eq_wave}). In vacuum, the allowed solutions can be obtained by first substituting the differential operators by their equivalents for a particular plane wave:
\begin{equation} \nabla \rightarrow -\ii\kk, \quad\quad\quad 
\frac{\partial} {\partial t} \rightarrow \ii\omega, \label{eq_difftok}\end{equation}
so the wave equation (\ref{eq_wave}) can be modified in the following way:
$$					\nabla (\nabla \cdot \E) - \nabla^2 \E				  =	-\mu_0 \varepsilon_0 \frac{\partial^{2} \E}{\partial t^{2}},  $$
$$				 -\ii\kk (-\ii\kk \cdot \E)  - (-\ii\kk \cdot -\ii\kk) \E = -\mu_0 \varepsilon_0 (\ii\omega)^2 \E, $$
$$   - \kk (\kk \cdot \E)      +          k^2 \E            = +\mu_0 \varepsilon_0 \omega^{2} \E,  $$
\begin{equation}  \TT\E            = \frac{\mu_0 \varepsilon_0 \omega^{2}}{k^2} \E.  \label{eq_wavek}\end{equation}
The linear operation $\TT$ of the left hand side can be geometrically interpreted as taking the transverse component of the field $\E$, that is, perpendicular to the wavevector $\kk$. 
\begin{equation} \TT\E :=  -\frac{\kk (\kk \cdot \E)}{k^2} + \E     \equiv     - \frac{\kk\times(\kk\times\E)}{k^2} \label{eq_transverse1}\end{equation}
Rewriting it explicitly, in this thesis we define the \textit{wavefront projection tensor}\index{tensor!of wavefront projection} that will be useful in the following chapters, too:
\begin{equation} \TT = \frac{1}{k^2} 
\left(\begin{array}{ccc} 
	k_y^2+k_z^2  	& -k_x k_y 		& -k_x k_z \\ 
	-k_y k_x 		& k_x^2+k_z^2	& -k_y k_z \\ 
	-k_z k_x 		& -k_z k_y		& k_x^2+k_y^2
	\end{array} \right) 
\text{ or equivalently, }
(\TT)_{ij} = - \frac{k_i k_j}{k^2} + \delta_{ij} .  \label{eq_transverse2}\end{equation}
The solutions of Eq. (\ref{eq_wavek}) for a harmonic plane wave in vacuum can be divided into two groups: 
\begin{enumerate}
 \item{\textit{Transverse electromagnetic waves}, with the electric field and wave vector being perpendicular, i.e. $(\kk \cdot \E) = 0$ and thus $\TT\E = \E$. Therefore, the \textit{dispersion relation}\index{dispersion!in vacuum} for a transverse plane wave in vacuum is linear:
\begin{equation} k~= \sqrt{\mu_0 \varepsilon_0}\; \omega = \frac{\omega}{c}, \label{eq_dispeq_vac}\end{equation}
	where we defined the \textit{light velocity}\index{velocity!of light in vacuum} $c := \frac{1}{\sqrt{\mu_0 \varepsilon_0}}$. The corresponding solution is plotted in Fig. \ref{fg_dcvac}, independent of the orientation of the wavevector $\kk$. As an unanimously accepted convention, the pseudovector $\HH$ is chosen to be oriented so that the perpendicular \textit{right-handed triplet}\index{right-handed triplet} $\E,\HH,\kk$ can be easily indicated with the thumb, index and middle finger of the right hand (cf. Fig. \ref{fg_separating_wave}). Vacuum is thus an example of the broad group of \textit{right-handed}\index{medium!right-handed} media.
} 
 \item{\textit{Longitudinal electromagnetic waves}, with $\kk\,|| \pm \E$ and $\TT\E = 0$, require the right side of Eq. (\ref{eq_wavek}) to be zero. In vacuum, there is no such a solution, except for a homogeneous static electric field ($k = 0, \omega = 0$), but they will be shown to exist in dispersive media.} 
 \end{enumerate}

 \begin{figure}[t] \caption{Left panel shows the permittivity $\epsrl'$ and permeability $\murl'$ of vacuum being equal to 1; as a result, the dispersion curve for a transverse wave in vacuum forms a straight line in the right panel. As in the case of similar plots later, both axes of the right panel are normalized to an arbitrary frequency of $\omega_0$, because the characteristic curve shapes are independent of the scale.} \label{fg_dcvac} \centering 
	\includegraphics[width=\textwidth]{img/dispersion_landau_lifshitz/dispersion_vacuum.pdf}
\end{figure}
%}}}

\subsection{Local response of media to the electromagnetic field} \label{loc_response_of_media}
\paragraph{Local response definition} \label{subsection_local_resp} %{{{
In the whole chapter, we expect the medium properties to be time-invariant, linear, and homogeneous (i.e. independent of time, field amplitude and position in space, respectively). 
In this section, we focus on the special case when the response of the medium in point $\rr$ is not influenced by the electric field in any other point $\brho \neq \rr$. The medium is then said to be \textit{local}\index{medium!local}. 
For most media found in nature, this approximation is very close to reality, except for optics on the nanometric scale \cite{wubs2013nonlocal}. Hence, most electrodynamics textbooks omit the discussion of nonlocal effects. However, using the local theory to a Bloch's wave propagating through periodic structures has a priori no justification and may lead to completely wrong results. The local theory presented here will serve as a basis for the nonlocal theory developed in following chapters.
% Applying the more familiar local electrodynamics to these structures is a mistake common to many papers. While the local theory is mathematically consistent, it predicts unusual spectral shapes and does not allow their further physical interpretation.

When an electric field $\E$ is applied, the medium responds by a change of the electric displacement $\D$ in a characteristic way. 
The immediate linear relation between $\E$ and $\D$ observed in vacuum, appearing at the right side of the constitutive equations (\ref{eq_ce}), remains unchanged, but the response of the matter complements it with a new term called \textit{electric polarisation}\index{electric polarisation}. The polarisation is not instantaneous, so it is generally expressed as a convolution of the \textit{response function}\index{dispersion!in dielectrics}\index{response function|see dispersion} $\chi_e$ with the values of electric field in the previous time $\tau$:
\begin{equation} \D(t,\rr) = \varepsilon_0 \E(t,\rr) + \varepsilon_0 \int_{-\infty}^{t}\chi_e(t-\tau)\, \E(\tau,\rr)\,\mbox{d}\tau. \label{eq_loc_chi_convol}\end{equation}
Let us assume that a harmonic plane wave propagates through the medium, so $\E(t, \rr) := \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr}$, as given by Eq. (\ref{eq_pw}). This can be inserted in the above equation:
\begin{equation} \D(t,\rr) = \varepsilon_0 \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr} + \varepsilon_0 \int_{-\infty}^{t} \chi_e(t-\tau) \, \E_0 \, e^{\ii\omega \tau - \ii\kk\cdot\rr} \,\mbox{d}\tau. \label{eq_chi_convol_harm}\end{equation}
	Substituting $T:=t-\tau$, the exponent can be separated into two parts: one of which factors out from the integral, and one which turns the convolution into a temporal Fourier transform\index{transform!Fourier} of the medium response:
$$				 \D(t,\rr) = \varepsilon_0 \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr} + \varepsilon_0 \int_{-\infty}^{0} \chi_e(T) \, \E_0 \, e^{\ii\omega (t - T) - \ii\kk\cdot\rr} \,\mbox{d}T,$$
$$				 \D(t,\rr) = \varepsilon_0 \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr} + \varepsilon_0 \left( \int_{-\infty}^{0} \chi_e(T)  \, e^{-\ii\omega T}\,\mbox{d}T  \right) \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr}.$$
This can be viewed as an application of the convolution theorem: convolution in the time domain is equivalent to multiplication in the frequency domain.  %% FIXME normalisation by 2pi?
Consequently, we may introduce the local \textit{relative permittivity}\index{permittivity!relative of a dielectric}, or also \textit{dielectric function}, $\epsrl(\omega)$ as a function of frequency. It is a property of the medium that determines how strongly it develops the electric displacement $\D$ in response to the electric field $\E$ of a harmonic wave.
From Eq. (\ref{eq_ce}) it is clear that in vacuum,  $\epsrl = 1$. In an isotropic medium, the permittivity may be obtained as % \E(t, \rr) ??=  1 + \chi_e(\omega) = 
\begin{equation}  \epsrl(\omega) = \left.\frac{D(t,\rr)}{\varepsilon_0 E(t,\rr)} \right|_{\E(t, \rr) := \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr}} = 1 + \int_{-\infty}^{0} e^{-\ii\omega T} \,\chi_e(T) \,\mbox{d}T. \label{eq_eps_loc}\end{equation}

In general, the function $\epsrl(\omega)$ is complex to account for possible phase delay between the harmonic driving field and the medium response. Phase delay different than $0$ (or integer multiple of $\pi$) corresponds to energy being dissipated in the medium. In anisotropic media, the division in Eq. (\ref{eq_eps_loc}), as well as in Eq. (\ref{eq_mu_loc}) has to be replaced by solving a set of linear equations, leading to a tensor of permittivity.
Let us note that the definition of the Fourier transform is never subject to the sign change, even when the convention of $e^{\ii \omega t}$ is exchanged for $e^{-\ii \omega t}$.
%}}}
\paragraph{Response of a harmonic oscillator} \label{chap_lorentzmedia} %{{{
\begin{figure}[t] \caption{\textbf{(a)} Illustration of how a simplified medium responds to an electric field impulse in the shape of Dirac delta function  $E(t) = \delta(t)$. The response is composed from an instantaneous part from vacuum, $\delta(t)$, and from a delayed ringdown of one damped harmonic oscillator, described by $\chi_e(t) := 2\pi \sin(2\pi t)\,e^{-x/2}$; \textbf{(b)} The corresponding local permittivity $\epsrl(\omega)$, computed by the Fourier transform of the response. Note that the imaginary part of permittivity is negative due to the fact that the material is lossy (it dissipates energy) and that the $e^{\ii\omega t}$ convention is used.} \label{fg_oscillator_spectrum} \centering 
	\includegraphics[width=\textwidth]{img/oscillator_spectrum.pdf}
\end{figure}
The response function $\chi_e(T, \mathbf{R})$ of usual media is composed of different phenomena.  Each of them may react on a different time scale, thus the response of the medium usually has a relatively complicated shape in the time domain.  However, to a reasonable degree of approximation, each of the contributions can be treated separately, as it is demonstrated in the following.

Linear physical systems with an inertial load, friction force and a restoring force are known as \textit{damped harmonic oscillators}\index{oscillator!damped harmonic}. Under weak fields, this theory applies well to the electrons elastically bound to an atomic nucleus, as well as to the atoms elastically bound to their equilibrium positions in the lattice. The molecular rotations can also be modelled as (possibly overdamped) harmonic oscillators. Even free electrons in conductive media can fit into the model of a harmonic oscillator provided the restoring force is set to nearly zero. The response of a harmonic oscillator is easy to describe both in time and frequency domains, even without the explicit use of the Fourier transform from Eq. (\ref{eq_eps_loc}). The harmonic oscillator model thus becomes a convenient starting point to approximate the response of materials.

A damped harmonic oscillator is described by a second-order differential equation:
\begin{equation} \alpha \frac{\partial^{2} x(t)}{\partial t^{2}} + \beta\frac{\partial x(t)}{\partial t} + \zeta x(t) = f(t). \label{eq_harm_osc}\end{equation}
Provided the driving term on the right hand side  is harmonic, $f(t) = e^{+\ii \omega t}$, the system response is also a harmonic function, $x(t) = \chi(\omega) e^{+\ii \omega t}$. The differential equation (\ref{eq_harm_osc}) can be easily solved to show that the complex amplitude of the driven oscillations, $\chi(\omega)$, depends on the angular frequency and on the positive real parameters $\alpha$, $\beta$ and $\zeta$ in the following way:
\begin{equation} \chi(\omega) \equiv \frac{x(t)}{f(t)} = \frac{1}{\zeta - \alpha\omega^{2} + \ii\omega\beta}  = \frac{\alpha^{-1}}{\frac{\zeta}{\alpha}-\omega^{2} + \ii\omega\frac{\beta}{\alpha}}. \label{eq_harm_osc_result}\end{equation}
The physical meaning of $\alpha$, $\beta$ and $\zeta$ is of little importance in this text, but without loss of generality the result of Eq. (\ref{eq_harm_osc_result}) can be rewritten into
\begin{equation} \chi(\omega) = \frac{F}{\omega_0^{2}-\omega^{2} + \ii\omega\gamma}, \label{eq_harm_osc_rewritten}\end{equation}
where the physical interpretation of the three (real and positive) parameters is as follows:
\begin{itemize}
 \item{$\omega_0 = \sqrt{\zeta/\alpha}$ is the angular \textit{resonant frequency}\index{frequency!of oscillator resonance}, at which the response is purely imaginary and usually its modulus $|\chi(\omega=\omega_0)|$ is near its maximum.} 
 \item{$\gamma = \beta/\alpha$ is the \textit{damping rate}\index{oscillator!damping rate}. In time domain, it determines the time constant of exponential amplitude decay. } 
 \item{$F = \alpha^{-1}$ is the \textit{oscillator strength}\index{oscillator!strength}, determining the amplitude of the response function.}
 \end{itemize}


\paragraph{Permittivity of Lorentz media} Within the approximation of relatively weak fields, the oscillators act independently of each other.
The response of usual media in frequency domain can thus be decomposed with acceptable precision into a sum of $Q$ independent harmonic oscillators, each $q$-th oscillator having the angular resonant frequency $\omega_{0q}$, damping rate $\gamma_q$ and strength $F_q$.
The permittivity function of the material is a solution of the differential equation of a damped harmonic oscillator, driven by a harmonic source:
\begin{equation} \epsrl(\omega) = 1 + \sum_{q=1}^Q \frac{F_q}{\omega_{0q}^2 - \omega^2 + \ii\omega\gamma_q}. \label{eq_lorentz_eps}\end{equation} 
Advancing from the general formulation in Eq. (\ref{eq_eps_loc}) to the Lorentz oscillator model in Eq. (\ref{eq_lorentz_eps}) is of great importance for theoretical interpretation of the material response, and it has also become a framework for description of periodic structures even in the presence of spatial dispersion.  %% clumsy - remove this sentence?
A sample time- and frequency-domain response of a medium with one harmonic oscillator is shown in Fig. \ref{fg_oscillator_spectrum}.
% Note this is also the way how one communicates the material definition to the numerical simulation software, as described later (in Chap. \ref{chap_fdtd}). 

One can see that each oscillator increases the real part of the low-frequency permittivity, but in the high-frequency limit the contribution of the oscillator vanishes. This can be understood intuitively as that at low frequencies $\omega \ll \omega_0$, the system reacts fast enough to simultaneously follow the driving force, whereas at very high frequencies $\omega \gg \omega_0$, the system cannot respond to the driving force.

The contribution of one oscillator to the low-frequency permittivity $\Delta\epsrl(0)$, is inversely proportional to the oscillator restoring force, which links it to the inverse square of the resonant frequency:
\begin{equation} \Delta \epsrl'(\omega\rightarrow0) = \frac{F}{\omega_0^{2}}.  \label{eq_delta_eps} \end{equation}
A more detailed treatment of the permittivity spectra $\epsrl(\omega)$ may be found in many textbooks, e.g. \cite[p. 454]{klingshirn2007semiconductor}, \cite{dresselhaus1966optical}. 

We will return to the Lorentz oscillator model also in the Chapter \ref{def_of_mat}, where it appears to be essential for realistic definition of materials for accurate numerical simulations. The chapter also describes how a \textit{Debye relaxator}\index{Debye relaxator} (e.g. from overdamped molecular rotation) and the \textit{Drude term}\textit{Drude term} (from unbound motion of free charges) can be represented using correct parameters of an oscillator.
%}}}
\paragraph{Permeability}  %{{{ 
In a manner very similar to the above derivation of the local permittivity, the \textit{local permeability}\index{permeability!local} can be introduced by means of a response of the (isotropic) medium to the magnetic field:
\begin{equation} \murl(\omega) = 1 + \chi_m(\omega) = \left.\frac{B(t,\rr)}{\mu_0 H(t,\rr)} \right|_{HH(t, \rr) := \HH_0 \, e^{\ii\omega t - \ii\kk\cdot\rr}} = 1 + \int_{-\infty}^{0} e^{-\ii\omega T} \,\chi_m(T) \,\mbox{d}T, \label{eq_mu_loc}\end{equation}
where $\chi_m(T)$ is the \textit{magnetic response function}\index{response function!magnetic} of the medium and $\HH_0$ is the amplitude of the magnetic field. This obviously results in an analogous expression for the local permeability in the frequency domain:
\begin{equation} \murl(\omega) = 1 + \sum_{q=1}^Q \frac{F_q}{\omega_{0q}^2 - \omega^2 + i\omega\gamma_q}, \label{eq_lorentz_mu}\end{equation} 
where formally the same notation was used as in  Eq. (\ref{eq_lorentz_eps}): $\omega_{0q}$, $\gamma_q$ and $F_q$ are here the magnetic oscillator's angular frequency, damping frequency and strength. Unlike the electric response, most ordinary media have either almost no response to the magnetic field or their response is limited to low frequencies.

%}}}
\paragraph{Kramers-Kronig relations in local media}%{{{
Causality prevents any medium from reacting to the future electric (or magnetic) field, so the integration in Eq. (\ref{eq_loc_chi_convol}) goes up to the current time only, $\tau \in (-\infty, t)$. The response of the medium to a real-valued field must moreover be also real, no matter that the computations are often done with complex field amplitude [Eq. (\ref{eq_pw})] for sake of convenience. 

Thus, the basic physical laws impose relatively strict constraints on the time-domain response function $f(t)$, which translate into other constraints for the possible shape of the response in frequency domain $F(\omega)$. The intuitive physical derivation is based on the fact that any time-domain response function can be trivially separated into its odd and even parts, as shown in Fig. \ref{fg_kk}. 
\begin{equation}f(t) = f_{odd}(t) + f_{even}(t) = -f_{odd}(-t) + f_{even}(-t) \label{eq_odd_even_decomp}\end{equation}
\begin{figure}[t] \caption{Illustration of how a real causal function $f(t)$ can be decomposed into the odd and even parts, which then yield pure imaginary and pure real functions in the spectrum, respectively. Mathematically this is expressed in Eqs. (\ref{eq_odd_even_decomp}--\ref{eq_kkresult}).} \label{fg_kk} \centering 
	\includegraphics[width=\textwidth]{img/Kramers_Kronig_plot/kk.pdf}
\end{figure}
The Fourier transform of a real odd function is an imaginary function:
%F_{odd}(\omega)=& \int_{-\infty}^{\infty} e^{-\ii \omega t} f(t) \mbox{d}t = \\
		 %=&  \frac{1}{2} \int_{0}^{\infty} e^{-\ii \omega   t } [ f_{odd}(t) + f_{even}(t)] \mbox{d}t 
		   %+ \frac{1}{2} \int_{0}^{\infty} e^{-\ii \omega (-t)} [-f_{odd}(t) + f_{even}(t)] \mbox{d}t  \\
		 %=&  \int_{0}^{\infty} \frac{e^{-\ii \omega t}-e^{+\ii \omega t}}{2} f_{odd}(t) \mbox{d}t 
		   %+ \int_{0}^{\infty} \frac{e^{-\ii \omega t}-e^{+\ii \omega t}}{2} f_{odd}(t) \mbox{d}t  = \\
		 %=&\int_{-\infty}^{\infty} e^{-\ii \omega t} f(t) \mbox{d}t  = \\
		 %=&\int_{-\infty}^{\infty} e^{-\ii \omega t} f(t) \mbox{d}t  = \\
\begin{equation} 
\begin{split} 
F_{odd}(\omega)=& \int_{-\infty}^{+\infty} e^{-\ii \omega t} f_{odd}(t) \,\mbox{d}t = 
		 \int_{-\infty}^{+\infty} \frac{e^{-\ii \omega   t }}{2} f_{odd}(t) \,\mbox{d}t 
		   +  \int_{-\infty}^{+\infty} \frac{e^{-\ii \omega (-t)}}{2} [-f_{odd}(-t)] \,\mbox{d}t  \\
		 =&   \int_{-\infty}^{+\infty} \frac{e^{-\ii \omega t}-e^{+\ii \omega t}}{2} f_{odd}(t) \,\mbox{d}t 
		 = -\ii \underbrace{\int_{-\infty}^{+\infty} \sin(\omega t) \, f_{odd}(t) \,\mbox{d}t}_{\mbox{$\in \mathbb{R}$}},
\end{split} 
\label{eq_kkF}\end{equation}
whereas the Fourier transform of a real even function yields a real function:
\begin{equation} 
\begin{split} 
F_{even}(\omega)= \ldots =   \int_{-\infty}^{+\infty} \frac{e^{-\ii \omega t}+e^{+\ii \omega t}}{2} f_{even}(t) \,\mbox{d}t 
		 = \underbrace{\int_{-\infty}^{+\infty} \cos(\omega t) \, f_{even}(t) \,\mbox{d}t}_{\mbox{$\in \mathbb{R}$}}.
\end{split} 
%% is the following is WRONG?? see http://www.thefouriertransform.com/pairs/step.php#signum,  ho3.pdf
\label{eq_kkFeven}\end{equation}
The odd and even components of the time-domain response function correspond to the imaginary and real part of the response spectrum, respectively:
\begin{equation} F_{odd}(\omega) + F_{even}(\omega) = F(\omega), \text{ where } F_{odd}(\omega) = F''(\omega) \text{ and } F_{even}(\omega) = F'(\omega).  \label{eq_kkFodd}\end{equation}
At the same time, $f_{even}(t)$ and $f_{odd}(t)$ are related to each other by having the opposite sign for $t<0$ and the same sign for $t>0$, that is 
\begin{equation} f_{even}(t) = \mbox{sign}(t)\,f_{odd}(t). \label{eq_kkfeven}\end{equation}
The multiplication in the time domain translates into a convolution in the frequency domain
\begin{equation} 
F'(\omega) = \int_{-\infty}^{+\infty}  \frac{-2\ii}{\omega - \Omega} F''(\omega) \,\mbox{d}\Omega  \equiv  \left[\frac{-2\ii}{\omega}\right]\,\ast\,F''(\omega),
\label{eq_kkresult}\end{equation} 
where we used the knowledge that the $-2\ii/\omega$ function is the Fourier transform of $\mbox{sign}(t)$. Convolution with this function is also known as the \textit{Hilbert transform}\index{Hilbert transform}. % tODO cite

Obviously, Eq. (\ref{eq_kkfeven}) can also be converted to $f_{odd}(t) = \mbox{sign}(t)\,f_{even}(t)$, thus the relation between the real and imaginary parts of the response spectrum also holds when $F_{odd}(\omega)$ and $F_{even}(\omega)$ are exchanged in Eq. (\ref{eq_kkresult}). 

A related mathematical proof of the Kramers-Kronig relations can be derived from the analyticity of the response function in the complex plane of frequency. \cite[p. 125]{klingshirn2007semiconductor}

% maybe:: why do K-K relations forbid to choose arbitrary phase difference over an unit cell? 
%     ... ie. Why the N_eff can be unambiguously extracted from spectra?
% maybe: show how it dictates the "constant area of imaginary part " when resonance Q is changed
%}}}

\subsection{Dispersion relations in local Lorentz media} \label{disp_rel_local_media}
\paragraph{Lower and upper polariton branches of transverse waves} %{{{
Returning to the derivation of dispersion relations, we start with modifying the constitutive relations (\ref{eq_ce}) to a plane wave propagating in a medium:
\begin{equation}		\D := \varepsilon_0\epsrl(\omega)	\E, \quad\quad\quad						\B := \mu_0\mu_r(\omega)		\HH,				 \label{eq_ce2}\end{equation}
with the relative permittivity $\epsrl(\omega)$ and permeability $\mu_r(\omega)$ being two dimensionless functions of frequency, defined in Eq. (\ref{eq_lorentz_eps}, \ref{eq_lorentz_mu}). 
The wave equation (\ref{eq_wavek}) then changes to
\begin{equation}  \TT\E = \varepsilon_0 \mu_0  \epsrl(\omega) \murl(\omega) \, \frac{ \omega^{2}}{k^2} \E.  \label{eq_wavek_disp}\end{equation}
	For transverse waves in isotropic media, the electric field is perpendicular to the wave vector ($\TT\E = \E$), and Eq. (\ref{eq_wavek_disp}) can be simplified to  
\begin{equation} k(\omega) = \sqrt{\varepsilon_0 \mu_0} \sqrt{\epsrl(\omega) \mu_{r}(\omega)}\;\omega = \sqrt{\epsrl(\omega) \mu_{r}(\omega)}\; \frac{\omega}{c}, \label{eq_dispeq_loc}\end{equation}
with the added frequency-dispersive term responsible for the deviation of the curve in Fig. \ref{fg_dcsimpleel} from the original straight light line in Fig. \ref{fg_dcvac}. 

In the simplest example of a single electric resonance with negligible losses, as shown in Fig. \ref{fg_dcsimpleel}, the curve is divided into two separate branches. The \textit{lower polariton branch}\index{polariton branches} is below the angular resonant frequency $\omega_0$ and it is characterized by the Lorentz oscillator being in phase with the electric field. Above $\omega_0$, the dipoles of the Lorentz oscillator can no more follow the electric field and point in the opposite direction. With further increase of the frequency, the permittivity crosses zero at the frequency $\omega_{\text{L}}$ where the \textit{upper polariton branch} starts. In case of a single (or well-isolated) Lorentz oscillator, the difference $\omega_{\text{L}} - \omega_0$ can be computed from the magnitude of the oscillator (using the Lydanne-Sachs-Teller relation) \cite{klingshirn2007semiconductor}.

The same behaviour is observed for a single resonance in the permeability $\mu_r(\omega)$, and will be typical also of the spectra of resonances of macroscopic structures described later.

Note that the formation of upper and lower polariton branches can be also interpreted \cite{landau1984electrodynamics} using the theory of coupled oscillators as the result of \textit{anticrossing}\index{anticrossing}, or also \textit{avoided crossing}\index{avoided crossing|see {anticrossing}}, between the oscillator at the frequency $\omega_0$ (forming a horizontal line) and the photon branch (forming a straight growing light-line). \label{anticrossing}

When losses are present, the lower and upper polariton branches (in Fig. \ref{fg_dcsimpleel}) are connected by a smooth line of \textit{anomalous}\index{dispersion!anomalous} dispersion and very high losses.  

%}}}
\paragraph{Longitudinal waves in dispersive media} %{{{
The wave equation in local dispersive media (\ref{eq_wavek_disp}) 
%\begin{equation} - \kk (\kk \cdot \E) + k^2 \E = \varepsilon_0 \epsrl(\omega) \mu_0 \murl(\omega) \omega^{2} \E, \tag{\ref{eq_wavek_disp} \again} \end{equation}
also allows the existence of longitudinal waves with electric field parallel to the wave vector $\E || \kk$. It was shown earlier that there is no solution for longitudinal waves in vacuum except for a static homogeneous field.

\begin{figure}[t] \caption{Influence of a single resonance in the real part of relative permittivity $\epsrl'(\omega)$ (magenta line in the left panel) to the shape of dispersion curves in the right panel (dashed green line, computed using Eq. (\ref{eq_dispeq_loc}). The lower and upper polariton branches are separated by a spectral region $\omega \in (0.3, 0.65)$, where the wave does not propagate on an appreciable distance.} \label{fg_dcsimpleel} \centering 
	\includegraphics[width=\textwidth]{img/dispersion_landau_lifshitz/dispersion_simple_el.pdf}
\end{figure}
If a local medium is assumed, and the wavenumber $k$ is nonzero, such waves can have a solution with nonzero $\E$ when $\epsrl(\omega)' = 0$, or with nonzero $\HH$ when $\mu_r(\omega)' = 0$. Therefore, the corresponding dispersion curve for a longitudinal wave is a horizontal line at $\omega = \omega_{\text{L}}$, independent of $k$. This would be equivalent to a standing oscillation that maintains the spatial amplitude envelope that was originally excited. 

Different physical phenomena can lead to $\epsrl(\omega)' = 0$, some of which introduce relatively low losses at the corresponding $\omega_{\text{L}}$; namely lattice vibrations in nonconductive crystals or electrons in inductive media (like metals and dilute plasma). 

At the interface of two media with differing signs of permittivity, another type of waves can be excited with an intermediate frequency $\omega < \omega_{\text{L}}$ that cannot propagate in either of the bulk medium. The dispersion curve of such waves is not flat, allowing them to propagate along the interface. Depending on the mechanism, they are known as  \textit{surface plasmons}\index{surface plasmons} or \textit{surface phonon-polaritons}\index{surface phonon-polaritons} \cite[p. 87]{klingshirn2007semiconductor}, respectively. Accordingly, \textit{surface magnons}\index{surface magnons} should be observed at interfaces where $\mu$ changes its sign.
%}}}
\paragraph{Anisotropy of permittivity} \label{par_anisotropy} %{{{
It shall be noted here that the permittivity $\epsrl$ was introduced as a scalar, assuming that the vector of electric field $\E$ and electric induction $\D$ are always parallel:
\begin{equation} \D = \varepsilon_0 \epsrl(\omega) \E \equiv \varepsilon_0  
	\left(\begin{array}{ccc} 
			\epsrl(\omega) & 0 & 0  \\
			0 & \epsrl(\omega) & 0  \\
			0 & 0 &\epsrl(\omega)  
	\end{array} \right) \cdot \E
	\label{eq_ed_orientation}
\end{equation}
In some media with a lower rotational symmetry (such as many crystals, or liquids under static electric field), the medium response depends on the electric or magnetic field direction, and the medium is denoted as \textit{anisotropic}\index{medium!anisotropic}. At the beginning of the chapter we assumed the fields to be relatively weak, which enables one to describe the medium as \textit{linear}\index{medium!linear}. Whatever the linear relation of $\D := \mathcal{L}(\E)$, it must obey the rule
$$\D_1 + \D_2 = \mathcal{L}(\E_1) + \mathcal{L}(\E_2) = \mathcal{L}(\E_1 + \E_2)$$
for any vectors $\E_1, \E_2$. Such a relation can be fully described by a \textit{tensor of permittivity}\index{tensor!of permittivity} 
\begin{equation} \D = \varepsilon_0 \epsrl(\omega) \E \equiv \varepsilon_0  
	\left(\begin{array}{ccc} 
	{\epsrl}_{xx}(\omega) & {\epsrl}_{xy}(\omega) & {\epsrl}_{xz}(\omega)  \\
	{\epsrl}_{yx}(\omega) & {\epsrl}_{yy}(\omega) & {\epsrl}_{yz}(\omega)  \\
	{\epsrl}_{zx}(\omega) & {\epsrl}_{zy}(\omega) & {\epsrl}_{zz}(\omega)  
	\end{array} \right) \cdot \E.
	\label{eq_epstensor}
\end{equation}
An elaborate discussion on all possible forms of this tensor and their physical interpretations can be found e.g. in \cite[pp. 678--686]{born1999book}. An analogous treatment can be applied to the magnetic permeability, though it is less often needed.

%}}}
\paragraph{Dispersion relations in anisotropic local media}  %{{{
If the medium response depends on the direction of the field, the dispersion relations cannot be directly obtained by substitution into the wave equation as in Eq. (\ref{eq_dispeq_loc}).
% a linear transformation of the vector $\E$, which can be expressed by the tensor
%$$ - \kk (\kk \cdot \E) \equiv 
	%-\left(\begin{array}{c} k_x \\ k_y \\ k_z \end{array} \right) \cdot
	%\left[\left(\begin{array}{ccc} k_x & k_y & k_z \end{array} \right) \cdot
	%\left(\begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right)\right] 
		%\equiv
	%\left(\begin{array}{ccc} -k_x^2 & -k_x k_y & -k_x k_z \\ -k_y k_x & -k_y^2 & -k_y k_z \\ -k_z k_x & -k_z k_y & -k_z^2 \end{array} \right) \cdot  
	%\left(\begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right).
	%$$
%$$ k^2 \E \equiv 
%	\left[\left(\begin{array}{ccc} k_x & k_y & k_z \end{array} \right) \cdot 
%	\left(\begin{array}{c} k_x \\ k_y \\ k_z \end{array} \right)\right] \cdot 
%	\left(\begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right),
%	$$
%$$ \mu_0 \varepsilon_0 \omega^{2} \E \equiv 
%	\left(\begin{array}{ccc} k_x^2+k_y^2+k_z^2 & 0 & 0 \\ 0 & k_x^2+k_y^2+k_z^2 & 0 \\ 0 & 0 & k_x^2+k_y^2+k_z^2 \end{array} \right) \cdot  
%	\left(\begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right),
%	$$
The dispersion relation can however still be solved \cite[pp. 667]{born1999book} as a set of three linear algebraic equations based on Eq. (\ref{eq_wavek_disp}).
%\begin{equation} \left|\begin{array}{ccc} 
	%k_y^2+k_z^2 - \mu_0 \varepsilon_0 \omega^{2} 		& -k_x k_y 		& -k_x k_z \\ 
	%-k_y k_x 		& k_x^2+k_z^2-\mu_0\varepsilon_0\omega^{2} 		& -k_y k_z \\ 
	%-k_z k_x 		& -k_z k_y		& k_x^2+k_y^2-\mu_0\varepsilon_0\omega^{2}
	%\end{array} \right| = 0, \label{eq_dispdet}\end{equation}
For simplicity, we assume here that the relative permeability $\murl = 1$; the extension to other cases is possible, too. %% TODO Do epsilon and mu commute, if they are both tensors??
A solution of Eq. (\ref{eq_wavek}) can exist with nonzero $\E$ if and only if the determinant of the set of three linear equations is zero:
\begin{equation} 
\det\left[
\TT -
	\frac{\mu_0 \varepsilon_0 \omega^2}{k^2}
	\left(\begin{array}{ccc} 
	{\epsrl}_{11}(\omega) & {\epsrl}_{21}(\omega) & {\epsrl}_{31}(\omega)  \\
	{\epsrl}_{12}(\omega) & {\epsrl}_{22}(\omega) & {\epsrl}_{32}(\omega)  \\
	{\epsrl}_{13}(\omega) & {\epsrl}_{23}(\omega) & {\epsrl}_{33}(\omega)  
	\end{array} \right) \right] = 0, \label{eq_dispdet}\end{equation}
The search for dispersion curves in an anisotropic medium is thus transformed into finding zeroes of this function of four scalar variables, $k_x$, $k_y$, $k_z$ and $\omega$.
In the most general case, it can be solved by means of numerical algebra software. 
%It would be nice to write about the number of solutions

%}}}
\paragraph{Isofrequency contours}  %{{{
%A resonance in local permittivity (or permeability) creates two polariton branches -- thus the dispersion curves cannot be expressed as one simple function $\omega(k)$ of the wavevector $k$. It would then appear adequate to treat dispersion curves as a function of frequency, $k(\omega)$, but as is shown below in the presence of nonlocal response, neither $k$ is necessarily a simple function of $\omega$. Also 
It is often important to describe the dispersion curves also for different wave angles, which is the best accomplished by plotting the % (real part of)
frequency $\omega$ as the function of \textit{wave vector}\index{wave vector} $\kk$. In three dimensions, this would require mapping a function of three independent variables, $\omega(\kk) = \omega(k_x, k_y, k_z)$. However, in most cases the projection of two selected components of $\kk$ is sufficient to understand all relevant phenomena, and naturally it is much easier to visualize. 
\begin{figure}
  \begin{minipage}[c]{0.5\textwidth}
\hfill
    \includegraphics[width=7cm]{img/dispersion_curves_to_ifc.pdf}
  \end{minipage}
  \begin{minipage}[c]{0.5\textwidth}
    \caption{Relation between a dispersion curve for one photonic branch and the corresponding isofrequency contours. At a selected frequency, all points in the $k_x$-$k_y$ plane are drawn for which a nonzero solution of Maxwell equations exists. For transverse waves in local media, this is equivalent to finding a solution to the dispersion equation (\ref{eq_dispeq_loc}).\\Multiple frequencies can be plotted to describe the frequency dependence. For illustration, an isotropic medium was used, thus all IFCs take the form of a circle. To save space, only one quarter of the circle was plotted here.} \label{fg_ifc_dc}
  \end{minipage}
\end{figure}

Such plots are known as \textit{isofrequency contours}\index{isofrequency contours} (IFC), or also \textit{equifrequency contours}\index{equifrequency contours|see {isofrequency contours}} (EFC), and they allow intuitive geometrical analysis \cite{lindell2001bw} of various phenomena such as light refraction, beam walk-off, total internal reflection etc. 
The relation between a dispersion curve for one photonic branch and the corresponding isofrequency contours in an isotropic medium is illustrated in Fig. \ref{fg_ifc_dc}.

The limitation is, however, that IFC plots do not show the imaginary part and thus are applicable to plot the dispersion in media with no or negligible losses only.  
Each photonic branch also has to be plotted in a separate plot to prevent the contours from overlapping (see the right panel of Fig. \ref{fg_dcsimpleel}). 
Moreover, in every single photonic band, Eq. (\ref{eq_dispdet}) yields two solutions for two possible polarisations of transverse waves \cite[p. 46]{klingshirn2007semiconductor}. The IFCs for these solutions are in general different in anisotropic media, but we always restrict the discussion only to one polarisation in the following figures for simplicity. 
% discuss this further? In uniaxial media, where two diagonal components ...

\begin{figure}[ht] \caption{Isofrequency contours for three different frequencies: \textbf{(a)} IFCs in the isotropic medium are circular and centered in the $k=0$ point. \textbf{(b)} An anisotropic medium with the optical axis perpendicular to the interface, where IFCs take the form of ellipses. \textbf{(c)} A similar case of another anisotropic medium, with the orientation of the optical axis (drawn as the dashed black line) that warrants introducing the index of refraction for the shown direction of the wave vector $\kk$. Note that for clarity, the plots (b) and (c) show the IFCs only for one wave polarisation. %%% todo what is the orientation of the group velocity??
%% FIXME \textbf{(d)} A dispersive \textit{hyperbolic medium}\index{hyperbolic medium}, with one component being a function of frequency and crossing zero.
} \label{fg_ifc} \centering
	\includegraphics[width=.8\textwidth]{img/ifc_freqdispersion.pdf} 
\end{figure}
IFCs are valuable for graphical prediction of wave refraction on interface of two media \cite[p. 118]{shalaev2010book}, \cite{boardman2005negative}. Starting, e.g., by an isotropic medium in Fig. \ref{fg_ifc}a, the wavevectors are known for given frequencies (as three coloured arrows).
The component of wave vector parallel to the interface (chosen as $k_x$ here and indicated by the vertical dashed lines) must be conserved upon refraction. This rule can be intuitively deduced from the continuity of the wave phase at the interface, as well as from the Noether theorem applied to the infinite translational symmetry of the interface. Transferring the vertical dashed lines to an IFC plot of another medium and finding the intersections with an IFC of the corresponding frequency, one can find the new wavevector. 
In Fig. \ref{fg_ifc}, we provide examples of IFCs for one isotropic and two anisotropic media.

%}}}
\paragraph{Index of refraction and its applicability}  %{{{
\label{indexofrefraction}
With the background of the theory developed above, the notion of the \textit{index of refraction}\index{index of refraction}\index{refraction!index of} $N$ can be properly introduced and discussed. 
In the strictest sense, the index of refraction is defined only for \textit{isotropic}\index{medium!isotropic} media. Then it is equivalent to the ratio of the wavenumber $k(\omega) \equiv |\kk|$ to the wavenumber in vacuum at the same frequency: 
\begin{equation} N(\omega) := k(\omega) \, \frac{c}{\omega} \equiv \sqrt{\epsrl(\omega)\murl(\omega)}, \label{eq_dispeqN}\end{equation}
as directly follows from Eq. (\ref{eq_dispeq_loc}).

Starting with a harmonic plane wave with the frequency $\omega$ and the wave vector $\kk^{(1)}$ refracting at an interface of two isotropic media, the projection of $\kk^{(1)}$ to the interface is given as 
$$ k_x = k^{(1)}\,\sin \alpha^{(1)}, $$
where $\alpha^{(1)}$ is the angle between the wavevector in the first medium $\kk^{(1)}$ and the  normal to the interface. In the second medium, a similar relation must apply. Therefore, the wavenumber $k^{(2)}$ may be different and as a result, the angle in the second medium is
\begin{equation} \alpha^{(2)} = \arcsin\left( \frac{k^{(1)}}{k^{(2)}} \,\sin \alpha^{(1)} \right) = \arcsin\left( \frac{N_1}{N_2} \,\sin \alpha^{(1)} \right), \label{eq_snell}\end{equation}
which is known as the \textit{Snell}\index{Snell law}\index{refraction!Snell law} (or also \textit{Snell-Descartes}) law.

The majority of (effective) media discussed in this thesis are, however, more or less anisotropic, and somewhat surprisingly, the notion of \textit{effective index of refraction}\index{index of refraction!effective}\index{refraction!index of, effective} $\Neff$ seems to be widely used for them in the literature anyway. The author thus feels there is a need for a conscientious extension of $\Neff$ for anisotropic cases, too. An extremely loose definition of $\Neff$ could be based on computing the ratio $\Neff(\omega, \kk) := kc/\omega$  for any medium. This could be formally done, but then the wavenumber $k$ would also depend on the direction $\alpha$, and the Snell law in Eq. (\ref{eq_snell}) would become an implicit equation, losing its original purpose of making the computation explicit and notably simple.

The author proposes instead to restrict the term \textit{index of refraction} to all cases where IFCs are perpendicular to $\kk$. This covers not only all isotropic media, but also those cases when the waves propagate close the optical axis of any anisotropic media. Fig. \ref{fg_ifc}c shows an example of such an anisotropic medium with its optical axis oriented parallel to the light wave vector, thus enabling one to approximate the IFC by an osculating circle and to use the Snell law to retrieve the correct angle of refraction. This approximation, however, can be used for a narrow range of angles only.
\begin{equation} \frac{\partial k}{\partial \alpha} \ll k,\label{eq_osculating}\end{equation}

We will show in the following that the prediction of the \textit{beam} propagation is more complex, because it is sensitive to the \textit{curvature of IFCs}\index{isofrequency contours!curvature of}. 

%}}}
%% \paragraph{Isofrequency contours of hyperbolic media}  %{{{ todo? add from Belov2003 
%%   The IFCs in Fig. \ref{fg_ifc}a,b,c were plotted for non-dispersive media, i. e. with $\epsrl$ and $\murl$ independent of frequency. The case of one polariton branch of a dispersive isotropic medium is already shown in Fig. \ref{fg_ifc_dc}, characteristic by non-equidistant IFCs which lead to the refraction angle being dependent on the frequency.
%%   
%%   A more complicated cases are the \textit{hyperbolic media}\index{hyperbolic media}, where, in a certain spectral region, one component of the anisotropic permittivity is positive, while another one is negative. The solution of Eq. (\ref{eq_dispdet}) then results in a hyperboloid, which is either of one sheet or two sheets, depending on the sign of the remaining diagonal component of permittivity tensor. At the frequencies where the dispersive component of permittivity is negative,  
%%   Note that a negative permittivity is always connected with the dispersive behaviour, thus \ref{fg_ifc}d is another example of a dispersive medium. 
%%   
%%   %}}}
\paragraph{Group velocity}  %{{{  
So far, only the propagation and refraction of a plane harmonic wave was discussed, and it was shown that it is determined by the shape of IFC at the given frequency.
Temporal modulation of the wave 
% can be expressed by a complex vector function $\E_0(t)$ that replaces the constant vector $\E_0$ in Eq. (\ref{eq_pw}). This 
is equivalent to the wave being formed by superposition of multiple frequency components in the frequency domain, with their respective amplitudes given by the Fourier transform of the field envelope. The temporal position of the envelope is determined by their \textit{mutual phase difference}, not by the absolute wave phase. 
%redund As a special case, when all components propagate with the same velocity, the envelope does not shift against the wave phase. 

A similar argument can be made with regard to the spatial modulation of the wave. Any wave shape other than the infinite plane wave can be expanded into a superposition of waves with different wavevectors, and the direction of propagation of its spatial envelope is given by mutual phase difference of the constituent waves. %

The velocity vector of the envelope propagation, denoted as the \textit{group velocity}\index{velocity!group v.} $\vg$, can be found \cite{mikki2009electromagnetic} as the gradient of frequency by the wavevector: 
\begin{equation} \vg := \frac{\partial \omega}{\partial \kk} \equiv 
\left(\begin{array}{c} 
	\partial \omega/\partial k_x \\ 
	\partial \omega/\partial k_y \\ 
	\partial \omega/\partial k_z \\ 
	\end{array} \right)
 \label{eq_vg}\end{equation}
In the IFC plot, the group velocity can be found visually as directing always perpendicular to the IFC, with a magnitude being proportional to the density of the IFCs.

If the group velocity is different from the phase velocity, the envelope $\E_0(t)$ is maintained in time, but it continuously shifts against the underlying wave. Thus the actual temporal shape of $\E(t)$ changes upon passing through a dispersive medium.
On the contrary, in vacuum or media with negligible dispersion, each frequency component of the wave acquires an additional phase strictly proportional to its frequency. Then the group velocity coincides with the phase velocity: 
$$\hspace{3cm}\vg = \kk \omega / k^2 \hspace{3cm} \text{(in nondispersive isotropic media)}.$$
%The orientation of $\vg$ and $\kk$ is usually the same, but later it will be shown that it can be otherwise, with peculiar physical consequences.

Usually, in spectral regions near a resonance, also the quadratic or even higher terms in the Taylor expansion of the $\omega(\kk)$ dependence shall be taken into account. This effect is known as the \textit{group velocity dispersion}\index{velocity!dispersion of group v.} as it is equivalent to the group velocity $\vg(\kk)$ being dependent on the wavevector (or, if reformulated, on frequency). It results in a temporal distortion of the wave envelope $\E_0(t)$.  %% TODO ... and what about the wave distortion in space?
%}}}
\paragraph{Beam propagation in anisotropic media}  %{{{  
The refraction of a beam is illustrated in Fig. \ref{fg_ifc2} by means of three slightly different wavevectors it is composed of. For simplicity, a monochromatic wave was used, so brown, black and violet were chosen for the three example wavevectors to prevent confusion with the rainbow-like frequency colour map used in Fig. \ref{fg_ifc}. 

In the first plot, \ref{fg_ifc2}a, the case of an isotropic medium is illustrated. Upon refraction into a general anisotropic medium in Fig. \ref{fg_ifc2}b, each component must maintain its wavevector projection to the interface, thus the wavevectors spread their directions. The resulting beam propagates along the group velocity that is different from the central wavevector; % The components of the resulting beam are   is narrower and does not propagate 
this is also known as the spatial \textit{walk-off}\index{spatial walk-off}.
\begin{figure}[ht] \caption{Examples of IFCs similar to Fig. \ref{fg_ifc}, but now at a single frequency and different wavevectors, corresponding to a monochromatic beam refracting on the interface. \textbf{(a)} The wavevectors in the isotropic medium lie on a circle.
		\textbf{(b)} Generally, anisotropic media exhibit different orientations of the wave vector $\kk$ and the beam propagation. \textbf{(c)} For the special cases of propagation along the optical axis, the beam refracts similarly as in isotropic media.
	\textbf{(d)} In anisotropic media with hyperbolic shape of the IFC, the beam refracts to the opposite direction than the wave vector \cite{lindell2001bw}. } \label{fg_ifc2} \centering 
	\includegraphics[width=\textwidth]{img/ifc_kdispersion_hyp.pdf} 
\end{figure}

Fig. \ref{fg_ifc2}c shows again the special case of the anisotropic medium, where the wave propagates near an optical axis, and thus where the Snell law can still be used to determine the refraction of a plane wave, based the generalized notion of the index of refraction.
However, there is a pitfall if one tries to apply the Snell law for prediction of how the beam will refract. The problem originates from the differential nature of the group velocity definition in Eq. (\ref{eq_vg}). 
%While the magnitude of the wavevector changes only negligibly, with the square of a 
While the phase velocity of a monochromatic wave at the frequency $\omega_0$ does not appreciably change upon a small deviation of the angle $\alpha$ from the optical axis  [cf.  Eq. (\ref{eq_osculating})], the group velocity direction does, because it has a nonzero linear component in its dependence on the wavevector:
$$ \frac{\partial \vg(\kk)}{\partial \alpha} \neq 0.$$
As a result, the group velocity in anisotropic media is always much more sensitive to the angle than the phase velocity, and Snell law does not predict it correctly. It can be regarded as a spatially-dependent manifestation of the group velocity dispersion.

In Fig. \ref{fg_ifc2}d, an extreme case of spatial walk-off is shown on the example of a \textit{hyperbolic}\index{medium!hyperbolic} medium with different signs of the permittivity along different axes. %% TODO study Belov2003 and possibly correct this if hyperbolas are wrong
 The normal to its IFC is nearly perpendicular to the wavevector, and accordingly, the beam refracts in opposite angle with regard to the incident wave. The $k_x$ component of the wave vector is however still maintained. Further geometrical treatise of the beam refraction is in Ref. \cite[p. 46]{klingshirn2007semiconductor}.

%If the amplitude and phase of the wave is not  temporal shape of a wave can be composed as a superposition of different frequency components with certain amplitude and phase. Similarly, any spatial shape of a beam can be composed of infinite waves of different wavevectors. The isofrequency contours, and their dependence on frequency, allow to determine how the temporal and spatial envelope will propagate. 
%Assuming a short pulse for simplicity

%}}}
\paragraph{The sign of the phase and group velocities}  %{{{
In the discussion of refraction both in Fig. \ref{fg_ifc} and \ref{fg_ifc2}, the solution of the vertical wave vector component $k_y$ pointing downwards was always selected without justification. In fact, both upwards and downwards pointing wave vectors provide a valid solution. The choice was made so that in the first medium represented by the leftmost plots (Fig. \ref{fg_ifc}a and \ref{fg_ifc2}a), the incident wave propagates downwards to the interface, and it was assumed that also the refracted wave will propagate downwards, from the interface.

A more complicated case occurs when the wave vector $\kk$ and group velocity $\vg$ point in opposite directions (or, more generally, when they have opposite projections on the normal to the interface).  This manifests itself by the IFC radius \textit{decreasing} with frequency growing, as shown in  Fig. \ref{fg_ifcnr}b,c.
Such a case can indeed occur in natural or artificial media, as described in more detail below. 
\begin{figure}[ht] \caption{Frequency-dependent IFCs of \textbf{(a)} an ordinary medium, \textbf{(b)} an isotropic medium with a negative index of refraction and \textbf{(c)} an anisotropic medium with a negative index of refraction. } \label{fg_ifcnr} \centering  %% TODO comment?
	\includegraphics[width=.9\textwidth]{img/ifc_negrefr.pdf} 
\end{figure}
\begin{figure}[ht] \caption{Wavevector-dependent IFCs for one frequency, for the same examples of media as in Fig. \ref{fg_ifcnr}. At an interface of ordinary medium \textbf{(a)} with an isotropic medium of opposite group and phase velocity \textbf{(b)}, both wavevector and beam direction refract under negative angle. The principle of refractive index can be applied to compute the wave vector of a wave propagating along the optical axis of anisotropic media \textbf{(c)}, but the beam refraction follows a different rule as discussed in the text.} \label{fg_ifcnrk} \centering  
	\includegraphics[width=.9\textwidth]{img/ifc_negrefrk.pdf} 
\end{figure}

Whenever the wave vector $\kk$ is nearly perpendicular to the isofrequency contour the Snell law is still applicable and we then speak of a \textit{negative index of refraction}\index{index of refraction!negative}. 
This is typically, but not exclusively, observed in isotropic media or in media where the wave vector points in a direction close to the optical axis.
The refraction between ordinary media and two examples involving negative-index media are plotted in Fig. \ref{fg_ifcnr}. It should be noted that the negative-index media always exhibits temporal dispersion. Therefore, the refraction of temporally short pulses disperses different frequencies into different angles, as can be seen from the wavevector orientation dependent on frequency in Fig. \ref{fg_ifcnr}b,c.
%% maybe?? the propagation velocities are more complicated 

%}}}
%}}} /section
\section{Electrodynamics of nonlocal homogeneous media} %{{{
\label{chap_nonloc}
\subsection{Nonlocal response} 
\paragraph{Definition of nonlocal media}%{{{
The previous two chapters that concerned local media can be generalized into the theory of \textit{nonlocal}\index{medium!nonlocal} (or, \textit{spatially dispersive}\index{medium!spatially dispersive}\index{dispersion!spatial}) media.  The downside of the spatial-dispersive model of media is that it is more complicated, leading e.g. to an implicit dispersion equation. Its great advantage is however that it provides a necessary level of generality for the description of periodic structures discussed below. 

Some phenomena observed in the frequency spectrum are in fact consequences of the spatial dispersion \cite[p. 359]{landau1984electrodynamics}, which is the case, for example, of the Doppler broadening of resonance lines in gases \cite{makarov2004permittivity}. These phenomena are primarily dependent on the wave vector $\kk$, and the apparent broadening on the frequency axis is only due to the dispersion curve defining a simple relation between the frequency broadening and wavevector broadening.  

In this section, a general class of media is discussed, where the medium response depends explicitly on the history of $\E(\tau, \brho)$ in previous time $\tau < t$ and in all surrounding points $\brho$, and therefore it is described by a spatio-temporal convolution:
\begin{equation} \D(t,\rr) = \varepsilon_0 \E(t,\rr) + \varepsilon_0\int_{V} \int_{-\infty}^{t} \chi_e(t-\tau, \rr-\brho) \, \E(\tau,\brho) \,\mbox{d}\tau \,\mbox{d}^3\brho. \label{eq_chi_convol_nonloc}\end{equation}
In a very similar manner as in the local theory above, we assume that a plane wave $\E(t, \rr) := \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr}$ propagates through the medium. This is without loss of generality, since it is possible to express any wave as a superposition of monochromatic plane waves.
\begin{equation} \D(t,\rr) = \varepsilon_0 \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr} + \varepsilon_0\int_{V} \int_{-\infty}^{t} \chi_e(t-\tau, \rr-\brho) \, \E_0 \, e^{\ii\omega \tau - \ii\kk\cdot\brho} \,\mbox{d}\tau \,\mbox{d}^3\brho. \label{eq_chi_convol_harm_nonloc}\end{equation}
After two  substitutions, $T:=t-\tau$, $\mathbf{R}:=\rr-\brho$, the exponent can again be separated into the original plane wave (which factors out), and a spatio-temporal Fourier transform of the medium response:
$$				 \D(t,\rr) = \varepsilon_0 \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr} + \varepsilon_0\int_{V} \int_{-\infty}^{0} \chi_e(T, \mathbf{R}) \, \E_0 \, e^{\ii\omega (t - T) - \ii\kk\cdot(\rr - \mathbf{R})} \,\mbox{d}T \,\mbox{d}^3\brho,$$
$$				 \D(t,\rr) = \varepsilon_0 \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr} + \varepsilon_0\left( \int_{V} \int_{-\infty}^{0} \chi_e(T, \mathbf{R})  \, e^{-\ii\omega T + \ii\kk\cdot \mathbf{R}}\,\mbox{d}T \,\mbox{d}^3 \mathbf{R} \right) \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr}.$$
The response of the medium to the electric field of any harmonic plane wave can now be expressed as a function of frequency $\omega$ and wave vector $\kk$. It is defined as the ratio between the electric displacement and the electric field:
\begin{equation} \epsrn(\omega, \kk) = \left.\frac{\D(t,\rr)}{\varepsilon_0 \E(t,\rr)} \right|_{\E(t, \rr) := \E_0 \, e^{\ii\omega t - \ii\kk\cdot\rr}}  = 1 + \int_{V} \int_{-\infty}^{0} \chi_e(T, \mathbf{R}) \,e^{-\ii\omega T + \ii\kk\cdot \mathbf{R}} \,\mbox{d}T \,\mbox{d}^3\mathbf{R} \label{eq_eps_nonloc}\end{equation}
Converting the problem from the spatio-temporal domain into the wavenumber-frequency domain allows to express the relation between $\D$ and $\E$ by the \textit{permittivity}\index{permittivity!as a function of frequency} function $\epsrn(\omega, \kk)$ and completely avoid the convolution from Eq. (\ref{eq_chi_convol_harm_nonloc}). Note that both the response function $\chi_e$ and the permittivity $\epsrn$ 
may be either scalar functions, or rank-2 tensor functions; the latter case accounts for possible anisotropy of the medium.

The terms of \textit{nonlocality}\index{nonlocality|see {spatial dispersion}} and of \textit{spatial dispersion} are used interchangeably in the literature. The difference seems to be related to the way one thinks about the medium -- while \textit{nonlocality} is obviously related to the description in the real space [cf. Eq. (\ref{eq_chi_convol_nonloc}), \textit{spatial dispersion}, in contrast, derives from that the response is not a constant function in the reciprocal $\kk$-space. The term \textit{spatial dispersion} is therefore of a slightly narrower meaning, as it implies that an infinite plane wave of a defined wavevector is being considered.

%}}}
\paragraph{Power expansion of the medium parameters} %{{{
Assuming the permittivity $\epsrn(\omega,\kk)$ or permeability $\murn(\omega,\kk)$ are smooth functions varying slowly with $\kk$, we can express them in general as power series \cite[p. 367]{landau1984electrodynamics}:
\begin{equation} 
\left.  \begin{array}{c}
\epsrn(\omega,\kk) = 1 + \chi_e(\omega) + \ii \gamma_{\text{e}}(\omega) \kk + [\alpha_{\text{e}}(\omega) \kk] \kk + \ldots, \medskip  \\
\murn(\omega,\kk) = 1 + \chi_m(\omega) + \ii \gamma_{\text{m}}(\omega) \kk + [\alpha_{\text{m}}(\omega) \kk] \kk + \ldots, 
\end{array} \quad \right\} \quad \text{(in any media)}
\label{eq_epsmusd1}\end{equation} %% maybe change gamma->eta, alpha->phi, ...->psi, ...-xi ? 
where $\chi_{e,m}(\omega)$ are second rank tensors, $\gamma_{e,m}(\omega)$ third rank tensors, $\alpha_{e,m}(\omega)$ fourth rank, and similarly for possible higher orders of expansion. After the corresponding number of matrix multiplication with $\kk$, they all yield rank-2 tensors that add up to form the tensor of permittivity or permeability.
%In this customary symmetric approximation of local media, both the electric permittivity and permeability in Eq. (\ref{eq_eps_loc}) are composed of two terms, one caused by the immediate response of vacuum, and another by the response of matter. In nonlocal media, both functions depend on the wavevector $\kk$, making the solution of Maxwell equations more complicated.

Note the response function for a local medium can be formally derived from its nonlocal formulation by replacing the spatial dependence in Eq. (\ref{eq_chi_convol_nonloc}) by a Dirac delta function as follows: %% is the sencence correct?
\begin{equation} \chi_e(t-\tau, \rr-\brho) \rightarrow \delta^{3}(\rr-\brho) \; \chi_e(t-\tau), \label{eq_loc_chi}\end{equation}
which allows to simplify Eq. (\ref{eq_chi_convol_nonloc}) so that in local media only the temporal convolution has to be computed. Then the higher order terms including $\gamma_{e,m}$ and $\alpha_{e,m}$ in Eqs. (\ref{eq_epsmusd1}) vanish and
\begin{equation} 
\left.  \begin{array}{c}
\epsrn(\omega,\kk) = 1 + \chi_e(\omega) = \epsrl(\omega), \medskip \\
\murn(\omega,\kk) = 1 + \chi_m(\omega) = \murl(\omega). 
\end{array} \quad \right\} \quad \text{(in local media)}
\label{eq_epsmusd2}\end{equation}

%}}}
\paragraph{Magnetic effects can be described by nonlocal permittivity} %{{{
Here, we will follow the approach of Landau and Lifshitz \cite{landau1984electrodynamics} to show that the magnetic response of any medium can be fully expressed by a certain form of permittivity dependence on $\kk$. (For more details, see Refs. \cite{agranovich2006spatial}, \cite[pp. 95-130]{krowne2007book}  and \cite[p. 19]{fietz2011homogenization}.) This leads to introducing new \textit{Landau-Lifshitz}\index{permittivity!Landau-Lifshitz form} permittivity $\epsLL$ and permeability $\muLL$, which are, in general, different from those used in the more customary symmetric model:
$$\epsLL(\omega, \kk) \not\equiv \epsrn(\omega, \kk),\quad\quad\quad 1 = \muLL(\omega, \kk) \not\equiv \murn(\omega, \kk).$$
The Maxwell equations (\ref{eq_me1}-\ref{eq_me4}) however still hold when these new parameters are substituted for the original ones. Requiring the relative permeability to be unity implies a trivial dependence of the magnetic field on the magnetic induction in this model:
$$ \mu_0 \HHsd = \B. $$

The advantage is that the relative magnetic permeability is defined as a mere constant of the magnetic response of vacuum $\muLL(\omega, \kk) := 1$, thus reducing the complexity of the computation compared to the symmetric spatial-dispersive model. %% [as in Eq. (\ref{eq_ce})]
Therefore, the \textit{Landau-Lifshitz} model developed below is also denoted as the $\E\D\B$-model, since it allows the substitution into the Maxwell equations to avoid explicit use of the magnetic field $\HH$.

%}}}
\paragraph{Local medium in the Landau-Lifshitz model} %{{{
An important step towards using the Landau-Lifshitz model is to derive how a classical, local medium with the magnetic response is represented. From the principle of correspondence, all local media can be expressed with this model without any change in the dispersion curves predicted. 

In the Landau-Lifshitz model, the new spatial-dispersive permittivity $\epsLL(\omega, \kk)$ consists of
\begin{enumerate}
 \item{the component caused by the local electric response of matter,} 
 \item{a new component fully accounting for the local \textit{magnetic} response of matter, thanks to a particular shape of its spatial dispersion.}
\end{enumerate}
Later, higher-order expansion terms can be easily added to describe all sorts of the nonlocal response.

Recalling the Maxwell equation (\ref{eq_me4}) that links the magnetic field $\HH$ with the electric induction $\D$, 
\begin{equation} \nabla \times \HH =  \frac{\partial \D} {\partial t}, \tag{\ref{eq_me4} \again} \end{equation}
it is clear that if one defines a new pair of vector fields
\begin{equation} \HHsd = \HH + \frac{\partial\mathbf{X}}{\partial t}, \label{eq_HHsd}\end{equation}
\begin{equation} \Dsd  = \D  + \nabla\times \mathbf{X}, \label{eq_Dsd}\end{equation}
then Eq. (\ref{eq_me4}) maintains exactly the same form with the new fields, for any differentiable vector field $\mathbf{X}$:
\begin{equation} \nabla \times \HHsd = \nabla \times \HH + \left(\nabla\times \frac{\partial\mathbf{X}}{\partial t}\right) = \frac{\partial \D}{\partial t}+ \frac{\partial(\nabla\times \mathbf{X})}{\partial t} =  \frac{\partial \Dsd} {\partial t}, \label{eq_me4sd} \end{equation}
because for reasonably shaped functions the temporal and spatial derivatives commute.

%% todo how this transform is denoted?
With the \textit{gauge freedom}\index{gauge freedom in electrodynamics} of choice of $\mathbf{X}$, we impose the above mentioned requirement that whole magnetic response of the matter is expressed by the constitutive equation for permittivity. Therefore in the spatially dispersive theory, the constitutive equation 
for magnetic induction is defined as the same as in vacuum, cf. Eq. (\ref{eq_ce}):
\begin{equation} \mu_0 \HHsd := \mu_0 \murl(\omega) \HH = \B. \label{eq_mu_sd}\end{equation}
When this equation is rearranged into the form similar to \ref{eq_HHsd}, we obtain a prescription for the sought $\mathbf{X}$: 
$$ \HHsd = \HH + (\murl(\omega) -1)\HH = \HH + \underbrace{\left(\frac{\murl(\omega)-1}{\mu_0\murl(\omega)}\right)\B}_{=:\,\partial\mathbf{X}/\partial t}$$
Without loss of generality, we again restrict the discussion to a plane wave (\ref{eq_pw}), thus the time derivative equals to multiplication by $\ii\omega$.
\begin{equation} \mathbf{X} = \frac{1}{\ii\omega}\left(\frac{\murl(\omega)-1}{\mu_0\murl(\omega)}\right)\B = \frac{1}{\ii\omega\mu_0}\left(1 - \frac{1}{\murl(\omega)}\right)\B. \label{eq_Xsd}\end{equation}
The new electric displacement $\Dsd$ in the  Landau-Lifshitz model, that also accounts for magnetic phenomena, is obtained by substitution of Eq. (\ref{eq_Xsd}) into Eq. (\ref{eq_Dsd}):
\begin{equation} \Dsd := \D - \ii\kk\times \mathbf{X} =  \D - \ii  \frac{1}{\ii\omega\mu_0}\left(1 - \frac{1}{\murl(\omega)}\right) \kk\times \B  \label{eq_Dsd2}\end{equation}
By means of the other Maxwell equation (\ref{eq_me3}), the magnetic induction $\B$ can be substituted by $\kk\times\E / \omega$ to obtain an expression that contains the electric quantities only.
\begin{equation} \Dsd = \D - \frac{1}{\omega^2 \mu_0}\left(1 - \frac{1}{\murl(\omega)}\right) \kk\times(\kk\times \E).  \label{eq_Dsd3}\end{equation} 
Double vector multiplication on the right hand side can be identified with the wavefront projection tensor $\TT$, cf. Eqs. (\ref{eq_rotrot}, \ref{eq_transverse1}):
\begin{equation} \Dsd = \D + \frac{k^2}{\omega^2 \mu_0}\left(1 - \frac{1}{\murl(\omega)}\right) \, \TT \E,  \label{eq_Dsd4}\end{equation} 

%}}}
\paragraph{Tensor form of the Landau-Lifshitz permittivity of local isotropic media}%{{{
%Continuing in the derivation outlined in , 
From Eq. (\ref{eq_Dsd3}) we can derive the tensor form of spatial-dispersive permittivity $\epsLL(\omega,\kk)$:
\begin{equation} 
\left.  \begin{array}{c}
\epsLL(\omega,\kk) = 1 \,+\, \chi_e(\omega) \;\;+\;\; \frac{k^2}{\omega^2 \mu_0}\left(1 - \frac{1}{\murl(\omega)}\right) \, \TT, \medskip \\
\muLL(\omega,\kk) = 1, 
\end{array} \quad \right\} \quad \text{(in local media)}
\label{eq_epsmusd3} \end{equation} 
where we make use of the wavefront projection tensor defined as
\begin{equation} \TT = \frac{1}{k^2} 
\left(\begin{array}{ccc} 
	k_y^2+k_z^2		& -k_x k_y		& -k_x k_z \\ 
	-k_y k_x		& k_x^2+k_z^2	& -k_y k_z \\ 
	-k_z k_x		& -k_z k_y		& k_x^2+k_y^2
	\end{array} \right) 
\text{ or equivalently, }
(\TT)_{ij} = - \frac{k_i k_j}{k^2} + \delta_{ij} . \tag{\ref{eq_transverse2} \again} \end{equation}

%}}}
\paragraph{Other formulations}%{{{
Let us note that the classical approach using symmetric parameters $\epsrl(\omega,\kk),\murl(\omega,\kk)$, and the Landau-Lifshitz approach of gathering all medium properties in the permittivity $\epsLL(\omega,\kk)$ are only two examples of all possible gauge transformations of medium parameters under which the Maxwell equations are invariant. 
Formally, one can distribute the medium electric and magnetic responses between these parameters arbitrarily -- the distribution weight may be even frequency dependent, leading to a physically realistic pair of parameters with custom spectral shape \cite{skaar2014diamagnetism}.
%}}}

\subsection{Dispersion relations in nonlocal homogeneous media} 
\paragraph{Dispersion relation as an implicit equation} %{{{
In the dispersion relation for local media, Eq. (\ref{eq_dispdet}), the wave vector could be found by solving a set of linear equations. In the particular case of \textit{isotropic} local media, the computation has an even simpler form of an explicit computation of one square root, see Eq. (\ref{eq_dispeq_loc}).

An intrinsic issue of spatial-dispersive media is that the permittivity $\epsLL(\omega,\kk)$ is a function of the wavevector $\kk$ on its own. The dispersion equation then becomes an implicit equation similar to Eq. (\ref{eq_dispdet}), but now with $\epsLL$ dependent on $\kk$: 
\begin{equation} 
\det\left[
\TT -
	\frac{\mu_0 \varepsilon_0 \omega^2}{k^2}
	\left(\begin{array}{ccc} 
	{\epsLL}_{11}(\omega,\kk) & {\epsLL}_{21}(\omega,\kk) & {\epsLL}_{31}(\omega,\kk)  \\
	{\epsLL}_{12}(\omega,\kk) & {\epsLL}_{22}(\omega,\kk) & {\epsLL}_{32}(\omega,\kk)  \\
	{\epsLL}_{13}(\omega,\kk) & {\epsLL}_{23}(\omega,\kk) & {\epsLL}_{33}(\omega,\kk)  
	\end{array} \right) \right] = 0. \label{eq_dispdetLL}\end{equation}
If all tensor elements, ${\epsLL}_{ij}(\omega,\kk)$, can be approximated in the form of a polynomial expansion to a low order in $\kk$, the solution can again be found as roots of a characteristic polynomial. The degree of this polynomial may be substantially higher than in Eq. (\ref{eq_dispdet}), the dispersion curves can nonetheless be found by a brute-force numerical search. %% ?? for any ${\epsLL}_{ij}(\omega,\kk)$.
\begin{figure}[t] \caption{Dispersion curves for three different types of local media, comparing the local and Landau-Lifshitz representations of constitutive parameters. Three rows show media with \textbf{(a)} a  resonance in permittivity, \textbf{(b)} in permeability, and \textbf{(c)} in both parameters simultaneously.\\
Left column: Local parameters $\epsrl(\omega)$, $\murl(\omega)$ as functions of frequency. Right column: Dispersion curves were computed either with the local parameters from Eq. (\ref{eq_epsmusd2}), plotted as dashed green lines, or with the Landau-Lifshitz parameters from Eq. (\ref{eq_epsmusd3}), plotted as black line. The Landau-Lifshitz permittivity $\epsLL(\omega,\kk)$ is colour shaded on the background. %% TODO make clearer?
} \label{fg_dcll} \centering  
\textbf{(a)}\\\includegraphics[width=1\textwidth]{img/dispersion_landau_lifshitz/dispersion_ll_el.pdf}    
\textbf{(b)}\\\includegraphics[width=1\textwidth]{img/dispersion_landau_lifshitz/dispersion_ll_mag.pdf}
\textbf{(c)}\\\includegraphics[width=1\textwidth]{img/dispersion_landau_lifshitz/dispersion_ll_elmag.pdf}
\end{figure}
\clearpage

%}}}
\paragraph{Dispersion in a local dielectric in the Landau-Lifshitz model} %{{{
The dispersion curves predicted by the classical and Landau-Lifshitz representations of the constitutive parameters for local media in Eq. (\ref{eq_epsmusd3}) are mathematically identical. This is illustrated in Fig. \ref{fg_dcll} on three examples: that of a medium with an electric resonance (\ref{fg_dcll}a), one with  a magnetic resonance (\ref{fg_dcll}b) and  with both resonances overlapping (\ref{fg_dcll}c). 

The plots in the left column of Fig. \ref{fg_dcll} show the \textit{local} permittivity $\epsrl(\omega)$ and permeability $\murl(\omega)$.
Its right column features a thin black contour connecting all points where the dispersion equation (\ref{eq_dispdetLL}) was found to hold by a numerical search.
This is complemented by the Landau-Lifshitz permittivity $\epsLL(\omega, k)$ as a background colour map with blue tone for negative values and red for positive ones. 
%As a confirmation of the compatibility of the Landau-Lifshitz model with the classical one,
Additionally, the original green dashed dispersion curve as in Fig. \ref{fg_dcsimpleel} is retained, as computed using the classical approach based on the Eq. (\ref{eq_dispeq_loc}). It can thus be seen that both models always predict exactly the same dispersion.

For a local isotropic medium with a single electric resonance in Fig. \ref{fg_dcll}a, the plotted curves are identical to Fig. \ref{fg_dcsimpleel}. On the right panel, the Landau-Lifshitz permittivity $\epsLL(\omega, k)$ follows a resonance curve in frequency, but is independent of the wavenumber $k$ (as long as the medium is local). With increasing frequency, the lower polariton branch bends towards higher $k$, as $\epsLL(\omega, k)$ increases towards the resonance at $\omega_0$, then a band of frequency follows where and no solution of the wave equation \cite{klingshirn2007semiconductor}
exists due to $\epsLL(\omega, k)$ being negative, and finally the upper polariton branch starts when $\epsLL(\omega, k)$ crosses zero and becomes positive again.

A local medium with a single \textit{magnetic} resonance, Fig. \ref{fg_dcll}b, is predicted by the symmetric model to exhibit similar dispersion curves. In the right panel of Fig. \ref{fg_dcll}b, the magnetic resonance is represented by the
contribution that grows proportionally to $k^2$: 
\begin{equation} \Dsd = \D + \underbrace{\frac{k^2}{\omega^2 \mu_0}\left(1 - \frac{1}{\murl(\omega)}\right) \, \TT \E}_{\text{magnetic contribution}},      \tag{\ref{eq_Dsd4} \again} \end{equation}
% TODO FK: Nějak to tam nevidím. To jako pro ω/ω0> 1.3? Chtělo by to popsat nějak názorněji.
and unlike the case of the electric resonance, the maximum magnitude of the magnetic contribution is located at the frequency $\omega_{mp}$ where $\murl(\omega_{mp}) \rightarrow 0$. In the plot \ref{fg_dcll}b, this is at $\omega/\omega_0 \approx 1.3$.
% TODO FK: Ani tady se mi nezdá zjevné, kde je tato podmínka splněna.
This shape of $\epsLL(\omega, k)$ causes the lower polariton branch to bend and approach a horizontal asymptote, which is again separated by a band of no allowed wave from the upper polariton branch, starting at $\omega_{mp}$.

Finally, in Fig. \ref{fg_dcll}c, both resonances are combined. The main difference from the cases of isolated resonances occurs in the frequency range where originally $\epsLL(\omega, k)$ was negative and no wave could propagate. %% TODO "Ja to v techobrazech nejak nevidim, nebo presne rerozumim, co tym myslis."
However, the magnetic resonance increases $\epsLL(\omega, k)$ by a term proportional to $k^2$, and consequently a new photonic branch is formed with $\mathrm{d}\omega/\mathrm{d}k < 0$, that is, with opposite group and phase velocities. IFCs for the new band then appear as sketched in Figs. \ref{fg_ifcnr}b and \ref{fg_ifcnrk}b.
% (electric quadrupoles..., [Merlin])
% todo Why magnetic resonances do not contribute to the static permeability?

%% XXX up to here revision of FK 2015-06-28 XXX
%% XXX from  here revision of FK 2015-09-28 XXX

%}}}
\paragraph{Dispersion in nonlocal media} %{{{ 
\begin{figure}[t] \caption{Dispersion curves for two nonlocal media. They differ by the value of the fourth-order expansion coefficient $\alpha(\omega)$, which was plotted with a thin blue line. The left and right columns of plots show the same information as in Fig. \ref{fg_dcll}. The Landau-Lifshitz permittivity $\epsLL(\omega,\kk)$ is colour shaded on the background.  } \label{fg_dcll_nl} \centering  
\textbf{(a)}\\	\includegraphics[width=1\textwidth]{img/dispersion_landau_lifshitz/dispersion_ll_quadrupp.pdf}
\textbf{(b)}\\	\includegraphics[width=1\textwidth]{img/dispersion_landau_lifshitz/dispersion_ll_quadrupn.pdf}
\end{figure}
Further terms in the permittivity $\epsLL(\omega,\kk)$ expansion in Eq. (\ref{eq_epsmusd3}) make the dispersion curves deviate from those predicted for local media. This corresponds to the black contours deviating from the green lines in Fig. \ref{fg_dcll_nl}.
As the simplest example we add one scalar term $\alpha(\omega) k^4$ to $\epsLL(\omega,\kk)$ in Eq. (\ref{eq_epsmusd1}). The shape of $\alpha(\omega)$  was chosen as a weak Lorentz oscillator at the resonant frequency of the ordinary electric response $\chi_e(\omega)$. This term is shown by adding a thin blue line in the left column of Fig. \ref{fg_dcll_nl}.

\label{chap_sd}
The choice of the same resonant frequency for $\chi_e(\omega)$ and $\alpha(\omega)$ is not arbitrary; it is assumed that both terms arise from the same resonance mode that has a field shape different from that of a simple dipole.  The relation of higher-order expansion terms of the Landau-Lifshitz permittivity  to the multipole expansion of the field shape is developed e.g. in Refs. \cite{agranovich2006spatial}, \cite{vinogradov2002form} and \cite{fietz2011homogenization}.

The choices of a positive amplitude (Fig. \ref{fg_dcll_nl}a) and a negative one (Fig. \ref{fg_dcll_nl}b) have typical impacts on the dispersion curves. In the former case, both polariton branches deviate one from another with growing frequency. In the latter case, they come closer to each other with growing $k$. Eventually, in the upper right corner of the right plot in Fig. \ref{fg_dcll_nl}b, they merge into one loop. The author, however, believes this merging may not be observed in nature, and that its occurrence is only due to unrealistic values of the $\alpha(\omega)$ coefficient or the absence of higher-order expansion terms.

%}}}
\paragraph{Existence of additional waves}   %{{{ 
For both cases shown in Fig. \ref{fg_dcll_nl}, it follows that the dispersion equations can allow multiple solutions with different wavenumbers $k$ at one frequency $\omega$, even when the waves have the same orientation and polarisation. (Note that, conversely, multiple solutions with different $\omega$ for a given wavevector $\kk$ are commonly present, as a usual consequence of frequency dispersion even in local media.)

The waves propagating with the higher wavenumber $k$ are denoted as \textit{additional}\index{wave!additional} \cite{agranovich2006spatial, agranovich2004linear, krowne2007book, agranovich1962crystal}. They were predicted by the works of Pekar et al. and also suggested by experimental data of dispersion near exciton levels, e.g. in cadmium sulphide \cite{pekar1975spatial}. 

The dispersion curves suggest the existence of additional waves around $\omega \sim 0.75 \omega_0$ in  Fig. \ref{fg_dcll_nl}a and $\omega \sim 1.75 \omega_0$ in  Fig. \ref{fg_dcll_nl}b. In both cases, one of the two solutions depicted has opposite signs of the wavevector and the group velocity ($|\vg| = \mathrm{d}\omega/\mathrm{d}k < 0$), predicting a negative refraction also in natural homogeneous media.

%}}}
\paragraph{Odd-power expansion terms and optical activity }   %{{{ 
Returning to the power expansion of $\epsLL(\omega, \kk)$ in terms of $k$ in Eq. (\ref{eq_epsmusd1}), we can identify the term constant in $k$ with the electric dipole moment $\chi_e(\omega)$, the term proportional to $k^2$ with the magnetic dipole moment $\chi_m(\omega)$ (or, also the electric quadrupole moment), and the recently discussed term proportional to $k^4$ with an electric octupole or magnetic quadrupole (\cite{agranovich2006spatial, agranovich2004linear, krowne2007book}).

The odd-power expansion terms were not discussed yet, although they have an important physical interpretation -- their nonzero values break the spatial inversion symmetry of the medium, and are thus related to optical activity \cite{bungay1993equivalency}. 
In media with nonzero odd-power terms, the corresponding eigenwaves are circularly polarized, and they propagate with different velocity depending on the direction. %, opposite signs of the odd expansion terms. 
Thus, the two plots in the right column of Fig. \ref{fg_dcllactivity} can also be viewed as the dispersion curves of the same medium, for the left and right circularly polarized waves.

\begin{figure}[t] \caption{Dispersion curves for two media with optical activity. The left and right columns of plots show the same information as in Figs. \ref{fg_dcll} and \ref{fg_dcll_nl}. The frequency dependence of the function $\gamma_e(\omega)$, which occurs in the term linearly proportional to $k$ in the expansion (\ref{eq_epsmusd1}), is plotted in the left column as a thin red line. } \label{fg_dcllactivity} \centering  
\textbf{(a)}\\\includegraphics[width=1\textwidth]{img/dispersion_landau_lifshitz/dispersion_ll_activep.pdf}
\textbf{(b)}\\\includegraphics[width=1\textwidth]{img/dispersion_landau_lifshitz/dispersion_ll_activen.pdf}
\end{figure}

%}}}

\subsection{Reflectance and transmittance at an interface of two local media}
\paragraph{Continuity requirements} %{{{ 
In the previous text, only the phase-related phenomena were discussed. The dependence of the dispersion curves and IFCs were computed as a result of the local or nonlocal response of the medium. Geometrical arguments were then used to infer the angle of refraction at the interface of two media, showing e.g. that a positive or negative refraction may occur. It was also shown that the beam may refract in a different direction than the wave vector, and that, in some cases, the notion of index of refraction can be used to simplify the problem to an application of the Snell law.

The conservation of the wave phase at the interface is, however, not the only constraint to the refraction and reflection problem. Assuming there is no surface current nor surface charge, and that both media are local, the components of the fields $\E$, $\HH$ \textit{parallel} with regard to the interface must be continuous and the \textit{perpendicular} component of the displacements $\D$, $\B$ is continuous, too. This rule can be derived from the Maxwell equations in Eqs. (\ref{eq_me1}--\ref{eq_me4}) by integrating them over infinitesimally thin loops or surfaces that are symmetrically placed on the interface \cite[pp. 26-29]{klingshirn2007semiconductor}. 

%}}}
\paragraph{Impedance} %{{{ 
The continuity requirement determines the \textit{amplitudes} of the waves reflected and transmitted at the interface. Corresponding complex \textit{Fresnel coefficients}\index{Fresnel coefficients} of \textit{reflectance}\index{reflectance!of interface} $r$ and \textit{transmittance}\index{transmittance!of interface} $t$ depend on the incidence angle and the polarisation of the wave \cite[p. 38]{born1999book}. They are derived in many textbooks with different levels of generalisation. 

In the case of isotropic media, the medium is sufficiently described by the ratio of permeability $\murl$ to permittivity $\epsrl$, whose square root is denoted as the \textit{impedance}\index{impedance of medium} of each medium 
\begin{equation} Z = \sqrt{\murl/\epsrl}. \label{eq_Z}\end{equation}
It shall be noted that most optics textbooks assume the media not to respond to the magnetic field, $\murl=1$, in which case $Z = 1/N = \sqrt{1/\epsrl}$ and the Fresnel coefficients then are be expressed as functions of $N$. This assumption is however not applicable to most structures studied in this thesis.

%}}}
\paragraph{Perpendicular incidence of two local media} %{{{ 
In the simplest case where the wave perpendicularly impinges a single interface of two isotropic local media, described by their respective impedances $Z_1, Z_2$, the reflectance $r$ and transmittance $t$ of the interface are
\begin{equation} r = \frac{Z_2 - Z_1}{Z_2+Z_1}, \quad t = \frac{2 Z_2}{Z_2 + Z_1}, \label{eq_reflection}\end{equation}

%}}}
\paragraph{Refraction on an interface with a nonlocal medium}   %{{{
The task to compute the amplitudes of reflectance and transmittance for nonlocal media is substantially more complicated than for the local media. The principal reason is related to the fact that the definition of the medium response (\ref{eq_chi_convol_nonloc}) contains a convolution over an infinite space, which naturally does not account for the interface. For the problem to be better formulated, the integral would have to be modified, either with the nonlocal response of medium to extend also behind the interface, or with the medium response sharply truncated at the interface. For a weak spatial dispersion, different approaches to the problem are discussed in \cite{golubkov1995boundary}; a similar problem for one example of periodic structures is numerically analyzed in \cite{lapine2012surface}.

Strong enough spatial dispersion can also allow the plane wave to refract into a superposition of two independent plane waves, which differ by their component of the wavevector perpendicular to the interface. The respective amplitudes of these components depend on the \textit{additional boundary conditions}\index{additional boundary conditions}\index{additional boundary conditions} \cite{agranovich2006spatial}.

%}}}

\subsection{Phase, group, energy and signal velocities}
\label{chap_vfvg}
\paragraph{Signs of the phase and group velocities}%{{{ 
The negative refraction in media shown in Figs. \ref{fg_ifcnr}b,c and \ref{fg_ifcnrk}b,c was a result of the requirement for the group velocity $\vg$ to conserve its component perpendicular to the interface, and for the wave vector $\kk$ to conserve its projection onto the plane of interface. 

It was assumed that the group velocity coincides with the energy or information propagation, and therefore that it should always be "positive", that is, it should propagate towards the interface in the first medium and outwards from it in the second one. In the following paragraphs, this assumption will be challenged for special cases of absorbing media, but it will be shown that in media with low losses it is correct and the term of \textit{negative phase}\index{velocity!negative phase velocity}  is used when it is opposite to the group velocity.
Note that since the phase velocity does not carry energy, it can be also higher than the speed of light in vacuum.

%}}}
\paragraph{Signs of the group and energy velocities}   %{{{
%  \mdf{
%  % TODO 0401 how can one tell apart negative refraction due to eps-mu and due to higher, purely nonlocal terms?
%  *doubly resonant local media -> negative phase velocity, positive group velocity, above resonant frequencies - theoretically arbitrarily small losses\\
%  *singly resonant local media -> positive phase velocity, negative group velocity, this happens near resonance so there are always high losses\\
%  *			  nonlocal media -> weird things happen,
%  }
In the discussion of wave refraction on an interface, it was assumed that the envelope (i.e. the modulation, carrying information) of the wave approaches the interface in the first medium, then refracts and propagates away \textit{from} the interface in the second medium. The envelope propagates with the group velocity $\vg$ as defined by Eq. (\ref{eq_vg}). There is another quantity, the \textit{Poynting vector}\index{Poynting vector} $\mathbf{S} := \E\times\HH$, describing the direction and density of the \textit{power} carried by the optical wave \cite[p. 16]{klingshirn2007semiconductor}. 

The group velocity usually points in the same direction as the Poynting vector, but this has not necessarily to be always true. A typical counterexample can be found near resonances in lossy (local) media. Such a behaviour can be traced back to Fig. \ref{fg_oscillator_spectrum}, where the permittivity drops from high values to negative ones. Since the medium is defined lossy in the figure, the curve of the wavenumber $k(\omega)$ is continuous and smooth, and as a result the magnitude of the group velocity $v_g = \mathrm{d}\omega / \mathrm{d}k$ is negative for a range of frequencies around the resonance for $\omega \sim 2\pi$, %% TODO "negative" means "opposite to kk"? But we are speaking of scalars...  which is also known as \textit{anomalous dispersion}\index{dispersion!anomalous}. 

It follows from this that the group velocity $\vg$ can also point in a direction opposite to the Poynting vector that represents propagation of the light beam energy, $\mathbf{S}$. 
% This fact cannot be easily deduced from the dispersion curves or IFCs. %% TODO  Or can it?
This can be observed in narrow parts of the spectrum only, around the resonant frequencies where the media have high losses, and also a strong group velocity dispersion.
Assuming passivity and absence of sources in the second medium, it is obviously the Poynting vector that has to conserve its perpendicular component and that the group velocity has inevitably to point towards the interface also in the second medium, which seems contradictory to causality.
This result 
%of the light envelope propagating with the group velocity \textit{from}\index{from} the second medium \textit{towards}\index{towards} the interface may seem unphysical, 
may however be explained by a strong deformation of the envelope shape on a short distance. %, which ensures the information carried by the wave modulation to propagate causally, out from the source, and always with speed lower or equal to the speed of light in vacuum $c$.

A negative group velocity has been experimentally observed in a thin sample. It manifested itself as a negative shift of the light envelope center-of-mass \cite{dolling2006simultaneous},  compared to the absence of the medium. Note that the negative group velocity is independent of the sign of the phase velocity, which may be both positive and negative \cite{mikki2009electromagnetic}. 

%}}}
\paragraph{Signal velocity}%{{{
The notion of \textit{signal} or \textit{information propagation velocity}\index{velocity!of information propagation} is sometimes % TODO cite
identified with the group velocity, $\vg$. However, this may be misleading, as the differential definition of the group velocity in Eq. (\ref{eq_vg}) enables one to define it for slow-enough modulation only, as long as  the span of frequencies is narrow and the second derivative, corresponding to the group velocity dispersion, can be neglected.

Assuming the information is carried by a wave modulation that is limited in time (e. g. presence or absence of an optical pulse), one leaves the comfortable approximation of narrow spectrum: from the convolution-multiplication theorem already used in Eq. (\ref{eq_kkresult}), it can be shown that any information carrying function with a compact support has to span over an infinite spectrum. The simplified version of the proof \cite{hill2013uncertainty} can be based on the fact that under multiplication by a well chosen compact-supported binary function, e.g., 
$$ f(t) \rightarrow f(t) \cdot \mathrm{sign}^2[f(t)], $$
the original function obviously does not change, yet the Fourier transform of this binary function is nonzero over almost all frequencies. 

A temporally limited optical pulse will be always more or less distorted upon propagation in a dispersive medium where the dispersion curves are not straight. Then the information velocity becomes a problematic term. If the information could be detected exactly in the first moment when the fields deviate from strict zero, one would come to the surprising conclusion that, even in dispersive media, the information would propagate with the speed of light. In the opposite example of a highly noisy transmission, the pulse may be reliably detected only after most of its energy already arrived, resulting in information propagation even slower than the group velocity. 

To conclude, the author is convinced that the notion of \textit{information velocity} is too vague unless its exact definition is provided first, and should not be directly associated with the group velocity. 

%}}}
%% TODO \paragraph{Speed of light and causality in nonlocal media}%{{{
%% \mdf{TODO REF Kramers-Kronig relations in nonlocal medium
%% \cite{skettrup1970kramers}
%% \cite{kirzhnitz1976}
%% \cite{melrose1977generalised}
%% \cite{sun1989kramers}
%% \cite{rozanov2003}
%% \cite{bruleanalysis}
%% \cite{makarov2013kramers}
%% }%}}}
%}}} /section
\section{Electromagnetic waves in periodic structures}%{{{
\subsection{Periodic structures and the Bloch's theorem}
\paragraph{Inhomogeneity}%{{{
The previous sections discussed infinite media, with the only deviation from homogenity at an interface of two media, where refraction of the waves is observed, and it was shown that the resulting orientation of the wave vector and of the group velocity could be easily deduced on a geometrical basis. The amplitudes of the reflected and refracted waves can be also easily computed for local media, whereas their computation for nonlocal media is much less straightforward. 

There is a broader class of shapes for which analytical or semianalytical methods have been developed to compute their interactions with electromagnetic waves, such as the Mie theory for scattering on dielectric or metallic spheres and cylinders, propagation through arbitrary stacks of parallel layers, diffraction on narrow apertures, resonances in orthogonal, cylindrical or spherical cavities, or wave guiding in high-symmetry waveguides or optical fibres. 

The interactions of electromagnetic waves with most of the possible shapes are too complex to be expressed analytically, and they can only be accessed by numerical methods, some of which are described in Chapter \ref{chapter_numerical}. However, when these elementary shapes are arranged into an infinite array, any resulting periodic structure behaves in a way typical for the periodicity and its most important traits can again be partially understood on an analytical basis.
This chapter focuses on these general properties shared by periodic structures, leaving the particular numerical simulations to the Results section.

%}}}
\paragraph{Periodicity}%{{{
Under the notion of \textit{periodicity}\index{periodicity} we understand the existence of discrete translational symmetries. In three dimensions, we can write
\begin{equation} \epsrl(\omega, \rr) = \epsrl(\omega, \rr + m_1 \ava + m_2 \avb + m_3 \avc) \label{eq_trsym}\end{equation}
where $m_1, m_2, m_3$ are integers and $\ava, \avb, \avc$ are three linearly independent vectors. The local permittivity, given by Eq. (\ref{eq_lorentz_eps}), was intuitively generalized to a function of the position $\rr$. The constituent media are described by \textit{local} quantities only, to avoid possible problems at the boundary when computing the spatial convolution  in nonlocal media, Eq. (\ref{eq_eps_nonloc}).  The same periodicity is also imposed on the permeability $\murl(\omega,\rr)$. 

The points generated by all combinations of possible translations of the unit cell center by $m_1 \ava + m_2 \avb + m_3 \avc$ form a periodic lattice.
The volume defined as a \textit{set of points closest to} one given point of the lattice will be denoted as a \textit{unit cell}\index{unit cell} (similar to the \textit{Wigner-Seitz cell}\index{Wigner-Seitz cell} in solid-state physics). Obviously, the permittivity or permeability in periodic structures needs only to be specified within one unit cell.

The choice of lattice vectors $\avabc$ limits the set of the rotation or mirror symmetries of the structure. Based on the allowed symmetries, all lattices in three dimensions can be classified into six \textit{crystallographic systems}\index{crystallographic systems}, namely \textit{cubic}, \textit{tetragonal}, \textit{ortorhombic}, \textit{monoclinic}, \textit{hexagonal-trigonal}, all periodic structures can be classified into 230 \textit{crystallographic space groups}. Numerous crystal optics textbooks (e.g. \cite[p. 678]{born1999book}) give more rigorous definitions. 

Unless stated otherwise, we will assume that the cubic lattice is used, which allows the highest possible symmetry. In the cubic lattice, the lattice vectors $\avabc$ are of the same magnitude $a$ and mutually orthogonal. % Note, however, that even the structure with a periodic modulation of permittivity is no more isotropic, i.e. posessing continuous rotational symmetry.

%}}}
\paragraph{The Bloch's theorem}%{{{
The \textit{Bloch's} (or \textit{Bloch's-Floquet's}) theorem states that while the harmonic wave is no more a solution for the Maxwell equations in a periodic structure that conforms to Eq. (\ref{eq_trsym}), a solution can always be found as a \textit{Bloch's wave}\index{Bloch's!wave} -- a product of a harmonic function and another periodic one:
\begin{equation} 
\E(t, \rr) = \mathbf{u_e}(\rr)\,\mathrm{e}^{\ii\omega t - \ii\KK\cdot\rr}, \text{ where } \mathbf{u_e}(\rr) = \mathbf{u_e}(\rr + m_1 \ava + m_2 \avc + m_3 \avc),
\label{eq_bloch}\end{equation} 
\begin{equation}
\HH(t, \rr) = \mathbf{u_m}(\rr)\,\mathrm{e}^{\ii\omega t - \ii\KK\cdot\rr}, \text{ where } \mathbf{u_m}(\rr) = \mathbf{u_m}(\rr + m_1 \ava + m_2 \avb + m_3 \avc).
\label{eq_blochh}\end{equation} 
As a rule, in all linear systems, any sum of Bloch's waves is also a proper solution, but for simplicity we will focus on one Bloch's wave at a time.

The functions $\mathbf{u_e}(\rr)$ and $\mathbf{u_m}(\rr)$ have the same periodicity as the structure, and will be denoted as the \textit{mode functions}\index{mode functions}. They are, in general, complex vector functions, so they not only alter the direction and magnitude of the electric and magnetic fields, but can also introduce a \textit{phase modulation} of the wave in each unit cell. 

The remaining term, $e^{\ii\omega t - \ii\KK\cdot\rr}$, is analogous to that of a plane wave
\begin{equation} \E(t, \rr) := \E_0\, e^{\ii\omega t - \ii\kk\cdot\rr}, \tag{\ref{eq_pw} \again} \end{equation}
except for the capital $\KK$ being used to distinguish the wave vector of the Bloch's wave envelope from the wave vector $\kk$ in homogeneous media. 

Note that the Bloch's theorem does not determine the shape of $\mathbf{u_e}(\rr)$, $\mathbf{u_m}(\rr)$, nor the direction and magnitude of $\KK$, it only states a solution in the form of Eqs. (\ref{eq_bloch}, \ref{eq_blochh}) can be found.
% TODO proof /home/filip/PhD/Sources_MM_theory/Bloch_Theorem_Proof.pdf
% TODO example illustration of Bloch's wave -> plotted in Python

%}}}
\paragraph{Proof of the Bloch's theorem in one dimension}%{{{
This theorem is essential for understanding the electromagnetic behaviour of periodic structures, and it deserves a proof. Originally, it was developed for the electron wave function $\psi$ in crystals on the basis of quantum mechanics. The outline of such a proof in one dimension is based on the following:
\label{blochproof}
\begin{enumerate}
 \item{We assume $\psi$ is an eigenfunction of the Hamiltonian: $\exists h\in \mathbb{C}: \hat H\psi = h\psi$.} 
 \item{We also assume that the Hamiltonian operator $\hat H$ commutes with the operator of discrete translation $\hat T$ by the inter-atomic distance: $\forall \psi: \hat H\hat T\psi = \hat T\hat H\psi$. } 
 \item{Then $(\hat T\psi)$ is an eigenfunction of $H$, because obviously $\hat H(\hat T\psi) \stackrel{1.}{=} \hat T\hat H\psi \stackrel{2.}{=} \hat Th\psi = h(\hat T\psi)$. }
 \item{From two eigenfunction relations, $\hat H\psi\stackrel{2.}{=} h\psi$ and $\hat H(\hat T\psi) \stackrel{3.}{=} h(\hat T\psi)$, it also follows that $\psi$ and $\hat T\psi$ must either represent the very same physical eigenstate of $\hat H$ that is uniquely related to its eigenvalue $h$, or otherwise that there must exist two or more different \textit{degenerate} states $\psi_1\neq \psi_2$ with the same eigenvalue $h$.\\ 
 The latter case of \textit{degeneracy} is proven in many textbooks not to change the conclusion that $\psi$ must be also an eigenfunction of the translation operator, that is,
$$\exists K\in \mathbb{C}: \hat T\psi = e^{-\ii Ka}\psi.$$ 
}
\end{enumerate}
 The physical consequence is that when the Hamiltonian operator is invariant to a discrete translation, its eigenfunctions are also \textit{almost unchanged} upon this particular discrete translation, since they may differ only by a phase shift of $-Ka$. Setting $a$ to be the unit cell size, $K$ becomes the wavenumber of the Bloch's wave envelope.
%
%Two commuting operators share the same eigenstates, which means that the electron wave function at a given energy (i.e. Hamiltonian eigenfunction) must also have a discrete translational symmetry provided its phase can change between atoms (i.e. it is an eigenfunction of the non-hermitian translation operator).

%}}}
\paragraph{The Bloch's theorem proof in the electromagnetic formulation}%{{{
\index{Bloch's!theorem, proof}
The steps can be reformulated replacing the abstract Hamiltonian with an operator derived from the Maxwell equations. In a periodic structure one can no longer assume the solution in the form of a plane wave (\ref{eq_pw}), but as long as the structure is time-invariant and linear, the monochromatic electric and magnetic fields can still be decomposed into a product of some complex function of space position, and a harmonic function of time:
\begin{equation} 
\begin{array}{cc}
\E(t, \rr) = \E(\rr) e^{\ii \omega t}, \\
\HH(t, \rr) = \HH(\rr) e^{\ii \omega t}. 
\end{array}
\label{eq_harmonic}\end{equation}
One thus only needs to prove the Bloch's theorem for the time-invariant parts of the fields, which will play the same role as the wavefunction $\psi$ in the quantum-mechanical formulation.

\begin{enumerate}
\item{
Assuming $\E(\rr)$ to be a  valid solution of Maxwell equations in the periodic structure at the angular frequency $\omega$,
it can be derived from the Maxwell Eqs. (\ref{eq_me3}, \ref{eq_me4}), taking into account the local medium response to the fields as defined by Eq. (\ref{eq_epstensor}) to the harmonic wave from Eq. (\ref{eq_harmonic}), that
\begin{equation} 
{\epsrl}^{-1}(\omega,\rr) \nabla\times \left[{\murl}^{-1 }(\omega,\rr) \nabla\times \E(\rr) \right] = \frac{\omega^2}{c^2}\E(\rr),   \label{eq_eigen_e}
\end{equation}
In analogy with the quantum-mechanical proof, the left hand side of Eq. (\ref{eq_eigen_e}) can be associated with the Hamiltonian $\hat H\psi$, and the right hand side with its eigenvalue $h\psi$ \cite{johnson2003introduction}. 
} 
\item{
The translation operator acts by substitution of the position vector $\rr$ in the argument. For example, the translated electric field is
\begin{equation} \hat T\psi \rightarrow \E(\rr+\mathbf{a_1}). \label{eq_blocht}\end{equation}
The commutation relation directly follows from the periodicity in Eq. (\ref{eq_trsym}), thus the right terms in Eqs. (\ref{eq_blochht}) and (\ref{eq_blochth}) are identical by definition:
\begin{equation} \hat H \hat T \psi \quad  \rightarrow \quad  {\epsrl}^{-1}(\omega,\rr    ) \nabla\times \left[{\murl}^{-1 }(\omega,\rr    ) \nabla\times \E(\rr+\ava)\right]  \label{eq_blochht}\end{equation}
\begin{equation} \hat T \hat H \psi \quad  \rightarrow \quad  {\epsrl}^{-1}(\omega,\rr+\ava) \nabla\times \left[{\murl}^{-1 }(\omega,\rr+\ava) \nabla\times \E(\rr+\ava)\right] .  \label{eq_blochth}\end{equation}
In three dimensions, this argument is valid for three different translation operators that correspond to the displacements by the different lattice vectors, $\ava$, $\avb$ and $\avc$.
} 
\item{
Then the wave translated by any of the lattice vectors is also a solution of the Maxwell equations:
\begin{equation}  {\epsrl}^{-1}(\omega,\rr) \nabla\times \left[{\murl}^{-1 }(\omega,\rr) \nabla\times   \E(\rr+\avabc)\right] =  \frac{\omega^2}{c^2}\E(\rr+\avabc).  \label{eq_blochme}\end{equation}
}
\item{In analogy with the fourth step in the Bloch's theorem proof, there exists at least one constant $K_1$ for which
\begin{equation}  \E(\rr+\ava) = e^{-\ii K_1 a_1} \E(\rr),  \label{eq_blochtt}\end{equation}
	where $K_1 a_1$ represents the phase shift between the adjacent cells along the direction of the lattice vector $\ava$. In a similar manner, constants $K_2$ and $K_3$ can be associated with the translations by the $\avb$ and $\avc$ vectors, respectively. 
}
 \end{enumerate}

%}}}
\paragraph{Bloch's theorem in three dimensions} %{{{ %% TODO unify 1-D 2-D and 3-D terminology
The three constants $K_{1,2,3}$ then define the \textit{Bloch's wave vector}\index{wave vector!of Bloch's wave}\index{Bloch's!wave vector} $\KK$ in three dimensions: 
\begin{equation} \KK := 
\frac{2\pi K_1(\avb\times\avc)}{\ava\cdot(\avb\times\avc)} +  
\frac{2\pi K_2(\avc\times\ava)}{\avb\cdot(\avc\times\ava)} +  
\frac{2\pi K_3(\ava\times\avb)}{\avc\cdot(\ava\times\avb)}.
\label{eq_reciprocalK}\end{equation}
Each of three addends in Eq. (\ref{eq_reciprocalK}) relates to one of the lattice vectors and is orthogonal to the two remaining lattice vectors. %In contrast, inner product of $\KK$ with any of the three elementary lattice vectors is equal to $2\pi$. % elementary translation by 
If the lattice vectors $\avabc$ form an orthogonal triplet, Eq. (\ref{eq_reciprocalK}) simplifies to
% $K_{1,2,3}$ directly define the projection of the $\KK$:
\begin{equation} \KK := 
	\frac{2\pi K_1 \ava}{a_1^{2}} +  
	\frac{2\pi K_2 \avb}{a_2^{2}} +  
	\frac{2\pi K_3 \avc}{a_3^{2}}.
\label{eq_reciprocalKorto}\end{equation}

Repeating the procedure for the magnetic field, a similar operator can be obtained, where the multiplication by ${\epsrl}^{-1}(\omega,\rr)$ and ${\murl}^{-1 }(\omega,\rr)$ occurs in the opposite order. Thus the same arguments, with the identical wave vector $\KK$, have to be valid also for the magnetic field $\HH(\rr)$, where the Hamiltonian is associated with an operator very similar to that of the electric field in Eq. (\ref{eq_eigen_e}):
\begin{equation}
{\murl}^{-1 }(\omega,\rr) \nabla\times \left[{\epsrl}^{-1}(\omega,\rr) \nabla\times \HH(\rr)\right] = \frac{\omega^2}{c^2}\HH(\rr).   \label{eq_eigen_h}
\end{equation}

%}}}
\paragraph{Virtual periodicity and ambiguity of the mode function}%{{{
Homogeneous media, such as vacuum, naturally fulfill the definition of periodicity in Eq. (\ref{eq_trsym}).
From the principle of correspondence, the Bloch's theorem must predict
the already known solution of a harmonic plane wave in vacuum
\begin{equation} \E(t, \rr) := \E_0\, e^{\ii\omega t - \ii\kk\cdot\rr}. \tag{\ref{eq_pw} \again} \end{equation}
The expected solution of the Bloch's wave in vacuum can be directly found as a formal modification of the dispersion relation,
\begin{equation}  
\E(t, \rr) = \mathbf{u_e}(\rr)\,\mathrm{e}^{\ii\omega t - \ii\KK\cdot\rr}, \text{ where } \mathbf{u_e}(\rr) := \E_0  \text{ and } \KK := \kk.
\label{eq_blochvac}
\end{equation}
However, this is not the only possible representation of a plane wave as the Bloch's wave, as the Eqs. (\ref{eq_bloch}, \ref{eq_blochh}) are mathematically ambiguous. 
For illustration, using the fact that any lattice of periodic unit cells may be imagined in vacuum, one can choose a cubic lattice with an arbitrary unit cell size $a$, % TODO represent by lattice vectors G?
for which the plane wave in vacuum can be simultaneously represented by any of infinitely many other combinations of
\begin{equation} 
\left.  \begin{array}{cc}
	\mathbf{u_{e}}(\rr) := \E_0 e^{2\pi\ii \left(\frac{m_1}{\ava\cdot\rr} + \frac{m_2}{\avb\cdot\rr} + \frac{m_3}{\avc\cdot\rr} \right)} \vspace{2mm} \\ %% /a^{2} TODO use mathbf here? or is this whole a rubbish?
%\KK := \kk - \frac{2\pi m_1}{a_1} - \frac{2\pi m_2}{a_2}- \frac{2\pi m_3}{a_3},
\end{array}\quad \right\} \quad \forall m_{1,2,3} \in \mathbb{Z}
\label{eq_blochvac2}
\end{equation}
%\begin{equation} \KK := 
	%\frac{2\pi K_1 \mathbf{a_1}}{a_1^{2}} +  
	%\frac{2\pi K_2 \mathbf{a_2}}{a_2^{2}} +  
	%\frac{2\pi K_3 \mathbf{a_3}}{a_3^{2}}.
%\label{eq_reciprocalKorto}\end{equation}
which all 
%maintain the requirement of the mode function $\mathbf{u_{e}}(\rr)$ being periodic and still 
give the exactly same resulting plane wave. 

\begin{figure}[ht] \caption{Folded and unfolded dispersion curves for free space and a periodic structure viewed along one of its axes. The phase difference across a unit cell can be expressed either by the wavenumber $K$, as it is in the unfolded plots \textbf{(a)} and \textbf{(b)}, or it can be predominantly absorbed into the periodic mode function $\mathbf{u}(\rr)$  as in the folded plots \textbf{(c)} and \textbf{(d)} with inaccessible $K$-values greyed out. \\
The lattice periodicity allows one to draw the dispersion curves of periodic structures using natural scale-invariant units, with the Bloch's wave number $K$ divided by the spatial frequency of the lattice $2\pi/a$, and the angular frequency $\omega$ multiplied by the time needed for the light to traverse the unit cell.
} \label{fg_phcfolding} \centering 
	\includegraphics[width=\textwidth]{img/PhC_folding_illustration.pdf} 
\end{figure}
\clearpage

%}}}

\subsection{Dispersion in periodic structures}
\paragraph{Folding of the dispersion curves} %{{{
All possible dispersion curves for a \textit{Bloch's wave in vacuum}\index{Bloch's!wave in vacuum} with the $\KK$ vector oriented along one lattice axis in vacuum are plotted in Fig. \ref{fg_phcfolding}a. The dispersion curve for the original forward wave from Eq. (\ref{eq_blochvac}) is plotted in blue; another solution exists for a wave propagating in the opposite direction which is plotted in red. All other solutions generated by Eq. (\ref{eq_blochvac2}) are plotted in gray, for both forward and backward waves. 
Note that the actual physical shape of the fields in space,  $\E(\rr)$ and $\HH(\rr)$, cannot be deduced from the dispersion curve. To fully determine the fields, the dispersion curves would have to be complemented by the mode functions $\mathbf{u_e}(\rr)$ and $\mathbf{u_m}(\rr)$, which depend on the frequency. 

\label{par_disp_curv_per}
The inherent ambiguity of the description on the basis of the Bloch's wave can be used to save space in the plot by showing the dispersion curves only for $K\in\langle-\pi/a, \pi/a\rangle$, as shown in Fig. \ref{fg_phcfolding}c. Plotting any other interval of equal width would be equivalent \cite[p. 177]{brillouin2003wave}. 

The space of the plot can be further halved, since all forward and backward propagating solutions are symmetrical with respect to the vertical $K=0$ axis, except for structures with optical activity (cf. Fig. \ref{fg_dcllactivity}) or those breaking the time-reversal symmetry (cf. Ref. \cite{vanwolleghem2009unidirectional}), which are not discussed in this thesis. 

In many papers on periodic structures, all dispersion curves are plotted as \textit{folded}\index{dispersion!curves,folding of} into the $K\in\langle -\pi/a, \pi/a\rangle$ region only \cite{obrien2002photonic, yannopapas2005negative, chakrabarti2012magnetic}, while some employ the symmetry to fold  all curves further into the $K\in\langle0, \pi/a\rangle$ range. % (e.g. \cite{}). 
In some other references, such as \cite{mortensen2010unambiguous} or \cite{yeh1977electromagnetic}, unfolded dispersion curves are used.

The mathematical interpretation of folding the dispersion curves in vacuum can be found from Eq. (\ref{eq_blochvac2}): 
\begin{enumerate}
\item{
For  $\omega\in\langle 0, \pi c /a\rangle$, the mode function is constant in space, $\mathbf{u_e}(\rr) := E_0$, and $K := \omega/c$. 
} 
\item{
At the first point of folding for $K=\pi/a$, the mode function changes to $\mathbf{u_e}(\rr) := E_0 e^{2\pi\ii r/a}$ and the formal solution of the \textit{backward} wave is used with $K := 2\pi r/a - \omega r/c$.  The phase increase across one unit cell remains positive, since the phase decrease of the backward-wave envelope is less than the phase contribution of the mode function.
} 
\item{
When the dispersion curve touches the $K=0$ axis at the frequency $\omega = 2\pi c/a$, the mode function is not changed, but the dispersion curve continues again along with the \textit{forward} wave $K := - 2\pi r/a + \omega r/c$. 
} 
 \end{enumerate}
In vacuum, this process repeats with higher orders of the mode function: $\mathbf{u_e}(\rr) := E_0 e^{4\pi\ii r/a}$, $\mathbf{u_e}(\rr) := E_0 e^{6\pi\ii r/a}$, etc.

We will show in the Results section on a particular example that some numerical approaches provide rules for the selection of one particular dispersion curve in the unfolded plot. The Results section also contains a physical interpretation of this finding.

%When the wavelength is less than double cell spacing, or equivalently when $2\pi c /K < 2 a$, the unit cells can scatter the forward propagating wave into the backward direction and vice versa. In an infinite periodic structure the energy transfer must be conserved, so the forward and backward waves must be of the same amplitude. (The equilibrium of forward and backward wave amplitudes may, however, differ in some cases even in infinite periodic structure. One example is when the constituent materials lose time-reversal symmetry due to static magnetic field.) We assume this also in this example with virtual periodicity. The \textit{Bloch's wave}\index{Bloch's wave} is a product of two functions, and it gives one degree of freedom to select which part of (\ref{eq_bloch}) expresses the phase shift across one unit cell. Two different conventions are used in the literature.
% \begin{enumerate}
%  \item{In the \textit{unfolded}\index{unfolded} plot, the phase shift is expressed entirely by the envelope $\mathrm{e}^{\mathrm{i}Kz}$. One can always find a line coming from the front side of a unit cell to its back side, along which the mode function does not change its phase significantly $\mathbf{u(\mathbf{r})}$. In other words, the volume of the unit cell is not divided by any \textit{nodal plane}\index{nodal plane} in the real part of electric field  $\mathbf{u(\mathbf{r})}$. (However, closed nodal planes within the unit cell may appear due to individual resonances.)
% 
% In vacuum, the mode function takes the simplest form, being constant $\mathbf{u(\rr)} = 1$ at all frequencies. The dispersion relation of vacuum forms a direct line as plotted in Fig. \ref{fg_phcfolding}a.} 
%  \item{The \textit{folded}\index{folded} plot requires $K/\left(\frac{2\pi}{a}\right) \in (-\frac{1}{2}, \frac{1}{2})$ and the remaining phase is mostly expressed by the mode function. 
% Further savings of space in the plot can be made using the symmetry of the forward and backward waves, so only the positive half for $K/\left(\frac{2\pi}{a}\right) \in (0, \frac{1}{2})$ has to be plotted to describe the mode structure.} 
% \end{enumerate}
% Although maybe less instructive, usually the \textit{folded}\index{folded} bands are plotted in the literature and they are used also in Fig. \ref{fg_1dbd}, \ref{fg_rodh} and \ref{fg_erod_radius11} below. The same approach can be used for periodic structures, where the bands are separated by band gaps, cf. Fig. \ref{fg_phcfolding}b and \ref{fg_phcfolding}d.

%}}}
\paragraph{Fourier expansion of the mode function}%{{{ 
If the periodic structure is not homogeneous, that is, $\epsrl(\omega,\rr)$ or $\murl(\omega,\rr)$ depends on $\rr$, the mode functions $\mathbf{u_{em}}(\rr)$ acquires more complex, anharmonic shape than assumed in Eq. (\ref{eq_blochvac2}). They however still show the periodicity of the lattice, as dictated by the Bloch's theorem in Eq. (\ref{eq_bloch}), and their exact shapes can thus be decomposed into Fourier series,
characteristic for the structure and the frequency of operation, defining nonzero amplitudes for all branches. The grey lines plotted in Fig. \ref{fg_phcfolding} then are no more hypothetical solutions to be chosen from, but they all have real amplitudes, and are present simultaneously. 
Note this does not contradict the ambiguity of choice of the Bloch's wavevector $\KK$, which stems from the  mathematical representation of the Bloch's wave.

%}}}
\paragraph{Brillouin zones in the reciprocal space} %{{{ 

The periodically repeating interval of $K\in\langle-\pi/a, \pi/a\rangle$ can be viewed as the unit cell in the space of wave vectors, i.e. in the \textit{reciprocal space}\index{reciprocal space}. This unit cell is then denoted as the first \textit{Brillouin zone}\index{Brillouin zone}. 

In two or three dimensions, any $m$-th Brillouin zone is defined as the set of all points for which the center point $K=0$ is the $m$-th closest point from the regular lattice of all points where $\KK = (2\pi m_1/a, 2\pi m_2/a, \ldots)$ for $m_1, n_2, \ldots \in \mathbb{Z}$.
All Brillouin zones have equal measure (i.e. length, area or volume), and moreover, they can be transformed into the first one by simple translation. %% TODO specify what this re-arrangement means

In the 1-D case, the definition of higher Brillouin zones trivially consists in dividing the axis symetrically into equal intervals, such as $$K\in\langle-2\pi/a, -\pi/a\rangle \cup \langle \pi/a, 2\pi/a\rangle$$ for the second one etc. 

In two or three dimensions the shapes of higher Brillouin zones become much more complex, and are unlike to each other \cite[pp. 134--135]{klingshirn2007semiconductor}: in two dimensions, all Brillouin zones are composed of polygons that are connected in their vertices to enclose the previous Brillouin zone \cite[p. 126]{brillouin2003wave}. Similarly, in three dimensions, they are composed of polyhedra connected by their edges, and their shape roughly resembles sphere shells. 

The above illustration of how the plot space can be saved by folding the dispersion curves into the first Brillouin zone also applies to two- or three-dimensional structures.
The exact shapes of the Brillouin zones are determined by the crystal family of the real-space lattice and its set of lattice vectors $\avabc$. 
As noted above, the space occupied by the plot can be further reduced thanks to the symmetry of the dispersion curves. This applies naturally to higher dimensions as well; the minimum part of the Brillouin zone that contains all necessary information is denoted as the \textit{irreducible Brillouin zone}\index{Brillouin zone!irreducible} and is determined by the crystallographic point group of the lattice. 

A more detailed listing and comparison of Brillouin zone shapes is beyond the scope of this thesis, but can be found in numerous solid-state textbooks, e.g. \cite[pp. 96-99]{klingshirn2007semiconductor}. 

%}}}
\paragraph{High-symmetry points of the Brillouin zones} %{{{ 
\begin{figure}[t] \caption{High-symmetry points in the reciprocal $\KK$-space and the Brillouin zones for \textbf{(a)} a two-dimensional square lattice and \textbf{(b)} a three-dimensional cubic lattice. For these lattices, the first Brillouin zones have square and cubic shapes, respectively. In each figure, one of the irreducible Brillouin zones is highlighted.  } \label{fg_phcbrillouin} \centering 
	\includegraphics[width=.9\textwidth]{img/PhC_high_symmetry_points.pdf}
\end{figure}
The fundamental difference between plotting one dispersion curve and the isofrequency contours (IFCs) is that IFCs provide information about the dispersion for a two-dimensional subspace of the wave vector $\KK$, whereas the dispersion curve is limited to a 1-D scan along some line in the $\KK$-space (cf. Fig. \ref{fg_ifc_dc}). The limitation is even more significant in the case of three dimensions, where it is virtually impossible to visualize the frequency $\omega(\KK)$ as a function of three components of $\KK$.

A hybrid approach can be taken, however, that maintains most of the physical information about the dispersion in the two- or three-dimensional structures. It consists in plotting the  dispersion curves around the boundary or edges of one selected irreducible Brillouin zone. It is expected that no important information is lost by leaving out its inner surface.

The notation of high-symmetry points in the $\KK$-space is standardized. For the square lattice in two dimensions, four high-symmetry points are named as follows: "$\mathbf{\Gamma}$" corresponds to the center, $K=0$. Provided that the first lattice vector $\ava$ is parallel to the $x$-axis of the coordinate system, the point where $\KK = (\pi/a, 0)$ is denoted as  "$\mathbf{X}$" and the point $\KK = (0, \pi/a)$ as "$\mathbf{Y}$". Finally, the diagonal point of $\KK = (\pi/a, \pi/a)$ is known as "$\mathbf{M}$".
In a two-dimensional case, the dispersion curves are typically plotted along the triangle encircling the irreducible Brillouin zone: $\mathbf{\Gamma}-\mathbf{X}-\mathbf{M}-\mathbf{\Gamma}$.

In three-dimensional cubic lattice, these points maintain their meaning in the plane perpendicular to the $z$-axis \cite[p. 99]{klingshirn2007semiconductor}. Additionally, "$\mathbf{R}$" denotes the spatial diagonal of $\KK = (\pi/a, \pi/a, \pi/a)$. Other lattices introduce more complex sets of the high-symmetry points and the paths along which the dispersion curves are plotted do not seem to be standardized.

%}}}
%% \paragraph{Dispersion curves for the reduced Brillouin zone} %{{{ 
%% %}}}
\paragraph{Tight-binding model for the dispersion curves} %{{{
The notation of the high-symmetry points in the Brillouin zone can outline a different way \cite{shi2007} of understanding the origin of the dispersion, which can be drawn in analogy with the \textit{tight-binding},   %% TODO cite?
or also \textit{hopping},  %% TODO cite?
model in solid-state physics. 

\label{tight-binding}
Analyzing one unit cell isolated in free space, one obtains its natural resonant frequency.  When the cell is surrounded by other cells in the periodic lattice, their mutual coupling alters the resonant frequency \cite[p. 75]{klingshirn2007semiconductor}; the actual sign of the coupling effect does not seem to be determined by any simple general rule, however.	

For a steady oscillatory solution to be obtained, the electromagnetic fields in the surrounding cells must have the same modulus of oscillation amplitude, but they do not need to share the phase. This is equivalent to a single Bloch's wave propagating in the structure. 

Identical phase in all cells corresponds to a zero Bloch's wave vector, $\KK = 0$, and to the point $\mathbf{\Gamma}$ in the center of the Brillouin zone. Under the basic assumption that the fields $\E, \HH$ are continuous and smooth vector functions, the requirement of an identical phase of adjacent cells translates into
\begin{equation} \frac{\partial \E}{\partial\mathbf{n}} = \frac{\partial \HH}{\partial\mathbf{n}} = 0 \quad\text{(Neumann boundary conditions).} \label{eq_bcond_neumann}\end{equation} 
Here, the operator $\partial/\partial \mathbf{n}$ denotes the derivative along the normal to the boundaries.

Other special cases are observed when $\KK$ corresponds to some of the high-symmetry points of the Brillouin zone. For the simple case of a two-dimensional square lattice, the $\mathbf{M}$-point corresponds to the adjacent cells having exactly opposite phase of the fields. This requires that the field amplitudes at the boundary are zero:
\begin{equation} \E = \HH = 0 \quad\text{(Dirichlet boundary conditions).} \label{eq_bcond_dirichlet}\end{equation} 

The $\mathbf{X}$ and $\mathbf{Y}$ points are combination of the Dirichlet and Neumann boundary conditions for two different directions.
As a generalization of the above mentioned criteria, the boundary conditions can be found for any point in the Brillouin zone, provided that a correct phase difference between all opposite faces is conserved.

The stronger the coupling between the neighbouring cells, the bigger the difference between the resonant frequencies for the $\mathbf{\Gamma}$ and $\mathbf{X}$ points. Returning to the definition of the group velocity in Eq. (\ref{eq_vg}), one comes to the already expected fact that the energy transfer is proportional to the strength of the inter-cell coupling.

%}}}
%%% \paragraph{Spatial dispersion in photonic structures} %{{{  TODO develop
%%% Near the $\mathbf{\Gamma}$, $\mathbf{\Gamma}$K~or $\mathbf{\Gamma}$M points, the spatial dispersion must be strong to account for the zero group velocity, and actual bending of the dispersion curve downwards in the adjacent zone. The effect can however not be directly interpreted as the negative group velocity, as the dispersion curve only corresponds to one term in the Fourier expansion of the field shape. 
%%% %}}}

\subsection{Band gaps} 
\paragraph{Properties of band gaps} %{{{
In the plots of dispersion curves, one can identify ranges of frequencies which are not associated with any dispersion curve. 
They can be observed for a local dielectric (Figs. \ref{fg_dcsimpleel} and \ref{fg_dcll}), for nonlocal dielectrics (Figs. \ref{fg_dcll_nl}, \ref{fg_dcllactivity}), as well as for periodic structures in Fig. \ref{fg_phcfolding}. Such regions are denoted as (electromagnetic) \textit{band gaps}\index{band gap}, or sometimes also \textit{stop-bands}\index{stop bands|see {band gap}}. They are a very common phenomenon both in dispersive homogeneous media and periodic structures. 
While the wave propagation received due attention in the above, also the wave "non-propagation" deserves to be commented.
%It is essential to point out the similarities and eventual differences between them.

The band-gap behaviour is the simplest to be described in isotropic media, analogously to the one-dimensional case: for a frequency within the band gap, no wavevector $\KK$ can be found with which the electromagnetic wave would propagate, or more precisely, there is no real wavevector for which the (Bloch's) wave would be a solution of Maxwell equations. 

In the more complicated cases of a harmonic plane wave in anisotropic homogeneous media, or of a Bloch's wave in any periodic structure, the band-gap properties depend on the field polarisation and on the direction in which the wave propagates. Many such media, for example, exhibit a band gap for one polarisation and only in some set of wave vector orientations.

%}}}
\paragraph{Reflection from an interface with a band-gap medium} %{{{
A familiar scenario is that a wave propagates through a transparent medium (e.g. glass) and impinges upon an interface with a band-gap medium (e.g. metal). In the latter, it is not allowed to propagate, and all its energy reflects back into the first medium. The Maxwell equations however require that the fields do not end abruptly at the interface; in particular, it follows from the linearity of the medium that a so-called \textit{evanescent}\index{wave!evanescent} wave amplitude decays in the second medium in an exponential manner with the distance below the interface. % TODO unclear: make a distinction between plasma-like medium, PBG and absorptive media

If losses are present, a part of the impinging wave energy may be also dissipated. Nonetheless, the exponential decay of the evanescent wave should not be confused with the exponential nature of the Lambert-Beer law in absorbing media. 
% todo this is not true--  The distinguishing property of the lossless evanescent wave is zero phase change in the direction perpendicular to the interface, corresponding pure imaginary projection of $\KK$ to the interface normal.
% todo discuss that in PBG of a PhC, there can be  also the phase shift...


%}}}
\paragraph{Radiating dipole inside a band-gap medium} %{{{
A less usual, but equally instructive, scenario is when a dipole source is embedded in an infinite medium and is forced by some external mechanism to oscillate at some frequency inside the band gap. Intuitively, one could expect such a forced dipole must inevitably radiate energy into the medium. In reality, the effect of the medium surrounding the dipole will make the dipole act, with regard to the driving force, as a purely elastic or inertial load. 

In particular, if the dipole is realized as an antenna driven by a high-frequency circuit, it will appear to the circuit as purely capacitive or inductive load. Thus, for continuous oscillations, all energy will be returned to the driving mechanism. This way, the environment can enhance or suppress the radiation from an excited atom or molecule, which is known as the \textit{Purcell effect}\index{Purcell effect}.

%To prevent the energy from being returned to the source, one can further imagine the energy is fed into the dipole in a very short pulse, and then the dipole completely disappears. However, the nature of the band gap cannot be overcome: from the convolution-multiplication theorem, such a pulsed source would have a broad enough spectrum to radiate the energy in the neighbouring allowed bands. 

% }}}
\paragraph{Band gaps stemming from individual resonances}%{{{
One common class of band gaps is based on \textit{individual resonances}\index{resonances!individual} in the medium or structure. These band gaps are often observed in realistic natural media, where the resonances arise from vibrations of electrons, atoms or molecules, as described above. 
Although realistic crystalline media are composed of a lattice of individual atoms, which would suggest that even here the wave has to be described in the Bloch's form of Eqs. (\ref{eq_bloch}) and (\ref{eq_blochh}), the optical wavelengths are roughly $10^{4}$ times larger than the inter-atomic distances, and in Eq. (\ref{eq_pw}) the electromagnetic wave can be satisfactorily approximated as a harmonic wave with a well-defined wavelength. To make this harmonic approximation, it is important that these particles be small and close to each other compared to the wavelength. % ; otherwise  % 

Periodic structures can similarly exhibit a resonance confined in the unit cell. 
To the author's knowledge,  %% XX?
every possible shape formed from a dielectric or metal has multiple resonance modes,  %% ... high-permittivity and...
characterized by a relatively high localisation of energy close to the structure. %% ... More specific rules will be given in the Results.
Arranging such resonant elements with a translational symmetry is a fundamental step in building a periodic structure with resonant behaviour. When such a structure is viewed as a homogeneous medium, its dispersion for the Bloch's wave follows a resonance curve similar (but not exactly the same) to that in Fig. \ref{fg_oscillator_spectrum}. Sometimes the behaviour of the Bloch's wave is described as \textit{macroscopic}, while the fields localized in the resonant elements as \textit{mesoscopic} \cite{felbacq2005theory} and finally the word \textit{microscopic} is left for the much finer field in the constituent atoms or molecules.

Another important characteristic of individual resonances consists in a \textit{qualitative} change in the field shape between the lower and upper edge of the band gap. This is true for natural media as well as for lattices of individual resonators. Stipulated by Eq. \ref{eq_harm_osc_result} for a harmonic oscillator, the dipoles are in a phase with the wave at the lower band-gap edge, whereas they are in opposite phase at the upper edge.  % TODO note also that these indiv resonances are typical by a resonance shape in the permitivity
The actual resonant frequency is located at some point inside the bandgap, as will be shown on the examples of simple individual resonators, e.g. in Fig. \ref{fg_CutWires_wireradius1u_cutwidth_comparison}c.

Yet another characteristic is whether the resonant-field shape in the structure unit cell can be approximated as a point dipole or multipole, or their combination. Further, the dipole moment may be electric, magnetic, zero or both simultaneously, as given by symmetry of the field. Most types of structures exhibit multiple resonances at different frequencies, each of them having usually either an electric or a magnetic dipole. In structures of lower symmetry, one resonant mode may have both nonzero dipole moments simultaneously; this leads to optical activity. 

Other resonances have both dipole moments zero, due to higher-order rotational symmetry of the fields. Their interaction with the electromagnetic waves is mediated by the electric or magnetic \textit{quadrupole}\index{quadrupole} (or even higher \textit{multipole}\index{multipole} \cite{merlin2009metamaterials}) moment. This interaction is usually weak, however it grows with the wavenumber $k$. Such multipole moments often complement dipole resonances, and are responsible for the spatial dispersion, see Chapter \ref{chap_sd}. %% TODO this should be moved elsewhere...
%% TODO how much should be "predicted" here, and nhow much should be left for the Results section?

%}}}
\paragraph{Bragg band gaps} %{{{ TODO
Band gaps of a different type, denoted as \textit{Bragg band gaps}\index{band gap!Bragg band gap} or also as \textit{photonic band gaps}\index{band gap!photonic} (PBG), are observed exclusively in periodic structures when an integer number of half-waves fits into the unit cell. 
In the above mentioned periodic lattice of isolated particles, the Bragg band gap can therefore be understood as a resonance of the electromagnetic wave in the space \textit{between} the particles, in contrast with the individual resonances where most of the energy is localised \textit{inside} the particles -- or at least very close to them.

The first Bragg band gap is observed in a periodic structure when the wavenumber $K$ reaches the $\mathbf X$, $\mathbf Y$ or $\mathbf Z$ points of the Brillouin zone. For a wave propagating parallel to any axis of a square or cubic lattice with unit cell size $a$, the rule for a band gap is particularly simple:
\begin{equation} K= \frac{\pi m}{a},\quad m\in\mathbb{Z}.  \label{eq_braggbg}\end{equation}
The Bragg band gaps exhibit a constant $m$ over the whole band gap, without any resonance-like shape of the $\Neff(f)$ %$= \KK(f)c/(2\pi f)$ 
curve observed for individual resonances. 

In analogy with the solid-state physics the Bragg band gaps can be called \textit{direct}\index{band gap!direct and indirect} if the dispersion curves at the lower and upper band gap edges have the same wavenumber $K$.  In contrast, the individual resonances are analogous to the \textit{indirect} band gaps, as the of phase increase over an unit cell is different  for frequencies below and above the individual resonance. 
% TODO 2 cos(K d) = Tr(M)    -- \cite Y. Fink, J.N. Winn, S. Fan, C. Chen, J. Michel, J.D. Joannopoulos, E.L.  Thomas, Science 282 (1998) 1679
%% Other types of band gaps
%% In the literature, also other types of band gaps were identified.  zero-average-index band gaps?
% a) the zero averaged refractive index, b) the zero average permeability and c) the zero permittivity gaps
% Interaction between non-Bragg band gaps in 1D metamaterial photonic crystals

The width of the band gap depends on the scattering strength of the particles; it can be  described by the transfer-matrix formalism and analytically computed for a one-dimensional layered system \cite{laktionov2008}. 
In the Results section, it is also shown that for a particular choice of parameters the width of the Bragg band gap can vanish. 

%}}}
%}}} /section
\section{Historical notes on metamaterials and photonic crystals}%{{{
\subsection{One history of three paradigms} %{{{
The scientific progress rarely follows a straightforward, efficient way. Making a brief historical review and following the development is therefore important to understand the current state of research, and possibly also the actual physical theory. 

Moreover, the field covered by this thesis is also burdened by confusion due to the used terminology: objects falling under the definition of \textit{periodic structures} have been denoted in the literature also as \textit{metamaterials}, \textit{photonic crystals}, \textit{nanostructured electromagnetic materials} or \textit{composites}, \textit{photonic band-gap materials}, \textit{artificial dielectrics} etc. Typically, these terms are used without a proper definition, and may and need not \cite{smith2006homogenization} overlap in their meaning. 

The split in terminology is apparently not only a formal issue, as it psychologically divides the scientific community and thus directs the future research.
Physics, and the field under discussion in particular, has developed rapidly in the recent decades. Some of the seemingly novel concepts were in fact conceived much longer ago, as indicated in the following. Sometimes a mere change in terminology seems to have lead the community to disregard, and later re-invent concepts that were already known \cite{sihvola2002electromagnetic},\cite[p. 5]{klingshirn2007semiconductor}. 


Following the author's personal view,
the historical notes on metamaterials and photonic crystals below are organized into three independent paradigms that developed in parallel, and only later they were unified. Each of them, on its own, once might seem to already have given their fruits, but their \textit{unification} appears to have renewed the scientific interest and pushed the research further. 

This chapter can in no way cover all relevant papers and ideas over the last century; much better resources on the history of the field are in \cite{shamonina2007metamaterials,eleftheriades2012transforming, paudyal2013left}, at the web pages of Moroz \cite{moroz} etc. Also the actual designs of the periodic structures are mentioned very briefly, since they will be the subject of the Results section.

% TODO incorporate TODO Note that the concept of metamaterials is not even limited to electrodynamics; macroscopic mechanical structures were demonstrated to exhibit comparable phenomena for acoustic waves. %CITE
% Photonic crystal slabs with band gaps for electromagnetic and acoustic waves are also investigated %CITE O. Painter's lab
% Staying in the electrodynamics field, metamaterial research is not limited to spatial structures only; so called \textit{metasurfaces} are developed to widely engineer the impedance of the surface, to enhance nonlinearity and to construct absorbers or high-contrast filters with tunable properties.

%}}}

\subsection{Photonic band-gap structures} 
\paragraph{One-dimensional band gaps and dielectric mirrors}%{{{ ======================== 1 cite "diel mirr"
In 1887, Lord Rayleigh \cite{johnson2003introduction} noticed the light reflecting in the volume of a transparent crystal composed of thin periodic layers of slightly different optical properties. X-ray scattering on crystals was observed by Max von Laue \cite{friedrich1913interferenzerscheinungen} and explained by W. H. and W. L. Braggs in the 1910s \cite{bragg1913reflection}. 
Both effects are enabled by the constructive interference of the waves reflected from the macroscopic interfaces in the Rayleigh's crystal, or directly from the atomic layers in the X-ray case. This phenomenon became known as the \textit{Bragg band gap}\index{band gap!Bragg band gap}.

The very same principle was later %% who was the first?
employed for the design and fabrication of \textit{dielectric mirrors}\index{dielectric mirrors}, which are made by stacking multiple layers of two alternating dielectrics with different wave impedances and quarter-wavelength thicknesses. They find a wide variety of applications, owing to their advantages over classical metallic mirrors such as higher reflectance and angular and spectral selectivity.  %% this is related to the \textit{structural colours}\index{structural colours} found in nature -- who and when did it explain?

Embedding one layer different from the others, a \textit{defect}, into the middle of a dielectric mirror can introduce a narrow transmittance window within its photonic band gap. The corresponding \textit{defect mode}\index{defect mode} is characterized by a strong concentration of the electric field $\E$ in the defect, enabling one to construct narrow-band photonic filters that can be tuned by the defect properties \cite{nemec2005highly}.

%}}}
\paragraph{Photonic crystals}%{{{ ============================= find 1 cite "phc"
\begin{figure}[t] \caption{Two approaches for efficient fabrication of a photonic crystal: \textbf{(a)} dielectric rods stacking into a \textit{woodpile}\index{woodpile structure} structure, \textbf{(b)} top view and side view of soft X-ray litography inscribing a \textit{yablonovite}-like hole lattice into a negative photoresist. Both drawings adapted from patent applications  \cite{ho1994periodic} and \cite{sweatt2005method}, respectively.} \label{fg_phc_patents} \centering \includegraphics[width=\textwidth]{img/patents/phc_patents.pdf} \end{figure}
The success in prohibiting the light to propagate along one direction inspired the design of two- and three-dimensional \textit{photonic crystals}\index{photonic crystal!history}, named in analogy with the band gaps for the electron wave in natural crystals \cite{joannopoulos2011photonic}. A nontrivial accomplishment of manufacturing the first three-dimensional structure with a delimited range of frequencies for which there exists no propagating wave in any direction, a \textit{complete band gap}\index{band gap!complete}, is attributed to Yablonovitch \cite{yablonovitch1987} in 1987.
%% phc term first used in ??

The structure, known as \textit{yablonovite}\index{yablonovite}, was designed for microwaves, since it was made by precise mechanical drilling into a plastic cube. Different technologies were developed in the 1990s in attempts to fabricate photonic crystals operating  efficiently in the optical range, e.g. a similar approach based on X-ray nanolitography (see Fig. \ref{fg_phc_patents}a), stacking of structured layers (Fig. \ref{fg_phc_patents}b), direct laser writing in a resin or self-assembly of dielectric spheres into opal-like lattice.

The demands for achieving a complete three-dimensional band gap can be reduced when a singly or doubly periodic modulation is etched into an index-guiding dielectric slab, as was developed by Zengerle in the 1970s and 1980s \cite{zengerle1987light}. Planar structures are also much easier to fabricate using existing technologies such as photolitography.  
% "Any dielectric material can be used to fabricate a photonic crystal with a sizeable and robust complete photonic bandgap (CPBG) in three dimensions, as long as small metal inclusions can be added." \cite{moroz2002metallo}

%}}}
\paragraph{Band gap used in fibres and waveguides} %{{{ ============= (cites: DBR)
Embedding a linear defect into a photonic crystal with a complete band gap forms a waveguide from which the wave cannot escape, even when it is bent under sharp angles.
Likewise, the \textit{planar band-gap waveguide}\index{band gap!in a planar waveguide} can be formed by a linear defect in a photonic crystal slab, which restricts the light by the total internal reflection horizontally and by a means of the two-dimensional band gap vertically. Such structures can form waveguides, beamsplitters and filters in integrated optics. 

The \textit{photonic-crystal fibre}\index{photonic crystal!fibre} (PCF), conceived by Russel in 1991 \cite{russell2007photonic}, differs from the standard index-guiding fibres in that the light is confined to its core by a photonic band gap in both lateral directions. 
Its special feature is the strong guiding of waves even if its core is hollow, that is, with refractive index lower than the cladding. Broad-band single mode operation can be achieved, and the high field concentration can lead to very strong nonlinear interaction.

As a complementary concept, a one-dimensional photonic band-gap structure can be inscribed into conventional index-guiding fibres, forming a \textit{distributed Bragg reflector}. Its filtering capabilities can be made extremely narrowband, which is employed in fibre lasers, telecommunication, sensing, etc. % TODO cite, when first?

Within the sole paradigm of band-gap engineering, new concepts emerge even in the last decades: Obviously, a minor shift of the band gap has a strong impact on whether the light can propagate. Such a shift can be introduced e.g. by static magnetic field in low symmetry PhCs, resulting in enhanced magneto-optic interaction for one-dimensional light propagation \cite{vanwolleghem2009unidirectional}.

%}}}

\subsection{Homogeneous media with uncommon parameters} 
\paragraph{Negative parameters}%{{{ 
By \textit{negative parameters}\index{permittivity!negative}\index{permeability!negative}\index{index of refraction!negative} we  understand either a negative real value of permittivity $\epsrl(\omega)$, permeability $\murl(\omega)$, 
or a negative index of refraction. In transparent media without spatial dispersion, negative index of refraction $\Neff'(\omega)$ requires $\epsrl(\omega)$ and $\murl(\omega)$ to be negative simultaneously, at least if its imaginary part is supposed to be small \cite{pazoutova2011dp}.  %% elaborate this a bit, using complex multiplication? Explain why negative branch is chosen, but this perhaps should be put above in the discussion of group velocity.
Note that all theoretical papers cited in this subsection used the electrodynamics of a \textit{homogeneous} medium, without imposing its microscopic periodicity nor any other technical way how such a medium should be obtained. 

Media with \textit{negative permittivity} are in fact very common in the nature, including metals, doped semiconductors and plasmas; additionally, negative permittivity also results from lattice or electronic vibrations, forming the so called \textit{reststrahlen}\index{reststrahlen bands} bands. Reflection from metals has been explained  % when? by whom?
as the opposite signs of $\epsrl(\omega)$ and $\murl(\omega)$ causing the wave to become \textit{evanescent}\index{wave!evanescent} and rapidly decaying under the surface. %A particular interest in the same phenomenon arose with 
The same mechanism causing the radio wave reflection from the dilute plasma in the Earth's ionosphere was predicted as early as in 1839 by Gauss; its application for over-the-horizon communication was suggested in 1901-1902 by Heaviside and Kennelly.

\textit{Negative permeability}\index{permeability!negative} manifests itself similarly to the electrodynamic point of view, but is observed much less often, typically above magnetic resonances, \textit{magnons}\index{magnons!bulk}.

Theoretical investigation of electrodynamics in media with a \textit{negative index of refraction}\index{index of refraction!negative, history of}, or slightly more generally, those with antiparallel group and phase velocities (cf. Chapter \ref{chap_vfvg}), can be traced back \cite{agranovich2006spatial} to early discussions of anomalous dispersion near resonances by Lamb or in a 1904 book from Schuster \cite{schuster1904introduction, boardman2005negative}. Mandelshtam's lectures on optics illustrated negative refraction of light in the 1940s \cite{mandelstam1971lectures} and were further complemented, among others, in 1957 by Sivukhin's discussion of energy concerns \cite{sivukhin1957energy} and Pafomov's notes about reversed \v{C}erenkov radiation \cite{pafomov1956cerenkov}. 

It was Veselago's review of the negative-index homogeneous media from 1968 that received later significantly wider attention in the literature \cite{veselago1968}. Veselago used % is he first?
the term of \textit{left-handed media}\index{media!left-handed} and speculated to what extent they might exist in the nature, but he did not discuss the way of realizing such media, nor the more general theory of the spatial dispersion which may also lead to natural negative refraction.
% add a note about negative index and spatial dispersion in homogeneous media? 
% (todo) The electrodynamics of spatially-dispersive media ...
% Or in an extra paragraph later, under 2nd unification?

%}}}
\paragraph{Terminology of media with unusual parameters}  %{{{
In the literature, different terms are used for a similar class of media. In most cases it appears \cite{lindell2001bw} that it is a result of disunited terminology rather than a need to distinguish fine nuances. 

\begin{enumerate}
\item{The most general term is \textit{negative-refraction} media. %% TODO accumulate literature
	However, negative refraction can occur due to a range of phenomena, for example, in the anisotropic crystal of calcite under well-chosen incidence angle \cite{agranovich2006spatial}.} 
\item{The aforementioned term is most probably supposed to mean \textit{negative-refractive-index}, or simply \textit{negative-index} media, which by definition behave as (approximately) isotropic at least in a limited range of incidence angles, and therefore their index of refraction can be defined as it is discussed in Sect. \ref{disp_rel_local_media}. 
	
Isotropy also implies the phase and group velocities are parallel or anti-parallel, which legitimates to call the latter group as \textit{backward-wave} \cite{lindell2001bw}.} 

\item{Yet a narrower group is defined by the terms of \textit{left-handed} \cite{veselago1968} or \textit{doubly-negative} media. Their etymology is based on the spatial orientation of the vector triplet $(\E, \HH, \KK)$. In ordinary media with $\epsrl > 0$, $\murl > 0$, the pseudovector $\HH$ is chosen so that the triplet can be associated to the thumb, index and middle finger of the right hand, in this order. 

By simultaneous reversal the sign of both $\epsrl$ and $\murl$, the electric and magnetic induction change their sign to be antiparallel to the fields. The wavevector $\KK$ is reversed, associating the triplet with the left hand. 

In the author's view, both terms can therefore be associated exclusively with a subset of the \textit{negative-index} media, for which the local \textit{effective permittivity} $\eeff(f)$ and \textit{effective permeability} $\meff(f)$ make physical sense. % todo definition of eeff and meff here?
}
\end{enumerate}
% TODO? add about "plasmonic" and epsilon-negative (ENG), or single-negative (SNG) MMs

%}}}
\paragraph{Prediction of focusing effects by a negative-index slab} %{{{  ================== OK? (readpaper, +img)
As also discussed by Veselago \cite[p. 511]{veselago1968} and a decade later by Silin \cite{silin1978possibility},  %: Silin in 1978  %% TODO READPAPER, does not discuss amplification of evanescent waves?
a plane-parallel slab with an isotropic negative index of refraction and sufficient thickness would refract waves emanated from a point source towards a new focus in its volume, and similar negative refraction forms a second focus behind the slab. It is thus often referred to as a \textit{negative-index lens}, although neither its focal point nor optical axis are defined as in classical lenses.

Clearly, a negative-index slab would be free of spherical aberrations, allowing a wide-angle and high numerical aperture for the imaging. To the knowledge of the author, the papers from this time did not discuss its application for overcoming the diffraction limit nor of amplification of evanescent waves.  % should the thin silver lens be noted here, or elsewhere?

%}}}
\paragraph{Media with parameters close to zero}%{{{ ================= (optional todos)
Media where $\epsrl(\omega)$ or $\murl(\omega)$ have very low values compared to vacuum form another class uncommon in nature. They are often denoted as \textit{epsilon-near-zero}\index{permittivity!near zero} (ENZ) and \textit{mu-near-zero}\index{permeability!near zero} (MNZ).
% and perhaps note superhigh opt activity?
The result of $\epsrl(\omega)\sim 0$ or $\murl(\omega)\sim 0$ is also the index of refraction being close to zero, $N(\omega) \sim 0$, particularly when both conditions are satisfied simultaneusly. 

Although e.g. the physics of the ionosphere involves regions where the permittivity transits from negative to positive values, the scientific interest in designing such media significantly grew
% with the paper of ...todo.... published  in ...
after the merging of the paradigm with metamaterials, and therefore such media are often denoted syncretically as \textit{zero-index metamaterials}\index{zero-index metamaterials}. This is a somewhat confusing term, as the papers focusing on their peculiar macroscopic properties \cite{basharin2013epsilon} 
usually do not discuss any structuration and are thus also applicable to fully homogeneous media of such properties.

A feature typical of this class of media is a very small wavevector $K\sim 0$, which can be approximated as the whole $N\sim 0$ region oscillating in phase. The wave entering from air into a $N\sim 0$ medium at near-normal incidence refracts under large angles from the normal; for angles higher than a small critical angle, the wave undergoes the \textit{total external reflection} \cite{schwartz2003total}. 
% note: this has been observed in thinned waveguides, too
% Epsilon-near-zero metamaterials and electromagnetic sources: Tailoring the radiation phase pattern

%the 2000s have brought also the \textit{near-zero}\index{near-zero} effective parameter structures into the limelight, which are also unusual in nature
% electromagnetic nihility?
%add some notable papers; quotes may be in: Epsilon-Mu Zero Meta-material with Zero Refractive Index An Electromagnetic Nihility, Shantanu Das

%}}}
\paragraph{Transformation optics and electromagnetic cloaks} %{{{ ================= OK (cites)
Slow spatial changes in the refractive index cause the electromagnetic wave to follow a curved path, which can be alternatively viewed as propagation through a particularly curved (transformed) space filled with homogeneous medium. This is the basis of \textit{transformation optics}\index{transformation optics}.

The \textit{electromagnetic cloak}\index{electromagnetic cloak} is a hollow structure that guides the electromagnetic waves around its core and reconstructs the wavefront behind it, so that it casts no shadow, independent of the wave direction and of an obstacle present inside the cloak. 
The curved path used to evade the obstacle is always longer than a line segment, however, so typically a medium with continuously varying $\Neff\in \langle0,1\rangle$ must be properly designed by transformation-optics calculations \cite{eleftheriades2012transforming}.
%todo \cite{alitalo2009electromagnetic}
Electromagnetic cloaking % by homogeneous media ???
was considered in 1961 by Dollin \cite{dollin1961possibility} and later by Kerker \cite{kerker1975invisible}.

The concept of cloaking was realized \cite{schurig2006metamaterial} in 2006 using a microwave metamaterial, and followed by many others.
All reported electromagnetic cloaks based on $N<1$ media are inevitably dispersive, which limits their band width; they are also usually polarisation sensitive. The presence of any macroscopic cloak made of known materials is also revealed by its losses, at least at the infrared or optical frequencies. Obviously a great conceptual breakthrough would be needed to overcome these issues. 

A simplified task is to construct the \textit{carpet cloak}\index{carpet cloak} which conceals the obstacle near a surface \cite{valentine2009optical}.
The carpet cloak can be built with $N<1$ media as in full-angle cloaks, however for a single light polarisation, it can also be made  of virtually loss-less anisotropic dielectrics with quite usual values of parameters \cite{wang2013homogeneous}.

%}}}

\subsection{Artificial dielectrics and metamaterials} 
\paragraph{Analysis of nonresonant sub-wavelength structures} %{{{ ================== 3x cite
The third independent paradigm was based on arranging sub-wavelength particles into an optically dense lattice to obtain some desired macroscopic behaviour of the light. 

Its original source of inspiration can be traced to the late 19th century, when the scientific community resolved the question of how the dipoles of individual atoms in matter affect the macroscopic permittivity and permeability. 
The atomic size is negligible compared to the optical wavelengths, so all optically transparent natural media can be easily approximated as homogeneous. The explanation of the diamagnetic behaviour of some materials was proposed in the form of  each atom behaving as a microscopic conductive loop in 1853 \cite{weber1852relationship}. However, the inter-atomic electric dipole interactions are not negligible, and except for gases they preclude simple dipole averaging, so the explanation of the permittivity was formulated two decades later using the Lorentz-Lorenz (or, Clausius-Mossoti) formula.

Theories of dielectric behaviour of composites having inhomogenities larger than the atomic scale have been developed throughout the 20th century, among which the Maxwell-Garnett (1904) and Bruggemann (1935) are most known among multiple approaches. 
Their common assumption is that the wavelength is much larger than the particles; furthermore, they suppose there are no internal resonances in the embedded particles. More elaborate theories are needed when the particle have unusual shapes, particularly near the so-called percolation threshold when they become connected.
%% "process of replacing a complex structure of subwavelength sized components with an “effective medium” with uni-
%% form properties. It is a fundamentally important notion which can be traced back to the earliest days of electro-
%% magnetic theory, to the Lorentz-Lorenz and Maxwell-Garnet effective medium models [1–3]
%			[1] J. C. M. Garnett, Phil Trans. R. Soc. A 203, 385 (1904).
%			[2] D. E. Aspnes, Thin Solid Films 89, 249 (1982).
%			[3] W. Cai and V. Shalaev, Optical Metamaterials: Fundamentals and Applications (Springer, New York, 2010).

%}}}
\paragraph{Synthesis of nonresonant sub-wavelength structures} %{{{ ====================== OK
\begin{figure}[h] \caption{Three examples of artificial dielectrics designs patented half a century ago: \textbf{(a)} sheets containing cut wires of alternating orientation and resonant frequency, enabling one to manipulate the microwave polarisation, \textbf{(b)} cross-section through a lightweight lens of 3 m diameter designed for 50-500 MHz radio waves, made of partially metallized plastic layers, \textbf{(c)} artificial dielectric made of non-resonant cut wires arranged into a nearly isotropic 3-D lattice. Drawings adapted from patent applications \cite{wickersham1960artificial},  \cite{anderson1967artificial} and \cite{hannan1966artificial}, respectively.} \label{fg_mm_patents} \centering \includegraphics[width=\textwidth]{img/patents/mm_patents.pdf} \end{figure}

The design of inhomogeneous structures that could be treated as homogeneous with some desired properties has also a surprisingly long history. 
Inspired by the Lorentz-Lorenz formula, Lord Rayleigh theoretically elaborated the wave propagation through a rectangular lattice of cylindrical or spherical particles in 1892. He also noted the analogy between electromagnetic and acoustic waves \cite[p. 498]{rayleigh1892}. His thoughtful analysis was however still limited to the \textit{static limit}\index{limit!static} of frequencies well below any resonant frequency of the particle, and to the \textit{long-wavelength limit}\index{limit!long-wavelength } well below the first photonic band gap. 

Subwavelength wire arrays have been experimentally used at least since 1898 as microwave polarisers. The concept of assembling a microwave polarisation rotator from organic or metallic helices comes from that era, too \cite{bose1898rotation, emerson}.

Periodic structures denoted as \textit{artificial dielectrics}\index{artificial dielectrics} \cite{brown1955properties} found their use since the 1940s to engineer the broad-band permittivity tensor for the applications in lightweight lenses for microwave frequencies \cite{kock1948metallic}. Initially, they were operated in the static limit below any resonance. Usually the particles consisted of hollow metal waveguides or spheres/disks which increased or reduced the wave phase velocity, respectively. 

%}}}

\subsection{First unification: metamaterials with uncommon parameters} 
\paragraph{Negative effective permittivity} %{{{ ============================= OK
The 1950s brought the first step towards achieving negative effective parameters from man-made periodic structures: The geometry of the wire antennas was chosen to put the fundamental resonance of the electric dipole \textit{below} the frequency of operation \cite{rotman1962plasma,boardman2005negative}. The spectrum of \textit{effective permittivity}\index{permittivity!effective} $\eeff(\omega)$ of this structure had the shape similar to the Lorentz oscillator shown in Fig. \ref{fg_oscillator_spectrum}. It formed a frequency range of values that were negative or at least lower than one. This feature lead to the occasional use of the term \textit{artificial plasma}\index{artificial plasma} \cite{merkel1973simulation}. 

In all cases of $\eeff(\omega)<1$, the effective refractive index $\Neff(\omega)$ was also lower than that of vacuum \cite{brown1953artificial}, but, as follows from
\begin{equation} N(\omega) := k(\omega) \, \frac{c}{\omega} \equiv \sqrt{\epsrl(\omega)\murl(\omega)}, \tag{\ref{eq_dispeqN} \again} \end{equation}
	it could not become negative, since the effective permeability $\meff'(\omega) > 0$ was still positive and its imaginary part $\meff''(\omega)$ was negligible. The unit cell size was kept much smaller than the wavelength, so the problem of \textit{homogenisation}\index{homogenisation} appeared as settled and the nonlocal effects apparently did not cause much concern.

%}}}
\paragraph{Negative effective permeability and negative refractive index}  %{{{ ================== 
The literature is somewhat ambiguous, but a thorough search shows that the term \textit{metamaterial}\index{metamaterial!history} (MM) was first used on conference by Walser in 1999. Although apparently redundant to \textit{artificial dielectrics}, which was used for decades before\footnote{An interesting search into the terminology usage by means of the \textit{Google Books Ngram Viewer} shows that the notion of  \textit{artificial dielectrics} has been the most used in the 1950s and early 1960s, and then its usage declined, to be greatly exceeded by the usage of \textit{metamaterial} in early 2000s. Simultaneously, the usage of \textit{photonic crystal} has rose an order of magnitude higher.},  the new term of \textit{metamaterial} can be viewed as more fortunate as it does not imply whether the structure is man made or not, nor whether its effective parameters mimic a dielectric or, e.g., a metal.

It remains unclear whether there was an objective cause, or a semantic bias from the term \textit{artificial dielectrics} only, that the periodic structures with magnetic effects received less attention than those with electric ones. Apparently the first theoretical studies thereof by Lewin \cite{lewin1947electrical} and Schelkunoff \cite{schelkunoff1952antennas} date back to the 1940s, but the first realisations of split-ring resonator arrays and other magnetically resonant structures seem to have been published no earlier than in the 1980s. Further details can be found in the Results chapter on pp. \pageref{negn_srr} and \pageref{negn_diel}.

% TODO add about history from SRRs and diel-spheres
The research of artificial dielectrics, or newly, metamaterials, gained on significant popularity after the papers from 1999 by Pendry \cite{pendry1999magnetism} and from 2000 by Smith \cite{smith2000composite}. The former outlined how a negative index of refraction can be obtained in a structure composed of metallic wires introducing $\eeff<0$ and of small metallic loops, known as \textit{split-ring resonators}\index{split-ring resonator!history} (SRR), introducing $\meff<0$. The latter paper described the experimental verification of negative refraction in the microwave range. Thousands of papers followed suite which covered different spectral ranges, demonstrated other MM designs, discussed theoretically when and how the metamaterial can be viewed as homogeneous, and sought for new applications. 

%}}}
\paragraph{Imaging beyond diffraction limit with metamaterials} %{{{ ===================== cite fn
Arguably the most intense excitement was caused by another Pendry's paper \cite{pendry2000negative} from 2000. It stated that a flat lens made of a lossless local medium with a negative refractive index can not only enable imaging, but it should do so with a resolution better than is allowed by the diffraction limit, theoretically down to infinitely small details. This device for \textit{subdiffraction imaging}, denoted as the \textit{superlens}\index{superlens}, acquires this extraordinary feature through the formation of high-amplitude evanescent waves at its rear surface. The resulting extremely sharp image in some distance behind the superlens would then be composed not only of negatively refracted propagating waves, but also of the evanescent waves.

If a metamaterial is used as a negative-index medium, subwavelength imaging requires the unit cells being \textit{finer} than the finest details to be imaged. The requirement is in fact even stricter, as also the spatial dispersion has to be negligible for all wavenumbers needed for the desired resolution. % which could be achieved either with sophisticated engineering of the dispersion, or with massive reduction of the unit cell size, so that all wave vectors are near the center of Brillouin zone. % refer to some Fig. of the PHC BZ
% Aside from spatial dispersion, a great challenge are the dissipative losses. 
%If the wave components with high spatial frequency are to be resolved, they also may not be significantly damped. 
For the formation of highly concentrated evanescent waves at the rear surface, the small energy flow-through and simultaneous high energy concentration require extremely low losses. In contrast, down-scaling of the unit cells inherently requires a higher energy localisation, and this is usually at the expense of the losses growing.

A related concept is the \textit{hyperlens}\index{hyperlens}, an inhomogeneous structure which employs anisotropic negative-index media to couple near fields to waves radiated in the air, thus enabling subdiffractive imaging in the far field. Again, the ultimate resolution of the hyperlens is determined by the spatial dispersion of the metamaterial used.

The problems of dissipative losses in resonant periodic structures were not sufficiently overcome yet. Some quantitative improvements can be easily done, namely thickening the metallic structures and increasing the inductance-to-capacitance ratios of the split-ring resonators (SRR) \cite{zhou2008efficient}. At the near-infrared and optical wavelengths, SRRs are overcome by perforated metallic sheets, so called \textit{fishnets}\index{metamaterial!fishnet}, % todo cite 
and it may be favorable to substitute metals with certain types of oxides with metallic-like permittivity in the optical range \cite{naik2011oxides}.
To the author's knowledge, all demonstrations of metamaterial sub-wavelength imaging were so far either only proof-of-concept, microscopic devices, or were limited to the easier accessible microwave region. % (verify there is no breakthrough in the IR lately?)
%The question remains whether the significant progress of technology required to fabricate suitable deeply sub-wavelength structures does not also eliminate their need, by simultaneously providing much better detectors operating in the near field.

%}}}
\paragraph{Metamaterials without unusual parameters and metasurfaces} %{{{ ============== 
In spite of the metamaterials being the most often associated with superlenses and cloaks, many MM applications \cite{zheludev2010road} do not involve negative or close-to-zero effective parameters \cite[p. 15]{sihvola2002electromagnetic}. Among these are novel antenna designs, sensors, light modulators or devices making use of the enhanced nonlinear interaction. Strong interaction of resonant structures with electromagnetic waves enables one to efficiently manipulate the surface impedance by two-dimensional periodic layers, so called \textit{metasurfaces}\index{metasurfaces}. 
%% Many MMs do not bring negative refraction, other (composites) do, but not by means of electromagnetic resonances: 
%% The tunable negative permittivity and negative permeability of percolative Fe/Al2O3 composites in radio frequency range with

A recently published paper \cite{zhao2009repulsive}, even if it is somewhat controversial \cite{silveirinha2010comment}, presents a metamaterial application where a strong chirality leads to reversing the Casimir force between two close matematerial surfaces, which then becomes repulsive.

\begin{figure}[h] \caption{Author's view of the development in the metamaterial and photonic crystal research over the past century} \label{fg_history} \centering 
	\begin{tikzpicture}[box/.style={draw,rounded corners,text width=4cm,align=center}, node distance=5mm, arrow/.style={rounded corners}]
	\node[box] (a)             {\textbf{Media with negative parameters}};
	\node[box, below=of a] (b) {\textbf{Artificial dielectrics and metamaterials}};
	\node[box, below=of b] (c) {\textbf{Photonic band-gap structures}};
	\node[box, right=of a] (d) {\textit{1950s\\ 1st unification}};
	\node[box, right=of d] (e) {\textit{2000s\\ 2nd unification}};
	\node[coordinate, right=of e] (f) {};
	\draw[->] (a) -- (d);
	\draw[->] (a) -- (d);
	\draw[->] (b.east) -| (d.south);
	\draw[->] (d) -- (e);
	\draw[->] (c.east) -| (e.south);
	\draw[->] (e.east) -- (f);
	\end{tikzpicture}
\end{figure}

%}}}

%% XXX up to here revision of FK 2015-10-01 XXX
\subsection{Second unification: metamaterials with photonic crystals} 
\paragraph{Experimental research}%{{{
Over the 20th century, the paradigms of metamaterials (MM) and photonic crystals (PhC) apparently developed independent of each other. The situation started to change in 2000s, when the metamaterial structures were realized to operate at near-infrared or optical frequencies. As already stated, increasing the operation frequency is intricate; not only is the contemporary three-dimensional submicroscopic fabrication technology somewhat limited in its resolution, but even more importantly, with any choice of available materials the dissipative losses become a major issue at the optical frequencies, particularly when the tight confinement of the field is required. 

For both reasons,  a compromise is usually made at the near-infrared and optical frequencies: the particle size is chosen only few times smaller than the wavelength of operation. 
Realized metamaterial structures can no more be viewed as deeply subwavelength \cite{paul2011reflection}, and classifies at the boundary of metamaterials and photonic crystals. This second unification of paradigms is sketched in Fig. \ref{fg_history}. 
%For example, this is the case of  % TODO 3 examples here!  \cite{yannopapas2005negative} \cite{shamonina2007metamaterials}
%where the unit cell size compared to the wavelength $2\pi a/K \approx $, %todo
%or %todo
%or %todo

%}}}
\paragraph{Theoretical investigation of the structures on the MM-PhC boundary}%{{{
Along with experiments, the interplay between the individual and Bragg-type resonances was also studied on a theoretical basis. The former type of resonances is typical for metamaterials, while the latter for photonic crystals. Typically, one either scans the relative cell spacing (this can be done in simulations \cite{shi2007, dominec2014transition} as well as in experiments \cite{rybin2015phase}),  % and \cite{antonoyiannakis1997mie}?
or tunes the individual resonant frequency keeping the lattice parameters unchanged (e.g. \cite{chakrabarti2012magnetic}). The frequencies of the individual and Bragg resonances have different sensitivities to these parameters, and can thus be identified in the spectra. One such study of a typical structure, the array of dielectric rods, is elaborated in the Results section of this thesis.

%PhCs with metallic/metamaterial inclusions
%Antonnoyiannakis1997-Mie_resonances_and_bonding_in_PhC.pdf
% --- weird, examine
% 

In more abstract studies, the structure of the unit cell with its full electrodynamic behaviour is substituted by a flat homogeneous layer characterized only by its effective parameters which may be negative at certain frequencies. % This is equivalent to examining a one-dimensional structure of alternating layers of positive and negative effective parameters.
This approach is only an approximation, since usually the unit cell behaviour involves its coupling with nonradiative fields, which can no more be taken into account when they are replaced by a layer. The advantage is however in the possibility to use analytic means for understanding of the factors forming the band diagram.

Such an unusual photonic crystal can exhibit a specific \textit{zero-order photonic band gap}\index{band gap!photonic, of zero order} \cite{li2003photonic, chen2006derivation, panoiu2006zero} when the phase advance through the layer stack is zero, as has been confirmed also on an experimental basis \cite{kocaman2011zero}. In contrast, the typical photonic band gaps in positive-index dielectric layers occur when the phase advance is $\pi M, M\in \mathbb{N}$ \cite{joannopoulos2011photonic}. 
The characteristic property of the \textit{zero-order} band gap lies in its insensitivity to the structure scaling. The band gap can thus be observed even when the spacing of the hypothetical PhC is deeply subwavelength.

Another type of a band gap with a different field pattern was predicted \cite{wang2004omnidirectional} for a case when layers with $\epsrl>0, \murl<0$ and $\epsrl<0, \murl>0$  are stacked. %It is assumed that the energy transfer is mediated through near-field coupling of surface plasmons.

%Li2003-PBG_from_stack_of_PIM_and_NIM.pdf
%	hypothetical nondispersive N<0 layers interlaced by air - zero-transmittance disp curve of adjacent bandgaps
%	"avg(n)=0" band gap blocks light even for deeply subwavelength layers
%	they denote the bandgap above N<0 band as "avg(n)=0" band - is it correct?
%	they write that in a bandgap, complex eps and mu are unphysical as there is no absorption
%Panoiu2006-ZeroN_PBG_in_PhC_superlattices.pdf
%	nondispersive layers -> peculiar dispersion curves with v_g->infty
%Kocaman_Panoiu-Zero_order_PBG_in_NIM_PhC_in_near-infrared.pdf
%	experimental zero-avg(N) bandgap, denoted also as zero-order photonic band gap
% 	this band gap is invariant to the wave angle
%Chen2006_Dispersion_of_defect_modes_in_PBG_from_DNG-DPS.pdf
%	defect modes in "avg(n)=0" bands, expression is given for dispersion curves
%Jiang2004-Props_of_1D_PhC_containing_SNG.pdf
%	analytic proof of the zero-avgN gap,  special by: the frequency is independent of the lattice	
% 	demonstrate a numeric example of an uncommon MM shapes
%	??? zero-gap width when 1) impedance is matched AND 2) (opt.?) thickness is equal
%LGWang2004_Omnidirectional_PBG_of_1D_PhC_with_SNG.pdf
% 	zero-avg(N) band gap is invariant to the wave angle
% 	LHM-RHM stacks, ENG-MNG stacks
% 	all assume a metamaterial purveys the negative eff params
%Liang2003-Unusual_narrow_bands_in_LHM_stack.pdf todo
%	
%YihangChen2009-Unusual_bands_of_1D_PhC_containing_SNG.pdf todo
%
%Depine2006-ENG_and_MNG_band_gaps_in_1D_PhC.pdf todo

The common property of structures operating at Bloch's wavelengths of comparable magnitude to the cell spacing (i.e. $2\pi a/K \gtrsim 0.2$) is the strong manifestation of the spatial dispersion. The effective index of refraction then forms a typical resonance curve clipped by the closest pair of Brillouin zone boundaries. In most such cases, the use of the local effective parameters $\eeff(\omega)$ and $\meff(\omega)$ hinders the possibility to obtain physical insight. On the contrary, it leads to a deviation from the expected Lorentz oscillator curves in $\eeff(\omega)$ and $\meff(\omega)$ and formation of \textit{antiresonances}\index{antiresonance} \cite{koschny2003resonant,wallen2011anti} and complex values of $\eeff$ and $\meff$ even in lossless media \cite{li2003photonic}. These problems, however, appear to be only artifacts of using the local theory outside its applicability. Instead it is necessary to adopt the full theory of spatially dispersive media, as sketched in Chapter \ref{chap_nonloc} of this thesis.

%}}}
%}}} /section
\section{The boundary between photonic crystals and metamaterials} %{{{
In the light of the previous %{{{
historical review, it appears that the notions of MMs and PhCs developed, for historical reasons, as different paradigms during the 20th century. It was however argued above that these paradigms and research comumnities are undergoing a process of unification, which raises the question whether there are any objective rules to tell apart MMs and PhCs.
%}}}
\paragraph{The wavelength criterion} %{{{
Originally, the determining parameter was surely the ratio of the unit cell size to the wavelength of the electromagnetic wave (at the given frequency of operation).
Even some recent papers \cite{tsukerman2011nonlocal} define MMs as 
\begin{displayquote}
\textit{(...) artificial periodic structures with features smaller than the vacuum wavelength}.
\end{displayquote}

Studies of MMs and PhCs were both inspired by a wave propagating through the lattice of a natural crystal. The difference was however that in the case of MMs, it was an \textit{optical} wave, whereas in the case of PhCs the inspiration came from the \textit{electron wavefunction}. The wavelength of the valence-band electrons is similar to the inter-atomic distance ($\approx 10^{-10}$ m), always requiring the description by the Bloch's theorem, whereas the wavelength of light is roughly four orders of magnitude larger, enabling one to approximate the crystal as a homogeneous medium.

Analogously, the earlier metamaterials (or, artificial dielectrics) were mostly view\-ed as deeply subwavelength, which provided a clear and easy distinction from all PhCs. It also allowed to neglect spatial dispersion for the wave. The situation however appears to change with realistic MMs scaling to the optical range, since the wavelength becomes of the same order of magnitude as the unit cell size and both regions start to overlap as illustrated in Fig. \ref{fg_mm_phc_diagram}. Besides, the PhC-MM boundary then becomes also frequency dependent.

\begin{figure}[h] \caption{Rough illustration of the ratio of the Bloch's wavelength ($2\pi/K$) to the unit cell size $a$ in metamaterials and photonic crystals} \label{fg_mm_phc_diagram} \centering \includegraphics[width=.4\textwidth]{img/mm-phc-diagram.pdf} \end{figure}
%}}}
\paragraph{The energy criterion} %{{{
The time-averaged energy of a plane wave in vacuum is constant, but any sort of interaction of MMs or PhCs with the wave is inherently connected to an uneven distribution of the electromagnetic energy in the structure. This may become the rationale for a new criterion: In the waves in MMs near a resonance, the energy is highly concentrated inside a part of the unit cell, typically close to sharp edges of metals or within a particle of high-permittivity dielectric. On the contrary, the resonances in typical PhCs do not concentrate the energy, and instead there is a moderately high intensity of the field in the majority of the space between the scattering structure elements.

As a result, one can easily understand that the Bragg band gaps typical of a PhC are more sensitive to the unit-cell spacing. 
By contrast, the high concentration of the energy within a small part of the structure reduces the sensitivity of the parameters to a moderate disorder in the cell position, which has also been viewed as a trait of metamaterials, \cite{peng2007}: 
\begin{displayquote}
	\textit{(the fact that) the randomized positions do not influence significantly the left-handed properties indicates that such composite is different from (...) photonic crystals}
\end{displayquote}
% This is the reason behind the aforementioned parametric scans \cite{shi2007,rybin2015phase,dominec2014transition} through the unit cell spacing, % (todo add refs from AdMonsoriu ...) since they recognize the nature of each single resonance mode, and justify even the problematic structures at the boundary to be viewed either as a MM or a PhC in different parts of the spectra.

%}}}
\paragraph{The phase criterion} %{{{
\label{phasecriterion}
The above criteria are somewhat vague, but there exists a well-defined, rigorous criterion dividing the periodic structures clearly into two groups at a given frequency: When a metamaterial is excited above its resonant frequency, and the resonant energy is well localized, the part of the unit cell volume contains field with the sign opposite to that in the excitation wave. This is similar to the behaviour of the Lorentz oscillator illustrated in Fig. \ref{fg_oscillator_spectrum}. %% TODO check Fig 2.2 
This opposite-field domain is delimited by a surface where either the electric or magnetic field is constantly zero, which will be denoted as the \textit{nodal plane}\index{nodal plane}. For the individual resonance, the nodal plane is closed and there exists some path through the unit cell which does not cross it. The Bragg band gap leads to open nodal surfaces delimiting the opposite-field domain which traverse the unit cell boundaries. % A limiting case is the band gap in a plane-parallel dielectric mirror.

The difference between the topologies of localized and delocalized pair of nodal planes can be interpreted as the boundary between MMs and PhCs. In practice, this well-defined difference manifests itself in the change of the phase along each unit cell, without any obvious impact on the propagation of energy. In more details this is demonstrated in the Chapter \ref{sect_diel_rods_el}  %TODO point to the diel array
by numerical studies of the transition cases in dielectric rod arrays.
%Under this criterion, the decision is frequency dependent. All structures fall into the PhC class in the high-frequency limit.

%}}}
\paragraph{Less fortunate criteria} %{{{
Some sources include the requirement of unusual parameters into the MM definition, but this is hardly acceptable since it would obviously exclude some established MMs from the definition; the previous historical chapter attempts to give several examples. It also mentions that even some natural materials can exhibit unusual parameters.

Also the demand for "artificiality" in some other definitions appears to be somewhat artificial on its own, provided that one is interested in the actual structure properties rather than in its origin:
\begin{displayquote}
	\textit{Metamaterial is an arrangement of artificial structural elements, designed to achieve advantageous and unusual electromagnetic properties} \cite{metamorphose-web} 
\end{displayquote}

Some confusion may arise from the fact that some authors seem to overly broaden the applicability of the popular term of metamaterials. For example, two-dimension\-al\-ly periodic structures, in which a wave interacts with a single layer of unit cells only, have substantially different behaviour and applications, and in the author's view should be denoted as \textit{metasurfaces}.

It is even more unfortunate to associate the concept of MMs with one iconic design of a unit cell. Under no thinkable definition can a metamaterial be made by attaching a single split-ring resonator to a microwave transmission line. 
%As one other of many lamentable examples, it also does not seem appropriate to focus a whole chapter of a metamaterial textbook % TODO CITE PAGES
%on a topic such as longitudinally modulated waveguide filters \cite{}.

Some of other definitions realize an important trait of metamaterials, which is \textit{emergence}\index{emergence in metamaterials} \cite{sihvola2002electromagnetic}. While the complexity of the electromagnetic field grows proportionally to the number of unit cells stacked together, new kinds of relatively simple and understandable behaviour can be observed for a large enough system. The introduction of effective parameters, a behaviour emergent from the structure of each unit cell, is one of the best examples.
\begin{displayquote}
\textit{A metamaterial is a (\ldots) substance whose properties depend on its inter-atomic structure rather than on the composition of atoms themselves,}  (whatis.com)
\end{displayquote}
%    (http://whatis.techtarget.com/definition/metamaterial):
or, in more words,
\begin{displayquote}
\textit{metamaterials are artificially structured materials used to control light, sound and many other physical phenomena. The properties of metamaterials are derived from the inherent properties of their constituent materials, as well as from the geometrical arrangement of these (\ldots)}  (website of the Duke university) % TODO cite
\end{displayquote}
%http://metamaterials.duke.edu/research/metamaterials
In the author's view, such quotes represent brilliant observations of emergence in physics, but they cannot be used as a distinction between metamaterials and photonic crystals. In fact, they apply to both classes of structures.

%}}}
\paragraph{Metamaterials as all structures subject to homogenisation} %{{{
The criteria of wavelength or energy localisation leave many structures on the indefinite boundary. The phase-advance criterion is rigorous, but of little practical use. There are other criteria proposed that appear to fail in some way. This leads the author to the personal view that the difference in one's decision on how each structure is characterized. 

The unifying concept of virtually all MM studies is to describe \textit{how} the wave propagates though them. Their overwhelmingly complicated electromagnetic interaction is approximated by \textit{homogenisation}\index{homogenisation}, i.e. by finding such properties of a homogeneous medium that would behave in similar way.

In contrast, the PhC research focuses on \textit{whether} the wave can propagate at all. When a mirror, filter, waveguide cladding or other structure with a band gap is designed, it is important to prevent or allow the light propagation rather than to achieve some desired effective parameters for the Bloch's wave. 

This also justifies denoting MMs and PhCs as two different \textit{paradigms} rather than two distinct \textit{classes} of periodic structures, and it is the reason why the neutral term of \textit{periodic structure} is used in this thesis. 

With this in mind, the author believes that there is no need for a strict rule for classification. Some structures are then more useful to be subject to homogenisation, and with some other it is more interesting to focus on their band and band-gap structure. 
%anything that is subject to homogenisation - i.e. a structure for which we can apply the Bloch's-Floquet theorem
%and determine the K-vector.

% emergent phenomena
% Wikipedia: https://en.wikipedia.org/wiki/Metamaterial
% "Metamaterials derive their properties not from the properties of the base materials, but from their newly designed structures."

%}}}
%}}} /section