Plasmonic resonator (DGTD)

This example uses the discontinuous Galerkin time-domain (DGTD) method to model and analyze a plasmonic structure which exhibits a Fano-like resonance. It demonstrates how to compute the extinction spectrum, how to identify the individual resonance frequencies and quality factors and how to obtain the fields of the individual modes.

download example

Overview

Understand the simulation workflow and key results

Plasmonic resonators confine and strongly enhance electromagnetic fields at the nanoscale. In addition, by designing structures which exhibit both radiative (bright) and non-radiative (dark) modes, it is possible to engineer Fano-like resonances which are characterized by relatively sharp and potentially asymmetric spectral features.

One example of such a resonator is the Dolmen structure shown below. It consists of three bars made of gold. The material parameters and the exact dimensions are taken from Ref. [1] and are shown below. In a simplified picture, we can consider the top bar as a dipole which couples well to a plane wave polarized along the x-axis (bright mode). The pair of vertical bars at the bottom can be viewed as an electrical quadrupole which poorly couples to a plane wave (dark mode). The coupling of the bright mode to the dark mode and the resulting interference causes the characteristic spectral feature.

Run and results

Instructions for running the model and discussion of key results

Step 1: Plane-wave excitation and cross section calculation

The script enables a plane wave source with propagation direction along the z-axis and with polarization \( \vec{E} _0=(2,1,0) \) which corresponds to a polarization angle of approximately 26.5 degrees. This angle is chosen to break any obvious symmetry and to allow coupling to all modes. The plane wave is injected using a rectangular total-field/scattered-field (TF/SF) region. A frequency domain monitor ‘flux_monitor’ is used to record the flux through the surface of the TF/SF box.

After the simulation is finished, the 'sigma_front' flux result corresponds to the scattering cross-section and the 'sigma_back' to the absorption cross-section. The sum of scattering and absorption cross section is the extinction cross section. The script will automatically plot all three cross sections as shown below

Step 2: Extraction of resonance frequencies/Q-factors and mode profiles

The script re-configures the simulation to exploit the symmetry of the structure by introducing a mirror plane, normal to the x-axis, through the center of the simulation domain. Instead of running a single simulation of the entire device, the script will run two separate simulations of half the structure. This allows to separate and more easily identify the individual modes. By setting the boundary conditions at the mirror plane to either “PEC” (first run) or “PMC” (second run), we obtain only the anti-symmetric or only the symmetric modes, respectively.

Since not all modes can be excited via plane waves, the script disables the plane wave source and instead uses a global excitation with a randomized but divergence-free spatial profile. This feature can be manually enabled by checking the “Use global excitation” group in the “Advanced“-tab of the “DGTD solver” object.

The “point_monitor” is used to record the electric field at a point close to the bottom of the vertical bar. This time trace can then be analyzed using the “findresonances” script command. As a result, we find three modes with the following parameters.

To find the corresponding mode profiles, the script uses time-data for the entire TF/SF region as recorded by the "volume_monitor". From the call of "findresonances", we already know the number of modes (N=3) and the respective parameters \( \omega \), \( \alpha_k \) and \( \phi_k \) for each mode. With this, we can then describe the time-dependent fields \( \vec E(\vec x,t) \) as

$$ \vec E(\vec x,t) \approx \sum_{k=1}^N \vec E_k(\vec x)cos(\omega_kt + \phi_k)e^{- \alpha_kt}  $$

This equation can then be solved for the individual modes \( \vec E_k(\vec x) \) in a least-squares sense by using singular value decomposition (SVD). The script automatically performs this analysis step and generates the dataset 'plasmonicModes' containing the three modes.

Below, we show the \( \vec E_z \) component in a plane normal to the z-axis and intersecting at z=15nm (i.e. 5nm above the structure). The sign of the \( \vec E_z \) component in a plane closely above (or below) the structure closely resembles the charge distribution on the surface of a planar particle (as indicated by the "+" and "-" symbols).

Important model settings

Description of important objects and settings used in this model

Materials

For this example, gold is modelled using a Drude/Plasma model. The parameters are taken from Ref. [1].

Global excitationWhen using the “use global excitation” option to excite the plasmonic modes in Step 2, it is important to specify the frequency/wavelength range of interest, which is 0.55μm - 2μm in this example.

Boundary conditions

For this example, basic absorbing boundary conditions (ABCs) are used instead of Perfectly Matched Layers (PMLs). While ABCs typically absorb less radiation than PMLs, they offer easier handling (no free parameters) and require less computation time. However, one needs to ensure that the outer surface of the simulation domain is sufficiently far away from the scatterer. In practice, a distance of at least  \(\frac{\lambda_{max}}{2} \) is often recommended, where \( \lambda_{max} \) is the largest wavelength of interest. For this example, a radius of 0.5μm for the spherical simulation domain gives good results.

Point monitor

The point monitor data is used to extract the resonance frequencies and quality factors of the modes. One therefor has to pick a location that has good overlap with all modes. In this example we use the point at \( \vec x=(30nm,-95nm,10nm) \) which is at the bottom of the right vertical bar. Other points such as between the vertical and the horizontal bar work equally well but they must be sufficiently close to the structure.

Updating the model with your parameters

Instructions for updating the model based on your device parameters

The general procedure for extracting modes is directly applicable to any type of plasmonic scatterer. However, depending on the specific structure, it will be necessary to adjust several settings:

• The separation of modes by symmetry used in Step 2 is only possible if the underlying structure has the corresponding mirror symmetry. For a structure that does not have obvious mirror planes, one must run a single simulation of the full structure.
• The location of the point monitor must be chosen such that it has good overlap with all modes. If the modes are not known in advance, it makes sense to add several point monitors, all relatively close to the structure.
• The frequency/wavelength range used for the excitation (of both the plane wave and the global excitation) must be chosen to include the resonances of the structure. If there is no prior knowledge about the spectral location of the resonances, it makes sense to first run an excitation with a very broad-band pulse.

Taking the model further

Information and tips for users that want to further customize the model

Resonator on a substrate

This example assumes that the resonator is embedded in a homogeneous background material. A number of modifications are required if the structure is on a substrate or membrane. The global excitation and mode extraction techniques can be applied without any modification, but it is recommended to use a rectangular simulation domain instead of the sphere. In addition, the TF/SF injection of a plane wave currently only supports homogeneous background materials, which makes step 1 difficult to replicate.

Array of resonators

This example uses an individual structure, but for experiments it is often easier to study large, periodic arrays. For periodic arrays, one can model only a single unit cell and then use periodic boundary conditions on the sides of the unit cell. Instead of calculating cross sections, one then studies the reflection and transmission by monitoring the flux through planes above and below the structure.

Convergence

The default settings for this example were chosen to give reasonably accurate results while keeping the simulation time short. The relative deviations from the resonance frequencies reported in Ref. [1] are below 1% while the relative deviations for the Q-factors are below 10%. To improve the accuracy further, one can simply increase the polynomial order in the DGTD solver object from 3 to 4. This roughly halves the maximum deviations for both the frequencies and the Q-factors.

Additional resources

Additional documentation, examples and training material

Related publications

1. R. Faggiani et al., “Modal analysis of the ultrafast dynamics of optical nanoresonators”, ACS Photonics 2017, 4, 4, 897-904

Related Ansys Innovation Courses