Near to far field transformation

See on github, run on colab, or just follow along with the output below.

This tutorial will show you how to solve for electromagnetic fields far away from your structure using field information stored on a nearby surface.

This technique is called a ‘near field to far field transformation’ and is very useful for reducing the simulation size needed for structures involving lots of empty space.

As an example, we will simulate a simple zone plate lens with a very thin domain size to get the transmitted fields measured just above the structure. Then, we’ll show how to use the Near2Far feature from tidy3D to extrapolate to the fields at the focal plane above the lens.

[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
import matplotlib.pyplot as plt
import sys

# import client side tidy3d
import tidy3d as td
from tidy3d import web

Problem Setup

Below is a rough sketch of the setup of a near field to far field transformation.

The transmitted near fields are measured just above the metalens on the blue line, and the near field to far field transformation is then used to project the fields to the focal plane above at the red line.

c7b4fd792be746b6b42d6069b54628cc

Define Simulation Parameters

As always, we first need to define our simulation parameters. As a reminder, all length units in tidy3D are specified in microns.

[2]:
# 1 nanometer in units of microns (for conversion)
nm = 1e-3

# free space central wavelength
wavelength = 1.0

# numerical aperture
NA = 0.8

# thickness of lens features
H = 200 * nm

# space between bottom PML and substrate (-z)
# and the space between lens structure and top pml (+z)
space_below_sub = 1.5 * wavelength

# thickness of substrate (um)
thickness_sub = wavelength / 2

# side length (xy plane) of entire metalens (um)
length_xy = 40 * wavelength

# Lens and substrate refractive index
n_TiO2 = 2.40
n_SiO2 = 1.46

# define material properties
air = td.Medium(epsilon=1.0)
SiO2 = td.Medium(epsilon=n_SiO2**2)
TiO2 = td.Medium(epsilon=n_TiO2**2)

# resolution of simulation (15 or more grids per wavelength is adequate)
grids_per_wavelength = 20

# Number of PML layers to use around edges of simulation, choose thickness of one wavelength to be safe
npml = grids_per_wavelength

Process Geometry

Next we perform some conversions based on these parameters to define the simulation.

[3]:
# grid size (um)
dl = wavelength / grids_per_wavelength

# because the wavelength is in microns, use builtin td.C_0 (um/s) to get frequency in Hz
f0 = td.C_0 / wavelength

# Define PML layers, for this application we surround the whole structure in PML to isolate the fields
pml_layers = [npml, npml, npml]

# domain size in z, note, we're just simulating a thin slice: (space -> substrate -> lens thickness -> space)
length_z = space_below_sub + thickness_sub + H + space_below_sub

# construct simulation size array
sim_size = np.array([length_xy, length_xy, length_z])

Create Geometry

Now we create the ring metalens programatically

[4]:
# define substrate
substrate = td.Box(
    center=[0, 0, -length_z/2 + space_below_sub + thickness_sub / 2.0],
    size=[td.inf, td.inf, thickness_sub],
    material=SiO2)

# create a running list of structures
geometry = [substrate]

# focal length
focal_length = length_xy / 2 / NA * np.sqrt(1 - NA**2)

# location from center for edge of the n-th inner ring, see https://en.wikipedia.org/wiki/Zone_plate
def edge(n):
    return np.sqrt(n * wavelength * focal_length + n**2 * wavelength**2 / 4)

# loop through the ring indeces until it's too big and add each to geometry list
n = 1
r = edge(n)

while r < 2 * length_xy:
    # progressively wider cylinders, material alternating between air and TiO2
    cyl = td.Cylinder(
        center = [0,0,-length_z/2  + space_below_sub + thickness_sub + H / 2],
        axis='z',
        radius=r,
        height=H,
        material=TiO2 if n % 2 == 0 else air,
        name=f'cylinder_n={n}'
    )
    geometry.append(cyl)
    n += 1
    r = edge(n)

# reverse geometry list so that inner, smaller rings are added last and therefore override larger rings.
geometry.reverse()

Create Source

Create a plane wave incident from below the metalens

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

# Gaussian source offset; the source peak is at time t = offset/fwidth
offset = 4.

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

source = td.PlaneWave(
            source_time=gaussian,
            injection_axis='+z',
            position=-length_z/2 + space_below_sub / 2, # halfway between PML and substrate
            polarization='x')

# Simulation run time
run_time = 40 / fwidth

Create Monitor

Create a near field monitor to measure the fields just above the metalens

[6]:
# place it halfway between top of lens and PML
monitor_near = td.FreqMonitor(
    center=[0., 0., -length_z/2 + space_below_sub + thickness_sub + H + space_below_sub / 2],
    size=[length_xy, length_xy, 0],
    freqs=[f0],
    name='near_field')

Create Simulation

Put everything together and define a simulation object

[7]:
sim = td.Simulation(size=sim_size,
                    mesh_step=[dl, dl, dl],
                    structures=geometry,
                    sources=[source],
                    monitors=[monitor_near],
                    run_time=run_time,
                    pml_layers=pml_layers)
Initializing simulation...
Mesh step (micron): [5.00e-02, 5.00e-02, 5.00e-02].
Simulation domain in number of grid points: [840, 840, 114].
Total number of computational grid points: 8.04e+07.
Total number of time steps: 15397.
Estimated data size (GB) of monitor near_field: 0.0307.

Visualize Geometry

Lets take a look and make sure everything is defined properly

[8]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(24, 8))

# Time the visualization of the 2D plane
sim.viz_eps_2D(normal='x', position=0.1, ax=ax1);
sim.viz_eps_2D(normal='z', position=-length_z/2 + space_below_sub + thickness_sub + H / 2, ax=ax2);
[8]:
<AxesSubplot:title={'center':'xy-plane at z=2.50e-01um'}, xlabel='x (um)', ylabel='y (um)'>
../_images/examples_Near2Far_ZonePlate_16_1.png

Run Simulation

Now we can run the simulation and download the results

[9]:
# Run simulation
project = web.new_project(sim.export(), task_name='near2far_docs')
web.monitor_project(project['taskId'])
Uploading the json file...
Project 'near2far_docs-16b6e151d8baf288' status: success...

[10]:
# 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 f:
     print(f.read())
Downloading results
Applying source normalization to all frequency monitors using source index 0.
Simulation domain Nx, Ny, Nz: [840, 840, 114]
Applied symmetries: [0, 0, 0]
Number of computational grid points: 8.0438e+07.
Using subpixel averaging: True
Number of time steps: 15397
Automatic shutoff factor: 1.00e-05
Time step (s): 8.6662e-17

Compute source modes time (s):     0.0573
Compute monitor modes time (s):    0.0699

Rest of setup time (s):            10.1887

Starting solver...
- Time step    245 / time 2.12e-14s (  1 % done), field decay: 1.00e+00
- Time step    615 / time 5.33e-14s (  4 % done), field decay: 6.31e-02
- Time step   1231 / time 1.07e-13s (  8 % done), field decay: 8.24e-03
- Time step   1847 / time 1.60e-13s ( 12 % done), field decay: 2.83e-03
- Time step   2463 / time 2.13e-13s ( 16 % done), field decay: 5.84e-04
- Time step   3079 / time 2.67e-13s ( 20 % done), field decay: 1.30e-04
- Time step   3695 / time 3.20e-13s ( 24 % done), field decay: 4.33e-05
- Time step   4311 / time 3.74e-13s ( 28 % done), field decay: 2.61e-05
- Time step   4927 / time 4.27e-13s ( 32 % done), field decay: 1.64e-05
- Time step   5542 / time 4.80e-13s ( 36 % done), field decay: 1.39e-05
- Time step   6158 / time 5.34e-13s ( 40 % done), field decay: 1.20e-05
- Time step   6774 / time 5.87e-13s ( 44 % done), field decay: 1.13e-05
- Time step   7390 / time 6.40e-13s ( 48 % done), field decay: 1.09e-05
- Time step   8006 / time 6.94e-13s ( 52 % done), field decay: 1.04e-05
- Time step   8622 / time 7.47e-13s ( 56 % done), field decay: 9.66e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   346.5106
Post-processing time (s):          0.3770

Visualization

Let’s inspect the near field using the Tidy3D builtin field visualization methods.
For more details see the documentation of viz_field_2D.
[11]:
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, val in zip(axes, ('re', 'abs', 'int')):
    im = sim.viz_field_2D(monitor_near, eps_alpha=0, comp='x', val=val, cbar=True, ax=ax)
plt.show()
../_images/examples_Near2Far_ZonePlate_21_0.png

Setting Up Near 2 Far

To set up near to far, we first need to grab the data from the nearfield monitor.

[12]:
# near field monitor data dictionary
monitor_data = sim.data(monitor_near)

# grab the raw data for plotting later
xs = monitor_data['xmesh']
ys =monitor_data['ymesh']
E_near = np.squeeze(monitor_data['E'])

Then, we create a td.Near2Far object using the monitor data dictionary as follows.

This object just stores near field data and provides various methods for looking at various far field quantities.

[13]:
# from near2far_tidy3d import Near2Far
n2f = td.Near2Far(monitor_data)

Getting Far Field Data

With the Near2Far object initialized, we just need to call one of it’s methods to get a far field quantity.

For this example, we use Near2Far.get_fields_cartesian(x,y,z) to get the fields at an x,y,z point relative to the monitor center.

Below, we scan through x and y points in a plane located at z=z0 and record the far fields.

[14]:
# points to project to
num_far = 40
xs_far = 4 * wavelength * np.linspace(-0.5, 0.5, num_far)
ys_far = 4 * wavelength * np.linspace(-0.5, 0.5, num_far)

# get a mesh in cartesian, convert to spherical
Nx, Ny = len(xs), len(ys)

# initialize the far field values
E_far = np.zeros((3, num_far, num_far), dtype=complex)
H_far = np.zeros((3, num_far, num_far), dtype=complex)

# loop through points in the output plane
for i in range(num_far):
    sys.stdout.write(" \rGetting far fields, %2d%% done"%(100*i/(num_far + 1)))
    sys.stdout.flush()
    x = xs_far[i]
    for j in range(num_far):
        y = ys_far[j]

        # compute and store the outputs from projection function at the focal plane
        E, H = n2f.get_fields_cartesian(x, y, focal_length)
        E_far[:, i, j] = E
        H_far[:, i, j] = H
sys.stdout.write("\nDone!")
Getting far fields, 95% done
Done!

Plot Results

Now we can plot the near and far fields together

[15]:
# plot everything
f, ((ax1, ax2, ax3),
    (ax4, ax5, ax6)) = plt.subplots(2, 3, tight_layout=True, figsize=(10, 5))

def pmesh(xs, ys, array, ax, cmap):
    im = ax.pcolormesh(xs, ys, array.T, cmap=cmap, shading='auto')
    return im

im1 = pmesh(xs, ys, np.real(E_near[0]), ax=ax1, cmap='RdBu')
im2 = pmesh(xs, ys, np.real(E_near[1]), ax=ax2, cmap='RdBu')
im3 = pmesh(xs, ys, np.real(E_near[2]), ax=ax3, cmap='RdBu')
im4 = pmesh(xs_far, ys_far, np.real(E_far[0]), ax=ax4,  cmap='RdBu')
im5 = pmesh(xs_far, ys_far, np.real(E_far[1]), ax=ax5, cmap='RdBu')
im6 = pmesh(xs_far, ys_far, np.real(E_far[2]), ax=ax6, cmap='RdBu')


ax1.set_title('near field $E_x(x,y)$')
ax2.set_title('near field $E_y(x,y)$')
ax3.set_title('near field $E_z(x,y)$')
ax4.set_title('far field $E_x(x,y)$')
ax5.set_title('far field $E_y(x,y)$')
ax6.set_title('far field $E_z(x,y)$')

plt.colorbar(im1, ax=ax1)
plt.colorbar(im2, ax=ax2)
plt.colorbar(im3, ax=ax3)
plt.colorbar(im4, ax=ax4)
plt.colorbar(im5, ax=ax5)
plt.colorbar(im6, ax=ax6)

plt.show()
../_images/examples_Near2Far_ZonePlate_29_0.png
[16]:
# we can also use the far field data and plot the field intensity to see the focusing effect

intensity_far = np.sum(np.square(np.abs(E_far)), axis=0)

_, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))

im1 = pmesh(xs_far, ys_far, intensity_far, ax=ax1, cmap='magma')
im2 = pmesh(xs_far, ys_far, np.sqrt(intensity_far), ax=ax2, cmap='magma')

ax1.set_title('$|E(x,y)|^2$')
ax2.set_title('$|E(x,y)|$')

plt.colorbar(im1, ax=ax1)
plt.colorbar(im2, ax=ax2)
plt.show()
../_images/examples_Near2Far_ZonePlate_30_0.png