Reproducing Water Defect Mode Plot in MATLAB (Only Getting 0.1748 Transmittance Instead of Sharp Peak of 1)

Hi everyone,

I'm trying to reproduce a plot of water defect modes in a one-dimensional photonic biosensor using MATLAB, but I'm not getting the sharp transmission peak that I expected based on a research paper. Instead, I'm seeing a peak of only around 0.1748, whereas the research paper's plot shows a sharp peak with a transmittance close to 1.

Background: The design I'm working on is for a photonic biosensor intended to detect several types of waterborne bacteria. The defect mode should show a sharp resonance peak when pure water is used as the defect layer in the photonic crystal (PhC) structure. I'm using the transfer matrix method to model the transmittance.

However, I seem to be missing something or making an error in my code that's preventing me from getting the correct results. I’ve tried adjusting parameters like the defect layer thickness and the design wavelength, but it still doesn’t match the expected result. I’m new to this, so any help would be appreciated!

clc;
clear;
close all;


n_Si = 3.465;    % Refractive index of Si
n_SiO2 = 1.442;  % Refractive index of SiO2
n_water = 1.333; % Refractive index of water (defect layer)

% Layer thicknesses based on the quarter-wave condition
lambda_design = 1700; % Design wavelength in nm
d_Si = lambda_design / (4 * n_Si);    % Thickness of Si layer
d_SiO2 = lambda_design / (4 * n_SiO2);% Thickness of SiO2 layer
d_defect = d_Si + d_SiO2; % Defect layer thickness (1D)

% Wavelength range for simulation (1100 nm to 2750 nm)
lambda_range = linspace(1100, 2750, 1000); % Wavelength range in nm

% Number of periods
N = 10; 

% Calculate transmittance for each wavelength
transmittance = arrayfun(@(lambda) transfer_matrix_method(lambda, N, n_Si, n_SiO2, n_water, d_Si, d_SiO2, d_defect), lambda_range);

% Plot the transmittance spectrum (Figure 2(a))
figure;
plot(lambda_range, transmittance, 'LineWidth', 2);
xlabel('Wavelength (nm)');
ylabel('Transmittance');
title('Transmittance Spectrum for N = 10, d_d = 1D with Pure Water Defect Layer');
grid on;

% Zoomed-in view around the defect mode (Figure 2(b))
figure;
plot(lambda_range, transmittance, 'LineWidth', 2);
xlim([1800 1900]); % Narrow range to zoom into the defect mode
xlabel('Wavelength (nm)');
ylabel('Transmittance');
title('Zoomed-in View of Water Defect Mode');
grid on;


function T = transfer_matrix_method(lambda, N, n_Si, n_SiO2, n_defect, d_Si, d_SiO2, d_defect)
    n_air = 1; % Refractive index of air
    n_substrate = 1.517; % Refractive index of substrate (glass)

    % Initial transmission matrix
    M_total = eye(2);
    
    for i = 1:N/2
        % Si layer matrix
        delta_Si = 2 * pi * n_Si * d_Si / lambda;
        M_Si = [cos(delta_Si), 1i * sin(delta_Si) / n_Si; 1i * n_Si * sin(delta_Si), cos(delta_Si)];
        
        % SiO2 layer matrix
        delta_SiO2 = 2 * pi * n_SiO2 * d_SiO2 / lambda;
        M_SiO2 = [cos(delta_SiO2), 1i * sin(delta_SiO2) / n_SiO2; 1i * n_SiO2 * sin(delta_SiO2), cos(delta_SiO2)];
        
        % Multiply matrices for Si/SiO2 bilayers
        M_total = M_Si * M_SiO2 * M_total;
    end
    
    % Defect layer matrix
    delta_defect = 2 * pi * n_defect * d_defect / lambda;
    M_defect = [cos(delta_defect), 1i * sin(delta_defect) / n_defect; 1i * n_defect * sin(delta_defect), cos(delta_defect)];
    
    % Complete structure (AB)^N/2 D (AB)^N/2
    M_total = M_total * M_defect * M_total;
    
   
    A = M_total(1,1);
    B = M_total(1,2);
    C = M_total(2,1);
    D = M_total(2,2);
   
    T = (4 * n_air * n_substrate) / abs(A * n_substrate + B + C * n_air * n_substrate + D * n_air)^2;
end

This is the original plot

This is what I am getting

Answer 1

[Alan Stevens]

The size of the peak is very sensitive to the increments in your lambda range:

n_Si = 3.465;    % Refractive index of Si
n_SiO2 = 1.442;  % Refractive index of SiO2
n_water = 1.333; % Refractive index of water (defect layer)

% Layer thicknesses based on the quarter-wave condition
lambda_design = 1700; % Design wavelength in nm
d_Si = lambda_design / (4 * n_Si);    % Thickness of Si layer
d_SiO2 = lambda_design / (4 * n_SiO2);% Thickness of SiO2 layer
d_defect = d_Si + d_SiO2; % Defect layer thickness (1D)

% Wavelength range for simulation (1100 nm to 2750 nm)
lambda_range = linspace(1100, 2750, 6603); % Wavelength range in nm %%%%%%%%%%%%%%%%%%%%%%
% Number of periods
N = 10; 

% Calculate transmittance for each wavelength
transmittance = arrayfun(@(lambda) transfer_matrix_method(lambda, N, n_Si, n_SiO2, n_water, d_Si, d_SiO2, d_defect), lambda_range);

% Plot the transmittance spectrum (Figure 2(a))
figure;
plot(lambda_range, transmittance, 'LineWidth', 2);
xlabel('Wavelength (nm)');
ylabel('Transmittance');
title('Transmittance Spectrum for N = 10, d_d = 1D with Pure Water Defect Layer');
grid on;

% Zoomed-in view around the defect mode (Figure 2(b))
% figure;
% plot(lambda_range, transmittance, 'LineWidth', 2);
% xlim([1800 1900]); % Narrow range to zoom into the defect mode
% xlabel('Wavelength (nm)');
% ylabel('Transmittance');
% title('Zoomed-in View of Water Defect Mode');
% grid on;
% 

function T = transfer_matrix_method(lambda, N, n_Si, n_SiO2, n_defect, d_Si, d_SiO2, d_defect)
    n_air = 1; % Refractive index of air
    n_substrate = 1.517; % Refractive index of substrate (glass)

    % Initial transmission matrix
    M_total = eye(2);
    
    for i = 1:N/2
        % Si layer matrix
        delta_Si = 2 * pi * n_Si * d_Si / lambda;
        M_Si = [cos(delta_Si), 1i * sin(delta_Si) / n_Si; 1i * n_Si * sin(delta_Si), cos(delta_Si)];
        
        % SiO2 layer matrix
        delta_SiO2 = 2 * pi * n_SiO2 * d_SiO2 / lambda;
        M_SiO2 = [cos(delta_SiO2), 1i * sin(delta_SiO2) / n_SiO2; 1i * n_SiO2 * sin(delta_SiO2), cos(delta_SiO2)];
        
        % Multiply matrices for Si/SiO2 bilayers
        M_total = M_Si * M_SiO2 * M_total;
    end
    
    % Defect layer matrix
    delta_defect = 2 * pi * n_defect * d_defect / lambda;
    M_defect = [cos(delta_defect), 1i * sin(delta_defect) / n_defect; 1i * n_defect * sin(delta_defect), cos(delta_defect)];
    
    % Complete structure (AB)^N/2 D (AB)^N/2
    M_total = M_total * M_defect * M_total;
    
   
    A = M_total(1,1);
    B = M_total(1,2);
    C = M_total(2,1);
    D = M_total(2,2);
   
    T = (4 * n_air * n_substrate) / abs(A * n_substrate + B + C * n_air * n_substrate + D * n_air)^2;
end