# Waveguide model for a subwavelength slit as a quarter-wave retarder # (Chapter 2 of my PhD dissertation, "Two-dimensional optics") Main file: simulation.py

simulation.py

import warnings
import numpy as N
import numpy.fft as FT
from scipy import integrate
from spp_generation import interface_calculation


class SlitSystem:
    """Simulation of a subwavelength slit in metal"""
    def __init__(self, slit_widths, wavelength, metal_epsilon, metal_thickness,
        incidence_index, slit_index, outcoupling_index, collection_index,
        numerical_aperture, angle_of_incidence,
        C1=1.0, C2=1.0, C3=1.0, caching=True):
        """
        Do the calculations for the subwavelength slit system.

        @slit_widths: An array of slit widths, in meters.
        @wavelength: Wavelength of the incident light, in meters.
        @metal_epsilon: Dielectric constant of the metal film.
        @metal_thickness: Thickness of the metal film, in meters.
        @incidence_index: Index of refraction of the medium from which the light
        is incident.
        @slit_index: Index of refraction of the medium inside the slit.
        @outcoupling_index: Index of refraction of the medium on the other side
        of the slit.
        @collection_index: Index of refraction of the medium in which the
        collection lens is situated.
        @numerical_aperture: numerical aperture of the collection lens.
        @angle_of_incidence: Angle under which the plane wave is incident on the
        slit.
        @C1, @C2, @C3: fitting parameter.
        @caching: whether to try to read the waveguide calculation from a file,
        and write it out if it is not there.
        """
        # Shorthand for parameters
        b = slit_widths
        wl = wavelength
        eps = metal_epsilon
        z = metal_thickness
        n1 = incidence_index
        n2 = slit_index
        n3 = outcoupling_index
        n4 = collection_index
        NA = numerical_aperture
        theta = angle_of_incidence

        k0 = 2 * N.pi / wl

        # P. Lalanne, J. P. Hugonin, J. C. Rodier. (2006). "Approximate model
        # for surface-plasmon generation at slit apertures." JOSA A 23 (7),
        # p. 1608.
        #
        # Model is valid for b <= wl, approximately.
        if N.any(b > wl):
            warnings.warn('Model is not valid for slit widths larger than the '
                + 'wavelength.')

        # eps_Ti = -5.649 + 25.787j  # Johnson & Christy
        _, _, alpha12, beta12 = interface_calculation(b, wl, eps, n1, n2, theta)
        _, _, alpha23, _ = interface_calculation(b, wl, eps, n3, n2, theta)

        # Snyder & Love, pp. 240-244.

        # Propagation constants of fundamental slit waveguide mode
        if caching:
            cache_filename = 'waveguide_eff_indices_{0}_{1}.txt'.format(int(wl * 1e9), b.size)
            try:
                beta_TE, beta_TM = N.loadtxt(cache_filename, unpack=True).view(complex)
            except IOError:
                from waveguide import effective_mode_indices
                beta_TE, beta_TM = effective_mode_indices(b, eps, n2, wl,
                    filename=cache_filename)
        else:
            from waveguide import effective_mode_indices
            beta_TE, beta_TM = effective_mode_indices(b, eps, n2, wl, filename=None)

        # Calculate collection efficiency of the optics.
        # First we calculate the mode profile at the exit of the slit, in a
        # one-dimensional half-space. (The modes are symmetric.)
        num = 2048
        num_padded = 32768
        lim = b.max() * 4
        space = N.linspace(0, lim, num=num).reshape(num, 1)

        # Extend the space for each slit width
        x = space.repeat(b.size, 1)

        # Mask for space inside the slit
        slit = x < b / 2

        # Electric field at end of slit
        root_TE = N.sqrt(k0 ** 2 * n2 ** 2 - beta_TE ** 2)
        root_TM = N.sqrt(k0 ** 2 * n2 ** 2 - beta_TM ** 2)
        root2_TE = N.sqrt(beta_TE ** 2 - k0 ** 2 * eps)
        root2_TM = N.sqrt(beta_TM ** 2 - k0 ** 2 * eps)

        Ey_TE_metal = N.exp(-(N.abs(x) - b / 2) * root2_TE)
        Ey_TE_metal[slit] = 0.0
        Ey_TE_slit = N.cos(x * root_TE) / N.cos(b * root_TE / 2)
        Ey_TE_slit[~slit] = 0.0
        Ey_TE = (Ey_TE_metal + Ey_TE_slit) * N.exp(1j * beta_TE * z)

        Ex_TM_metal = (n2 ** 2 / eps) * N.exp(-(N.abs(x) - b / 2) * root2_TM)
        Ex_TM_metal[slit] = 0.0
        Ex_TM_slit = N.cos(x * root_TM) / N.cos(b * root_TM / 2)
        Ex_TM_slit[~slit] = 0.0
        Ex_TM = (Ex_TM_metal + Ex_TM_slit) * N.exp(1j * beta_TM * z)

        self.field_Ey_TE = N.r_[Ey_TE[::-1], Ey_TE]
        self.field_Ex_TM = N.r_[Ex_TM[::-1], Ex_TM]
        self.field_x = N.r_[-x[::-1], x]

        # Fourier transform of electric field
        self.Fy_TE = FT.fft(Ey_TE, n=num_padded, axis=0)[:num] / num
        self.Fx_TM = FT.fft(Ex_TM, n=num_padded, axis=0)[:num] / num

        # Convert from spatial frequencies to angles
        spatial_frequencies = FT.fftfreq(num_padded, d=lim / num)[:num]
        max_spatial_frequency = spatial_frequencies[spatial_frequencies < (n3 / wl)].argmax() + 1
        # spatial_frequencies is in cycles / unit, so no extra factor of 2pi
        angles = N.append(N.arcsin(wl * spatial_frequencies[:max_spatial_frequency - 1] / n3), [N.pi / 2])
        self.spatial_frequencies = spatial_frequencies

        self.Fy_TE2 = N.abs(self.Fy_TE[:max_spatial_frequency, :]) ** 2
        self.Fx_TM2 = N.abs(self.Fx_TM[:max_spatial_frequency, :]) ** 2
        self.angles = angles

        # Calculate Fresnel reflection coefficient at n3-n4 interface
        if n3 == n4:
            self.R_TE, self.R_TM = N.zeros_like(angles), N.zeros_like(angles)
        else:
            critical_angle = N.arcsin(n4 / n3)
            mask = angles < critical_angle
            # For numerical computation considerations, we set R=1 above the
            # critical angle; this
            term1, term2 = N.ones_like(angles), N.zeros_like(angles)
            term1[mask] = n3 * N.cos(angles[mask])
            term2[mask] = n4 * N.sqrt(1 - ((n3 / n4) * N.sin(angles[mask])) ** 2)
            self.R_TE = ((term1 - term2) / (term1 + term2)) ** 2

            term1, term2 = N.ones_like(angles), N.zeros_like(angles)
            term1[mask] = n3 * N.sqrt(1 - ((n3 / n4) * N.sin(angles[mask])) ** 2)
            term2[mask] = n4 * N.cos(angles[mask])
            self.R_TM = ((term1 - term2) / (term1 + term2)) ** 2

        # This is what the objective sees
        self.profile_TE = (1 - self.R_TE)[:, N.newaxis] * self.Fy_TE2
        self.profile_TM = (1 - self.R_TM)[:, N.newaxis] * self.Fx_TM2

        # Integrate the intensity picked up by the objective normalized by the
        # total intensity
        max_collection_angle = N.arcsin(NA / n3)
        NA_cutoff = angles < max_collection_angle
        self.collection_eff_TE = \
            integrate.trapz(self.profile_TE[NA_cutoff], angles[NA_cutoff], axis=0) \
            / integrate.trapz(self.Fy_TE2, angles, axis=0)
        self.collection_eff_TM = \
            integrate.trapz(self.profile_TM[NA_cutoff], angles[NA_cutoff], axis=0) \
            / integrate.trapz(self.Fx_TM2, angles, axis=0)

        # Fresnel reflection and transmission coefficients to and from the slit
        # mode
        self.r21_TM = (k0 * n1 - beta_TM) / (k0 * n1 + beta_TM)
        t12_TM = (k0 * n1 / beta_TM) * (1 - self.r21_TM)
        self.r23_TM = (k0 * n3 - beta_TM) / (k0 * n3 + beta_TM)
        t23_TM = (beta_TM / (k0 * n3)) * (1 + self.r23_TM)
        r21_TE = (beta_TE - k0 * n1) / (beta_TE + k0 * n1)
        t12_TE = 1 - r21_TE
        r23_TE = (beta_TE - k0 * n3) / (beta_TE + k0 * n3)
        t23_TE = 1 + r23_TE

        # Fabry-Perot transmission coefficients
        self.t123_TM = t12_TM * t23_TM * N.exp(1j * beta_TM * z) \
            / (1 - self.r23_TM * self.r21_TM * N.exp(2j * beta_TM * z))
        self.t123_TE = t12_TE * t23_TE * N.exp(1j * beta_TE * z) \
            / (1 - r23_TE * r21_TE * N.exp(2j * beta_TE * z))

        # Surface plasmon generation coefficients
        self.s12 = beta12 + (t12_TM * self.r23_TM * alpha12
            * N.exp(2j * beta_TM * z)
            / (1 - self.r23_TM * self.r21_TM * N.exp(2j * beta_TM * z)))
        self.s23 = (t12_TM * alpha23 * N.exp(1j * beta_TM * z)
            / (1 - self.r23_TM * self.r21_TM * N.exp(2j * beta_TM * z)))

        # Total transmission (taking into account loss to surface plasmons on
        # both sides)
        self.T_TM = C1 * self.collection_eff_TM * (
            (n3 / n1) * N.abs(self.t123_TM) ** 2
            - 2 * N.abs(self.s12) ** 2
            - 2 * C2 * N.abs(self.s23) ** 2)
        self.T_TE = self.collection_eff_TE * (n3 / n1) ** 2 * N.abs(self.t123_TE) ** 2

        # Phase lag between TM and TE modes
        self.dphase = (N.angle(self.t123_TM) - N.angle(self.t123_TE)) \
            % (2 * N.pi)

        # Output results as object members
        self.beta_TE, self.beta_TM = beta_TE, beta_TM

if __name__ == '__main__':
    b = N.arange(10, 830, 2) * 1e-9
    b2 = N.arange(10, 830, 2) * 1e-9
    calc = SlitSystem(b,
        wavelength=830e-9,
        metal_epsilon=-29.321 + 2.051j,
        metal_thickness=200e-9,
        incidence_index=1.00027487,
        slit_index=1.00027487,
        outcoupling_index=1.5102,
        collection_index=1.00027487,
        numerical_aperture=0.40,
        angle_of_incidence=0.0)
    calc2 = SlitSystem(b2,
        wavelength=830e-9,
        metal_epsilon=-29.321 + 2.051j,
        metal_thickness=200e-9,
        incidence_index=1.5102,
        slit_index=1.00027487,
        outcoupling_index=1.00027487,
        collection_index=1.00027487,
        numerical_aperture=0.40,
        angle_of_incidence=0.0)

    from matplotlib import pyplot as P
    fig = P.figure()
    ax = fig.gca()
    ax.plot(b * 1e9, calc.collection_eff_TM, label='forward TM')
    ax.plot(b * 1e9, calc.collection_eff_TE, label='forward TE')
    ax.plot(b2 * 1e9, calc2.collection_eff_TM, label='reverse TM')
    ax.plot(b2 * 1e9, calc2.collection_eff_TE, label='reverse TE')
    ax.set_xlabel('Slit width (nm)')
    ax.set_ylabel('Fraction of total radiated energy collected by objective')
    ax.legend(loc='best')
    ax.set_xlim(0, 500)
    P.show()

spp_generation.py

"""Plasmon generation efficiency from fundamental slit mode
P. Lalanne, J. P. Hugonin, J. C. Rodier. (2006). "Approximate model for
surface-plasmon generation at slit apertures." JOSA A 23 (7), p. 1608."""

import numpy as N
from scipy import special


def I0_integrand(u, w):
    return N.sinc(N.outer(u, w)) ** 2


def I1_integrand(u, w, gamma):
    # The full integrand (weighted by sqrt(1 - u^2)) is singular at -1 and +1.
    # Adaptive QUADPACK quadrature and AdaptSim (Matlab's quad() function)
    # both fail on it for |x| > 1, presumably because there are poles close to
    # the real axis.
    # Gaussian quadrature succeeds. A fixed order of 200 seems to produce the
    # most accurate results balanced by computing time. However, the weighted
    # Gaussian quadrature from Lalanne's code works much faster.
    gamma_u = N.sqrt(N.asarray(1 - u ** 2, dtype=complex))
    wu = N.outer(u, w)
    return N.sinc(wu) * N.exp(-1j * N.pi * wu) / (gamma_u + gamma)[:, N.newaxis]


def lalanne_integral(func, order=20, *args):
    """Gaussian quadrature of a function, adapted from Lalanne's code.

    Computes integral of func(x, *args) / sqrt(1 - x^2) from -inf to +inf.
    Does this by splitting it up into two parts, [0, 1) and (1, inf].
    """
    points, weights = special.orthogonal.p_roots(order)
    # Points and weights for the first integral
    x1, w1 = (N. pi / 4) * (points + 1), (N.pi / 4) * weights
    # Points and weights for the second integral
    x2, w2 = (points + 1) / 2, weights / 2

    # for [0, 1), u = cos(x1)
    gauss1 = func(N.cos(x1), *args) + func(-N.cos(x1), *args)
    # for (1, inf], uu = sqrt(1 - uuu^2), uuu = 1 / x2 - 1
    uuu = 1 / x2 - 1
    uu = N.sqrt(uuu ** 2 + 1)
    gauss2 = (func(uu, *args) + func(-uu, *args)) \
        * (((uuu + 1) ** 2) / uu)[:, N.newaxis]

    return (w1[:, N.newaxis] * gauss1).sum(axis=0) \
        - 1j * (w2[:, N.newaxis] * gauss2).sum(axis=0)


def calc_I0(w_norm):
    """Calculate the I0 integral for normalized slit width w'."""
    return lalanne_integral(I0_integrand, 20, w_norm)


def calc_I1(w_norm, epsilon, n):
    """Calculate the I1 integral for normalized slit width w'.

    @w_norm: normalized slit width
    @epsilon: dielectric constant of metal
    @n: index of material on outside of slit"""
    gamma_SP = -N.sqrt(n ** 2 / (n ** 2 + epsilon))
    return lalanne_integral(I1_integrand, 120, w_norm, gamma_SP)


def interface_calculation(w, wl, epsilon, n1, n2, theta):
    """Surface-plasmon generation dynamics at a slit aperture.

    Returns r0, t0, alpha, beta
    @r0: Reflection coefficient from inside the slit to outside
    @t0: Transmission coefficient from outside the slit to inside
    @alpha: SP generation coefficient for illumination from the fundamental mode
      inside
    @beta: SP generation coefficient for illumation from a plane wave outside

    @w: slit width
    @wl: wavelength
    @n1: index of medium outside the slit
    @n2: index of medium inside the slit
    @theta: angle of incidence of plane wave"""
    w_norm = w * n1 / wl
    temp = (n1 / n2) * w_norm * calc_I0(w_norm)
    r0 = (temp - 1) / (temp + 1)

    alpha = (-1j * (1 - r0) * calc_I1(w_norm, epsilon, n1)
        * N.sqrt((w_norm * n1 ** 2 / (n2 * N.pi))
        * (N.sqrt(N.abs(epsilon)) / (-epsilon - n1 ** 2))))

    t0 = ((1 - r0) * N.sinc(w_norm * N.sin(theta))
        ) * N.sqrt(n1 / (n2 * N.cos(theta)))

    beta = -alpha * t0 / (1 - r0)

    return r0, t0, alpha, beta


def test_integrals():
    """Test our numerical integration by comparing it to Table 1 in the paper,
    which tabulates certain values"""
    import time

    bn = N.array([0.1, 0.3, 0.5, 0.7])
    epsilon = -26.27 + 1.85j

    I0_time = -time.clock()
    I0 = calc_I0(bn)
    I0_time += time.clock()
    I1_time = -time.clock()
    I1 = calc_I1(bn, epsilon, 1.0)
    I1_time += time.clock()

    print 'I0 took {0:.4} s'.format(I0_time)
    print 'I1 took {0:.4} s'.format(I1_time)

    I0_ref = N.array([3.09 - 4.09j, 2.72 - 1.68j, 2.13 - .63j, 1.54 - .18j])
    I1_ref = N.array([.53 - 2.93j,  1.75 - 1.8j,  1.79 - .4j,  1.01 + .35j])
    dI0 = N.abs(I0_ref - I0) / N.abs(I0_ref)
    dI1 = N.abs(I1_ref - I1) / N.abs(I1_ref)

    print 'Table 1. I1 and I0 for Gold at 800 nm (epsilon=-26.27+1.85i)'
    print " w' {0[0]:10.2} {0[1]:10.2} {0[2]:10.2} {0[3]:10.2}".format(bn)
    print ' I0 {0[0]:10.2f} {0[1]:10.2f} {0[2]:10.2f} {0[3]:10.2f}'.format(I0)
    print 'dI0 {0[0]:10.0%} {0[1]:10.0%} {0[2]:10.0%} {0[3]:10.0%}'.format(dI0)
    print ' I1 {0[0]:10.2f} {0[1]:10.2f} {0[2]:10.2f} {0[3]:10.2f}'.format(I1)
    print 'dI1 {0[0]:10.0%} {0[1]:10.0%} {0[2]:10.0%} {0[3]:10.0%}'.format(dI1)

if __name__ == '__main__':
    test_integrals()

waveguide.py

import numpy as N
from scipy.optimize import fsolve

# Snyder & Love, pp. 240-244


# Trick to get scipy to solve complex equations
# http://mail.scipy.org/pipermail/scipy-user/2006-November/009943.html
def pack_into_real(cy):
    '''Pack a complex array into a real array.'''
    # Too clever by half, but it works.
    return N.column_stack([cy.real, cy.imag]).ravel()


def unpack_into_complex(y):
    '''Unpack a complex array from a real array.'''
    return y[::2] + 1j * y[1::2]


# Eigenvalue equations
def ev_TE(beta_packed, epsilon, n2, b, wl):
    beta2 = unpack_into_complex(beta_packed) ** 2
    k02 = (2 * N.pi / wl) ** 2
    W1 = N.sqrt(beta2 - k02 * epsilon)
    U1 = N.sqrt(k02 * n2 ** 2 - beta2)
    return pack_into_real(W1 / U1 - N.tan(0.5 * b * U1))


def ev_TM(beta_packed, epsilon, n2, b, wl):
    beta2 = unpack_into_complex(beta_packed) ** 2
    k02 = (2 * N.pi / wl) ** 2
    n22 = n2 ** 2
    W1 = N.sqrt(beta2 - k02 * epsilon)
    U1 = N.sqrt(k02 * n22 - beta2)
    return pack_into_real(n22 * W1 / (epsilon * U1) - N.tan(0.5 * b * U1))


# Initial guesses for solver
def TE_guess(b, wl):
    """Initial TE guess"""
    ny_TE_guess = N.array((0.439 * wl) / b, dtype=complex)
    beta_TE_guess = 2 * N.pi * N.sqrt(1 - ny_TE_guess ** 2) / wl
    return beta_TE_guess


def TM_guess(b):
    """Initial TM guess"""
    ny_TM_guess = 8e6 + 1e5j * N.ones_like(b)
    return ny_TM_guess


def _solve_beta(b, epsilon, n2, wl, beta_TE_guess, beta_TM_guess, verbose=False):
    """Solve the eigenvalue equations numerically"""

    # Output arrays
    beta_TE = N.zeros(b.shape, dtype=complex)
    beta_TM = N.zeros_like(beta_TE)

    for ix, width in enumerate(b):
        if verbose:
            print ix,
            if ix % 10 == 9:
                print ' '

        # first TE
        beta_packed, _, ier, mesg = fsolve(ev_TE,
            pack_into_real(beta_TE_guess[ix]),
            args=(epsilon, n2, width, wl),
            full_output=True)
        if ier != 1:
            print "For width={0} nm, no TE solution found: {1}".format(width * 1e9, mesg)
            beta_TE[ix] = N.nan
        else:
            beta_TE[ix] = unpack_into_complex(beta_packed)[0]

        # then TM
        beta_packed, _, ier, mesg = fsolve(ev_TM,
            pack_into_real(beta_TM_guess[ix]),
            args=(epsilon, n2, width, wl),
            full_output=True)
        if ier != 1:
            print "For width={0} nm, no TM solution found: {1}".format(width * 1e9, mesg)
            beta_TM[ix] = N.nan
        else:
            beta_TM[ix] = unpack_into_complex(beta_packed)[0]

    return beta_TE, beta_TM


def effective_mode_indices(b, epsilon, n2, wl, filename='waveguide_eff_indices.txt', verbose=True):
    beta_TE_guess = TE_guess(b, wl)
    beta_TM_guess = TM_guess(b)

    beta_TE, beta_TM = _solve_beta(b, epsilon, n2, wl, beta_TE_guess, beta_TM_guess, verbose)

    if verbose:
        print ' '
    if filename is not None:
        N.savetxt(filename, N.column_stack([
            beta_TE.view(float).reshape(-1, 2),
            beta_TM.view(float).reshape(-1, 2)
        ]))

    return beta_TE, beta_TM

if __name__ == '__main__':
    import matplotlib.pyplot as P

    b = N.arange(10, 1000, 2) * 1e-9
    wl = 830e-9
    epsilon = -29.321 + 2.051j
    n2 = 1.0
    teguess = TE_guess(b, wl)
    tmguess = TM_guess(b)
    te, tm = _solve_beta(b, epsilon, n2, wl, teguess, tmguess)

    fig = P.figure()
    ax = fig.gca()
    ax.plot(b * 1e9, te.real, 'b-')
    ax.plot(b * 1e9, teguess.real, 'b--')
    ax.plot(b * 1e9, te.imag, 'r-')
    ax.plot(b * 1e9, teguess.imag, 'r--')
    ax.plot(b * 1e9, tm.real, 'k-')
    ax.plot(b * 1e9, tmguess.real, 'k--')
    ax.plot(b * 1e9, tm.imag, 'g-')
    ax.plot(b * 1e9, tmguess.imag, 'g--')
    ax.set_ylim(-1e7, 10e7)
    P.show()