10
Techniques for exciting surface plasmons

Figures in this chapter that were generated by Mathematica routines contained in the Code section are labelled "Mathematica simulation" in the figure caption.  The reader can regenerate all of these graphs automatically by selecting "Evaluate Notebook" under the "Evaluation" menu bar in Mathematica.  If Mathematica asks "Do you want to automatically evaluate all the initialization cells in the notebook?" it is sufficient to answer "No."  If the reader is interested in generating similar graphs for different parameter ranges or different materials, it is generally straightforward to do this as well.  The Code section should first be expanded by double clicking on the down arrow of the rightmost bracket.  This will expand the notebook to display "Setup" and "Code for generating figures" subsubsections.  If the notebook has not yet been evaluated, the "Setup" subsubsection should be executed first by selecting (clicking once on) its rightmost bracket and entering "Shift-Return" on a standard keyboard or just "Enter" on an extended keyboard.  Then the "Code for generating figures" subsection should be expanded by double clicking on the downward arrow on the rightmost bracket of that subsubsection.  This subsubsection can be further expanded, if necessary, by again double clicking on the right bracket with the downward pointing arrow.  After expanding the subsubsections, eventually a list of figure captions is obtained.  By double-clicking on the up arrow of the figure caption, a text section and a code section for that figure are displayed.  The parameters for generating the graph are described in the text section.  The reader simply needs to modify these parameters in the code section, which is right below the text section, and then replot the figure by selecting (clicking once on) the rightmost bracket that encloses the code section for the figure and hitting "Shift-Return" on the keyboard or just "Enter" on an extended keyboard.  The figure should then be replotted using the new parameters.  The third line in the Setup subsection of the Code section is normally "option := optC" for generating color figures, but the line can be changed to "option := optB" for generating black and white figures instead if desired.
This notebook also contains a section with all of the figures in the text, a section describing the Chandezon vector diffraction technique, a section with the exercises and code for generating relevant figures, and a section with all the chapter references.  Each section can, as usual, be opened by double clicking on the down arrow of the right bracket.

Code

Many of the figures in this chapter were generated by Mathematica routines contained in this Code section.  The reader can regenerate these graphs automatically by selecting "Evaluate Notebook" under the "Evaluation" menu bar in Mathematica.  If the reader is interested in generating similar graphs for different parameter ranges or different materials, it is generally straightforward to do this as well.  If the notebook has not yet been evaluated, the "Setup" subsection below should be executed first by selecting the rightmost bracket of that section and entering "Shift-Return" on a standard keyboard or just "Enter" on an extended keyboard.  Then the figure in the chapter which is to be replotted should be located in the "Figures" subsection below.  The subsection can be expanded to show the figures by double clicking on the right bracket with the downward pointing arrow.  The parameters for generating the graph are described before each figure.  The reader simply needs to modify the parameters the parameters and replot the figure by selecting the rightmost bracket that encloses the input code for the figure and hitting "Shift-Return" on the keyboard or just "Enter" on an extended keyboard.  The figure should then be replotted with the new parameters.  For color figures, edit the third line in the "Setup" subsection to read "option := optC"

Setup

ClearAll["Global`*"];
Needs["PlotLegends`"];
option :=optC;   (*  change this option to "optC" for color plots or "optB" for black and white  *)
base="C:\Users\Bill\Desktop\Surface Plasmons\Final figures\Ch 10 launching\jpg";

Units

yocto=(10.^(-24));
zepto=(10.^(-21));
atto=(10.^(-18));
femto=(10.^(-15));
pico=(10.^(-12));
nano=(10.^(-9));
micro=(10.^(-6));
milli=(10.^(-3));
kilo=(10.^3);
mega=(10.^6);
giga=(10.^9);
tera=(10.^12);
peta=(10.^15);
exa=(10.^18);
zetta=(10.^21);
yotta=(10.^24);

Physical constants

NSub={
cc->2.99792458 (10^8) ,
e->1.6021773 (10^(-19)),
h->6.6260755 (10^(-34)),
ℏ->1.05457 (10^(-34)),
Kb->1.380658 (10^(-23)),
me->9.109389 (10^(-31)),
μ0->4 π (10^(-7)),
ϵ0->8.85419 (10^(-12)),
Lorentz->24.6974 nano
};

Plotting Examples

Graphics Styles

fs=12;(*FontSize for LabelStyle*)
is={350,200};(*ImageSize for one figure*)
gr={400,150};(*ImageSize for GraphicsRow*)
gg={400,300};(*ImageSize for GraphicsGrid*)
gs=-30; (*GraphicsSpacing for two figures*)
ggs={-20,20};(*GraphicsSpacing for four figures*)
thinthick=Thickness[0.005];
thick=Thickness[0.007]; (*Thickness of lines*)

Graphics color option: optC

optC={
LabelStyle->{FontSize -> fs},
GridLines->Automatic,
Frame->True,
ImageSize->is,
PlotRange->All,
PlotStyle->{
{Red,Thick},
{Green,Dashing[{0.06,0.02}],Thick},{Blue,Dashing[{0.02,0.01}],Thick},
{Magenta,DotDashed,Thick},{Purple,Thick,Dashing[{0.04,0.01}]},{Orange,Thick,Dashing[{0.01,0.01}]}}};

Graphics black option: optB

optB={
GridLines -> Automatic,
Frame -> True,
ImageSize->is,
LabelStyle -> {FontSize -> fs},
PlotRange->{All,All},
PlotStyle->{
{Black,thick},
{Black,thick,Dashing[{0.06,0.02}]},
{Black,thick,Dashing[{0.02,0.01}]},
{Black,thick,DotDashed},
{Black, thick, Dashing[{0.04, 0.01}]},
{Black,thick,Dashing[{0.04, 0.03}]}}};

General two plots using optX

xmin=-π;xmax=π;
ymin=-1;ymax=1;
Fig01=Plot[
{Sin[x],Sin[2x],Sin[3x],Sin[4x]},
{x,xmin,xmax},
PlotRange->{ymin,ymax},
FrameLabel -> {{"x",""},{"sinx", ""}},
Evaluate[optC]]

xmin=-π;xmax=π;
ymin=-1;ymax=1;
Fig02=Plot[
{Sin[x],Sin[2x]},
{x,xmin,xmax},Filling -> Axis,
FillingStyle -> Opacity[0.2],
PlotRange->{ymin,ymax},
FrameLabel -> {{"x",""},{"sinx", ""}},
Evaluate[optC]]

GraphicsRow[{Fig01,Fig02},
ImageSize->gr,
Frame->False,
Spacings->gs]

General four plots using optX

GraphicsGrid[{
  {Fig01,Fig02},
  {Fig01,Fig02}},
 ImageSize -> gg,
 Frame -> False,
 Spacings -> gs]

Manipulate plots

amin=1;
amax=4;
Manipulate[
Plot[
 Sin[a θ],
{θ,0,6},
Frame->True,
FrameLabel->{"θ (deg)","sin(a θ)","Function",StringJoin["a = ",ToString[.01Ceiling[100.a]],""]},ImageSize->.75is,
PlotRange->{-1,1},
GridLines->Automatic,
PlotStyle->{
{Red,thinthick},
{Green,thinthick,Dashing[{0.06,0.02}]},{Blue,thinthick,Dashing[{0.02,0.01}]},
{Magenta,thinthick,DotDashed},
{Purple,thinthick,Dashing[{0.04,0.01}]}}],
{{a,amin,"Phase"},amin,amax}]

Tabulated metallic refractive indices

A note on refractive indices:  There are several sources of tabulated refractive indices for silver, gold and other metals in the literature and on the web.  The Handbook of Optical Constants of Solids edited by Palik [Academic, San Diego, 1998] is probably the most frequently cited source in recent years and is one of the more comprehensive sources.  For the noble metals, the data of Johnson and Christy [Phys. Rev. B, 6 (1972) 4370-4379] are also frequently cited, especially in literature prior to the publication of the Handbook.  Indices from both of these sources have been included for use in simulations in this chapter.  However, the experimental precision in the Johnson and Christy data is somewhat limited and the data in the Handbook are frequently stitched together from several sources over an extended wavelength range.  The result of this is that calculations of extinction coefficients and peak field intensities often exhibit anomalous peaks at certain wavelengths corresponding to small changes in the refractive index from the limited experimental precision or a data stitch.  In order to avoid these unphysical effects in the modeling, it is generally much better to use an analytical model of the refractive indices that is smooth and continuous and has no precision problems.  Therefore, we have also included the anayltical models developed by Rakic et al. [Appl. Opt. 37 (1998) 5271-5283] and generally recommend their use for all modeling.
	Refractive indices for gold are taken from D. W. Lynch and W. R. Hunter, "Gold" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 293-295.  Refractive indices for silver are taken from D. W. Lynch and W. R. Hunter, "Silver" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 355-357.  Refractive indices for copper are taken from D. W. Lynch and W. R. Hunter, "Copper" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 284-285.  Refractive indices for aluminum are taken from D. Y. Smith, E. Shiles and M. Inokuti, "The optical properties of metallic aluminum" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 397-400.  Refractive indices for iridium are taken from D. W. Lynch and W. R. Hunter, "Iridium" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 302-303.  Refractive indices for rhodium are taken from D. W. Lynch and W. R. Hunter, "Rhodium" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 348-349.

GoldIndices = {{0.2883,1.742,1.9},{0.2952,1.776,1.918},{0.3024,1.812,1.92},{0.310, 1.83, 1.916}, {0.3179, 1.84, 1.904}, {0.3263, 1.824, 1.878}, {0.3351, 1.798, 1.86}, {0.3444, 1.766, 1.846}, {0.3542, 1.74, 1.848}, {0.3647, 1.716, 1.862}, {0.3757, 1.696, 1.906}, {0.3875, 1.674, 1.936}, {0.400, 1.658, 1.956}, {0.4133, 1.636, 1.958}, {0.4275, 1.616, 1.94}, {0.4428, 1.562, 1.904}, {0.4592, 1.426, 1.846}, {0.4769, 1.242, 1.796}, {0.4959, 0.916, 1.84}, {0.5166, 0.608, 2.12}, {0.5391, 0.402, 2.54}, {0.5636, 0.306, 2.88}, {0.5904, 0.236, 2.96}, {0.6199, 0.194, 3.05}, {0.6526, 0.166, 3.15}, {0.6888, 0.16, 3.8}, {0.7293, 0.164, 4.35},{0.7749, 0.174, 4.86}, {0.8266, 0.188, 5.39}, {0.8856, 0.21, 5.88}, {0.9537, 0.236, 6.47}, {1.033, 0.272, 7.07}, {1.127, 0.312, 7.93}, {1.240, 0.372, 8.77}, {1.265, 0.389, 8.09}, {1.291, 0.403, 8.25}, {1.319, 0.419, 8.42}, {1.378, 0.454, 8.77}, {1.409, 0.473, 8.96}, {1.442, 0.493, 9.15}, {1.476, 0.515, 9.36}, {1.512, 0.537, 9.58}, {1.550, 0.559, 9.81}, {1.590, 0.583, 10.1}, {1.631, 0.609, 10.3}, {1.675, 0.636, 10.6}, {1.722, 0.665, 10.9}, {1.771, 0.696, 11.2}, {1.823, 0.73, 11.5}, {1.879, 0.767, 11.9}, {1.937, 0.807, 12.2}, {2.000, 0.85, 12.6}, {2.066, 0.896, 13}, {2.138, 0.947, 13.4}, {2.214, 1.002, 13.9}, {2.296, 1.063, 14.4}, {2.384, 1.13, 14.9}, {2.480, 1.205, 15.5}, {2.583, 1.287, 16.1}, {2.695, 1.379, 16.8}, {2.818, 1.482, 17.5}, {2.952, 1.598, 18.3}};

SilverIndices = {{0.200,1.072,1.24},{0.20327869,1.098,1.26},{0.20666667,1.125,1.27},{0.2137931,1.173,1.29},{0.22142857,1.208,1.3},{0.22962963,1.238,1.31},{0.23846154,1.265,1.33},{0.248,1.298,1.35},{0.25306122,1.32,1.35},{0.25833333,1.343,1.35},{0.26382979,1.372,1.35},{0.26956522,1.404,1.33},{0.27555556,1.441,1.31},{0.28181818,1.476,1.26},{0.28837209,1.502,1.19},{0.2952381,1.519,1.08},{0.29879518,1.522,0.992},{0.30243902,1.496,0.882},{0.30617284,1.432,0.766},{0.310,1.323,0.647},{0.31155779,1.246,0.586},{0.31392405,1.149,0.54},{0.31552163,1.044,0.514},{0.31794872,0.932,0.504},{0.31958763,0.815,0.526},{0.32207792,0.708,0.565},{0.32375979,0.616,0.609},{0.32631579,0.526,0.663},{0.33066667,0.371,0.813},{0.33243968,0.321,0.902},{0.33513514,0.294,0.986},{0.33972603,0.259,1.12},{0.34444444,0.238,1.24},{0.35428571,0.209,1.44},{0.36470588,0.186,1.61},{0.37575758,0.2,1.67},{0.3875,0.192,1.81},{0.400,0.173,1.95},{0.41333333,0.173,2.11},{0.42758621,0.16,2.26},{0.44285714,0.157,2.4},{0.45925926,0.144,2.56},{0.47692308,0.132,2.72},{0.496,0.13,2.88},{0.51666667,0.13,3.07},{0.53913043,0.129,3.25},{0.56363636,0.12,3.45},{0.59047619,0.121,3.66},{0.620,0.131,3.88},{0.65263158,0.14,4.15},{0.68888889,0.14,4.44},{0.72941176,0.148,4.74},{0.775,0.143,5.09},{0.82666667,0.145,5.5},{0.88571429,0.163,5.95},{0.95384615,0.198,6.43},{1.03333333,0.226,6.99},{1.12727273,0.251,7.67},{1.240,0.329,8.49},{1.26530612,0.375,7.78},{1.29166667,0.383,7.92},{1.31914894,0.392,8.06},{1.34782609,0.401,8.21},{1.37777778,0.411,8.37},{1.40909091,0.421,8.37},{1.44186047,0.431,8.7},{1.47619048,0.442,8.88},{1.51219512,0.455,9.08},{1.550,0.469,9.32},{1.58974359,0.485,9.57},{1.63157895,0.501,9.84},{1.67567568,0.519,10.1},{1.72222222,0.537,10.4},{1.77142857,0.557,10.7},{1.87878788,0.6,11.4},{1.9375,0.624,11.8},{2.000,0.65,12.2}};

CopperIndices={{0.2952,1.42,1.64},{0.3099,1.34,1.72},{0.3263,1.34,1.81},{0.3444,1.31,1.87},{0.3646,1.27,1.95},{0.3874,1.18,2.21},{0.4133,1.18,2.21},{0.4428,1.17,2.36},{0.4768,1.15,2.5},{0.5166,1.12,2.6},{0.539,1.04,2.59},{0.5635,0.826,2.6},{0.5904,0.468,2.81},{0.6199,0.272,3.24},{0.6525,0.214,3.67},{0.6702,0.215,3.86},{0.688,0.213,4.05},{0.7084,0.214,4.24},{0.7293,0.223,4.43},{0.8265,0.26,5.26},{1.240,0.433,8.46},{1.265,0.496,6.78},{1.291,0.505,6.92},{1.319,0.515,7.06},{1.348,0.525,7.21},{1.378,0.536,7.36},{1.409,0.547,7.53},{1.442,0.559,7.7},{1.476,0.572,7.88},{1.512,0.586,8.06},{1.550,0.606,8.26},{1.589,0.627,8.47},{1.631,0.649,8.69},{1.675,0.672,8.92},{1.722,0.697,9.16},{1.771,0.723,9.41},{1.823,0.752,9.68},{1.878,0.782,9.97},{1.937,0.815,10.3},{2.000,0.85,10.6},{2.066,0.89,11},{2.138,0.933,11.3},{2.214,0.98,11.8},{2.296,1.03,12.2},{2.384,1.09,12.7},{2.480,1.15,13.2},{2.583,1.22,13.7}};

AluminumIndices={{0.300,0.276,3.61},{0.310,0.294,3.74},{0.3179,0.31,3.84},{0.3263,0.326,3.95},{0.3351,0.344,4.06},{0.3444,0.364,4.17},{0.350,0.375,4.24},{0.3542,0.385,4.3},{0.3647,0.407,4.43},{0.3757,0.432,4.56},{0.3875,0.46,4.71},{0.400,0.49,4.86},{0.4133,0.523,5.02},{0.4275,0.558,5.2},{0.4428,0.598,5.38},{0.450,0.618,5.47},{0.4592,0.644,5.58},{0.4769,0.695,5.8},{0.4959,0.755,6.03},{0.500,0.769,6.08},{0.5061,0.789,6.15},{0.5166,0.826,6.28},{0.5276,0.867,6.42},{0.5391,0.912,6.55},{0.550,0.958,6.69},{0.551,0.963,6.7},{0.5636,1.02,6.85},{0.5767,1.08,7},{0.5904,1.15,7.15},{0.600,1.2,7.26},{0.6048,1.22,7.31},{0.6199,1.3,7.48},{0.6358,1.39,7.65},{0.650,1.47,7.79},{0.6526,1.49,7.82},{0.6702,1.6,8.01},{0.6888,1.74,8.21},{0.700,1.83,8.31},{0.7085,1.91,8.39},{0.7293,2.14,8.57},{0.750,2.4,8.62},{0.7514,2.41,8.62},{0.7749,2.63,8.6},{0.7999,2.8,8.45},{0.825,2.75,8.31},{0.8266,2.74,8.31},{0.850,2.61,8.22},{0.8551,2.58,8.21},{0.875,2.38,8.18},{0.8856,2.24,8.21},{0.900,2.06,8.3},{0.9184,1.86,8.44},{0.925,1.77,8.49},{0.950,1.49,8.88},{0.9537,1.47,8.95},{0.9919,1.37,9.49},{1.000,1.35,9.58},{1.033,1.26,10},{1.078,1.21,10.6},{1.127,1.2,11.2},{1.181,1.21,11.8},{1.200,1.21,12},{1.240,1.21,12.5},{1.305,1.23,13.2},{1.378,1.26,14},{1.459,1.33,14.9},{1.500,1.38,15.4},{1.550,1.44,16},{1.653,1.59,17.1},{1.771,1.77,18.3},{1.907,1.99,19.8},{2.000,2.15,20.7}};

IridiumIndices={{0.2,0.94,1.68},{0.2033,0.94,1.73},{0.2066,0.95,1.78},{0.2101,0.97,1.82},{0.2138,0.96,1.86},{0.2175,0.97,1.92},{0.2214,0.98,1.98},{0.2254,0.99,2.03},{0.2296,1.0,2.09},{0.2339,1.02,2.15},{0.2384,1.04,2.22},{0.2431,1.07,2.29},{0.248,1.1,2.35},{0.253,1.13,2.42},{0.2583,1.18,2.49},{0.2638,1.24,2.56},{0.2695,1.31,2.6},{0.2755,1.37,2.65},{0.2818,1.45,2.68},{0.2883,1.51,2.7},{0.2952,1.58,2.71},{0.3024,1.62,2.7},{0.31,1.64,2.68},{0.3179,1.64,2.67},{0.3263,1.61,2.69},{0.3351,1.57,2.72},{0.3444,1.52,2.81},{0.3542,1.5,2.93},{0.3647,1.53,3.05},{0.3757,1.57,3.15},{0.3875,1.62,3.26},{0.4,1.68,3.35},{0.4133,1.73,3.43},{0.4275,1.77,3.51},{0.4428,1.81,3.61},{0.4592,1.85,3.73},{0.4769,1.91,3.86},{0.4959,1.98,4.00},{0.5166,2.07,4.14},{0.5391,2.18,4.26},{0.5636,2.29,4.38},{0.5904,2.4,4.48},{0.6199,2.5,4.57},{0.6526,2.57,4.68},{0.6888,2.64,4.81},{0.7293,2.69,4.92},{0.7749,2.68,5.08},{0.8266,2.65,5.39},{0.8856,2.72,5.74},{0.9537,2.85,6.07},{1.033,2.96,6.41},{1.127,3.04,6.84},{1.24,3.15,7.31},{1.378,3.19,7.88},{1.55,3.14,8.61},{1.771,2.93,9.78},{2.066,2.79,11.6}};

RhodiumIndices={{0.2000,0.78,1.85},{0.2066,0.76,1.93},{0.2138,0.79,2.00},{0.2214,0.8,2.06},{0.2296,0.8,2.14},{0.2384,0.79,2.23},{0.248,0.79,2.34},{0.253,0.79,2.39},{0.2583,0.79,2.46},{0.2638,0.79,2.52},{0.2695,0.78,2.6},{0.2755,0.79,2.68},{0.2818,0.8,2.76},{0.2883,0.83,2.85},{0.2952,0.84,2.94},{0.3024,0.84,3.03},{0.31,0.86,3.12},{0.3179,0.88,3.23},{0.3263,0.91,3.34},{0.3351,0.95,3.45},{0.3444,0.99,3.58},{0.3542,1.04,3.71},{0.3647,1.11,3.84},{0.3757,1.2,3.97},{0.3875,1.3,4.09},{0.4000,1.41,4.2},{0.4133,1.53,4.29},{0.4275,1.63,4.36},{0.4592,1.8,4.49},{0.4769,1.85,4.55},{0.4959,1.88,4.65},{0.5166,1.9,4.78},{0.5391,1.94,4.94},{0.5636,2.0,5.11},{0.5904,2.05,5.3},{0.6199,2.12,5.51},{0.6526,2.2,5.76},{0.6888,2.3,6.02},{0.7293,2.42,6.33},{0.7749,2.6,6.64},{0.8266,2.78,6.97},{0.8856,3.01,7.31},{0.9537,3.26,7.63},{0.9919,3.38,7.8},{1.033,3.51,7.94},{1.078,3.6,8.09},{1.127,3.67,8.26},{1.181,3.72,8.44},{1.240,3.71,8.67},{1.305,3.69,8.94},{1.378,3.62,9.36},{1.459,3.64,9.81},{1.550,3.63,10.3},{1.653,3.63,11.0},{1.771,3.67,11.7},{1.907,3.77,12.6},{2.066,3.87,13.5}};

GoldIndices2 and SilverIndices2 are from Johnson and Christy, "Optical constants of the noble metals," Phys. Rev. B, 6 (1972) 4370-4379.  Note that at short wavelengths there is a substantial difference between the refractive indices of Ag in this list and that from the Handbook of Optical Constants of Solids.  When comparing to examples in the literature it is important to use the same source for the optical constants!

GoldIndices2 = {{0.18788, 1.28, 1.188}, {0.19165, 1.32, 1.203}, {0.19528, 1.34, 1.226}, {0.19936, 1.33, 1.251}, {0.20328, 1.33, 1.277}, {0.20736, 1.3, 1.304}, {0.21197, 1.3, 1.35}, {0.2164, 1.3, 1.387}, {0.22143, 1.3, 1.427}, {0.22628, 1.31, 1.46}, {0.23134, 1.3, 1.497}, {0.23709, 1.32, 1.536}, {0.24266, 1.32, 1.577}, {0.2490, 1.33, 1.631}, {0.25514, 1.33, 1.688}, {0.2616, 1.35, 1.749}, {0.26898, 1.38, 1.803}, {0.27617, 1.43, 1.847}, {0.2844, 1.47, 1.869}, {0.29245, 1.49, 1.878}, {0.30097, 1.53, 1.889}, {0.31078, 1.53, 1.893}, {0.32041, 1.54, 1.898}, {0.33155, 1.48, 1.883}, {0.34254, 1.48, 1.871}, {0.35429, 1.5, 1.866}, {0.36795, 1.48, 1.895}, {0.38154, 1.46, 1.933}, {0.39744, 1.47, 1.952}, {0.41333, 1.46, 1.958}, {0.43056, 1.45, 1.948}, {0.45091, 1.38, 1.914}, {0.47148, 1.31, 1.849}, {0.4960, 1.04, 1.833}, {0.52101, 0.62, 2.081}, {0.54867, 0.43, 2.455}, {0.58216, 0.29, 2.863}, {0.61692, 0.21, 3.272}, {0.65957, 0.14, 3.697}, {0.70455, 0.13, 4.103}, {0.7561, 0.14, 4.542}, {0.82119, 0.16, 5.083}, {0.89209, 0.17, 5.663}, {0.98413, 0.22, 6.35}, {1.08772, 0.27, 7.15}, {1.21569, 0.35, 8.145}, {1.39326, 0.43, 9.519}, {1.61039, 0.56, 11.21}, {1.9375, 0.92, 13.78}};

SilverIndices2={{0.18788,1.07,1.212},{0.191.65,1.1,1.232},{0.19528,1.12,1.255},{0.19936,1.14,1.277},{0.20328,1.15,1.296},{0.20736,1.18,1.312},{0.21197,1.2,1.325},{0.21640,1.22,1.336},{0.22143,1.25,1.342},{0.22628,1.26,1.344},{0.23134,1.28,1.357},{0.23709,1.28,1.367},{0.24266,1.3,1.378},{0.24900,1.31,1.389},{0.25514,1.33,1.393},{0.26160,1.35,1.387},{0.26898,1.38,1.372},{0.27617,1.41,1.331},{0.28440,1.41,1.264},{0.29245,1.39,1.161},{0.30097,1.34,0.964},{0.31078,1.13,0.616},{0.32041,0.81,0.392},{0.33155,0.17,0.829},{0.34254,0.14,1.142},{0.35429,0.1,1.419},{0.36795,0.07,1.657},{0.38154,0.05,1.864},{0.39744,0.05,2.07},{0.41333,0.05,2.275},{0.43056,0.04,2.462},{0.45091,0.04,2.657},{0.47148,0.05,2.869},{0.49600,0.05,3.093},{0.52101,0.05,3.324},{0.54867,0.06,3.586},{0.58216,0.05,3.858},{0.61692,0.06,4.152},{0.65957,0.05,4.483},{0.70455,0.04,4.838},{0.75610,0.03,5.242},{0.82119,0.04,5.727},{0.89209,0.04,6.312},{0.98413,0.04,6.992},{1.08772,0.04,7.795},{1.21569,0.09,8.828},{1.39326,0.13,10.1},{1.61039,0.15,11.85},{1.9375,0.24,14.08}};

IronIndices2 and CobaltIndices2 are from Johnson and Christy, "Optical constants of transition metals: Ti, V, Cr, Mn, Fe, Co, Ni, and Pd," Phys. Rev. B, 9 (1974) 5056 - 5070.

IronIndices2 = {{0.18788,1.29,1.35},{0.19165,1.35,1.37},{0.19528,1.42,1.39},{0.19936,1.45,1.4},{0.20328,1.47,1.4},{0.20736,1.49,1.41},{0.21197,1.47,1.43},{0.2164,1.47,1.44},{0.22143,1.47,1.47},{0.22628,1.47,1.49},{0.23134,1.48,1.53},{0.23709,1.48,1.57},{0.24266,1.5,1.61},{0.249,1.51,1.66},{0.25514,1.53,1.7},{0.2616,1.56,1.75},{0.26898,1.59,1.79},{0.27617,1.62,1.84},{0.2844,1.64,1.88},{0.29245,1.65,1.94},{0.30097,1.67,2},{0.31078,1.69,2.06},{0.32041,1.74,2.12},{0.33155,1.78,2.19},{0.34254,1.85,2.27},{0.35429,1.93,2.35},{0.36795,2.02,2.43},{0.38154,2.12,2.5},{0.39744,2.24,2.58},{0.41333,2.35,2.65},{0.43056,2.48,2.71},{0.45091,2.59,2.77},{0.47148,2.67,2.82},{0.496,2.74,2.88},{0.52101,2.86,2.91},{0.54867,2.95,2.93},{0.58216,2.94,2.99},{0.61692,2.88,3.05},{0.65957,2.92,3.1},{0.70455,2.86,3.19},{0.7561,2.87,3.28},{0.82119,2.94,3.39},{0.89209,2.96,3.56},{0.98413,2.92,3.79},{1.08772,2.97,4.06},{1.21569,3.03,4.39},{1.39326,3.09,4.83},{1.61039,3.11,5.39},{1.9375,3.17,6.12}};

CobaltIndices2={{0.18788,1.16,1.59},{0.19165,1.21,1.63},{0.19528,1.26,1.67},{0.19936,1.29,1.71},{0.20328,1.32,1.75},{0.20736,1.36,1.78},{0.21197,1.38,1.82},{0.2164,1.41,1.84},{0.22143,1.43,1.85},{0.22628,1.45,1.86},{0.23134,1.47,1.87},{0.23709,1.47,1.89},{0.24266,1.46,1.91},{0.249,1.45,1.93},{0.25514,1.45,1.97},{0.2616,1.44,2.01},{0.26898,1.44,2.04},{0.27617,1.44,2.09},{0.2844,1.44,2.14},{0.29245,1.44,2.19},{0.30097,1.44,2.25},{0.31078,1.44,2.31},{0.32041,1.45,2.37},{0.33155,1.46,2.44},{0.34254,1.48,2.52},{0.35429,1.49,2.61},{0.36795,1.5,2.71},{0.38154,1.53,2.82},{0.39744,1.57,2.93},{0.41333,1.61,3.05},{0.43056,1.67,3.17},{0.45091,1.74,3.28},{0.47148,1.81,3.41},{0.496,1.88,3.55},{0.52101,1.97,3.68},{0.54867,2.05,3.82},{0.58216,2.13,3.96},{0.61692,2.19,4.11},{0.65957,2.25,4.27},{0.70455,2.31,4.45},{0.7561,2.4,4.64},{0.82119,2.53,4.88},{0.89209,2.65,5.16},{0.98413,2.78,5.5},{1.08772,2.94,5.88},{1.21569,3.17,6.31},{1.39326,3.42,6.77},{1.61039,3.61,7.26},{1.9375,3.87,7.79}};

Routines to compute refractive indices of various metals

Drude[λ_]:=Module[{n,ω,ϵ,c,h,ωp,γ,ϵ∞},
c=2.99792456×(10^14);          (*  speed of light in μm/s  *)
h=4.1359×(10^(-15));                (*  Plank's constant in eV s  *)
ϵ∞=4;                                     (*  background dielectric constant  *)
ωp=15;                                       (*  ωp in eV  *)
γ = 0.1;                               (*  damping coefficient in eV  *)
ω = 7.7905/λ;                     (*  ang. frequency in eV when λ is in μm  *)
ϵ=ϵ∞-(ωp^2)/(ω^2 + I γ ω);
n = (Sqrt[ϵ]);
(*  Print["λ=",λ," n=",n];  *)
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeAg[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.065,0.124,0.011,0.840,5.646};
γ={3.886,0.452,0.065,0.916,2.419};
ωω={0.816,4.481,8.185,9.083,20.29};
ωp=9.01;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.845(ωp^2)/(ω^2+I 0.048ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,5}];
n = (Sqrt[ϵ]);
Return[n]];


(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeAu[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.024,0.010,0.071,0.601,4.384};
γ={0.241,0.345,0.870,2.494,2.214};
ωω={0.415,0.830,2.969,4.304,13.32};
ωp=9.03;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.760(ωp^2)/(ω^2+I 0.053ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,5}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeCu[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.061,0.104,0.723,0.638};
γ={0.378,1.056,3.213,4.305};
ωω={0.291,2.957,5.300,11.18};
ωp=10.83;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.575(ωp^2)/(ω^2+I 0.030ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeAl[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.227,0.050,0.166,0.030};
γ={0.333,0.312,1.351,3.382};
ωω={0.162,1.544,1.808,3.473};
ωp=14.98;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.523(ωp^2)/(ω^2+I 0.047ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeBe[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.031,0.140,0.530,0.130};
γ={1.664,3.395,4.454,1.802};
ωω={0.100,1.032,3.183,4.604};
ωp=18.51;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.084(ωp^2)/(ω^2+I 0.035ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeCr[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.151,0.150,1.149,0.825};
γ={3.175,1.305,2.676,1.335};
ωω={0.121,0.543,1.970,8.775};
ωp=10.75;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.168(ωp^2)/(ω^2+I 0.047ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeNi[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.100,0.135,0.106,0.729};
γ={4.511,1.334,2.178,6.292};
ωω={0.174,0.582,1.597,6.089};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.096(ωp^2)/(ω^2+I 0.048ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudePd[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.649,0.121,0.638,0.453};
γ={2.950,1.555,4.621,3.236};
ωω={0.336,0.501,1.659,5.715};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.330(ωp^2)/(ω^2+I 0.008ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudePt[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.191,0.659,0.547,3.576};
γ={0.517,1.838,3.668,8.517};
ωω={0.780,1.314,3.141,9.249};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.333(ωp^2)/(ω^2+I 0.080ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeTi[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.899,0.393,0.187,0.001};
γ={2.276,2.518,1.663,1.762};
ωω={0.777,1.545,2.509,19.43};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.148(ωp^2)/(ω^2+I 0.082ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
LorentzDrudeW[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.054,0.166,0.706,2.590};
γ={0.530,1.281,3.332,5.836};
ωω={1.004,1.917,3.580,7.498};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.206(ωp^2)/(ω^2+I 0.064ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(*  from D.W.Lynch and W.R.Hunter,"Gold" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 293-295  *)
Gold[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<GoldIndices[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>GoldIndices[[Length[GoldIndices],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=GoldIndices[[i,1]];
n2=GoldIndices[[i,2]];
K2=GoldIndices[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=GoldIndices[[i-1,1]];
n1=GoldIndices[[i-1,2]];
K1=GoldIndices[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[GoldIndices]}]];
Return[n]];

(*  from D.W.Lynch and W.R.Hunter "Iridium" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 302-303  *)
Iridium[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<IridiumIndices[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>IridiumIndices[[Length[IridiumIndices],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=IridiumIndices[[i,1]];
n2=IridiumIndices[[i,2]];
K2=IridiumIndices[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=IridiumIndices[[i-1,1]];
n1=IridiumIndices[[i-1,2]];
K1=IridiumIndices[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[IridiumIndices]}]];
Return[n]];

(*  from D.W.Lynch and W.R.Hunter "Rhodium" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 348-349  *)
Rhodium[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<RhodiumIndices[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>RhodiumIndices[[Length[RhodiumIndices],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=RhodiumIndices[[i,1]];
n2=RhodiumIndices[[i,2]];
K2=RhodiumIndices[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=RhodiumIndices[[i-1,1]];
n1=RhodiumIndices[[i-1,2]];
K1=RhodiumIndices[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[RhodiumIndices]}]];
Return[n]];

(*  from Johnson and Christy,"Optical constants of the noble metals,"Phys. Rev. B,6(1972) 4370-4379  *)
Gold2[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<GoldIndices2[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>GoldIndices2[[Length[GoldIndices2],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=GoldIndices2[[i,1]];
n2=GoldIndices2[[i,2]];
K2=GoldIndices2[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=GoldIndices2[[i-1,1]];
n1=GoldIndices2[[i-1,2]];
K1=GoldIndices2[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[GoldIndices2]}]];
Return[n]];

(*  from D.W.Lynch and W.R.Hunter "Silver" in Handbook of Optical Constants of Solids ed.E.D.Palik (Academic,San Diego,1998) 355-357  *)
Silver[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<SilverIndices[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>SilverIndices[[Length[SilverIndices],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=SilverIndices[[i,1]];
n2=SilverIndices[[i,2]];
K2=SilverIndices[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=SilverIndices[[i-1,1]];
n1=SilverIndices[[i-1,2]];
K1=SilverIndices[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[SilverIndices]}]];
Return[n]];

(*  from Johnson and Christy,"Optical constants of the noble metals,"Phys. Rev. B,6(1972) 4370-4379  *)
Silver2[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<SilverIndices2[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>SilverIndices2[[Length[SilverIndices2],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=SilverIndices2[[i,1]];
n2=SilverIndices2[[i,2]];
K2=SilverIndices2[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=SilverIndices2[[i-1,1]];
n1=SilverIndices2[[i-1,2]];
K1=SilverIndices2[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[SilverIndices2]}]];
Return[n]];

(*  from Johnson and Christy, Phys. Rev. B,9(1974) 5056-5070  *)
Iron2[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<IronIndices2[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>IronIndices2[[Length[IronIndices2],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=IronIndices2[[i,1]];
n2=IronIndices2[[i,2]];
K2=IronIndices2[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=IronIndices2[[i-1,1]];
n1=IronIndices2[[i-1,2]];
K1=IronIndices2[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[IronIndices2]}]];
Return[n]];

(*  from Johnson and Christy, Phys. Rev. B,9(1974) 5056-5070  *)
Cobalt2[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<CobaltIndices2[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>CobaltIndices2[[Length[CobaltIndices2],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=CobaltIndices2[[i,1]];
n2=CobaltIndices2[[i,2]];
K2=CobaltIndices2[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=CobaltIndices2[[i-1,1]];
n1=CobaltIndices2[[i-1,2]];
K1=CobaltIndices2[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[CobaltIndices2]}]];
Return[n]];

(*  from D.W.Lynch and W.R.Hunter "Copper" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 284-285  *)
Copper[λ_] := Module[{n, n1, n2, K1, K2, λ1, λ2},
   If[λ < CopperIndices[[1, 1]], Print[" Input wavelength is too short."]; Abort[]];
   If[λ > CopperIndices[[Length[CopperIndices], 1]], Print[" Input wavelength is too long."]; Abort[]];
   Catch[Do[λ2 = CopperIndices[[i, 1]];
     n2 = CopperIndices[[i, 2]];
     K2 = CopperIndices[[i, 3]];
     If[λ == λ2, n = n2 + I K2; Throw[n],
      If[λ < λ2, λ1 = CopperIndices[[i - 1, 1]];
       n1 = CopperIndices[[i - 1, 2]];
       K1 = CopperIndices[[i - 1, 3]];
       n = n1 + (n2 - n1) (λ - λ1)/(λ2 - λ1) + I (K1 + (K2 - K1) (λ - λ1)/(λ2 - λ1));
       Throw[n]]], {i, Length[CopperIndices]}]];
   Return[n]];

(*  from D.Y.Smith E.Shiles and M.Inokuti "The optical properties of metallic aluminum" in Handbook of Optical Constants of Solids ed.E.D.Palik (Academic,San Diego,1998) 397-400  *)
Aluminum[λ_] := Module[{n, n1, n2, K1, K2, λ1, λ2},
   If[λ < AluminumIndices[[1, 1]], Print[" Input wavelength is too short."]; Abort[]];
   If[λ > AluminumIndices[[Length[AluminumIndices], 1]], Print[" Input wavelength is too long."]; Abort[]];
   Catch[Do[λ2 = AluminumIndices[[i, 1]];
     n2 = AluminumIndices[[i, 2]];
     K2 = AluminumIndices[[i, 3]];
     If[λ == λ2, n = n2 + I K2; Throw[n],
      If[λ < λ2, λ1 = AluminumIndices[[i - 1, 1]];
       n1 = AluminumIndices[[i - 1, 2]];
       K1 = AluminumIndices[[i - 1, 3]];
       n = n1 + (n2 - n1) (λ - λ1)/(λ2 - λ1) + I (K1 + (K2 - K1) (λ - λ1)/(λ2 - λ1));
       Throw[n]]], {i, Length[AluminumIndices]}]];
   Return[n]];

Various thin film routines

The dielectric function is computed from the refractive index in this routine.

ϵ[λ_,n_]:=Module[{nmetal,ϵm},
nmetal=Which[n==-1,Gold[λ],n==-2,Silver[λ],n==-3,Copper[λ],n==-4,Aluminum[λ],n==-5,Drude[λ],n==-6,LorentzDrudeAu[λ],n==-7,LorentzDrudeAg[λ],n==-8,LorentzDrudeCu[λ],n==-9,LorentzDrudeAl[λ],n==-10,Gold2[λ],n==-11,Silver2[λ],n==-12,LorentzDrudeBe[λ],n==-13,LorentzDrudeCr[λ],n==-14,LorentzDrudeNi[λ],n==-15,LorentzDrudePd[λ],n==-16,LorentzDrudePt[λ],n==-17,LorentzDrudeTi[λ],n==-18,LorentzDrudeW[λ],n==-19,Iron2[λ],n==-20,Cobalt2[λ],n==-21,Rhodium[λ],n==-22,Iridium[λ],True,n];
ϵm=(nmetal^2);
Return[ϵm]];

The reflectivity in the Otto or Kretschmann configurations is calculated in these routines.  For the Otto configuration, the input parameters are the wavelength λ, the angle of incidence θ of the beam on the bottom surface of the prism, the index of the prism (Subscript[n, 0]), the index of the dielectric in the gap (Subscript[n, 1]), and the index of the metal (Subscript[n, 2]), and the thickness of the dielectric film d.  For the Kretschmann configuration the input parameters are the wavelength λ, the angle of incidence θ, the index of the incident medium (Subscript[n, 0]), the metal index (Subscript[n, 1]), the index of the transmitting medium (Subscript[n, 2]), and the thickness of the metal film d.  RTM computes the reflectivity for incident light polarized in the plane of incidence (TM or p polarization), while RTE computes the reflectivity for incident light polarized perpendicular to the plane of incidence (TE or s polarization).

RTM[λ_,θ_,n0_,n1_,n2_,d_]:=Module[{nn1,nn2,k0,cos0,cos1,cos2,sin1,sin2,a1,a2,b1,b2,exp1,refl,tran,R,T},
k0=2 π/λ;
nn1=Which[n1==-1,Gold[λ],n1==-2,Silver[λ],n1==-3,Copper[λ],n1==-4,Aluminum[λ],n1==-5,Drude[λ],n1==-6,LorentzDrudeAu[λ],n1==-7,LorentzDrudeAg[λ],n1==-8,LorentzDrudeCu[λ],n1==-9,LorentzDrudeAl[λ],n1==-10,Gold2[λ],n1==-11,Silver2[λ],n1==-12,LorentzDrudeBe[λ],n1==-13,LorentzDrudeCr[λ],n1==-14,LorentzDrudeNi[λ],n1==-15,LorentzDrudePd[λ],n1==-16,LorentzDrudePt[λ],n1==-17,LorentzDrudeTi[λ],n1==-18,LorentzDrudeW[λ],n1==-19,Iron2[λ],n1==-20,Cobalt2[λ],n1==-21,Rhodium[λ],n1==-22,Iridium[λ],True,n1];
nn2=Which[n2==-1,Gold[λ],n2==-2,Silver[λ],n2==-3,Copper[λ],n2==-4,Aluminum[λ],n2==-5,Drude[λ],n2==-6,LorentzDrudeAu[λ],n2==-7,LorentzDrudeAg[λ],n2==-8,LorentzDrudeCu[λ],n2==-9,LorentzDrudeAl[λ],n2==-10,Gold2[λ],n2==-11,Silver2[λ],n2==-12,LorentzDrudeBe[λ],n2==-13,LorentzDrudeCr[λ],n2==-14,LorentzDrudeNi[λ],n2==-15,LorentzDrudePd[λ],n2==-16,LorentzDrudePt[λ],n2==-17,LorentzDrudeTi[λ],n2==-18,LorentzDrudeW[λ],n2==-19,Iron2[λ],n2==-20,Cobalt2[λ],n2==-21,Rhodium[λ],n2==-22,Iridium[λ],True,n2];
cos0=Cos[2 π θ/360];
sin1=n0 Sin[2 π θ/360]/nn1;
cos1=Sqrt[1-(sin1^2)];
sin2=n0 Sin[2 π θ/360]/nn2;
cos2=Sqrt[1-(sin2^2)];
a1=(n0 cos1-nn1 cos0)/(nn2 cos1-nn1 cos2);
a2=(nn1 cos0+n0 cos1)/(nn2 cos1+nn1 cos2);
b1=(n0 cos1-nn1 cos0)/(nn2 cos1+nn1 cos2);
b2=(nn1 cos0+n0 cos1)/(nn2 cos1-nn1 cos2);
exp1=Exp[2 I nn1 k0 d cos1];
refl=(a1-a2 exp1)/(b1 exp1-b2);
R=(Abs[refl])^2;
tran=(a2+b1 refl)Exp[ I nn1 k0 d cos1];
T= nn2  Re[cos2] (Abs[tran])^2/(n0 cos0);
Return[R]];

RTE[λ_,θ_,n0_,n1_,n2_,d_]:=Module[{nn1,nn2,k0,cos0,cos1,cos2,sin1,sin2,a1,a2,b1,b2,exp1,refl,tran,R,T},
k0=2 π/λ;
nn1=Which[n1==-1,Gold[λ],n1==-2,Silver[λ],n1==-3,Copper[λ],n1==-4,Aluminum[λ],n1==-5,Drude[λ],n1==-6,LorentzDrudeAu[λ],n1==-7,LorentzDrudeAg[λ],n1==-8,LorentzDrudeCu[λ],n1==-9,LorentzDrudeAl[λ],n1==-10,Gold2[λ],n1==-11,Silver2[λ],n1==-12,LorentzDrudeBe[λ],n1==-13,LorentzDrudeCr[λ],n1==-14,LorentzDrudeNi[λ],n1==-15,LorentzDrudePd[λ],n1==-16,LorentzDrudePt[λ],n1==-17,LorentzDrudeTi[λ],n1==-18,LorentzDrudeW[λ],n1==-19,Iron2[λ],n1==-20,Cobalt2[λ],n1==-21,Rhodium[λ],n1==-22,Iridium[λ],True,n1];
nn2=Which[n2==-1,Gold[λ],n2==-2,Silver[λ],n2==-3,Copper[λ],n2==-4,Aluminum[λ],n2==-5,Drude[λ],n2==-6,LorentzDrudeAu[λ],n2==-7,LorentzDrudeAg[λ],n2==-8,LorentzDrudeCu[λ],n2==-9,LorentzDrudeAl[λ],n2==-10,Gold2[λ],n2==-11,Silver2[λ],n2==-12,LorentzDrudeBe[λ],n2==-13,LorentzDrudeCr[λ],n2==-14,LorentzDrudeNi[λ],n2==-15,LorentzDrudePd[λ],n2==-16,LorentzDrudePt[λ],n2==-17,LorentzDrudeTi[λ],n2==-18,LorentzDrudeW[λ],n2==-19,Iron2[λ],n2==-20,Cobalt2[λ],n2==-21,Rhodium[λ],n2==-22,Iridium[λ],True,n2];
cos0=Cos[2 π θ/360];
sin1=n0 Sin[2 π θ/360]/nn1;
cos1=Sqrt[1-(sin1^2)];
sin2=n0 Sin[2 π θ/360]/nn2;
cos2=Sqrt[1-(sin2^2)];
a1=(nn1 cos1-n0 cos0)/(nn1 cos1+nn2 cos2);
a2=(nn1 cos1+n0 cos0)/(nn1 cos1-nn2 cos2);
b1=(nn1 cos1-n0 cos0)/(nn1 cos1-nn2 cos2);
b2=(nn1 cos1+n0 cos0)/(nn1 cos1+nn2 cos2);
exp1=Exp[2 I nn1 k0 d cos1];
refl=(a1 exp1-a2 )/(b1-b2 exp1);
R=(Abs[refl])^2;
tran=(a1+b2 refl)Exp[ I nn1 k0 d cos1];
T= nn2  Re[cos2] (Abs[tran])^2/(n0 cos0);
Return[R]];

The field in the medium outside of the prism in the Kretschmann configuration is calculated in this routine.  The input parameters are the wavelength λ, the angle of incidence θ, the index of the incident medium (Subscript[n, 0]), the metal index (Subscript[n, 1]), the index of the transmitting medium (Subscript[n, 2]), the thickness of the metal film d, and the distance below the surface z.  The output is the list {|(Subscript[E, x])|,|(Subscript[E, z])|,|E|}.

field[λ_,θ_,n0_,n1_,n2_,d_,z_]:=Module[{nmetal,k0,cos0,cos1,cos2,sin0,sin1,sin2,a1,a2,b1,b2,c1,c2,exp1,refl,tran,R,T,Ex,Ez,E,E1,E2},
k0=2 π/λ;
nmetal=Which[n1==-1,Gold[λ],n1==-2,Silver[λ],n1==-3,Copper[λ],n1==-4,Aluminum[λ],n1==-5,Drude[λ],n1==-6,LorentzDrudeAu[λ],n1==-7,LorentzDrudeAg[λ],n1==-8,LorentzDrudeCu[λ],n1==-9,LorentzDrudeAl[λ],n1==-10,Gold2[λ],n1==-11,Silver2[λ],n1==-12,LorentzDrudeBe[λ],n1==-13,LorentzDrudeCr[λ],n1==-14,LorentzDrudeNi[λ],n1==-15,LorentzDrudePd[λ],n1==-16,LorentzDrudePt[λ],n1==-17,LorentzDrudeTi[λ],n1==-18,LorentzDrudeW[λ],n1==-19,Iron2[λ],n1==-20,Cobalt2[λ],n1==-21,Rhodium[λ],n1==-22,Iridium[λ],True,n1];
sin0=Sin[2 π θ/360];
cos0=Cos[2 π θ/360];
sin1=n0 Sin[2 π θ/360]/nmetal;
cos1=Sqrt[1-(sin1^2)];
sin2=n0 Sin[2 π θ/360]/n2;
cos2=Sqrt[1-(sin2^2)];
a1=(n0 cos1-nmetal cos0)/(n2 cos1-nmetal cos2);
a2=(nmetal cos0+n0 cos1)/(n2 cos1+nmetal cos2);
b1=(n0 cos1-nmetal cos0)/(n2 cos1+nmetal cos2);
b2=(nmetal cos0+n0 cos1)/(n2 cos1-nmetal cos2);
exp1=Exp[2 I nmetal k0 d cos1];
refl=(a1-a2 exp1)/(b1 exp1-b2);
R=(Abs[refl])^2;
tran=(a2+b1 refl)Exp[ I nmetal k0 d cos1];
T= n2  Re[cos2] (Abs[tran])^2/(n0 cos0);
If[z>0.0,
Ex=cos2 tran Exp[I n2 k0 cos2 z];
Ez=sin2 tran Exp[I n2 k0 cos2 z],
If[z>-d,
c1=(nmetal cos2+n2 cos1)/(2 nmetal cos1);
c2=(n2 cos1-nmetal cos2)/(2 nmetal cos1);
E1=c1 Exp[-I nmetal k0 d cos1] tran;
E2=c2 Exp[+I nmetal k0 d cos1] tran;
Ex=E1 cos1 Exp[+I nmetal k0 cos1 (z+d)]-E2 cos1 Exp[-I nmetal k0 cos1 (z+d)];
Ez=E1 sin1 Exp[+I nmetal k0 cos1 (z+d)]+E2 sin1 Exp[-I nmetal k0 cos1 (z+d)],
Ex=cos0 Exp[+I n0 k0 cos0 (z+d)]-refl cos0 Exp[-I n0 k0 cos0 (z+d)];
Ez=sin0 Exp[+I n0 k0 cos0 (z+d)]+refl sin0 Exp[-I n0 k0 cos0 (z+d)]]];
E=Sqrt[Abs[Ex]^2 + Abs[Ez]^2];
Return[{{z,Abs[Ex]},{z,Abs[Ez]},{z,E}}]];

Calculation of the complex SP wavevector for a Drude metal with damping assuming the SP frequency is real.

Ksp:=Ω Sqrt[(Ω^2 -1+I Γ Ω)/(2 Ω^2 -1 + 2 I Γ Ω)];

Diffraction grating routines

Resonance angle

Routine "resonance" computes the approximate polar resonance angle θ as a function of the wavelength λ, the azimuthal angle ϕ, the refractive index of the metal (Subscript[n, m]), the refractive index of the surrounding dielectric medium (Subscript[n, d]), the grating period d, and the order m.

((((Subscript[n, d]) sinθ))^2)+((2m (Subscript[n, d]) λ)/d) sinθ sinϕ+(((((m λ)/d)))^2)=(((Subscript[ϵ, d])(Subscript[ϵ, m]))/((Subscript[ϵ, d])+(Subscript[ϵ, m])))

resonance[λ_, ϕ_, nm_, nd_, d_, m_] := Module[{nmetal, ϵm, ϵd, a, b, c, sinθ1, sinθ2, θ1, θ2},
   nmetal = Which[nm == -1, Gold[λ], nm == -2, Silver[λ], nm == -3, Copper[λ], nm == -4, Aluminum[λ], nm == -5, Drude[λ], nm == -6, LorentzDrudeAu[λ], nm == -7, LorentzDrudeAg[λ], nm == -8, LorentzDrudeCu[λ], nm == -9, LorentzDrudeAl[λ], nm == -10, Gold2[λ], nm == -11, Silver2[λ], nm == -12, LorentzDrudeBe[λ], nm == -13, LorentzDrudeCr[λ], nm == -14, LorentzDrudeNi[λ], nm == -15, LorentzDrudePd[λ], nm == -16, LorentzDrudePt[λ], nm == -17, LorentzDrudeTi[λ], nm == -18, LorentzDrudeW[λ], nm == -19, Iron2[λ], nm == -20, Cobalt2[λ], nm == -21, Rhodium[λ], nm == -22, Iridium[λ], True, nm];
   ϵm = nmetal^2;
   ϵd = nd^2;
   a = ϵd;
   b = (2 m nd λ/d) Sin[2 π ϕ/360];
   c = (m λ/d)^2 - ϵd Re[ϵm]/(ϵd + Re[ϵm]);
   sinθ1 = (-b + Sqrt[b^2 - 4 a c])/(2 a);
   sinθ2 = (-b - Sqrt[b^2 - 4 a c])/(2 a);
   θ1 = ArcSin[sinθ1];
   θ2 = ArcSin[sinθ2];
   Return[{180 θ1/π, 180 θ2/π}]];

Vector diffraction program setup

Parameters and functions used in the grating program.

initialize:=(
Clear[x,y,slope,offset,C0,D0,a,b,c,d];
Clear[A1,B1c,B1g,B1s,C1c,C1g,C1s,D1];
Clear[cp,dp,α,βc,βg,βs];
Clear[Ltop,Lbot,Ztop1,Zbot1];

ηc[1]=ηg[1]=ηs[1]=1.;
K:=2π/period;           (* grating vector *)
Nsegments=50;            (* # line segments for groove profile *)
N1 := norders+1;
N2 := 2 norders + 1;
N4 := 4 norders + 2;
z0=(Sqrt[μ0/ϵ0])/.NSub;    (* impedance of free space *)

(* x and y axis coordinates of piecewise linear groove profile *)
x[n_]:=n period/Nsegments;         
y[n_]:=0.5 amp(1-Cos[n 2π/Nsegments]); 
slope[n_]:=((y[n]-y[n-1])/(x[n]-x[n-1]));
offset[n_]:=y[n]-slope[n] x[n];
C0[n_]:=(1/(1+((slope[n])^2)));
D0[n_]:=((slope[n])/(1+((slope[n])^2)));
cp[p_]:=If[p==0,(∑_(n=1)^Nsegments)C0[n]((x[n]-x[n-1])/period),(∑_(n=1)^Nsegments)((I C0[n])/(2π p))((E^(-I p K x[n] ))-(E^(-I p K x[n-1] )))];
dp[p_]:=If[p==0,(∑_(n=1)^Nsegments)D0[n]((x[n]-x[n-1])/period),(∑_(n=1)^Nsegments)((I D0[n])/(2π p))((E^(-I p K x[n] ))-(E^(-I p K x[n-1] )))];
α[m_]:=nc  k Sin[π θi/180.0]+ m K;
βc[m_]:=(Sqrt[(nc^2) (k^2)-((α[m])^2)]); 
βg[m_,ng_]:=(Sqrt[(ng^2) (k^2)-((α[m])^2)]);
βs[m_,ns_]:=(Sqrt[(ns^2) (k^2)-((α[m])^2)]);
a[T_]:=Select[T,Re[#]>0&&Abs[Im[#]]<micro& ];
b[T_]:=Select[T,Im[#]>micro&];
c[T_]:=Select[T,Re[#]<0&&Abs[Im[#]]<micro&];
d[T_]:=Select[T,Im[#]<-micro&]; 
Ltop[m_,t_]:=(1/period)(∑_(n=1)^Nsegments)If[(slope[n]==0&&m==0)||(t==0&&m==0),
x[n]-x[n-1],(E^(I t offset[n]))(((E^(I (t slope[n]-m K)x[n]))-(E^(I (t slope[n]-m K)x[n-1])))/(I (t slope[n]-m K)))];
Lbot[m_,t_]:=(1/period)(∑_(n=1)^Nsegments)If[(slope[n]==0&&m==0)||(t==0&&m==0),
x[n]-x[n-1],(E^(I t offset[n]))(((E^(I (t slope[n]-m K)x[n]))-(E^(I (t slope[n]-m K)x[n-1])))/(I (t slope[n]-m K)))];
Ztop1:=Table[Ltop[m,-βc[0]],{m,-norders,norders}];
Zbot1:=Table[ηc[mode](-βc[0]+m K ((α[0])/(βc[0])))Lbot[m,-βc[0]],{m,-norders,norders}];)

Main vector diffraction program

This is the main routine "Grating" for calculating the vector diffraction amplitudes of reflectivity and transmissivity for each order.  The loop iterates over angle of incidence or wavelength.  The output is a list, {R, T}, where R is a list of {wavelength/angle, R} and T is a list of {wavelength/angle, T}.  If mode = 1, R and T are calculated for TE incidence.  If mode = 2, R and T are calculated for TM incidence.

Grating[θ1_,θ2_,θinc_,λ1_,λ2_,λinc_,nlayers_,thickness_,amp_,period_,norders_,nc_,ng_,ns_]:=Module[{Test,ng1,ns1,B1g,C1g,Tg1,Tg2,Tg,Tgrq,Tgmm,MMg,Φg,Tstmp,Ntleft,θt,Ntmp1,Ntmp2,Ntop,Nbot,Ntright,T0,TTEPoint,TTMPoint},
initialize;
If [θ1≠θ2,NumInc = Abs[(θ2 - θ1)/θinc] + 1,NumInc = Abs[(λ2 - λ1)/λinc] + 1];
Do[Do[{
If [θ1≠θ2,{θi=θ1 + (index-1) θinc,λ=λ1},{θi=θ1,λ=λ1 + (index-1) λinc}],
ng1=Which[ng==-1,Gold[λ],ng==-2,Silver[λ],ng==-3,Copper[λ],ng==-4,Aluminum[λ],ng==-5,Drude[λ],ng==-6,LorentzDrudeAu[λ],ng==-7,LorentzDrudeAg[λ],ng==-8,LorentzDrudeCu[λ],ng==-9,LorentzDrudeAl[λ],ng==-10,Gold2[λ],ng==-11,Silver2[λ],ng==-12,LorentzDrudeBe[λ],ng==-13,LorentzDrudeCr[λ],ng==-14,LorentzDrudeNi[λ],ng==-15,LorentzDrudePd[λ],ng==-16,LorentzDrudePt[λ],ng==-17,LorentzDrudeTi[λ],ng==-18,LorentzDrudeW[λ],ng==-19,Iron2[λ],ng==-20,Cobalt2[λ],ng==-21,Rhodium[λ],ng==-22,Iridium[λ],True,ng];
ns1=Which[ns==-1,Gold[λ],ns==-2,Silver[λ],ns==-3,Copper[λ],ns==-4,Aluminum[λ],ns==-5,Drude[λ],ns==-6,LorentzDrudeAu[λ],ns==-7,LorentzDrudeAg[λ],ns==-8,LorentzDrudeCu[λ],ns==-9,LorentzDrudeAl[λ],ns==-10,Gold2[λ],ns==-11,Silver2[λ],ns==-12,LorentzDrudeBe[λ],ns==-13,LorentzDrudeCr[λ],ns==-14,LorentzDrudeNi[λ],ns==-15,LorentzDrudePd[λ],ns==-16,LorentzDrudePt[λ],ns==-17,LorentzDrudeTi[λ],ns==-18,LorentzDrudeW[λ],ns==-19,Iron2[λ],ns==-20,Cobalt2[λ],ns==-21,Rhodium[λ],ns==-22,Iridium[λ],True,ns];
ηc[2]=1/(nc^2);ηg[2]=1/(ng1^2);ηs[2]=1/(ns1^2);
k=2π/λ,
γ1=k 1.z0,
γ2c=k (nc^2)/z0,
γ2g=k (ng1^2)/z0,
γ2s=k (ns1^2)/z0,

A1=Table[α[n-N1] dp[m-n],{m,1,N2},{n,1,N2}];
B1=Table[cp[m-n],{m,1,N2},{n,1,N2}];
B1c=B1/ηc[mode];
B1s=B1/ηs[mode];
C1=Table[-α[m-N1] α[n-N1] cp[m-n],{m,1,N2},{n,1,N2}];
C1c=Table[If[m==n,ηc[mode](γ1 γ2c+C1[[m,n]]),ηc[mode]C1[[m,n]]],{m,1,N2},{n,1,N2}];
C1s=Table[If[m==n,ηs[mode](γ1 γ2s+C1[[m,n]]),ηs[mode]C1[[m,n]]],{m,1,N2},{n,1,N2}];
If[nlayers==1,{B1g=B1/ηg[mode],C1g=Table[If[m==n,ηg[mode](γ1 γ2g+C1[[m,n]]),ηg[mode]C1[[m,n]]],{m,1,N2},{n,1,N2}]}];
D1=Table[α[m-N1] dp[m-n],{m,1,N2},{n,1,N2}];

Tc1 = Partition[Flatten[Riffle[A1,B1c]],N4];
Tc2 = Partition[Flatten[Riffle[C1c,D1]],N4];
Tc = Join[Tc1,Tc2];
Tcrq=Eigenvalues[Tc];
Tcmm=Eigenvectors[Tc];
(* Print["Tcrq=",MatrixForm[Chop[Tcrq]]]; *)

If[nlayers==1,{
Tg1 = Partition[Flatten[Riffle[A1,B1g]],N4],
Tg2 = Partition[Flatten[Riffle[C1g,D1]],N4],
Tg = Join[Tg1,Tg2],
Tgrq=Eigenvalues[Tg],
Tgmm=Eigenvectors[Tg]}];
(* Print["Tgrq=",MatrixForm[Chop[Tgrq]]]; *)

Ts1 = Partition[Flatten[Riffle[A1,B1s]],N4];
Ts2 = Partition[Flatten[Riffle[C1s,D1]],N4];
Ts = Join[Ts1,Ts2];
Tsrq=Eigenvalues[Ts];
Tsmm=Eigenvectors[Ts];
(* Print["Tsrq=",MatrixForm[Chop[Tsrq]]]; *)

abcd=Flatten[Append[Flatten[Append[Flatten[Append[a[Tcrq],b[Tcrq]]],c[Tcrq]]],d[Tcrq]]];
MMc=Transpose[Table[Tcmm[[Position[Tcrq,abcd[[i]]][[1,1]]]],{i,1,N4}]];
Tcrq = abcd;        (* Replace Tcrq with the reordered vector. *)
Φc=Table[If[i==j,1,0],{i,1,N4},{j,1,N4}];

If[nlayers==1,{
abcd=Flatten[Append[Flatten[Append[Flatten[Append[a[Tgrq],b[Tgrq]]],c[Tgrq]]],d[Tgrq]]],
MMg=Transpose[Table[Tgmm[[Position[Tgrq,abcd[[i]]][[1,1]]]],{i,1,N4}]],
Tgrq = abcd,
Φg=Table[If[i==j,(E^(I thickness abcd[[i]])),0],{i,1,N4},{j,1,N4}]}];

abcd=Flatten[Append[Flatten[Append[Flatten[Append[a[Tsrq],b[Tsrq]]],c[Tsrq]]],d[Tsrq]]];
MMs=Transpose[Table[Tsmm[[Position[Tsrq,abcd[[i]]][[1,1]]]],{i,1,N4}]];
Tsrq = abcd;
Φs=Table[If[i==j,1,0],{i,1,N4},{j,1,N4}];

(*** Set up Minc matrix for incident & reflected propagating/evanescent fields *****)
Clear[Tctmp];
Ztop=Transpose[Table[If[n==1,Ztop1,Table[0,{m,1,N2}]],{n,1,N2}]];
Zbot=Transpose[Table[If[n==1,Zbot1,Table[0,{m,1,N2}]],{n,1,N2}]];
Z = Join[Ztop,Zbot];

Pr=0;Per = 0;
Do[{If[Abs[Im[βc[n]]] < micro,
  {Pr=Pr+1,θr[Pr]=180.0 ArcTan[βc[n],α[n]]/π,Tctmp[Pr]=βc[n],Mtmp1[Pr]=Table[Ltop[m-n,βc[n]],{m,-norders,norders}],
Mtmp2[Pr]=Table[ηc[mode](βc[n] - (m - n) K α[n]/βc[n])Ltop[m-n,βc[n]],{m,-norders,norders}]}]},{n,-norders,norders}];
Mtop = Transpose[Table[Mtmp1[m],{m,1,Pr}]];
Mbot = Transpose[Table[Mtmp2[m],{m,1,Pr}]];
MMc1 = Join[Mtop,Mbot];
Do[If[Im[Tcrq[[n]]] > micro,{
Per=Per+1,Tctmp[Pr+Per]=Tcrq[n],MMc1=MapThread[Append,{MMc1,Transpose[MMc][[n]]}]}],{n,1,N4}];

If[Pr+Per≠N2,
{Do[{If[Abs[Im[Tcrq[[n]]]]<micro,
If[Re[Tcrq[[n]]]>0,
{Test=True,
Do[If[Abs[Tcrq[[n]]-βc[m]]<0.1,Test=False],{m,-norders,norders}],
If[Test,{Per=Per+1,Tctmp[Pr+Per]=Tcrq[[n]],
MMc1=MapThread[Append,{MMc1,Transpose[MMc][[n]]}]}]}]],
If[Pr+Per==N2,Break[]]},{n,1,N4}]}];
Minc = Transpose[Join[Transpose[MMc1],Transpose[Z]]];
Do[Tcrq[[n]]=Tctmp[n],{n,1,N4}];

(******** Set up Nt matrix for transmitted propagating/evanescent fields ********)
Clear[Tstmp];
Pt = 0;Ped = 0;
If[Im[ns1]>0,Nt=MMs,
{Ntleft=Transpose[Table[Transpose[MMs][[n]],{n,1,N2}]],Do[If[Abs[Im[βs[n,ns1]]] < micro,
  {Pt=Pt+1,θt[Pt]=180.0 ArcTan[βs[n,ns1],α[n]]/π,Tstmp[Pt]=βs[n,ns1],Ntmp1[Pt]=Table[Ltop[m-n,-βs[n,ns1]],{m,-norders,norders}],
Ntmp2[Pt]=Table[ηs[mode](-βs[n,ns1] + (m - n) K α[n]/βs[n,ns1])Ltop[m-n,-βs[n,ns1]],{m,-norders,norders}]}],{n,-norders,norders}],
Ntop=Transpose[Table[Ntmp1[m],{m,1,Pt}]],
Nbot=Transpose[Table[Ntmp2[m],{m,1,Pt}]],
Ntright=Join[Ntop,Nbot],
Do[If[Im[Tsrq[[n]]] < -micro,{Ped=Ped+1,Tstmp[Pt+Ped]=Tsrq[[n]],Ntright=MapThread[Append,{Ntright,Transpose[MMs][[n]]}]}],{n,1,N4}],
Nt=Transpose[Join[Transpose[Ntleft],Transpose[Ntright]]]}];

If[Pt+Ped≠N2,
{Do[{If[Abs[Im[Tsrq[[n]]]]<micro,
If[Re[Tsrq[[n]]]<0,
{Test=True,
Do[If[Abs[Tsrq[[n]]-βs[m,ns1]]<0.1,Test=False],{m,-norders,norders}],
If[Test,{Ped=Ped+1,Tstmp[Pt+Ped]=Tsrq[[n]],
Nt=MapThread[Append,{Nt,Transpose[MMs][[n]]}]}]}]],If[Pt+Ped==N2,Break[]]},{n,1,N4}]}];
Do[Tsrq[[n]]=Tstmp[n],{n,1,N4}];

(*******   Generate transfer matrices.   *******)
Mtemp=If[nlayers==0,Minc,MMg];
W[0] = Inverse[Nt].Mtemp;
H = IdentityMatrix[N4];
Htemp1 = MMg.Φg;
H = Htemp1.Inverse[MMg].H;
W[1] = Inverse[Htemp1].Minc;
S11[0] = IdentityMatrix[N2];
S12[0] = Table[0,{i,1,N2},{j,1,N2}];
S21[0] = S12[0];
S22[0] = S11[0];

Do[{W11[i-1] = W[i-1][[1;;N2,1;;N2]],
W12[i-1] = W[i-1][[1;;N2,N2+1;;N4]],
W21[i-1] = W[i-1][[N2+1;;N4,1;;N2]],
W22[i-1] = W[i-1][[N2+1;;N4,N2+1;;N4]],
Tinv = Inverse[W11[i-1] - S12[i-1].W21[i-1]],
S11[i] = Tinv.S11[i-1],
S12[i] = Tinv.(S12[i-1].W22[i-1] - W12[i-1]),
S21[i] = S21[i-1] + S22[i-1].W21[i-1].S11[i],
S22[i] = S22[i-1].W22[i-1] + S22[i-1].W21[i-1].S12[i]},{i,1,nlayers+1}];
Do[If[Abs[θr[m]-θi]<0.01,m0=m],{m,1,Pr}];
If[θ1≠θ2,R0={θr[m0],(((Abs[S12[nlayers+1][[m0,1]]]))^2)},R0={1000 λ,(((Abs[S12[nlayers+1][[m0,1]]]))^2)}];
If[mode==1,{If[index==1,RTEPoint=R0,RTEPoint=Join[RTEPoint,R0]],RTEg=Partition[RTEPoint,2]},{If[index==1,RTMPoint=R0,RTMPoint=Join[RTMPoint,R0]],RTMg=Partition[RTMPoint,2]}];
If[Im[ns1]<micro,
{Do[If[Abs[ns1 Sin[π θt[m]/180.0]-nc Sin[π θi/180.0]]<0.01,m0=m],{m,1,Pt}],
If[θ1≠θ2,T0={θt[m0],((ns1 Cos[π θt[m0]/180.0])/(nc Cos[π θi/180.0]))(((Abs[S22[nlayers+1][[m0,1]]]))^2)},T0={1000 λ,((ns1 Cos[π θt[m0]/180.0])/(nc Cos[π θi/180.0]))(((Abs[S22[nlayers+1][[m0,1]]]))^2)}],
If[mode==1,{If[index==1,TTEPoint=T0,TTEPoint=Join[TTEPoint,T0]],TTEg=Partition[TTEPoint,2]},{If[index==1,TTMPoint=T0,TTMPoint=Join[TTMPoint,T0]],TTMg=Partition[TTMPoint,2]}]}]},
{index,1,NumInc}],{mode,1,2}];
Return[{RTEg,RTMg,TTEg,TTMg}]];

Output

(* Dummy values to suppress error message for first point in graphs: *)
RTEg={0,0};RTMg={0,0};TTEg={0,0};TTMg={0,0};

Dynamic[ListLinePlot[{RTEg,RTMg},PlotLegend->{"TE","TM"},LegendPosition-> {-0.25,-1.1},LegendSize->{0.3,0.25},LegendShadow->{.01,-.01},FrameLabel->{"Angle (°) or Wavelength (nm)","Reflectivity"},GridLines->Automatic,Frame->True,PlotRange->All,PlotStyle->Thick,PlotLabel->"Zero Order Reflectivity"]]

Dynamic[ListLinePlot[{TTEg,TTMg},PlotLegend->{"TE","TM"},LegendPosition-> {-0.25,-1.1},LegendSize->{0.3,0.25},LegendShadow->{.01,-.01},FrameLabel->{"Angle (°) or Wavelength (nm)","Transmissivity"},GridLines->Automatic,Frame->True,PlotRange->All,PlotStyle->Thick,PlotLabel->"Zero Order Transmissivity"]]

Figures

10.2  Otto technique

The routines for computing the reflectivity in the Otto configuration are RTM and RTE.  Reflectivity is computed for the field in the plane of incidence for RTM and perpendicular to the plane of incidence for RTE.  The input parameters are the wavelength λ, the angle of incidence θ of the beam on the bottom surface of the prism, the index of the prism (Subscript[n, 0]), the index of the dielectric in the gap (Subscript[n, 1]), and the index of the metal (Subscript[n, 2]), and the thickness of the dielectric film d.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

θ1 = 40;      (* initial angle of incidence within prism *)
θ2 = 45;      (* final angle of incidence within prism *)
θinc=0.01;    (* increment of angle  *)
λ=0.8;            (* wavelength in μm  *)
num = Round[Abs[θ2 - θ1]/θinc] + 1;
n0=1.5;   (*  index of glass prism  *)
n1=1.0;    (*  index of dielectric in gap between prism and metal surface  *)
n2=-6;   (*  index of gold vs. wavelength  *)

d=0.6;   (*  thickness of dielectric in gap in μm  *)
d1=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=1.0;   (*  thickness of dielectric in gap in μm  *)
d2=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=1.4;   (*  thickness of dielectric in gap in μm  *)
d3=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=1.0;   (*  thickness of dielectric in gap in μm  *)
d4=Table[{θ=θ1+(i-1)θinc,RTE[λ,θ,n0,n1,n2,d]},{i,num}];

Fig3=ListLinePlot[{d1,d2,d3,d4}, LegendPosition -> {-0.7, -0.35},LegendLabel->"Gap", PlotLegend->{"0.6 μm, TM","1.0 μm, TM","1.4 μm, TM","1.0 μm, TE"},LegendSize -> {0.5,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[option]] 
Fig3BW=ListLinePlot[{d1,d2,d3,d4}, LegendPosition -> {-0.7, -0.35},LegendLabel->"Gap", PlotLegend->{"0.6 μm, TM","1.0 μm, TM","1.4 μm, TM","1.0 μm, TE"},LegendSize -> {0.5,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.3.JPG"}],Fig3BW,ImageResolution->600];
Clear[θ,θ1,θ2,θinc,λ,num,n0,n1,n2,d1,d2,d3,d4];

The function for computing the dispersion of a lossy Drude metal with a real frequency is (Subscript[K, sp]).  The parameters that get passed to the routine are the normalized frequency Ω and the normalized damping coefficient Γ.

Γ=0;
Ω1 = 0.0;      (* initial normalized wavevector *)
Ω2 = 1/Sqrt[2]-0.01;      (* final normalized wavelength *)
Ωinc=0.01;    (* increment of wavevector  *)
num = Round[Abs[Ω2 - Ω1]/Ωinc] + 1;
ReΩ1=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];
Vacuum={{0.0,0.0},{10,10}};
Prism={{0.0,0.0},{15,10}};

Fig4=ListLinePlot[{ReΩ1,Vacuum,Prism}, LegendPosition -> {0.3, -0.2}, PlotLegend->{"SP","Dielectric","Prism"},LegendSize -> {0.45,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"K", "Ω"}, LabelStyle -> Bold, GridLines -> Automatic, Frame -> True, PlotRange->{{0,3},{0,1}},Evaluate[option]]
Fig4BW=ListLinePlot[{ReΩ1,Vacuum,Prism}, LegendPosition -> {0.3, -0.2}, PlotLegend->{"SP","Dielectric","Prism"},LegendSize -> {0.45,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"K", "Ω"}, LabelStyle -> Bold, GridLines -> Automatic, Frame -> True, PlotRange->{{0,3},{0,1}},Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.4.JPG"}],Fig4BW,ImageResolution->600];
Clear[Ω,Ω1,Ω2,Ωinc,Γ,num,ReΩ1,Vacuum,Prism];

10.3  Kretschmann technique

The routines for computing the reflectivity in the Otto configuration are RTM and RTE.  Reflectivity is computed for the field in the plane of incidence for RTM and perpendicular to the plane of incidence for RTE.  The input parameters are the wavelength λ, the angle of incidence θ of the beam on the bottom surface of the prism, the index of the prism (Subscript[n, 0]), the index of the dielectric in the gap (Subscript[n, 1]), and the index of the metal (Subscript[n, 2]), and the thickness of the dielectric film d.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

θ1 = 40;      (* initial angle of incidence within prism *)
θ2 = 50;      (* final angle of incidence within prism *)
θinc=0.01;    (* increment of angle  *)
λ=0.8;            (* wavelength in μm  *)
num = Round[Abs[θ2 - θ1]/θinc] + 1;
n0=1.5;   (*  index of glass prism  *)
n1=-6;    (*  index of gold film  *)
n2=1.0;   (*  index of dielectric below gold film  *)

d=0.02;   (*  thickness of dielectric in gap in μm  *)
d1=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=0.03;   (*  thickness of dielectric in gap in μm  *)
d2=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=0.05;   (*  thickness of dielectric in gap in μm  *)
d3=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=0.07;   (*  thickness of dielectric in gap in μm  *)
d4=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];

Fig6=ListLinePlot[{d1,d2,d3,d4}, LegendPosition -> {0.3, -0.3},LegendLabel->"Thickness", PlotLegend->{"0.02 μm","0.03 μm","0.05 μm","0.07 μm"},LegendSize -> {0.4,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[option]]
Fig6BW=ListLinePlot[{d1,d2,d3,d4}, LegendPosition -> {0.3, -0.3},LegendLabel->"Thickness", PlotLegend->{"0.02 μm","0.03 μm","0.05 μm","0.07 μm"},LegendSize -> {0.4,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.6.JPG"}],Fig6BW,ImageResolution->600];
Clear[θ,θ1,θ2,θinc,λ,num,n0,n1,n2,d1,d2,d3,d4];

10.4 Diffraction gratings and Wood's anomalies

The parameters that are set before a calculation are as follows:
	(Subscript[θ, 1]): initial angle (°)
	(Subscript[θ, 2]): final angle (°)
	(Subscript[θ, inc]): increment of angle (°)
	(Subscript[λ, 1]): initial wavelength (μm)
	(Subscript[λ, 2]): final wavelength (μm)
	(Subscript[λ, inc]): increment of wavelength (μm)
	nlayers: # thin film layers (either 0 or 1)
	thickness: thickness of thin film at x = 0 in μm
	amp: peak-to-peak amplitude of sinusoidal groove profile in μm
	period: period of grating in μm
	norders: # of diffraction orders to retain in the calculation
	nc: refractive index of cover (upper dielectric)
	ng: refractive index of guide layer
	ns: refractive index of substrate

The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

After setting up the desired parameters, the Grating routine is executed.  Because it can take many minutes to generate a complete graph, the output is plotted in real time in the Code/Setup/Output section where the reader can watch the graph as it is being generated.

θ1 = 17.;         (* starting angle of incidence *)
θ2 = 23.;                (* final angle of incidence *)
(*θinc = 0.01;              (* increment of angle *)*)
θinc = 0.02;              (* increment of angle *)
λ1 = .700 ;                  (* initial free space wavelength in microns *)
λ2 = .700;                  (* final free space wavelength in microns *)
λinc = 0.0;                  (* increment in wavelength in microns *)
nlayers =0;               (* # of thin films ≤ 1*)
thickness = 0.0;      (* thickness of guide thin film in microns *)
amp = 0.050;                (* peak-to-peak amplitude of sinusoidal grating in microns *)
period = 1.0 ;                 (* grating period length in microns*)
norders = 2;                 (* # orders to include in calculation *)
nc = 1.0;                             (* index of cover layer *)
ng = 1.0;                          (* index of guide thin film *)
ns = -6;                          (* Drude-Lorentz fit to refractive index of gold *)

(* Dummy values to suppress error message for first point in graphs: *)
RTEg={0,0};RTMg={0,0};TTEg={0,0};TTMg={0,0};  {RTEg,RTMg,TTEg,TTMg}=Grating[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng,ns];
Fig11=ListLinePlot[{RTEg, RTMg}, PlotLegend -> {" TE", " TM"}, LegendPosition -> {0.4, -.2}, LegendSize -> {0.25,0.25}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"}, Evaluate[option]]
Fig11BW=ListLinePlot[{RTEg, RTMg}, PlotLegend -> {" TE", " TM"}, LegendPosition -> {0.2, -.1}, LegendSize -> {0.25,0.25}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"}, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.11.JPG"}],Fig11BW,ImageResolution->600];
Clear[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng,ns];

The parameters that are set before a calculation are as follows:
	(Subscript[θ, 1]): initial angle (°)
	(Subscript[θ, 2]): final angle (°)
	(Subscript[θ, inc]): increment of angle (°)
	(Subscript[λ, 1]): initial wavelength (μm)
	(Subscript[λ, 2]): final wavelength (μm)
	(Subscript[λ, inc]): increment of wavelength (μm)
	nlayers: # thin film layers (either 0 or 1)
	thickness: thickness of thin film at x = 0 in μm
	amp: peak-to-peak amplitude of sinusoidal groove profile in μm
	period: period of grating in μm
	norders: # of diffraction orders to retain in the calculation
	nc: refractive index of cover (upper dielectric)
	ng: refractive index of guide layer
	ns: refractive index of substrate
	mode: 1 for TE incidence, 2 for TM incidence

θ1 = 19.8;         (* starting angle of incidence *)
θ2 = 19.8;                (* final angle of incidence *)
θinc = 0.00;              (* increment of angle *)
λ1 = .650 ;                  (* initial free space wavelength in microns *)
λ2 = .740;                  (* final free space wavelength in microns *)
λinc = 0.001;                  (* increment in wavelength in microns *)
nlayers =0;               (* # of thin films ≤ 1*)
thickness = 0.0;      (* thickness of guide thin film in microns *)
amp = 0.050;                (* peak-to-peak amplitude of sinusoidal grating in microns *)
period = 1.0 ;                 (* grating period length in microns*)
norders = 2;                 (* # orders to include in calculation *)
nc = 1.0;                             (* index of cover layer *)
ng = 1.0;                          (* index of guide thin film *)
ns=-6;(* index of substrate *)

{RTEg,RTMg,TTEg,TTMg}=Grating[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng,ns];
Fig12=ListLinePlot[{RTEg, RTMg},PlotLegend -> {" TE", " TM"}, LegendPosition -> {0.3, -.3}, LegendSize -> 0.3, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (nm)", "R"}, Evaluate[option]]
Fig12BW=ListLinePlot[{RTEg, RTMg},PlotLegend -> {" TE", " TM"}, LegendPosition -> {0.3, -.2}, LegendSize -> 0.3, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (nm)", "R"}, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.12.JPG"}],Fig12BW,ImageResolution->600];
Clear[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng,ns];

Routine "resonance" computes the approximate polar resonance angle θ as a function of the wavelength λ, the azimuthal angle ϕ, the refractive index of the metal (Subscript[n, m]), the refractive index of the surrounding dielectric medium (Subscript[n, d]), the grating period d, and the order m.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

ϕ1 = 0;      (* initial azimuthal angle of incidence *)
ϕ2 = 90;      (* final angle of incidence within prism *)
ϕinc=5;    (* increment of angle  *)
λ=0.7;            (* wavelength in μm  *)
nm=-6;         (*  refractive index choice for gold  *)
nd=1;           (*  refractive index for surrounding medium  *)
d=1.0;          (*  grating period in μm *)
m=1;             (*  order  *)
num = Round[Abs[ϕ2 - ϕ1]/ϕinc] + 1;

θvd={{0,50.3},{5,45.2},{10,41.1},{15,37.5},{20,34.5},{25,31.9},{30,29.7},{35,27.8},{40,26.2},{45,25.0},{50,23.8},{55,22.9},{60,22.1},{65,21.4},{70,20.9},{75,20.5},{80,20.2},{85,20.0},{90,20.0}};
data=Table[resonance[0.7,ϕϕ,.257+I 3.717,1,1,1][[1]], {ϕϕ, 0, 90}];
Fig18=ListPlot[{data,θvd}, LegendPosition -> {-0.1, 0.3}, PlotLegend->{"approximate","vector diffraction"},LegendSize -> {0.7,0.2}, LegendShadow -> {.01, -.01},Frame->True,PlotStyle->{Thick,Directive[Red,PointSize[0.02]]},Joined->{True,False},FrameLabel -> {"ϕ (°)", "θ (°)"},Evaluate[option]]
Fig18BW=ListPlot[{data,θvd}, LegendPosition -> {-0.1, 0.3}, PlotLegend->{"approximate","vector diffraction"},LegendSize -> {0.7,0.2}, LegendShadow -> {.01, -.01},Frame->True,PlotStyle->{Directive[Black,Thick],Directive[Black,PointSize[0.02]]},Joined->{True,False},FrameLabel -> {"ϕ (°)", "θ (°)"},Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.18.JPG"}],Fig18BW,ImageResolution->600];
Clear[θ,ϕ1,ϕ2,ϕinc,λ,num,nm,nd,d,m,data,θvd];

rmin = {{10, 0.883}, {15, 0.828}, {20,0.756}, {25,0.672}, {30, 0.580}, {35, 0.487}, {40, 0.396}, {45, 0.313}, {50,0.240}, {55, 0.178}, {60,0.119}, {65, 0.074}, {70, 0.043}, {75, 0.025}, {80,0.007}, {85, 0.001}, {90, 0.004},{95,0.006},{100,0.017},{105,0.026},{110,0.040},{115,0.051},{120,0.063},{125,0.076},{130,0.086},{135,0.096},{140,0.106},{145,0.115},{150,0.125}};
Fig19=ListPlot[rmin, PlotStyle -> Blue,FrameLabel -> {"d (nm)", "(Subscript[R, min])"}, Joined->True,PlotMarkers->Automatic,Evaluate[option]]
Fig19BW=ListPlot[rmin,FrameLabel -> {"d (nm)", "(Subscript[R, min])"}, Joined->True,PlotMarkers->Automatic,Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.19.JPG"}],Fig19BW,ImageResolution->600];
Clear[rmin];

The parameters that are set before a calculation are as follows:
	(Subscript[θ, 1]): initial angle (°)
	(Subscript[θ, 2]): final angle (°)
	(Subscript[θ, inc]): increment of angle (°)
	(Subscript[λ, 1]): initial wavelength (μm)
	(Subscript[λ, 2]): final wavelength (μm)
	(Subscript[λ, inc]): increment of wavelength (μm)
	nlayers: # thin film layers (either 0 or 1)
	thickness: thickness of thin film at x = 0 in μm
	amp: peak-to-peak amplitude of sinusoidal groove profile in μm
	period: period of grating in μm
	norders: # of diffraction orders to retain in the calculation
	nc: refractive index of cover (upper dielectric)
	ng: refractive index of guide layer
	ns: refractive index of substrate

The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

θ1 = 17.;             (* starting angle of incidence *)
θ2 = 23.;                (* final angle of incidence *)
θinc =0.02;              (* increment of angle *)
λ1 = .700 ;                  (* initial free space wavelength in microns *)
λ2 = .700;                  (* final free space wavelength in microns *)
λinc = 0.0;                  (* increment in wavelength in microns *)
nlayers =0;               (* # of thin films ≤ 1*)
thickness = 0.0;      (* thickness of guide thin film in microns *)
amp = 0.050;                (* peak-to-peak amplitude of sinusoidal grating in microns *)
period = 1.0 ;                 (* grating period length in microns*)
norders = 2;                 (* # orders to include in calculation *)
nc = 1.0;                             (* index of cover layer *)
ng = 1.0;                          (* index of guide thin film *)
 
ns = -6;                          (* Drude-Lorentz fit to refractive index of gold *)
{RTEAu,RTMAu,TTEAu,TTMAu}=Grating[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng,ns];
ns = -7;                          (* Drude-Lorentz fit to refractive index of gold *)
{RTEAg,RTMAg,TTEAg,TTMAg}=Grating[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng,ns];
Fig21=ListLinePlot[{RTMAu, RTMAg}, PlotLegend -> {" Au", " Ag"}, LegendPosition -> {0.3, -.3}, LegendSize -> {0.25,0.25}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[option]]
Fig21BW=ListLinePlot[{RTMAu, RTMAg}, PlotLegend -> {" Au", " Ag"}, LegendPosition -> {0.15, -.1}, LegendSize -> {0.25,0.25}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.21.JPG"}],Fig21BW,ImageResolution->600];
Clear[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng];

The parameters that are set before a calculation are as follows:
	(Subscript[θ, 1]): initial angle (°)
	(Subscript[θ, 2]): final angle (°)
	(Subscript[θ, inc]): increment of angle (°)
	(Subscript[λ, 1]): initial wavelength (μm)
	(Subscript[λ, 2]): final wavelength (μm)
	(Subscript[λ, inc]): increment of wavelength (μm)
	nlayers: # thin film layers (either 0 or 1)
	thickness: thickness of thin film at x = 0 in μm
	amp: peak-to-peak amplitude of sinusoidal groove profile in μm
	period: period of grating in μm
	norders: # of diffraction orders to retain in the calculation
	nc: refractive index of cover (upper dielectric)
	ng: refractive index of guide layer
	ns: refractive index of substrate

(*************   Parameters to be set by user *************)
θ1 = 17.;                  (* starting angle of incidence *)
θ2 = 35.;                (* final angle of incidence *)
θinc = 0.02;              (* increment of angle *)
λ1 = .700 ;                  (* initial free space wavelength in microns *)
λ2 = .700;                  (* final free space wavelength in microns *)
λinc = 0.0;                  (* increment in wavelength in microns *)
nlayers =1;               (* # of thin films ≤ 1*)
thickness = 0.025;      (* thickness of guide thin film in microns *)
amp = 0.110;                (* peak-to-peak amplitude of sinusoidal grating in microns *)
period = 1.1 ;                 (* grating period length in microns*)
norders = 2;                 (* # orders to include in calculation *)
nc = 1.0;                             (* index of cover layer *)
ng = 0.14+I(4.523);      (* index of guide thin film *)
ns=1.0;                            (* index of substrate *)

RTEg={0,0};RTMg={0,0};TTEg={0,0};TTMg={0,0};  
{RTEg,RTMg,TTEg,TTMg}=Grating[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng,ns];
Fig23=ListLinePlot[{RTMg, TTMg},PlotLegend -> {" Refl", " Tran"}, LegendPosition -> {0.45, -.1}, LegendSize -> {0.25,0.25}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R & T"},Evaluate[option]]
Fig23BW=ListLinePlot[{RTMg, TTMg},PlotLegend -> {" Refl", " Tran"}, LegendPosition -> {0.2, 0}, LegendSize -> {0.25,0.25}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R & T"},Evaluate[optB]];

Export[FileNameJoin[{base,"Fig10.23.JPG"}],Fig23BW,ImageResolution->600];
Clear[θ1,θ2,θinc,λ1,λ2,λinc,nlayers,thickness,amp,period,norders,nc,ng,ns];

10.11  Exercises

θ1 = 40;      (* initial angle of incidence within prism *)
θ2 = 50;      (* final angle of incidence within prism *)
θinc = 0.01;    (* increment of angle  *)
λ = 0.9;            (* wavelength in μm  *)
num = Round[Abs[θ2 - θ1]/θinc] + 1;
n0 = 1.5;   (*  index of glass prism  *)
n1 = -9;    (*  index of aluminum film  *)
n2 = 1.0;   (*  index of dielectric below gold film  *)

d = 0.006;   (*  thickness of dielectric in gap in μm  *)
d1 = Table[{θ = θ1 + (i - 1) θinc, RTM[λ, θ, n0, n1, n2, d]}, {i, num}];
d = 0.008;   (*  thickness of dielectric in gap in μm  *)
d2 = Table[{θ = θ1 + (i - 1) θinc, RTM[λ, θ, n0, n1, n2, d]}, {i, num}];
d = 0.010;   (*  thickness of dielectric in gap in μm  *)
d3 = Table[{θ = θ1 + (i - 1) θinc, RTM[λ, θ, n0, n1, n2, d]}, {i, num}];
d = 0.012;   (*  thickness of dielectric in gap in μm  *)
d4 = Table[{θ = θ1 + (i - 1) θinc, RTM[λ, θ, n0, n1, n2, d]}, {i, num}];

ex2=ListLinePlot[{d1, d2, d3, d4}, LegendPosition -> {0.2, 0.1}, LegendLabel -> "Thickness", PlotLegend -> {"0.006 μm", "0.008 μm", "0.010 μm", "0.012 μm"}, LegendSize -> {0.5, 0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[option]]
ex2BW=ListLinePlot[{d1, d2, d3, d4}, LegendPosition -> {0.2, 0.1}, LegendLabel -> "Thickness", PlotLegend -> {"0.006 μm", "0.008 μm", "0.010 μm", "0.012 μm"}, LegendSize -> {0.5, 0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[optB]];

Export[FileNameJoin[{base,"ex2.JPG"}],ex2BW,ImageResolution->600];
Clear[θ, θ1, θ2, θinc, λ, num, n0, n1, n2,d, d1,d2,d3,d4];

θ1 = 30;      (* initial angle of incidence within prism *)
θ2 = 60;      (* final angle of incidence within prism *)
θinc=0.01;    (* increment of angle  *)
λ=0.32;            (* wavelength in μm  *)
num = Round[Abs[θ2 - θ1]/θinc] + 1;
n0=1.5;   (*  index of glass prism  *)
n1=-7;    (*  index of silver film  *)
n2=1.0;   (*  index of dielectric below gold film  *)

d=0.01;   (*  thickness of dielectric in gap in μm  *)
d1=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=0.02;   (*  thickness of dielectric in gap in μm  *)
d2=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=0.03;   (*  thickness of dielectric in gap in μm  *)
d3=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d=0.04;   (*  thickness of dielectric in gap in μm  *)
d4=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];

ex3a=ListLinePlot[{d1,d2,d3,d4}, LegendPosition -> {-0.6,0.2},LegendLabel->"Thickness", PlotLegend->{"0.02 μm","0.03 μm","0.04 μm","0.05 μm"},LegendSize -> {0.4,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[option]] 
ex3aBW=ListLinePlot[{d1,d2,d3,d4}, LegendPosition -> {-0.6,0.2},LegendLabel->"Thickness", PlotLegend->{"0.02 μm","0.03 μm","0.04 μm","0.05 μm"},LegendSize -> {0.4,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"},Evaluate[optB]] ;

Export[FileNameJoin[{base,"ex3a.JPG"}],ex3aBW,ImageResolution->600];
Clear[θ,θ1,θ2,θinc,λ,num,n0,n1,n2,d,d1,d2,d3,d4];

z1 = -0.1;      (* initial distance below metal film in microns *)
z2 = 0.5;      (* final distance below metal film in microns *)
zinc=0.001;    (* increment of distance  *)
n0=1.5;         (*  index of incident medium - glass prism  *)
n1=-7;            (*  index of silver  *)
n2=1.0;          (*  index of transmitting medium - air  *)
λ=0.32;            (*  wavelength in microns  *)

num = Round[Abs[z2 - z1]/zinc] + 1;
n1=-7;             (*  silver  *)
d=0.03;        (*  thickness of metal film in microns  *)
θ=45.0;         (*  angle of incidence  *)
Ag=Table[z=z1+(i-1)zinc;field[λ,θ,n0,n1,n2,d,z],{i,num}];

ex3b=ListLinePlot[{Ag[[All,1]],Ag[[All,2]],Ag[[All,3]]}, LegendPosition -> {.4, 0}, PlotLegend->{"|(Subscript[E, x])|","|(Subscript[E, z])|","|E|"},LegendSize ->{0.2,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"z (μm)", "|Field Amplitudes|"},Evaluate[option]]
ex3bBW=ListLinePlot[{Ag[[All,1]],Ag[[All,2]],Ag[[All,3]]}, LegendPosition -> {.4, 0}, PlotLegend->{"|(Subscript[E, x])|","|(Subscript[E, z])|","|E|"},LegendSize ->{0.2,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"z (μm)", "|Field Amplitudes|"},Evaluate[optB]];

Export[FileNameJoin[{base,"ex3b.JPG"}],ex3bBW,ImageResolution->600];
Clear[λ,θ,n0,n1,n2,d,z1,z2,zinc,num,Ag];

The routines for computing the reflectivity in the Kretschmann configuration are RTM and RTE.  Reflectivity is computed for the field in the plane of incidence for RTM and perpendicular to the plane of incidence for RTE.  The input parameters are the wavelength λ, the angle of incidence θ of the beam on the bottom surface of the prism, the index of the prism (Subscript[n, 0]), the index of the metal film (Subscript[n, 1]), and the index of the dielectric below the metal (Subscript[n, 2]), and the thickness of the metal film d.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

λ1 = 0.3;      (* initial wavelength in μm *)
λ2 = 1.2;      (* final wavelength μm *)
λinc=0.001;    (* increment of wavelength  *)

num = Round[Abs[λ2 - λ1]/λinc] + 1;
n0=1.5;   (*  index of glass prism  *)
n2=1.0;   (*  index of dielectric below metal film  *)
d=0.05;   (*  thickness of metal film in μm  *)
n1=-7;    (*  index of silver film  *)

θ=43;            (* angle of incidence in deg. *)
d1=Table[{λ=λ1+(i-1)λinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
θ=45;            (* angle of incidence in deg. *)
d2=Table[{λ=λ1+(i-1)λinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
ex4=ListLinePlot[{d1,d2}, LegendPosition -> {0.45, -0.3}, PlotLegend->{"43°","45°"},LegendSize -> {0.25,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "R"}, Evaluate[option]] 
ex4BW=ListLinePlot[{d1,d2}, LegendPosition -> {0.4, -0.25}, PlotLegend->{"43°","45°"},LegendSize -> {0.25,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "R"}, Evaluate[optB]] ;

Export[FileNameJoin[{base,"ex4.JPG"}],ex4BW,ImageResolution->600];
Clear[θ,θ1,θ2,θinc,λ,num,n0,n1,n2,d1,d2,d3,d4];

The routine for computing the reflectivity in the Otto configuration is r.  The parameters that get pased to the routine are the normalized wavevector K and the normalized damping coefficient Γ.  The input parameters are the wavelength λ, the angle of incidence θ of the beam on the bottom surface of the prism, the index of the prism (Subscript[n, 0]), the index of the dielectric in the gap (Subscript[n, 1]), and the index of the metal (Subscript[n, 2]), and the thickness of the dielectric film d.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

θ1 = 68;      (* initial angle of incidence within prism *)
θ2 = 80;      (* final angle of incidence within prism *)
θinc=0.01;    (* increment of angle  *)
λ=0.8;            (* wavelength in μm  *)
num = Round[Abs[θ2 - θ1]/θinc] + 1;
n0=1.5;   (*  index of glass prism  *)
n1=1.38;    (*  index of dielectric in gap between prism and metal surface  *)
n2=-6;   (*  index of gold vs. wavelength  *)
d=0.55;   (*  thickness of dielectric in gap in μm  *)
d1=Table[{θ=θ1+(i-1)θinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
d2=Table[{θ=θ1+(i-1)θinc,RTE[λ,θ,n0,n1,n2,d]},{i,num}];

ex5=ListLinePlot[{d1,d2}, LegendPosition -> {-0.6, -0.3}, PlotLegend->{"TM","TE"},LegendSize -> {0.25,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"}, Evaluate[option]] 
ex5BW=ListLinePlot[{d1,d2}, LegendPosition -> {-0.6, -0.3}, PlotLegend->{"TM","TE"},LegendSize -> {0.25,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"θ (°)", "R"}, Evaluate[optB]] ;

Export[FileNameJoin[{base,"ex5.JPG"}],ex5BW,ImageResolution->600];
Clear[θ,θ1,θ2,θinc,λ,num,n0,n1,n2,d1,d2];

The routines for computing the reflectivity in the Kretschmann configuration are RTM and RTE.  Reflectivity is computed for the field in the plane of incidence for RTM and perpendicular to the plane of incidence for RTE.  The input parameters are the wavelength λ, the angle of incidence θ of the beam on the bottom surface of the prism, the index of the prism (Subscript[n, 0]), the index of the metal film (Subscript[n, 1]), and the index of the dielectric below the metal (Subscript[n, 2]), and the thickness of the metal film d.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

λ1 = 0.3;      (* initial wavelength in μm *)
λ2 = 1.2;      (* final wavelength μm *)
λinc=0.001;    (* increment of wavelength  *)
θ=43;            (* angle of incidence in deg. *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
n0=1.5;   (*  index of glass prism  *)
n2=1.0;   (*  index of dielectric below metal film  *)
d=0.05;   (*  thickness of metal film in μm  *)

n1=-6;    (*  index of gold film  *)
d1=Table[{λ=λ1+(i-1)λinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
n1=-7;    (*  index of silver film  *)
d2=Table[{λ=λ1+(i-1)λinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
n1=-8;    (*  index of copper film  *)
d3=Table[{λ=λ1+(i-1)λinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];
n1=-9;    (*  index of aluminum film  *)
d=0.018;   (*  thickness of metal film in μm  *)
d4=Table[{λ=λ1+(i-1)λinc,RTM[λ,θ,n0,n1,n2,d]},{i,num}];

ex6=ListLinePlot[{d1,d2,d3,d4}, LegendPosition -> {0.45, -0.3}, PlotLegend->{"Au","Ag","Cu","Al"},LegendSize -> {0.25,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "R"}, Evaluate[option]]
ex6BW=ListLinePlot[{d1,d2,d3,d4}, LegendPosition -> {0.45, -0.3}, PlotLegend->{"Au","Ag","Cu","Al"},LegendSize -> {0.25,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "R"}, Evaluate[optB]];

Export[FileNameJoin[{base,"ex6.JPG"}],ex6BW,ImageResolution->600];
Clear[θ,θ1,θ2,θinc,λ,num,n0,n1,n2,d1,d2,d3,d4];

Black and white figures

Fig3BW

Fig4BW

Fig6BW

Fig11BW

Fig12BW

Fig18BW

Fig19BW

Fig21BW

Fig23BW

ex2BW

ex3aBW

ex3bBW

ex4BW

ex5BW

ex6BW

Figures

10.1  Introduction

10.2  Otto configuration

10.3  Kretschmann configuration

10.4  Diffraction gratings and Wood's anomalies

10.6  End fire coupling

10.7  Near field launching

Appendix:  description of grating code

In this section we discuss the Chandezon technique [31] for solving the vector diffraction of a plane wave incident upon a multilayer conformal thin film stack on a grating.  As originally described, this technique was derived for conformal multilayer gratings, i.e., multilayer gratings of thin films in which the groove profile was conformal in each layer.  A very useful enhancement of the technique was described by Priest et al. [32] for multilayer films with periodic grooves which are not conformal from layer to layer, and we include this enhancement in the discussion.  Finally, the transfer matrix technique as originally used by Chandezon et al., has been modified into a scattering matrix technique by Cotter et al. [33] in order to handle films which are exceptionally thick or grooves which are exceptionally deep without exceeding the numerical precision of conventional PC's.  This option is also described and implemented.

We should also point out that rigorous vector diffraction theory has undergone an immense development in the last two decades by many research groups.  Both the Chandezon technique and the other popular approach, the Fourier modal method, have been greatly improved and optimized for efficient computation of the diffracted fields.  The following discussion is not meant to incorporate all of these improvements and advances, but to simply provide an introduction to the Chandezon technique so that the interested reader may more easily understand the literature.

10.5.1  Definition of groove profile

The basic geometry of the problem is shown in Fig. 10.31.  As before, the grooves are parallel to the z axis and perpendicular to the x axis.  The y axis is normal to the surface of the grating.  The groove profile in each layer has a period d.  The profile itself is allowed to vary from layer to layer, however.   For the bottom of layer j the functional dependence of the groove profile is

(Subscript[y, j])=(Subscript[d, j]) + (Subscript[a, j])(x).

We choose a(x=0) = 0 for each layer, so

(Subscript[d, j])=(Subscript[e, 1]) + (Subscript[e, 2])+ ··· + (Subscript[e, j-1])

with

(Subscript[d, 0])=(Subscript[d, 1]) =0.

Subscript 0 denotes the substrate and subscript Q+1 represents the upper dielectric for Q layers.

For the purposes of this implementation, we assume the groove profile in each layer is composed of piecewise continuous line segments with the same number of line segments and the same x coordinate end points as shown in Fig. 10.31.  We should point out that the original Chandezon technique described here is more suited for a continuously differentiable groove profile.  The piecewise linear groove profile generally requires a larger number of terms in the Fourier expansions to converge to an accurate result.  This limitation has been removed by Plumey et al., [34] but their work is beyond the scope of this discussion.  Nevertheless, in practice even the groove profile shown in Fig. 10.31 is handled well by the Chandezon technique for the problems in which we are interested.

By defining the groove profile in piecewise continuous line segments, the groove profile in each layer is specified by the set of L+1 line segment end points, {(0,0), (Subscript[((Subscript[x, 1]), (Subscript[y, 1])),....,((Subscript[x, L]), (Subscript[y, L]))}, j]) for the (j^th) layer.  The thickness of the (j^th) layer at x = 0 is (Subscript[e, j]).  The (n^th) line segment for layer j is given by

(Subscript[a, lj])(x)=(((x-(Subscript[x, l-1]^j))/((Subscript[x, l]^j)-(Subscript[x, l-1]^j)))) ((Subscript[y, l]^j)-(Subscript[y, l-1]^j))+(Subscript[y, l-1]^j).

The slope of this line segment is

((∂(Subscript[a, lj])(x))/(∂x))=(Subscript[a, lj]^′)=((((Subscript[y, l]^j)-(Subscript[y, l-1]^j))/((Subscript[x, l]^j)-(Subscript[x, l-1]^j)))).

The change in vertical distance between line segments is, in general, a function of x with period d,

(Subscript[Δ, lj])(x)≡ (Subscript[a, l,j+1])(x) - (Subscript[a, l,j])(x).

If the films are conformal then (Subscript[Δ, j]) = 0.  We also define the derivative of this function,

(Subscript[Δ, lj]^′)= (∂/(∂x))((Subscript[a, l,j+1]) - (Subscript[a, lj]))=(Subscript[a, l,j+1]^′)-(Subscript[a, lj]^′).

Then

(Subscript[Δ, lj])(x)=(Subscript[Δ, lj]^′)x+[-(Subscript[Δ, lj]^′) ((Subscript[x, l-1]^(j+1))-(Subscript[x, l-1]^j))+((Subscript[y, l-1]^(j+1))-(Subscript[y, l-1]^j))]≡(Subscript[Δ, lj]^′)x + (Subscript[f, lj])

where we have explicitly exhibited the x dependence of (Subscript[Δ, lj]).
The function exp(i (Subscript[r, q]^j) (Subscript[Δ, j])) also has period d and may be expanded in a Fourier series,

exp[i (Subscript[r, q]^j) (Subscript[Δ, j])(x)]=(Underscript[(∑^∞), m=-∞])(Subscript[L, m]^j)((Subscript[r, q]))exp(i m K x)

where the grating vector is

K=((2 π)/d),

and

(∫_0^d)exp[i (Subscript[r, q]^j) (Subscript[Δ, j])(x)] exp(-i n K x)dx=(∫_0^d)(∑_m)(Subscript[L, m]^j)((Subscript[r, q])) exp(i m K x)exp(-i n K x)dx

=(∑_m)d (Subscript[L, m]^j)((Subscript[r, q])) (Subscript[δ, mn])=d (Subscript[L, n]^j)((Subscript[r, q]))

or

(Subscript[L, n]^j)((Subscript[r, q]))=(Piecewise[{{1
0, for (Subscript[Δ, j])=0and n = 0
for (Subscript[Δ, j])=0and n ≠ 0}, {(1/d)(∫_0^d)exp[i (Subscript[r, q]^j) (Subscript[Δ, j])(x)-i n K x] dx,      otherwise}}]).

Because (Subscript[Δ, j]) is segmented along the groove in our scheme, the integration in Eq. (10.13) is also a summation over each segment,

(Subscript[L, n]^j)((Subscript[r, q]))=(1/d)(∑_(l=1)^L)(Piecewise[{{(Subscript[x, l])-(Subscript[x, l-1]), for (Subscript[Δ, lj])=0and n = 0}, {(∫_(Subscript[x, l-1])^(Subscript[x, l]))exp[i (Subscript[r, q]^j) (Subscript[Δ, lj])(x)-i n K x] dx, otherwise}}]).

The integration in Eq. (10.14) is easily carried out with Eq. (10.8).
For future reference, we define two additional functions for each layer, (Subscript[C, j])(x) and (Subscript[D, j])(x), which are also constant within each line segment of the groove profile,

(Subscript[C, nj])=(1/(1+((((Subscript[a, nj]^′)))^2)))

and

(Subscript[D, nj])=((Subscript[a, nj]^′)/(1+((((Subscript[a, nj]^′)))^2))).

The functions (Subscript[C, j])(x) and (Subscript[D, j])(x) are periodic with the groove profile and can be expanded in Fourier series themselves.  Specifically,

(Subscript[C, j])(x)=(∑_(p=-N)^N)(Subscript[c, pj]) exp(i p K x)

and

(Subscript[D, j])(x)=(∑_(p=-N)^N)(Subscript[d, pj]) exp(i p K x)

in the limit N -> ∞, where

(Subscript[c, pj])=(∑_(n=1)^L)(Piecewise[{{(((Subscript[C, nj]) ((Subscript[x, n])-(Subscript[x, n-1])))/d), p = 0}, {((Subscript[iC, nj])/(2 π p))[exp(-i p K (Subscript[x, n]) )-exp(-i p K (Subscript[x, n-1]) )],      otherwise}}]).

and

(Subscript[d, pj])=(∑_(n=1)^L)(Piecewise[{{(((Subscript[D, nj]) ((Subscript[x, n])-(Subscript[x, n-1])))/d), p = 0}, {((Subscript[iD, nj])/(2 π p))[exp(-i p K (Subscript[x, n]) )-exp(-i p K (Subscript[x, n-1]) )],      otherwise}}]).

Eqs. (10.19) and (10.20) are easily derived from the Fourier integral.  For example, Eq. (10.11) follows from Eq. (10.10) by

(∫_0^d)(Subscript[C, j])(x) exp(-i q K x )dx=(∑_(p=-N)^N)(∫_0^d)(Subscript[c, pj]) exp(i p K x) exp(-i q K x)dx

=(∑_(p=-N)^N)(d (Subscript[δ, pq]) )(Subscript[c, pj])=d (Subscript[c, qj]).

By Eq. (10.13), (Subscript[C, j])(x) is piecewise constant, so Eq. (10.17) can be rewritten as

d (Subscript[c, qj])=(∑_(n=1)^L)(Subscript[C, nj]) (∫_(Subscript[x, n-1])^(Subscript[x, n]))exp(-i q K x )dx

or, with the substitution q->p we arrive at Eq. (10.19), and similarly for (Subscript[d, pj]).
To recap, in order to compute the Fourier coefficients (Subscript[c, pj]) and (Subscript[d, pj]) of the functions (Subscript[C, j]) and (Subscript[D, j]) we follow the following procedure.  First, for each layer j, we compute the slopes, (Subscript[a, nj]^′), for each line segment in the groove profile according to Eq. (10.5).  From Eqs. (10.15) and (10.16) we compute the coefficients (Subscript[C, nj]) and (Subscript[D, nj]) for each line segment n within layer j.  Finally, from Eqs. (10.19) and (10.20), the Fourier coefficients (Subscript[c, pj]) and (Subscript[d, pj]) of the functions  (Subscript[C, j]) and (Subscript[D, j]), respectively, are obtained.
We also note that for sinusoidal groove profiles, analytical expressions can be obtained for (Subscript[c, pj]), (Subscript[d, pj]), and (Subscript[L, n]^j)((Subscript[r, q])) as described by L. Li. [35]

10.5.2  Coordinate system transformation

The next step is to define the fields everywhere.  One of the clever aspects of Chandezon's technique is the use of coordinate transforms so that the groove profiles are planar surfaces in the new coordinate systems.  Another clever idea is to represent the reflected propagating (diffracted) fields in the upper dielectric medium as a Fourier plane wave expansion everywhere above the top interface, even inside the grooves where such an expansion is not physically correct, and to then define the total field in the upper medium as a sum of the propagating modes and evanescent modes, the latter of which are chosen to provide the total true field even within the grooves.
The incident field is a plane wave with wavevector (Subscript[k, i]) and polar angle of incidence θ in the xy plane as shown in Fig. 10.31,

(F^i)=exp((Subscript[ik, i])x sin θ - (Subscript[ik, i])y cos θ).

The reflected, plane wave propagating modes are given by the expansion

(F^r)=(∑_(n∈(Subscript[U, r])))(Subscript[B, n]^r)exp(i (Subscript[α, n]) x + i (Subscript[β, n])y).

We will assume that there are (Subscript[P, r]) of these modes, each with unknown amplitude (Subscript[B, n]^r).  The transmitted plane wave propagating modes are given by the expansion

(F^t)=(∑_(n∈(Subscript[U, t])))(Subscript[B, n]^t)exp(i (Subscript[α, n]) x - i (Subscript[β, n])y)

and we assume that there are (Subscript[P, t]) of these modes with unknown amplitudes (Subscript[B, n]^t).  The (Subscript[α, n]) are the wavevector components along the x-direction of the various diffracted orders,

(Subscript[α, n])= (Subscript[k, i]) sin θ+n K,

and the (Subscript[β, n]) are the wavevector components along the y-direction of the various diffracted orders,

(Subscript[β, n]^j)=(Sqrt[(Subscript[k, j]^2)-(Subscript[α, n]^2)]) = (Sqrt[(Subscript[ϵ, j])(Subscript[k, 0]^2)-(Subscript[α, n]^2)])        for      (Subscript[k, i]) ≥ |(Subscript[α, n])|, i.e., n∈ (Subscript[U, r,t]).

(Subscript[ϵ, j]) is the relative permittivity (dielectric constant) of layer j, and (Subscript[k, 0]) is the vacuum wavevector previously defined.  (Subscript[α, n]) is the same for each layer of the film stack as well as the upper dielectric and the substrate by virtue of Snell's law.  However, (Subscript[β, n]^j) is a function of the index of refraction and in general will be different for each layer, the upper dielectric, and the substate.  For that reason it is given the superscript j.  When n ∉ (Subscript[U, r,t]), (Subscript[k, i])< |(Subscript[α, n])|, then (Subscript[β, n]) is imaginary and the fields in the expansions, Eqs. (10.25) and (10.26), are evanescent.  When the total fields in the upper dielectric and/or the substrate are expressed by this expansion of real and evanescent orders, this is the well known Raleigh approximation for the diffracted fields.  Unfortunately, this expansion does not correctly reproduce the fields within the groove region and, therefore, fails when the grooves have any appreciable depth. [36]
In the Chandezon technique the evanescent fields in the upper dielectric and the substrate are computed in the transformed coordinate system to properly represent the total field in this region.  There is a separate coordinate system for each layer in the film stack, but for now we will not include the layer label j and consider a generic transformed coordinate system.  Following the standard notation in the literature, the axes of the transformed coordinate system are labelled v, u, and w where

v=x,                    (Subscript[u, j])=y-(Subscript[a, j])(x)                 and             w=z.

In the transformed coordinate system it is convenient to define two different sets of basis vectors (which are not all unit vectors) in terms of the standard Cartesian coordinate system basis vectors.  One set is represented by subscripts,

(Subscript[(OverVector[e]), 1])=(OverHat[i])+a'(OverHat[j ])             (Subscript[(OverVector[e]), 2])=(OverHat[j]),          and                (Subscript[(OverVector[e]), 3])=(OverHat[k])

and the other by superscripts,

((OverVector[e])^1)=(OverHat[i]) ,                ((OverVector[e])^2)=(OverHat[j])-a'(OverHat[i ]),     and               ((OverVector[e])^3)=(OverHat[k]).

The basis vector (Subscript[(OverHat[e]), 1]) is parallel to the surface of the groove profile and is, therefore, useful for computing the tangential components of the field vectors at the interfaces, while the ((OverVector[e])^2) basis vector is normal to the surface of the groove profile and is useful for computing the normal components of the field vectors at the interfaces.  There are several useful relationships between the two sets of basis vectors,

((OverVector[e])^i)·(Subscript[(OverVector[e]), j])=(Subscript[δ, ij]),       ((OverHat[e])^i)×(Subscript[(OverVector[e]), j])=(Subscript[e, ijk])(Subscript[(OverVector[e]), k])       and         (Subscript[(OverVector[e]), i])×(Subscript[(OverVector[e]), j])=(Subscript[e, ijk])((OverVector[e])^k)

where (Subscript[δ, ij]) is the Kronecker delta function and (Subscript[e, ijk]) is Levi-Civita symbol or permutation symbol defined by

(Subscript[e, ijk])=(Piecewise[{{     0                    , for i=j, i=k, or j=k                   
}, {+1                 , for (1,2,3), (2,3,1), or (3,1,2) 
}, {-1                 , for (1,3,2), (3,2,1), or (2,1,3).}}])

Any vector (OverVector[A]) can be written

(OverVector[A])=(Subscript[A, v])(Subscript[(OverVector[e]), 1])+(Subscript[A, u])(Subscript[(OverVector[e]), 2])+(Subscript[A, w])(Subscript[(OverVector[e]), 3])=(A^1)(Subscript[(OverVector[e]), 1])+(A^2)(Subscript[(OverVector[e]), 2])+(A^3)(Subscript[(OverVector[e]), 3]).

The tangential components of the vector (OverVector[A]) at an interface are

(Subscript[(OverHat[e]), 1])·(OverVector[A])=((Subscript[(OverVector[e]), 1])/(Sqrt[1+(((a'))^2)]))·(OverVector[A])=((A^1)/(Sqrt[1+(((a'))^2)]))

and

(Subscript[(OverHat[e]), 3])·(OverVector[A])=(A^3).

The normal component of the vector (OverVector[A]) at an interface is

((OverHat[e])^2)·(OverVector[A])=(((OverVector[e])^2)/(Sqrt[1+(((a'))^2)]))·(OverVector[A])=((Subscript[A, u])/(Sqrt[1+(((a'))^2)])).

The vector component (A^2) can be represented in terms of the subscript basis vectors by

(A^2)=((((OverVector[A])-(A^1)((OverVector[e])^1)-(A^3)((OverVector[e])^3))· (Subscript[(OverVector[e]), 2]))/(1+(((a'))^2)))=((((Subscript[A, v])(Subscript[(OverVector[e]), 1])+(Subscript[A, u])(Subscript[(OverVector[e]), 2])+(Subscript[A, w])(Subscript[(OverVector[e]), 3]))·((OverVector[e])^2))/(1+(((a'))^2)))-(((A^1)((OverVector[e])^1)· ((OverVector[e])^2))/(1+(((a'))^2)))

=((Subscript[A, u])/(1+(((a'))^2)))+((a'(A^1))/(1+(((a'))^2)))=C (Subscript[A, u]) + D (A^1)

where C and D are the functions previously defined in Eqs. (10.15) and (10.18).  Similarly,

(Subscript[A, v])=((((OverVector[A])-(Subscript[A, u])(Subscript[(OverVector[e]), 2])-(Subscript[A, w])(Subscript[(OverVector[e]), 3]))· (Subscript[(OverVector[e]), 1]))/(1+(((a'))^2)))=((((A^1)((OverVector[e])^1)+(A^2)((OverVector[e])^2)+(A^3)((OverVector[e])^3))·(Subscript[(OverVector[e]), 1]))/(1+(((a'))^2)))=((A^1)/(1+(((a'))^2)))-((a'(Subscript[A, u]))/(1+(((a'))^2)))=C (A^1)-D (Subscript[A, u])

and

(Subscript[A, w])=((OverVector[A])-(Subscript[A, v])(Subscript[(OverVector[e]), 1])-(Subscript[A, u])(Subscript[(OverVector[e]), 2]))· (Subscript[(OverVector[e]), 3])=((A^1)((OverVector[e])^1)+(A^2)((OverVector[e])^2)+(A^3)((OverVector[e])^3))·(Subscript[(OverVector[e]), 3])=(A^3).

The differential operator is

∇=((∂/(∂v)))((OverVector[e])^1)+((∂/(∂u)))((OverVector[e])^2)+((∂/(∂w)))((OverVector[e])^3)

so the curl of a vector is

∇×(OverVector[A])=(((∂(A^3))/(∂u))-((∂(A^2))/(∂w)))(Subscript[(OverVector[e]), 1])+(((∂(A^1))/(∂w))-((∂(A^3))/(∂v)))(Subscript[(OverVector[e]), 2])+(((∂(A^2))/(∂v))-((∂(A^1))/(∂u)))(Subscript[(OverVector[e]), 3]).

By combining Eq. (10.39) with Eq. (10.43), we can separate the component parts of the Maxwell-like equation,

∇×(OverVector[A])=i γ (OverVector[R])

into

(((∂(A^3))/(∂u))-((∂(A^2))/(∂w)))=((∂(A^3))/(∂u))-((∂(C (Subscript[A, u]) + D (A^1)))/(∂w))=i γ (Subscript[R, v])=i γ (C (R^1)-D (Subscript[R, u])) ,

(((∂(A^1))/(∂w))-((∂(A^3))/(∂v)))=i γ (Subscript[R, u]),

and

(((∂(A^2))/(∂v))-((∂(A^1))/(∂u)))=((∂(C (Subscript[A, u]) + D (A^1)))/(∂v))-((∂(A^1))/(∂u))=i γ (Subscript[R, w])=i γ (R^3).

10.5.3  Solution of Maxwell's equations within each layer

Maxwell's equations are

∇×(OverVector[E])=i (Subscript[γ, 1]) (OverVector[H])

and

∇×(OverVector[H])=-i (Subscript[γ, 2]) (OverVector[E])

where

(Subscript[γ, 1])=((ω (Subscript[Z, 0]))/c) (Subscript[μ, r])= (Subscript[k, 0])(Subscript[Z, 0]) (Subscript[μ, r]),

(Subscript[γ, 2])=((ω )/((Subscript[Z, 0])c)) (Subscript[ϵ, r]) = ((Subscript[k, 0])/(Subscript[Z, 0])) (Subscript[ϵ, r]),

and

(Subscript[Z, 0])≡(Sqrt[((Subscript[μ, 0])/(Subscript[ϵ, 0]))]).

(Subscript[μ, r]) and (Subscript[ϵ, r]) are the relative permeability and permittivity of the material in a given layer, respectively, and (Subscript[μ, 0]) and (Subscript[ϵ, 0]) are the permeability and permittivity of free space, respectively.
When the incident light is in the xy plane, Maxwell's equations can be reduced to three coupled equations for the TE and TM modes separately.  Specifically, we define new fields (OverVector[F]) and (OverVector[G]) where

(OverVector[F])≡(OverVector[E])          and         (OverVector[G])≡(Subscript[γ, 1]) (OverVector[H])     for the TE mode, and

(OverVector[F])≡(Subscript[Z, 0])(OverVector[H])          and         (OverVector[G])≡-(Subscript[Z, 0])(Subscript[γ, 2]) (OverVector[E])     for the TM mode.

Therefore, the F field has components along the w (z) axis, so (Subscript[F, v]), (Subscript[F, u]), (F^1), and (F^2) are all zero, while the G field has components in the uv (xy) plane of incidence, so (Subscript[G, w]) and (G^3) are zero.  Then Maxwell's equations for both modes can be expressed by two equations,

∇×(OverVector[F])=i (OverVector[G])

and

∇×(OverVector[G])=-i (Subscript[γ, 1])(Subscript[γ, 2]) (OverVector[F.])

Due to the symmetry of the problem, there is no w dependence of the fields.  From Eqs. (10.45) - (10.47), we have

((∂(F^3))/(∂u))=i (C (G^1)-D (Subscript[G, u])) ,

i((∂(F^3))/(∂v))=(Subscript[G, u]),

and

((∂(C (Subscript[G, u]) + D (G^1)))/(∂v))-((∂(G^1))/(∂u))=-i (Subscript[γ, 1])(Subscript[γ, 2]) (F^3).

We eliminate (Subscript[G, u]) from Eqs. (10.57) and (10.59) using Eq. (10.58), and then simplify the notation by letting (F^3)->F and (G^1)->G.  Maxwell's equations have now been reduced to two linear differential equations which define the fields within each layer of the film stack,

((∂F)/(∂u))=i C G+D((∂F)/(∂v))

and

((∂G)/(∂u))=i (Subscript[γ, 1])(Subscript[γ, 2]) F+(∂/(∂v))(i C((∂F)/(∂v)))+((∂(D G))/(∂v)).

To solve these equations, we recall that the C and D functions have been expanded in a Fourier series in Eqs. (10.17) and (10.18) in the variable v = x.  We rewrite those equations and their partial derivatives with respect to v.

(Subscript[C, j])(v)=(∑_p)(Subscript[c, pj]) exp(i p K v)

((∂(Subscript[C, j])(v))/(∂v))=(∑_p)i p K(Subscript[c, pj]) exp(i p K v)

(Subscript[D, j])(v)=(∑_p)(Subscript[d, pj]) exp(i p K v)

and

((∂(Subscript[D, j])(v))/(∂v))=(∑_p)i p K(Subscript[d, pj]) exp(i p K v).

The fields F and G and their partial derivatives are also expanded in Fourier series in v.  The layer variable j is explicitly exhibited.

(F^j)(u,v)=(∑_m)(Subscript[F, m]^j)(u) exp(i (Subscript[α, m]) v)

((∂(F^j)(u,v))/(∂v))=(∑_m)i (Subscript[α, m])(Subscript[F, m]^j)(u) exp(i (Subscript[α, m]) v)

(G^j)(u,v)=(∑_m)(Subscript[G, m]^j)(u) exp(i (Subscript[α, m]) v)

((∂(G^j)(u,v))/(∂v))=(∑_m)i (Subscript[α, m]) (Subscript[G, m]^j)(u) exp(i (Subscript[α, m]) v)

The Fourier coefficients (Subscript[F, m]) and (Subscript[G, m]) are the field amplitudes in the transformed coordinate system within each layer.  Entering Eqs. (10.62) - (10.69) into Eq. (10.60) gives

((∂(F^j)(u,v))/(∂u))=i (∑_(m,p))((Subscript[c, pj])(Subscript[G, m]^j)  +(Subscript[α, m]) (Subscript[d, pj]) (Subscript[F, m]^j)) exp(i (Subscript[α, m+p]) v)

since

(Subscript[α, m]) + p K = ((Subscript[k, i]) sin θ + m K) + p K = (Subscript[α, m+p]).

Substituting m -> m - p, we can eliminate the v dependence and obtain an equation for each Fourier component of the field F,

((∂(Subscript[F, m]^j)(u))/(∂u))=i (∑_p)((Subscript[c, pj])(Subscript[G, m-p]^j)  +(Subscript[α, m-p]) (Subscript[d, pj]) (Subscript[F, m-p]^j)).

Finally, we let p -> m - p,

((∂(Subscript[F, m]^j)(u))/(∂u))=i (∑_p)((Subscript[c, m-p,j])(Subscript[G, p]^j)  +(Subscript[α, p]) (Subscript[d, m-p,j]) (Subscript[F, p]^j)).

With Eqs. (10.64), (10.65), (10.68) and (10.69) we find

((∂(D G))/(∂v))=((∂D)/(∂v))G+D((∂G)/(∂v))

=i (∑_p)i p K(Subscript[d, pj]) exp(i p K v)·(∑_m)(Subscript[G, m]^j)(u) exp(i (Subscript[α, m]) v)+i (∑_p)(Subscript[d, pj]) exp(i p K v)·(∑_m)(Subscript[α, m]) (Subscript[G, m]^j)(u) exp(i (Subscript[α, m]) v)

=i(∑_(m,p))(p K + (Subscript[α, m]))(Subscript[d, pj]) (Subscript[G, m]^j)  exp(i (Subscript[α, m+p]) v)

=i(∑_(m,p))(Subscript[α, m+p])(Subscript[d, pj]) (Subscript[G, m]^j)  exp(i (Subscript[α, m+p]) v)

=i(∑_(m,p))(Subscript[α, m])(Subscript[d, pj]) (Subscript[G, m-p]^j)  exp(i (Subscript[α, m]) v)

where in the last step we have again made the substitution m -> m - p.
Also, Eqs. (10.62), (10.63), (10.66), and (10.67) give

(∂/(∂v))(i C((∂F)/(∂v)))=i ((∂C)/(∂v))·((∂F)/(∂v))+i C (((∂^2)F)/(∂(v^2)))

=i (∑_p)i p K(Subscript[c, pj]) exp(i p K v)·(∑_m)i (Subscript[α, m])(Subscript[F, m]^j)(u) exp(i (Subscript[α, m]) v) -

i (∑_p)(Subscript[c, pj]) exp(i p K v)·(∑_m)(Subscript[α, m]^2)(Subscript[F, m]^j)(u) exp(i (Subscript[α, m]) v)

=-i (∑_(m,p))(Subscript[α, m]) (p K + (Subscript[α, m]))(Subscript[c, pj]) (Subscript[F, m]^j) exp(i (Subscript[α, m+p]) v)

=-i (∑_(m,p))(Subscript[α, m]) (Subscript[α, m+p])(Subscript[c, pj]) (Subscript[F, m]^j) exp(i (Subscript[α, m+p]) v)

=-i (∑_(m,p))(Subscript[α, m-p]) (Subscript[α, m])(Subscript[c, pj]) (Subscript[F, m-p]^j) exp(i (Subscript[α, m]) v)

Using Eqs. (10.78) and (10.84) in Eq. (10.61) gives

((∂(G^j))/(∂u))=i ((∑_(m,p))[(Subscript[γ, 1]^j)(Subscript[γ, 2]^j)(Subscript[δ, mp])(Subscript[F, p]^j)]+(∑_(m,p))[-(Subscript[α, m-p]) (Subscript[α, m])(Subscript[c, pj]) (Subscript[F, m-p]^j)+(Subscript[α, m])(Subscript[d, pj]) (Subscript[G, m-p]^j)])exp(i (Subscript[α, m]) v).

In the second summation, we make the substitution p -> m - p,

((∂(G^j)(u,v))/(∂u))=i ((∑_(m,p))[(Subscript[γ, 1]^j)(Subscript[γ, 2]^j)(Subscript[δ, mp])(Subscript[F, p]^j)]+(∑_(m,p))[-(Subscript[α, p]) (Subscript[α, m])(Subscript[c, m-p,j]) (Subscript[F, p]^j)+(Subscript[α, m])(Subscript[d, m-p,j]) (Subscript[G, p]^j)])exp(i (Subscript[α, m]) v),

and then combine the summations and pick out the terms of the individual Fourier components of (G^j),

((∂(Subscript[G, m]^j)(u))/(∂u))=i (∑_p)[((Subscript[γ, 1]^j)(Subscript[γ, 2]^j)(Subscript[δ, mp])-(Subscript[α, p]) (Subscript[α, m])(Subscript[c, m-p,j]))(Subscript[F, p]^j) +(Subscript[α, m])(Subscript[d, m-p,j]) (Subscript[G, p]^j)].

Eqs. (10.73) and (10.87) are an infinite set of linear differential equations in the field components.  To satisfy these equations, the u dependence of the field components (Subscript[F, m]) and (Subscript[G, m]) must clearly be exponential.  These two sets of equations can be put into the form of a linear partial differential vector equation.  Specifically, we let the (Subscript[F, m]) and (Subscript[G, m]) field amplitudes of the various Fourier components in layer j, where -N ≤ m ≤ +N, be represented by the vector ((OverVector[ζ])^j) with 4N+2 components,

((OverVector[ζ])^j)(u)≡{(Subscript[F, -N]^j),....,(Subscript[F, +N]^j),(Subscript[G, -N]^j),....,(Subscript[G, +N]^j)}

where in principle N -> ∞, but in practice we truncate N at a sufficiently large number to achieve the desired accuracy in the calculation.  Eqs. (10.73) and (10.87) can be written as a matrix equation,

-i((∂((OverVector[ζ])^j))/(∂u))=(({{   ((OverTilde[A])^j)   ,   ((OverTilde[B])^j)  
}, {((OverTilde[C])^j),  ((OverTilde[D])^j)}}))((OverVector[ζ])^j)

where (OverTilde[A]), (OverTilde[B]), (OverTilde[C]), and (OverTilde[D]) are rank 2N+1 square matrices with components given by,

(Subscript[A, mn]^j)=(Subscript[α, n])(Subscript[d, m-n,j])

(Subscript[B, mn]^j)= (Subscript[c, m-n,j])

(Subscript[C, mn]^j)= (Subscript[γ, 1]^j)(Subscript[γ, 2]^j) (Subscript[δ, mn]) - (Subscript[α, m])(Subscript[α, n]) (Subscript[c, m-n,j])

and

(Subscript[D, mn]^j)= (Subscript[α, m]) (Subscript[d, m-n,j]).

The indices m and n range between -N and +N.
When we consider the boundary conditions at the interfaces between layers in the next section, we will find that it is easier work with modified fields within each layer, (OverHat[F]) and (OverHat[G.])  First of all, we define the new variable (Subscript[η, j]),

(Subscript[η, j]) ≡ (Piecewise[{{(1/(Subscript[μ, r])), for the TE mode
}, {(1/(Subscript[ϵ, r])),    for the TM mode.}}])

Then

(Subscript[(OverHat[F]), m])≡ (Subscript[F, m])

and

(Subscript[(OverHat[G]), m])≡(Subscript[η, j]) (Subscript[G, m]^1).

In terms of the new fields, the vector

((Overscript[ζ, (OverVector[^])])^j)(u)≡{(Subscript[(OverHat[F]), -N]),....,(Subscript[(OverHat[F]), +N]),(Subscript[(OverHat[G]), -N]),....,(Subscript[(OverHat[G]), +N])}

is a solution of the linear partial differential equation

-i((∂((Overscript[ζ, (OverVector[^])])^j))/(∂u))=(({{  ((OverTilde[A])^j)   ,   (((OverTilde[B])^j)/(Subscript[η, j]))  
}, {(Subscript[η, j])((OverTilde[C])^j),    ((OverTilde[D])^j)}}))((Overscript[ζ, (OverVector[^])])^j)≡(OverTilde[T]) ((Overscript[ζ, (OverVector[^])])^j).

(The caret on top of this vector is the standard notation used in the literature, but does not indicate that this is a unit vector, which it is not).  (OverTilde[T]) is a rank 4N+2 square matrix, and there are 4N+2 independent solutions to Eq. (10.98), indexed by q, of the form

((OverVector[(OverHat[ζ])])^j)(u)=(Subscript[(OverVector[(OverHat[ζ])]), q]^j) exp(i (Subscript[r, q]^j) u).

The vector amplitudes (Subscript[(OverVector[(OverHat[ζ])]), q]^j) of these independent solutions of the fields which can exist in each layer are fixed by matching boundary conditions.  Once the boundary conditions have been applied, the total field within each layer is a sum of these q solutions with appropriate weighting amplitudes.   Plugging Eq. (10.99) into Eq. (10.98) gives

-i((∂((OverVector[(OverHat[ζ])])^j)(u))/(∂u))=(Subscript[r, q]^j) (Subscript[(OverVector[(OverHat[ζ])]), q]^j) exp(i (Subscript[r, q]^j) u)=(Subscript[r, q]^j)((OverVector[(OverHat[ζ])])^j)=(OverTilde[(T^j)]) ((OverVector[(OverHat[ζ])])^j),

so we see that the constants (Subscript[r, q]) are simply the eigenvalues of the matrix (OverTilde[T]), and the (Subscript[(OverVector[(OverHat[ζ])]), q]^j) are vectors proportional to the eigenvectors of (OverTilde[T]).  The eigenvalues, (Subscript[r, q]), determine whether or not the modes corresponding to each vector (Subscript[(OverVector[(OverHat[ζ])]), q]) are evanescent or propagating.  If the former, then the (Subscript[r, q]) value determines if the field decays exponentially in the upward or downward direction.  If the latter, the (Subscript[r, q]) value determines if the field is propagating in the upward or downward direction.  There will be 2N+1 eigenvalues for which Re((Subscript[r, q]))>0 (some propagating and some evanescent modes), and 2N+1 eigenvalues for which Re((Subscript[r, q]))<0.  For propagating modes, Im((Subscript[r, q])) = 0.  If, in addition, Re((Subscript[r, q]))>0, then the modes are propagating in the upward direction (we assume an exp(-iωt) time dependence for the fields).  Likewise, Re((Subscript[r, q])) < 0 with Im((Subscript[r, q]))=0 corresponds to a downward propagating mode.  If Im((Subscript[r, q])) ≠0, then the modes either grow or decay exponentially in the vertical direction.  If Im((Subscript[r, q])) > 0, then the modes decay exponentially as u (or y) -> ∞.  If Im((Subscript[r, q])) < 0, then the modes decay exponentially as u (or y) -> -∞.
Let (OverHat[V]) be the complete set of normalized eigenvectors of (OverTilde[T]).  Then

(Subscript[(OverVector[(OverHat[ζ])]), q]^j)=(Subscript[b, q]^j)(Subscript[(OverHat[V]), q]^j)

where the (Subscript[b, q]^j) are the  4N+2 (as yet unknown) constants of proportionality.  The complete solution to Eq. (10.100) is

((OverVector[(OverHat[ζ])])^j)(u)=(∑_q)(Subscript[(OverVector[(OverHat[ζ])]), q]^j) exp(i (Subscript[r, q]^j) u)=(∑_q)(Subscript[b, q]^j)(Subscript[(OverHat[V]), q]^j)exp(i (Subscript[r, q]^j) u).

Each separate diffraction order within a layer is given by one component of the vector sum ((OverVector[(OverHat[ζ])])^j)(u) on the left hand side of this equation.  For example, if N = 2 and we are interested in the amplitude of +(1^st) diffracted order of the TE mode, then from Eqs. (10.55) and (10.95) - (10.97),

((Overscript[ζ, (OverVector[^])])^j)(u)≡{(Subscript[E, -2]),(Subscript[E, -1]),(Subscript[E, 0]),(Subscript[E, +1]),(Subscript[E, +2]),(Subscript[η, j])(Subscript[γ, 1])(Subscript[H, -2]),(Subscript[η, j])(Subscript[γ, 1])(Subscript[H, -1]),(Subscript[η, j])(Subscript[γ, 1])(Subscript[H, 0]),(Subscript[η, j])(Subscript[γ, 1])(Subscript[H, +1]),(Subscript[η, j])(Subscript[γ, 1])(Subscript[H, +2])}

so we need to solve for the fourth component of ((Overscript[ζ, (OverVector[^])])^j) which is itself a sum of the fourth components of the 4N+2 eigenvectors (Subscript[(OverHat[V]), q]^j) weighted by (Subscript[b, q]^j) exp(i (Subscript[r, q]^j) u).  For convenience we follow the procedure of Priest et al. [32] and multiply each term in the sum by the phase factor exp(-i (Subscript[r, q]^j) (Subscript[d, j])) and include the reciprocal factor in the values for (Subscript[b, q]^j) in order to fix the phase to zero at the bottom interface of layer j.
If we now define a new square matrix ((OverTilde[M])^j) of rank 4N+2 whose columns are composed of the eigenvectors (Subscript[(OverHat[V]), q]^j), and a new vector ((OverVector[b])^j) of the corresponding eigenvalues (Subscript[b, q]^j) then Eq. (10.102) can be written as a matrix equation,

((OverVector[(OverHat[ζ])])^j)(u)=((OverTilde[M])^j)((OverTilde[Φ])^j)(u) ((OverVector[b])^j)

or

(Subscript[(OverHat[ζ]), p]^j)(u)=(∑_q)(Subscript[M, pq]^j)exp[i (Subscript[r, q]^j) (u-(Subscript[d, j]))](Subscript[b, q]^j)

where ((OverTilde[Φ])^j) is a diagonal matrix incorporating the phase factors,

(Subscript[Φ, mq]^j)(u)=(Subscript[δ, mq]) exp[i (Subscript[r, q]) (u-(Subscript[d, j]))]=(Subscript[δ, mq]) exp(i (Subscript[r, q])(Subscript[e, j]) )exp(i (Subscript[r, q])(Subscript[Δ, j]) )≡(Subscript[Γ, mq]^j)exp(i (Subscript[r, q])(Subscript[Δ, j]) ).

The unknown (Subscript[b, q]^j) eigenvalue amplitudes are determined by matching boundary conditions at the interfaces between layers.  At the bottom interface of layer j, ((OverTilde[Φ])^j)((Subscript[d, j]))=(OverTilde[I]), the identity matrix.  Also, for the fields within the substrate (j = 0), the phase matrix at the interface with the bottom of the thin film stack is the identity matrix by definition.  For the top interface of layer j,

y=(Subscript[d, j+1]) + (Subscript[a, j+1])(x)

or

(Subscript[u, j]) =y-(Subscript[a, j])(x)= (Subscript[d, j+1]) + (Subscript[Δ, j])(x),

so

(Subscript[Φ, mq]^j)( (Subscript[d, j+1]) + (Subscript[Δ, j])(x))=(Subscript[δ, mq]) exp(i (Subscript[r, q]) (Subscript[e, j])) exp(i (Subscript[r, q]) (Subscript[Δ, j])(x)).

We note that in general (Subscript[Δ, j]) is a periodic function of x.  When the films are conformal, Δ = 0 and

(Subscript[Φ, mq]^j)( (Subscript[d, j+1]) + (Subscript[Δ, j]))=(Subscript[δ, mq]) exp(i (Subscript[r, q]) (Subscript[e, j])) .

10.5.4  Fields in the incident medium and substrate for conformal thin films

10.5.4.1  Incident field

The field amplitudes in Eq. (10.105) are defined in the transformed coordinate system.  The amplitude of the incident field is given by Eq. (10.24), the field for the reflected propagating waves within the incident medium by Eq. (10.25), and the field for the transmitted propagating waves within the substrate by Eq. (10.26), in terms of the original Cartesian coordinate system.  We now derive the form of  these fields in the transformed coordinate system.  The F component of the incident field is

(F^i)=exp(i (Subscript[α, 0])x) exp( - i (Subscript[β, 0]^(Q+1)) y)

=exp(i (Subscript[α, 0])v) exp( - i (Subscript[β, 0]^(Q+1)) u) exp(- i (Subscript[β, 0]^(Q+1))(Subscript[a, Q+1])(v))

where we substitute for y using Eq. (10.29).  The final exponential on the right hand side of Eq. (10.102) is periodic in v and can be expanded in a Fourier series.  Substituting the generic t for  (Subscript[β, 0]^(Q+1)),

exp(- i t (Subscript[a, Q+1])(v))=(∑_m)(Subscript[L, m]^(Q+1))(t) exp(i m K v)

where

(Subscript[L, m]^(Q+1))(t)=(Piecewise[{{1

0
, if a or t= 0, and m = 0

if a or t= 0, and m ≠ 0}, {(1/d)(∫_0^d)exp(- i t (Subscript[a, Q+1])(v))·exp(-i m K v)dv, otherwise}}]).

We note that this is the same function as that given in Eq. (10.14) except that (Subscript[a, Q+1])(v) is substituted for (Subscript[Δ, lj]).  The integral in Eq. (10.114) becomes a summation over the line segments in the groove profile.  Using Eq. (10.113) in Eq. (10.112) gives

(F^i)=(∑_m)(Subscript[L, m]^(Q+1))((Subscript[β, 0]^(Q+1))) exp(i (Subscript[α, m])v- i (Subscript[β, 0]^(Q+1)) u)=(Piecewise[{{exp(i (Subscript[α, 0]) ν), for (Subscript[β, 0]^(Q+1))=0
}, {(∑_m)(Subscript[F, m]^i)((Subscript[β, 0]^(Q+1))) exp(i (Subscript[α, m])v), otherwise}}])

in the transformed coordinate system, where

(Subscript[F, m]^i)((Subscript[β, 0]^(Q+1)))≡(Subscript[L, m]^(Q+1))((Subscript[β, 0]^(Q+1))) exp(- i (Subscript[β, 0]^(Q+1)) u).

The G component of the incident field is obtained from Eq. (10.55), remembering that the only nonzero component of F is (F^3),

i ((OverVector[G])^i)=∇×((OverVector[F])^i)=(((∂(F^i3))/(∂u)))(Subscript[(OverVector[e]), 1])-(((∂(F^i3))/(∂v)))(Subscript[(OverVector[e]), 2])

or, remembering that G = (G^1) and using Eqs. (10.38) and (10.43),

i (G^i)=[1+((((Subscript[a, Q+1]^′)))^2)](((∂(F^i3))/(∂u)))-(Subscript[a, Q+1]^′)(((∂(F^i3))/(∂v))).

From Eq. (10.123),

((∂(F^i))/(∂u))= - i (Subscript[β, 0]^(Q+1))[exp(i (Subscript[α, 0])v) exp( - i (Subscript[β, 0]^(Q+1)) u) exp(- i (Subscript[β, 0]^(Q+1))(Subscript[a, Q+1])(v))]=- i (Subscript[β, 0]^(Q+1))(F^i)

and

((∂(F^i))/(∂v))=i (Subscript[α, 0]) [exp(i (Subscript[α, 0])v) exp( - i (Subscript[β, 0]^(Q+1)) u) exp(- i (Subscript[β, 0]^(Q+1))(Subscript[a, Q+1])(v))] - i (Subscript[β, 0]^(Q+1))(Subscript[a, Q+1]^′)(v)[exp(i (Subscript[α, 0])v) exp( - i (Subscript[β, 0]^(Q+1)) u) exp(- i (Subscript[β, 0]^(Q+1))(Subscript[a, Q+1])(v))]

=[i (Subscript[α, 0])-i (Subscript[β, 0]^(Q+1))(Subscript[a, Q+1]^′)(v)](F^i).

Plugging Eqs. (10.127) and (10.129) into (10.126) gives

(G^i)=-[ (Subscript[β, 0]^(Q+1))+ (Subscript[α, 0])(Subscript[a, Q+1]^′)(v)](F^i).

From Eq. (10.121),

(∂/(∂v))exp(- i t (Subscript[a, Q+1])(v))=(∂/(∂v))(∑_m)(Subscript[L, m])(t) exp(i m K v)

or

- i t (Subscript[a, Q+1]^′)(v)exp(- i t (Subscript[a, Q+1])(v))=(∑_m)i m K(Subscript[L, m])(t) exp(i m K v).

so

(Subscript[a, Q+1]^′)(v)exp(- i t (Subscript[a, Q+1])(v))=(∑_m)-(((m K)/( t)))(Subscript[L, m])(t) exp(i m K v).

From Eqs. (10.120) and (10.133),

(Subscript[a, Q+1]^′)(v)(F^i)=[(Subscript[a, Q+1]^′)(v) exp(- i (Subscript[β, 0]^(Q+1))(Subscript[a, Q+1])(v))]exp(i (Subscript[α, 0])x) exp( - i (Subscript[β, 0]^(Q+1)) u)

=-(∑_m)(((m K)/(Subscript[β, 0]^(Q+1))))(Subscript[L, m])(t) exp(i (Subscript[α, m])x)exp( - i (Subscript[β, 0]^(Q+1)) u)

so we can substitute the expansions (10.123) and (10.135) into (10.130), giving

(G^i)=-(∑_m)(Subscript[L, m]^(Q+1))((Subscript[β, 0]^(Q+1))) [(Subscript[β, 0]^(Q+1))-((( (Subscript[α, 0])m K)/(Subscript[β, 0]^(Q+1))))]exp(i (Subscript[α, m])v- i (Subscript[β, 0]^(Q+1)) u)

=(∑_m)(Subscript[G, m]^i)((Subscript[β, 0]^(Q+1))) exp(i (Subscript[α, m])v)

where

(Subscript[G, m]^i)((Subscript[β, 0]^(Q+1)))≡(∑_m)(Subscript[L, m]^(Q+1))((Subscript[β, 0]^(Q+1))) [-(Subscript[β, 0]^(Q+1))+m K((( (Subscript[α, 0]))/(Subscript[β, 0]^(Q+1))))]exp(- i (Subscript[β, 0]^(Q+1)) u).

The total incident field has now been expressed in Eqs. (10.115) and (10.130) as a Fourier expansion in the grating period.  In vector form,

((OverVector[ζ])^i)=(({{(OverVector[l])'}, {(OverVector[l])''}}))

where the 2N+1 (OverHat[F]) components are taken from Eq. (10.116),

(Subscript[ζ, m]^′)=(Subscript[L, m]^(Q+1))(-(Subscript[β, 0]^(Q+1))) exp(- i (Subscript[β, 0]^(Q+1)) u)

and the 2N+1 (OverHat[G]) components are taken from Eq. (10.130),

(Subscript[ζ, m]^(′′))=(Subscript[η, j]) [-(Subscript[β, 0]^(Q+1))+m K((( (Subscript[α, 0]))/(Subscript[β, 0]^(Q+1))))](Subscript[L, m]^(Q+1))(-(Subscript[β, 0]^(Q+1)))exp(- i (Subscript[β, 0]^(Q+1)) u)

where the (Subscript[η, j]) factor comes from Eq. (10.96).  When matching boundary conditions at the top surface of the grating (u = 0), we note that the factor exp(- i (Subscript[β, 0]^(Q+1)) u) in Eqs. (10.132) and (10.133) is equal to 1.

10.5.4.2  Reflected propagating field

As previously stated, the reflected field consists of (Subscript[P, r]) upward propagating plane wave modes in the original Cartesian coordinate system.  In the transformed coordinate system, from Eqs. (10.25) and (10.29),

(F^r)=(∑_(n∈(Subscript[U, r])))(Subscript[B, n]^r)exp(i (Subscript[α, n]) v)exp(i (Subscript[β, n]^(Q+1))u)exp[i (Subscript[β, n]^(Q+1))(Subscript[a, Q+1])(v)].

Again, the final factor on the right hand side is periodic and can be expanded in a Fourier series using Eq. (10.113),

exp( i (Subscript[β, n]^(Q+1))(Subscript[a, Q+1])(v))=(∑_m)(Subscript[L, m]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i m K v)

so

(F^r)=(∑_(n∈(Subscript[U, r])))(Subscript[B, n]^r)(∑_m)(Subscript[L, m]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[α, n+m]) v+i (Subscript[β, n]^(Q+1))u).

Letting m -> m - n,

(F^r)=(∑_(n∈(Subscript[U, r])))(Subscript[B, n]^r)(∑_m)(Subscript[L, m-n]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[α, m]) v+i (Subscript[β, n]^(Q+1))u).

The (OverHat[G]) field components are also obtained from the analog to Eq. (10.122),

i (G^r)=[1+((((Subscript[a, Q+1]^′)))^2)](((∂(F^r3))/(∂u)))-(Subscript[a, Q+1]^′)(((∂(F^r3))/(∂v))).

where

((∂(Subscript[F, n]^r))/(∂u))=i (Subscript[β, n]^(Q+1))(Subscript[F, n]^r).

((∂(Subscript[F, n]^r))/(∂v))=i (Subscript[α, n])(Subscript[F, n]^r)+i (Subscript[β, n]^(Q+1))(Subscript[a, Q+1]^′)(v)(Subscript[F, n]^r).

so

(Subscript[G, n]^r)=((Subscript[β, n]^(Q+1))-(Subscript[α, n]) (Subscript[a, Q+1]^′))(Subscript[F, n]^r).

Taking the derivative of Eq. (10.143),

((1/(i (Subscript[β, n]^(Q+1)))))(∂/(∂v))exp(i (Subscript[β, n]^(Q+1))(Subscript[a, j])(v))=((1/(i (Subscript[β, n]^(Q+1)))))(∂/(∂v))(∑_m)(Subscript[L, m]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i m K v)

for (Subscript[β, n]^(Q+1))≠0.  So

(Subscript[a, Q+1]^′)(v)exp( i (Subscript[β, n]^(Q+1))(Subscript[a, j])(v))=(∑_m)((( m K)/( (Subscript[β, n]^(Q+1)))))(Subscript[L, m]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i m K v).

Combining Eqs. (10.143) and (10.151), the periodic function (Subscript[a, Q+1]^′)(v)(F^r) has the expansion

(Subscript[a, Q+1]^′)(v)(F^r)=(∑_(n∈(Subscript[U, r])))(Subscript[B, n]^r)exp(i (Subscript[α, n]) v)[(Subscript[a, Q+1]^′)(v)exp(i (Subscript[β, n]^(Q+1))(Subscript[a, Q+1])(v))]exp(i (Subscript[β, n]^(Q+1))u)

=(∑_(n∈(Subscript[U, r])))(Subscript[B, n]^r)(∑_m)((( m K)/( (Subscript[β, n]^(Q+1)))))(Subscript[L, m]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[α, n+m]) v)exp(i (Subscript[β, n]^(Q+1))u)

for (Subscript[β, n]^(Q+1))≠0 and letting m -> m - n,

=(∑_(n∈(Subscript[U, r])))(Subscript[B, n]^r)(∑_m)((( (m-n) K)/( (Subscript[β, n]^(Q+1)))))(Subscript[L, m-n]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[α, m]) v)exp(i (Subscript[β, n]^(Q+1))u).

Inserting Eqs. (10.145) and (10.154) into Eq. (10.149) gives

(Subscript[G, n]^r)=(Subscript[B, n]^r)(∑_m)[(Subscript[β, n]^(Q+1))-(m-n) K((( (Subscript[α, n]))/( (Subscript[β, n]^(Q+1)))))](Subscript[L, m-n]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[α, m]) v+i (Subscript[β, n]^(Q+1))u)

for (Subscript[β, n]^(Q+1))≠0.  It is entirely possible that a diffracted order will have a reflected angle of 0° for which (Subscript[β, n]^(Q+1))= 0.  Because we expect the field amplitude to change continuously as the diffracted angle scans through 0°, we can determine the value for (Subscript[G, n]^r) at (Subscript[β, n]^(Q+1))=0 by taking the limit,

(Subscript[G, n]^r)((Subscript[β, n]^(Q+1))=0)=(Underscript[lim, (Subscript[β, n]^(Q+1))->0])((Subscript[B, n]^r)(∑_m)[(Subscript[β, n]^(Q+1))-(m-n) K((( (Subscript[α, n]))/( (Subscript[β, n]^(Q+1)))))](Subscript[L, m-n]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[α, m]) v+i (Subscript[β, n]^(Q+1))u))

=(Underscript[lim, (Subscript[β, n]^(Q+1))->0])(-(Subscript[B, n]^r)(∑_m)[(m-n) K((( (Subscript[α, n]))/( (Subscript[β, n]^(Q+1)))))](Subscript[L, m-n]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[α, m]) v))

But for m ≠ 0,

(Underscript[lim, t->0])(Subscript[L, m]^(Q+1))(t) =(1/d)(∫_0^d)(1- i t (Subscript[a, Q+1])(v))·exp(-i m K v)dv

=-((i t)/d)(∫_0^d) (Subscript[a, Q+1])(v)·exp(-i m K v)dv

so Eq. (10.149) becomes

(Subscript[G, n]^r)((Subscript[β, n]^(Q+1))=0)=(Subscript[B, n]^r)(∑_m)[(m-n) K (Subscript[α, n])]((∫_0^d) ((i (Subscript[a, Q+1])(v))/d)·exp(-i (m-n) K v)dv) exp(i (Subscript[α, m]) v)

for (Subscript[β, n]^(Q+1))=0.  Of course, the integral becomes a summation of integrals for our piecewise continuous groove profile.  From Eq. (10.4),

(∫_0^d) (Subscript[a, Q+1])(v)·exp(-i (m-n) K v)dv=(∑_(l=1)^L)(∫_(Subscript[ν, l-1])^(Subscript[ν, l]))((Subscript[a, l,Q+1]^′)ν +(Subscript[b, l,Q+1]))exp(-i (m-n) K v)dv

where

(Subscript[b, l,Q+1])≡(Subscript[y, l-1]^(Q+1))-(((Subscript[ν, l-1]^(Q+1))((Subscript[y, l]^(Q+1))-(Subscript[y, l-1]^(Q+1))))/((Subscript[ν, l]^(Q+1))-(Subscript[ν, l-1]^(Q+1)))).

The total reflected propagating field can now be represented by a matrix equation,

((OverVector[ζ])^r)=(OverTilde[M])'((OverVector[B])^r)

where ((OverVector[ζ])^r) is a 4N+2 component field vector, (OverTilde[M])' is a matrix with 4N+2 rows and (Subscript[P, r]) columns, and ((OverVector[B])^r) is a vector with (Subscript[P, r]) components of the unknown amplitudes of the reflected propagating plane waves.  (OverTilde[M])' can be split into upper and lower matrices,

(OverTilde[M])'≡(({{(OverTilde[M])''}, {(OverTilde[M])'''}}))

where according to Eqs. (10.137) and (10.147),

(Subscript[M, mn]^(′′))=(Piecewise[{{(Subscript[L, m-n]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[β, n]^(Q+1))u),     for (Subscript[β, n]^(Q+1))≠0
}, {0,    for (Subscript[β, n]^(Q+1))=0}}])

and

(Subscript[M, mn]^(′′′))=(Piecewise[{{(Subscript[η, j])[(Subscript[β, n]^(Q+1))-(m-n) K((( (Subscript[α, n]))/( (Subscript[β, n]^(Q+1)))))](Subscript[L, m-n]^(Q+1))(+(Subscript[β, n]^(Q+1))) exp(i (Subscript[β, n]^(Q+1))u),      for (Subscript[β, n]^(Q+1))≠0
}, {(Subscript[η, j])[(m-n) K (Subscript[α, n])]((∫_0^d) ((i (Subscript[a, Q+1])(v))/d)·exp(-i (m-n) K v)dv) ,      for (Subscript[β, n]^(Q+1))=0}}]).

In Eqs. (10.157) and (10.158), the m index runs between -N and +N while the n index corresponds to the (Subscript[P, r]) modes that are propagating in the upward direction, i.e., for which

Im((Subscript[β, n]^(Q+1))) = 0             and              Re((Subscript[β, n]^(Q+1))) > 0.

Once again, at the top surface the factor exp(i (Subscript[β, n]^(Q+1))u) is equal to 1.

10.5.4.3  Reflected evanescent field

There are only (Subscript[P, r])  reflected propagating modes given by Eq. (10.155), but there are a total of 2N+1 modes with Re((Subscript[β, n]^(Q+1))) > 0 in the upper dielectric medium.  The remaining 2N+1-(Subscript[P, r]) modes which are not propagating modes must be evanescent modes.  The evanescent modes can be represented by the eigenvectors from Eq. (10.101),

((OverVector[ζ])^er)=(∑_(q∉U))(Subscript[b, q]^(Q+1))(Subscript[(OverTilde[V]), q]^(Q+1))exp(i (Subscript[r, q]^j) u)

such that

Im((Subscript[r, q]^(Q+1))) > 0.

In terms of the M matrix,

((OverVector[ζ])^er)=((OverTilde[M])^(Q+1))((OverTilde[Φ])^(Q+1))((Subscript[u, Q+1])) ((OverVector[b])^(Q+1))

where the ((OverHat[M])^(Q+1)) matrix in this case is a rectangular matrix with 4N+2 rows and only those 2N+1-(Subscript[P, r]) columns for which the eigenvector corresponds to one of the eigenvalues satisfying inequality (10.161).   The total field in the incident medium is thus

((OverVector[ζ])^R)=((OverVector[ζ])^i)+((OverVector[ζ])^r)+((OverVector[ζ])^er)

in which there are 2N+1 total unknown amplitudes.

Transmitted propagating field

By Eqs. (10.26) and (10.29) the transmitted propagating modes are given by the expansion

(F^t)=(∑_(n∈(Subscript[U, t])))(Subscript[B, n]^t)exp(i (Subscript[α, n]) v - i (Subscript[β, n]^0)(u+(Subscript[a, 1])(v))).

exp(- i (Subscript[β, n]^0) (Subscript[a, 1])(v))=(∑_m)(Subscript[L, m]^1)(-(Subscript[β, n]^0)) exp(i m K v)

(F^t)=(∑_(n∈(Subscript[U, t])))(Subscript[B, n]^t)exp(i (Subscript[α, n]) v) exp(-i (Subscript[β, n]^0)u)exp(-i (Subscript[β, n]^0)(Subscript[a, 1])(v))

=(∑_(n∈(Subscript[U, t])))(Subscript[B, n]^t) exp(-i (Subscript[β, n]^0)u)(∑_m)(Subscript[L, m]^1)(-(Subscript[β, n]^0)) exp(i (Subscript[α, n+m]) v)

=(∑_(n∈(Subscript[U, t])))(Subscript[B, n]^t) (∑_m)(Subscript[L, m-n]^1)(-(Subscript[β, n]^0)) exp(i (Subscript[α, m]) v-i (Subscript[β, n]^0)u).

The (OverHat[G]) field components are once again obtained from the analog to Eq. (10.122),

i (G^t)=[1+((((Subscript[a, 0]^′)))^2)](((∂(F^t3))/(∂u)))-(Subscript[a, 0]^′)(((∂(F^t3))/(∂v)))

where

((∂(Subscript[F, n]^t))/(∂u))=-i (Subscript[β, n]^0)(Subscript[F, n]^t).

((∂(Subscript[F, n]^t))/(∂v))=i (Subscript[α, n])(Subscript[F, n]^t)-i (Subscript[β, n]^0)(Subscript[a, 0]^′)(v)(Subscript[F, n]^t).

so

(Subscript[G, n]^t)=-[ (Subscript[β, n]^0)+(Subscript[α, n])(Subscript[a, 0]^′)(v)](Subscript[F, n]^t).

From Eq. (10.165),

(∂/(∂v))exp(- i (Subscript[β, n]^0) (Subscript[a, 1])(v))=(∂/(∂v))(∑_m)(Subscript[L, m]^1)(-(Subscript[β, n]^0)) exp(i m K v)

(Subscript[a, 1]^′)(v)exp(- i (Subscript[β, n]^0) (Subscript[a, 1])(v))=-(∑_m)((m K)/((Subscript[β, n]^0) ))(Subscript[L, m]^1)(-(Subscript[β, n]^0)) exp(i m K v)

and

(Subscript[a, 1]^′)(v)(F^t)=(∑_(n∈(Subscript[U, t])))(Subscript[B, n]^t)exp(i (Subscript[α, n]) v) [(Subscript[a, 1]^′)(v)exp(-i (Subscript[β, n]^0)u)]exp(-i (Subscript[β, n]^0)(Subscript[a, 1])(v))

=-(∑_(n∈(Subscript[U, t])))(Subscript[B, n]^t) (∑_m)((m K)/((Subscript[β, n]^0) ))(Subscript[L, m]^1)(-(Subscript[β, n]^0)) exp(i (Subscript[α, n+m]) v)exp(-i (Subscript[β, n]^0)(Subscript[a, 1])(v))

=-(∑_(n∈(Subscript[U, t])))(Subscript[B, n]^t) (∑_m)(((m-n) K)/((Subscript[β, n]^0) ))(Subscript[L, m-n]^1)(-(Subscript[β, n]^0)) exp(i (Subscript[α, m]) v)exp(-i (Subscript[β, n]^0)(Subscript[a, 1])(v))

Then substituting Eqs. (10.168) and (10.177) into Eq. (10.172) gives

(Subscript[G, n]^t)= (∑_m)[-(Subscript[β, n]^0)+(m-n) K((Subscript[α, n])/((Subscript[β, n]^0) ))](Subscript[L, m-n]^1)(-(Subscript[β, n]^0)) exp(i (Subscript[α, m]) v-i (Subscript[β, n]^0)u)(Subscript[B, n]^t).

The (Subscript[P, t]) transmitted propagating modes are represented by the matrix equation

((OverVector[ζ])^t)=(OverTilde[N])'((OverVector[B])^t)

where ((OverVector[ζ])^t) is a 4N+2 component transmitted field vector, ((OverVector[B])^t)is a vector of the (Subscript[P, t]) unknown transmitted field amplitudes, and (OverTilde[N])' is a rectangular matrix with 4N+2 rows and (Subscript[P, t]) columns which can be divided into an upper and lower matrix, each with 2N+1 rows,

(OverTilde[N])'=(({{(OverTilde[N])''}, {(OverTilde[N])'''}})),

defined according to Eqs. (10.168) and (10.178),

(Subscript[N, mn]^(′′))=(Subscript[L, m-n]^1)(-(Subscript[β, n]^0)) exp(-i (Subscript[β, n]^0)u)

and

(Subscript[N, mn]^(′′′))=(Subscript[η, j]) [-(Subscript[β, n]^0)+(m-n) K((Subscript[α, n])/((Subscript[β, n]^0) ))](Subscript[L, m-n]^1)(-(Subscript[β, n]^0))  exp(-i (Subscript[β, n]^0)u).

As with the corresponding equation for the reflected field, the m index in these matrices ranges from -N to +N, while the n index covers all (Subscript[P, t]) values for which

Im((Subscript[β, n]^0)) = 0           and            Re((Subscript[β, n]^0)) < 0.

At the bottom surface of the grating, the factor exp(-i (Subscript[β, n]^0)u) in Eqs. (10.181) and (10.182) is equal to 1.

10.5.4.4  Transmitted evanescent field

Finally, there are 2N+1-(Subscript[P, t]) evanescent modes in the substrate and we can use the previously derived field expression in the transformed coordinate system, Eq. (10.105), to express these components,

((OverVector[ζ])^et)=(Subscript[(OverTilde[M]), ev]^0) (Subscript[(OverVector[b]), ev]^0),

where ((OverVector[ζ])^et) is the evanescent field vector in the substrate with 4N+2 components, (Subscript[(OverTilde[M]), ev]^0) is a rectangular matrix with 4N+2 rows and 2N+1-(Subscript[P, t]) columns, and (Subscript[(OverVector[b]), ev]^0) is a vector with 2N+1-(Subscript[P, t]) unknown evanescent field amplitudes.  The columns in (Subscript[(OverTilde[M]), ev]^0) correspond to the columns in ((OverTilde[M])^0) which represent these evanescent modes, i.e., for which the eigenvalues satisfy the condition

Im((Subscript[r, q]^0)) < 0.

The total transmitted field is then

((OverVector[ζ])^T)=((OverVector[ζ])^t)+((OverVector[ζ])^et)=(OverTilde[N])'((OverVector[B])^t)+(Subscript[(OverTilde[M]), ev]^0) (Subscript[(OverVector[b]), ev]^0)≡(Subscript[(OverTilde[M]), s]) (({{((OverVector[B])^t)}, {(Subscript[(OverVector[b]), ev]^0)}}))

where we have defined the new square matrix (Subscript[(OverTilde[M]), s]) with rank 2N+1 that is obtained by combining (OverTilde[N])' and (Subscript[(OverTilde[M]), ev]^0).

10.5.5  Boundary conditions for conformal thin films

The unknown field amplitudes in (OverHat[F]) and (OverHat[G]) within each layer are determined by satisfying the boundary conditions.  We know that the tangential components of the (OverVector[E]) and (OverVector[H]) fields are continuous across each interface.  We now start with the simplest case, that for which all layers have the identical groove profile.  For TE polarization, the E field lies along the w (z) axis and is tangential to each interface.  The H field is orthogonal to the E field and (H^1) is the component of the H field which is tangential to the interface, because it is parallel to the (Subscript[(OverHat[e]), 1]) basis vector,

(Subscript[(OverVector[e]), 1])·(OverVector[H])=(H^1).

By Eq. (10.53),

(OverHat[F])≡ F =(Subscript[E, w])

and

(Subscript[(OverHat[G]), m])=(Subscript[η, j]) (Subscript[G, m]^1)= ((Subscript[G, m]^1)/(Subscript[μ, r]))=(((Subscript[γ, 1])(H^1))/(Subscript[μ, r]))=(Subscript[k, 0])(Subscript[Z, 0])(H^1).

Therefore, (OverHat[F]) and (OverHat[G]) as defined in Eqs. (10.95) and (10.96) are proportional to (Subscript[E, w]) and (H^1), respectively, and must also be continuous across each interface for TE polarization.
For TM polarization, the (OverVector[H]) field is along the w (z)-direction and is tangential to each interface.  The (E^1) field component is orthogonal to the (OverVector[H]) field and is tangential to the interface.  Using Eq. (10.54) and the following definitions for (OverHat[F]) and (OverHat[G]),

(OverHat[F])≡ F=(Subscript[Z, 0])(Subscript[H, w])

and

(OverHat[G])=(Subscript[η, j]) (Subscript[G, m]^1)=((-(Subscript[Z, 0])(Subscript[γ, 2]) (E^1))/(Subscript[ϵ, r]))=-(Subscript[k, 0]) (E^1),

we find that (OverHat[F]) and (OverHat[G]) are again both proportional to field components that are continuous across each interface.  Therefore, for the conformal multilayer film stack, the vector (OverVector[(OverHat[ζ])]) defined in Eq. (10.97) is continuous across each interface for both TE and TM polarizations.
 
The boundary condition at the top of layer j is

((Overscript[ζ, (OverVector[^])])^(j+1))(u=(Subscript[d, j+1]))=((Overscript[ζ, (OverVector[^])])^j)(u=(Subscript[d, j+1]))

By Eq. (10.105), and remembering that the phase factor is defined to be zero at the bottom interface of a given layer,

((OverTilde[M])^(j+1))((OverTilde[Φ])^(j+1))((Subscript[d, j+1])) ((OverVector[b])^(j+1))=((OverTilde[M])^(j+1)) ((OverVector[b])^(j+1))=((OverTilde[M])^j)((OverTilde[Φ])^j)((Subscript[d, j+1])) ((OverVector[b])^j).

10.5.6  Transfer matrix for conformal thin films

Eq. (10.193) provides the basic mechanism for generating the transfer matrix that connects the fields at the top surface of the films to the fields within the substrate.  For multiple thin film layers,

((OverTilde[M])^(Q+1)) ((OverVector[b])^(Q+1))=((OverTilde[M])^Q)((OverTilde[Φ])^Q)((Subscript[d, Q+1])) ((OverVector[b])^Q)=((OverTilde[M])^Q)((OverTilde[Φ])^Q) ((Subscript[d, Q+1]))(((((OverTilde[M])^Q)))^(-1)) ((OverTilde[M])^(Q-1))((OverTilde[Φ])^(Q-1))((Subscript[d, Q])) ((OverVector[b])^(Q-1))

=[((OverTilde[M])^Q)((OverTilde[Φ])^Q) (((((OverTilde[M])^Q)))^(-1))][ ((OverTilde[M])^(Q-1))((((OverTilde[Φ])^(Q-1))(((OverTilde[M])^(Q-1))))^(-1))]···[((OverTilde[M])^1)((OverTilde[Φ])^1) (((((OverTilde[M])^1)))^(-1))]  ((OverTilde[M])^0)((OverTilde[Φ])^0) ((OverVector[b])^0)

=[((OverTilde[H])^Q) ((OverTilde[H])^(Q-1))··· ((OverTilde[H])^1)] ((OverTilde[M])^0)((OverTilde[Φ])^0)((Subscript[d, 0])) ((OverVector[b])^0)=(OverTilde[H])(((OverTilde[M])^0) ((OverVector[b])^0))

where

((OverTilde[H])^j)≡((OverTilde[M])^j)[((OverTilde[Φ])^j) ((Subscript[d, j+1]))](((((OverTilde[M])^j)))^(-1))

and

(OverTilde[H])≡((OverTilde[H])^Q) ((OverTilde[H])^(Q-1))··· ((OverTilde[H])^1).

(OverTilde[H]) is the transfer matrix, a square matrix of rank 4N+2.  Eqs. (10.194)-(10.196) connect the unknown field amplitudes inside the substrate, ((OverVector[b])^0), to the unknown field amplitudes in the incident medium, ((OverVector[b])^(Q+1)), via the transfer matrix (OverTilde[H]).
((OverTilde[M])^(Q+1))((OverVector[b])^(Q+1)) in Eq. (10.194) is the total field in the incident medium at the top surface, which was derived in Eq. (10.163).  ((OverTilde[M])^0) ((OverVector[b])^0) is the total transmitted field at the bottom interface, which was derived in Eq. (10.186).  Therefore, we can rewrite Eqs. (10.194)-(10.196) as

((OverVector[ζ])^R)((Subscript[u, Q+1]))=(OverTilde[H])((OverVector[ζ])^T)((Subscript[u, 0])).

This matrix equation may be separated into 4N+2 linear equations with 4N+2 unknown amplitudes, and thus it provides the solution for all the unknown (b^0) and (b^(Q+1)) amplitudes.  Once these amplitudes are known, the amplitudes for the field within any other layer are found simply by truncating Eq. (10.194) at the appropriate layer.
It may not be imediately obvious, however, how to solve Eq. (10.199).  The procedure is as follows.  First we expand Eq. (10.199) into all of its components.

((OverVector[(OverHat[ζ])])^i)+M'((OverVector[B])^r)+((OverTilde[M])^(Q+1))((OverTilde[Φ])^(Q+1))((Subscript[u, Q+1])) ((OverVector[b])^(Q+1))=(OverTilde[H])((OverTilde[N])'((OverVector[B])^t)+((OverTilde[M])^0)((OverVector[b])^0))=((OverTilde[H])(OverTilde[N])')((OverVector[B])^t)+((OverTilde[H]) ((OverTilde[M])^0))((OverVector[b])^0).

Solving for ((OverVector[(OverHat[ζ])])^i) and remembering that the definition of the phase factor means that ((OverTilde[Φ])^(Q+1)) is the identity matrix,

-M'((OverVector[B])^r)-((OverTilde[M])^(Q+1)) ((OverVector[b])^(Q+1))+((OverTilde[H])(OverTilde[N])')((OverVector[B])^t)+((OverTilde[H]) ((OverTilde[M])^0))((OverVector[b])^0)=((OverVector[(OverHat[ζ])])^i).

Now, the matrices on the left hand side of this equation all have 4N+2 rows, and the total number of columns is also 4N+2, so they can be combined into a single square matrix, (OverTilde[Z]), with rank 4N+2,

(OverTilde[Z])=[({{-M',         -((OverTilde[M])^(Q+1)),         (OverTilde[H])(OverTilde[N])',        (OverTilde[H]) ((OverTilde[M])^0)}})].

We also define the new vector of the 4N+2 unknown amplitudes,

(OverVector[B])=(({{((OverVector[B])^r)}, {((OverVector[b])^(Q+1))}, {((OverVector[B])^t)}, {((OverVector[b])^0)}})).

Then

(OverTilde[Z])(OverVector[B])=((OverVector[(OverHat[ζ])])^i)

or

(OverVector[B])=((((OverTilde[Z])))^(-1))((OverVector[(OverHat[ζ])])^i).

The quantities on the right hand side of this equation are all known or can be computed from the groove profile and the optical properties of the thin film stack.  The vector on the left hand side includes all the unknown field amplitudes within the incident medium and the substrate.
A simple example makes this much more obvious.  Let us consider the case where N = 1, and (Subscript[P, r]) and (Subscript[P, t]) are both equal 1 as well.  Thus, there will be two evanescent modes (2N+1-P) in both incident and transmitted media.  Let ((OverTilde[M])^0) be so arranged that the column vectors in ((OverTilde[M])^0) corresponding to the evanescent modes satisfying inequality (10.185) are the last two columns of this 6 × 6 matrix, and let ((OverTilde[M])^(Q+1)) be arranged so that the column vectors satisfying inequality (10.161) are the second and third columns.  Then Eq. (10.199) can be written

(({{(Subscript[F, -1]^i)}, {(Subscript[F, 0]^i)}, {(Subscript[F, +1]^i)}, {(Subscript[G, -1]^i)  }, {(Subscript[G, 0]^i)}, {(Subscript[G, +1]^i)}}))+(({{(Subscript[M, 11]^″)}, {(Subscript[M, 21]^″)}, {(Subscript[M, 31]^″)}, {(Subscript[M, 41]^(′′′))}, {(Subscript[M, 51]^(′′′))}, {(Subscript[M, 61]^(′′′))}}))(Subscript[b, 1]^r)+(({{(Subscript[M, 15]^(Q+1)), (Subscript[M, 16]^(Q+1))}, {(Subscript[M, 25]^(Q+1)), (Subscript[M, 26]^(Q+1))}, {(Subscript[M, 35]^(Q+1)), (Subscript[M, 36]^(Q+1))}, {(Subscript[M, 45]^(Q+1)), (Subscript[M, 46]^(Q+1))}, {(Subscript[M, 55]^(Q+1)), (Subscript[M, 56]^(Q+1))}, {(Subscript[M, 65]^(Q+1)), (Subscript[M, 66]^(Q+1))}}))(({{(Subscript[b, 5]^(Q+1))}, {(Subscript[b, 6]^(Q+1))}}))=(({{(Subscript[H, 11]), (Subscript[H, 12]), (Subscript[H, 13]), (Subscript[H, 14]), (Subscript[H, 15]), (Subscript[H, 16])}, {(Subscript[H, 21]), (Subscript[H, 22]), (Subscript[H, 23]), (Subscript[H, 24]), (Subscript[H, 25]), (Subscript[H, 26])}, {(Subscript[H, 31]), (Subscript[H, 32]), (Subscript[H, 33]), (Subscript[H, 34]), (Subscript[H, 35]), (Subscript[H, 36])}, {(Subscript[H, 41]), (Subscript[H, 42]), (Subscript[H, 43]), (Subscript[H, 44]), (Subscript[H, 45]), (Subscript[H, 46])}, {(Subscript[H, 51]), (Subscript[H, 52]), (Subscript[H, 53]), (Subscript[H, 54]), (Subscript[H, 55]), (Subscript[H, 56])}, {(Subscript[H, 61]), (Subscript[H, 62]), (Subscript[H, 64]), (Subscript[H, 64]), (Subscript[H, 65]), (Subscript[H, 61])}}))((({{(Subscript[N, 11]^″)}, {(Subscript[N, 21]^″)}, {(Subscript[N, 31]^″)}, {(Subscript[N, 41]^(′′′))}, {(Subscript[N, 51]^(′′′))}, {(Subscript[N, 61]^(′′′))}}))(Subscript[b, 1]^t)+(({{(Subscript[M, 12]^0)+(Subscript[M, 13]^0)}, {(Subscript[M, 22]^0)+(Subscript[M, 23]^0)}, {(Subscript[M, 32]^0) +(Subscript[M, 33]^0)}, {(Subscript[M, 42]^0)+(Subscript[M, 43]^0)}, {(Subscript[M, 52]^0) +(Subscript[M, 53]^0)}, {(Subscript[M, 62]^0)+(Subscript[M, 64]^0)}}))(({{(Subscript[b, 2]^0)}, {(Subscript[b, 3]^0)}})))

(({{(Subscript[F, -1]^i)}, {(Subscript[F, 0]^i)}, {(Subscript[F, +1]^i)}, {(Subscript[G, -1]^i)  }, {(Subscript[G, 0]^i)}, {(Subscript[G, +1]^i)}}))+(({{(Subscript[M, 11]^″), (Subscript[M, 15]^(Q+1)), (Subscript[M, 16]^(Q+1))}, {(Subscript[M, 21]^″), (Subscript[M, 25]^(Q+1)), (Subscript[M, 26]^(Q+1))}, {(Subscript[M, 31]^″), (Subscript[M, 35]^(Q+1)), (Subscript[M, 36]^(Q+1))}, {(Subscript[M, 41]^(′′′)), (Subscript[M, 45]^(Q+1)), (Subscript[M, 46]^(Q+1))}, {(Subscript[M, 51]^(′′′)), (Subscript[M, 55]^(Q+1)), (Subscript[M, 56]^(Q+1))}, {(Subscript[M, 61]^(′′′)), (Subscript[M, 65]^(Q+1)), (Subscript[M, 66]^(Q+1))}}))(({{(Subscript[b, 1]^r)}, {(Subscript[b, 5]^(Q+1))}, {(Subscript[b, 6]^(Q+1))}}))=(({{(Subscript[H, 11]^′), (Subscript[H, 12]^′), (Subscript[H, 13]^′)}, {(Subscript[H, 21]^′), (Subscript[H, 22]^′), (Subscript[H, 23]^′)}, {(Subscript[H, 31]^′), (Subscript[H, 32]^′), (Subscript[H, 33]^′)}, {(Subscript[H, 41]^′), (Subscript[H, 42]^′), (Subscript[H, 43]^′)}, {(Subscript[H, 51]^′), (Subscript[H, 52]^′), (Subscript[H, 53]^′)}, {(Subscript[H, 61]^′), (Subscript[H, 62]^′), (Subscript[H, 64]^′)}}))(({{(Subscript[b, 1]^t)}, {(Subscript[b, 2]^0)}, {(Subscript[b, 3]^0)}}))

where the (OverTilde[H]) matrix multiplication has been performed to obtain the new matrix (OverTilde[H']) and the right hand side of the equation has been simplified into a single matrix product.  Then the matrices can be rearranged and combined,

(({{-(Subscript[M, 11]^″), -(Subscript[M, 15]^(Q+1)), -(Subscript[M, 16]^(Q+1)),   (Subscript[H, 11]^′), (Subscript[H, 12]^′), (Subscript[H, 13]^′)}, {-(Subscript[M, 21]^″), -(Subscript[M, 25]^(Q+1)), -(Subscript[M, 26]^(Q+1)),   (Subscript[H, 21]^′), (Subscript[H, 22]^′), (Subscript[H, 23]^′)}, {-(Subscript[M, 31]^″), -(Subscript[M, 35]^(Q+1)), -(Subscript[M, 36]^(Q+1)),   (Subscript[H, 31]^′), (Subscript[H, 32]^′), (Subscript[H, 33]^′)}, {-(Subscript[M, 41]^(′′′)), -(Subscript[M, 45]^(Q+1)), -(Subscript[M, 46]^(Q+1)),   (Subscript[H, 41]^′), (Subscript[H, 42]^′), (Subscript[H, 43]^′)}, {-(Subscript[M, 51]^(′′′)), -(Subscript[M, 55]^(Q+1)), -(Subscript[M, 56]^(Q+1)),   (Subscript[H, 51]^′), (Subscript[H, 52]^′), (Subscript[H, 53]^′)}, {-(Subscript[M, 61]^(′′′)), -(Subscript[M, 65]^(Q+1)), -(Subscript[M, 66]^(Q+1)),   (Subscript[H, 61]^′), (Subscript[H, 62]^′), (Subscript[H, 64]^′)}}))(({{(Subscript[b, 1]^r)}, {(Subscript[b, 5]^(Q+1))}, {(Subscript[b, 6]^(Q+1))}, {(Subscript[b, 1]^t)}, {(Subscript[b, 2]^0)}, {(Subscript[b, 3]^0)}}))=(({{(Subscript[F, -1]^i)}, {(Subscript[F, 0]^i)}, {(Subscript[F, +1]^i)}, {(Subscript[G, -1]^i)}, {(Subscript[G, 0]^i)}, {(Subscript[G, +1]^i)}})).

or

(({{(Subscript[b, 1]^r)}, {(Subscript[b, 5]^(Q+1))}, {(Subscript[b, 6]^(Q+1))}, {(Subscript[b, 1]^t)}, {(Subscript[b, 2]^0)}, {(Subscript[b, 3]^0)}}))=(((({{-(Subscript[M, 11]^″), -(Subscript[M, 15]^(Q+1)), -(Subscript[M, 16]^(Q+1)), (Subscript[H, 11]^′), (Subscript[H, 12]^′), (Subscript[H, 13]^′)}, {-(Subscript[M, 21]^″), -(Subscript[M, 25]^(Q+1)), -(Subscript[M, 26]^(Q+1)), (Subscript[H, 21]^′), (Subscript[H, 22]^′), (Subscript[H, 23]^′)}, {-(Subscript[M, 31]^″), -(Subscript[M, 35]^(Q+1)), -(Subscript[M, 36]^(Q+1)), (Subscript[H, 31]^′), (Subscript[H, 32]^′), (Subscript[H, 33]^′)}, {-(Subscript[M, 41]^(′′′)), -(Subscript[M, 45]^(Q+1)), -(Subscript[M, 46]^(Q+1)), (Subscript[H, 41]^′), (Subscript[H, 42]^′), (Subscript[H, 43]^′)}, {-(Subscript[M, 51]^(′′′)), -(Subscript[M, 55]^(Q+1)), -(Subscript[M, 56]^(Q+1)), (Subscript[H, 51]^′), (Subscript[H, 52]^′), (Subscript[H, 53]^′)}, {-(Subscript[M, 61]^(′′′)), -(Subscript[M, 65]^(Q+1)), -(Subscript[M, 66]^(Q+1)), (Subscript[H, 61]^′), (Subscript[H, 62]^′), (Subscript[H, 64]^′)}})))^(-1))(({{(Subscript[F, -1]^i)}, {(Subscript[F, 0]^i)}, {(Subscript[F, +1]^i)}, {(Subscript[G, -1]^i)}, {(Subscript[G, 0]^i)}, {(Subscript[G, +1]^i)}})).

10.5.7  Boundary conditions for nonconformal thin films

If the groove profiles are not conformal, as shown in Fig. 10.32, then the situation is more complicated.  The (Subscript[E, w]) and (Subscript[H, w]) field components along the w (z)-direction are continuous as before, so the boundary condition for (OverHat[F]) remains the same,

((OverHat[F])^(j+1))((Subscript[u, j+1])=(Subscript[d, j+1]))=((OverHat[F])^j)((Subscript[u, j])=(Subscript[d, j+1])+(Subscript[Δ, j])(ν)).

(Subscript[Δ, j])(ν) is a periodic function of ν.
The (G^1) field component in layer j is no longer tangential to the j+1 interface, so (Subscript[η, j]) (G^1) is no longer continuous across the interface.  (G^1) and (Subscript[G, u]) in the j'th layer must be appropriately combined to obtain the field component that is tangential to the j+1 interface.

Referring to Fig. 10.32, it is clear that the field component of G in the j'th layer which is tangential to the j+1 interface is

((OverHat[G])^j)=(Subscript[(OverHat[G]), t]^j) cos Θ + (Subscript[(OverHat[G]), n]^j) sin Θ.

where

(Subscript[(OverHat[G]), t]^j)=(Subscript[(OverHat[e]), 1])·(OverVector[(OverHat[G])])=((Subscript[(OverVector[e]), 1])/(Sqrt[1+((((Subscript[a, j]^′)))^2)]))·(OverVector[(OverHat[G])])=(((OverHat[G])^(1,j))/(Sqrt[1+((((Subscript[a, j]^′)))^2)]))

and

(Subscript[(OverHat[G]), n]^j)=((OverHat[e])^2)·(OverVector[(OverHat[G])])=(((OverVector[e])^2)/(Sqrt[1+((((Subscript[a, j]^′)))^2)]))·(OverVector[(OverHat[G])])=((Subscript[(OverHat[G]), u]^j)/(Sqrt[1+((((Subscript[a, j]^′)))^2)])).

Θ is the angle between the line segment of the j+1 interface and the j interface, which can be written in terms of the slopes of these line segments as follows.

Θ=(Subscript[Θ, 2]) - (Subscript[Θ, 1])

where

tan (Subscript[Θ, 2]) = (Subscript[a, j+1]^′)       and       tan (Subscript[Θ, 1]) = (Subscript[a, j])' .

With a little algebraic manipulation and making use of standard trigonometric identities, one finds

cos Θ = ((1+(Subscript[a, j])'·(Subscript[a, j+1]^′))/((Sqrt[1+((((Subscript[a, j]^′)))^2)]) (Sqrt[1+((((Subscript[a, j+1]^′)))^2)])))

and

sin Θ = (((Subscript[a, j+1]^′)-(Subscript[a, j]^′))/((Sqrt[1+((((Subscript[a, j]^′)))^2)]) (Sqrt[1+((((Subscript[a, j+1]^′)))^2)]))).

Combining Eqs. (10.219) - (10.221), (10.224) and (10.225) gives

((OverHat[G])^j)=(((1+(Subscript[a, j]^′)·(Subscript[a, j+1]^′)) ((OverHat[G])^(1,j)) +((Subscript[a, j+1]^′)-(Subscript[a, j]^′))(Subscript[(OverHat[G]), u]^j))/([1+((((Subscript[a, j]^′)))^2)] (Sqrt[1+((((Subscript[a, j+1]^′)))^2)]))) .

This field is tangential to the j+1 interface and is, therefore, equal to the tangential component of the field in the j+1 medium,

(Subscript[(OverHat[G]), t]^(j+1))=(((OverHat[G])^(1,j+1))/(Sqrt[1+((((Subscript[a, j+1]^′)))^2)]))=(((1+(Subscript[a, j]^′)·(Subscript[a, j+1]^′)) (((OverHat[G])^(1,))^j) +((Subscript[a, j+1]^′)-(Subscript[a, j]^′))(Subscript[(OverHat[G]), u]^j))/([1+((((Subscript[a, j]^′)))^2) ](Sqrt[(((1+(Subscript[a, j+1]^′)))^2)])))

or

((OverHat[G])^(j+1))=(((1+(Subscript[a, j]^′)·(Subscript[a, j+1]^′))  )/(1+((((Subscript[a, j]^′)))^2) ))((OverHat[G])^j)+((((Subscript[a, j+1]^′)-(Subscript[a, j]^′)))/(1+((((Subscript[a, j]^′)))^2)))(Subscript[(OverHat[G]), u]^j)

where we have dropped the "1" superscript as was done previously following Eq. (10.66).  This result can be further simplified by making use of the definitions in Eqs. (10.23) and (10.24) for (Subscript[C, j]) and (Subscript[D, j]), as well as Eq. (10.15) for (Subscript[Δ, j]^′).

(Subscript[Δ, j]^′) (Subscript[D, j]) + 1 = ((Subscript[a, j+1]^′)-(Subscript[a, j]^′))(((Subscript[a, j]^′)/(1+((((Subscript[a, j]^′)))^2)))) + 1=((1+(Subscript[a, j]^′)·(Subscript[a, j+1]^′))/(1+((((Subscript[a, j]^′)))^2) ))

and

((OverHat[G])^(j+1))=((Subscript[Δ, j]^′) (Subscript[D, j]) + 1) ((OverHat[G])^j)+(Subscript[Δ, j]^′)(Subscript[C, j])(Subscript[(OverHat[G]), u]^j)

within each of the n line segments of the groove profile.  We now want to eliminate (Subscript[(OverHat[G]), u]) from the boundary condition in Eq. (10.228) and express this equation in terms of the (OverHat[G]) and (OverHat[F]) fields only.  From Eq. (10.63) we have

((OverVector[e])^2)·(∇×(OverVector[F]))=((OverVector[e])^2)·i (OverVector[G])

or

((∂(F^1))/(∂w))-((∂(F^3))/(∂v))=i (Subscript[G, u]).

There is no field dependence along the w-direction, so the first term is zero.  Indeed, (F^3) is the only nonzero component of F, so we previously dropped the superscript "3."  So we can replace (Subscript[G, u]) in Eq. (10.230) by a function of F from Eq. (10.232),

((OverHat[G])^(j+1))=((Subscript[Δ, j]^′)(ν) (Subscript[D, j]) + 1) ((OverHat[G])^j)+i (Subscript[Δ, j]^′)(ν)(Subscript[η, j])(Subscript[C, j])((∂(OverHat[F]))/(∂v)).

(Subscript[Δ, j]^′)(ν) is also a periodic function of ν.
Eqs. (10.210) and (10.225) are the boundary conditions for the nonconformal multilayer film stack.  For the conformal multilayer film stack, Δ'=0, and Eq. (10.225) reduces to

((OverHat[G])^(j+1))= ((OverHat[G])^j)

as expected.
The alert reader at this point may be uncomfortable with the manner in which we have matched fields at the interface between nonconformal grooves.  We are essentially making a one-to-one matching of the Fourier components of the field amplitudes across the interface, albeit making sure that we are matching tangential components.  This approach is known as the Rayleigh approximation.  It has been shown that this approximation fails when the the groove depth becomes an appreciable fraction of a wavelength.  The reason for this failure is that the fields in the region above the grating can be expanded in a Fourier series, but the fields within the grooves cannot, due to the evanescent components.  Similarly, when the grooves are not conformal, the upper interface of layer j is a grating within the transformed coordinate system of layer j.  Therefore, to properly handle the field within layer j at the upper interface, we should really transform each field component within layer j into the coordinate system of layer j+1 and then apply the matching boundary conditions.  Essentially each diffracted component in the Fourier series expansion of the field within layer j must be expanded in another Fourier series at the upper interface of this layer to properly account for all the field amplitudes within the groove region at this interface.  L. Li was the first to point out this issue. [37]

10.5.8  Transfer matrix for nonconformal thin films

Now that the boundary conditions have been obtained, we can apply them explicitly to the fields on both sides of each interface.  We consider first the (OverHat[F]) field and then the (OverHat[G]) field.  It is convenient to expand Eq. (10.104) as

((OverVector[ζ])^j)((Subscript[u, j]) )=(({{((OverVector[(OverHat[F])])^j)}, {((OverVector[(OverHat[G])])^j)}}))=((OverTilde[M])^j)((OverTilde[Φ])^j)((Subscript[u, j]) ) ((OverVector[b])^j)=(({{((OverTilde[M])^j)(+)}, {((OverTilde[M])^j)(-)}}))((OverTilde[Φ])^j)((Subscript[u, j]) ) ((OverVector[b])^j)

where ((OverTilde[M])^j)(+) is the upper half of matrix ((OverTilde[M])^j), and ((OverTilde[M])^j)(-) is the lower half.  Both are rectangular matrices of dimension 2N+1 rows by 4N+2 columns.  So

((OverVector[(OverHat[F])])^j)((Subscript[u, j]) )=((OverTilde[M])^j)(+)((OverTilde[Φ])^j)((Subscript[u, j]) )((OverVector[b])^j)

or from Eq. (10.109),

(Subscript[(OverHat[F]), p]^j)((Subscript[u, j]) )=(∑_q)(Subscript[M, pq]^j)(+)exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[r, q]^j) (Subscript[Δ, j]))(Subscript[b, q]^j),

and

((OverVector[(OverHat[G])])^j)((Subscript[u, j]) )=((OverTilde[M])^j)(-)((OverTilde[Φ])^j)((Subscript[u, j]) ) ((OverVector[b])^j)

or

(Subscript[(OverHat[G]), p]^j)((Subscript[u, j]) )=(∑_q)(Subscript[M, pq]^j)(-)exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[r, q]^j) (Subscript[Δ, j]))(Subscript[b, q]^j).

Matching the (OverHat[F]) field at the interface using Eqs. (10.66) and (10.210),

((OverHat[F])^(j+1))((Subscript[u, j+1]),v)=(∑_m)(Subscript[(OverHat[F]), m]^(j+1))((Subscript[u, j+1])) exp(i (Subscript[α, m]) v)=(∑_m)(Subscript[(OverHat[F]), m]^j)((Subscript[u, j])) exp(i (Subscript[α, m]) v)=((OverHat[F])^j)((Subscript[u, j]),v).

where (Subscript[u, j+1])= (Subscript[d, j+1]) and (Subscript[u, j])= (Subscript[d, j+1]) + (Subscript[Δ, j]) (and (Subscript[Δ, j]) is defined for each line segment of the groove profile in Eq. (10.6) and is zero if the films are conformal).  Substituting Eq. (10.229) into this equation gives

(∑_(m,q))(Subscript[M, mq]^(j+1))(+) exp(i (Subscript[α, m]) v)(Subscript[b, q]^(j+1))=(∑_(m,q))(Subscript[M, mq]^j)(+) exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[r, q]^j) (Subscript[Δ, j]))exp(i (Subscript[α, m]) v)(Subscript[b, q]^j).

The factor exp(i (Subscript[r, q]^j) (Subscript[Δ, j])) on the right hand side of Eq. (10.233), not present with conformal films, is periodic with the grating and was previously expanded in a Fourier series in Eq. (10.9).  Substituting this into Eq. (10.233) gives

(∑_mq)(Subscript[M, mq]^(j+1))(+) exp(i (Subscript[α, m]) v)(Subscript[b, q]^(j+1))=(∑_(m,q))(Subscript[M, mq]^j)(+) (∑_n)(Subscript[L, n]^j)((Subscript[r, q]))exp(i n K v)exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, m]) v)(Subscript[b, q]^j)

=(∑_mnq)(Subscript[M, mq]^j)(+)(Subscript[L, n]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, m+n]) v)(Subscript[b, q]^j).

Replacing n with n - m,

=(∑_mnq)(Subscript[M, mq]^j)(+)(Subscript[L, n-m]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, n]) v)(Subscript[b, q]^j)

and then swapping m and n,

=(∑_mnq)(Subscript[M, nq]^j)(+)(Subscript[L, m-n]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, m]) v)(Subscript[b, q]^j).

We can now equate term-by-term on both sides of the equality,

(∑_q)(Subscript[M, mq]^(j+1))(+) (Subscript[b, q]^(j+1))=(∑_nq)(Subscript[M, nq]^j)(+)(Subscript[L, m-n]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))(Subscript[b, q]^j).

We can define the new rectangular matrix with 2N+1 rows and 4N+2 columns, ((OverTilde[X])^j)(+), where

(Subscript[X, nq]^j)(+)=(∑_m)(Subscript[M, mq]^j)(+)(Subscript[L, n-m]^j)((Subscript[r, q]^j)).

and the square diagonal phase matrix of rank 4N+2,

(Subscript[Γ, mq]^j)=exp(i (Subscript[r, q]^j) (Subscript[e, j])) (Subscript[δ, mq])

as defined in Eq. (10.106).  Then Eq. (10.238) can be written

((OverTilde[M])^(j+1))(+) ((OverVector[b])^(j+1))=((OverTilde[X])^j)(+)((OverTilde[Γ])^j) ((OverVector[b])^j).

This is the result for nonconformal films that corresponds to the top half of the matrix equation for conformal films, Eq. (10.193).  We see that the top half of matrix ((OverTilde[M])^j) on the right hand side of Eq. (10.193) has been replaced by ((OverTilde[X])^j)(+) in Eq. (10.241).  This makes sense since as seen in Eq. (10.21), the function (Subscript[L, n-m]^j) is zero when the films are conformal except for n - m = 0, when (Subscript[L, n-m]^j) = 1.
A similar relation must be obtained for the (OverHat[G]) field amplitudes.  Substituting Eqs. (10.17), (10.18), (10.67), and (10.68) into the boundary condition Eq. (10.225),

(∑_m)(Subscript[(OverHat[G]), m]^(j+1))(u) exp(i (Subscript[α, m]) v)=((Subscript[Δ, j]^′) (∑_p)(Subscript[d, pj]) exp(i p K x) + 1) (∑_m)(Subscript[(OverHat[G]), m]^j)(u) exp(i (Subscript[α, m]) v)+i (Subscript[Δ, j]^′)(Subscript[η, j])(∑_p)(Subscript[c, pj]) exp(i p K x)(∑_m)i (Subscript[α, m])(Subscript[F, m]^j)(u) exp(i (Subscript[α, m]) v).

Substituting Eqs. (10.237) and (10.239) into (10.250) then gives

(∑_mq)(Subscript[M, mq]^(j+1))(-) exp(i (Subscript[α, m]) v)(Subscript[b, q]^(j+1))=(Subscript[Δ, j]^′) (∑_mpq)(Subscript[d, pj]) exp(i p K x)(Subscript[M, mq]^j)(-)exp(i (Subscript[r, q]^j) (Subscript[e, j])) exp(i (Subscript[r, q]^j) (Subscript[Δ, j])) exp(i (Subscript[α, m]) v)(Subscript[b, q]^j) +(∑_mq)(Subscript[M, mq]^j)(-)exp(i (Subscript[r, q]^j) (Subscript[Δ, j])) exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, m]) v)(Subscript[b, q]^j) -(Subscript[Δ, j]^′)(Subscript[η, j])(∑_mpq)(Subscript[α, m])(Subscript[c, pj]) (Subscript[M, mq]^j)(+)exp(i p K x)exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[r, q]^j) (Subscript[Δ, j])) exp(i (Subscript[α, m]) v)(Subscript[b, q]^j).

Now, from Eq. (10.9),

exp[i (Subscript[r, q]^j) (Subscript[Δ, j])(v)]exp(i (Subscript[α, m])v)=(∑_s)(Subscript[L, s]^j)((Subscript[r, q]))exp(i (Subscript[α, m+s])v)

=(∑_s)(Subscript[L, s-m]^j)((Subscript[r, q]))exp(i (Subscript[α, s])v)

Differentiating both sides with respect to v, and substituting Eq. (10.245) back in for one term,

i (Subscript[r, q]^j) (Subscript[Δ, j]^′)exp[i (Subscript[r, q]^j) (Subscript[Δ, j])(v)]exp(i (Subscript[α, m])v)+i (Subscript[α, m])(∑_s)(Subscript[L, s-m]^j)((Subscript[r, q]))exp(i (Subscript[α, s])v)=(∑_s)i (Subscript[α, s])(Subscript[L, s-m]^j)((Subscript[r, q]))exp(i (Subscript[α, s])v)

or

(Subscript[Δ, j]^′)exp[i (Subscript[r, q]^j) (Subscript[Δ, j])(v)]exp(i (Subscript[α, m])v)=(1/(Subscript[r, q]^j))(∑_s)((Subscript[α, s])-(Subscript[α, m]))(Subscript[L, s-m]^j)((Subscript[r, q]))exp(i (Subscript[α, s])v)

=(1/(Subscript[r, q]^j))(∑_s)(s-m)K (Subscript[L, s-m]^j)((Subscript[r, q]))exp(i (Subscript[α, s])v).

We can use Eq. (10.248) to eliminate (Subscript[Δ, j]^′) in Eq. (10.243), and Eq. (10.9) to eliminate (Subscript[Δ, j]),

(∑_mq)(Subscript[M, mq]^(j+1))(-) exp(i (Subscript[α, m]) v)(Subscript[b, q]^(j+1))=(∑_p)(Subscript[d, pj])exp(i p K v)(Underscript[ ∑, mqs])(1/(Subscript[r, q]^j)) (Subscript[M, mq]^j)(-)(s-m)K (Subscript[L, s-m]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, s])v)(Subscript[b, q]^j) +(∑_mq)(Subscript[M, mq]^j)(-)exp(i (Subscript[r, q]^j) (Subscript[e, j])) (∑_n)(Subscript[L, n]^j)((Subscript[r, q])) exp(i (Subscript[α, m+n]) v)(Subscript[b, q]^j) -(Subscript[η, j])(∑_p)(Subscript[c, pj])exp(i p K v)(∑_mqs)((Subscript[α, m])/(Subscript[r, q]^j)) (Subscript[M, mq]^j)(+)(s-m)K (Subscript[L, s-m]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, s])v)(Subscript[b, q]^j).

In the first and third terms on the right hand side, let s -> n - p, and in the second term let n -> n - m,

=(∑_p)(Subscript[d, pj])exp(i p K v)(∑_mnq)(1/(Subscript[r, q]^j)) (Subscript[M, mq]^j)(-)(n-m-p)K (Subscript[L, n-m-p]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, n-p])v)(Subscript[b, q]^j) +(∑_mnq)(Subscript[M, mq]^j)(-)(Subscript[L, n-m]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, n]) v)(Subscript[b, q]^j) -(Subscript[η, j])(∑_p)(Subscript[c, pj]) exp(i p K v)(∑_mnq)((Subscript[α, m])/(Subscript[r, q]^j))(Subscript[M, mq]^j)(+)(n-m-p)K (Subscript[L, n-m-p]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, n-p])v)(Subscript[b, q]^j)

and swap m and n in all the terms on the right hand side, and combine the exponential factors,

=(∑_p)(Subscript[d, pj])(∑_mnq)(1/(Subscript[r, q]^j)) (Subscript[M, nq]^j)(-)(m-n-p)K (Subscript[L, m-n-p]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, m])v)(Subscript[b, q]^j) +(∑_mnq)(Subscript[M, nq]^j)(-) (Subscript[L, m-n]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, m]) v)(Subscript[b, q]^j) -(Subscript[η, j])(∑_p)(Subscript[c, pj]) (∑_mnq)((Subscript[α, n])/(Subscript[r, q]^j))(Subscript[M, nq]^j)(+)(m-n-p)K (Subscript[L, m-n-p]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))exp(i (Subscript[α, m])v)(Subscript[b, q]^j).

We can now set the exp(i (Subscript[α, m])v) Fourier components equal to each other on right and left sides term-by-term, and we again swap m and n in the terms on the right hand side,

(∑_q)(Subscript[M, mq]^(j+1))(-)(Subscript[b, q]^(j+1))=K(∑_pmq)((Subscript[d, pj])/(Subscript[r, q]^j)) (Subscript[M, mq]^j)(-)(n-m-p) (Subscript[L, n-m-p]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))(Subscript[b, q]^j) +(∑_q)(Subscript[M, mq]^j)(-) (Subscript[L, n-m]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))(Subscript[b, q]^j) -K (Subscript[η, j])(∑_pmq)(Subscript[c, pj]) ((Subscript[α, m])/(Subscript[r, q]^j))(Subscript[M, mq]^j)(+)(n-m-p)(Subscript[L, n-m-p]^j)((Subscript[r, q]))exp(i (Subscript[r, q]^j) (Subscript[e, j]))(Subscript[b, q]^j).

We now define the rectangular matrix ((OverTilde[X])^j)(-) with 2N+1 rows and 4N+2 columns,

(Subscript[X, nq]^j)(-)=(∑_m)(Subscript[M, mq]^j)(-)(Subscript[L, n-m]^j)((Subscript[r, q]^j))+(K/(Subscript[r, q]^j))(∑_mp)(n-m-p)((Subscript[d, pj])(Subscript[M, mq]^j)(-)-(Subscript[η, j]) (Subscript[α, m])(Subscript[c, pj])(Subscript[M, mq]^j)(+)) (Subscript[L, n-m-p]^j)((Subscript[r, q]^j)).

Eq. (10.252) simplifies to

((OverTilde[M])^(j+1))(-)((OverVector[b])^(j+1))=((OverTilde[X])^j)(-)(Subscript[(OverTilde[Γ]), q]^j)((Subscript[e, j]))((OverVector[b])^j).

This is the result for nonconformal films that corresponds to the bottom half of the matrix equation for conformal films, Eq. (10.253).  The bottom half of matrix ((OverTilde[M])^j) on the right hand side of Eq. (10.253) has been replaced by ((OverTilde[X])^j)(-) in Eq. (10.254).  Because the function (Subscript[L, n]^j) is zero when the films are conformal except for n  = 0, when (Subscript[L, n]^j) = 1, Eq. (10.253) is quickly seen to reduce to the conformal case.
Eqs. (10.241) and (10.254) can be combined into a single equation,

((OverTilde[M])^(j+1))((OverVector[b])^(j+1))=((OverTilde[X])^j)((OverTilde[Γ])^j)((OverVector[b])^j)

where (OverTilde[X]) is a square matrix of rank 4N+2,

((OverTilde[X])^j)=(({{((OverTilde[X])^j)(+)}, {((OverTilde[X])^j)(-)}}))

and 0 ≤ j ≤ Q.  Eq. (10.255) for nonconformal films is the counterpart to Eq. (10.193) for conformal films.
We can use Eq. (10.255) to generate the transfer matrix.

((OverTilde[M])^(Q+1))((OverVector[b])^(Q+1))=((OverTilde[X])^Q)((OverTilde[Γ])^Q)((OverVector[b])^Q)=((OverTilde[X])^Q)((((OverTilde[Γ])^Q)(((OverTilde[M])^Q) ))^(-1))((OverTilde[X])^(Q-1))((OverTilde[Γ])^(Q-1))((OverVector[b])^(Q-1))

=((OverTilde[X])^Q)((((OverTilde[Γ])^Q)(((OverTilde[M])^Q) ))^(-1))((OverTilde[X])^(Q-1))((OverTilde[Γ])^(Q-1))···(((((OverTilde[M])^2) ))^(-1))((OverTilde[X])^1)((((OverTilde[Γ])^1)(((OverTilde[M])^1) ))^(-1))((OverTilde[X])^0)((OverTilde[Γ])^0)((OverVector[b])^0).

We define the new matrices,

((OverTilde[H])^j)=((OverTilde[X])^j) ((OverTilde[Γ])^j) (((((OverTilde[M])^j) ))^(-1))

so Eq. (10.258) becomes

(OverTilde[H])(((OverTilde[X])^0)((OverVector[b])^0))=((OverTilde[H])^Q)((OverTilde[H])^(Q-1))···((OverTilde[H])^1)((OverTilde[X])^0)((OverVector[b])^0)=((OverTilde[M])^(Q+1))((OverVector[b])^(Q+1))

where

(OverTilde[H])≡((OverTilde[H])^Q)((OverTilde[H])^(Q-1))···((OverTilde[H])^2)(H^1).

Eq. (10.260) is the transfer matrix equation that expresses the field amplitudes inside the substrate, ((OverVector[b])^0), in terms of the field amplitudes in the incident medium, ((OverVector[b])^(Q+1)).  At the bottom interface, (Subscript[Δ, 0])=0, and

(Subscript[L, m-n]^0)((Subscript[r, q]))={({{1,    for  m=n}, {0,    for  m≠n}}).

Then

(Subscript[X, nq]^0)(+)=(∑_m)(Subscript[M, mq]^0)(+)(Subscript[L, m-n]^0)((Subscript[r, q]))=(Subscript[M, nq]^0)(+)

and

(Subscript[X, nq]^0)(-)=(∑_m)(Subscript[M, mq]^0)(-)(Subscript[L, m-n]^0)((Subscript[r, q]))+(K/(Subscript[r, q]^j))(∑_mp)(n-m-p)((Subscript[d, pj])(Subscript[M, mq]^0)(-)-(Subscript[η, j]) (Subscript[α, m])(Subscript[c, pj])(Subscript[M, mq]^0)(+)) (Subscript[L, n-m-p]^0)((Subscript[r, q]))

=(Subscript[M, nq]^0)(-).

Therefore, within the substrate,

((OverTilde[X])^0)((OverVector[b])^0)=((OverTilde[M])^0)((OverVector[b])^0),

which is the total field in the substrate.  The transfer matrix Eq. (10.260) is solved in exactly the same manner as was discussed for conformal thin films with the matrix definition in Eq. (10.259) rather than Eq. (10.197).
A simpler expression for Eq. (10.253) has been given by Li. [37]

10.5.9  Scattering matrix

Although we have now derived the complete technique for solving the vector diffraction problem of a multilayer thin film stack with nonformal gratings at each layer, one finds that these equations exhibit numerical instability when the films become too thick or the grooves too deep due to the exponential dependence of the evanescent field amplitudes and the finite precision of the machine computation.  It is relatively straightforward to convert the transfer matrix formalism, in which the fields on one side of an interface are related to the fields on the other side of the interface, into a scattering matrix formalism in which the fields incident upon the interface from above and below are related to the fields reflected from the interface into the upper and lower media.  It turns out that this approach eliminates the problems in loss of precision caused by the exponentials in the calculation.  We follow the procedure described by Cotter et al. [38]
With reference to Fig. 10.33 we can denote the 2N+1 field amplitudes for which Re((Subscript[r, q]))>0, corresponding to the "upward" propagating and evanescent fields, with a "+" subscript.  Likewise, the 2N+1 field amplitudes for which Re((Subscript[r, q]))<0, corresponding to the "downward" propagating and evanscent fields, are denoted with a "-" subscript.

Let us now define a scattering matrix ((OverTilde[S])^j) which connects the 2N+1 field amplitudes which are incident upon the bottom interface, (Subscript[(OverVector[b]), +]^0), and the 2N+1 field amplitudes which are incident upon j'th interface, (Subscript[(OverVector[b]), -]^j), to the 2N+1 field amplitudes which are reflected from the bottom interface, (Subscript[(OverVector[b]), -]^0), and the 2N+1 field amplitudes which are reflected from the j'th interface, (Subscript[(OverVector[b]), +]^j),

(({{(Subscript[(OverVector[b]), +]^j)}, {(Subscript[(OverVector[b]), -]^0)}}))=(({{(Subscript[(OverTilde[S]), 11]^j), (Subscript[(OverTilde[S]), 12]^j)}, {(Subscript[(OverTilde[S]), 21]^j), (Subscript[(OverTilde[S]), 22]^j)}}))(({{(Subscript[(OverVector[b]), +]^0)}, {(Subscript[(OverVector[b]), -]^j)}})).

The eigenvalues (Subscript[r, q]) associated with each column in the transfer matrix, Eq. (10.263), for the 2N+1 "upward" modes satisfy either

Im((Subscript[r, q]^j)) = 0             and              Re((Subscript[r, q]^j)) > 0

for the propagating modes, or

Im((Subscript[r, q]^(Q+1))) > 0

for the evanescent modes.  We rearrange columns in ((OverTilde[X])^j) so that the upward propagating modes are the left-most columns of the matrix.  The remaining columns in the left half of the transfer matrix are associated with the evanescent modes.  Similarly, there are 2N+1 "downward" modes.  The eigenvalues associated with the column vectors in the transfer matrix for the downward modes satisfy

Im((Subscript[r, q]^j)) = 0             and              Re((Subscript[r, q]^j)) < 0

for the propagating modes, or

Im((Subscript[r, q]^(Q+1))) < 0

for the evanescent modes.  The column vectors associated with the downward propagating modes are placed on the left side of the right half of the transfer matrix.  Finally, the remaining column vectors in the right side of the transfer matrix are associated with the downward evanescent modes.  We can then write Eq. (10.263) as

((OverTilde[M])^(j+1))(({{(Subscript[(OverVector[b]), +]^(j+1))}, {(Subscript[(OverVector[b]), -]^(j+1))}}))=((OverTilde[X])^j)((OverTilde[Γ])^j)(({{(Subscript[(OverVector[b]), +]^j)}, {(Subscript[(OverVector[b]), -]^j)}}))

where the matrix ((OverTilde[Φ])^j) has also been appropriately computed according to Eq. (10.248).  Solving for the amplitude vector on the right hand side of this equation,

(({{(Subscript[(OverVector[b]), +]^j)}, {(Subscript[(OverVector[b]), -]^j)}})) = (({{(Subscript[(OverTilde[W]), 11]^j), (Subscript[(OverTilde[W]), 12]^j)}, {(Subscript[(OverTilde[W]), 21]^j), (Subscript[(OverTilde[W]), 22]^j)}}))(({{(Subscript[(OverVector[b]), +]^(j+1))}, {(Subscript[(OverVector[b]), -]^(j+1))}}))

where we have defined the new transfer matrix,

((OverTilde[W])^j)≡(((((OverTilde[X])^j)((OverTilde[Γ])^j)))^(-1))((OverTilde[M])^(j+1))=(((((OverTilde[Γ])^j)))^(-1))(((((OverTilde[X])^j)))^(-1))((OverTilde[M])^(j+1)).

We also note that at the bottom interface with the substrate, by virtue of Eqs. (10.186) and (10.266), we can reduce Eq. (10.272) to

((OverTilde[M])^1)(({{(Subscript[(OverVector[b]), +]^1)}, {(Subscript[(OverVector[b]), -]^1)}}))=((OverTilde[M])^0)(({{(Subscript[(OverVector[b]), +]^0)}, {(Subscript[(OverVector[b]), -]^0)}})) =(Subscript[(OverTilde[M]), s])(({{(Subscript[(OverVector[b]), +]^0)}, {(Subscript[(OverVector[b]), -]^0)}}))

so

(({{(Subscript[(OverVector[b]), +]^0)}, {(Subscript[(OverVector[b]), -]^0)}})) = (({{(Subscript[(OverTilde[W]), 11]^0), (Subscript[(OverTilde[W]), 12]^0)}, {(Subscript[(OverTilde[W]), 21]^0), (Subscript[(OverTilde[W]), 22]^0)}}))(({{(Subscript[(OverVector[b]), +]^1)}, {(Subscript[(OverVector[b]), -]^1)}}))

where

((OverTilde[W])^0)≡((((Subscript[(OverTilde[M]), s])))^(-1))((OverTilde[M])^1).

We now want to find the components of the scattering matrix for the (((j+1))^st) interface, ((OverTilde[S])^(j+1)), in terms of ((OverTilde[S])^j) and ((OverTilde[W])^j).  Eqs. (10.275) and (10.276) can be expanded into the four equations

(Subscript[(OverVector[b]), +]^j)=(Subscript[(OverTilde[S]), 11]^j)(Subscript[(OverVector[b]), +]^0) + (Subscript[(OverTilde[S]), 12]^j) (Subscript[(OverVector[b]), -]^j)

(Subscript[(OverVector[b]), -]^0)=(Subscript[(OverTilde[S]), 21]^j)(Subscript[(OverVector[b]), +]^0) + (Subscript[(OverTilde[S]), 22]^j) (Subscript[(OverVector[b]), -]^j)

(Subscript[(OverVector[b]), +]^j)=(Subscript[(OverTilde[W]), 11]^j)(Subscript[(OverVector[b]), +]^(j+1)) + (Subscript[(OverTilde[W]), 12]^j) (Subscript[(OverVector[b]), -]^(j+1))

(Subscript[(OverVector[b]), -]^j)=(Subscript[(OverTilde[W]), 21]^j)(Subscript[(OverVector[b]), +]^(j+1)) + (Subscript[(OverTilde[W]), 22]^j) (Subscript[(OverVector[b]), -]^(j+1))

We eliminate the field amplitudes with the superscript j from these equations by equating Eqs. (10.278) and (10.280),

(Subscript[(OverTilde[S]), 11]^j)(Subscript[(OverVector[b]), +]^0) + (Subscript[(OverTilde[S]), 12]^j) (Subscript[(OverVector[b]), -]^j)=(Subscript[(OverTilde[W]), 11]^j)(Subscript[(OverVector[b]), +]^(j+1)) + (Subscript[(OverTilde[W]), 12]^j) (Subscript[(OverVector[b]), -]^(j+1))

and by substituting (Subscript[(OverVector[b]), -]^j) from Eq. (10.281) into Eqs. (10.279) and (10.282),

(Subscript[(OverVector[b]), -]^0)=(Subscript[(OverTilde[S]), 21]^j)(Subscript[(OverVector[b]), +]^0) + (Subscript[(OverTilde[S]), 22]^j) ((Subscript[(OverTilde[W]), 21]^j)(Subscript[(OverVector[b]), +]^(j+1)) + (Subscript[(OverTilde[W]), 22]^j) (Subscript[(OverVector[b]), -]^(j+1)))

(Subscript[(OverTilde[S]), 11]^j)(Subscript[(OverVector[b]), +]^0) + (Subscript[(OverTilde[S]), 12]^j) ((Subscript[(OverTilde[W]), 21]^j)(Subscript[(OverVector[b]), +]^(j+1)) + (Subscript[(OverTilde[W]), 22]^j) (Subscript[(OverVector[b]), -]^(j+1)))=(Subscript[(OverTilde[W]), 11]^j)(Subscript[(OverVector[b]), +]^(j+1)) + (Subscript[(OverTilde[W]), 12]^j) (Subscript[(OverVector[b]), -]^(j+1)).

Rearranging terms gives

(Subscript[(OverVector[b]), -]^0)=(Subscript[(OverTilde[S]), 21]^j)(Subscript[(OverVector[b]), +]^0) + ((Subscript[(OverTilde[S]), 22]^j)(Subscript[(OverTilde[W]), 21]^j))(Subscript[(OverVector[b]), +]^(j+1))+ ( (Subscript[(OverTilde[S]), 22]^j) (Subscript[(OverTilde[W]), 22]^j) )(Subscript[(OverVector[b]), -]^(j+1))

and

((Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)-(Subscript[(OverTilde[W]), 11]^j))(Subscript[(OverVector[b]), +]^(j+1))=((Subscript[(OverTilde[W]), 12]^j)-(Subscript[(OverTilde[S]), 12]^j) (Subscript[(OverTilde[W]), 22]^j)) (Subscript[(OverVector[b]), -]^(j+1))- (Subscript[(OverTilde[S]), 11]^j)(Subscript[(OverVector[b]), +]^0) .

This last equation is rewritten

(Subscript[(OverVector[b]), +]^(j+1))= [((((Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)-(Subscript[(OverTilde[W]), 11]^j)))^(-1))((Subscript[(OverTilde[W]), 12]^j)-(Subscript[(OverTilde[S]), 12]^j) (Subscript[(OverTilde[W]), 22]^j))] (Subscript[(OverVector[b]), -]^(j+1))-  [((((Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)-(Subscript[(OverTilde[W]), 11]^j)))^(-1))(Subscript[(OverTilde[S]), 11]^j)](Subscript[(OverVector[b]), +]^0) .

Now we can eliminate the (Subscript[(OverVector[b]), +]^(j+1))term in Eq. (10.285) using Eq. (10.287),

(Subscript[(OverVector[b]), -]^0)=[(Subscript[(OverTilde[S]), 21]^j)-((Subscript[(OverTilde[S]), 22]^j)(Subscript[(OverTilde[W]), 21]^j))((((Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)-(Subscript[(OverTilde[W]), 11]^j)))^(-1))(Subscript[(OverTilde[S]), 11]^j)](Subscript[(OverVector[b]), +]^0) + [((Subscript[(OverTilde[S]), 22]^j)(Subscript[(OverTilde[W]), 21]^j))((((Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)-(Subscript[(OverTilde[W]), 11]^j)))^(-1))((Subscript[(OverTilde[W]), 12]^j)-(Subscript[(OverTilde[S]), 12]^j) (Subscript[(OverTilde[W]), 22]^j))+( (Subscript[(OverTilde[S]), 22]^j) (Subscript[(OverTilde[W]), 22]^j) ) ](Subscript[(OverVector[b]), -]^(j+1))

Now, the scattering matrix for layer j+1 provides two additional matrix equations,

(Subscript[(OverVector[b]), +]^(j+1))=(Subscript[(OverTilde[S]), 11]^(j+1))(Subscript[(OverVector[b]), +]^0) + (Subscript[(OverTilde[S]), 12]^(j+1)) (Subscript[(OverVector[b]), -]^(j+1))

and

(Subscript[(OverVector[b]), -]^0)=(Subscript[(OverTilde[S]), 21]^(j+1))(Subscript[(OverVector[b]), +]^0) + (Subscript[(OverTilde[S]), 22]^(j+1)) (Subscript[(OverVector[b]), -]^(j+1)).

Comparison of Eq. (10.287) with Eq. (10.289) immediately gives

(Subscript[(OverTilde[S]), 11]^(j+1))=((((Subscript[(OverTilde[W]), 11]^j)-(Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)))^(-1))(Subscript[(OverTilde[S]), 11]^j)

and

(Subscript[(OverTilde[S]), 12]^(j+1))=((((Subscript[(OverTilde[W]), 11]^j)-(Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)))^(-1))((Subscript[(OverTilde[S]), 12]^j) (Subscript[(OverTilde[W]), 22]^j)-(Subscript[(OverTilde[W]), 12]^j)).

Comparison of Eq. (10.290) with Eq. (10.288), and making use of Eqs. (10.291) and (10.292) to simplify the notation, we find

(Subscript[(OverTilde[S]), 21]^(j+1))=(Subscript[(OverTilde[S]), 21]^j)-((Subscript[(OverTilde[S]), 22]^j)(Subscript[(OverTilde[W]), 21]^j))((((Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)-(Subscript[(OverTilde[W]), 11]^j)))^(-1))(Subscript[(OverTilde[S]), 11]^j)=(Subscript[(OverTilde[S]), 21]^j)+(Subscript[(OverTilde[S]), 22]^j)(Subscript[(OverTilde[W]), 21]^j)(Subscript[(OverTilde[S]), 11]^(j+1))

and

(Subscript[(OverTilde[S]), 22]^(j+1))=((Subscript[(OverTilde[S]), 22]^j)(Subscript[(OverTilde[W]), 21]^j))((((Subscript[(OverTilde[S]), 12]^j)(Subscript[(OverTilde[W]), 21]^j)-(Subscript[(OverTilde[W]), 11]^j)))^(-1))((Subscript[(OverTilde[W]), 12]^j)-(Subscript[(OverTilde[S]), 12]^j) (Subscript[(OverTilde[W]), 22]^j))+( (Subscript[(OverTilde[S]), 22]^j) (Subscript[(OverTilde[W]), 22]^j) )=(Subscript[(OverTilde[S]), 22]^j) (Subscript[(OverTilde[W]), 22]^j)+(Subscript[(OverTilde[S]), 22]^j)(Subscript[(OverTilde[W]), 21]^j)(Subscript[(OverTilde[S]), 12]^(j+1)).

Eqs. (10.291) - (10.294) determine the (((j+1))^st) scattering matrix in terms of the (j^th) scattering and transfer matrices.  Starting with ((OverTilde[S])^0), the identity matrix, once the transfer matrices have been determined for each layer, we can successively determine the scattering matrix for each layer until we reach j= Q+1, the upper dielectric medium.  This final matrix connects the incident field amplitude (Subscript[(OverVector[b]), -]^(Q+1)) to the reflected field amplitudes, (Subscript[(OverVector[b]), +]^(Q+1)), and transmitted field amplitudes, (Subscript[(OverVector[b]), -]^0).

(({{(Subscript[(OverVector[b]), +]^(Q+1))}, {(Subscript[(OverVector[b]), -]^0)}}))=(({{(Subscript[(OverTilde[S]), 11]^(Q+1)), (Subscript[(OverTilde[S]), 12]^(Q+1))}, {(Subscript[(OverTilde[S]), 21]^(Q+1)), (Subscript[(OverTilde[S]), 22]^(Q+1))}}))(({{(Subscript[(OverVector[b]), +]^0)}, {(Subscript[(OverVector[b]), -]^(Q+1))}})).

We note that the amplitudes denoted by the vector (Subscript[(OverVector[b]), +]^0) are all identically zero - there are no incident fields from below the substrate onto the thin film stack.  Therefore, the reflected field amplitudes are

(Subscript[(OverVector[b]), +]^(Q+1))=(Subscript[(OverTilde[S]), 12]^(Q+1))(Subscript[(OverVector[b]), -]^(Q+1))

and the transmitted field amplitudes are

(Subscript[(OverVector[b]), -]^0)=(Subscript[(OverTilde[S]), 22]^(Q+1))(Subscript[(OverVector[b]), -]^(Q+1)).

In the paper by Priest et al., [32] there is an additional somewhat complex procedure described for converting these results into reflectivity and transmissivity amplitudes.  However, because we have already employed the procedure of Chandezon to rewrite the reflected and transmitted field amplitudes as sums of propagating plane waves and evanescent modes, Eqs. (10.296) and (10.297) are in fact the final result.  In particular, the first component of (Subscript[(OverVector[b]), -]^(Q+1)) is unity for the amplitude of the incident plane wave, and all other components are zero.  Therefore, the first column of (Subscript[(OverTilde[S]), 12]^(Q+1))is composed of the amplitude components for all the reflected orders, and the first column of (Subscript[(OverTilde[S]), 22]^(Q+1)) is composed of the amplitude components for all the transmitted orders.
When computing reflected and transmitted intensities for each diffracted order, it is useful to consider the component of the Poynting vector along the y-direction.  This component will be proportional to n·|E((|)^2) cos((Subscript[θ, m])) or |H((|)^2) cos((Subscript[θ, m]))/n where (Subscript[θ, m]) is the polar angle of the diffracted order and n is the refractive index of the medium.  For TE mode we are solving for the electric field amplitude, so we use the first of these functions to obtain the field intensity for both the diffracted order and the incident plane wave, while for TM mode we are solving for the magnetic field amplitude so we use the second of these functions to do the same thing.  We then normalize the diffracted intensity by the incident intensity to obtain the reflectivity and/or transmissivity of the grating film stack for each order.

The scattering matrices can also be used to find the field amplitudes within each layer.  Eqs. (10.289) and (10.290) with (Subscript[(OverVector[b]), +]^0) = 0 become

(Subscript[(OverVector[b]), +]^j)=(Subscript[(OverTilde[S]), 12]^j) (Subscript[(OverVector[b]), -]^j)

and

(Subscript[(OverVector[b]), -]^0)=(Subscript[(OverTilde[S]), 22]^j) (Subscript[(OverVector[b]), -]^j)             or            (Subscript[(OverVector[b]), -]^j)=((((Subscript[(OverTilde[S]), 22]^j)))^(-1)) (Subscript[(OverVector[b]), -]^0) .

Therefore, the total amplitude vector for layer j is a function of the scattering matrix for that layer and the transmitted amplitude vector,

((OverVector[b])^j)=(({{(Subscript[(OverTilde[S]), 12]^j) ((((Subscript[(OverTilde[S]), 22]^j)))^(-1)) 
}, {((((Subscript[(OverTilde[S]), 22]^j)))^(-1)) }}))·(({{(Subscript[(OverVector[b]), -]^0)}, {(Subscript[(OverVector[b]), -]^0)}})).

10.5.10  Field profile

Once the components of the (OverVector[b]) scattering amplitude vector have been computed for each layer, the field amplitude can be computed through the thickness of the film stack.  We consider the F component of the vector (OverVector[ζ]), which corresponds to the (OverVector[E]) field for TE polarization and to (Subscript[γ, 1])(OverVector[H]) for TM polarization.  For some arbitrary position on the x axis we can make use of the expansion in Eq. (10.66) to determine the field profile along the y-direction.  However, for simplicity we will consider the field profile at x = 0.  In that case,

(F^j)(u)=(∑_m)(Subscript[F, m]^j)(u)

where the different components (Subscript[F, m]) are given by the top half of the (OverVector[ζ]) vector.  The (OverVector[ζ]) vector is itself a summation of eigenvectors according to Eq. (10.101),

((OverVector[ζ])^j)=(∑_q)(Subscript[b, q]^j) (Subscript[(OverHat[V]), q]^j) (e^(i (Subscript[r, q]^j)( u-(Subscript[d, j]))))

so, in a given layer,

(({{(Subscript[F, -m])}, {(Subscript[F, -m+1])}, {:}, {(Subscript[F, m-1])}, {(Subscript[F, m])}}))=(({{(Subscript[b, 1])(Subscript[M, 11]) (e^(i (Subscript[r, 1]) u)) + (Subscript[b, 2])(Subscript[M, 12]) (e^(i (Subscript[r, 2]) u)) + .... + (Subscript[b, 4N+2])(Subscript[M, 1,4N+2]) (e^(i (Subscript[r, 4N+2]) u))}, {(Subscript[b, 1])(Subscript[M, 21]) (e^(i (Subscript[r, 1]) u)) + (Subscript[b, 2])(Subscript[M, 22]) (e^(i (Subscript[r, 2]) u)) + .... + (Subscript[b, 4N+2])(Subscript[M, 2,4N+2]) (e^(i (Subscript[r, 4N+2]) u))}, {:}, {(Subscript[b, 1])(Subscript[M, 2N,1]) (e^(i (Subscript[r, 1]) u)) + (Subscript[b, 2])(Subscript[M, 2N,2]) (e^(i (Subscript[r, 2]) u)) + .... + (Subscript[b, 4N+2])(Subscript[M, 2N,4N+2]) (e^(i (Subscript[r, 4N+2]) u))}, {(Subscript[b, 1])(Subscript[M, 2N+1,1]) (e^(i (Subscript[r, 1]) u)) + (Subscript[b, 2N+1,])(Subscript[M, 12]) (e^(i (Subscript[r, 2]) u)) + .... + (Subscript[b, 4N+2])(Subscript[M, 2N+1,4N+2]) (e^(i (Subscript[r, 4N+2]) u))}})).

Combining Eq. (10.301) with Eq. (10.303) gives the total field in layer j at x = 0,

(F^j)=(Underscript[(∑^(4N+2)), n=1])(Subscript[b, n]^j) (e^(i (Subscript[r, n]^j) (u-(Subscript[d, j]))))(Underscript[ (∑^(2N+1)), m=1])(Subscript[M, mn]^j).

Exercises

1.  Recognizing that (Subscript[ϵ, m]^′) generally becomes more negative as the wavelength increases, what do you expect to happen to the SP wavevector with increasing wavelength? (see Eq. 10.2)  What is the limiting value of (Subscript[k, sp]) as (Subscript[ϵ, m]^′)->-∞?

2. For what thickness of aluminum film is the Kretschmann configuration optimized at a incident wavelength of 900 nm?

ex2

3. At a wavelength of 320 nm, the refractive index for silver is 0.765 + i 0.611.  Is the reflectivity dip at this wavelength in Fig. 10.8 due to a SP resonance?  At what silver thickness is this dip minimized?

ex3a

ex3b

4. Is the component of the wavevector of an incident plane wave which lies in the plane of the surface larger or smaller for a larger angle of incidence?  In light of the answer to question 1, if the angle of incidence is changed from 43° to 45° what should happen to the SP resonance wavelength in the Kretschmann configuration?  Replot Fig. 10.8 for this new angle of incidence to verify your answer.

ex4

5.  It was stated in the text that a low index dielectric like (Subscript[MgF, 2]) can be used in place of the air spacer in the Otto configuration.  For (Subscript[MgF, 2]) with a refractive index of 1.38, what thickness of the dielectric film would give the largest dip in reflectivity for gold and at what angle of incidence is the SP excited?

ex6

6.  Explore the wavelength dependence of the coupling efficiency for gold, silver, aluminum and copper over the wavelength range of 300 nm to 1200 nm.  What thickness of aluminum film gives optimal coupling?

ex6

References

[1] A. Otto.  Excitation of nonradiative surface plasma waves in silver by the method of frustrated total reflection.  Z. Phys. 216 (1968) 398.
[2] E. Kretschmann and H. Raether.  Radiative decay of nonradiative surface plasmons excited by light.  Z. Naturforsch. 23a (1968) 2135.
[3] H. Kano.  Excitation of surface plasmon polaritons by a focused laser beam.  In Near-Field Optics and Surface Plasmon Polaritons, ed. S. Kawata (New York:  Springer, 2001) 189.
[4] R. W. Wood.  On the remarkable case of uneven distribution of light in a diffraction grating spectrum.  Philos. Mag. 4 (1902) 396.
[5] R. W. Wood.  Anomalous diffraction gratings.  Phys. Rev. 48 (1935) 928.
[6] Lord Rayleigh.  Note on the remarkable case of diffraction spectra described by Prof. Wood.  Philos. Mag. 14 (1907) 60.
[7] U. Fano.  The theory of anomalous diffraction gratings and of quasi-stationary waves on metallic surfaces (Sommerfeld’s waves).  J. Opt. Soc. Am. 31 (1941) 213.
[8] A. D. Rakic, A. B. Djurišic, J. M. Elazar and M. L. Majewski.  Optical properties of metallic films for vertical-cavity optoelectronic devices.  Appl. Opt. 37 (1998) 5271.
[9] S. J. Elston, G. P. Bryan-Brown, T. W. Priest and J. R. Sambles.  Surface-resonance polarization conversion mediated by broken surface symmetry.  Phys. Rev. B  44 (1991) 3483.
[10] S. J. Elston, G. P. Bryan-Brown and J. R. Sambles.  Polarization conversion from diffraction gratings.  Phys. Rev. B 44 (1991) 6393.
[11] E. A. Stern.  Plasma radiation by rough surfaces.  Phys. Rev. Lett. 19 (1967) 1321.
[12] E. Kretschmann.  The angular dependence and the polarisation of light emitted by surface plasmons on metals due to roughness.  Opt. Commun. 5 (1972) 331
[13] E. Kretschmann.  Die bestimmung der oberflächenrauhigkeit dünner schichten durch messung der winkelabhängigkeit der streustrahlung von oberflächenplasmaschwingungen.  Opt. Commun. 10 (1974) 353. 
[14] H. Raether, Surface Plasmons on Smooth and Rough Surfaces and on Gratings (Berlin: Springer, Berlin, 1988). 
[15] K. Arya, Z. B. Su and J. L. Birman.  Localization of the surface plasmon polariton caused by random roughness and its role in surface-enhanced optical phenomena.  Phys. Rev. Lett. 54 (1985) 1559
[16] S. I. Bozhevolnyi, B. Vohnsen, I. I. Smolyaninov and A. V. Zayats.  Direct observation of surface polariton localizaton caused by surface roughness.  Opt. Commun. 117 (1995) 417.
[17] S. I. Bozhevolnyi.  Localization phenomena in elastic surface-polariton scattering caused by surface roughness.  Phys. Rev. B 54 (1996) 8177
[18] S. I. Bozhelvolnyi,  I. I. Smolyaninov and A. V. Zayats.  Near-field microscopy of surface-plasmon polaritons:  Localization and internal interface imaging.  Phys. Rev. B 51 (1995) 17916.
[19] G. I. Stegeman, R. F. Wallis and A. A. Maradudin.  Excitation of surface polaritons by end-fire coupling.  Opt. Lett. 8 (1983) 386.
[20] J. J. Burke, G. I. Stegeman and T. Tamir.  Surface-polariton-like waves guided by thin, lossy metal films.  Phys. Rev. B 33 (1986) 5186.
[21] R. Charbonneau, N. Lahoud and P. Berini.  Demonstration of integrated optics elements based on long-range surface plasmon polaritons.  Opt. Exp. 13 (2005) 977.
[22] Z. Sun and D. Zeng.  Coupling of surface plasmon waves in metal/dielectric gap waveguides and single interface waveguides.  J. Opt. Soc. Am. B 24 (2007) 2883.
[23] A. Degiron, S-Y. Cho, C. Harrison, N. M. Jokerst, C. Dellagiacoma, O. J. F. Martin and D. R. Smith.  Experimental comparison between conventional and hybrid long-range surface plasmon waveguide bends.  Phys. Rev. A 77 (2008) 021804.
[24] B. Hecht, H. Bielefeldt, L. Novotny, Y. Inouye, and D. W. Pohl.  Local excitation, scattering, and interference of surface plasmons.  Phys. Rev. Lett. 77 (1996) 1889.
[25] H. Ditlbacher, J. R. Krenn, G. Schider, A. Leitner and F. R. Aussenegg.  Two-dimensional optics with surface plasmon polaritons.  Appl. Phys. Lett. 81 (2002) 1762.
[26] S. Kim, Y. Lim, H. Kim, J. Park, and B. Lee.  Optical beam focusing by a single subwavelength metal slit surrounded by chirped dielectric surface gratings.  Appl. Phys. Lett. 92 (2008) 013103.
[27] Y. S. Jung, J. Wuenschell, T. Schmidt, and H. K. Kim.  Near- to far-field imaging of free-space and surface-bound waves emanating from a metal nanoslit.  Appl. Phys. Lett. 92 (2008) 023104.
[28] H. W. Kihm, K. G. Lee, D. S. Kim, J. H. Kang and Q-H. Park.  Control of surface plasmon generation efficiency by slit-width tuning.  Appl. Phys. Lett. 92 (2008) 051115.
[29] T. Xu, Y. Zhao, D. Gan, C. Wang, C. Du and X. Luo.  Directional excitation of surface plasmons with subwavelength slits.  Appl. Phys. Lett. 92 (2008) 101501.
[30] L. Yin, V. K. Vlasko-Vlasov, A. Rydh, J. Pearson, U. Welp, S.-H. Chang, S. K. Gray, G. C. Schatz, D. B. Brown and C. W. Kimball.  Surface plasmons at single nanoholes in Au films.  Appl. Phys. Lett. 85 (2004) 467.
[31]  J. Chandezon, M. T. Dupuis, G. Cornet and D. Mayster.  Multicoated gratings: a differential formalism applicable in the entire optical region.  J. Opt. Soc. Am. 72 (1982) 839.
[32] T. W. Priest, N. P. K. Cotter and J. R. Sambles.  Periodic multilayer gratings of arbitrary shape.  J. Opt. Soc. Am. A 12 (1995) 1740.
[33] N. P. K. Cotter, T. W. Priest and J. R. Sambles.  Scattering matrix approach to multilayer diffraction.  J. Opt. Soc. Am. A 12 (1995) 1097.
[34] J. P. Plumey, B. Guizal and J. Chandezon.  Coordinate transformation method as applied to asymmetric gratings with vertical facets.  J. Opt. Soc. Am. A 14 (1997) 610.
[35]  L.Li.  Multilayer-coated diffraction gratings: differential method of Chandezon et al. revisited.  J Opt. Soc. Am. A 11 (1994) 2816.
[36]  D. Maystre.  Rigorous vector theories of diffraction gratings.  In Progress in Optics, ed. E Wolf  (Amsterdam:  Elsevier, 1984) 21, 1.
[37] L. Li.  Periodic multilayer gratings of arbitrary shape:  Comment.  J. Opt. Soc. Am. A 13 (1996) 1475
[38] N. P. K. Cotter, T. W. Priest and J. R. Sambles.  Scattering-matrix approach to multilayer diffraction.  J. Opt. Soc. Am. A 12 (1995) 1097.