Plasmonic cavity resonator consisting of a gold nanorod array#

Plasmonic structures have been under intense research interest in the past few decades, primarily because of their unique ability to confine light to subwavelength dimensions, significantly enhancing electromagnetic fields at their surfaces. This phenomenon has opened up new avenues in various scientific fields, including photovoltaics, sensing, and nanoscale imaging, by enabling the manipulation of light at the nanometer scale.

In this notebook, we demonstrate the modeling of a hexagonal lattice of gold nanorods, which supports localized surface plasmon resonances. When excited by a plane wave, these resonances manifest as distinct dips in the reflectance spectrum, indicative of the efficient transfer of incident power into the plasmon modes. Through the examination of field profiles at these resonance points, it becomes possible to discern the standing waves correlating to various harmonics. The design is adopted from David P. Lyvers, Jeong-Mi Moon, Alexander V. Kildishev, Vladimir M. Shalaev, and Alexander Wei, ACS Nano 2008, 2, 12, 2569–2576 DOI: 10.1021/nn8006477.

Schematic of the nanorod array

This notebook mainly aims to demonstrate how to model a hexagonal lattice within the FDTD simulation domain, which is rectangular. In addition, the plane wave source in Tidy3D by default has an incident power of 1 W. In this example, we also demonstrate how to perform field normalization to obtain a meaningful relative field enhancement factor. For more discussion on source normalization, refer to the FAQ entry.

Periodic or quasi-periodic structures are common in photonic devices. In Tidy3D’s example library, we have demonstrated a gradient metasurface reflector, a metalens at the visible frequency, and a graphene metamaterial absorber. In addition, we also investigated a frequency selected surface for the microwave frequency range.

[1]:

import numpy as np
import matplotlib.pyplot as plt

import tidy3d as td
import tidy3d.web as web

Simulation Setup#

Define the simulation wavelength (frequency) range.

[2]:

lda0 = 1.15  # central wavelength
freq0 = td.C_0 / lda0  # central frequency
ldas = np.linspace(0.4, 1.9, 201)  # wavelength range
freqs = td.C_0 / ldas  # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs))  # width of the source frequency range

Define the geometric parameters for the nanorod array.

[3]:

h = 0.649  # height of the nanorods
a = 0.098  # spacing between the nanorods
d = 0.0636  # diameter of the nanorods
t_base = 0.05  # thickness of the base plate

Define the materials used in the simulation. Tidy3D’s material library provides several options for gold. In this particular case, we will use the data from the widely used paper by Johnson and Christy in 1972.

The Al\(_2\)O\(_3\) is defined by a constant refractive index material in the wavelength range of interest.

[4]:

# define gold material from the material library
gold = td.material_library["Au"]["JohnsonChristy1972"]

# define Al2O3 material as a constant refractive index
n_eff = 1.61  # refractive index of Al2O3
al2o3 = td.Medium(permittivity=n_eff**2)

Next we define the structures, which include the dielectric filling, two nanorods, and a base plate.

The nanorods are arranged in a hexagonal lattice. In FDTD, the simulation domain needs to be rectangular. For a hexagonal lattice, we can identify a rectangular unit cell as indicated by the blue dashed box below. For a plane wave excitation, we can apply periodic boundary conditions for normal incidence or Bloch boundary conditions for oblique incidence.

In addition, we notice that the unit cell structure has two reflection symmetries. Therefore, in the case of normal incidence of a plane wave linearly polarized either in the horizontal or vertical directions, we can utilize the symmetry to further reduce the simulation domain size by a factor of 4. For more details on symmetry, please refer to the tutorial here. Next we are going to take this approach and only model the unit cell.

Schematic of the unit cells in a haxagonal lattice

[5]:

# define the Al2O3 filling structure
dielectric = td.Structure(
    geometry=td.Box(center=(0, 0, h / 2), size=(td.inf, td.inf, h)), medium=al2o3
)

# define the first nanorod structure
rod_1 = td.Structure(
    geometry=td.Cylinder(center=(0, 0, h / 2), radius=d / 2, length=h), medium=gold
)

# define the second nanorod structure
rod_2 = td.Structure(
    geometry=td.Cylinder(center=(a / 2, np.sqrt(3) * a / 2, h / 2), radius=d / 2, length=h),
    medium=gold,
)

# define the third nanorod structure
rod_3 = td.Structure(
    geometry=td.Cylinder(center=(a / 2, -np.sqrt(3) * a / 2, h / 2), radius=d / 2, length=h),
    medium=gold,
)

# define the forth nanorod structure
rod_4 = td.Structure(
    geometry=td.Cylinder(center=(-a / 2, np.sqrt(3) * a / 2, h / 2), radius=d / 2, length=h),
    medium=gold,
)

# define the fifth nanorod structure
rod_5 = td.Structure(
    geometry=td.Cylinder(center=(-a / 2, -np.sqrt(3) * a / 2, h / 2), radius=d / 2, length=h),
    medium=gold,
)

# define the base plate structure
base = td.Structure(
    geometry=td.Box(center=(0, 0, -t_base / 2), size=(td.inf, td.inf, t_base)), medium=gold
)

Next we define the source and monitors. A PlaneWave polarized in the \(y\) direction is injected from the top of the nanorods. By default the plane wave is polarized in the \(x\) direction. To make it \(y\)-polarized, we set pol_angle to \(\pi\)/2.

To measure reflection, we define a FluxMonitor above the PlaneWave. In addition, we define a FieldMonitor to help visualize the field profiles at the resonances. In principle, one needs to run the simulation first without the FieldMonitor to determine the resonance wavelengths (frequencies) and then rerun it with the FieldMonitor set to record the resonance frequencies. For presentation purposes, the resonance wavelengths (frequencies) are pre-determined here.

[6]:

# define a plane wave incident on the nanorod array from the top
plane_wave = td.PlaneWave(
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    size=(td.inf, td.inf, 0),
    center=(0, 0, h + 0.4 * lda0),
    direction="-",
    pol_angle=np.pi / 2,
)

# define a flux monitor above the plane wave to measure reflection
flux_monitor = td.FluxMonitor(
    center=[0, 0, h + 0.5 * lda0], size=[td.inf, td.inf, 0], freqs=freqs, name="R"
)

# pre-determined resonance wavelengths
res_ldas = np.array([1.6, 1.25, 0.93, 0.765, 0.675])
# corresponding resonance frequencies
res_freqs = td.C_0 / res_ldas

# define a field monitor at the resonance frequencies to visualize the excited modes
monitor_field = td.FieldMonitor(
    center=[a / 4, a * np.sqrt(3) / 4, h / 2],
    size=[0, a * np.sqrt(3) / 2, 2 * h],
    freqs=res_freqs,
    name="field",
)

With everything in place, we are ready to define the Simulation. Since the simulation domain is very small, we can afford to use a very fine grid without increasing the cost. In this case, we use an automatic nonuniform grid with a minimum steps per wavelength of 80.

For the boundary condition, we apply the periodic boundary condition in the \(x\) and \(y\) directions and PML in the \(z\) direction. symmetry=(1,-1,0) is used since the incident light is polarized in the \(y\) direction.

The default automatic shutoff level is 1e-5 but we decreased it to 1e-7 to further enhance the simulation accuracy, which is typically not required. It can be needed in simulations like plasmonics where strong fields are concentrated in relatively small volumes.

[7]:

run_time = 2e-13  # simulation run time

# define simulation domain size
sim_size = (a, np.sqrt(3) * a, 1.5 * lda0 + h)

# set up simulation
sim = td.Simulation(
    center=(0, 0, h / 2),
    size=sim_size,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=80, wavelength=lda0),
    structures=[dielectric, base, rod_1, rod_2, rod_3, rod_4, rod_5],
    sources=[plane_wave],
    monitors=[flux_monitor, monitor_field],
    run_time=run_time,
    boundary_spec=td.BoundarySpec(
        x=td.Boundary.periodic(), y=td.Boundary.periodic(), z=td.Boundary.pml()
    ),
    symmetry=(1, -1, 0),
    shutoff=1e-7,  # reducing the default shutoff level
)

After the simulation is set up, we can inspect the structure and grid to ensure the setup is correct and the grid is sufficiently fine.

[8]:

ax = sim.plot(z=h / 2)
sim.plot_grid(z=h / 2, ax=ax)
plt.show()

../_images/notebooks_PlasmonicNanorodArray_17_0.png

Also plot the 3D visualization.

[9]:

sim.plot_3d()

Run Simulation#

Now we are ready to submit the simulation to the server.

[10]:

sim_data = web.run(simulation=sim, task_name="nanorod_array")

10:13:08 Eastern Standard Time Created task 'nanorod_array' with task_id
                               'fdve-a32680fd-2fdc-47f2-b738-8854615b5f43' and
                               task_type 'FDTD'.

                               View task using web UI at
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-a32680fd-2fdc-47f2-b738-8854615b5f43'.

10:13:10 Eastern Standard Time status = queued

10:13:26 Eastern Standard Time status = preprocess

10:13:30 Eastern Standard Time Maximum FlexCredit cost: 0.031. Use
                               'web.real_cost(task_id)' to get the billed
                               FlexCredit cost after a simulation run.

                               starting up solver

                               running solver

                               To cancel the simulation, use
                               'web.abort(task_id)' or 'web.delete(task_id)' or
                               abort/delete the task in the web UI. Terminating
                               the Python script will not stop the job running
                               on the cloud.

10:13:52 Eastern Standard Time early shutoff detected at 64%, exiting.

                               status = postprocess

10:13:57 Eastern Standard Time status = success

                               View simulation result at
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-a32680fd-2fdc-47f2-b738-8854615b5f43'.

10:13:59 Eastern Standard Time loading simulation from simulation_data.hdf5

Postprocessing and Result Visualization#

Since the simulation domain is quite small, the simulation took less than a minute to finish. First, we plot the reflectance spectrum. Vertical lines are added to indicate the resonance wavelengths. At resonance, the incident wave excites the plasmon modes and thus the reflection is suppressed.

[11]:

R = sim_data["R"].flux  # extract reflectance

# plotting
plt.plot(ldas, R, c="black", linewidth=2)

# adding vertical lines at the resonances
for res_lda in res_ldas:
    plt.axvline(x=res_lda, linestyle="--", c="red")

plt.xlim(min(ldas), max(ldas))
plt.ylim(0, 1)
plt.xlabel("Wavelength($\mu m$)", fontsize=15)
plt.ylabel("Reflectance", fontsize=15)
plt.yticks(fontsize=15)
plt.xticks(fontsize=15)
plt.show()

../_images/notebooks_PlasmonicNanorodArray_25_0.png

Finally, plot the field profiles at the resonance wavelengths. In particular, we plot the absolute value of \(E_y\). Here we can clearly observe the standing wave patterns with increasing orders between the rods, consistent with the result presented in Figure 4 of the original publication.

[12]:

fig, ax = plt.subplots(1, len(res_ldas), figsize=(9, 4), tight_layout=True)
for i, res_freq in enumerate(res_freqs):
    E = np.sqrt(sim_data["field"].intensity.sel(f=res_freq))  # field strength

    # plotting
    E.plot(x="y", y="z", ax=ax[i], vmin=0, vmax=1e3, cmap="jet")
    ax[i].set_title(f"$\lambda$={res_ldas[i]} $\mu m$")
    ax[i].set_ylim(-0.2, h + 0.2)
plt.show()

../_images/notebooks_PlasmonicNanorodArray_27_0.png

Plasmonic cavity resonator consisting of a gold nanorod array

Contents

Plasmonic cavity resonator consisting of a gold nanorod array#

Simulation Setup#

Run Simulation#

Postprocessing and Result Visualization#

Normalization of the Field#