Scattering of a plasmonic nanoparticle¶
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.
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
# 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.
# 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.
# 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.
# 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]
Create Simulation¶
Now we can put everything together and define the simulation.
# 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.
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
sim.plot(y=0, ax=ax1)
sim.plot_grid(y=0, ax=ax1)
sim.plot(y=0, ax=ax2, monitor_alpha=0.2)
sim.plot_grid(y=0, ax=ax2)
zoom = radius * 2
ax2.axis([-zoom, zoom, -zoom, zoom])
plt.show()
Run Simulation¶
# Run simulation
sim_data = web.run(
sim,
task_name="plasmonic_nanoparticle",
path="data/plasmonic_nanoparticle.hdf5",
verbose=True,
)
Output()
Output()
Output()
Output()
Output()
Extract the results¶
# 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.
# 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()
