Login Logout Register Discuss Discuss Edit Edit
Jump to: navigation, search

Optical simulation of multilayer stack of two dimensional (2D) metallic photonic crystals

Problem Definition

Goal of this simulation model repository is to simulate far field optical scattering from a metalic two dimensional (2D) photonic crystal consisting of square metal nano particles arranged in a square lattice as shown in figure below. . Structure is also known as a 2D square grating. We are interested in simulating oblique (conical) incidence of unpolarized optical plane waves on the gratings. We want to know how much light is expected to be reflected, how much is expected to be transmitted through the entire structure and how much gets absorbed by the structure. We are also interested in the breakup of reflected and transmitted power components into different Bragg components. For example, we want to know how much light is reflected at different Bragg scattering angles.

Simulator Used

We use the Rigorous Electromagnetic Simulator (REMS) developed at Stanford University. REMS is a rigorous and general purpose simulator for simulating periodic or non-periodic stratified (multilayer) structures. Simulator is technical rigorous when each layer in the structure consists of two dimensional (2D) or ine dimensional (1D) periodic photonic cyrstal or same periodicity in each layer.

REMS is theoretically based on the Rigorous Coupled Wave Analysis, Fourier Modal Method and Scattering Matrix method.

Model Repository

This simulation model is a publicly shared repository. Click here to access the repository or fork/clone the repository on the Kogence compute platform. Anybody can fork or clone this simulation model on Kogence platform by clicking on any of the file links below. One can then modify the project as needed. User can then decide to publish the modified project as a new public project or keep private to the user or to a chosen group of collaborators.

Model File Descriptions

REMS provides a Matlab syntax programming interface (API) to the main simulation engine. Users create and simulate new structures in REMS by providing three main types of files written in Matlab syntax: 1. User Defined Scripts ( the *.uds files), 2. Problem Definition File (the *.prb file) and 3. Command File (the *.cmd file). These files are "scripts", meaning they cannot contain or call user defined custom functions/subroutines. They can call other user defined custom scripts (i.e. other *.uds files). They can also call most common built-in Matlab subroutines/functions -- this becomes very handy as users can post process or plot their simulation results using common Matlab commands within REMS *.uds scripts as described below.

*.UDS Files

File extension uds stands for user defined script. As the name suggests, user can write as many scripts and use them as they wish by calling them from any of model files described in this section. Following three are the most common purpose for which UDS files are used in REMS.

Geometry UDS File

Geometry uds file is a script written in Matab syntax. We use this script to define and parameterize the geometry. Main objective of the script is to divide a given structure into several layers, discretize each layer into same number of in-plane grid points and then generate a 3D matrix of consituent material in each grid point and in each layer. One can achieve this goal by using any Matlab language construct one likes.

In the simulation model presented here, we defined 1D line gratings that have rectangular cross sections both in transverse and in-plane directions. Geometry is parameterized in terms of periodicity (PeriodX and PeriodY), duty cycles (FilFacX and FilFacY). number of layers, heights of each layer, materials of ridges and groves in each layer and material of substrate and the ambient medium.

These paramers are then assigned desired real values in the problem definition file (the *.PRB file) as described below.

Normal Vector Field UDS File

REMS provides several different types of core simulation engines. One of the simulation engines uses what is called a Normal Vector Field (NVF). NVF is a user defined arbitrary vector field defined in the plane of each layer of the structure such that it point normal to any interface between two materials. In the simulation model presented in this repository, we do not use the NVF simulation engine. So we actually do not need this *.uds file.

Post Processing UDS File

Users can write any Matlab syntax scripts in a *.uds file or processing and plotting the simulation results as they want. In this simulation repository we have used three files: StructureDefinition.uds, PlotAndSave.uds and PlotAndSaveLoop.uds. These post processing scripts are called in the *.cmd files as described below. PlotAndSave.uds is used to plot and save results when optical scattering from gratings of one fixed geometry is simulated. Once we have validated the results we scan/sweep various geometric parameters of the grating in order to find optimal design with respect to a given target (maximum absorption in this case). Such parameter sweep rules are defined in the *.cmd file and in this multiple simulation case we call a different PlotAndSaveLoop.uds file post processing and plotting main trends from the simulations of full ensemble of grating structures.

*.PRB File

The Problem Definition File (*.prb) is the main Matlab syntax script that defines one instance of the full problem that we want to simulate. It defines the geometry of the gratings, materials, wavelengths, angle of incidence, and polarization of the light. It also defines several other parameters that decides what is the desired outcome of the simulation -- do we want to simulate far field scattering or do we want to generate near-field optical electromagnetics field profiles etc. This files also defines which core simulation engine we want to use for simulation and defines some parameters that trades-off computation time versus the computation accuracy.

REMS Matlab API provides users an access to several variables, datastructures and functions that can be used to define the geometry and the problem that user wants to simulate. For example, REMS exposes a datastructure called REMS_PV (PV stands of Path Variables). Variables inside this datastructure REMS_PV.* help us define a project_name, path to InputFiles folder that contain the material properties of all the materials used in the geometry, path to GeometryDefDirectory that contains the user defined scripts (*.uds) to define the geometry that user wants to simulate among some other path variables as shown below.

% ---- File Names
%%
% project_name would be used as suffix for all files that are saved. For
% example data_project_name would have all the matlab workspace.
% struct_project_name would have the data that helps you visualize the
% structure or re-simulate the structure.

% ProjectDirectory is the full path to the root directory where you want
% project data files to be saved. Program would create Data, Figures,
% Log etc. directories under the ProjectDirectory location you specified.

% InputDirectory is the full path to the directory from where simulator
% would pick the required input files like refractive index data for
% various materials.

% GeometryDefDirectory is the full path to the directory from where
% simulator would pick the user defined geometry definition files.

% NVFDefDirectory is the full path to the directory from where
% simulator would pick the user defined geometry definition files.

% SimulatorRootPath is the full path to the directory inside which
% InternalFiles directory is installed.

REMS_PV.project_name='sim3_11';
REMS_PV.SimulatorRootPath='/opt/software/REMS/Version20/src';
REMS_PV.ProjectDirectory='/tmp/user1/job0';
REMS_PV.InputDirectory=[REMS_PV.ProjectDirectory '/InputFiles'];
REMS_PV.GeometryDefDirectory=[REMS_PV.ProjectDirectory '/GeometryDefinitions'];
REMS_PV.NVFDefDirectory=[REMS_PV.ProjectDirectory '/NVFDefinitions'];

In the next section of the prb file shown below we select the core simulation engine using another REMS_flags datastructure that the Matlab API of REMS exposes to the users. REMS provide several core simulation engines -- a heavily optimized version of the original RCWA algorithm published by Moharam and Gaylord, another engine based on the FMM formulation published by Lifeng Li, another engine that utilizes c2v mirror symmetries on top of the Lifeng Li's FMM formulation, an alternative engine that uses the Normal Vector Fields and lastly an engine optimized by Dr. Mukul Agrawal of Stanford University which is particularly well suited for simulating non-perioidic (aperiodic/random) structures.

Depending upon what type of scattering property user wants to simulate, REMS can further optimize the computation engine. Fastest computations happen when user is only interested in computing absolute power intensities of various scattered reflected and transmitted modes. We can also choose to calculate how much light is absorbed in each of the layers of the multi-layer structure. On the other hand, if user wants to calculate how much light is absorbed in each constituent materials then the compute engine has to first calculate full near-field optical electromagnetic field profiles and then perform numerical volume integration over each material. This selection consumes most computation time. REMS uses Newton-Cotes integration to perform these volume integration. API also exposes some flags by which user can select the order of Newton-Cotes integration.

% --- Various flags to choose the core simulation engine
%%

% REMS_flags.simulator_type:- 
% 1=MG RCWA 
% 2=LL RCWA without symmetry 
% 3=LL RCWA with symmetry 
% 4=NVF RCWA 
% 5=MA RCWA
REMS_flags.simulator_type=3;

% Only one of the following flags should be set to 1.
REMS_flags.only_abs_ref_trans_flag=1;
REMS_flags.layer_abs_identification_flag=0;
REMS_flags.material_region_identification_flag=0;
REMS_flags.field_profile_identification_flag=0;

% Following two parameters are needed only when
% material_region_identification_flag is switched ON.
% REMS_SP.type_of_integration= 1 , 2 and 3 are Newton-Cotes integration; 
% 4,5 and 6 are rectangular integration. 1, 4 keep all fields in memory; 
% 2, 5 keep fields within one layer in memory; 3, 6 keep fields at one Z 
% point in memory. 
% Its basically a memory versus computation time trade-off.
% REMS_SP.NCorder defines the order of Newton-Cotes integration. Orders 2 
% to 11 are programmed currently.

REMS_flags.type_of_integration=5;
REMS_flags.NCorder=5;

API also exposes some general purpose flags by which user can decide if data should be saved and/or overwritten. User can also choose to generate a 3D visulation of the structure that is being simulated.

% Multiple of the following flags can be set to 1.
REMS_flags.save_files_flag=1;
% If this flag is switched ON, it would save the final data file, the
% strcuture definitions, the logs and the figures (see the next two flags
% as well).
REMS_flags.overwrite_files_flag=1;
% If this flag is switched ON, existing files would be overwritten.
REMS_flags.logging_option_flag=1;
% If this flag is switched ON, then program would output logs on screen.
% Logs are also saved in the log directory if save_files_flag is switched
% ON.
REMS_flags.generate_new_structure_flag=1;
% If this flag is OFF, then GeometryGenerator program is skipped, instead a
% saved strcuture data file is uploaded. Also
% material_region_identification_flag and field_profile_identification_flag
% and grid parameters would be uploaded. You can still analyze the same
% structure for different wavelengths, angles, polarizations or different
% number of harmonics but with same old grid.
REMS_flags.structure_visualization_flag=0;
% If this flag is ON then several plots would be plotted for
% structure visualization and debugging purposes. 
%%

Next we define the geometry that we want to simulate. REMS Matlab API exposes a datastructure called structure parameters (REMS_SP). This datastructure is available for use in this prb script as well as inside the geometry definition uds script as discussed below. User can introduce new variables inside the REMS_SP.* datastructure inside the geometry definition uds script and then assign an actual value to those variables inside this prb file as shown below.

% --- Following paramters control the geometry you want to simulate
%%
% First five of the following Structure Parameters (SP) are required.
% Rest, user can define any variable inside the structure variable REMS_SP 
% that would help user define geometry in the User Defined Script (*.uds) 
% file that is placed inside the REMS_PV.GeometryDefDirectory directory.
% Only the variables defined inside REMS_SP and inside Grid Prameters
% structure variable REMS_GP would be available for use inside the User 
% Defined Scripts placed under the REMS_PV.GeometryDefDirectory directory.


REMS_SP.Structure_Type='rect_rect';
REMS_SP.PeriodX=500.0e-9;              
% absolute period along X in m
REMS_SP.PeriodY=500.0e-9;              
% absolute period along Y in m
REMS_SP.Cover_material={'air'};        
% materials here can not be numbers -- have to have string name with a file associated
REMS_SP.Substrate_material={'air'};     
% materials here can not be numbers -- have to have string name with a file associated
REMS_SP.groove_materials={'air'};      
% materials here can not be numbers -- have to have string name with a file associated
REMS_SP.ridge_materials={'ag_sopra'};        
% materials here can not be numbers -- have to have string name with a file associated
REMS_SP.DispX=[0];                     
% displacement of center of ridges from origin along X axis given as fraction of period along X and its absolute value can not be greater than 1.0
REMS_SP.DispY=[0];                     
% displacement of center of ridges from origin along Y axis given as fraction of period along Y and its absolute value can not be greater than 1.0
REMS_SP.FilFacX=[0.5];   
% length of ridges along X axis given as fraction of period along X and it can not be less than 0 or greater than 1.0
REMS_SP.FilFacY=[0.5];                 
% length of ridges along Y axis given as fraction of period along Y and it can not be less than 0 or greater than 1.0
% a homegeneous layer of groove_material can be generated by specifyng the FilFacX=FilFacY=0
% a homegeneous layer of ridge_material can be generated by specifyng the FilFacX=FilFacY=1
REMS_SP.depth=[50].*1.0e-9;            
% absolute depth of grating layer in m
%%

Next, we need to define some grid parameters (REMS_GP.*) and some harmonic parameters (REMS_HP.*). REMS computes the optical fields in reciprocal space. Therefore, REMS_HP.* parameters are the main knobs that trade-offs accuracy versus computation time. REMS_GP.* paramters help define the geometry on a grid that is then Fourier transformed to calculate Fourier components for the compute engine.

% --- Following paramters control the grid in real and Fourier spaces 
%%
% Z grid (i.e. Zp) is needed only if material_region_identification_flag or
% field_profile_identification_flag is switched ON.
% Start_Z and End_Z is is needed only if material_region_identification_flag or
% field_profile_identification_flag or layer_abs_identification_flag is switched ON.
% Variables defined under REMS_GP would also be available for use inside 
% the User Defined Scripts placed under the 
% REMS_PV.GeometryDefDirectory directory.

REMS_GP.Xp=1001;
REMS_GP.Yp=1001;
REMS_GP.Zp=1;
REMS_GP.Start_Z=25.0e-9;%0;
REMS_GP.End_Z=25.0e-9;%sum(REMS_SP.depth);
%%

% --- Following two paramters selects spatial frequencies
%%
% Following two parameters are needed to select which spatial frequencies
% are retained in the problem analysis. 
% For REMS_flags.simulator_type=5, REMS_HP.total_no_of_harmonics most 
% important frequencies are selected by Fourier analyzing the layer 
% REMS_HP.representative_layer_no.
% For all other simulator types, lowest REMS_HP.MX spatial frequencies are 
% selected along X and lowest REMS_HP.NY spatial frequencies are selected 
% along Y.

REMS_HP.MX=11;
REMS_HP.NY=11;
REMS_HP.total_no_of_harmonics=101;
REMS_HP.representative_layer_no=1;

Next we define the wavelengths, angle of incidence, and polarization of the light using the another external variables datastructure that API exposes (REMS_EV).

%%
% ---- Incident radiation
%%
% Polarization parameters

% psi is the polarization angle in degrees. psi=0 => TE polarization,
% psi=90 => TM polarization.Polarization angle is defined as the direction
% of incident E field with respect to plane of incidence (a uniquely
% defined entity even in case of non-symmetrical 2D diffraction).  For
% normal incidence, TE polarization is defined to mean incident E field
% oscillating along Y axis. 

% delta is the retardation phase angle from TE to TM in degrees. This
% parameter can be used to mix two polarizations (together with a non-TE
% and non-TM value of psi) differing in temporal phase of oscillations.


REMS_EV.psi	= 0;            
REMS_EV.delta	= 0; 		

% Wavelengths

% If wavelength becomes exactly equal to an integral multiple periodicity,
% end results may not be accurate. To avoid this, add a small irrational
% part on wavelength grid.

REMS_EV.start_lambda=(380.0+pi/100)*1.0e-9;
REMS_EV.end_lambda=(880.0+pi/100)*1.0e-9;
REMS_EV.no_of_lambda_points=51;      


% Polar angles

% Following parameters define the polar angles of incidence in cover medium
% in degree.

REMS_EV.start_angle_degree=0;
REMS_EV.end_angle_degree=0;
REMS_EV.no_of_angle_points=1;    


% Azimuthal angles

% Following parameters define the azimuth angles of incidence in cover
% medium in degree.

REMS_EV.start_azimuth_degree=0;
REMS_EV.end_azimuth_degree=0;
REMS_EV.no_of_azimuth_points=1; 

%%

*.CMD File

The Command File (*.cmd) is another Matlab syntax script that allows us to sweep variety of problem parameters and simulate multiple structures either simultaneously or one by one. Below is the simplest cmd file that is used to simulate only one structure.
% Following script when execute simulates the problem as defined by the user in the ProblemDefinitionFile (*.prb).
% This script helps in designing multitude of experiments and batch processing of such experiments.
% One can use entire Matlab scripting language in this script and make "for loops", "if-else" conditions etc.
% For example, one can change geometry variables, incidence conditions etc. and resimulate by repeating some of the commands below.


SimulatorInitializer;
% Needs to be called only once

GeometryProcessor;
% Needs to be called again whenever geometry variables (REMS_SP.* in *.prb file) or the grid variables (REMS_GP.* in *.prb file) are changed.

RefractiveIndexInterpolator;
% Needs to be called again whenever wavelength or materials in the geometry are changed.

IncidenceConditionProcessor;
% Needs to be called again whenever optical incidence conditions (REMS_EV.* in *.prb file) or the cover materials (REMS_SP.Cover_material in *.prb file) are changed.

StructureVisualization.uds;
% This is a UserDefinedScript (*.uds) and needs to be executed only when user wants.

Simulator;
% Needs to be called whenever a new simulation is required.

SimulatorFinisher;
% Needs to be called only once

PlotAndSave.uds; 
% This is a UserDefinedScript (*.uds) and needs to be executed only when user wants.

REMS is designed as a set of component modules. Command file calls various component modules of the REMS one by one. Main benefit of this approach is that when user wants to sweep certain problem parameter, user does not need to repeatedly call all the modules. This brings significant benefit in net computation time requires to solve a large design or optimization problem.

Following is an example of the cmd file that is used to loop over several problem parameters.

% Following script when execute simulates the problem as defined by the
% user in the ProblemDefinitionFile (*.prb). This script helps in designing
% multitude of experiments and batch processing of such experiments. One
% can use entire Matlab scripting language in this script and make "for
% loops", "if-else" conditions etc. For example, one can change geometry
% variables, incidence conditions etc. and resimulate by repeating some of
% the commands below.


SimulatorInitializer;
% Needs to be called only once

% thickness of first layer (closer to cover)
REMS_UDO.t1=10.0e-9:10.0e-9:100.0e-9;
% thickness of second layer (closer to substrate)
REMS_UDO.t2=10.0e-9:10.0e-9:160.0e-9;

for ii=1:1:length(REMS_UDO.t1)
    
    
    for jj=1:1:length(REMS_UDO.t2)
        
        REMS_SP.depth=[REMS_UDO.t1(ii),REMS_UDO.t2(jj)];
        
        GeometryProcessor;
        % Needs to be called again whenever geometry variables (REMS_SP.* in *.prb file) or the grid variables (REMS_GP.* in *.prb file) are changed.
        
        RefractiveIndexInterpolator;
        % Needs to be called again whenever wavelength or materials in the geometry are changed.
        
        IncidenceConditionProcessor;
        % Needs to be called again whenever optical incidence conditions (REMS_EV.* in *.prb file) or the cover materials (REMS_SP.Cover_material in *.prb file) are changed.
        
        %StructureVisualization.uds;
        % This is a UserDefinedScript (*.uds) and needs to be executed only when user wants.
        
        Simulator;
        % Needs to be called whenever a new simulation is required.
        
        
        
        cd(REMS_PV.InputDirectory);
        load('am15')
        cd(REMS_PV.ProjectDirectory);
        
        
        am15_wavelength=am15(:,1).*1.0e-9;                                                   % m
        am15_irradiance=am15(:,5).*1.0e6;                                                    % W/m^2/m; normalized to 1000 W/m^2
        input_irradiance=interp1(am15_wavelength,am15_irradiance,REMS_EV.wavelength,'cubic');            % W/m^2/m
        input_photon_count=abs(REMS_EV.wavelength(2)-REMS_EV.wavelength(1)).*input_irradiance.*REMS_EV.wavelength./REMS_PC.hbar./2./pi./REMS_PC.cc;            % number/s/m^2
        
        am15_weighted_reflected_current=sum(input_photon_count.*REMS_EV.reflectance_z_total.'./100)*REMS_PC.ee*1.0e3/1.0e4;          
        %mA/cm^2
        total_incident_current=sum(input_photon_count)*REMS_PC.ee*1.0e3/1.0e4; 
        %mA/cm^2
        REMS_UDO.am15_weighted_reflectivity(ii,jj)=am15_weighted_reflected_current*100/total_incident_current;      % in percentage
        
        
    end
    
    
end

SimulatorFinisher;
% Needs to be called only once

PlotAndSave.uds;
% This is a UserDefinedScript (*.uds) and needs to be executed only when user wants.

Results

Fig 2. Total optical power absorption (all modes), total transmitance and total reflectance is shown.
We used an 8 Intel cores with 32GB of RAM to simulate above described structure. Spectral distribution of absorbed power, reflected power and transmitted power is shown in the Fig. 2 on the right. Towards the shorter wavelengths, we see that the light not only scatters specularly they are some bragg scattering modes are also available in the far field scattering pattern. Figure 3 below shows that reflectance and transmittance scattered components. Left panel shows the components of reflection and right panel shows the components of transmission. In the left panel we see that the total reflected power is much higher than the specularly scattered power. As the wavelength of incident light increases the becomes smaller than the recreciprocal lattice vector . Therefore, even the first bragg scaattering mode is cut-off. Meaning, the component of the first scattered bragg mode is imaginary and therefore such a mode only exist in the near fields and not in far field patterns. Such modes also cannot carry any optical power. This is clear seen in the parameter known as haze. Haze is the ration of scattered power to the total power. We see that the reflectance haze in short wavelengths is very high and at 500nm wavelength it drops abruptly to zero. This is when the bragg scattering mode is cut off.

List of Other REMS Simulation Models on Kogence

Following are some more recently added public REMS models available for cloning/copying by anyone:


See the full list of public REMS models available for cloning/copying here: Category:REMS. Alternatively, you could also click on the category links at the bottom of this page to navigate to other simulation models in similar subject are or similar computational methodology.


References

Model New Results