COMSOL 6.4 - Simulation of Metal–Air Surface Plasmon Polariton Propagation and Dispersion

Simulation of Metal–Air Surface Plasmon Polariton Propagation and Dispersion

Electromagnetic waves that are confined to propagate along a surface, such as surface plasmon polaritons (SPPs), are of great research interest due to their potential applications in nanoscale manipulation of light. This model will discuss how to set up a simulation to visualize the propagation as well as the frequency versus propagation constant dispersion relationship of SPPs.

Propagating electromagnetic waves can exist in a variety of forms — such as plane waves, spherical waves, Gaussian beams, and so on. There are also propagating electromagnetic waves that are confined in space, such as waveguide modes propagating in metallic or dielectric waveguides.

In addition, there is another special type of electromagnetic wave that is confined to a planar surface. This type of wave propagates tangentially to the surface with exponentially decaying fields in the perpendicular direction. Its wavelength is often smaller than the free space wavelength at the same frequency. Thus, this type of wave provides promises for nanoscale control and manipulation of photons, which is desirable in many applications, ranging from optical communication and information processing to solar energy harvesting and digital displays.

The existence of this wave type was found at metal–dielectric interfaces and is known as surface plasmon polariton. The term plasmon refers to the collective oscillation of charges in a metal, whereas the term polariton generally refers to a quasi-particle that is the result of the coupling between a photon with a polar excitation in a material. The SPPs are the result of coupling of surface plasmons and light.

Model Definition

This model demonstrates SPPs in the simplest system that supports it, namely the bulk metal–dielectric interface. Imagine that the metal–dielectric interface is represented by the xz plane at y = 0. The dielectric region is above the plane, y > 0, and the metal region is below the plane, y < 0. The surface wave propagates in the x direction.

Figure 1: A metal–dielectric interface at y = 0. In this model, the SPP propagates in the x direction and exponentially decays in the y direction.

Here, it is of interest to study a TM mode surface wave. The expressions for the electric and magnetic fields in the dielectric and the metal are then given by

(1)

(2)

(3)

(4)

where the d and m superscripts denote quantities in dielectric and metallic domains, respectively. β is the complex SPP propagation constant. Both γd and γm are positive real numbers that describe the decay of the field away from the metal–dielectric interface.

The tangential components of the electric and magnetic fields, and the normal component of the electric displacement field, are continuous across the metal–dielectric boundary. Therefore,

(5)

(6)

(7)

Since there is no external charge and the permittivity is constant inside the respective metallic and dielectric domains, ∇⋅E = 0 must hold inside the two regions. Combining this condition with Equation 2 and Equation 4, gives the following two equations

(8)

(9)

Subsequently, if those two equations are combined, you get the simple relation

(10)

From this relation, it is clear why SPPs only exist between a dielectric, having a positive real part of the permittivity, and a metal, having a negative real part of the permittivity. To have the field decay in the y direction, both γd and γm must be positive, which means that the permittivities εd and εm must have opposite signs.

To derive the expression for the propagation constant β, the Helmholtz equation is used

(11)

where k0 = ω/c is the free space wave number.Plugging Equation 2 and Equation 4 into the Helmholtz equation leads to

(12)

(13)

Finally, combining Equation 10, Equation 12, and Equation 13, the following expression for the SPP propagation constant is derived

(14)

The real part of β is related to the SPP wavelength by λSPP = 2π/Re(β), while the imaginary part describes the SPP propagation loss. Generally, the permittivities εd and εm are frequency dependent, so β is also frequency-dependent. The relationship between β and the frequency is often what needs to be known to characterize the SPPs in a system.

Remember that the above discussion is purely based on the assumption that the SPP is a TM wave. For the possibility of the TE wave, one can simply follow the previous derivation steps and show that all the field amplitudes must be zero. This means that SPPs only exist as a TM wave.

This model simulates SPPs propagating at the interface between silver (metal) and air (dielectric), using material properties for silver from the built-in Optical Material Library. Numeric ports are used on the left and right boundaries of the model. The left port, with excitation turned on, will launch the SPP, while the right port, with excitation turned off, will absorb the SPP without reflection. To get the mode field for both ports, two Boundary Mode Analysis study steps are added. Finally, a Frequency Domain study step is used when solving for the domain field.

Results and Discussion

Figure 2 shows the wave as it propagates in the x direction. The decay in the y direction is faster on the metal side due to the strong absorption.

Figure 2: The y-component of the electric field for a wavelength of 370 nm (3.35 eV). The wave gets weaker, due to absorption, as it propagates from left to right. The arrows indicate the electric field strength and direction.

Figure 3 confirms that there is strong absorption for short wavelengths, whereas for longer wavelengths the transmission increases.

Figure 3: A plot of the transmittance and absorptance versus wavelength.

The mode field plot in Figure 4 clearly shows that the SPP decays much faster in the silver domain than in the air domain.

Figure 4: The norm of the electric mode field for Port 1.

To capture the quantitative relationship between the frequency (photon energy) and β, a plot of freq*h_const on the y-axis versus the propagation constant ewfd.beta_1 on the x-axis (shown by the circle markers) is shown in Figure 5.

When studying SPPs, it is customary to define a figure of merit, often referred to as the Q factor, as the ratio of the real and the imaginary part of β. When β has a smaller imaginary part (equivalently, a larger Q factor), the SPP can propagate a long distance relative to its wavelength before it is attenuated. A larger Q factor is usually desirable for practical applications. The Q factor can conveniently be plotted using a color expression added to the dispersion curve. Here, a brighter color represents higher Q factors and a darker color represents lower Q factors. In addition, a dashed line representing f = c/λ, often called the light line, is added. The light line is the dispersion relation of free space photons. Lastly, the analytical expression from Equation 13 is plotted as a solid curve. The simulated dispersion and the analytical expression show good agreement.

The dispersion plot in Figure 5 is very representative for the SPP dispersion in noble metals. This plot is useful for gaining insight of the characteristics of SPPs. More importantly, it shows that the dispersion curve of SPPs always lies on the right side of the light line. The implication of this is that the SPP wavelength is always smaller than the wavelength of free space waves. This is why SPPs can be used as a way to compress the wavelength of light to achieve a higher field concentration.

Furthermore, the mismatch between the free-space wave number and the SPP propagation constant means that it is not possible to excite SPPs simply by shining light onto the metal surface. Some external mechanism to achieve phase matching is required. The excitation of SPPs is often done using total internal reflection from a prism or using diffraction from a grating. The goal of with these techniques is to prepare the electromagnetic field, so that its propagation constant matches that of the SPP at the same frequency.

Figure 5: A plot of the simulated dispersion curve of the SPP at the interface between silver and air. The simulated result (circles) is consistent with the analytical calculation (solid line). The free space light dispersion, or light line, is represented by a dashed line. The color represents the Q factor of SPP.

For more SPP simulation examples, see the blog post in Ref. 1.

Reference

1. X. Chen, “Modeling Surface Plasmon Polaritons in COMSOL®,” COMSOL Blog, 12 Oct. 2022; www.comsol.com/blogs/modeling-surface-plasmon-polaritons-in-comsol/.

Wave_Optics_Module/Waveguides/metal_air_surface_plasmon_polariton

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
L	2[um]	2E-6 m	Simulation domain length
H	3[um]	3E-6 m	Simulation domain height
lda0	400[nm]	4E-7 m	Wavelength
f0	c_const/lda0	7.4948E14 1/s	Frequency

The constant c_const above represents the speed of light.

Geometry 1

The geometry is very simple, consisting of only one rectangle that is split into two halves.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type L.

In the text field, type H.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (nm)
Layer 1	H/2

Click

Materials

Pick silver from the Optical material library.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > > > .

Click the button in the window toolbar.

Materials

Air

In the window, under right-click and choose .

In the window for , type Air in the text field.

Select Domain 2 only.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Refractive index, real part	n_iso ; nii = n_iso, nij = 0	1	1	Refractive index

In the toolbar, click

to close the window.

Electromagnetic Waves, Frequency Domain (ewfd)

Port 1

In the toolbar, click

and choose .

Select Boundaries 1 and 3 only.

In the window for , locate the section.

From the list, choose .

Port 2

In the toolbar, click

and choose .

Select Boundaries 6 and 7 only.

In the window for , locate the section.

From the list, choose .

Study 1

Step 1: Boundary Mode Analysis

In the window, under click .

In the window for , locate the section.

In the text field, type f0.

Select the checkbox. In the associated text field, type 5.

Step 3: Boundary Mode Analysis 1

Right-click > and choose .

Right-click and choose .

In the window for , locate the section.

In the text field, type 2.

Step 3: Frequency Domain

In the window, click .

In the window for , locate the section.

In the text field, type f0.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
lda0 (Wavelength)		nm

Click

In the dialog, type 340 in the text field.

In the text field, type 10.

In the text field, type 600.

Click .

Definitions

Analytic 1 (an1)

In the toolbar, click

In the window for , type epsilonr in the text field.

Locate the section. In the text field, type (mat1.rfi.nr(x)-1i*mat1.rfi.ni(x))^2.

Click to expand the section. Select the checkbox.

This analytical function calculates the relative permittivity of silver.

Mesh 1

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

Size 2

Right-click > and choose .

Drag and drop below .

In the window for , locate the section.

From the list, choose .

Select Boundary 4 only.

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type (lda0/real(sqrt(epsilonr(lda0/1[m])/(1+epsilonr(lda0/1[m])))))/12.

The maximum element size should be much smaller than SPP wavelength.

Click

The mesh at the metal-air interface is finer than the one built using Physics-controlled mesh.

Study 1

In the toolbar, click

Results

Electric Field (ewfd)

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Surface 1

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type ewfd.Ey.

Locate the section. From the list, choose .

From the list, choose .

Arrow Surface 1

In the window, right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose > > > .

Locate the section. Find the subsection. In the text field, type 50.

Find the subsection. In the text field, type 50.

Locate the section. From the list, choose .

In the toolbar, click

Select the checkbox. In the associated text field, type 1e-3.

Transmittance and Absorptance (ewfd)

Before creating the frequency versus propagation constant dispersion plot, it is good to inspect some of the default plots that are automatically generated when performing a simulation with Numeric ports.

The reflectance is small and not of interest here. So, curves related to the reflectance will be removed from the plot.

In the window, under click .

In the window for , type Transmittance and Absorptance (ewfd) in the text field.

Locate the section. In the text field, type Transmittance and absorptance (1).

Locate the section. From the list, choose .

Global 1

In the window, expand the node, then click .

In the window for , locate the section.

Ctrl-click to select table rows 1 and 3.

Click

In the toolbar, click

This figure shows that the absorption is large for short wavelengths and falls off to low values for longer wavelengths.

Electric Mode Field, Port 1 (ewfd)

In the window, under click .

This figure shows the electric mode field norm for . It is clear that the wave decays rapidly in the metal and is extending more into the dielectric domain.

SPP Dispersion

In the toolbar, click

In the window for , type SPP Dispersion in the text field.

Locate the section. From the list, choose .

Click to expand the section. From the list, choose .

In the text area, type SPP dispersion of bulk silver.

Locate the section.

Select the checkbox. In the associated text field, type Propagation constant (1/m).

Locate the section. From the list, choose .

Global 1

Right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Expression	Unit	Description
freq*h_const	eV	Energy

Locate the section. From the list, choose .

From the list, choose .

In the text field, type ewfd.beta_1.

Click to expand the section. Find the subsection. From the list, choose .

From the list, choose .

Find the subsection. From the list, choose .

Click to expand the section. From the list, choose .

In the table, enter the following settings:


Legends
Simulated SPP dispersion

Global 2

Right-click and choose .

In the window for , locate the section.

In the text field, type ewfd.k0.

Locate the section. Find the subsection. From the list, choose .

From the list, choose .

Find the subsection. From the list, choose .

Locate the section. In the table, enter the following settings:


Legends
Light line

Color Expression 1

In the window, right-click and choose .

In the window for , locate the section.

In the text field, type -real(ewfd.beta_1)/imag(ewfd.beta_1).

Locate the section. From the list, choose .

Click to expand the section. Select the checkbox.

In the text field, type 0.

In the text field, type 800.

Global 3

Right-click and choose .

In the window for , locate the section.

In the text field, type ewfd.omega/c_const*sqrt(epsilonr(lda0)/(epsilonr(lda0)+1)).

Locate the section. Find the subsection. From the list, choose .

From the list, choose .

Find the subsection. From the list, choose .

Locate the section. In the table, enter the following settings:


Legends
Analytic SPP dispersion

Color Expression 1

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type real((ewfd.omega/c_const)*sqrt(epsilonr(lda0)/(epsilonr(lda0)+1)))/imag(-(ewfd.omega/c_const)*sqrt(epsilonr(lda0)/(epsilonr(lda0)+1))).

Locate the section. Clear the checkbox.

In the toolbar, click