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


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)\).

\[ \newcommand{\intinfty}{\int_{-\infty}^{\infty}} \newcommand{\me}{{\mathrm{e}}} \newcommand{\mi}{{\mathrm{i}}} \newcommand{\md}{\,\mathrm{d}} \begin{align} E(x,z) &= \intinfty \tilde{E}(k_x)\, \me^{\mi k_z z} \me^{\mi k_x x} \md k_x \end{align} \]

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


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

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

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


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))\
    out(i) = quadgk(f,k0*sqrt(epsilon1)*sin(theta-spread),k0*sqrt(epsilon1)*sin(theta+spread),\
    'RelTol',reltol,'AbsTol',abstol, 'MaxIntervalCount',1e8,'Waypoints',kxi);
out = out.^2;