Electromagnetic metamaterials and metasurfaces are artificially engineered structures made of subwavelength resonating unit cells. By utilizing various design principles, they have enabled fascinating electromagnetic phenomena and capabilities such as negative refraction, cloaking, perfect absorption, high-NA focusing, and so on. In Tidy3D's example library, we have demonstrated a dielectric metasurface absorber, a metalens at the visible frequency, and a graphene metamaterial absorber.
This notebook demonstrates a near-infrared metasurface reflector consisting of plasmonic antennas. First, we model the unit cell and extract the relationship between a geometric parameter and the corresponding complex reflection coefficient. Then, a super cell consisting of ten antennas is designed to exhibit a linear reflection phase such that the incident light can be diffracted to a particular angle efficiently. By a similar procedure, other beam shaping devices and metalenses can be designed according to different phase profiles.
This example is based on the landmark work Sun, S.; Yang, K.-Y.; Wang, C.-M.; Juan, T.-K.; Chen, W.T.; Liao, C.Y.; He, Q.; Xiao, S.; Kung, W.-T.; Guo, G.-Y.; et al. High-Efficiency Broadband Anomalous Reflection by Gradient Meta-Surfaces. Nano Lett. 2012, 12, 6223–6229.
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
import tidy3d.web as web
Unit Cell Simulation¶
Simulation Setup¶
In the first part of the notebook, we simulate a unit cell that consists of a rectangular gold antenna, a dielectric spacer, and a gold film on glass substrate. The thickness of the gold film is d3=130 nm. The thickness of the spacer is d2 = 50 nm. The thickness and width of the antenna are d1=30 nm and w=90 nm. The only parameter we can vary is the length of the antenna L. The periodicities in x and y directions are Lx=120 nm and Ly=300 nm, respectively.
First, define the simulation frequency and wavelength ranges.
lda0 = 0.85 # central wavelength
freq0 = td.C_0 / lda0 # central frequency
Define geometric parameters.
Lx = 0.12 # unit cell size in x direction
Ly = 0.3 # unit cell size in y direction
d1 = 0.03 # antenna thickness
d2 = 0.05 # spacer thickness
d3 = 0.13 # gold layer thickness
w = 0.09 # antenna width
inf_eff = 1e2 # effective infinity
Three materials will be defined in this model: gold for the antennas and the back plate, magnesium fluoride for the spacer, and silica for the substrate. Magnesium fluoride and silica will be modeled as lossless and non-dispersive dielectrics while gold will be modeled as a dispersive medium.
More specifically, magnesium fluoride has a permittivity of 1.892 and silica has a permittivity of 2.25. For gold, we use one option from the Tidy3D's built-in material library. More specifically, we use the evaporated gold data from Olmon et al. There are other options for gold, as well as other common materials, in the material library.
# define MgF2 material for the spacer layer
mgf2 = td.Medium(permittivity=1.892)
# using material library gold refractive index
au = td.material_library["Au"]["Olmon2012evaporated"]
# define SiO2 material for the substrate
sio2 = td.Medium(permittivity=2.25)
Construct the substrate, gold film, spacer, and antenna structures.
# define SiO2 substrate
sub = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-inf_eff, -inf_eff, -inf_eff), rmax=(inf_eff, inf_eff, -d3)
),
medium=sio2,
)
# define gold layer
gold_layer = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-inf_eff, -inf_eff, -d3), rmax=(inf_eff, inf_eff, 0)
),
medium=au,
)
# define MgF2 spacer layer
spacer = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-inf_eff, -inf_eff, 0), rmax=(inf_eff, inf_eff, d2)
),
medium=mgf2,
)
Define a linearly polarized PlaneWave as the excitation source. The polarization is chosen to be in the y direction. Also define a DiffractionMonitor on the reflection side to extract the complex reflection coefficient.
fwidth = freq0 / 10 # width of the gaussian pulse
# define a plane wave excitation source
plane_wave = td.PlaneWave(
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
size=(td.inf, td.inf, 0),
center=(0, 0, 0.3 * lda0),
direction="-",
pol_angle=np.pi / 2, # polarization is set to y direction
)
# define a diffraction monitor to calculate the reflection coefficient
monitor_r = td.DiffractionMonitor(
center=[0, 0, 0.6 * lda0], size=[td.inf, td.inf, 0], freqs=[freq0], name="R"
)
To perform a parameter sweep over antenna length, it is convenient to define a function that builds the simulation with a given antenna length L. Later, this function will be called repeatedly to construct a simulation batch for parameter sweep.
Lz = 1.5 * lda0 # simulation domain size in z direction
sim_size = [Lx, Ly, Lz]
run_time = 3e-13 # simulation run time
boundary_spec = td.BoundarySpec(
x=td.Boundary.periodic(),
y=td.Boundary.periodic(),
z=td.Boundary(minus=td.PML(), plus=td.PML()),
)
# define a function to build simulation given antenna length L
def make_sim(L):
# define the gold antenna
antenna = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-w / 2, -L / 2, d2), rmax=(w / 2, L / 2, d2 + d1)
),
medium=au,
)
unit_cell = [sub, gold_layer, spacer, antenna]
# set up simulation
sim = td.Simulation(
size=sim_size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda0),
structures=unit_cell,
sources=[plane_wave],
monitors=[monitor_r],
run_time=run_time,
boundary_spec=boundary_spec, # pml is applied to z direction. x and y directions are periodic
)
return sim
To ensure the simulation setup is correct, let's test a specific case where the antenna length is 150 nm and visualize the simulation. From the plot, the geometry, source, and monitor all seem to be correctly defined.
# test a case where the antenna length is 150 nm
sim = make_sim(0.15)
# visualize the simulation
sim.plot(x=0)
plt.show()
To have a better visualization, we can also plot the simulation in 3D.
sim.plot_3d()
Now we are ready to perform the parameter sweep. In this particular case, we investigate antenna length from 40 nm to 280 nm. This is broken down into 14 simulations. A simulation batch is constructed by calling the make_sim function with different antenna lengths.
Ls = np.linspace(0.04, 0.28, 14) # antenna lengths for parameter sweep
sims = {f"L={L:.2f}": make_sim(L) for L in Ls} # construct simulation batch
Submit the simulation batch to the server.
# submit simulation batch to the server
batch = web.Batch(simulations=sims, verbose=True)
batch_results = batch.run(path_dir="data")
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- e9cb8b5f-8eaf-4a8d-853c-1b1f79bccd01v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- ea215c6c-27f8-4dae-b193-3c90dbdec55cv1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- c65925e3-464b-42bf-8323-3040b5a315b5v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- a546a4e8-3069-418b-990f-3dd273f2650bv1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- fa7a4d07-6cd1-4789-8751-162ecd3a0434v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- f82f063f-f7c8-4887-8fe5-f3b9cdf19db0v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 34f68355-13a7-4ca5-9814-d5e6dc0d3a00v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 72490fcf-240e-480c-80e4-a27a55ddf850v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 74c8edad-027a-4e99-8cb6-8a8bca309036v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- cc069453-3324-4783-824a-b3ba3e41d4fcv1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 2dec033c-1279-4529-ae28-01f965d10b6cv1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- a7150ddc-2dbb-4e87-933f-37dde47ae60av1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 1e6b582f-08e4-4c36-86b9-b8f2d4a7aff7v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 25341407-28a4-4923-bb2d-9fa9e1886344v1'.
[15:26:06] Started working on Batch. container.py:475
[15:26:11] Maximum FlexCredit cost: 0.350 for the whole batch. container.py:479 Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
[15:26:48] Batch complete. container.py:522
Result Analysis¶
After the simulations are complete, we are ready to extract the reflection coefficients. The reflection coefficient is simply the amplitude of the zero-th diffraction order.
# extract the reflection coefficient at each antenna length
r = np.zeros(len(Ls), dtype="complex")
for i, L in enumerate(Ls):
sim_data = batch_results[f"L={L:.2f}"]
r[i] = np.array(sim_data["R"].amps.sel(f=freq0, polarization="s"))[0][0]
To visualize the results, we plot the phase of the reflection coefficient and the reflectivity as a function of antenna length. Since the absolute value of the phase does not carry physical significance, we define 0 phase as the phase when the antenna is 40 nm long. The reflectivity decreases slightly when the antenna is 100 nm to 200 nm long but overall the reflectivity is above 80%.
# plot the reflection phase
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
theta = np.unwrap(np.angle(r)) * 180 / np.pi
theta = theta - theta[0]
ax1.plot(Ls, theta)
ax1.set_xlabel("L ($\mu m$)")
ax1.set_ylabel("Reflection phase (degree)")
# plot the reflectivity
ax2.plot(Ls, np.abs(r) ** 2)
ax2.set_xlabel("L ($\mu m$)")
ax2.set_ylabel("Reflectivity")
plt.show()
Super Cell Simulation¶
Simulation Setup¶
The above phase relation provides a look-up table such that a gradient metasurface with an arbitrary phase profile can be designed. For example, when designing a metalens, a parabolic phase profile can be used. In this example, we aim to construct a beam steering reflector. Therefore, a linear phase profile is used.
Here we adapt the idea from the referenced paper and design a super cell that consists of 10 antennas of 5 different lengths. The lengths are 40 nm, 106 nm, 128 nm, 150 nm, and 260 nm. The plot below confirms that this design indeed yields a nearly linearly varying phase across the super cell.
N_sample = 5 # number of antenna lengths in a super cell
N_repeat = 2 # number of repeated antenna lengths
Ls_design = np.array([0.04, 0.106, 0.128, 0.15, 0.26]) # designed antenna lengths
theta_design = np.interp(
Ls_design, Ls, theta
) # extract the reflection phase by interpolation the above simulation result
# plot phase across the super cell
plt.scatter(np.linspace(1, N_sample, N_sample), theta_design)
plt.xlabel("Site")
plt.ylabel("Reflection phase (degree)")
plt.show()
To define the geometry of the super cell, we first generate the antennas systematically.
# construct ten antennas given the designed antenna lengths
antennas = []
xs = np.linspace(
Lx / 2, (N_sample * N_repeat - 0.5) * Lx, N_sample * N_repeat
) # x coordinates of each antenna
# systematically generate antennas in the super cell
for i in range(N_sample):
for j in range(N_repeat):
antennas.append(
td.Structure(
geometry=td.Box.from_bounds(
rmin=(-w / 2 + xs[2 * i + j], -Ls_design[i] / 2, d2),
rmax=(w / 2 + xs[2 * i + j], Ls_design[i] / 2, d2 + d1),
),
medium=au,
)
)
super_cell = antennas + [spacer, sub, gold_layer]
To visualize the reflected field, we add a FieldMonitor in the xz plane. The simulation domain is lengthened in the z direction such that we can better visualize the reflected field.
# define a field monitor to visualize the reflected field
field_monitor = td.FieldMonitor(
center=(0, 0, 0),
size=(td.inf, 0, td.inf),
freqs=[freq0],
name="field",
)
# define the super cell simulation
sim = td.Simulation(
size=(Lx * N_sample * N_repeat, Ly, 3 * Lz),
center=(Lx * N_sample * N_repeat / 2, 0, 0),
grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda0),
structures=super_cell,
sources=[plane_wave],
monitors=[monitor_r, field_monitor],
run_time=run_time,
boundary_spec=boundary_spec,
)
[15:26:57] 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.
Visualize the generated antenna structures.
# visualize the antennas in the super cell
sim.plot(z=d2 + d1)
plt.show()
Visualize the simulation setup from the xz plane.
# visualize the xz plane
sim.plot(y=0)
plt.show()
To have a better visualization, we can also plot the simulation in 3D.
sim.plot_3d()
Submit the simulation to the server.
job = web.Job(simulation=sim, task_name="beam_steering_metasurface", verbose=True)
sim_data = job.run(path="data/simulation_data.hdf5")
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- f493f3f2-5a6e-4204-8ce4-b70611e4702av1'.
Result Analysis¶
After the simulation is complete, we first extract the reflected angles and power from the DiffractionMonitor. Then, plot the power and angle of each diffraction order as a scatter plot in polar coordinate.
theta = np.array(sim_data["R"].angles[0]) # diffraction angle theta
phi = np.array(sim_data["R"].angles[1]) # diffraction angle phi
power = np.array(sim_data["R"].power) # diffraction power of each order
# plot the power and angle in polar coordinate
fig, ax = plt.subplots(subplot_kw={"projection": "polar"})
plt.scatter(phi, theta * 180 / np.pi, s=60, c=power, vmin=0, vmax=1, cmap="bwr")
ax.set_rlim(0, 90)
plt.colorbar()
plt.show()
To see the power values more clearly, we plot the relationship between the power and the diffraction order.
plt.scatter(sim_data["R"].orders_x, sim_data["R"].power.values[:, 0])
plt.xlabel("Diffraction order")
plt.ylabel("Power")
plt.show()
The Above result shows that about 80% of the power is diffracted to about 45 degree. This is consistent with our design since we have 2π/Lx≈0.71k0, where Lx = 1200 nm is the super cell size in the x direction. This gives a diffraction angle of sin−1(0.71)≈45∘ at the first diffraction order. Note that another 5% of the power that is also diffracted to 45 degree corresponds to the -1 diffraction order.
Lastly, visualize the field distribution at the xz plane. The field above the super cell resembles a plane wave propagating at about 45∘. Clearly it is not a perfect plane wave. The distortion is due to the fact that 5% of the power is diffracted to the -1 order and about 1% to the 0 order. This design can potentially be further optimized.
# plot Ey distribution
fig, ax = plt.subplots()
Ey = sim_data["field"].Ey.sel(f=freq0).real
Ey.plot(x="x", y="z", ax=ax, vmin=-100, vmax=100, cmap="bwr")
ax.set_aspect("equal")
