Fail to obtain reasonable results for a deep 2D dieletric resonant grating with GratingGSM

[Samson-Ming]

Hi Professor Shcherbakov, 

These days I have been trying to use the GratingGSM to simulate a deep 2D dieletric resonant grating. But unfortunatelly,  I find it seems quite strange that I can't get the resonable results even with hard endaver to tune the control parameters, nevertherless it works well with non-resonant subwavelength gratings with proper tuning.

Could you please help me find out where I may get wrong? Thank you so much for your kind instruction.

The parameters of the grating are listed below, and I also attach the program file for your reference:

lambda=1.55 um;
TE, Normal incident from substrate;
Px=Py=0.65 um;
Cylinder radius r0=0.2515 um;

Cylinder hegith H=0.8 um;
n_sub=1.45;
n_cyl=3.361;
n_sup=1;

And the reference R and T by RCWA (also veified by FEM) are
R=0.0823775905
T=0.9176244730
Besides, the grating has a resonant dip around 1.519 um

[Samson-Ming]

The schematic of the structure:

My main file is also pasted:

clc;
format long;
%%
% fill the method structure
gsm_method.no = [30, 30]; % number of Fourier harmonics
gsm_method.nsl = 100; % number of slices
gsm_method.eps_b = 1; % basis permittivity of the grating layer
gsm_method.restart = []; % GMRes parameter
gsm_method.maxit = 500; % maximum number of iterations for the GMRes
gsm_method.tol = 1e-6; % tolerance for the GMRes

NO = gsm_method.no(1) * gsm_method.no(2); % total number of harmonics
ic = (ceil(gsm_method.no(1)/2)-1) * ...
gsm_method.no(2)+ceil(gsm_method.no(2)/2); % central harmonic

% grating containing multiple cylindrical pitches on a period

ncyl = [1, 1]; % numbers of cylinders along X and Y dimensions
Ncyl = ncyl(1) * ncyl(2); % total number of cylinders
% coordinates of centers relative to the period:
% cx = -0.5 + 0.5/ncyl(1) + (1/ncyl(1))[0:(ncyl(1)-1)];
% cy = -0.5 + 0.5/ncyl(2) + (1/ncyl(2))[0:(ncyl(2)-1)];
cx=0;
cy=0;
[CX, CY] = meshgrid(cx, cy);

grating.layer_type = "2D grating";
grating.type = "cylinder";
grating.period = 0.65 * ncyl;
grating.depth = 0.8;
grating.radius = 0.2515;%0.15 + 0.02*(rand(1, Ncyl)-0.5); %Note this random!
grating.center_x = CX;
grating.center_y = CY;
grating.eps = 3.361^2 * ones(1, Ncyl); % permittivities of cylinders
grating.eps_m = 1; % permittivity of the medium surrounding the cylinders
grating.eps_sub = 1.45^2; % substrate permittivity
grating.eps_sup = 1; % superstrate permittivity

% incidence angles:

phi_inc = 0;
theta_inc = 1e-10; %use 1e-10 for theta=0
% incidence wavevector projections:
k_inc = [sin(pitheta_inc/180)cos(piphi_inc/180); sin(pitheta_inc/180)sin(piphi_inc/180)];
Vinc = zeros(2*NO, 2); % incident amplitude vector
Vinc(ic, 1) = 1; % TE plane wave incidence from the substrate side
%Vinc:1st column:substrate; 2nd: superstrate; 1st half of column: TE; 2nd:TM
%Vinc(NO+ic, 2) = 1; %TM plane wave incidence from the superstrate side

wavelength = 1.55;
wavenumber = 2*pi/wavelength;
kg = [wavelength/grating.period(1), wavelength/grating.period(2)];

incidence.wavelength = wavelength;
incidence.k_inc = k_inc;
incidence.Vinc = Vinc;

% calculate wavevector projections and moduli for the reciprocal lattice:

[kz, ~, kx, ky, kxy] = kxyz(gsm_method.no, k_inc, kg, [gsm_method.eps_b, gsm_method.eps_b]);
gsm_method.kx = kx;
gsm_method.ky = ky;
gsm_method.kxy = kxy;
gsm_method.kz = kz;
gsm_method.kh = wavenumber * grating.depth;
% centers of slices along vertical direction:
gsm_method.zs = -0.5gsm_method.kh + ((1 : gsm_method.nsl)-0.5)(gsm_method.kh / gsm_method.nsl);

[Samson-Ming]

Later I tried to change the parameter gsm_method.eps_b to [(3.361^2+1)/2] or [3.361^2pi.2515^2/.65^2+1*(1-pi*.2515^2/.65^2)], the results improve a little but still far from convergence.

I check the effect of each geometry parameter and material permittivity, and I find that when

                                                  (grating.eps-grating.eps_m)*grating.depth

is too large (threshold?), the GratingGSM tends to converge quite slow. Unfortunately, due to the limitted resources of my PC, the code would report "Out of memery" if I set control parameters larger than gsm_method.no = [40, 40] and gsm_method.nsl = 400 before I get a convergent result.

I'm not sure if I make a reasonable assumption, or do you have some good advice to handle it other than simply enlarge the memory? Thank you so much.

Xianshun

[Samson-Ming]

Dear Professor Shcherbakov,

Sure, thank you so much for your instruction.

Best regards

Xianshun

[aashcher]

Dear Professor Shcherbakov,

Sure, thank you so much for your instruction.

Best regards

Xianshun

Indeed the factorization rules are implemented. Did you run the method gsm_diffraction with bnorm=1 ?

[Samson-Ming]

Yeah, I ran the method with "[Vout, Veff, balance, rv] = gsm_diffraction(incidence, gsm_method, grating, 1, 0);", it seems that bnorm has already be set to 1 (gsm_diffraction(incidence, method, grating, bnorm, bini)).

[Samson-Ming]

Dear Professor Shcherbakov,

Do you have some more advice on this specific case? Thank you so much.

The best results I obtained with many endeavor with GratingGSM are
R{0,0} = 0.132342;
T{0,0} = 0.867658;

which are still a little far from the reference results of
R=0.0823775905;
T=0.9176244730.

Xianshun

[Samson-Ming]

Dear Professor Shcherbakov,
Sure, thank you so much for your instruction.
Best regards
Xianshun

Dear Professor Shcherbakov,

It seems that the factorization rules are not implemented in the code (feps_cyl.m & fnxy_cyl.m) where only the analytical
2D Fourier transformation is adopted.

Best
Xianshun

[Samson-Ming]

Hi Professor Shcherbakov,

I have tried hard with GSM for this specific example even with the S vector method equipped and NVM technique used, but still I failed to obtain the correct answer, the best I acquire is

R{0,0} = 0.132342;
T{0,0} = 0.867658;

with NO=60×60 and total number of slices=500 （5×100） considering the limitation of CPU and memory.

So I am thinking that if GSM is suitable for grating layer with effective relative thickness d_eff>1? Here it is d_eff=(3.361-1)*0.8/1.55≈1.22. I have tried to reduce the d_eff to less than 1 by decreasing n_g、n_m、d_g or increasing λ, and all can output consistent result.

I wonder if you could try this simple case to check what the problem is or push the limit. Best regards and looking forward to your sooner reply.