RCWA - simple 1D binary gratings

Hello,

I am trying to model a simple 1D binary, metallic grating using the RCWA method and am having trouble. I'm trying to use the formulation of Moharam (JOSA A, 12 pp. 1068). When I run my code with a dielectric grating, I recover the results presented in this work, but when I change my index of refraction to a complex value I get a reflections > 1, which is obviously incorrect. I'll post my matlab code below, hoping that there is an simple mistake that I am missing. Any help is appreciated. Thanks:



function [R, T] = binaryGratingDiffractionTE(nI, nGR, nRD, nII, Lambda, f, d, lambda, theta, N)

% Returns the reflection and transmission efficiency of a binary grating
% using the rigurous coupled wave analysis technique given in the paper by
% Moharam et. al. in the paper, "Formulation for stable and efficient
% implementation of rigorous coupled-wave analysis of binary gratings,"
% JOSAA 12, 1068 (1995).
%
% Inputs: nI - refractive index in incident plane
% nGR - refractive index in ridge
% nRD - refractive index in groove
% nII - refractive index in output plane
% Lambda - period of grating
% f - fraction of period occupied by ridge
% d - height of grating
% lambda - wavelength of incident light (in a vacuum)
% theta - angle down from normal (degrees)
% N - trunction order (odd)
% Outputs: R - reflection efficiencies
% T - transmission efficiencies
% [order, efficiency]

% Convert angle to radians
theta = theta * pi / 180;

% Set order array
if (mod(N,2) == 0)
N = N + 1;
display('Order not odd: Increased by 1 for calculation');
end
n = ((1-N):(N-1))';
m = (-(N-1)/2:(N-1)/2)';
m0 = (N-1)/2 + 1;

% Set up wave-vectors: Eq. 6, 7
k0 = 2 * pi / lambda;
kxi = k0 .* (nI * sin(theta) - m .* (lambda/Lambda));
kIzi = conj(k0 .* sqrt(nI^2 - (kxi/k0).^2));
kIIzi = conj(k0 .* sqrt(nII^2 - (kxi/k0).^2));

% Calculatet Fourier harmonics of dielectric in grating region
% Note: Matlab calculates sin(n*pi) = (+/-)3.6739E-16 for integer n instead
% of 0... could be a source of some error.
epsilon = (nRD^2-nGR^2) .* sin(pi * f .* n) ./ (pi .* n);
epsilon(N) = nRD^2 * f + nGR^2*(1-f);

% Set up the A matrix: Eq. 16
Kx = diag(kxi/k0);
E = zeros(N);
for j = 1:N
E(j, :) = epsilon(j:j+N-1);
end
E = fliplr(E);
A = Kx * Kx - E;

% Find eigenvalues and eigenvectors of A matrix
[W, Q] = eig(A, 'nobalance');
Q = sqrt(Q);
q = diag(Q);

% Set up matrix V
V = W * Q;

% Set up matrix equation for cp and cm: Eq. 19, 20, 22, and 23
YI = diag(kIzi/k0);
YII = diag(kIIzi/k0);
X = diag(exp(-k0*q*d));
L = [i*YI*W + V, (i*YI*W - V) * X; (V - i*YII*W)*X, -V - i*YII*W];
R = [zeros(m0-1, 1); i*(kIzi(m0)/k0 + nI * cos(theta)); zeros(N+m0-1, 1)];

% Solve equations
C = L\R;
% C = pinv(L) * R;

% Find Ri: Eq. 21
LR1 = [eye(N); -i*YI];
RR1 = [W, W * X; V, -V*X] * C - [zeros(m0-1, 1); 1; zeros(m0-1,1); zeros(m0-1,1); i*cos(theta)* nI; zeros(m0-1,1)];
RI = LR1\RR1;
DEri = RI .* conj(RI) .* real(kIzi / (k0 * nI * cos(theta)));
R = [m, DEri];

% Find Ti: Eq. 24
LT1 = [eye(N); i*YII];
RT1 = [W*X, W; V*X, -V] * C;
TI = LT1\RT1;
DEti = TI .* conj(TI) .* real(kIIzi / (k0 * nI * cos(theta)));
T = [m, DEti];

I am working on a similar problem (however, I am intersted in TM mode). Did you ever resolve the problem with using complex refractive indices?

You said you where able to reproduce results from the Moharam paper. What values did you use? I gives values for nII,nRD,nI, and theta i Fig 2. What values did you use for nGR? Any guidance you can give is greatly appreciated. Thanks.

Yeah... it has to do with the way matlab takes square roots; basically for a complex nII you need to get rid of the conj on kIIzi.

However, I found that the convergence for this is truly horrible. For a gold in air grating that is 20 nm, working with visible wavelengths, I needed to keep something like N=401 for convergence. That took way to long to make this usable.

Still, as long as you don't mind working with dielectrics...

could you post the actual matlab function file, some of the symbols are not present.

also, can you post an example of how you called this function using the complex permittivities? Not the results, as you said it takes a long time to converge, but just the input parameters :nI, nGR, nRD, nII, Lambda, f, d, lambda, theta, N

Thanks!

I think it should substitute nI^2 for nI*conj(nI).So the matrix E may be still symmetrical with real elements.

Did you solve your problem?

I modified your code to TM wave and it was working fine. But I tried to extend it to conical diffraction and the program was not producing the same figure as the paper. Have you done the conical diffraction code?