Vortex beam generation with metasurfaces leverages ultra-thin, planar structures to induce helical phase fronts in light, forming beams with orbital angular momentum (OAM). These metasurfaces are composed of precisely designed arrays of nanostructures that manipulate light at subwavelength scales, enabling the creation of vortex beams with high efficiency and accuracy. By engineering the phase profile of the metasurface, vortex beams with specific topological charges can be generated, providing compact and versatile solutions for applications in optical communication, quantum information processing, and particle manipulation.

In this notebook, we will use Tidy3D to simulate the plasmonic interface fabricated as described in the paper titled: Yu, N., Genevet, P., Kats, M. A., Aieta, F., Tetienne, J.-P., Capasso, F., & Gaburro, Z. (n.d.). Light Propagation with Phase Discontinuities: Generalized Laws of Reflection and Refraction Doi: 10.1126/science.1210713.

The metasurface is composed of 8 different unit cells, each distributed in one of the 8 octants of the metasurface. This arrangement creates a phase discontinuity in steps of $\pi/4$, inducing a spiral-like phase shift in an incident wave and generating a vortex beam.

For more examples involving metamaterials and metasurfaces, please refer to our example library, were you can find interesting case studies such as Gradient metasurface reflector, Dielectric metasurface absorber, Tunable chiral metasurface based on phase change material and a inverse-design-optimized diffractive metasurface.

If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.

FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.

Simulation Setup¶

First, we will import the necessary libraries.

[1]:

# importing libraries

from warnings import warn

import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
from tidy3d import web

The simulation will be parameterized by the number of unit cells in each direction (variable N), the distance of the source from the structure and from the PML, and the lattice constant.

The model consists of a semi-infinite silicon substrate in the -z direction, the metasurface, and an air region. An x-polarized Gaussian beam is incident from the -z to the +z direction. A field monitor is placed to store the field profile in the frequency domain and project it to the far-field.

Simulation Parameters¶

It is always good practice in FDTD simulations to parameterize the setup. This way, changing the structure only requires modifying a single variable.

[2]:

N = 20

# even distribution is used to avoid octant ambiguities
if N % 2 != 0:
    N += 1
    warn("Odd numbers of unit cells are not allowed. Adding +1 to the NN variable.")

# target wavelength
Wl = 8

# lattice constant
lattice_constant = 1.5

# defining the central frequency and linewidth of the source
fcen = td.C_0 / Wl
fwidth = 0.2 * fcen

# simulation run time
run_time = 1e-12


# parameters of the fabricated antennas
width = 0.22
thickness = 0.05

# distance between structures and pml along z
pml_buffer = 4

Now, we will define each parameter needed to build the simulation as a function of the parameters defined above:

[3]:

# the size of the metasurface plane as a function of the number of unit cells and the lattice constant
#   a 1.5 lattice constant is added on each side to avoid structures touching the PML

sx = (N + 3) * lattice_constant
sy = sx

# defining the z-size of the simulation
sz = 2 * pml_buffer + thickness

# defining the Z position of the metasurface
structure_z_position = 0

# defining monitor Z position
monitor_z_position = sz / 2 - pml_buffer / 2

# defining the center position and size of the source
center_source = (0, 0, -sz / 2 + pml_buffer / 2)
size_source = (td.inf, td.inf, 0)

Here, we will define an auxiliary function to return the antenna structure. The antenna is modeled as two Box arms that are tilted apart from each other by the angle delta. The antennas are then rotated by an angle theta around an axis that passes through their center, defined as the junction of the arms. To make the structures smoother and more similar to the SEM images, a cylinder is added at the arms' junction.

The vAntennaBlock function returns a Transformed object. For more information on how to apply geometry transformations, refer to this notebook.

[4]:

def vAntennaBlock(center, delta, theta, length):
    """
    Create a V-antenna structure.

    Parameters:
        center (tuple): Center coordinates of the antenna.
        delta (float): Angle in between the antenna arms (deg).
        theta (float): Tilt angle of the antenna arms (deg).
        length (float): Length of the antenna arms.

    Returns:
        Tidy3D structure representing the V-antenna.
    """

    # Create arms and cylinder
    arm1 = td.Box(size=(length, width, thickness), center=(0, 0, 0))
    arm2 = arm1.copy()
    cylinder1 = td.Cylinder(radius=width / 2, center=(0, 0, 0), axis=2, length=thickness)

    # rotate each arm around its left corner in opposite directions
    arm1 = arm1.translated(length / 2, 0, 0).rotated((delta / 2) * np.pi / 180, 2)
    arm2 = arm2.translated(length / 2, 0, 0).rotated((-delta / 2) * np.pi / 180, 2)

    antenna = cylinder1 + arm1 + arm2

    # rotate and translate the whole antenna
    antenna = antenna.rotated(theta * np.pi / 180, 2).translated(*center)

    return antenna


# create a dictionary to define the parameters of the anthenas in each octante
# to create the appropriate phase change
dicNumerator = {
    0: (90, -45, 1.1),
    1: (90 + 45, -45, 0.9),
    2: (180, -45, 0.8),
    3: (45, 45, 1.3),
    4: (-90, 45, 1.1),
    5: (-90 - 45, 45, 0.9),
    6: (180, 45, 0.8),
    7: (45, -45, 1.3),
}

[5]:

# plot each antenna
fig, ax = plt.subplots(1, 8)
for i in dicNumerator.keys():
    vAntennaBlock((0, 0, 0), *dicNumerator[i]).plot(
        z=0,
        ax=ax[i],
        linewidth=0,
        facecolor="navy",  # patches kwargs
    )
    ax[i].set_axis_off()
    ax[i].set_title("")
fig.tight_layout()
plt.show()

No description has been provided for this image

Building the Metasurface¶

To properly position the unit cells and build the metasurface, we defined the function classifyOctant, which takes the x and y position of a unit cell and returns the octant where the structure belongs. The function is adapted so that the first octant matches the zero-phase antenna.

[6]:

def classifyOctant(x, y):
    angle = 180 + 180 * (np.arctan2(y, x) / np.pi)
    octant = angle // 45.00
    return octant, angle

[7]:

metasurface = None
# loops iterating in both dimensions
for x in range(N):
    for y in range(N):
        # defining the position of the unit cell
        posX = ((N / 2) - x) * lattice_constant - lattice_constant * 0.5
        posY = ((N / 2) - y) * lattice_constant - lattice_constant * 0.5

        # defining angle and octant
        octant, angle = classifyOctant(posX, posY)

        # condition to remove the 45 degree row in the lower half of the metasurface, for it is more similar
        # to what was done in the experiment
        if (abs(angle - 45) >= 0.1) and not (posX == 0 and posY > 0):
            # antenna parameters for a given octant
            armArgs = dicNumerator[octant]

            # adding to the structure
            if metasurface is None:
                metasurface = vAntennaBlock((posX, posY, structure_z_position), *armArgs)
            else:
                metasurface += vAntennaBlock((posX, posY, structure_z_position), *armArgs)

# define the medium e create the structure object
medium = td.material_library["Au"]["Olmon2012crystal"]  # built-in gold model for the
# frequency range of interest
metasurface = td.Structure(geometry=metasurface, medium=medium)

metasurface.plot(z=structure_z_position, linewidth=0)
plt.show()

Creating Simulation Object¶

Since the structures are much smaller than the target wavelength, it is necessary to create a mesh override region around the metasurface to properly resolve the structures, while using a coarser mesh for free space propagation.

We will also define a Gaussain beam source source with a gaussian time profile.

Additionally, we will include the substrate region.

For the monitors, we will add a far-field monitor, to project the fields on the server side, and a field monitor that will store the near-field information for projection locally after the simulation.

Note that we set colocate = False for the FieldMonitor, which means the raw field data will be recorded on the Yee grid. This option avoids double interpolation, leading to more accurate results for far-field projection on the user side.

[8]:

mesh_override = td.MeshOverrideStructure(
    geometry=td.Box(center=(0, 0, structure_z_position), size=(sx, sy, thickness)),
    dl=(0.02,) * 3,
)

source_time = td.GaussianPulse(freq0=fcen, fwidth=fwidth)

source = td.GaussianBeam(
    center=center_source,
    size=size_source,
    direction="+",
    pol_angle=0,
    source_time=source_time,
    waist_radius=16,
    waist_distance=center_source[2],
)

substrate = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-2 * sx, -2 * sy, -100),
        rmax=(2 * sx, 2 * sy, -thickness / 2),
    ),
    medium=td.Medium(permittivity=3.47**2),
)

far_field_monitor = td.FieldProjectionCartesianMonitor(
    center=(0, 0, monitor_z_position),
    size=(td.inf, td.inf, 0),
    normal_dir="+",
    freqs=[fcen],
    x=300 * np.linspace(-0.5, 0.5, 100),
    y=300 * np.linspace(-0.5, 0.5, 100),
    proj_axis=2,
    proj_distance=20 * Wl,
    far_field_approx=False,
    name="farFieldMon",
)

field_monitor = td.FieldMonitor(
    center=far_field_monitor.center,
    size=far_field_monitor.size,
    freqs=far_field_monitor.freqs,
    name="fieldMon",
    colocate=False,
)

# creating the simulation object

sim = td.Simulation(
    size=(sx, sy, sz),
    grid_spec=td.GridSpec.auto(
        min_steps_per_wvl=20,  # grid specification and mesh override
        override_structures=[mesh_override],
    ),
    structures=[substrate] + [metasurface],
    sources=[source],
    monitors=[far_field_monitor, field_monitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec(  # all boundaries are pmls in order to simulate a finite structure
        x=td.Boundary.pml(), y=td.Boundary.pml(), z=td.Boundary.pml()
    ),
)

22:46:37 -03 WARNING: Structure at structures[1] was detected as being less than
             half of a central wavelength from a PML on side x-min. To avoid    
             inaccurate results or divergence, please increase gap between any  
             structures and PML or fully extend structure through the pml.

             WARNING: Suppressed 3 WARNING messages.

We notice a warning that the structures are positioned close to the PML. However, since we're using a Gaussian beam source focused on the center of the structures, there will be no evanescent fields near the PML. Therefore, we can disregard this warning.

To make sure that everything is correct, we will plot the y-z and x-y planes of the model.The green line represents the source, the yellow line represents the monitors, the blue area indicates the substrate, and the red area shows the cross-section of the antennas. the override mesh region is depicted in gray.

[9]:

fig, Ax = plt.subplots(1, 2, tight_layout=True, figsize=(8, 3.8))

sim.plot(x=1, ax=Ax[0])
mesh_override.plot(x=1, ax=Ax[0], alpha=0.2)
Ax[0].set_xlim(-5, 5)
Ax[0].set_ylim(-5, 5)

sim.plot(z=structure_z_position, ax=Ax[1])
plt.show()

[10]:

# estimating the credit cost of the simulation
task_id = web.upload(sim, task_name="vortexMetasurface")
web.estimate_cost(task_id)

22:46:40 -03 Created task 'vortexMetasurface' with task_id                      
             'fdve-b8071be0-3346-4d48-bedb-679569aef286' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-b8071be0-334
             6-4d48-bedb-679569aef286'.

Output()

22:46:45 -03 Maximum FlexCredit cost: 2.465. Minimum cost depends on task       
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

[10]:

2.464539667330026

[11]:

# running the simulation
results = web.run(
    simulation=sim,
    task_name="vortexMetasurface",
    folder_name="data",
    path="data/1.hdf5",
    verbose="True",
)

             Created task 'vortexMetasurface' with task_id                      
             'fdve-a52c8e27-8112-4fb4-92ae-b7508879851d' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-a52c8e27-811
             2-4fb4-92ae-b7508879851d'.

Output()

22:46:50 -03 status = queued

             To cancel the simulation, use 'web.abort(task_id)' or              
             'web.delete(task_id)' or abort/delete the task in the web UI.      
             Terminating the Python script will not stop the job running on the 
             cloud.

Output()

22:46:56 -03 status = preprocess

22:47:09 -03 Maximum FlexCredit cost: 2.465. Use 'web.real_cost(task_id)' to get
             the billed FlexCredit cost after a simulation run.

             starting up solver

             running solver

Output()

22:48:17 -03 early shutoff detected at 64%, exiting.

             status = postprocess

Output()

22:48:26 -03 status = success

             View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-a52c8e27-811
             2-4fb4-92ae-b7508879851d'.

Output()

22:48:46 -03 loading simulation from data/1.hdf5

22:48:47 -03 WARNING: Structure at structures[1] was detected as being less than
             half of a central wavelength from a PML on side x-min. To avoid    
             inaccurate results or divergence, please increase gap between any  
             structures and PML or fully extend structure through the pml.

             WARNING: Suppressed 3 WARNING messages.

             WARNING: Warning messages were found in the solver log. For more   
             information, check 'SimulationData.log' or use                     
             'web.download_log(task_id)'.

Creating Far-Field Projection¶

Now we plot the far-field projections: one done in the cloud using the FarFieldMonitor and the other done locally. It is interesting to note the difference between them, as the projection done in the cloud has the window_size argument set to the default (0,0). Without applying the window filter, high-frequency artifacts appear in the projected field as they need to be truncated. On the other hand, applying the spatial filter leads to a more smooth far-field profile, more in agreement with the experimental results. More information about field projection can be found here. Since we are using a Gaussian beam source, the differences between the projections are minimal.

[12]:

# radial distance away from the origin at which to project fields
r_proj = 60 * Wl

# here we add the 'window_size' argument to ensure that the near fields decay to zero at the monitor boundaries
far_field_local_projection = td.FieldProjectionCartesianMonitor(
    center=(0, 0, monitor_z_position),
    size=(td.inf, td.inf, 0),
    normal_dir="+",
    freqs=[fcen],
    x=1000 * np.linspace(-0.5, 0.5, 100),
    y=1000 * np.linspace(-0.5, 0.5, 100),
    proj_axis=2,
    proj_distance=r_proj,
    far_field_approx=False,
    window_size=(0.5, 0.5),
    name="farFieldLocalProjection",
)

projector = td.FieldProjector.from_near_field_monitors(
    sim_data=results,
    near_monitors=[field_monitor],
    normal_dirs=["+"],
    pts_per_wavelength=10,
)
projected = projector.project_fields(far_field_local_projection)
far_field = projected.fields_cartesian

Output()

[13]:

# also plot the projected fields from the FarFieldMonitor
fig, Ax = plt.subplots(1, 2, tight_layout=True, figsize=(8, 3.8))

abs(results["farFieldMon"].fields_cartesian["Ey"]).plot(cmap="inferno", ax=Ax[0])
Ax[0].set_title("window_size = (0,0)")


# plotting only the cross-polarized component (Ey). In the experiments, a polarizer was used to block the
# scattered light with the same polarization as the incident wave.

# visualizing the electric field
abs(far_field["Ey"]).plot(cmap="inferno", ax=Ax[1])
Ax[1].set_title("window_size = (0.5,0.5)")
plt.show()

One can note that the results are consistent with those reported in the original paper, specifically in Figure 5(c)

Both the simulation and the experiment reveal the formation of a phase singularity at the center of the beam, leading to a region of zero field.

To obtain the phase profile, you can interfere the beam with a collinear Gaussian beam. To reproduce this, use the analytical expression for the Gaussian beam to create an interference pattern with the simulated beam.

[14]:

# gaussiam beam function
def E(r, z, lda=1, E0=1, w0=1, n=1):
    k = 2 * np.pi / lda
    zr = np.pi * w0**2 * n / lda
    wz = w0 * np.sqrt(1 + (z / zr) ** 2)
    Rz = z * (1 + (zr / z) ** 2)
    phi = np.arctan(z / zr)

    return (
        E0
        * (w0 / wz)
        * np.exp(-(r**2) / (wz**2))
        * np.exp(-1j * (k * z + k * r**2 / (2 * Rz) - phi))
    )


# creating the data points
x = np.linspace(-40, 40, 100)
y = x.copy()

X, Y = np.meshgrid(x, y)
r = np.sqrt(X**2 + Y**2)

# creating the gaussiam beam field
gb = E(r, 1, 8, w0=20)

# plotting the interference pattern
fig, ax = plt.subplots()
ax.imshow(
    abs(
        results["farFieldMon"].fields_cartesian["Ey"].squeeze().T
        / results["farFieldMon"].fields_cartesian["Ey"].max()
        + gb
    ),
    interpolation="bilinear",
    cmap="inferno",
)
ax.axis("off")
plt.show()