This page provides details on the numerical calculation carried out in the paper Interference of Conically Scattered Light in Surface Plasmon Resonance

## Overview

When a focussed beam is used to excite surface plasmon polaritions (SPPs) in a Kretchmann Attenuated Total Reflection (ATR) type setup, the reflected light in both the specular and conically scattered directions exhibits a unique one-sided spatial interference pattern that evolves with propagation into the far field.

This 2D spatial profile as a function of $$x$$ and $$z$$, $$E(x,z)$$ can be computed by taking the Fourier transform of the field on the interface, $$\tilde{E}(k_x,k_z)$$ multiplied by a free space transfer function $$\exp\left(\mathrm{i} k_z z\right)$$.


In practice the field in k-space can be represented by the Fresnel n-layer reflection multiplied by a (for example) Gaussian beam.

## Calculation

Numerically we present two ways (using GNU Octave) to compute this integral

• by using a numerical integration routine such as Gaussian quadrature, or
• with a fast Fourier transform (FFT)

The files may be downloaded here. This archive contains the following files

• fftwiggles.m - main routine for calculating the interference pattern
• nlayerfresnel.m - implements the n-layer Fresnel reflectivity for TM polarization
• showwidth.m - in case we oversample, this function rescales the output domain

To use these files, first set the variables as you please in fftwiggles.m. Then, call the function like so:

[x,E_spec,E_cone,E_gauss] = fftwiggles(500);

where ‘500’ in the above is an example propagation distance of 500 microns. You can plot the output with

plot(x,E_spec,x,E_cone,x,E_gauss)

Note the scaling which is applied in fftwiggles.m.

The example uses a three layer system, but the analysis works for an arbitrary number of layers. Just make sure that your variables d and epsilon contain the proper ordered list of permittivities.

If you wish to compute the integral using Gaussian quadrature, the following octave code will work.

out = zeros(1,N);
for i=1:N
f = @(kx) 1/sqrt(2*pi)*exp(1.0i.*kx.*x(i))\
.*exp(1.0i*sqrt(k0.*k0.*epsilon1-kx.*kx).*z)\
.*gausskx(kx).*nlayerfrensl(kx,epsilon,d);
out(i) = quadgk(f,k0*sqrt(epsilon1)*sin(theta-spread),k0*sqrt(epsilon1)*sin(theta+spread),\
'RelTol',reltol,'AbsTol',abstol, 'MaxIntervalCount',1e8,'Waypoints',kxi);
endfor
out = out.^2;