Plasmonic nanoparticles can exhibit interesting electromagnetic properties at certain frequencies, such as a negative real part of the relative permittivty. Due to their high electrical conductivity, gold nanoparticles can be challenging to model: the rapid field variations inside and near the particles require a fine local discretization. Therefore, an intelligent non-uniform meshing scheme is essential to make sure that the particle is well resolved, while ensuring that the empty space around it does not lead to wasted simulation effort.

Scattering from a 10 nm gold sphere is modeled in this example, and results are compared to the analytical Mie series. This example uses the total-field-scattered-field (TFSF) source to excite the particle, and illustrates how to compute various quantities like the total absorbed power, total scattering cross-section, and forward cross-section. A non-uniform mesh is carefully designed to make sure the sphere is well resolved without a significant sacrifice in efficiency. For a more detailed demonstration on TFSF, please refer to the TFSF tutorial.

For more simulation examples, please visit our examples page. If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials. FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.

[1]:

# standard python imports
import matplotlib.pyplot as plt
import numpy as np

# tidy3d imports
import tidy3d as td
import tidy3d.web as web

Define the structure and boundary conditions¶

Note the special treatment in creating the mesh: we need to make sure that the mesh is sufficiently fine within the sphere, but we can also make use of Tidy3D's non-uniform meshing algorithm to have a coarser grid outside the sphere, for better efficiency. Furthermore, the TFSF source works best if placed in a region with uniform mesh. Because of this, we will add a mesh override structure which is slightly larger than the sphere itself.

[2]:

# radius and location of the nanoparticle
radius = 5e-3
center = [0, 0, 0]

# nanoparticle material
medium = td.material_library["Au"]["RakicLorentzDrude1998"]

# free space central wavelength of the pulse excitation
wavelength = 530e-3
f0 = td.C_0 / wavelength

# bandwidth in Hz
fwidth = f0 / 5.0
fmin = f0 - fwidth
fmax = f0 + fwidth
wavelength_max = td.C_0 / fmin
wavelength_min = td.C_0 / fmax

# distance between particle and the boundary of the tfsf box
buffer_tfsf = 0.3 * radius
tfsf_size = (radius + buffer_tfsf) * 2

# distance between the particle and the scattered field region monitors (should be larger than buffer_tfsf)
buffer_out = 0.4 * radius
out_size = (radius + buffer_out) * 2

# distance between the particle and the total field region monitor (should be smaller than buffer_tfsf)
buffer_in = 0.2 * radius
in_size = (radius + buffer_in) * 2

# distance between particle and the mesh override region
buffer_override = 0.8 * radius
override_size = (radius + buffer_override) * 2

# The nanoparticle is very electrically small in the frequency range considered here, and meshing
# based on a standard 10-30 points per wavelength would lead to a grid too coarse to resolve the
# curvature of the sphere. Instead, we just define how many cells we want within the sphere diameter
# by creating a mesh override region
cells_per_wavelength = 50  # for meshing outside of the particle
cells_in_particle = 60  # for the mesh override region
mesh_override = td.MeshOverrideStructure(
    geometry=td.Box(center=center, size=(override_size,) * 3),
    dl=(2 * radius / cells_in_particle,) * 3,
)

# create the sphere
sphere = td.Structure(geometry=td.Sphere(center=center, radius=radius), medium=medium)
geometry = [sphere]

# distance between the surface of the sphere and the start of the PML layers along each cartesian direction
buffer_pml = wavelength_max / 2

# set the full simulation size along x, y, and z
sim_size = [(radius + buffer_pml) * 2] * 3

# define PML layers on all sides
boundary_spec = td.BoundarySpec.all_sides(boundary=td.PML())

Create Source¶

For our incident field, we use a TFSF source to inject a plane wave incident from below the sphere polarized in the x direction.

[3]:

# time dependence of source
gaussian = td.GaussianPulse(freq0=f0, fwidth=fwidth)

# the tfsf source is defined as a box around the particle
source = td.TFSF(
    center=center,
    size=(tfsf_size,) * 3,
    source_time=gaussian,
    injection_axis=2,  # inject along the z axis...
    direction="+",  # ...in the positive direction, i.e. along z+
    name="tfsf",
    pol_angle=0,
)

# Simulation run time
run_time = 10 / fwidth

Create Monitors¶

Next, we define a number of monitors to analyze the simulation.

A 3D FluxMonitor inside the TFSF source can be used to compute the power absorbed in the particle. If no power is absorbed, the recorded flux would be 0, since there are no power sources or sinks inside the monitor box. Any flux recorded in the FluxMonitor is then due to the imbalance between the incoming power and the outgoing one, which has been reduced by the particle absorption.
A 3D FluxMonitor outside the TFSF region can be used to compute the total cross-section, meaning the total power of scattered radiation in all directions.
A FieldProjectionAngleMonitor can be used to compute the radiated cross-section at a particular angle.
Finally, we also store the fields in a planar cross-section for visualization purposes.

[4]:

# set the list of frequencies at which to compute quantities
num_freqs = 100
freqs = np.linspace(f0 - fwidth, f0 + fwidth, 100)

monitor_flux_out = td.FluxMonitor(
    size=[out_size] * 3,
    freqs=freqs,
    name="flux_out",
)
monitor_flux_in = td.FluxMonitor(
    size=[in_size] * 3,
    freqs=freqs,
    name="flux_in",
)

# Create near field to far field projection monitor for forward RCS
monitor_n2f = td.FieldProjectionAngleMonitor(
    center=center,
    size=[out_size] * 3,
    freqs=freqs,
    name="n2f",
    phi=[0],
    theta=[0],
)

# Near field monitor for visualization
monitor_near = td.FieldMonitor(
    center=[0, 0, 0],
    size=[out_size, 0, out_size],
    freqs=[f0],
    name="near",
)

monitors = [monitor_flux_out, monitor_flux_in, monitor_n2f, monitor_near]

[15:21:31] WARNING: Default value for the field monitor           monitor.py:261
           'colocate' setting has changed to 'True' in Tidy3D                   
           2.4.0. All field components will be colocated to the                 
           grid boundaries. Set to 'False' to get the raw fields                
           on the Yee grid instead.

Create Simulation¶

Now we can put everything together and define the simulation.

[5]:

# since the mesh is extremely fine near the particle and transitions to a coarser mesh farther out,
# we can use the 'max_scale' field to set the largest allowed increase in mesh size from one cell to
# the next, to obtain a gentler gradient; the default is 1.4, and we'll set it to 1.2 here
grid_spec = td.GridSpec.auto(
    min_steps_per_wvl=cells_per_wavelength,
    override_structures=[mesh_override],
    max_scale=1.2,
)

sim = td.Simulation(
    size=sim_size,
    grid_spec=grid_spec,
    structures=geometry,
    sources=[source],
    monitors=monitors,
    run_time=run_time,
    boundary_spec=boundary_spec,
    shutoff=1e-8,
)

Visualize Geometry and Mesh¶

Let's take a look at the mesh and geometry and make sure everything is defined properly in the simulation. The nanoparticle is very small compared to the simulation and is not visible in the full plot. For better visibility, we also show a zoomed-in plot, which shows how the nanoparticle is meshed. Note that the mesh inside the nanoparticle is quite fine and may not render correctly in the notebook, possibly giving the false impression of being somewhat coarse.

[6]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
zoom = radius * 2
sim.plot(y=0, ax=ax1)
sim.plot_grid(y=0, ax=ax1)
sim.plot(y=0, ax=ax2, hlim=[-zoom, zoom], vlim=[-zoom, zoom], monitor_alpha=0.2)
sim.plot_grid(y=0, ax=ax2, hlim=[-zoom, zoom], vlim=[-zoom, zoom])
plt.show()

No description has been provided for this image

Run Simulation¶

[7]:

# Run simulation
sim_data = web.run(
    sim,
    task_name="plasmonic_nanoparticle",
    path="data/plasmonic_nanoparticle.hdf5",
    verbose=True,
)

[15:21:34] Created task 'plasmonic_nanoparticle' with task_id      webapi.py:188
           'fdve-3c1702e5-fda0-478d-a403-0c684108529ev1'.

           View task using web UI at                               webapi.py:190
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-              
           3c1702e5-fda0-478d-a403-0c684108529ev1'.

[15:21:36] status = queued                                         webapi.py:361

[15:21:46] status = preprocess                                     webapi.py:355

[15:21:50] Maximum FlexCredit cost: 1.253. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.

           starting up solver                                      webapi.py:377

           running solver                                          webapi.py:386

           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           '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.

[15:27:02] early shutoff detected, exiting.                        webapi.py:404

           status = postprocess                                    webapi.py:419

[15:27:08] status = success                                        webapi.py:426

[15:27:09] loading SimulationData from                             webapi.py:590
           data/plasmonic_nanoparticle.hdf5

Extract the results¶

[8]:

# Total absorbed power
absorbed = np.abs(sim_data["flux_in"].flux)

# Total scattered power
scattered = sim_data["flux_out"].flux

# Power scattered in the forward direction
RCS = np.squeeze(np.real(sim_data["n2f"].radar_cross_section.values))

Plot Results and Compare with Mie Series¶

The final results are compared against the analytical Mie series below, and very good agreement is observed. The small deviations can be reduced with a further refinement of the grid. Since the sphere's material is dispersive, no subpixel averaging scheme is applied, so the simulation approximates the curved permittivity profile in a staircase-like manner, which contributes to the error.

The txt files containing the Mie results can be downloaded from our documentation repo.

[9]:

# load Mie series data
savefile_flux_abs = "./misc/mie_plasmonic_flux_abs.txt"
savefile_flux_scat = "./misc/mie_plasmonic_flux_scat.txt"
savefile_RCS = "./misc/mie_plasmonic_RCS.txt"

flux_abs_mie = np.loadtxt(savefile_flux_abs, delimiter="\t", skiprows=1)[:, 1]
flux_scat_mie = np.loadtxt(savefile_flux_scat, delimiter="\t", skiprows=1)[:, 1]
RCS_mie = np.loadtxt(savefile_RCS, delimiter="\t", skiprows=1)[:, 1]


def to_db(val):
    return val
    return 10.0 * np.log10(val)


fig, ax = plt.subplots(1, 3, figsize=(15, 3))
sim_data.plot_field(field_monitor_name="near", field_name="Ex", val="abs", f=f0, ax=ax[0])
sim_data.plot_field(field_monitor_name="near", field_name="Ey", val="abs", f=f0, ax=ax[1])
sim_data.plot_field(field_monitor_name="near", field_name="Ez", val="abs", f=f0, ax=ax[2])

fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(td.C_0 / freqs * 1e3, to_db(RCS_mie), "-k", label="Mie")
ax.plot(td.C_0 / freqs * 1e3, to_db(RCS), "--r", label="Tidy3D", mfc="None")
ax.set(
    xlabel="Wavelength (nm)",
    ylabel="Forward scattering (dBs$\\mu$m)",
    yscale="linear",
    xscale="linear",
)
ax.legend()
ax.grid(visible=True, which="both", axis="both", linewidth=0.4)
plt.tight_layout()

fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(td.C_0 / freqs * 1e3, flux_abs_mie, "-k", label="Mie")
ax.plot(td.C_0 / freqs * 1e3, absorbed, "--r", label="Tidy3D", mfc="None")
ax.set(
    xlabel="Wavelength (nm)",
    ylabel="Total absorbed power",
    yscale="linear",
    xscale="linear",
)
ax.legend()
ax.grid(visible=True, which="both", axis="both", linewidth=0.4)
plt.tight_layout()

fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(td.C_0 / freqs * 1e3, flux_scat_mie, "-k", label="Mie")
ax.plot(td.C_0 / freqs * 1e3, scattered, "--r", label="Tidy3D", mfc="None")
ax.set(
    xlabel="Wavelength (nm)",
    ylabel="Total scattered power",
    yscale="linear",
    xscale="linear",
)
ax.legend()
ax.grid(visible=True, which="both", axis="both", linewidth=0.4)
plt.tight_layout()
plt.show()

[ ]:

TIDY3D
LEARNING CENTER

Scattering of a plasmonic nanoparticle¶

Define the structure and boundary conditions¶

Create Source¶

Create Monitors¶

Create Simulation¶

Visualize Geometry and Mesh¶

Run Simulation¶

Extract the results¶

Plot Results and Compare with Mie Series¶

TIDY3D LEARNING CENTER

Scattering of a plasmonic nanoparticle¶

Define the structure and boundary conditions¶

Create Source¶

Create Monitors¶

Create Simulation¶

Visualize Geometry and Mesh¶

Run Simulation¶

Extract the results¶

Plot Results and Compare with Mie Series¶

TIDY3D
LEARNING CENTER