Large-area metalens#

Run this notebook in your browser using Binder. View this project in Tidy3D Web App.

Introduction#

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!

First, we’ll show a simple approach for computing fields on the focal plane by simply placing a FieldMonitor at the focal plane to record the fields. However, the focal plane is several wavelengths away from the structure, which means that the above setup will involve simulating fields propagating through lots of empty space. Therefore, we will also show how one can instead record the near fields on a large surface just above the structure, and then project them to the focal plane without actually simulating the fields in the empty region in between, leading to a further improvement in the already-impressive simulation speed.

Setup#

We first perform basic imports of the packages needed.

[1]:
# standard python imports
import numpy as np
from numpy import random
import matplotlib.pyplot as plt

import tidy3d as td
from tidy3d import web
[18:50:12] WARNING  This version of Tidy3D was pip installed from the 'tidy3d-beta' repository on   __init__.py:102
                    PyPI. Future releases will be uploaded to the 'tidy3d' repository. From now on,                
                    please use 'pip install tidy3d' instead.                                                       
           INFO     Using client version: 1.8.2                                                     __init__.py:120

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.

diagram

[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(permittivity=1.0)
SiO2 = td.Medium(permittivity=n_SiO2**2)
TiO2 = td.Medium(permittivity=n_TiO2**2)
for diameter of 10.0 um, have 23 cells per side
full metalens has area of 100.0 um^2 and 529 total cells
[3]:
# 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

# 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 = (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.Structure(
    geometry=td.Box(
        center=[0, 0, -length_z/2 + space_below_sub + thickness_sub / 2.0],
        size=[td.inf, td.inf, thickness_sub],
    ),
    medium=SiO2
)

# 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

# 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]])


xs, ys = np.meshgrid(centers_x, centers_y, indexing='ij')
xs = xs.flatten()
ys = ys.flatten()

angles = theta(xs, ys)

# 2x2 rotation matrix angle `angle` with respect to x axis
rotation_matrix = np.array([
    [+np.cos(angles), -np.sin(angles)],
    [+np.sin(angles), +np.cos(angles)]
])

# rotate the origin vertices by this angle
vertices_rotated = np.einsum('ij, jkn -> nik', vertices_origin, rotation_matrix)

# shift the rotated vertices to be centered at (xs, ys)
vertices_shifted = vertices_rotated + np.stack([xs, ys], axis=-1)[:, None, :]

metalens_geometry = []
for vertices in vertices_shifted:
    # create a tidy3D PolySlab with these rotated and shifted vertices and thickness `H`
    metalens_geometry.append(
        td.PolySlab(
            vertices=vertices.tolist(),
            slab_bounds=(center_z-H/2, center_z+H/2),
            axis=2,
        ),
    )

metalens = td.Structure(
    geometry=td.GeometryGroup(geometries=metalens_geometry),
    medium=TiO2
)

Define grid specification#

We define the grid based on the properties of the geometry. The metalens is quasi-periodic in x and y, in that we have clearly defined unit cells, but each is slightly modified from its neighbors. Such structures are best resolved with a grid that matches the periodicity, which is why we use a uniform grid in x and y. In z, we use the automatic nonuniform grid that will place a higher grid density around the metalens region, and a lower one in the air region away from the metalens. To speed up the auto meshing in the region with the pillars, we put an override box in the grid specification.

[5]:
# steps per unit cell along x and y
grids_per_unit_length = 20
# uniform mesh in x and y
grid_x = td.UniformGrid(dl=S / grids_per_unit_length)
grid_y = td.UniformGrid(dl=S / grids_per_unit_length)
# in z, use an automatic nonuniform mesh with the wavelength being the "unit length"
grid_z = td.AutoGrid(min_steps_per_wvl=grids_per_unit_length)
# we need to supply the wavelength because of the automatic mesh in z
grid_spec = td.GridSpec(wavelength=wavelength, grid_x=grid_x, grid_y=grid_y, grid_z=grid_z)
# put an override box over the pillars to avoid parsing a large amount of structures in the mesher
grid_spec = grid_spec.copy(update=dict(override_structures = [
    td.Structure(
        geometry=td.Box.from_bounds(
            rmin=(-td.inf, -td.inf, -length_z/2 + space_below_sub),
            rmax=(td.inf, td.inf, center_z + H/2)
        ),
        medium=TiO2,
    )
]))

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

[6]:
# Bandwidth in Hz
fwidth = f0 / 10.0

# time dependence of source
gaussian = td.GaussianPulse(freq0=f0, fwidth=fwidth, phase=0)

source = td.PlaneWave(
    source_time=gaussian,
    size=(td.inf, td.inf, 0),
    center=(0,0,-length_z/2 + space_below_sub / 10.0),
    direction='+',
    pol_angle=0)

run_time = 50 / 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 field monitor documentation.

[7]:
# To decrease the amount of data stored, only store the E field in 2D monitors
fields = ["Ex", "Ey", "Ez"]

# get fields along x=y=0 axis
monitor_center = td.FieldMonitor(
    center=[0., 0., 0],
    size=[0, 0, td.inf],
    freqs=[f0],
    name='center'
)

# get the fields at a few cross-sectional planes
monitor_xz = td.FieldMonitor(
    center=[0., 0., 0.],
    size=[td.inf, 0., td.inf],
    freqs=[f0],
    name='xz',
    fields=fields,
)

monitor_yz = td.FieldMonitor(
    center=[0., 0., 0.],
    size=[0., td.inf, td.inf],
    freqs=[f0],
    name='yz',
    fields=fields,
)

monitor_xy = td.FieldMonitor(
    center=[0., 0., center_z + H/2 + f],
    size=[td.inf, td.inf, 0],
    freqs=[f0],
    name='focal_plane',
    fields=fields,
)

# 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 object to be run.

We get a number of warnings about structures being too close to the PML. In FDTD simulations, this can result in instability, as PML are absorbing for propagating fields, but can be amplifying for evanescent fields. This particular simulation runs without any issues even with PML on the sides, but it is best to heed these warnings to avoid problems. There are two ways that we can fix the simulation: one is to just put some space between the last of the metalens boxes and the PML. The other is to use adiabatic absorbers on the sides, which are always stable. The only downside of the absorbers is that they are slightly thicker than the PML, making the overall simulation size slighlty larger. This is why we only put them along x and y, while we leave the PML in z.

[8]:
sim = td.Simulation(
    size=sim_size,
    grid_spec=grid_spec,
    structures=[substrate, metalens],
    sources=[source],
    monitors=monitors,
    run_time=run_time,
    boundary_spec=td.BoundarySpec(
        x=td.Boundary.absorber(),
        y=td.Boundary.absorber(),
        z=td.Boundary.pml()
    )
)

Visualize Geometry#

Lets take a look and make sure everything is defined properly

[9]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 6))

sim.plot(x=0.1, ax=ax1);
sim.plot(y=0.1, ax=ax2);
sim.plot(z=-length_z/2 + space_below_sub + thickness_sub + H / 2, ax=ax3);
           WARNING  Override structures take no effect along x-axis. If intending to apply         grid_spec.py:519
                    override structures to this axis, use 'AutoGrid'.                                              
           WARNING  Override structures take no effect along y-axis. If intending to apply         grid_spec.py:519
                    override structures to this axis, use 'AutoGrid'.                                              
[9]:
<AxesSubplot: title={'center': 'cross section at z=-6.93'}, xlabel='x', ylabel='y'>
../_images/notebooks_Metalens_18_3.png

Run Simulation#

Now we can run the simulation over time and measure the results to plot

[10]:
job = web.Job(simulation=sim, task_name='metalens')
sim_data = job.run(path='data/simulation_data.hdf5')
[18:50:13] INFO     Using Tidy3D credentials from stored file.                                           auth.py:70
[18:50:17] INFO     Authentication successful.                                                           auth.py:30
[18:50:18] INFO     Created task 'metalens' with task_id 'b08f5061-d7f9-4913-af19-7130dc512887'.      webapi.py:120
[18:50:24] INFO     Maximum FlexUnit cost: 1.631                                                      webapi.py:252
           INFO     status = queued                                                                   webapi.py:261
[18:50:26] INFO     status = preprocess                                                               webapi.py:273
[18:50:35] INFO     starting up solver                                                                webapi.py:277
[18:50:46] INFO     running solver                                                                    webapi.py:283
[18:51:06] INFO     early shutoff detected, exiting.                                                  webapi.py:294
           INFO     status = postprocess                                                              webapi.py:300
[18:51:11] INFO     status = success                                                                  webapi.py:306
           INFO     Billed FlexUnit cost: 0.527                                                       webapi.py:310
           INFO     downloading file "output/monitor_data.hdf5" to "data/simulation_data.hdf5"        webapi.py:592
[18:51:15] INFO     loading SimulationData from data/simulation_data.hdf5                             webapi.py:414
[11]:
print(sim_data.log)
Simulation domain Nx, Ny, Nz: [541, 541, 589]
Applied symmetries: (0, 0, 0)
Number of computational grid points: 1.7426e+08.
Using subpixel averaging: True
Number of time steps: 3.1633e+04
Automatic shutoff factor: 1.00e-05
Time step (s): 3.1636e-17


Compute source modes time (s):     0.8775
Compute monitor modes time (s):    0.0056
Rest of setup time (s):            4.4001

Running solver for 31633 time steps...
- Time step    503 / time 1.59e-14s (  1 % done), field decay: 1.00e+00
- Time step   1265 / time 4.00e-14s (  4 % done), field decay: 1.00e+00
- Time step   2530 / time 8.00e-14s (  8 % done), field decay: 1.25e-01
- Time step   3795 / time 1.20e-13s ( 12 % done), field decay: 5.74e-03
- Time step   5061 / time 1.60e-13s ( 16 % done), field decay: 9.60e-04
- Time step   6326 / time 2.00e-13s ( 20 % done), field decay: 1.84e-04
- Time step   7591 / time 2.40e-13s ( 24 % done), field decay: 4.23e-05
- Time step   8857 / time 2.80e-13s ( 28 % done), field decay: 1.31e-05
- Time step  10122 / time 3.20e-13s ( 32 % done), field decay: 5.96e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   14.4930


As we can see from the logs, the total time to solve this problem (not including data transfer and pre/post-processing) was on the order of a minute!

For reference, the same problem run with FDTD on a single CPU core FDTD is projected to take 40 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

[12]:
focal_z = center_z+H/2+f
data_center_line = sim_data.at_centers('center')
I = abs(data_center_line.Ex)**2 + abs(data_center_line.Ey)**2 + abs(data_center_line.Ez)**2
I.plot()
plt.title('intensity(z)')
plt.show()
           WARNING  Override structures take no effect along x-axis. If intending to apply         grid_spec.py:519
                    override structures to this axis, use 'AutoGrid'.                                              
           WARNING  Override structures take no effect along y-axis. If intending to apply         grid_spec.py:519
                    override structures to this axis, use 'AutoGrid'.                                              
../_images/notebooks_Metalens_24_2.png

We now can inspect the field patterns on the area monitors using the Tidy3D visualization nethods.

[13]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 4))
sim_data.plot_field('xz', 'Ex', val='real', f=f0, y=0, ax=ax1)
sim_data.plot_field('yz', 'Ex', val='real', f=f0, x=0, ax=ax2)
sim_data.plot_field('focal_plane', 'Ex', val='real', f=f0, z=focal_z, ax=ax3)
ax1.set_title('x-z plane')
ax2.set_title('y-z plane')
ax3.set_title('x-y (focal) plane')
plt.show()
../_images/notebooks_Metalens_26_0.png
[14]:
# plot absolute value for good measure
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 4))
sim_data.plot_field('xz', 'Ex', val='abs', f=f0, y=0, ax=ax1)
sim_data.plot_field('yz', 'Ex', val='abs', f=f0, x=0, ax=ax2)
sim_data.plot_field('focal_plane', 'Ex', val='abs', f=f0, z=focal_z, ax=ax3)
ax1.set_title('x-z plane')
ax2.set_title('y-z plane')
ax3.set_title('x-y (focal) plane')
plt.show()
../_images/notebooks_Metalens_27_0.png
[15]:
# and let's plot the intensites as well
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 4))
sim_data.plot_field('xz', 'int', f=f0, y=0, ax=ax1)
sim_data.plot_field('yz', 'int', f=f0, x=0, ax=ax2)
sim_data.plot_field('focal_plane', 'int', f=f0, z=focal_z, ax=ax3)
ax1.set_title('x-z plane')
ax2.set_title('y-z plane')
ax3.set_title('x-y (focal) plane')
plt.show()
../_images/notebooks_Metalens_28_0.png

Field Projection to the Focal Plane#

Here we’ll show how to avoid simulating lots of empty space above the structure by using Tidy3D’s field projection functionality.

[16]:
# create a field projection monitor in Cartesian coordinates which records near fields just above the strucure,
# and projects them to points on the focal plane

# number of focal plane sampling points in the x and y directions
num_far = 200

# define the focal plane sample points at which to project fields
xs_far = np.linspace(-sim_size[0] / 2, sim_size[0] / 2, num_far)
ys_far = np.linspace(-sim_size[1] / 2, sim_size[1] / 2, num_far)

pos_monitor_z = -6
proj_distance = monitor_xy.center[2] - pos_monitor_z
monitor_proj = td.FieldProjectionCartesianMonitor(
    center=[0., 0., pos_monitor_z],
    size=[td.inf, td.inf, 0],
    freqs=[f0],
    name='focal_plane_proj',
    proj_axis=2,           # axis along which to project, in this case z
    proj_distance=proj_distance, # distance from this monitor to where fields are projected
    x=xs_far,
    y=ys_far,
    custom_origin=[0., 0., pos_monitor_z],
    far_field_approx=False, # the distance to the focal plane is comparable to the size of the structure, so
                            # turn off the far-field approximation and use an exact Green's function to
                            # project the fields
)

# create a simulation as before, but this time there's no need to include the large amount of
# empty space up to the focal plane, and include the projection monitor

# total domain size in z: (space -> substrate -> unit cell -> space)
length_z_new = space_below_sub + thickness_sub + 2 * H + space_below_sub
sim_size = (length_xy, length_xy, length_z_new)
sim_center = (0, 0, -length_z / 2 + length_z_new / 2)
sim_new = td.Simulation(
    size=sim_size,
    center=sim_center,
    grid_spec=grid_spec,
    structures=[substrate, metalens],
    sources=[source],
    monitors=[monitor_proj],
    run_time=run_time,
    boundary_spec=td.BoundarySpec(
        x=td.Boundary.absorber(),
        y=td.Boundary.absorber(),
        z=td.Boundary.pml()
    )
)

# visualize to make sure everything looks okay
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 6))
sim_new.plot(x=0.1, ax=ax1);
sim_new.plot(y=0.1, ax=ax2);
sim_new.plot(z=-length_z/2 + space_below_sub + thickness_sub + H / 2, ax=ax3);
[18:51:22] WARNING  Override structures take no effect along x-axis. If intending to apply         grid_spec.py:519
                    override structures to this axis, use 'AutoGrid'.                                              
           WARNING  Override structures take no effect along y-axis. If intending to apply         grid_spec.py:519
                    override structures to this axis, use 'AutoGrid'.                                              
[16]:
<AxesSubplot: title={'center': 'cross section at z=-6.93'}, xlabel='x', ylabel='y'>
../_images/notebooks_Metalens_30_3.png

Run#

Run the new simulation and note the simulation time recorded in the log.

[17]:
job = web.Job(simulation=sim_new, task_name='metalens')
sim_data_new = job.run(path='data/simulation_data_new.hdf5')
print(sim_data_new.log)
[18:51:23] INFO     Created task 'metalens' with task_id '6c1d97f5-20fd-4e2d-ace6-3375575fc991'.      webapi.py:120
[18:51:26] INFO     Maximum FlexUnit cost: 0.415                                                      webapi.py:252
           INFO     status = queued                                                                   webapi.py:261
[18:51:30] INFO     status = preprocess                                                               webapi.py:273
[18:51:36] INFO     starting up solver                                                                webapi.py:277
[18:51:47] INFO     running solver                                                                    webapi.py:283
[18:52:05] INFO     early shutoff detected, exiting.                                                  webapi.py:294
           INFO     status = postprocess                                                              webapi.py:300
[18:52:07] INFO     status = success                                                                  webapi.py:306
           INFO     Billed FlexUnit cost: 0.122                                                       webapi.py:310
           INFO     downloading file "output/monitor_data.hdf5" to "data/simulation_data_new.hdf5"    webapi.py:592
[18:52:09] INFO     loading SimulationData from data/simulation_data_new.hdf5                         webapi.py:414
Simulation domain Nx, Ny, Nz: [541, 541, 144]
Applied symmetries: (0, 0, 0)
Number of computational grid points: 4.3048e+07.
Using subpixel averaging: True
Number of time steps: 3.1633e+04
Automatic shutoff factor: 1.00e-05
Time step (s): 3.1636e-17


Compute source modes time (s):     1.0552
Compute monitor modes time (s):    0.0061
Rest of setup time (s):            4.2741

Running solver for 31633 time steps...
- Time step    503 / time 1.59e-14s (  1 % done), field decay: 1.00e+00
- Time step   1265 / time 4.00e-14s (  4 % done), field decay: 1.30e-01
- Time step   2530 / time 8.00e-14s (  8 % done), field decay: 1.26e-02
- Time step   3795 / time 1.20e-13s ( 12 % done), field decay: 2.13e-03
- Time step   5061 / time 1.60e-13s ( 16 % done), field decay: 4.33e-04
- Time step   6326 / time 2.00e-13s ( 20 % done), field decay: 8.90e-05
- Time step   7591 / time 2.40e-13s ( 24 % done), field decay: 2.04e-05
- Time step   8857 / time 2.80e-13s ( 28 % done), field decay: 7.28e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   17.6998
Field projection time (s):         8.1576

The solver time was reduced by a factor of about 3, and about 10x less data needs to be downloaded compared to the original simulation! Furthermore, only a few seconds (a small fraction of the total time) were spent on projecting the near fields to the focal plane.

Plot Fields#

Finally, let’s plot the fields on the focal plane and compare those recorded previously in the large simulation to those computed via field projections just now. Comparing the actual simulated fields to the projected ones, we notice an excellent match, with a several-fold reduction in compute time.

[18]:
# plot the focal plane electric field components recorded by directly placing a monitor in the first simulation,
# which we can use as a reference
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14.3, 3.6))
focal_field_data = sim_data['focal_plane']
Ex = focal_field_data.Ex.sel(f=f0, z=focal_z)
Ey = focal_field_data.Ey.sel(f=f0, z=focal_z)
Ez = focal_field_data.Ez.sel(f=f0, z=focal_z)
im1 = ax1.pcolormesh(Ex.y, Ex.x, np.real(Ex), cmap='RdBu', shading='auto')
im2 = ax2.pcolormesh(Ey.y, Ey.x, np.real(Ey), cmap='RdBu', shading='auto')
im3 = ax3.pcolormesh(Ez.y, Ez.x, np.real(Ez), cmap='RdBu', shading='auto')
fig.colorbar(im1, ax=ax1)
fig.colorbar(im2, ax=ax2)
fig.colorbar(im3, ax=ax3)

# now plot the projected fields computed via the second simulation
proj_fields = sim_data_new['focal_plane_proj'].fields_cartesian.sel(f=f0, z=monitor_proj.proj_distance)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14.3, 3.6))
im1 = ax1.pcolormesh(ys_far, xs_far, np.real(proj_fields.Ex), cmap='RdBu', shading='auto')
im2 = ax2.pcolormesh(ys_far, xs_far, np.real(proj_fields.Ey), cmap='RdBu', shading='auto')
im3 = ax3.pcolormesh(ys_far, xs_far, np.real(proj_fields.Ez), cmap='RdBu', shading='auto')
fig.colorbar(im1, ax=ax1)
fig.colorbar(im2, ax=ax2)
fig.colorbar(im3, ax=ax3)

[18]:
<matplotlib.colorbar.Colorbar at 0x12360f370>
../_images/notebooks_Metalens_35_1.png
../_images/notebooks_Metalens_35_2.png
[ ]: