In this example, we compute the reflection matrix of an open system with the input and output bases being Gaussian beams focused to different spatial locations. The magnitude of the diagonal elements of this matrix corresponds to what is measured in confocal microscopy. We first build the list of inputs B and list of outputs C, and then use mesti() to compute the reflection matrix.
In this example, we consider a test system of a dielectric cylindrical scatterer located at (x0, y0) in air.
clear
% System parameters
syst.length_unit = 'µm';
syst.wavelength = 1.0; % wavelength (µm)
syst.dx = syst.wavelength/20; % discretization grid size (µm)
nPML = 20; % number of PML pixels
W = 20; % width of simulation domain (including PML) (µm)
L = 10; % length of simulation domain (including PML) (µm)
r_0 = 0.5; % cylinder radius (µm)
n_bg = 1.0; % refractive index of the background
n_scat = 1.5; % refractive index of the cylinder
% Build the relative permittivity profile
nx = round(L/syst.dx);
ny = round(W/syst.dx);
x = (0.5:nx)*syst.dx;
y = (0.5:ny)*syst.dx;
x_0 = L/2; % location of the cylinder
y_0 = W/2; % location of the cylinder
[X, Y] = meshgrid(x, y);
epsilon = n_bg^2*ones(ny, nx);
epsilon((X-x_0).^2+(Y-y_0).^2 < r_0^2) = n_scat^2;
% Plot the relative permittivity profile
figure
imagesc(x, y, epsilon)
set(gca,'YDir','normal')
colorbar
colormap(flipud(gray))
title('$\varepsilon_{\rm r}(x,y)$','Interpreter','latex')
xlabel('x (µm)')
ylabel('y (µm)')
axis image
set(gca, 'FontSize', 18)
We consider inputs being Gaussian beams focused at (xf, yf). In this example, we fix the focal depth at xf = x0 (i.e., the depth of the scatterer), and scan the transverse coordinate yf of the focus.
Perfect Gaussian beams can be generated with the total-field/scattered-field (TF/SF) method. But since the cross section of the beam decays exponentially in y, we can generate Gaussian beams to a high accuracy simply by placing line sources at a cross section on the left, which is what we do here. We place the line sources at x = xsource, just in front of the PML.
To determine the required line sources, we (1) take the field profile of the desired incident Gaussian beam at the focal plane, Ein(xf, y) = E0exp(-(y - yf)2/w2), (2) project it onto the propagating channels (i.e., ignoring evanescent contributions) of free space, (3) back propagate it to the source plane to determine Ein(xsource, y) in the propagating-channel basis, and (4) determine the line source necessary to generate such Ein(xsource, y).
% Parameters of the input Gaussian beams
NA = 0.5; % numerical aperture
x_f = x_0; % location of the focus in x (fixed)
y_f_start = 0.3*W; % starting location of the focus in y
y_f_end = 0.7*W; % ending location of the focus in y
dy_f = syst.wavelength/(10*NA); % spacing of focal spots in y
% Parameters of the line sources
n_source = nPML + 1; % index of the source plane
x_source = x(n_source); % location of the source plane
w_0 = syst.wavelength/(pi*NA); % beam radius at x = x_f
y_f = y_f_start:dy_f:y_f_end; % list of focal positions in y
M_in = length(y_f); % number of inputs
% Step 1: Generate the list of E^in(x=x_f, y).
% Here, y.' is an ny-by-1 column vector, and y_f is a 1-by-M_in row vector.
% So, y.' - y_f is an ny-by-M_in matrix by implicit expansion.
% Then, E_f is an ny-by-M_in matrix whose m-th column is the cross section
% of the m-th Gaussian beam at x = x_f.
E_f = exp(-(y.' - y_f).^2/(w_0^2)); % size(E_f) = [ny, M_in]
% Get properties of propagating channels in the free space.
% We use PEC as the boundary condition for such channels since the default
% boundary condition in mesti() is PEC for TM waves, but the choice has
% little effect since E^in should be exponentially small at the boundary of
% the simulation domain.
channels = mesti_build_channels(ny,'TM','PEC',(2*pi/syst.wavelength)*syst.dx,n_bg^2);
% Transverse profiles of the propagating channels. Each column of u is
% one transverse profile. Different columns are orthonormal.
u = channels.fun_u(channels.kydx_prop); % size(u) = [ny, N_prop]
% Step 2: Project E^in(x_f, y) onto the propagating channels.
E_f_prop = u'*E_f; % size(E_f_prop) = [N_prop, M_in]
% Step 3: Back propagate from x = x_f to x = x_source.
% This step assumes a PEC boundary in y, so it is not exact with PML in y,
% but it is sufficiently accurate since E^in decays exponentially in y.
% Note we use implicit expansion here.
kx = reshape(channels.kxdx_prop/syst.dx, [], 1); % list of wave numbers
E_s_prop = exp(1i*kx*(x_source-x_f)).*E_f_prop; % size(E_s_prop) = [N_prop, M_in]
% Step 4: Determine the line sources.
% In a closed geometry with no PML in y, a line source of
% -2i*nu(a)*u(:,a) generates outgoing waves with transverse profile
% u(:,a). With PML in y, this is not strictly true but is sufficiently
% accurate since E^in(x=x_source,y) decays exponentially in y.
% Note we use implicit expansion here.
nu = reshape(channels.sqrt_nu_prop, [], 1).^2; % nu = sin(kx*dx)
B_L = u*(nu.*E_s_prop); % size(B_L) = [ny, M_in]
% We take the -2i prefactor out, to be multiplied at the end. The reason
% will be clear when we handle C below.
opts.prefactor = -2i;
% In mesti(), B_struct.pos = [m1, n1, h, w] specifies the position of a
% block source, where (m1, n1) is the index of the smaller-(y,x) corner,
% and (h, w) is the height and width of the block. Here, we put line
% sources (w=1) at n1 = n_source that spans the whole width of the
% simulation domain (m1=1, h=ny).
B_struct.pos = [1, n_source, ny, 1];
% B_struct.data specifies the source profiles inside such block, with
% B_struct.data(:, a) being the a-th source profile.
B_struct.data = B_L;
% We check that the input sources are sufficiently localized with little
% penetration into the PML; otherwise the Gaussian beams will not be
% accurately generated.
figure
imagesc([1, M_in], y, abs(B_L))
set(gca,'YDir','normal')
colorbar
colormap(flipud(gray))
title('|B_L|')
xlabel('Input index')
ylabel('y (µm)')
set(gca, 'FontSize', 18)
We consider output projections onto the same set of Gaussian beams focused at (xf, yf), with the projection done at the same plane as the source plane (x = xsource).
When the system has a closed boundary in y, as is the case in mesti2s(), the set of transverse modes form a complete and orthonormal basis, so it is clear what the output projection should be. But the Gaussian beams here are not orthogonal to each other, are not normalized, and do not form a complete basis. So, it is not obvious how our output projection should be defined.
What we do here is to convert everything onto the complete and orthonormal basis of transverse modes, and do the projection in such basis while accounting for the flux. Specifically, we (1) project the total field at the source plane, Etot(xsource, y) = Ein(xsource, y) + Esca(xsource, y), onto the propagating channels (i.e., ignoring evanescent contributions) of free space; the incident contribution will be subtracted later (2) back propagate such reflection to the focal plane at x = xf since the Esca(xsource, y) component supposedly comes from reflection, (3) take the previously computed Gaussian beams at the focal plane projected onto propagating channels of free space, and (4) take the inner product between the two while accounting for the longitudinal flux of the different propagating channels.
Above, the incident field Ein(x, y) was not subtracted. Contribution from the incident field will be subtracted using matrix D in the next step.
% We perform the output projection on the same plane as the line source
C_struct.pos = B_struct.pos;
% Step 1: Project E(x_source, y) onto the propagating channels.
% The projection will be C_L*E(:,n_source)
C_L = u'; % size(C_L) = [N_prop, ny]
% Step 2: Back propagate from x = x_source to x = x_f
C_L = exp(-1i*kx*(x_f-x_source)).*C_L; % size(C_L) = [N_prop, ny]
% Step 3: Project Gaussian beams at the focal plane onto the propagating
% channels. No need to repeat since this was already done earlier.
% E_f_prop = u'*E_f; % size(E_f_prop) = [N_prop, M_in]
% Step 4: Take the inner product between the two
% The longitudinal flux of a propagating channels is proportional to nu, so
% we weight the inner product with nu to account for flux dependence.
% Note we use implicit expansion here.
C_L = (E_f_prop') * (nu.*C_L); % size(C_L) = [M_in, ny]
% Normally, the next step would be
% C_struct.data = C_L.';
% However, we can see that C_L equals transpose(B_L)
fprintf('max(|C_L-transpose(B_L)|) = %g\n', max(abs(C_L - B_L.'), [], 'all'));max(|C_L-transpose(B_L)|) = 1.38778e-16
% That means we will have C = transpose(B). So, we can save some computing
% time and memory usage by specifying C = transpose(B).
% This is expected by reciprocity -- when the set of inputs equals the set
% of outputs, we typically have C = transpose(B) or its permutation.
clear C_struct
C = 'transpose(B)';The scattering matrix is given by S = C*inv(A)*B - D, with D = C*inv(A0)*B - S0 where A0 is a reference system for which its scattering matrix S0 is known. We consider A0 to be a homogeneous space with no scatterers, for which the reflection matrix S0 is zero.
syst.PML.npixels = nPML; % Put PML on all four sides
% For a homogeneous space, the length of the simulation domain doesn't
% matter, so we choose a minimal thickness of nx_temp = n_source + nPML
syst.epsilon = n_bg^2*ones(ny, n_source + nPML);
D = mesti(syst, B_struct, C, [], opts);System size: ny = 400, nx = 41; Ez polarization
UPML on -x +x -y +y sides; xBC = PEC; yBC = PEC
Building B,C... elapsed time: 0.001 secs
Building A ... elapsed time: 0.016 secs
< Method: APF using MUMPS with AMD ordering (symmetric K) >
Building K ... elapsed time: 0.005 secs
Analyzing ... elapsed time: 0.012 secs
Factorizing ... elapsed time: 0.035 secs
Total elapsed time: 0.080 secs
% Compute the reflection matrix
syst.epsilon = epsilon;
r = mesti(syst, B_struct, C, D, opts);System size: ny = 400, nx = 200; Ez polarization
UPML on -x +x -y +y sides; xBC = PEC; yBC = PEC
Building B,C... elapsed time: 0.000 secs
Building A ... elapsed time: 0.040 secs
< Method: APF using MUMPS with AMD ordering (symmetric K) >
Building K ... elapsed time: 0.017 secs
Analyzing ... elapsed time: 0.048 secs
Factorizing ... elapsed time: 0.242 secs
Total elapsed time: 0.374 secs
For most applications, it is not necessary to compute the full field profile, since most experiments measure properties in the far field. Here, we compute the full field profile for the purpose of visualizing the system as the incident Gaussian beams are scanned across y.
% Exclude the PML pixels from the returned field profiles
opts.exclude_PML_in_field_profiles = true;
x_interior = x((nPML+1):(nx-nPML));
y_interior = y((nPML+1):(ny-nPML));
field_profiles = mesti(syst, B_struct, [], [], opts);System size: ny = 400, nx = 200; Ez polarization
UPML on -x +x -y +y sides; xBC = PEC; yBC = PEC
Building B,C... elapsed time: 0.000 secs
Building A ... elapsed time: 0.046 secs
< Method: factorize_and_solve using MUMPS with AMD ordering >
Analyzing ... elapsed time: 0.086 secs
Factorizing ... elapsed time: 0.227 secs
Solving ... elapsed time: 0.341 secs
Total elapsed time: 0.732 secs
theta = linspace(0, 2*pi, 100);
circ_x = cos(theta);
circ_y = sin(theta);
cmap_bluered = colorcet_D09(); % use a blue-white-red colormap from ColorCET
% Loop through Gaussian beams focused at different locations
figure
for ii = 1:M_in
% Plot the total field profile; exclude PML
clf
ax1 = subplot(1,2,1);
imagesc(x_interior, y_interior, real(field_profiles(:, :, ii)));
set(gca,'YDir','normal')
caxis([-1,1]) % center the colorbar axis so zero = white
colormap(ax1, cmap_bluered);
hold on
plot(x_0+r_0*circ_x, y_0+r_0*circ_y, 'k-', 'linewidth', 0.8);
axis image off
% Plot the reflection matrix
ax2 = subplot(1,2,2);
imagesc(y_f, y_f, abs(r).^2)
set(gca,'YDir','normal')
colormap(ax2, flipud(gray))
xlabel('Input position')
ylabel('Output position')
hold on
% Show the current focal position
plot(y_f(ii)*[1,1], get(gca,'YLim'), 'b-')
axis image
xticks([])
yticks([])
set(gca, 'FontSize', 18)
drawnow
end