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

- 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;
```

`fftwiggles.tar.gz`

- the complete archive`fftwiggles.m`

- main routine`nlayerfresnel.m`

- n-layer Fresnel reflectivity for TM polarization`showwidth.m`

- rescales the output