---
title: "Dispersion relation of surface plasmons near photonic band gaps: influence of the interaction with light"
authors: ["V.N. Konopsky", "E.V. Alieva"]
affiliation: "Institute of Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow region, 142190, Russia"
journal: "Journal of Modern Optics"
year: 2001
volume: "48"
issue: "10"
pages: "1597-1615"
doi: "10.1080/09500340110056420"
type: journal-article
site_group: "Photonic band gaps and surface plasmons"
url_abstract: ""
url_pdf: "https://valery.konopsky.com/additional_pdf/Dispersion relation of surface plasmons near photonic band gaps_influence of the interaction with light.pdf"
language: en
source_tex: ""
source_pdf: ""
---
## Abstract

We present an experimental and theoretical study of the photonic band
gap in the propagation of surface plasmons (SPs) on periodically
corrugated surfaces. Our main purpose is to investigate the case where
the band gap width is larger than the energy distance between the SP
dispersion curve for a flat surface and the light line. We introduce a
physical model of the interaction of light waves with SPs and derive an
analytical expression for the SP wavevector near band gaps based on the
coupled-mode approach involving three interacting modes (two of them are
SP modes and one is a light mode). By using the interferometric
measurement we have studied for the first time the SP propagation
parameters in the vicinity of the photonic band gap ($10~\mu$m
wavelength region). The predictions of our theory are in good agreement
with the experimental data.

*Keywords: Surface plasmons, Photonic band gap, Coupled-mode theory.*

\

> **Abstract.** We present an experimental and theoretical study of the
> photonic band gap in the propagation of surface plasmons (SPs) on
> periodically corrugated surfaces. Our main purpose is to investigate
> the case where the band gap width is larger than the energy distance
> between the SP dispersion curve for a flat surface and the light line.
> We introduce a physical model of the interaction of light waves with
> SPs and derive an analytical expression for the SP wavevector near
> band gaps based on the coupled-mode approach involving three
> interacting modes (two of them are SP modes and one is a light mode).
> By using the interferometric measurement we have studied for the first
> time the SP propagation parameters in the vicinity of the photonic
> band gap ($10~\mu$m wavelength region). The predictions of our
> theory are in good agreement with the experimental data.

</div>

# Introduction

In recent years the optical properties of materials that possess a
periodic modulation of their refraction index on the scale of the
wavelength of light have received much attention . Such materials can
exhibit photonic band gaps that are very much like the electronic band
gaps for electron waves travelling in the periodic potential of a
crystal. In both cases frequency intervals exist where wave propagation
is forbidden. Materials with band gaps for propagation of bulk light
waves are called "photonic crystals". Surfaces that produce such a band
gap for propagation of surface modes are called "photonic surfaces" . We
will concentrate our discussion on the type of surface wave known as the
surface plasmon (SP). The SP is a nonradiative transverse magnetic mode
at a metal-dielectric interface . The dispersion relation (the
wavevector of the SP as a function of its frequency) at the metal-air
interface is:
``` math
\begin{equation}
k_{\mathrm sp}= {\omega\over c}\left({
\varepsilon_{\mathrm M}
\over \varepsilon_{\mathrm M}+1}\right)^{1/2}
\; ,

\end{equation}
```
where $\varepsilon_{\mathrm M}= \varepsilon'+i\varepsilon''$ is the
complex dielectric constant of the metal.

The simplest "photonic surface" is an ordinary diffraction grating
having a period $\Lambda$ close to half of the wavelength $\lambda$
($\Lambda\sim\lambda/2$). When the projection of the SP wavevector
${\mathbf k}_{\mathrm sp}$ on the normal to the grooves is equal to
half the Bragg vector of the grating ${\mathbf g}$, that is if
$k_{\mathrm
 sp}\cos(\varphi)=g/2$ ($\varphi$ is the angle between the vectors
${\mathbf k}_{\mathrm sp}$ and ${\mathbf g}$), then Bragg scattering
of SPs takes place. The two counter-propagating SP modes set up a
standing wave and, owing to the different surface charges and field
distributions on the grating surface, associated with the two
standing-wave solution, a bandgap in the dispersion of the mode opens
up.

There are several theoretical approaches to describe the SP band gap. A
short review of these approaches can be found in . We will use the
theory of Mills  (see also ) as a starting point in our consideration.
This theory gives an analytical expression for the SP dispersion
relation in the vicinity of the gaps opened by interaction with the
periodically corrugated surface. Figure 1 shows the SP dispersion curve
for a flat surface (the dashed line) and for corrugated surfaces.

According to Mills the band gap width $\delta\omega_G$ for gratings
with shallow grooves $h$ ($h<\Lambda\equiv 2\pi/g$) is
``` math
\begin{equation}
{\delta\omega_G\over\omega_0}=
{gh\over\sqrt{|\varepsilon'|}}\cos(\varphi)\; ,

\end{equation}
```
and the form of the SP dispersion curve is schematically shown by dotted
line in Figure 1.

But the theory of Mills (and other theories) does not describe the
situation when the band gap width $\delta\omega_G$ becomes larger than
the distance between the SP dispersion curve for a flat surface and the
light line (see Figure 1, the solid line). The energy distance is given
by
``` math
\begin{equation}
{\omega-\omega_0\over\omega_0}=
{ck-
ck
\sqrt{\varepsilon'+1\over
\varepsilon'}
\over  ck\sqrt{\varepsilon'+1\over
\varepsilon' }
}\approx
-{1\over 2\varepsilon'}=
{1\over 2|\varepsilon'|}\; .

\end{equation}
```

Comparing (<a href="#1" data-reference-type="ref" data-reference="1">[1]</a>)
and (<a href="#2" data-reference-type="ref" data-reference="2">[2]</a>)
one can see that the influence of the light line (that is the influence
of the interaction between the SP and light) on the SP dispersion curve
becomes important when $h>(g\cos(\varphi)\sqrt{\varepsilon'})^{-1}$.
If
``` math
\begin{equation}
h>{\Lambda\over 2\pi\cos(\varphi)
\sqrt{|\varepsilon'|}}\simeq
{\lambda\over 4\pi\cos^2(\varphi)
\sqrt{|\varepsilon'|}}\; ,

\end{equation}
```
then we have to modify the analytic expression for the SP dispersion
curve near the band gap. The
inequality (<a href="#3" data-reference-type="ref" data-reference="3">[3]</a>)
must be considered only as a first-order estimation, inasmuch as
increasing the strength of modulation of a corrugated surface not only
results in the band gap becoming wider, but the central frequency is
also reduced (see  for detail). The
inequality (<a href="#3" data-reference-type="ref" data-reference="3">[3]</a>)
is fulfilled even for not very large $h$. For instance, in the visible
region $\varepsilon'_{\mathrm Ag}\simeq -18$ ($\lambda=0.63~\mu$m)
and, therefore, a modification of the theory is necessary for
$h>\lambda/50$ (for $\varphi=0$).

This modification is especially needed at the middle and the
far-infrared regions, where the energy distance between the SP
dispersion curve and the light line is vanishingly small. For example,
at $\lambda=10.6~\mu$m,
$\varepsilon'_{\mathrm Ag,Au,Cu}\simeq -4500$ and
condition (<a href="#3" data-reference-type="ref" data-reference="3">[3]</a>),
for $\varphi=0$ becomes $h>\lambda/850$, so the Mills theory is not
applicable in this case. A knowledge of the real and imaginary parts of
the SP wavevector near the photonic bandgap is of more than only
theoretical interest. It is of vital importance for estimating the
decreasing of the SP intensity on the "photonic surface". If the
photonic surface is used as a mirror for high-power infrared laser
radiation, then the decreasing of the intensity of SPs can increase the
laser damage threshold of such a mirror .

The main objective of this paper is to investigate both theoretically
and experimentally the behaviour of the real and imaginary parts of the
SP wavevector near the band gap under
condition (<a href="#3" data-reference-type="ref" data-reference="3">[3]</a>).

We introduce a physical model of the interaction of light waves with SPs
and use the coupled-mode formalism to obtain the SP dispersion relation
near the gap in an analytical form, with the emphasis on the physical
explanation of the result. The method of interferometric SP spectroscopy
is used for the first time to obtain experimentally the real and
imaginary part of the SP wavevector near the bandgap. We used the
wavelength region $\lambda\simeq 10~\mu$m where a discrepancy between
the theory of Mills and experiment is present even for a very small
$h$.

The plan of this paper is as follows: section 2 presents the theory. We
will start with the theory of Mills and then give our expressions for
the SP wavevector near the band gap. Section 3 describes our
experimental approach and presents experimental results. In section 4
our experimental data will be discussed and compared with the
theoretical results from section 2. We offer some conclusions in
section 5.

# Theory

## Equation for ${\mathbf k_{\mathrm sp}}$ near band gaps from the theory of Mills

The dispersion relation near the center of the gap
($k_0$,$\omega_0$) obtained by Mills  in our notation has the
representation
``` math
\begin{equation}
\left({\delta}\omega-{\omega_0{\delta}k\over
k_0l}\right) \left({\delta}\omega+{\omega_0
{\delta}k\cos(2\varphi)\over k_0l}\right)=
{|\varepsilon'|^2\omega_0^2 D\over l^2
(\varepsilon'^2-1)^2} \; ,

\end{equation}
```
where $\varphi$ is the angle between the SP wavevector and the Bragg
vector of the grating,
``` math
\begin{equation}
l=l(\omega)=\left\{1+{\omega\over
2|\varepsilon'|(|\varepsilon'|-1)}
{{\mathrm d}\varepsilon'(\omega)
\over{\mathrm d}\omega} \right\}_{\omega=\omega_0} \; ,

\end{equation}
```
where ${{\mathrm d}\varepsilon'(\omega)/
{\mathrm d}\omega}$ is the derivative of the real part of the
dielectric constant with respect to the frequency and
``` math
\begin{equation}
D= {k_0^2\over|\varepsilon'|} h^2
\left|(\varepsilon'-1)\cos^2(\varphi) \right|^2\; .

\end{equation}
```

Usually
equation (<a href="#4" data-reference-type="ref" data-reference="4">[4]</a>)
is solved for ${\delta}\omega$ in terms of ${\delta}k$. Our
experimental results are represented as
$k_{\mathrm sp}=k_{\mathrm sp}(\varphi)$ at some fixed $\omega$ and
we have to
solve (<a href="#4" data-reference-type="ref" data-reference="4">[4]</a>)
in these terms:
``` math
\begin{equation}
{\delta k_{\pm}
\over k_0}=-{\delta\omega l\sin^2(\varphi)
\over\omega_0\cos(2\varphi)} \pm
\left(
{l^2\delta\omega^2\over\cos^2(2\varphi)\omega_0^2}
-
{|\varepsilon'|h^2k_0^2\over\cos(2\varphi)
(\varepsilon'+1)^2}
\right)^{1/2}\cos^2(\varphi) \; .

\end{equation}
```
We denote:
``` math
\begin{eqnarray}
{\delta}\omega&=&\omega-\omega_0
=\omega-{gc\over
2\cos(\varphi)}\sqrt{\varepsilon'+1\over\varepsilon'}\; ,
\nonumber\\
{\delta}k&=&k-k_0
=k-{g\over
2\cos(\varphi)}\; .

\end{eqnarray}
```
The solution of
equation (<a href="#4" data-reference-type="ref" data-reference="4">[4]</a>)
in these terms is
``` math
\begin{equation}
k_\pm(\varphi , \omega)=
{g\over 2\cos{\varphi}}-
{l\Delta \sin^2(\varphi)\over 2\cos(2\varphi) \cos{\varphi}}
\pm
{\sqrt{\varrho}\over 2
\cos(2\varphi)}
\; ,

\end{equation}
```
where
``` math
\begin{eqnarray}
\varrho&=&
-4{\left( \varepsilon'
          \over\varepsilon'+1
   \right)^2}
   \kappa^2\cos(2\varphi)
   +
l^2 \Delta^2 \cos^2(\varphi)
\; ,
\nonumber\\
\Delta&=&{2\omega\over c }
\sqrt{\varepsilon'\over\varepsilon'+1}
\cos(\varphi)
-g
\; ,
\nonumber\\
\kappa&=&{g^2h\over 4\sqrt{|\varepsilon'|}}
\; .

\end{eqnarray}
```

Equation (<a href="#8" data-reference-type="ref" data-reference="8">[8]</a>)
can be simplified in the limit $l\simeq 1$, $\varepsilon'+1\simeq
\varepsilon'$ and $\varphi = 0$:
``` math
\begin{equation}
k_\pm(0 , \omega)\simeq
{
g\over 2}
\pm    {i\over 2}
\sqrt{
4\left(
  {
g^2 h \over
4 \sqrt{|\varepsilon'|}
  }
\right)^2
-
\left( {2\omega\over c} - g
\right)^2
}
\; .

\end{equation}
```

Now we will show that a similar dispersion relation can be derived using
the coupled-mode formalism for two interacting SP modes. Later, in
subsection (<a href="#three" data-reference-type="ref"
data-reference="three">2.3</a>), we will extend this formalism to the
case of three interacting modes (two of them are SP modes and one is a
light mode). This enables us to take into account the interaction
between SPs and light waves.

## Coupled-mode approach: two interacting modes

At first, following Yariv , we will obtain expressions for the behaviour
of the wavevector near the bandgap in the case of two interacting modes.
Let us consider two electromagnetic modes with complex amplitudes $A$
and $B$ and wavevectors $k_a$ and $k_b$. These are taken as the
eigenmodes of the unperturbed medium so that they represent propagating
disturbances
``` math
\begin{eqnarray}
a&=&Ae^{i(\omega t+k_a z)}

\; ,
\\
b&=&Be^{i(\omega t-k_b z)}
\nonumber
\; ,
\end{eqnarray}
```
where $A$ and $B$ are constant. Mode $a$ correspond to a left
($-z$) travelling wave while $b$ travels to the right. In the
presence of a surface corrugation power is exchanged between modes $a$
and $b$. The complex amplitudes $A$ and $B$ in this case are no
longer constant but will depend on $z$. They obey the following
relations 
``` math
\begin{equation}
\left\{
       \begin{array}{lll}
{
\displaystyle
{\mathrm d}A\over
\displaystyle
{\mathrm d}z}&=&
\kappa_{ab}Be^{-i\Delta z}
\nonumber\$$\medskipamount]
{
\displaystyle
{\mathrm d}B\over
\displaystyle
{\mathrm d}z}&=&
\kappa_{ba}Ae^{+i\Delta z}
\; ,
       \end{array}
\right.

\end{equation}
```
where the phase-mismatch constant $\Delta$ depends on the wavevectors
$k_a$ and $k_b$ as well as on the Bragg vector ${\mathbf g}$ of
the coupling grating. The coupling coefficients $\kappa_{ba}$ and
$\kappa_{ab}$ are determined by the physical situation under
consideration.

It is extremely convenient to define $A$ and $B$ in such way that
$|A(z)|^2$ and $|B(z)|^2$ correspond to the power carried by modes
$a$ and $b$, respectively. The conservation of the total power is
thus expressed as
``` math
\begin{equation}
{{\mathrm d}\over {\mathrm d}z}
(|B(z)|^2 - |A(z)|^2) =0

\end{equation}
```
which,
using (<a href="#11" data-reference-type="ref" data-reference="11">[11]</a>),
is satisfied when
``` math
\begin{equation}
\kappa_{ab}=
\kappa_{ba}^*\; .

\end{equation}
```

On reducing
system (<a href="#11" data-reference-type="ref" data-reference="11">[11]</a>)
to one differential equation, we remove the exponential dependence on
$z$ and obtain a differential equation of the second order with
constant coefficients
``` math
\begin{equation}
{{\mathrm d}^2 B(z)\over {\mathrm d}z^2}-
i\Delta
{{\mathrm d} B(z)\over {\mathrm d}z}-
|\kappa_{ba}|^2 B(z)=0
\; .

\end{equation}
```
Seeking the solution in the form
``` math
\begin{equation}
B(z)={\mathrm const}\cdot e^{\mu z}
\; ,

\end{equation}
```
we obtain the quadratic equation for eigenvalues $\mu$
``` math
\begin{equation}
\mu^2-i\Delta\mu-|\kappa_{ba}|^2=0
\; .

\end{equation}
```

The general solution
of (<a href="#11" data-reference-type="ref" data-reference="11">[11]</a>),
(<a href="#14" data-reference-type="ref" data-reference="14">[14]</a>)
has the form (excepting the case when $\mu_1=\mu_2$)
``` math
\begin{equation}
B(z)=b_1e^{\mu_1 z} + b_2e^{\mu_2 z}
\; ,

\end{equation}
```
where
``` math
\begin{equation}
\mu_{1,2}={1\over 2}i\Delta \pm {1\over
2}S
\; ,

\end{equation}
```
``` math
\begin{equation}
S=\sqrt{4|\kappa_{ba}|^2-\Delta^2}
\; ,

\end{equation}
```
and the constants $b_1$ and $b_2$ are defined from boundary
conditions.

We take the mode $b$ to be incident at $z=0$ on the perturbed region
which occupies the space between $z=0$ and $z=L$. Since mode $a$
is generated by the perturbation we have
``` math
\begin{equation}
b(0)=b_0 \; , \quad
a(L)=0
\; .

\end{equation}
```
With these boundary conditions we have
``` math
\begin{eqnarray}
b_1
&=&
{
b_0\left(S-i\Delta\right)
\over
e^{SL}\left(i\Delta
+S\right)-i\Delta+
S
}
\; ,
\nonumber\\
b_2
&=&
{b_0\left(S+i\Delta\right)
e^{SL}
\over
e^{SL}\left(i\Delta
+S\right)-i\Delta+
S }
\; .

\end{eqnarray}
```
The wavevector $k'_b$ of the perturbed mode receives an addition
$\delta k_b$ to the wavevector $k_b$ of the unperturbed mode $b$:
``` math
\begin{equation}
k'_b=k_b+\delta k_b\; ,

\end{equation}
```
where $\delta k_b$, in the general case is given by
``` math
\begin{equation}
\delta k_b =-{
{\mathrm d}\over {\mathrm d} z
              }
\left[
\mbox{argument}
\left(
b_1e^{\mu_1 z} + b_2e^{\mu_2 z}
\right)
\right]
\; .

\end{equation}
```
But inasmuch as $|b_1|\simeq 1$, $|b_2|\simeq 0$ for
$\Delta<-2|\kappa_{ba}|$ and $|b_1|\simeq 0$, $|b_2|\simeq 1$ for
$\Delta>-2|\kappa_{ba}|$ we can write
``` math
\begin{eqnarray}
\delta k_b
&\simeq& i\mu_1 \quad \mbox{at} \quad
\Delta<-2|\kappa_{ba}|
\; ,
\nonumber\\
\delta k_b
&\simeq& i\mu_2 \quad \mbox{at} \quad
\Delta>-2|\kappa_{ba}|
\; .

\end{eqnarray}
```

We have obtained the solution in a general form without specifying the
nature of the contradirectional coupling between two modes. Comparing
equations (<a href="#20a" data-reference-type="ref" data-reference="20a">[20a]</a>),
(<a href="#22" data-reference-type="ref" data-reference="22">[22]</a>),
(<a href="#18" data-reference-type="ref" data-reference="18">[18]</a>)
and
(<a href="#18a" data-reference-type="ref" data-reference="18a">[18a]</a>)
with (<a href="#0-8" data-reference-type="ref" data-reference="0-8">[0-8]</a>),
it is apparent that for SP modes, coupled through the grating, the
coupling coefficient $\kappa_{ba}$ is
``` math
\begin{equation}
|\kappa_{ba}|={g^2h\over 4\sqrt{|\varepsilon'|}}
\; .

\end{equation}
```

From equations 
(<a href="#18" data-reference-type="ref" data-reference="18">[18]</a>),
(<a href="#18a" data-reference-type="ref" data-reference="18a">[18a]</a>)
and
(<a href="#22" data-reference-type="ref" data-reference="22">[22]</a>)
one can see that the maximum negative value of the imaginary part of the
SP wavevector is
``` math
\begin{equation}
-{\mathrm Im}(k_2) =  |\kappa_{ba}|
\; ,

\end{equation}
```
and it takes place at $\Delta = 0$. The dissipation in the mode $b$
is due to an energy transfer to the counter-propagating mode $a$. From
the same equations one can see that the maximum addition to the real
part of the SP wavevector occurs at
``` math
\begin{equation}
\Delta=\pm 2 |\kappa_{ba}|
\; ,

\end{equation}
```
(at these points ${\mathrm Im}(k_2)$ becomes zero) and is given by
``` math
\begin{equation}
{\mathrm Re}(k_{1,2}) = \mp |\kappa_{ba}|
\; .

\end{equation}
```
Therefore
``` math
\begin{equation}
{\delta k_G\over k_0}=
{\delta k_G\over (g/2)}=
{2 |\kappa_{ba}| \over (g/2)} =
{gh\over\sqrt{|\varepsilon'|}}\; .

\end{equation}
```

Here it may be noted that $\delta k_G$ it is the difference between
maximum and minimum values of $k$, which occur at the edges of the
*frequency* band gap. Therefore there are no needs to use a term
"$k$-gap" here, because it is really "$\omega$-gap" and $k$
becomes complex value into the gap. We believe that the term "$k$-gap"
must be reserved for conditions such as "a medium with modulated gain",
where real "$k$-gap" takes place, and frequency $\omega$ becomes
complex into the "$k$-gap" .

## Coupled-mode approach: three interacting modes

Now we introduce the third interacting mode — the light mode $c$ —
into the
system (<a href="#11" data-reference-type="ref" data-reference="11">[11]</a>).
The necessity of the introducing of the third mode is explained in
Figure 2. When $|{\mathbf k}_{\mathrm sc}|=|{\mathbf k}_{\mathrm sp}-
{\mathbf g}| \leq \omega / c$ the energy of SP mode $b$ is leaking
out to the light mode $c$, which propagates at a small angle
$\alpha= \pi/2 - \theta$ to the surface ($k_{\mathrm
sc}=(\omega/c)\cos(\alpha)$).

Furthermore, for $\alpha<\delta\theta \sim h/\Lambda$ an *inter*action
between light mode $c$ and the SP mode $b$ takes place. This occurs
because, during propagation along the surface at the near-grazing angle,
the light mode $c$ actually interacts with the modulated medium and
can obtain the momentum ${\mathbf - g}$ directed along the surface and
returns to its original state — to the SP mode $b$.

So, we have three interacting modes $a$, $b$, $c$:
``` math
\begin{eqnarray}
a&=&Ae^{i(\omega t+k_a z)}
\; ,
\nonumber\\
b&=&Be^{i(\omega t-k_b z)}

\; ,
\\
c&=&Ce^{i(\omega t+k_c z)}
\nonumber
\; .
\end{eqnarray}
```
They obey the following relations
``` math
\begin{eqnarray}
\left\{
       \begin{array}{lcc}
\frac{
\displaystyle
{\mathrm d}A}{
\displaystyle
{\mathrm d}z}&=&
\kappa_{ab}Be^{-i\Delta_{ba} z}
\\[\medskipamount]
{
\displaystyle
{\mathrm d}B\over
\displaystyle
{\mathrm d}z}&=&
\kappa_{ba}Ae^{i\Delta_{ba} z}
+\kappa_{bc}Ce^{i\Delta_{bc} z}
\\[\medskipamount]
{
\displaystyle
{\mathrm d}C\over
\displaystyle
{\mathrm d}z}&=&
\kappa_{cb}Be^{-i\Delta_{bc} z}
\; ,
         \end{array}
\right.

\end{eqnarray}
```
where $\Delta_{ba}$ and $\Delta_{bc}$ are the phase-mismatch
constants between the corresponding modes.

The equation for the conservation of the total power now becomes
``` math
\begin{equation}
{{\mathrm d}\over {\mathrm d}z}
(|B(z)|^2 - |A(z)|^2 - |C(z)|^2) \simeq 0

\end{equation}
```
which, using (<a href="#3-11" data-reference-type="ref"
data-reference="3-11">[3-11]</a>), is satisfied when
``` math
\begin{equation}
\kappa_{ab}\simeq
\kappa_{ba}^*  \quad  \mbox{and} \quad
\kappa_{cb}\simeq
\kappa_{bc}^*
\; .

\end{equation}
```

Reducing system (<a href="#3-11" data-reference-type="ref"
data-reference="3-11">[3-11]</a>) to one differential equation, we
remove the exponential dependence on $z$ and obtain a differential
equation of the third order with constant coefficients. Seeking the
solution in the form $B(z)={\mathrm const}\cdot e^{\mu z}$ we obtain a
cubic equation for the eigenvalues $\mu$
``` math
\begin{equation}
\mu^3+\alpha\mu^2+\beta\mu+\gamma =0
\; ,

\end{equation}
```
where the coefficients $\alpha$, $\beta$ and $\gamma$ are
``` math
\begin{eqnarray}
\alpha
&=&
-i(\Delta_{ba}+\Delta_{bc})
\; ,
\nonumber\\
\beta
&=&
-|\kappa_{ba}|^2-|\kappa_{bc}|^2-\Delta_{ba}
\Delta_{bc}
\; ,
 \\
\gamma
&=&
i\left(|\kappa_{ba}|^2\Delta_{bc} +
|\kappa_{bc}|^2\Delta_{ba}\right)
\nonumber
\; .
\end{eqnarray}
```

The general solution of (<a href="#3-11" data-reference-type="ref"
data-reference="3-11">[3-11]</a>) has the form (excepting the cases when
$\mu_i=\mu_j$)
``` math
\begin{equation}
B(z)=b_1e^{\mu_1 z} + b_2e^{\mu_2 z}
+ b_3e^{\mu_3 z}
\; ,

\end{equation}
```
and $\mu_{1,2,3}$ are the solutions of the cubic
equation (<a href="#3-16" data-reference-type="ref"
data-reference="3-16">[3-16]</a>)
``` math
\begin{eqnarray}
\mu_1
&=&
{1\over 6}
Q^{1/3}
-2{p\over
Q^{1/3}}
-{\alpha\over 3}
\; ,
\nonumber\\
\mu_{2,3}
&=&
-{1\over 12}
Q^{1/3}
+{p\over
Q^{1/3}} \pm
{\sqrt{3}\over 2}i
\left(
{1\over 6}
Q^{1/3}
+2{p\over
Q^{1/3}}
\right)
-{\alpha\over 3}
\; ,

\end{eqnarray}
```
where
``` math
\begin{eqnarray}
Q
&=&
-108q+12\sqrt{12p^3+81q^2}
\; ,
\nonumber\\
p
&=&
\beta
-{1\over 3}\alpha^2

\; ,
\\
q
&=&
\gamma-{1\over 3} \beta \alpha
+{2\over 27}\alpha^3
\nonumber
\; .
\end{eqnarray}
```

The constants $b_1$, $b_2$ and $b_3$ are defined from the boundary
conditions. If the boundary conditions are such that the mode $b$ is
incident at $z=0$ on the perturbed region which occupies the space
between $z=0$ and $z=L$, we have
``` math
\begin{equation}
b(0)=b_0
\; , \quad
a(L)=0
\; , \quad
c(L)=0
\; .

\end{equation}
```
In this case
``` math
\begin{eqnarray}
b_1
&=&
{
b_0 \over \Sigma}
(\mu_2-\mu_3)(\Delta_{ba}+i\mu_1)(\Delta_{bc}+i\mu_1)
\exp{[(\mu_2+\mu_3-i\Delta_{ba} -i\Delta_{bc})L]}
\; ,
\nonumber\\
b_2
&=&
{
b_0 \over \Sigma}
(\mu_3-\mu_1)(\Delta_{ba}+i\mu_2)(\Delta_{bc}+i\mu_2)
\exp{[(\mu_1+\mu_3-i\Delta_{ba} -i\Delta_{bc})L]}

\; ,
\\
b_3
&=&
{
b_0 \over \Sigma}
(\mu_1-\mu_2)(\Delta_{ba}+i\mu_3)(\Delta_{bc}+i\mu_3)
\exp{[(\mu_1+\mu_2-i\Delta_{ba} -i\Delta_{bc})L]}
\nonumber
\; ,
\end{eqnarray}
```
where
``` math
\begin{eqnarray}
\Sigma
=&
(\mu_2-\mu_3)(\Delta_{ba}+i\mu_1)(\Delta_{bc}+i\mu_1)
\exp{[(\mu_2+\mu_3-i\Delta_{ba} -i\Delta_{bc})L]}
\nonumber\\
+&
(\mu_3-\mu_1)(\Delta_{ba}+i\mu_2)(\Delta_{bc}+i\mu_2)
\exp{[(\mu_1+\mu_3-i\Delta_{ba} -i\Delta_{bc})L]}
 \\
+&
(\mu_1-\mu_2)(\Delta_{ba}+i\mu_3)(\Delta_{bc}+i\mu_3)
\exp{[(\mu_1+\mu_2-i\Delta_{ba} -i\Delta_{bc})L]}
\nonumber
\; .
\end{eqnarray}
```

The wavevector $k'_b$ of the perturbed mode receives an addition
$\delta k_b$ to wavevector $k_b$ of unperturbed mode $b$ – see
equation (<a href="#20a" data-reference-type="ref" data-reference="20a">[20a]</a>),
where $\delta k_b$, in the general case, is given by
``` math
\begin{equation}
\delta k_b =-{ {\mathrm
d}\over {\mathrm d} z } \left[ \mbox{argument} \left(
b_1e^{\mu_1 z}
+ b_2e^{\mu_2 z}
+ b_3e^{\mu_3 z}
\right)
\right]
\; .

\end{equation}
```

Equations (<a href="#3-17" data-reference-type="ref"
data-reference="3-17">[3-17]</a>),
(<a href="#3-18" data-reference-type="ref"
data-reference="3-18">[3-18]</a>),
(<a href="#3-16a" data-reference-type="ref"
data-reference="3-16a">[3-16a]</a>) and
(<a href="#3-20" data-reference-type="ref"
data-reference="3-20">[3-20]</a>) give the exact solution of the problem
but, instead of these rather unwieldy expressions, it is convenient to
have approximate expressions for some important quantities.

In the approximation
``` math
\begin{equation}
\Delta_{ba} -
\Delta_{bc} =
\Delta_{ac}
\approx
{\omega \over c}
{1 \over |\varepsilon'|} < |\kappa_{ba}|, |\kappa_{bc}|

\end{equation}
```
(this condition is equivalent to
(<a href="#3" data-reference-type="ref" data-reference="3">[3]</a>)) the
maximum negative value of the imaginary part of the SP wavevector at
$\Delta_{ba}=0$ is
``` math
\begin{equation}
\mbox{max}[-{\mathrm Im}(k_2),
-{\mathrm Im}(k_3)]\simeq
\sqrt{
|\kappa_{ba}|^2+
|\kappa_{bc}|^2
}         -{1\over 8}
{
|\kappa_{bc}|^2
(|\kappa_{bc}|^2 + 4
|\kappa_{ba}|^2)
\Delta_{ac}^2
\over
(|\kappa_{ba}|^2 +
|\kappa_{bc}|^2)^{5/2}
}
\; ,

\end{equation}
```
or if $\Delta_{ac}
\sim 0$
``` math
\begin{equation}
\mbox{max}[-{\mathrm Im}(k_2),
-{\mathrm Im}(k_3)]\approx
\sqrt{
|\kappa_{ba}|^2+
|\kappa_{bc}|^2
}

\end{equation}
```
(compare this with
(<a href="#23" data-reference-type="ref" data-reference="23">[23]</a>)).

The real part of $k_2$ in this approximation has the maximum addition
at
``` math
\begin{equation}
\Delta_{ba}\simeq
-2\sqrt{
|\kappa_{ba}|^2+
|\kappa_{bc}|^2
}

\end{equation}
```
(at this point ${\mathrm
Im}(k_2)$ becomes zero) — compare this with
(<a href="#24" data-reference-type="ref" data-reference="24">[24]</a>) —
and takes the form
``` math
\begin{equation}
{\mathrm Re}(k_2)
\simeq
+\sqrt{
|\kappa_{ba}|^2+
|\kappa_{bc}|^2
}  \; .

\end{equation}
```

Inasmuch as at $\Delta_{ba}\simeq
+2\sqrt{
|\kappa_{ba}|^2+
|\kappa_{bc}|^2
}$ the value ${\mathrm Re}(k_3)\approx -\Delta_{ac}/2$ ($\sim 0$)
we have
``` math
\begin{equation}
{\delta k_G \over k_0}=
{\delta k_G \over (g/2)} \approx
{
\sqrt{
|\kappa_{ba}|^2+
|\kappa_{bc}|^2}
\over (g/2)} +{1\over 2|\varepsilon'|}
\;

\end{equation}
```
(compare this with
(<a href="#26" data-reference-type="ref" data-reference="26">[26]</a>)).

Up to this point we elaborate the coupled-mode approach on the
assumption that the SP wavevector $k_b$ is collinear with the grating
Bragg vector $g$ ($\varphi=0$). If $\varphi\not=
0$ then instead of the
equation (<a href="#20a" data-reference-type="ref" data-reference="20a">[20a]</a>)
one should use
``` math
\begin{equation}
|k'_b|\simeq k_b+\delta k_b\cos(\varphi)

\end{equation}
```
(it is true at $\delta k_b<<k_b$), and use the following expressions
for $\Delta_{ba}$ and $\Delta_{bc}$ (see vector diagram on
Figure 3):
``` math
\begin{eqnarray}
\Delta_{ba}
&=&
2k_{\mathrm sp}
\cos(\varphi)-g
\; ,
\nonumber\\
\Delta_{bc}
&=&
k_{\mathrm sp}
\cos(\varphi)
+\sqrt{
\left[{\omega\over c}\sin(\theta_{\mathrm max})\right]^2
-\left[k_{\mathrm sp}
\sin(\varphi)\right]^2} -g

\end{eqnarray}
```
(about $\theta_{\mathrm max}$ see below,
equation (<a href="#54" data-reference-type="ref" data-reference="54">[54]</a>)).

Let us now determine the other values in the above expressions. The
value $\kappa_{ba}$ has been already determined by
equations (<a href="#22a" data-reference-type="ref" data-reference="22a">[22a]</a>).
The coupling coefficient $\kappa_{bc}$ can be derived from the
following consideration: from the
system (<a href="#3-11" data-reference-type="ref"
data-reference="3-11">[3-11]</a>), ignoring the interaction between
modes $b$ and $a$, one can find that if $\Delta_{bc}\sim 0$ (this
takes place at $0<\alpha<\delta\theta$) the complex amplitude of the
mode $b$ changes as follows $B(z)\sim
\exp(-|\kappa_{bc}|z)$. Therefore
``` math
\begin{equation}
|\kappa_{bc}|=\gamma_r/2\; ,

\end{equation}
```
where $\gamma_r$ is the constant of the radiative damping of the SP on
the grating. For $\gamma_r$ we used the expression from :
``` math
\begin{equation}
\gamma_r  \approx
{|a_p|^2\over c}=\left|
\left({|Z|\omega\cos(\theta) \over 2c}\right)^{1\over 2}
hg{\cos(\varphi)\over \cos(\theta)+Z}
         \right|^2
\; ,

\end{equation}
```
where $Z=Z'+iZ''=1/\sqrt\varepsilon_{\mathrm M}$ is the surface
impedance of the metal. $\gamma_r$ may be written in the
form (<a href="#50" data-reference-type="ref" data-reference="50">[50]</a>)
if angle $\varphi$ is not very large
($\sin(\varphi)<\cos(\varphi)$). In order words, this means that a
transformation of the SP to s-polarized light is neglected. Hereafter we
will work at this approximation.

Considering that
``` math
\begin{eqnarray}
Z'&=&{1\over\sqrt{2}}{\sqrt{
\sqrt{\varepsilon'^2_{ }+ \varepsilon''^2_{ }}
+\varepsilon'}\over
\sqrt{\varepsilon'^2_{ }+ \varepsilon''^2_{ }}}
\; ,
\\
Z''&=&-{1\over\sqrt{2}}{\sqrt{
\sqrt{\varepsilon'^2_{ }+ \varepsilon''^2_{ }}
-\varepsilon'}\over
\sqrt{\varepsilon'^2_{ }+ \varepsilon''^2_{ }}}
\nonumber
\; ,
\end{eqnarray}
```
the
expression (<a href="#50" data-reference-type="ref" data-reference="50">[50]</a>)
can be written in the form
``` math
\begin{equation}
\gamma_r  =
{1\over 2c}
{(\varepsilon'^2 +
\varepsilon''^2)^{1/4}\cos^2(\varphi)g^2h^2\cos(\theta)\omega
\over
(\cos^2(\theta)
(\varepsilon'^2 + \varepsilon''^2)^{1/2}
+   \sqrt{2}
\cos(\theta)
((\varepsilon'^2 + \varepsilon''^2)^{1/2} +\varepsilon')^{1/2}
+1)}
\; .

\end{equation}
```
In approximation $\varepsilon''<<-\varepsilon'$
``` math
\begin{equation}
\gamma_r  =
-{1\over 2}
{\sqrt{-\varepsilon'}\cos^2(\varphi)g^2h^2\cos(\theta)\omega
\over
c(\cos^2(\theta)
\varepsilon'
-1)}
\; .

\end{equation}
```
The function
(<a href="#53" data-reference-type="ref" data-reference="53">[53]</a>)
has a sharp maximum at
``` math
\begin{equation}
\theta=
\theta_{\mathrm
max}=\arccos(|Z|) \simeq
\arctan
\left(\sqrt{-\varepsilon'-1}
\right)\approx
{\pi\over 2}- {1\over \sqrt{-\varepsilon'}}
\; .

\end{equation}
```

In order to appreciate a physical interpretation of this maximum one
must recall the Fresnel’s formula for p-polarized light. In the
Leontovich approximation (see, for example, ) it is
``` math
\begin{equation}
R_p=\left|
{\cos(\theta)-Z
\over
\cos(\theta)+Z
}
\right|^2
\; .

\end{equation}
```
One can see that at $\theta=\arccos(|Z|)$
equation (<a href="#55a" data-reference-type="ref" data-reference="55a">[55a]</a>)
has a minimum (this is an analog of the Brewster angle). That is the
minimum in the reflection corresponds to the maximum in the plasmon
emission.

At this point
``` math
\begin{equation}
\gamma_r(\theta_{\mathrm max})  =
\gamma_0  =
{\omega\over 4c}
\cos^2(\varphi)
g^2h^2
\; .

\end{equation}
```
The substitution $\omega/c\simeq g/(2\cos(\varphi))$ yields
``` math
\begin{equation}
\gamma_0  =
{1\over 8}
g^3h^2
\cos(\varphi)
\; .

\end{equation}
```

We will proceed as follows: for $\theta >\theta_{\mathrm max}$, i.e.
for
``` math
\begin{equation}
\varphi<\varphi_0=
\arccos{\left(k_{\mathrm
sp}^2+g^2-(\omega\sin{(\theta_{\mathrm
max})}/c)^2\over 2k_{\mathrm sp}g \right)}
\; ,

\end{equation}
```
(see Figure 3) we will use
equation (<a href="#53" data-reference-type="ref" data-reference="53">[53]</a>)
to determine $\kappa_{bc}$ and put $\Delta_{bc}=0$; but for
$\theta \leq\theta_{\mathrm max}$, i.e. for $\varphi\geq\varphi_0$,
we will use the
equation (<a href="#56" data-reference-type="ref" data-reference="56">[56]</a>)
to determine $\kappa_{bc}$ and will use $\Delta_{bc}$ given by the
equation (<a href="#48a" data-reference-type="ref" data-reference="48a">[48a]</a>)
(such behaviour of the phase-mismatch constant may be described by the
following single analytical expressions:
$1/2\Delta_{bc}-1/2\sqrt{\Delta_{bc}^2+d}$, where $\Delta_{bc}$ was
taken from
equation (<a href="#48a" data-reference-type="ref" data-reference="48a">[48a]</a>)
and $d\rightarrow 0$).

When
condition (<a href="#3" data-reference-type="ref" data-reference="3">[3]</a>)
is satisfied one must use for estimates
equations (<a href="#47" data-reference-type="ref" data-reference="47">[47]</a>),
(<a href="#48" data-reference-type="ref" data-reference="48">[48]</a>)
and
(<a href="#44" data-reference-type="ref" data-reference="44">[44]</a>)
instead
of (<a href="#25" data-reference-type="ref" data-reference="25">[25]</a>),
(<a href="#26" data-reference-type="ref" data-reference="26">[26]</a>)
and
(<a href="#23" data-reference-type="ref" data-reference="23">[23]</a>),
respectively. The physical reason for the appearance of the coupling
coefficient between SPs and light modes ($\kappa_{bc}$) in new
equations is as follows: it is well known that the dominant contribution
to the change of the SP wavevector near the gaps comes from a process in
which the SP with wavevector $\mathbf{k}_{\mathrm sp}$ is scattered
into an intermediate state with wavevector
$\mathbf{k}_{\mathrm sc}=\mathbf{k}_{\mathrm sp}- \mathbf{g}$ and then
backscattered into the initial state by a second interaction with the
periodic surface structure. This process causes a change of the SP phase
velocity and thus also of the SP wavevector. When
condition (<a href="#3" data-reference-type="ref" data-reference="3">[3]</a>)
is satisfied, not only the SP modes, but also the light modes (which
propagate at near-grazing angles
$0<\alpha<\delta\theta \sim h/\Lambda$ to the surface) act as
intermediate states for SPs.

One can see that $|\kappa_{bc}|\sim h^2$, while
$|\kappa_{ba}|\sim h$ and therefore at $h\rightarrow 0$ the
coupled-mode solution for three modes transforms to the coupled-mode
solution for two modes (i.e. becomes similar to the theory of Mills).

# Experiment

## Grating

The diffraction grating was made on a silicon wafer by the method of
photolithography followed by etching, and, was coated by a 380 nm layer
of silver. The period of the grating is equal to
$\Lambda=(5.3008 \pm  0.0005)~\mu$m, the groove depth:
$H=0.24~\mu$m, the groove width: $\Lambda -s=1.41~\mu$m, the grating
dimension is 10$\times$<!-- -->10 mm. The profile of our grating was
recorded by a commercial scanning probe microscope "Solver P-47" of the
"NT-MDT" firm  (using the needles with radius of curvature
$\sim 10$ nm) and it is presented in Figure 4. The value of the
grating period with the accuracy mentioned above was obtained from
diffraction measurements of He-Ne laser radiation.

In our theoretical approach we consider the sinusoidal shape of grating,
whereas the present grating is rectangular. We will proceed as follows:
we resolve the grating in several Fourier components
$y(z)=H s/\Lambda+\sum_n h_n\cos(n g z)$ with amplitudes given by
``` math
\begin{equation}
 h_{\mathrm
n}={2H\over n\pi}\sin{\left(n\pi {s \over\Lambda}\right)}
\; .

\end{equation}
```
The amplitude of the first Fourier component $h=h_1$
($h_1\simeq0.113~\mu$m in our case) is used to find a splitting of the
SP dispersion curve in our theory. The presence of the remaining Fourier
components leads to a uniform shift of this curve. The value of this
shift can be estimated using a theory developed by E. Kretschmann and
E. Kröger :
``` math
\begin{equation}
\Delta k_{\mathrm sp} = \sum_{\mathrm n} {1\over 2}h_{\mathrm
n}^2{\omega_0^3\over c^3}\oint
g({\mathbf k}-{\mathbf k}_0)
A({\mathbf k},{\mathbf k}_0){\mathrm d}^2{\mathbf k}
\; .

\end{equation}
```
where $g({\mathbf k}-{\mathbf k}_0)=0.5(\delta({\mathbf
k}-{\mathbf k}_0-{\mathbf g}_{\mathrm
n})+\delta({\mathbf k}-{\mathbf k}_0+{\mathbf g}_{\mathrm
n}))$ (here $\delta$ is the Dirac function) and $A({\mathbf
k},{\mathbf k}_0)$ given in . In our numerical calculations the first
thirty Fourier components ($n=30$) were taken into account (it is
enough to approximate the real shape of the grating, so one can see from
Figure 4 that radius of curvature of the grating edges $\geq 20$ nm,
while $h_{30}\sim 5$ nm). In this case
$\Delta k_{\mathrm sp}\approx 1.38$ rad/cm.

The direction of SP propagation on the surface of grating was changed by
rotation of the grating. This allowed us to change the angle $\varphi$
(see Figure 3) between the SP wavevector and the Bragg vector of the
grating. The grating holder was mounted on a rotating table capable of
$1/60^\circ$ resolution inside the unit for SP excitation.

## Optical system and measurements

The interferometric phase SP spectroscopy with aperture excitation was
used to measure the propagation parameters of SP on the grating
surface . The scheme of this experiment is shown in Figure 5.

To launch the SP in the $10~\mu$m spectral range on the grating
surface CO$_2$ laser radiation corresponding to spectral lines P(10) –
P(28) and R(10) – R(26) is used. The corresponding wavelengths of the
radiation are in the spectral ranges 937 – 953 cm$^{-1}$ and 969 –
980 cm$^{-1}$. The CO$_2$ laser radiation was focused by a
cylindrical lens on the aperture gap between the sample surface and the
razor blade (see Figure 5). A part of the radiation was transformed into
SP which propagates along the grating surface, another part travels in
the form of a spreading packet of bulk radiation. To convert the SP
passing through the grating into bulk radiation the surface behind the
grating was partially covered with a thin dielectric overlayer (a
5$~\mu$m thick film of polystyrene). At this impedance step the SP was
transformed into bulk radiation  and interfered with the bulk radiation
diffracted at the aperture. The resulting interference pattern contains
information about the SP propagation parameters (real and imaginary
parts of SP wavevector) on the surface under study. The interference
pattern is registered by a detector moving in the direction
perpendicular to the sample surface. Joint processing of several
interferograms obtained at a fixed frequency for various values of
distances covered by the SP allows to determine the SP wavevector in the
given direction of the SP propagation. The SP damping was obtained from
the modulation depth of the interference pattern .

## Results

We have measured the dependence of the propagation parameters of SP on
the rotational angle of the grating near the SP band gap. The dispersion
curve of SP in the IR region differs from the light line by only
10$^{-4}$ – 10$^{-3}$. Figure 6 displays the experimentally measured
difference between the real part of SP wavevector and the wavevector of
bulk radiation with the wavelength 974.62 cm$^{-1}$ as a function of
the rotational angle of the grating. It is seen that
Re$(k_{\mathrm sp})$ exceeds the wavevector of bulk radiation
${\omega/ c}$ and changes from 1.00012${\omega/ c}$ to
1.0019${\omega/ c}$. The resonance angle of Bragg scattering at the
given frequency for the grating under study is indicated by the arrow in
Figure 6 and Figure 7. The wavevector of SP propagating on the smooth
part of the same silver covered silicon surface was measured also by the
method of interferometric phase SP spectroscopy: Re$(k_{\mathrm
sp})=(1.000102 \pm$<!-- -->0.000005)$\omega/c$ (at the wavelength
974.62 cm$^{-1}$).

The imaginary part of the SP wavevector determines the SP damping along
the surface. The modulation depth of the measured interference pattern
depends on the SP damping on the grating and was used to obtain the
angular dependence of the imaginary part of the SP wavevector presented
in Figure 7. So we examined experimentally how the real and imaginary
parts of the SP wavevector near the band gap depend on the relative
orientation of the Bragg vector of the metal grating and the propagation
direction of the SP.

# Discussion

In Figures 6 and 7 our experimental data and theoretical curves
calculated from the theory of Mills
(equation (<a href="#8" data-reference-type="ref" data-reference="8">[8]</a>))
– the dashed line, and from coupled-mode approach involving three
interacting modes (Section <a href="#three" data-reference-type="ref"
data-reference="three">2.3</a>) – the solid line, are presented
together. From these figures one can see that there is no satisfactory
agreement between Mills’ theory and the experimental points in the case
under consideration, whereas the presented theory gives much better fit
to the experimental data. However, from Figure 7 one can see that
coupled-mode theory with three interacting mode gives a slightly higher
values for imaginary part of the SP wavevector in the center of the gap.
In order to provide a more satisfactory description the model with four
interacting modes – the dotted line, was used (see Appendix for
details). From the results presented it can be seen that for gratings
usually used in the IR region (where the
inequality (<a href="#3" data-reference-type="ref" data-reference="3">[3]</a>),
as a rule, is always fulfilled) Mills’ theory is inadequate, while the
suggested theory is in good agreement with experiment.

The presented theoretical results can be represented not only in the
form $k_{\mathrm sp}=k_{\mathrm sp}(\varphi)$ at some fixed
$\omega$, but also as $k_{\mathrm sp}=k_{\mathrm sp}(\omega)$ at
some fixed $\varphi$ (i.e. as the ordinary SP dispersion curve). To
illustrate it in Figure 8 we present the theoretically calculated
dispersion curves of SP on the grating at $\varphi=0$ with different
corrugation amplitudes $h$ of the first Fourier component. The dashed
line shows the SP dispersion curve for a flat metal surface, the dotted
line that for $h=0.045~\mu$m and the solid line that for
$h=0.113~\mu$m (it is the real amplitude of our grating). From this
figure one can see that the shape of dispersion curve is perceptibly
asymmetric (i.e. differs from the Mills theory) even at $h=0.045~\mu$m
(i.e at $\lambda/235$).

# Conclusions

To describe the dispersion relation of SPs near photonic bandgaps a
physical model based on the coupled-mode approach involving three
interacting modes (two of them are SP modes and one is a light mode) is
developed. The inclusion of the interaction between bulk light waves and
SPs allows us to describe the experimental data for SP propagation
parameters in the IR region. The analytical expression obtained on the
basis of this model can be used to determine the form of the SP
dispersion curve near band gap for gratings with actually used
amplitudes of corrugations.
Equations (<a href="#48" data-reference-type="ref" data-reference="48">[48]</a>)
and (<a href="#45" data-reference-type="ref" data-reference="45">[45]</a>)
(using (<a href="#22a" data-reference-type="ref" data-reference="22a">[22a]</a>),
(<a href="#50a" data-reference-type="ref" data-reference="50a">[50a]</a>)
and (<a href="#56" data-reference-type="ref" data-reference="56">[56]</a>))
can be used for estimations of the dispersion curve splitting and SP
damping on gratings in the case under consideration.

#  Acknowledgments

Authors thank Prof. Yakovlev V.A. for useful discussions and
Dr. Zhukov A.A. for supplying at our disposal the silicon gratings. The
present research was supported by the Russian Foundation for Basic
Research and by program "Fundamental spectroscopy" of Russian Ministry
of Science.

## Appendix: four interacting modes

The fourth mode — the light one $d$ — is symmetrical to the light mode
$c$ with respect to the coordinate axis $y$, (see Figure 2) i.e. it
propagates to the right at a near-grazing angle to the surface. In an
"idealized" case (if an initial source excites only one undamped plasmon
wave $b$ on an infinite grating) the intensity of this mode is
vanishingly small. But in real experiments the intensity of the light
mode $d$ is nonzero, and this mode should be taken into account to
describe some fine effects, especially in recording the transmission and
reflection of SPs on Bragg gratings .

If all interactions between all four modes are taken into account, the
corresponding system of differential equations cannot be solved in an
analytical form. For this reason we take into account only the main (the
strongest) interactions between four modes. In this case we have four
interacting modes $a$, $b$, $c$, $d$:
``` math
\begin{eqnarray}
a&=&Ae^{i(\omega
t+k_a z)} \; , \nonumber\\ b&=&Be^{i(\omega t-k_b z)}
\; ,
\\
c&=&Ce^{i(\omega t+k_c z)}
\; ,
\nonumber\\
d&=&De^{i(\omega t-k_d z)}
\nonumber
\; ,
\end{eqnarray}
```
which obey relations of the type
``` math
\begin{eqnarray}
\left\{
       \begin{array}{lcl}
\frac{
\displaystyle
{\mathrm d}A}{
\displaystyle
{\mathrm d}z}&=&
\kappa_{ab}Be^{-i\Delta_{ba} z}
+\kappa_{ad}De^{-i\Delta_{da} z}
\\[\medskipamount]
{
\displaystyle
{\mathrm d}B\over
\displaystyle
{\mathrm d}z}&=&
\kappa_{ba}Ae^{i\Delta_{ba} z}
+\kappa_{bc}Ce^{i\Delta_{bc} z}
\\[\medskipamount]
{
\displaystyle
{\mathrm d}C\over
\displaystyle
{\mathrm d}z}&=&
\kappa_{cb}Be^{-i\Delta_{bc} z}
\\[\medskipamount]
{
\displaystyle
{\mathrm d}D\over
\displaystyle
{\mathrm d}z}&=&
\kappa_{da}Ae^{i\Delta_{da} z}
\; ,
         \end{array}
\right.

\end{eqnarray}
```
where $\Delta_{ba}$, $\Delta_{bc}$ and $\Delta_{da}$ are the
phase-mismatch constants between the corresponding modes.

From the equation for the conservation of the total power one can obtain
``` math
\begin{equation}
\kappa_{ab}\simeq
\kappa_{ba}^* \; ; \quad
\kappa_{cb}\simeq
\kappa_{bc}^*
\; ; \quad
\kappa_{ad}\simeq
\kappa_{da}^*
\; .

\end{equation}
```

Reducing
system (<a href="#a2" data-reference-type="ref" data-reference="a2">[a2]</a>)
to one differential equation, we remove the exponential dependence on
$z$ and obtain a differential equation of the fourth order with
constant coefficients. Seeking the solution in the form
$B(z)={\mathrm const}\cdot e^{\mu z}$ we obtain an equation for the
eigenvalues $\mu$ of the fourth order
``` math
\begin{equation}
\mu^4+\alpha\mu^3+\beta\mu^2+\gamma\mu+\epsilon =0
\; ,

\end{equation}
```
where the coefficients are
``` math
\begin{eqnarray}
\alpha
&=&
-i(\Delta_{bc}+2\Delta_{ba}-\Delta_{da})
\; ,
\nonumber\\
\beta
&=&
-|\kappa_{ba}|^2
-|\kappa_{da}|^2
-|\kappa_{bc}|^2
-2\Delta_{ba}\Delta_{bc}
+\Delta_{da}\Delta_{ba}
+\Delta_{da}\Delta_{bc}
-\Delta_{ba}^2
\; ,
  \\
\gamma
&=&
i\left(
|\kappa_{ba}|^2(\Delta_{bc} +
\Delta_{ba}
-\Delta_{da})
+|\kappa_{bc}|^2(2\Delta_{ba}
-\Delta_{da})
+|\kappa_{da}|^2\Delta_{bc}
+\Delta_{ba}\Delta_{bc}
(\Delta_{ba}-\Delta_{da})
\right)
\; ,
\nonumber\\
\epsilon
&=&
(|\kappa_{bc}|^2\Delta_{ba}
+|\kappa_{ba}|^2\Delta_{bc})(\Delta_{ba}-\Delta_{da})
+|\kappa_{da}|^2|\kappa_{bc}|^2
\nonumber
\; .
\end{eqnarray}
```

The general solution
of (<a href="#a2" data-reference-type="ref" data-reference="a2">[a2]</a>)
has the form (excepting the cases when $\mu_i=\mu_j$)
``` math
\begin{equation}
B(z)=b_1e^{\mu_1 z} + b_2e^{\mu_2 z}
+ b_3e^{\mu_3 z}
+ b_4e^{\mu_4 z}
\; ,

\end{equation}
```
where $\mu_{1,2,3,4}$ are solutions of
equation (<a href="#a4" data-reference-type="ref" data-reference="a4">[a4]</a>)
``` math
\begin{eqnarray}
\mu_{{{1\atop 2}\atop 3}\atop 4}
&=&
{
 {{{+\atop +}\atop -}\atop -}(2\xi)^{3/2}
 {{{+\atop -}\atop +}\atop -}2\sqrt{-2\xi^3-2p\xi^2
 {{{-\atop -}\atop +}\atop +}\sqrt{2}\xi^{3/2}q}
\over 4\xi}
 -{1\over 4}\alpha
\; ,

\end{eqnarray}
```
with
``` math
\begin{eqnarray}
\xi
&=&
S^{1/3}
-{R\over
3S^{1/3}}
-{p\over 3}
\; ,
\nonumber\\
S
&=&
-{1\over2}Q+{1\over 18}\sqrt{12R^3+81Q^2}
\; ,
\nonumber\\
Q
&=&
-{1\over 108}p^3-{1\over 8}q^2+{1\over 3}pr
\; ,
\nonumber\\
R
&=&
-{1\over 12}p^2-r
\; ,
\\
p
&=&
\beta
-{3\over 8}\alpha^2
\; ,
\nonumber\\
q
&=&
\gamma-{1\over 2} \beta \alpha
+{1\over 8}\alpha^3
\; ,
\nonumber\\
r
&=&
\epsilon-{1\over 4}\gamma\alpha+{1\over 16} \beta \alpha^2
-{3\over 256}\alpha^4
\nonumber
\; .
\end{eqnarray}
```

The constants $b_1$, $b_2$, $b_3$ and $b_4$ are determined from
the boundary conditions. Equations for these constants are cumbersome
and are not written here.

The wavevector $k'_b$ of the perturbed mode receives an addition
$\delta k_b$ to the wavevector $k_b$ of the unperturbed mode $b$ –
see
equation (<a href="#20b" data-reference-type="ref" data-reference="20b">[20b]</a>),
where $\delta k_b$, in general case, is given by
``` math
\begin{equation}
\delta k_b =-{ {\mathrm
d}\over {\mathrm d} z } \left[ \mbox{argument} \left(
b_1e^{\mu_1 z}
+ b_2e^{\mu_2 z}
+ b_3e^{\mu_3 z}
+ b_4e^{\mu_4 z}
\right)
\right]
\; ,

\end{equation}
```
or, what is the same
``` math
\begin{equation}
\delta k_b =i{{ {\mathrm
d}\over {\mathrm d} z }
\left(
b_1e^{\mu_1 z}
+ b_2e^{\mu_2 z}
+ b_3e^{\mu_3 z}
+ b_4e^{\mu_4 z}
\right)
\over
\left(
b_1e^{\mu_1 z}
+ b_2e^{\mu_2 z}
+ b_3e^{\mu_3 z}
+ b_4e^{\mu_4 z}
\right)
}
\; .

\end{equation}
```

For the new quantities in the presented equations the obvious relations
$|\kappa_{da}|\simeq|\kappa_{bc}|$ and $\Delta_{da}=-\Delta_{bc}$
can be used. The solution obtained has been used to derive the imaginary
part of the SP wavevector (the dotted line on Figure 7).

## References

[1] Yablonovitch, E., 1993, *J. opt. Soc. Am.* B, **10**, 283.

[2] Barnes, W.L., Kitson, S.C., Preist, T.W., and Sambles, J.R., 1997, *J.
opt. Soc. Am.* A, **14**, 1654.

[3] Raether, H., 1977, *Physics of thin films*, (New York: Acad.press),
**9**, 145.

[4] Barnes, W.L., Preist, T.W., Kitson, S.C., and Sambles, J.R., 1996,
*Phys. Rev.* B, **54**, 6227.

[5] Mills, D.L., 1977, *Phys. Rev.* B, **15**, 3097.

[6] Maradudin, A.A., 1982, *Surface polaritons*, edited by V.M. Agranovich
and D.L. Mills (Amsterdam: North-Holland).

[7] Konopsky, V.N., 2000, *Opt. Laser Technol.,* **32**, 15.

[8] Yariv, A., 1973, *IEEE J. quantum Electron.,* **QE9**, 919.

[9] Kogelnik, H., 1975, *Integrated Optics,* edited by T. Tamir (Berlin:
Springer-Verlag).

[10] Kogelnik, H., and Shank, C.V., 1972, *J. Appl. Phys.,* **43**, 2327.

[11] Gandelman, G.M., and Kondratenko, P.S., 1983, *Sov. Phys. – JETP
Letters*, **38**, 291. \[1983 *Pis’ma v ZhETF*, **38**, 246 (in
Russian)$$.

[12] Senior, T.B.A., 1960, *Appl. Sci. Res.* Sec.B, **8**, 418.

[13] http://www.ntmdt.ru

[14] Kretschmann, E., and Kröger, E., 1976, *Phys. stat. sol. (b),* **76**,
515.

[15] Alieva, E.V., Beitel, G., Kuzik, L.A., Sigarev, A.A., Yakovlev, V.A.,
Zhizhin, G.N., van der Meer, A.F.G., and van der Wiel, M.J., 1997,
*Appl. Spectroscopy,* **51**, 584.

[16] Schlesinger, Z., and Sievers, A.J., 1980, *Appl. Phys. Lett.,* **36**,
409.

[17] Alieva, E.V., and Konopsky, V.N. (to be published).

[18]