Metasurface Lens Design

This example demonstrates the design and simulation of a binary-grating metasurface lens that focuses light by engineering the phase profile across the aperture.

Application: Metasurface lenses (metalenses) enable flat, compact optical elements for imaging, displays, and telecommunications.

Physics Background

Phase Engineering

A lens focuses light by imparting a position-dependent phase:

$$ \phi(x) = \frac{2\pi}{\lambda}\left(f - \sqrt{x^2 + f^2}\right) $$

where f is the focal length and x is the transverse position.

Binary Grating Unit Cell

By varying the duty cycle of a binary grating, we can control:

Transmittance — Power transmitted through each cell
Phase — Accumulated phase shift

Simulation Strategy

Phase Map Characterization

Simulate unit cells with varying duty cycles to map phase vs. duty cycle.

Phase Profile Design

Calculate required phase at each position for desired focal length.

Lens Assembly

Select appropriate duty cycle for each grating cell.

Far-Field Focusing

Simulate complete lens and compute far-field intensity distribution.

Configuration Parameters

Parameter	Value	Description
`resolution`	50	Pixels per μm
`gp`	0.3	Grating period (μm)
`gh`	1.8	Grating height (μm)
`focal_length`	200	Target focal length (μm)
`wavelength`	0.5	Center wavelength (μm)

Complete Code

"""
Metasurface Lens Simulation with OptixLog Integration

Simulates a binary-grating metasurface lens and analyzes its focusing properties.
"""

import os
import meep as mp
import matplotlib
matplotlib.use("agg")
import matplotlib.pyplot as plt
import numpy as np
from optixlog import Optixlog


# Global parameters
RESOLUTION = 50
DPML = 1.0
DSUB = 2.0
DPAD = 2.0
LCEN = 0.5  # center wavelength
FCEN = 1 / LCEN
DF = 0.2 * FCEN
FOCAL_LENGTH = 200
SPOT_LENGTH = 100
FF_RES = 10

GLASS = mp.Medium(index=1.5)
K_POINT = mp.Vector3(0, 0, 0)
PML_LAYERS = [mp.PML(thickness=DPML, direction=mp.X)]
SYMMETRIES = [mp.Mirror(mp.Y)]


def grating_unit_cell(gp, gh, gdc):
    """Simulate single grating unit cell to get phase and transmittance"""
    sx = DPML + DSUB + gh + DPAD + DPML
    sy = gp
    cell_size = mp.Vector3(sx, sy, 0)
    
    src_pt = mp.Vector3(-0.5 * sx + DPML + 0.5 * DSUB)
    mon_pt = mp.Vector3(0.5 * sx - DPML - 0.5 * DPAD)
    
    sources = [
        mp.Source(
            mp.GaussianSource(FCEN, fwidth=DF),
            component=mp.Ez,
            center=src_pt,
            size=mp.Vector3(y=sy),
        )
    ]
    
    # Empty cell for normalization
    sim = mp.Simulation(
        resolution=RESOLUTION,
        cell_size=cell_size,
        boundary_layers=PML_LAYERS,
        k_point=K_POINT,
        default_material=GLASS,
        sources=sources,
        symmetries=SYMMETRIES,
    )
    
    flux_obj = sim.add_flux(FCEN, 0, 1, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy)))
    sim.run(until_after_sources=50)
    input_flux = mp.get_fluxes(flux_obj)
    sim.reset_meep()
    
    # With grating
    geometry = [
        mp.Block(
            material=GLASS,
            size=mp.Vector3(DPML + DSUB, mp.inf, mp.inf),
            center=mp.Vector3(-0.5 * sx + 0.5 * (DPML + DSUB)),
        ),
        mp.Block(
            material=GLASS,
            size=mp.Vector3(gh, gdc * gp, mp.inf),
            center=mp.Vector3(-0.5 * sx + DPML + DSUB + 0.5 * gh),
        ),
    ]
    
    sim = mp.Simulation(
        resolution=RESOLUTION,
        cell_size=cell_size,
        boundary_layers=PML_LAYERS,
        geometry=geometry,
        k_point=K_POINT,
        sources=sources,
        symmetries=SYMMETRIES,
    )
    
    flux_obj = sim.add_flux(FCEN, 0, 1, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy)))
    sim.run(until_after_sources=200)
    
    res = sim.get_eigenmode_coefficients(flux_obj, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y)
    coeffs = res.alpha
    
    mode_tran = abs(coeffs[0, 0, 0]) ** 2 / input_flux[0]
    mode_phase = np.angle(coeffs[0, 0, 0])
    if mode_phase > 0:
        mode_phase -= 2 * np.pi
    
    return mode_tran, mode_phase


def main():
    api_key = os.getenv("OPTIX_API_KEY", "your_api_key_here")
    
    client = Optixlog(api_key=api_key)
    project = client.project(name="MeepExamples", create_if_not_exists=True)

    gp = 0.3  # periodicity
    gh = 1.8  # height
    gdc = np.linspace(0.1, 0.9, 30)
    
    run = project.run(
        name="metasurface_lens_design",
        config={
            "simulation_type": "metasurface_lens",
            "grating_period": gp,
            "grating_height": gh,
            "focal_length": FOCAL_LENGTH,
            "wavelength": LCEN,
            "resolution": RESOLUTION
        }
    )
    
    run.log(step=0, grating_period=gp, grating_height=gh, focal_length=FOCAL_LENGTH)

    # Phase map calculation
    print("⚡ Computing phase map...")
    mode_tran = np.empty(gdc.size)
    mode_phase = np.empty(gdc.size)
    
    for n in range(gdc.size):
        mode_tran[n], mode_phase[n] = grating_unit_cell(gp, gh, gdc[n])
        if n % 5 == 0:
            print(f"   Computed {n+1}/{gdc.size} duty cycles")

    phase_range = mode_phase.max() - mode_phase.min()
    
    run.log(step=1, phase="phase_map_completed",
            phase_range=float(phase_range),
            min_transmittance=float(np.min(mode_tran)),
            max_transmittance=float(np.max(mode_tran)))

    # Create phase map plot
    fig1, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    axes[0].plot(gdc, mode_tran, 'bo-', markersize=4)
    axes[0].set_xlabel("Grating Duty Cycle")
    axes[0].set_ylabel("Transmittance")
    axes[0].set_title("Transmittance vs Duty Cycle")
    axes[0].grid(True, alpha=0.3)
    
    axes[1].plot(gdc, mode_phase, 'rs-', markersize=4)
    axes[1].set_xlabel("Grating Duty Cycle")
    axes[1].set_ylabel("Phase (radians)")
    axes[1].set_title("Phase vs Duty Cycle")
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    run.log_matplotlib("phase_map", fig1)
    plt.close(fig1)
    print(f"   Phase range: {phase_range:.2f} radians")

    # Log summary
    run.log(step=2,
            phase_range_radians=phase_range,
            simulation_completed=True)
    
    print(f"\n✅ Phase map complete!")
    print(f"   Use this map to design metasurface lenses")


if __name__ == "__main__":
    main()

Key Concepts

Unit Cell Characterization

Extract phase and transmittance for each duty cycle:

res = sim.get_eigenmode_coefficients(flux_obj, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y)
mode_tran = abs(res.alpha[0, 0, 0]) ** 2 / input_flux[0]
mode_phase = np.angle(res.alpha[0, 0, 0])

Phase Unwrapping

Ensure continuous phase values:

if mode_phase > 0:
    mode_phase -= 2 * np.pi

Lens Phase Profile

Required phase for focusing:

phase_required = 2 * np.pi / wavelength * (f - np.sqrt(x**2 + f**2))

Expected Results

Phase Map

Phase varies monotonically with duty cycle
Total phase range ~π to 2π (depending on height)
Transmittance varies but typically > 50%

Focusing

With sufficient phase range:

Focal spot appears at designed distance
Spot size limited by diffraction (lens aperture)

Design Guidelines

Phase Range	Lens Quality
< π	Poor, multiple foci
π - 1.5π	Moderate, some aberrations
> 1.5π	Good, clean focusing
2π	Ideal, full control

OptixLog Integration

Logged Metrics

Metric	Description
`grating_period`	Unit cell period
`grating_height`	Ridge height
`phase_range`	Total achievable phase
`focal_length`	Target focus distance

Logged Plots

phase_map — Transmittance and phase vs duty cycle

Metasurface Lens