Author: Amir Targholizadeh, Syracuse University
In this notebook, we simulate a metalens focusing light onto superconducting nanowires for single-photon applications, as presented in the paper by Amir Targholizadeh, Grigoriy Y. Nikulin, Pankaj K. Jha, "Cryogenic Dielectric Metasurface-Integrated Superconducting Nanowire Single-Photon Detectors", (2025). DOI: https://doi.org/10.48550/arXiv.2512.17115.
The model consists of a plane wave source incident on the metalens, with the nanowire positioned above it. We place a FieldMonitor around the nanowire region and use the field information to calculate the absorbed power.
[1]:
import tidy3d as td
import matplotlib.pyplot as plt
import numpy as np
from tidy3d import web
Simulation Setup¶
Defining global parameters.
[2]:
NA_list = np.array([0.451])
distance = np.array([8])
z_designed_focus_list = np.array([17.33])
z_focus_list = np.array([13.10])
f_offset = z_designed_focus_list - z_focus_list
d = 8 # Distance
lda = 1.55
h_antenna = 1.96
L = 40
NA = 0.451
f00 = (L / 2) / np.tan(np.arcsin(NA)) # Designed focal length
nanowire_t = 0.008
cavity_h = 0.233
reflector_h = 0.0
sub_t = 0.1
n_sio2 = 1.457
n_silicon = 3.4757
nanowire_w = 0.1
nanowire_l = 5
nanowire_n = 25
f = f00 - h_antenna - f_offset # Actual focal length
L_h = 2 * lda / 3 + sub_t + h_antenna + f + cavity_h
sub_z = -L_h / 2 + 2 * lda / 3 + sub_t / 2
cavity_z = sub_z + sub_t / 2 + h_antenna + f + cavity_h / 2
f0 = td.C_0 / lda
reflector_z = cavity_z + cavity_h / 2 + reflector_h / 2
Defining metalens and nanowires.
[3]:
# ---------------------------------------------------------------------------------------
# Materials
sio2 = td.Medium(
name="sio2",
permittivity=2.1228490000000004,
)
silicon = td.Medium(
name="silicon",
permittivity=n_silicon**2,
)
# -----------------------------------------------------------------------------------------
# Structures
# metalens substrate
metalens_sub = td.Structure(
geometry=td.Box(center=[0, 0, sub_z], size=[td.inf, td.inf, sub_t]),
name="metalens_sub",
medium=sio2,
)
# metalens
cell_size = 0.5 # µm
grid_size = L # µm
num_cells = int(grid_size / cell_size)
# Generate center positions for each cell
x_centers = np.linspace(
-grid_size / 2 + cell_size / 2, grid_size / 2 - cell_size / 2, num_cells
)
y_centers = np.linspace(
-grid_size / 2 + cell_size / 2, grid_size / 2 - cell_size / 2, num_cells
)
# Nanoantenna dimensions
lx = np.array(
[
0,
0.238,
0.220,
0.220,
0.210,
0.206,
0.206,
0.273,
0.248,
0.248,
0.248,
0.234,
0.227,
0.227,
0.213,
0.298,
0.266,
0.266,
0.259,
0.345,
0.425,
0.302,
0.386,
0.309,
0.287,
0.298,
0.284,
0.386,
0.266,
0.450,
0.263,
0.355,
0.316,
0.252,
0.248,
0.425,
0.241,
0.280,
0.351,
0.294,
0.355,
0.337,
0.100,
0.100,
0.446,
0.294,
0.153,
0.128,
0.118,
0.160,
0.153,
0.153,
0.149,
0.142,
0.210,
0.408,
0.174,
0.195,
0.188,
0.188,
0.178,
0.171,
0.263,
0.224,
0.210,
]
)
ly = np.array(
[
0,
0.238,
0.259,
0.270,
0.302,
0.323,
0.337,
0.160,
0.195,
0.227,
0.248,
0.256,
0.294,
0.280,
0.333,
0.153,
0.192,
0.224,
0.234,
0.345,
0.234,
0.425,
0.100,
0.153,
0.189,
0.210,
0.231,
0.323,
0.266,
0.227,
0.298,
0.139,
0.178,
0.369,
0.390,
0.302,
0.439,
0.280,
0.418,
0.245,
0.171,
0.206,
0.369,
0.386,
0.319,
0.273,
0.146,
0.220,
0.263,
0.270,
0.302,
0.309,
0.316,
0.386,
0.132,
0.404,
0.224,
0.245,
0.291,
0.273,
0.319,
0.355,
0.118,
0.174,
0.210,
]
)
metalens_geometry = []
for x in x_centers:
for y in y_centers:
phi_x = np.mod(
(-2 * np.pi / lda * (np.sqrt((x - d) ** 2 + y**2 + f00**2) - f00))
* 180
/ np.pi,
360,
) # Phase profile for x_polarized light
phi_y = np.mod(
(-2 * np.pi / lda * (np.sqrt((x + d) ** 2 + y**2 + f00**2)) - f00)
* 180
/ np.pi,
360,
) # Phase profile for x_polarized light
if 0 < phi_x < 45:
if 0 < phi_y < 45:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[1], ly[1], h_antenna],
)
)
elif 45 < phi_y < 90:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[2], ly[2], h_antenna],
)
)
elif 90 < phi_y < 135:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[3], ly[3], h_antenna],
)
)
elif 135 < phi_y < 180:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[4], ly[4], h_antenna],
)
)
elif 180 < phi_y < 225:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[5], ly[5], h_antenna],
)
)
elif 225 < phi_y < 270:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[6], ly[6], h_antenna],
)
)
elif 270 < phi_y < 315:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[7], ly[7], h_antenna],
)
)
elif 315 < phi_y < 360:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[8], ly[8], h_antenna],
)
)
elif 45 < phi_x < 90:
if 0 < phi_y < 45:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[9], ly[9], h_antenna],
)
)
elif 45 < phi_y < 90:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[10], ly[10], h_antenna],
)
)
elif 90 < phi_y < 135:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[11], ly[11], h_antenna],
)
)
elif 135 < phi_y < 180:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[12], ly[12], h_antenna],
)
)
elif 180 < phi_y < 225:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[13], ly[13], h_antenna],
)
)
elif 225 < phi_y < 270:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[14], ly[14], h_antenna],
)
)
elif 270 < phi_y < 315:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[15], ly[15], h_antenna],
)
)
elif 315 < phi_y < 360:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[16], ly[16], h_antenna],
)
)
elif 90 < phi_x < 135:
if 0 < phi_y < 45:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[17], ly[17], h_antenna],
)
)
elif 45 < phi_y < 90:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[18], ly[18], h_antenna],
)
)
elif 90 < phi_y < 135:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[19], ly[19], h_antenna],
)
)
elif 135 < phi_y < 180:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[20], ly[20], h_antenna],
)
)
elif 180 < phi_y < 225:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[21], ly[21], h_antenna],
)
)
elif 225 < phi_y < 270:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[22], ly[22], h_antenna],
)
)
elif 270 < phi_y < 315:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[23], ly[23], h_antenna],
)
)
elif 315 < phi_y < 360:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[24], ly[24], h_antenna],
)
)
elif 135 < phi_x < 180:
if 0 < phi_y < 45:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[25], ly[25], h_antenna],
)
)
elif 45 < phi_y < 90:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[26], ly[26], h_antenna],
)
)
elif 90 < phi_y < 135:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[27], ly[27], h_antenna],
)
)
elif 135 < phi_y < 180:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[28], ly[28], h_antenna],
)
)
elif 180 < phi_y < 225:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[29], ly[29], h_antenna],
)
)
elif 225 < phi_y < 270:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[30], ly[30], h_antenna],
)
)
elif 270 < phi_y < 315:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[31], ly[31], h_antenna],
)
)
elif 315 < phi_y < 360:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[32], ly[32], h_antenna],
)
)
elif 180 < phi_x < 225:
if 0 < phi_y < 45:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[33], ly[33], h_antenna],
)
)
elif 45 < phi_y < 90:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[34], ly[34], h_antenna],
)
)
elif 90 < phi_y < 135:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[35], ly[35], h_antenna],
)
)
elif 135 < phi_y < 180:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[36], ly[36], h_antenna],
)
)
elif 180 < phi_y < 225:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[37], ly[37], h_antenna],
)
)
elif 225 < phi_y < 270:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[38], ly[38], h_antenna],
)
)
elif 270 < phi_y < 315:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[39], ly[39], h_antenna],
)
)
elif 315 < phi_y < 360:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[40], ly[40], h_antenna],
)
)
elif 225 < phi_x < 270:
if 0 < phi_y < 45:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[41], ly[41], h_antenna],
)
)
elif 45 < phi_y < 90:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[42], ly[42], h_antenna],
)
)
elif 90 < phi_y < 135:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[43], ly[43], h_antenna],
)
)
elif 135 < phi_y < 180:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[44], ly[44], h_antenna],
)
)
elif 180 < phi_y < 225:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[45], ly[45], h_antenna],
)
)
elif 225 < phi_y < 270:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[46], ly[46], h_antenna],
)
)
elif 270 < phi_y < 315:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[47], ly[47], h_antenna],
)
)
elif 315 < phi_y < 360:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[48], ly[48], h_antenna],
)
)
elif 270 < phi_x < 315:
if 0 < phi_y < 45:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[49], ly[49], h_antenna],
)
)
elif 45 < phi_y < 90:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[50], ly[50], h_antenna],
)
)
elif 90 < phi_y < 135:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[51], ly[51], h_antenna],
)
)
elif 135 < phi_y < 180:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[52], ly[52], h_antenna],
)
)
elif 180 < phi_y < 225:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[53], ly[53], h_antenna],
)
)
elif 225 < phi_y < 270:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[54], ly[54], h_antenna],
)
)
elif 270 < phi_y < 315:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[55], ly[55], h_antenna],
)
)
elif 315 < phi_y < 360:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[56], ly[56], h_antenna],
)
)
elif 315 < phi_x < 360:
if 0 < phi_y < 45:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[57], ly[57], h_antenna],
)
)
elif 45 < phi_y < 90:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[58], ly[58], h_antenna],
)
)
elif 90 < phi_y < 135:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[59], ly[59], h_antenna],
)
)
elif 135 < phi_y < 180:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[60], ly[60], h_antenna],
)
)
elif 180 < phi_y < 225:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[61], ly[61], h_antenna],
)
)
elif 225 < phi_y < 270:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[62], ly[62], h_antenna],
)
)
elif 270 < phi_y < 315:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[63], ly[63], h_antenna],
)
)
elif 315 < phi_y < 360:
metalens_geometry.append(
td.Box(
center=[x, y, sub_z + sub_t / 2 + h_antenna / 2],
size=[lx[64], ly[64], h_antenna],
)
)
metalens = td.Structure(
geometry=td.GeometryGroup(geometries=metalens_geometry), medium=silicon
)
# Cavity
cavity = td.Structure(
# geometry = td.Box(center = [0, 0, cavity_z], size = [td.inf, td.inf, cavity_h]),
geometry=td.Box.from_bounds(
rmin=[-1e5, -1e5, cavity_z - cavity_h / 2], rmax=[1e5, 1e5, 1e5]
),
name="cavity",
medium=sio2,
)
# nanowires_x
nanowire_z = cavity_z - cavity_h / 2 - nanowire_t / 2
nanowire_y_x = -(2 * nanowire_n - 1) * nanowire_w / 2 + nanowire_w / 2
NbTiN = td.Medium(
name="NbTiN", permittivity=2.817499999999999, conductivity=0.8897831601063977
)
xnanowire = []
for n in range(0, nanowire_n):
xnanowire.append(
td.Box(
center=[d, nanowire_y_x + 2 * nanowire_w * n, nanowire_z],
size=[nanowire_l, nanowire_w, nanowire_t],
)
)
x_meander = td.Structure(geometry=td.GeometryGroup(geometries=xnanowire), medium=NbTiN)
# nanowires_y
nanowire_y_y = -d - (2 * nanowire_n - 1) * nanowire_w / 2 + nanowire_w / 2
ynanowire = []
for n in range(0, nanowire_n):
ynanowire.append(
td.Box(
center=[nanowire_y_y + 2 * nanowire_w * n, 0, nanowire_z],
size=[nanowire_w, nanowire_l, nanowire_t],
)
)
y_meander = td.Structure(geometry=td.GeometryGroup(geometries=ynanowire), medium=NbTiN)
Source and monitors.
[4]:
# ---------------------------------------------------------------------------------------------
# Sources
x_pol = td.PlaneWave(
name="x_pol",
center=[0, 0, sub_z - sub_t],
size=[td.inf, td.inf, 0],
source_time=td.GaussianPulse(
freq0=f0,
fwidth=f0 / 5,
),
num_freqs=1,
direction="+",
)
y_pol = td.PlaneWave(
name="y_pol",
center=[0, 0, sub_z - sub_t],
size=[td.inf, td.inf, 0],
source_time=td.GaussianPulse(
phase=3.141592653589793,
freq0=f0,
fwidth=f0 / 5,
),
num_freqs=1,
direction="+",
pol_angle=np.pi / 2,
)
# ---------------------------------------------------------------------------------------------
# Monitors
fieldmonitor_0 = td.FieldMonitor(
name="fieldmonitor_0",
size=[td.inf, 0, td.inf],
freqs=td.C_0 / 1.5500000000000003,
fields=["Ex", "Ey", "Ez"],
)
fieldmonitor_1 = td.FieldMonitor(
name="fieldmonitor_1",
center=[0, 0, nanowire_z],
size=[td.inf, td.inf, nanowire_t],
freqs=td.C_0 / 1.5500000000000003,
)
fieldmonitor_2 = td.FieldMonitor(
name="fieldmonitor_2",
center=[0, 0, sub_z - 0.8 * sub_t],
size=[td.inf, td.inf, 0],
freqs=td.C_0 / 1.5500000000000003,
)
# ---------------------------------------------------------
# Meshing for nanowires
refine_box_x = td.MeshOverrideStructure(
geometry=td.Box(
center=[d, 0, nanowire_z],
size=[nanowire_l, (2 * nanowire_n - 1) * nanowire_w, nanowire_t],
),
dl=[None, 0.02, 0.002],
)
refine_box_y = td.MeshOverrideStructure(
geometry=td.Box(
center=[-d, 0, nanowire_z],
size=[nanowire_l, (2 * nanowire_n - 1) * nanowire_w, nanowire_t],
),
dl=[0.02, None, 0.002],
)
Defining the simulation object. To prevent evanescent fields from reaching the PML, which can cause divergence, we add an extra 2 µm to the xy dimensions.
[5]:
# -----------------------------------------------------------------------
# Boundary condition
boundary_z = td.Boundary(minus=td.PML(), plus=td.PECBoundary())
sim = td.Simulation(
size=(L + 2, L + 2, L_h),
grid_spec=td.GridSpec.auto(
wavelength=1.55,
min_steps_per_wvl=6,
override_structures=[refine_box_x, refine_box_y],
),
structures=[metalens, metalens_sub, cavity, x_meander, y_meander],
boundary_spec=td.BoundarySpec(
x=td.Boundary.pml(), z=boundary_z, y=td.Boundary.pml()
),
sources=[x_pol],
monitors=[fieldmonitor_1, fieldmonitor_0, fieldmonitor_2],
run_time=1 * (f00 - f_offset) * 1e-14,
)
Visualizing the simulation.
[6]:
# Plot grid around nanowires
ax = sim.plot(z=ynanowire[0].center[2])
sim.plot_grid(ax=ax, z=ynanowire[0].center[2])
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
plt.show()
[7]:
# Visualize metalens structure
ax = sim.plot(z=metalens_geometry[0].center[2])
# sim.plot_grid(ax=ax,z=metalens_geometry[0].center[2])
plt.show()
[8]:
sim.plot_3d()
[9]:
# Run simulation
sim_data = web.run(sim, task_name="L 30 um", path="data/data.hdf5", verbose=True)
10:38:21 -03 Created task 'L 30 um' with resource_id 'fdve-eaea4719-e0f6-4b52-8d5b-b575ae49a8bc' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-eaea4719-e0f 6-4b52-8d5b-b575ae49a8bc'.
Task folder: 'default'.
Output()
10:38:27 -03 Estimated FlexCredit cost: 2.647. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
10:38:29 -03 status = success
Output()
10:38:46 -03 Loading simulation from data/data.hdf5
10:38:47 -03 WARNING: Simulation final field decay value of 0.102 is greater than the simulation shutoff threshold of 1e-05. Consider running the simulation again with a larger 'run_time' duration for more accurate results.
We can note that the simulation fields did not decay completely by the end of the simulation. However, the runtime was chosen so that the fields decay at the nanowire positions, while the reflected fields do not reach the -z boundaries.
Now, we will double-check the field normalization by integrating the approximate power flux, estimated from the electric field magnitude as
$P = \iint_S \frac{|\mathbf{E}|^2}{2\eta_0} \, dA$,
which corresponds to the plane-wave approximation of the Poynting vector magnitude,
$|\mathbf{S}| = \frac{|\mathbf{E}|^2}{2\eta_0}$.
[10]:
# Calculating input power
Ex = sim_data["fieldmonitor_2"].Ex.abs.isel(f=0).squeeze(drop=True)
Ey = sim_data["fieldmonitor_2"].Ey.abs.isel(f=0).squeeze(drop=True)
Ez = sim_data["fieldmonitor_2"].Ez.abs.isel(f=0).squeeze(drop=True)
I = Ex**2 + Ey**2 + Ez**2
impdance = 377
power = I / (2 * impdance)
x = power["x"].values
y = power["y"].values
domain = power.sel(x=slice(-L / 2, L / 2), y=slice(-L / 2, L / 2))
Pt = domain.integrate(coord="x").integrate(coord="y") # total input power
print(f"Normalization factor: {Pt:.2f}")
Normalization factor: 0.74
Next, we calculate the total absorbed power at the nanowires using:
$P_{abs} = \tfrac{1}{2} \sigma |E|^2$
[11]:
# --------------------------------------------------------------------------------------
# loss calculation in x meander
Ex = sim_data["fieldmonitor_1"].Ex.abs.isel(f=0).squeeze(drop=True)
Ey = sim_data["fieldmonitor_1"].Ey.abs.isel(f=0).squeeze(drop=True)
Ez = sim_data["fieldmonitor_1"].Ez.abs.isel(f=0).squeeze(drop=True)
I = Ex**2 + Ey**2 + Ez**2
conductivity = 0.8897831601063977
loss = I * conductivity / 2
x = loss["x"].values
y = loss["y"].values
z = loss["z"].values
offset = 0.001
Px = 0
Px_vals = []
for n in range(0, 25):
x = loss["x"].values
y = loss["y"].values
z = loss["z"].values
x_max, x_min = d + nanowire_l / 2, d - nanowire_l / 2
y_max, y_min = (
nanowire_y_x + nanowire_w / 2 + 2 * nanowire_w * n + offset,
nanowire_y_x - nanowire_w / 2 + 2 * nanowire_w * n - offset,
)
z_max, z_min = (
nanowire_z + nanowire_t / 2 + offset,
nanowire_z - nanowire_t / 2 - offset,
)
domain = loss.sel(
x=slice(x_min, x_max),
y=slice(y_min, y_max),
z=slice(z_min.item(), z_max.item()),
)
Px = domain.integrate(coord="x").integrate(coord="y").integrate(coord="z") + Px
print(f"Fraction of power absorbed in x meander: {Px / Pt:.2f}")
P = Px / Pt
Px_vals.append(P)
# -------------------------------------------------------------------------------------------
# loss calculation in y meander
Py_vals = []
Py = 0
for n in range(0, nanowire_n):
x = loss["x"].values
y = loss["y"].values
z = loss["z"].values
x_max, x_min = (
nanowire_y_y + nanowire_w / 2 + 2 * nanowire_w * n + offset,
nanowire_y_y - nanowire_w / 2 + 2 * nanowire_w * n - offset,
)
y_max, y_min = nanowire_l / 2, -nanowire_l / 2
z_max, z_min = (
nanowire_z + nanowire_t / 2 + offset,
nanowire_z - nanowire_t / 2 - offset,
)
domain = loss.sel(
x=slice(x_min, x_max),
y=slice(y_min, y_max),
z=slice(z_min.item(), z_max.item()),
)
Py = domain.integrate(coord="x").integrate(coord="y").integrate(coord="z") + Py
print(f"Fraction of power absorbed in y meander: {Py / Pt:.2f}")
P = Py / Pt
Py_vals.append(P)
Fraction of power absorbed in x meander: 0.47 Fraction of power absorbed in y meander: 0.01
Plotting the results.
[12]:
sim_data.plot_field("fieldmonitor_0", "E", val="abs^2", vmin=0, vmax=100)
plt.show()
[13]:
s = sim_data._get_scalar_field(
field_monitor_name="fieldmonitor_1", field_name="E", val="abs"
)
s2 = s**2
s2d = s2.sel(z=0, method="nearest").sel(f=193414489032258.03) # take a z-slice at z=0
plt.figure()
plt.pcolormesh(
s2d["x"], s2d["y"], s2d.T, cmap="magma", vmin=0, vmax=100, shading="auto"
)
plt.colorbar(label="|E|^2")
plt.xlabel("x (µm)")
plt.ylabel("y (µm)")
plt.tight_layout() # avoid label cutoff
plt.show()
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: