I need help on how to plot Surface Plasmon Polaritons (SPPs) for Multilayer Systems

I have derived the surface plasmon polariton (SPP) relation for the system under consideration. The system consists of three layers: semi-infinite layers of InSb as the upper layer and SiC (lower layer), separated by a vacuum gap of 10 nm. I have written a code to calculate the parallel wavevector (G) and plot it against frequency (ω). However, the results seem incorrect, as I only obtain numerical values for some frequency points. Could you help me identify and correct the error? Here is my code:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

clc

% Physical parameters

hbar = 1.054571596e-34;

h = 2 * pi * hbar;

global c0;

c0 = 2.99792458e+8;

kb = 1.3806503e-23;

u0 = 4 * pi * 10^-7;

ep0 = 8.854 * 10^-12;

qe = 1.602176462e-19;

T_h0 = [350];

T_l = [300];

T_SiC = T_h0;

Tl_SiC = T_l;

T_InSb = T_h0;

Tl_InSb = T_l;

d = 10e-9;

N = 200;

% Initialization of different parameters

w1 = 0.01e14;

w2 = 5e14;

w = linspace(w1, w2, N);

lambda = (2 * pi * c0) ./ w;

w_ev = w / 1.5193E15;

ep1 = 1+0*w; %gap

ep2 = epsilon_SiC1w_T(w, Tl_SiC);

ep3 = epsilon_InSbw_T(w, T_InSb);

% Initialize arrays to store solutions

G_solutions = zeros(N, 1);

abs_G_solutions = zeros(N, 1);

for i = 1:N

syms G

k0 = w(i) / c0; % vacuum

k1=k0; %gap

k2 = sqrt(G^2 - ep2(i) * k0.^2); % SiC material

k3 = sqrt(G^2 - ep3(i) * k0.^2); % InSb material

% Dispersion equation for three-layer system: InSb/gap/SiC

equation = exp(-2 * k1 * d) - (((k2 * ep1 + k1 * ep2(i)) .* (k3 * ep1 + k1 * ep3(i))) ./ ...

((k2 * ep1 - k1 * ep2(i)) .* (k3 * ep1 - k1 * ep3(i)))) == 0;

% Solve the equation

solutions = solve(equation, G);

numerical_solutions = vpa(solutions, 10);

% Check if numerical_solutions is empty

if ~isempty(numerical_solutions)

G_solutions(i) = double(numerical_solutions(1)); % Assuming first solution is desired

abs_G_solutions(i) = abs(G_solutions(i));

else

G_solutions(i) = NaN; % Assign NaN if no solution found

abs_G_solutions(i) = NaN;

end

end

% Display the numerical solutions

disp('Numerical solutions for G at different frequencies:');

disp(G_solutions);

disp('Absolute values of G at different frequencies:');

disp(abs_G_solutions);

% Plot the scatter plot

figure(1)

scatter(w, abs_G_solutions, 'filled')

xlabel('Frequency (w)')

ylabel('Absolute value of G (|G|)')

title('Scatter plot of |G| vs Frequency (w)')

grid on

figure(2)

set(gcf, 'Units', 'pixels', 'Position', [100, 100, 500, 300]);

scatter(w,abs_G_solutions,'blue','filled','LineWidth',3);

xlabel('Energy (eV)','FontSize',14);

ylabel('Absolute value of G (|G|)')

set(gca,'linewidth',3);

set(gca,'FontSize',14);

set(0, 'DefaultAxesFontName', 'Times New Roman');

Answer 1

[Ronit]

Hello Ambali,

To address the issues, you're encountering with plotting Surface Plasmon Polaritons (SPPs) for a multilayer system consisting of ‘InSb’ and ‘SiC’ layers separated by a vacuum gap, it seems that your MATLAB code has several key areas for improvement.

% Updated syms declaration to specify G as positive for physical solutions
syms G positive

% Updated vpasolve call with a specified range for G
solutions = vpasolve(equation, G, [0, inf]); % Limiting G to positive values

% Store the first solution if it exists
if ~isempty(solutions)
    G_solutions(i) = double(solutions(1));
    abs_G_solutions(i) = abs(G_solutions(i));
else
    G_solutions(i) = NaN; % Assign NaN if no solution is found
    abs_G_solutions(i) = NaN;
end

Explanation for the above changes:

• Positive ‘G’ Declaration: By declaring ‘G’ as positive in the symbolic variable declaration, the solver is directed to look for positive solutions, aligning with the physical expectation that the wavevector component should be real and positive for propagating modes.
• Specifying Range in ‘vpasolve’: The updated ‘vpasolve’ call includes a range [0, inf] for ‘G’. This helps the solver by providing a physical range of interest, potentially improving the accuracy and efficiency of the solution process.
• Handling No Solution Cases: The conditional check for empty solutions ‘(if ~isempty(solutions))’ ensures that the code can gracefully handle scenarios where no solution is found by the solver, assigning ‘NaN’ to the ‘G_solutions’ and ‘abs_G_solutions’ arrays for those frequencies. This prevents the code from crashing or throwing errors when solutions are not found.

Hope this helps!