Gold nanoparticle validation
Contents
Gold nanoparticle validation#
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 both near and far fields are compared to the analytical Mie series.
To obtain the scattered field, two simulations will be run: * one with the sphere to compute the total field near the sphere, and * one without the sphere, to compute just the incident field and subtract it from the above to get the scattered field. To this end, we ensure that * the far field projection grid is identical in the two simulations, so that the total and incident far fields are sampled at the same location; * the FDTD mesh is identical in the two simulations, so that they exhibit similar grid dispersion and other numerical behavior.
Note that in the future, this process will be simplified with the introduction of a total-field-scattered-field source, requiring only one simulation. Stay tuned!
The far fields are computed by a near-to-far transformation of the near fields recorded on a closed surface around the sphere. A non-uniform mesh is carefully designed to make sure the sphere is well resolved without a significant sacrifice in efficiency.
[1]:
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
# tidy3d imports
import tidy3d as td
import tidy3d.web as web
[22:25:56] WARNING This version of Tidy3D was pip installed from the 'tidy3d-beta' repository on __init__.py:103 PyPI. Future releases will be uploaded to the 'tidy3d' repository. From now on, please use 'pip install tidy3d' instead.
INFO Using client version: 1.9.0rc1 __init__.py:121
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.
[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
# 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, so to ensure a sufficiently fine mesh in the particle, a phantom sphere
# is created with an extremely large permittivity, which is chosen so that there are at least
# ``cells_in_particle`` cells along the diameter of the sphere
cells_per_wavelength = 60
cells_in_particle = 60
phantom_electrical_diameter = cells_in_particle / cells_per_wavelength
wavelength_in_phantom = 2 * radius * phantom_electrical_diameter
wavelength_ratio = wavelength_in_phantom / wavelength_max
epsr_phantom = (1.0 / wavelength_ratio) ** 2
phantom = td.Structure(
geometry=td.Sphere(center=center, radius=radius),
medium=td.Medium(permittivity=epsr_phantom),
)
# 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 = [buffer_PML + 2 * radius + buffer_PML] * 3
# define PML layers on all sides
boundary_spec = td.BoundarySpec.all_sides(boundary=td.PML())
Create Source#
For our incident field, we create a plane wave incident from below the sphere polarized in the x direction.
[3]:
# Gaussian source offset; the source peak is at time t = offset/fwidth
offset = 4.0
# time dependence of source
gaussian = td.GaussianPulse(freq0=f0, fwidth=fwidth, offset=offset)
# place the source below the sphere, propagating in the +z direction
source = td.PlaneWave(
center=(0, 0, -(radius + 3 * buffer_PML / 4)),
size=(td.inf, td.inf, 0),
source_time=gaussian,
direction="+",
pol_angle=0,
)
# Simulation run time past the source decay (around t=2*offset/fwidth)
run_time = 10 / fwidth
Create Monitor#
Next, we define the monitors that will capture the near field data.
First, we create a FieldProjectionAngleMonitor which allows performing the far field transformation on GPUs during the solver run.
Also, we create a point FieldMonitor to record the field a point just above the sphere.
[4]:
# distance between the sphere and the near field monitor along each cartesian direction
buffer_mon = 3 * radius
# create a volume monitor around the sphere
mon_size = 2 * radius + 2 * buffer_mon
# set the observation angles
phis = [np.pi / 2]
thetas = [0]
# set the list of frequencies at which to compute far fields
num_freqs = 100
freqs = np.linspace(f0 - fwidth, f0 + fwidth, 100)
# create the near-to-far field monitor
monitor_n2f = td.FieldProjectionAngleMonitor(
center=center,
size=[mon_size] * 3,
freqs=list(freqs),
name="n2f",
phi=phis,
theta=thetas,
)
monitor_near = td.FieldMonitor(
center=[0, 0, 2 * radius],
size=[0, 0, 0],
freqs=list(freqs),
name="near",
)
monitor_flux = td.FluxMonitor(
center=[0, 0, 1.5 * radius],
size=[10 * radius, 10 * radius, 0],
freqs=list(freqs),
name="flux",
)
monitor_thru = td.FieldMonitor(
center=[0, 0, 0],
size=[td.inf, td.inf, 0],
freqs=list(freqs),
name="thru",
)
monitors = [monitor_n2f, monitor_near, monitor_flux, monitor_thru]
Create Simulation#
Now we can put everything together and define the two simulation classes: with the sphere to get the total field, and without the sphere to get the incident field. Since the incident field must be subtracted from the total field, we need to make sure the FDTD grid is identical in the two simulations.
[5]:
grid_spec = td.GridSpec.auto(
min_steps_per_wvl=cells_per_wavelength,
override_structures=[phantom],
)
grid_spec_empty = td.GridSpec.auto(
min_steps_per_wvl=cells_per_wavelength,
override_structures=[phantom],
)
# create the simulation with the particle
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-9,
)
# create the reference simulation without the particle
sim_empty = td.Simulation(
size=sim_size,
grid_spec=grid_spec_empty,
structures=[],
sources=[source],
monitors=monitors,
run_time=run_time,
boundary_spec=boundary_spec,
shutoff=1e-9,
)
Visualize Geometry and Mesh#
Let’s take a look at the mesh and geometry and make sure everything is defined properly in both simulations.
The first column shows the setup with the nanoparticle, while the second one shows the empty simulation but with an identical mesh.
[6]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim.plot(y=0, ax=ax1)
sim.plot_grid(y=0, ax=ax1)
sim_empty.plot(y=0, ax=ax2)
sim_empty.plot_grid(y=0, ax=ax2)
INFO Auto meshing using wavelength 0.5300 defined from sources. grid_spec.py:510
INFO Auto meshing using wavelength 0.5300 defined from sources. grid_spec.py:510
[6]:
<AxesSubplot: title={'center': 'cross section at y=0.00'}, xlabel='x', ylabel='z'>
The nanoparticle is very small compared to the simulation and is not visible in the plots above. For better visibility, we show a zoomed-in plot below, which shows how the nanoparticle is meshed.
[7]:
fig, ax = plt.subplots(1, 1, figsize=(4, 3))
sphere.plot(y=0, ax=ax)
sim.plot_grid(y=0, ax=ax)
ax.axis([-0.02, 0.02, -0.02, 0.02])
[7]:
(-0.02, 0.02, -0.02, 0.02)
Run Simulations#
Now we can run both simulations over time and measure the results
[8]:
# Run simulation
sim_data = web.run(
sim, task_name="plasmonic_nanoparticle", path="data/plasmonic_nanoparticle.hdf5"
)
sim_empty_data = web.run(
sim_empty,
task_name="plasmonic_nanoparticle_empty",
path="data/plasmonic_nanoparticle_empty.hdf5",
)
Compute Scattered Far Fields#
We now subtract the far fields associated with the empty simulation, which represents just the incident field, from the far fields associated with the full simulation, which represents total fields, to get scattered far fields.
[9]:
# compute scattered far fields by subtracting out the incident fields from the empty simulation
n2f_data = sim_data[monitor_n2f.name]
n2f_empty_data = sim_empty_data[monitor_n2f.name]
for field, field_empty in zip(
n2f_data.field_components.values(), n2f_empty_data.field_components.values()
):
field.values -= field_empty.values
Extract the results#
[10]:
# grab the far field data and compute the RCS to compare against Mie series
RCS = np.squeeze(n2f_data.radar_cross_section.values)
# compute the scattered near fields which we'll also compare against Mie series
Ex_scat = (
sim_data[monitor_near.name].Ex.values - sim_empty_data[monitor_near.name].Ex.values
)
Ey_scat = (
sim_data[monitor_near.name].Ey.values - sim_empty_data[monitor_near.name].Ey.values
)
Ez_scat = (
sim_data[monitor_near.name].Ez.values - sim_empty_data[monitor_near.name].Ez.values
)
Hx_scat = (
sim_data[monitor_near.name].Hx.values - sim_empty_data[monitor_near.name].Hx.values
)
Hy_scat = (
sim_data[monitor_near.name].Hy.values - sim_empty_data[monitor_near.name].Hy.values
)
Hz_scat = (
sim_data[monitor_near.name].Hz.values - sim_empty_data[monitor_near.name].Hz.values
)
E = np.sqrt(
np.squeeze(Ex_scat) ** 2 + np.squeeze(Ey_scat) ** 2 + np.squeeze(Ez_scat) ** 2
)
H = np.sqrt(
np.squeeze(Hx_scat) ** 2 + np.squeeze(Hy_scat) ** 2 + np.squeeze(Hz_scat) ** 2
)
flux = sim_data[monitor_flux.name].flux - sim_empty_data[monitor_flux.name].flux
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.
[11]:
plt.rcParams.update({"font.size": 14})
# load Mie series data
savefile_E = "./data/mie_plasmonic_Eabs.txt"
savefile_H = "./data/mie_plasmonic_Habs.txt"
savefile_RCS = "./data/mie_plasmonic_RCS.txt"
Eabs_mie = np.loadtxt(savefile_E, delimiter="\t", skiprows=1)[:, 1]
Habs_mie = np.loadtxt(savefile_H, delimiter="\t", skiprows=1)[:, 1]
RCS_mie = np.loadtxt(savefile_RCS, delimiter="\t", skiprows=1)[:, 1]
def to_db(val):
val = np.abs(val)
val = val / np.max(np.abs(val))
return 10.0 * np.log10(val)
def normalize(val):
return np.abs(val) / np.max(np.abs(val))
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
sim_data.plot_field(
field_monitor_name="thru", field_name="Ex", val="abs", f=f0, ax=ax[0]
)
sim_data.plot_field(
field_monitor_name="thru", field_name="Ey", val="abs", f=f0, ax=ax[1]
)
sim_data.plot_field(
field_monitor_name="thru", field_name="Ez", val="abs", f=f0, ax=ax[2]
)
for _ax in ax:
_ax.set_xlim([-20e-3, 20e-3])
_ax.set_ylim([-20e-3, 20e-3])
fig, ax = plt.subplots(figsize=(7.5, 5))
ax.plot(td.C_0 / freqs * 1e3, normalize(RCS_mie), "-k", label="Mie")
ax.plot(td.C_0 / freqs * 1e3, normalize(RCS), "--r", label="Tidy3D", mfc="None")
ax.set(
xlabel="Wavelength (nm)",
ylabel="Normalized radar cross section",
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=(7.5, 5))
ax.plot(td.C_0 / freqs * 1e3, normalize(Eabs_mie), "-k", label="Mie")
ax.plot(td.C_0 / freqs * 1e3, normalize(E), "--r", label="Tidy3D", mfc="None")
ax.set(
xlabel="Wavelength (nm)",
ylabel="Normalized near field |$\\vec{E}$|",
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=(7.5, 5))
ax.plot(td.C_0 / freqs * 1e3, normalize(Habs_mie), "-k", label="Mie")
ax.plot(td.C_0 / freqs * 1e3, normalize(H), "--r", label="Tidy3D", mfc="None")
ax.set(
xlabel="Wavelength (nm)",
ylabel="Normalized near field |$\\vec{H}$|",
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=(7.5, 5))
# ax.plot(td.C_0 / freqs * 1e3, flux_mie, '-k', label="Mie")
ax.plot(td.C_0 / freqs * 1e3, flux, "--r", label="Tidy3D", mfc="None")
ax.set(
xlabel="Wavelength (nm)",
ylabel="Flux",
yscale="linear",
xscale="linear",
)
ax.legend()
ax.grid(visible=True, which="both", axis="both", linewidth=0.4)
plt.tight_layout()
INFO Auto meshing using wavelength 0.5300 defined from sources. grid_spec.py:510
[ ]: