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.

# 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.

# 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.

# 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

# 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

# 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.

# 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

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

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);

<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

# 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

# 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

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.

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()