Large-area metalens¶
Introduction¶
See on github, run on colab, or just follow along with the output below.
Here we use Tidy3D to simulate a very large dielectric metalens. We base this example off of the recent paper from Khorasaninejad et al. titled Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging, which was recently published in Science. In this paper, a 2-dimensional array of dielectric structures is used as a lens to focus transmitted light to a single position directly above the device.
Typically, these structures are simulated by simulating each dielectric unit cell individually to compute a phase and amplitude transmittance for each cell. While this approach makes for an approximation of the overall device performance, it would be useful to be able to simulate the entire device as a whole to capture the entire physics. However, a simulation of this scale requires several hours or days to do with a conventional CPU-based FDTD. With Tidy3D, we are able to complete the entire simulation in about 1 minute!
Setup¶
We first perform basic imports of the packages needed.
[1]:
# get the most recent version of tidy3d
!pip install -q --upgrade tidy3d
# make sure notebook plots inline
%matplotlib inline
# standard python imports
import numpy as np
from numpy import random
import matplotlib.pyplot as plt
# tidy3D import
import tidy3d as td
from tidy3d import web
Define Simulation Parameters¶
We now set the basic parameters that define the metalens. The following image (taken from the original paper) define the variables describing the unit cell of the metalens. The angle of each cell (θ) is chosen for maximum focusing effect. Note that microns are the default spatial unit in Tidy3D.
[2]:
# 1 nanometer in units of microns (for conversion)
nm = 1e-3
# free space central wavelength
wavelength = 600 * nm
# desired numerical aperture
NA = 0.5
# shape parameters of metalens unit cell (um) (refer to image above and see paper for details)
W = 85 * nm
L = 410 * nm
H = 600 * nm
S = 430 * nm
# space between bottom PML and substrate (-z)
space_below_sub = 1 * wavelength
# thickness of substrate
thickness_sub = 100 * nm
# side length of entire metalens (um)
side_length = 10
# Number of unit cells in each x and y direction (NxN grid)
N = int(side_length / S)
print(f'for diameter of {side_length:.1f} um, have {N} cells per side')
print(f'full metalens has area of {side_length**2:.1f} um^2 and {N*N} total cells')
# Define material properties at 600 nm
n_TiO2 = 2.40
n_SiO2 = 1.46
air = td.Medium(epsilon=1.0)
SiO2 = td.Medium(epsilon=n_SiO2**2)
TiO2 = td.Medium(epsilon=n_TiO2**2)
# resolution control
grids_per_wavelength = 25
# Number of PML layers to use along z direction
npml = 10
for diameter of 10.0 um, have 23 cells per side
full metalens has area of 100.0 um^2 and 529 total cells
Process Geometry¶
Next we need to do conversions to get the problem parameters ready to define the simulation.
[3]:
# grid size (um)
dl = wavelength / grids_per_wavelength
# using the wavelength in microns, one can use td.C_0 (um/s) to get frequency in Hz
# wavelength_meters = wavelength * meters
f0 = td.C_0 / wavelength
# Define PML layers, for this we have no PML in x, y but `npml` cells in z
pml_layers = [0, 0, npml]
# Compute the domain size in x, y (note: round down from side_length)
length_xy = N * S
# focal length given diameter and numerical aperture
f = length_xy / 2 / NA * np.sqrt(1 - NA**2)
# Function describing the theoretical best angle of each box at position (x,y). see paper for details
def theta(x, y):
return np.pi / wavelength * (f - np.sqrt(x ** 2 + y ** 2 + f ** 2))
# total domain size in z: (space -> substrate -> unit cell -> 1.7 focal lengths)
length_z = space_below_sub + thickness_sub + H + 1.7 * f
# construct simulation size array
sim_size = np.array([length_xy, length_xy, length_z])
Create Metalens Geometry¶
Now we can automatically generate a large metalens structure using these parameters.
We will first create the substrate as a td.Box
Then, we will loop through the x and y coordinates of the lens and create each unit cell as a td.PolySlab
[4]:
# define substrate
substrate = td.Box(
center=[0, 0, -length_z/2 + space_below_sub + thickness_sub / 2.0],
size=[length_xy, length_xy, thickness_sub],
material=SiO2)
# create a running list of structures
geometry = [substrate]
# define coordinates of each unit cell
centers_x = S * np.arange(N) - length_xy / 2.0 + S / 2.0
centers_y = S * np.arange(N) - length_xy / 2.0 + S / 2.0
center_z = -length_z/2 + space_below_sub + thickness_sub + H / 2.0
# convenience function to make an angled box at each x,y location using polyslab.
# For more details see, https://simulation.cloud/docs/html/generated/tidy3d.PolySlab.html
def angled_box(x, y, angle):
""" make a box of size (L, W, H) centered at `(x, y)` at `angle` from x axis"""
# x, y vertices of box of size (L, W) centered at the origin
vertices_origin = np.array([[+L/2, +W/2],
[-L/2, +W/2],
[-L/2, -W/2],
[+L/2, -W/2]])
# 2x2 rotation matrix angle `angle` with respect to x axis
rotation_matrix = np.array([[+np.cos(angle), -np.sin(angle)],
[+np.sin(angle), +np.cos(angle)]])
# rotate the origin vertices by this angle
vertices_rotated = vertices_origin @ rotation_matrix
# shift the rotated vertices to be centered at (x, y)
vertices = vertices_rotated + np.array([x, y])
# create a tidy3D PolySlab with these rotated and shifted vertices and thickness `H`
return td.PolySlab(
vertices=vertices,
z_cent=center_z,
z_size=H,
material=TiO2
)
# loop through the coordinates and add all unit cells to geometry list
for i, x in enumerate(centers_x):
for j, y in enumerate(centers_y):
angle = theta(x, y)
geometry.append(angled_box(x, y, angle))
Define Sources¶
Now we define the incident fields. We simply use a normally incident plane wave with Guassian time dependence centered at our central frequency. For more details, see the plane wave source documentation and the gaussian source documentation
[5]:
# Bandwidth in Hz
fwidth = f0 / 10.0
# time dependence of source
gaussian = td.GaussianPulse(f0, fwidth, phase=0)
source = td.PlaneWave(
source_time=gaussian,
injection_axis='+z',
position=-length_z/2 + space_below_sub / 10.0, # just next to PML
polarization='x')
# Simulation run time past the source decay (around t=2*offset/fwidth)
run_time = 40 / fwidth
Define Monitors¶
Now we define the monitors that measure field output from the FDTD simulation. For simplicity, we use measure the fields at the central frequency. We’ll get the fields at the following locations:
The y=0 cross section
The x=0 cross section
The z=f cross section at the focal plane
The central axis, along x=y=0
For more details on defining monitors, see the frequency-domain monitor documentation.
[6]:
# get fields along x=y=0 axis
monitor_center = td.FreqMonitor(
center=[0., 0., 0],
size=[0, 0, length_z],
freqs=[f0],
store=['E'],
name='central_line')
# get the fields at a few cross-sectional planes
monitor_xz = td.FreqMonitor(
center=[0., 0., 0.],
size=[length_xy, 0., length_z],
freqs=[f0],
store=['E'],
name='xz_plane')
monitor_yz = td.FreqMonitor(
center=[0., 0., 0.],
size=[0., length_xy, length_z],
freqs=[f0],
store=['E'],
name='yz_plane')
monitor_xy = td.FreqMonitor(
center=[0., 0., center_z + H/2 + f],
size=[length_xy, length_xy, 0],
freqs=[f0],
store=['E'],
name='focal_plane')
# put them into a single list
monitors=[monitor_center, monitor_xz, monitor_yz, monitor_xy]
Create Simulation¶
Now we can put everything together and define a simulation class to be run
[7]:
sim = td.Simulation(
size=sim_size,
mesh_step=[dl, dl, dl],
structures=geometry,
sources=[source],
monitors=monitors,
run_time=run_time,
pml_layers=pml_layers)
Initializing simulation...
Mesh step (micron): [2.40e-02, 2.40e-02, 2.40e-02].
Simulation domain in number of grid points: [412, 412, 680].
Total number of computational grid points: 1.15e+08.
Total number of time steps: 19246.
Estimated data size (GB) of monitor central_line: 0.0000.
Estimated data size (GB) of monitor xz_plane: 0.0065.
Estimated data size (GB) of monitor yz_plane: 0.0065.
Estimated data size (GB) of monitor focal_plane: 0.0041.
Visualize Geometry¶
Lets take a look and make sure everything is defined properly
[8]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 6))
# Time the visualization of the 2D plane
sim.viz_eps_2D(normal='x', position=0.1, ax=ax1);
sim.viz_eps_2D(normal='y', position=0.1, ax=ax2);
sim.viz_eps_2D(normal='z', position=-length_z/2 + space_below_sub + thickness_sub + H / 2, ax=ax3);
[8]:
<AxesSubplot:title={'center':'xy-plane at z=-6.93e+00um'}, xlabel='x (um)', ylabel='y (um)'>

Run Simulation¶
Now we can run the simulation over time and measure the results to plot
[9]:
# Run simulation
project = web.new_project(sim.export(), task_name='metalens')
web.monitor_project(project['taskId'])
# download and load the results
print('Downloading results')
web.download_results(project['taskId'], target_folder='output')
sim.load_results('output/monitor_data.hdf5')
# print stats from the logs
with open("output/tidy3d.log") as file:
print(file.read())
Uploading the json file...
Project 'metalens-16b875e7361b4fb0-16b875e785c961f0' status: success...
Downloading results
Applying source normalization to all frequency monitors using source index 0.
Simulation domain Nx, Ny, Nz: [412, 412, 680]
Applied symmetries: [0, 0, 0]
Number of computational grid points: 1.1543e+08.
Using subpixel averaging: True
Number of time steps: 19246
Automatic shutoff factor: 1.00e-05
Time step (s): 4.1598e-17
Compute source modes time (s): 0.2846
Compute monitor modes time (s): 0.2866
Rest of setup time (s): 0.3859
Starting solver...
- Time step 383 / time 1.59e-14s ( 1 % done), field decay: 1.00e+00
- Time step 769 / time 3.20e-14s ( 4 % done), field decay: 1.00e+00
- Time step 1539 / time 6.40e-14s ( 8 % done), field decay: 9.49e-01
- Time step 2309 / time 9.60e-14s ( 12 % done), field decay: 1.35e-01
- Time step 3079 / time 1.28e-13s ( 16 % done), field decay: 7.74e-02
- Time step 3849 / time 1.60e-13s ( 20 % done), field decay: 4.49e-02
- Time step 4619 / time 1.92e-13s ( 24 % done), field decay: 2.74e-02
- Time step 5388 / time 2.24e-13s ( 28 % done), field decay: 1.70e-02
- Time step 6158 / time 2.56e-13s ( 32 % done), field decay: 1.11e-02
- Time step 6928 / time 2.88e-13s ( 36 % done), field decay: 7.65e-03
- Time step 7698 / time 3.20e-13s ( 40 % done), field decay: 5.43e-03
- Time step 8468 / time 3.52e-13s ( 44 % done), field decay: 4.05e-03
- Time step 9238 / time 3.84e-13s ( 48 % done), field decay: 3.09e-03
- Time step 10007 / time 4.16e-13s ( 52 % done), field decay: 2.44e-03
- Time step 10777 / time 4.48e-13s ( 56 % done), field decay: 1.97e-03
- Time step 11547 / time 4.80e-13s ( 60 % done), field decay: 1.59e-03
- Time step 12317 / time 5.12e-13s ( 64 % done), field decay: 1.32e-03
- Time step 13087 / time 5.44e-13s ( 68 % done), field decay: 1.11e-03
- Time step 13857 / time 5.76e-13s ( 72 % done), field decay: 9.43e-04
- Time step 14626 / time 6.08e-13s ( 76 % done), field decay: 8.10e-04
- Time step 15396 / time 6.40e-13s ( 80 % done), field decay: 6.93e-04
- Time step 16166 / time 6.72e-13s ( 84 % done), field decay: 6.03e-04
- Time step 16936 / time 7.05e-13s ( 88 % done), field decay: 5.21e-04
- Time step 17706 / time 7.37e-13s ( 92 % done), field decay: 4.58e-04
- Time step 18476 / time 7.69e-13s ( 96 % done), field decay: 4.08e-04
- Time step 19245 / time 8.01e-13s (100 % done), field decay: 3.51e-04
Solver time (s): 49.9682
Post-processing time (s): 0.2859
As we can see from the logs, the total time to solve this problem (not including data transfer and pre/post-processing) was about 1 minute!
For reference, the same problem run with FDTD on a single CPU core FDTD is projected to take 11 hours!
Visualize Fields¶
Let’s see the results of the simulation as captured by our field monitors.
First, we look at the field intensity along the axis of the lens
[10]:
# split monitors by those that have area and those that are on a line
line_monitor, *area_monitors = monitors
# get the data from the line monitor and plot it
data = sim.data(line_monitor)
E = data['E']
zs = data['zmesh']
I = np.squeeze(np.sum(np.square(np.abs(E)), axis=0))
plt.plot(zs / wavelength, I / np.max(np.abs(I)))
plt.plot(np.ones(2) * (center_z + H/4 + f) / wavelength, np.array([0, 0.9]), 'k--', alpha=0.5 )
plt.xlabel('position along z axis ($\lambda_0$)')
plt.ylabel('field intensity')
plt.title('$|E(z)|^2$ along axis of metalens')
plt.show()

We now can inspect the field patterns on the area monitors using the Tidy3D built in field visualization methods. For more details see the documentation of viz_field_2D
[11]:
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
for ax, monitor, comp in zip(axes, area_monitors, ('x', 'y', 'x')):
im = sim.viz_field_2D(monitor, eps_alpha=0.99, comp=comp, val='abs', cbar=True, ax=ax)

Or you can use your own plotting functions using the raw data, similar to how we dealt with the line monitor.
[12]:
fig, axes = plt.subplots(1, len(area_monitors), figsize=(14, 4))
for ax, monitor in zip(axes, area_monitors):
data = sim.data(monitor)
E = data['E']
I = np.squeeze(np.sum(np.square(np.abs(E)), axis=0))
I = np.flipud(I.T)
im = ax.imshow(I / np.max(np.abs(I)), cmap='BuPu')
plt.colorbar(im, ax=ax)
ax.set_title(monitor.name)
plt.show()

[ ]: