TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / Metamaterials, Gratings, and Other Periodic Structures / Optical skyrmion generation using metasurfaces

Optical skyrmion generation using metasurfaces

Sign up now to simulate .ipynb

Optical Stokes skyrmions are a class of structured light beams whose polarization state, mapped point-by-point across the beam’s transverse plane, sweeps out the entire surface of the Poincaré sphere. The local polarization state is captured by the normalized Stokes vector \(\mathbf{s} = (\hat{x}S_1+\hat{y}S_2+\hat{z}S_3)/S_0\), with

\[ S_0 = |\psi_R|^2 + |\psi_L|^2, \\ S_1 = 2 \mathrm{Re}(\psi_L^* \psi_R), \\ S_2 = -2 \mathrm{Im}(\psi_L^* \psi_R), \\ S_3 = |\psi_R|^2 - |\psi_L|^2,\]

where \(\psi_L\) and \(\psi_R\) are the complex amplitudes of the left- and right-handed circular polarization components at each transverse position. The skyrmion number is the topological charge derived from the resulting polarization texture:

\[N_{\mathrm{sk}} \equiv \frac{1}{4\pi}\iint_\Omega \mathbf{s}\cdot\left(\frac{\partial \mathbf{s}}{\partial x}\times\frac{\partial \mathbf{s}}{\partial y}\right)\mathrm{d}x\,\mathrm{d}y.\]

Mata-Cervera et al. generate propagation-invariant skyrmions using a single dielectric metasurface made of elliptical a-Si pillars on a SiO2 substrate. Each pillar acts as a birefringent, spin-orbit-coupled meta-atom whose two principal diameters \((D_1, D_2)\) set the co-polarized (CoP) and cross-polarized (CrP) transmission amplitude and phase. Sweeping \((D_1, D_2)\) along an isophase curve \(D_1 = f(D_2)\) tunes the CoP/CrP conversion ratio while keeping the dynamic phase constant, so the generated skyrmion’s topology stays invariant during propagation.

This notebook reproduces select results from the paper through the following workflow:

  1. Unit-cell design: batch-simulate the unit cell over \((D_1, D_2)\) to map its transmission response and extract the isophase curve.
  2. Metasurface simulation: assemble a finite-aperture metasurface along the isophase curve and simulate it alongside a non-isophase reference design.
  3. Results: reconstruct the polarization texture of the generated beam and evaluate its topological charge \(N_\mathrm{sk}\).

Reference: N. Mata-Cervera et al., Nanophotonics 14(23), 4069–4077 (2025) DOI:10.1515/nanoph-2024-0736.

Note: This notebook costs ~15 FlexCredits for unit-cell batch simulations and ~25 FlexCredits for full metasurface simulations.

schematic

Simulation Setup

# Core numerical, plotting, and Tidy3D imports
import math

import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
from matplotlib.patches import Circle
from scipy.special import assoc_laguerre
from tidy3d import web

# Silence solver log output
td.config.logging.level = "ERROR"

Define the base parameters used throughout the notebook: materials, source wavelength, unit-cell geometry (period, pillar height, and isophase curve coefficients \(\alpha,\beta,\gamma\)), metasurface aperture diameter, and the beam waists of the incident and target modes.

# Source wavelength and frequency
wl0 = 1.55
f0 = td.C_0 / wl0
fwidth = f0 / 4

# Material refractive indices and media
n_aSi = 3.505
n_SiO2 = 1.444
aSi = td.Medium(permittivity=n_aSi**2)
SiO2 = td.Medium(permittivity=n_SiO2**2)

# Unit cell geometry
nm = 1e-3
P = 650 * nm  # unit cell period
h = 800 * nm  # pillar height

# 1D parameter sweep along the iso-phase curve: D1 = alpha/D2^beta + gamma
alpha = 0.069683
beta = 1.185457
gamma = 0.078194

# Metasurface aperture diameter (scaled down from the paper's 200 um) and pillar-grid half-count
D_MS = 60.0  # 200 in the paper.
N = int(D_MS / 2 / P)

# Beam waist of the incident mode and of the two target output modes
w0 = D_MS / 4
w1 = 0.8 * w0
w2 = 0.76 * w0

Unit-Cell Design

Define helper functions to build an elliptical pillar geometry and a periodic unit-cell simulation with a normal-incidence plane wave source and transmission monitor.

# Substrate thickness, buffer layers, and overall simulation domain size
t_sub = 2 * wl0 / n_SiO2
z_buffer1 = wl0 / n_SiO2 / 2
z_buffer2 = wl0 / 2

sim_size = (P, P, h + z_buffer1 + z_buffer2)
sim_center_z = -z_buffer1 + sim_size[2] / 2


# Build an elliptical (optionally rotated) pillar, extruded to height h
def make_pillar(D1, D2, center=(0, 0), N_point=360, name="pillar", rotate=0):
    dtheta = 2 * np.pi / N_point
    theta = dtheta * np.arange(N_point)
    x = D1 / 2 * np.cos(theta)
    y = D2 / 2 * np.sin(theta)
    x_rot = np.cos(rotate) * x - np.sin(rotate) * y
    y_rot = np.sin(rotate) * x + np.cos(rotate) * y
    vertices = np.column_stack([x_rot + center[0], y_rot + center[1]])

    pillar = td.Structure(
        geometry=td.PolySlab(vertices=vertices, slab_bounds=(0, h), axis=2),
        medium=aSi,
        name=name,
    )
    return pillar


# SiO2 substrate slab
substrate = td.Structure(
    geometry=td.Box(center=(0, 0, -t_sub / 2), size=(td.inf, td.inf, t_sub)),
    medium=SiO2,
    name="substrate",
)


# Assemble a periodic unit-cell simulation for a single pillar
def make_unitcell_sim(D1, D2, aux_mnt=False):
    pillar = make_pillar(D1, D2)

    # Normally incident, x-polarized plane wave source
    pw_source = td.PlaneWave(
        center=(0, 0, -z_buffer1 / 4),
        size=(td.inf, td.inf, 0),
        source_time=td.GaussianPulse(freq0=f0, fwidth=fwidth),
        direction="+",
    )
    # Transmission monitor above the pillar
    mnt_trans = td.FieldMonitor(
        center=(0, 0, h + z_buffer2 / 2),
        size=(td.inf, td.inf, 0),
        freqs=[f0],
        name="trans",
    )
    # Optional reflection and field-profile monitors, for visualization only
    mnts_aux = (
        [
            td.FieldMonitor(
                center=(0, 0, -z_buffer1 / 2),
                size=(td.inf, td.inf, 0),
                freqs=[f0],
                name="refl",
            ),
            td.FieldMonitor(
                center=(0, 0, 0), size=(td.inf, 0, td.inf), freqs=[f0], name="field_xz"
            ),
            td.FieldMonitor(
                center=(0, 0, 0), size=(0, td.inf, td.inf), freqs=[f0], name="field_yz"
            ),
        ]
        if aux_mnt
        else []
    )

    # Local mesh refinement around the pillar
    dl = wl0 / n_aSi / 40
    mos = td.MeshOverrideStructure(
        geometry=pillar.geometry.bounding_box, dl=(dl, dl, dl)
    )
    grid_spec = td.GridSpec.auto(
        wavelength=wl0,
        min_steps_per_wvl=20,
        override_structures=[mos],
        snapping_points=[(0, 0, pw_source.center[2]), (0, 0, mnt_trans.center[2])],
    )
    # Periodic boundaries in-plane, PML along propagation
    boundary_spec = td.BoundarySpec(
        x=td.Boundary(minus=td.Periodic(), plus=td.Periodic()),
        y=td.Boundary(minus=td.Periodic(), plus=td.Periodic()),
        z=td.Boundary(minus=td.PML(), plus=td.PML()),
    )
    # Assemble and return the unit-cell simulation
    return td.Simulation(
        size=sim_size,
        center=(0, 0, sim_center_z),
        structures=[substrate, pillar],
        sources=[pw_source],
        monitors=[mnt_trans] + mnts_aux,
        grid_spec=grid_spec,
        boundary_spec=boundary_spec,
        run_time=td.RunTimeSpec(quality_factor=50),
        symmetry=(-1, 1, 0),  # x-pol incidence
    )

Plot a cross section of an example unit-cell simulation.

# Example unit cell for visualization
sim = make_unitcell_sim(0.5, 0.25, aux_mnt=True)

# Top-down cross section at the pillar mid-height
fig, ax = plt.subplots(1, 2)
sim.plot(z=h / 2, ax=ax[0])
sim.plot_grid(z=h / 2, ax=ax[0])

# Side cross section through the pillar
sim.plot(y=0 + 1e-4, ax=ax[1])
sim.plot_grid(y=0 + 1e-4, ax=ax[1])
ax[1].set(ylim=(-0.4, 1.4));

Run Batch Simulations

Define a dictionary of simulations over \((D_1, D_2)\): a 2D sweep on a grid as well as 1D sweep along the isophase curve.

# 2D parameter sweep on grid
D_list = np.arange(210, 541, 15) * nm
sims_grid = {
    f"unit_{n1}_{n2}": make_unitcell_sim(d1, d2)
    for n1, d1 in enumerate(D_list)
    for n2, d2 in enumerate(D_list)
}

# Iso-phase curve sampled in D2 and mapped to D1
d2_iso = np.arange(210, 333.9, 3) * nm


def d2_to_d1(d2):
    return alpha / d2**beta + gamma


def d1_to_d2(d1):
    return (alpha / (d1 - gamma)) ** (1 / beta)


d1_iso = d2_to_d1(d2_iso)

# Unit-cell sims along the iso-phase curve, (D1, D2)
sims_isophase = {
    f"unit_isophase_{s}": make_unitcell_sim(d1, d2)
    for s, (d1, d2) in enumerate(zip(d1_iso, d2_iso))
}
# Transposed sims, (D2, D1), to obtain t_yy along the same curve
sims_isophase_t = {
    f"unit_isophase_{s}_t": make_unitcell_sim(d2, d1)
    for s, (d1, d2) in enumerate(zip(d1_iso, d2_iso))
}

# Combine all sweeps into one batch
sims = sims_grid | sims_isophase | sims_isophase_t

Submit batch simulations and download results.

# Submit all unit-cell simulations as one batch job
batch_data = web.run_async(
    simulations=sims, folder_name="Skyrmion", path_dir="./data/unitcell", verbose=False
)

Transmission Coefficients

Extract \(t_{xx}\) from each symmetry-enforced unit-cell simulation, and obtain \(t_{yy}\) by transposing \(D_1 \leftrightarrow D_2\). In the circular polarization basis \(\hat{u}_\pm = (\hat{x} \pm i\hat{y})/\sqrt{2}\), the Jones matrix for LCP (\(\hat{u}_+\)) and RCP (\(\hat{u}_-\)) is

\[ \frac{1}{2} \begin{bmatrix} t_{xx}+t_{yy} & t_{xx}-t_{yy} \\ t_{xx}-t_{yy} & t_{xx}+t_{yy} \end{bmatrix}.\]

Thus, co- and cross-polarization transmission coefficients are $ t_ = (t_{xx} + t_{yy})/2$ and \(t_\mathrm{CrP} = (t_{xx} - t_{yy})/2\), respectively.

def unitcell_transmission(field_data, P):
    """Normalized complex transmission coefficient from a unit-cell Ex field monitor."""
    Ex_avg = field_data.Ex.squeeze().integrate(coord=["x", "y"]) / P**2
    return np.sqrt(0.5 * np.abs(Ex_avg) ** 2 / td.ETA_0 * (P**2)) * np.exp(
        1j * np.angle(Ex_avg)
    )


# Collect t_xx over the 2D (D1, D2) grid
t_xx_mat = np.array(
    [
        [
            unitcell_transmission(batch_data[f"unit_{n1}_{n2}"]["trans"], P)
            for n2 in range(len(D_list))
        ]
        for n1 in range(len(D_list))
    ]
)

# Derive t_yy by transpose (D1 <-> D2), then form the CoP/CrP combinations
t_yy_mat = t_xx_mat.T
t_cop_mat = 0.5 * (t_xx_mat + t_yy_mat)
t_crp_mat = 0.5 * (t_xx_mat - t_yy_mat)

# Collect t_xx and t_yy along the iso-phase curve
n_iso = len(d1_iso)
t_xx_iso = np.array(
    [
        unitcell_transmission(batch_data[f"unit_isophase_{s}"]["trans"], P)
        for s in range(n_iso)
    ]
)
t_yy_iso = np.array(
    [
        unitcell_transmission(batch_data[f"unit_isophase_{s}_t"]["trans"], P)
        for s in range(n_iso)
    ]
)

# CoP/CrP transmission along the iso-phase curve
t_cop_iso = 0.5 * (t_xx_iso + t_yy_iso)
t_crp_iso = 0.5 * (t_xx_iso - t_yy_iso)

For LCP incidence, the position-dependent output field becomes \[\mathbf{E}_\mathrm{out}(\mathbf{r}) = e^{i\xi(\mathbf{r})} \left[\hat{u}_+ \cos\frac{\Delta\phi(\mathbf{r})}{2} + i\hat{u}_- e^{i2\theta(\mathbf{r})}\sin\frac{\Delta\phi(\mathbf{r})}{2}\right],\] where \(\xi \equiv (\phi_x + \phi_y)/2\) is a global phase factor, \(\Delta\phi = \phi_x - \phi_y\) is a birefringence factor, and \(\theta\) is a geometrical phase corresponding to the pillar rotation, for near-unitary transmission \(t_{xx} = e^{i\phi_x}\) and \(t_{yy} = e^{i\phi_y}\). Keeping \(\xi\) constant over the metasurface is thus crucial to suppress unwanted wavefront perturbation while using \(\Delta\phi\) and \(\theta\) as amplitude and phase degrees of freedom. Plot the isophase (solid) and non-isophase (dashed) paths in the \((D_1, D_2)\) parameter space.

# Meshgrid of the (D1, D2) sweep, for plotting
D1_grid, D2_grid = np.meshgrid(D_list, D_list, indexing="ij")
fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, tight_layout=True)

# Non-isophase reference path: D2 fixed at its minimum value
d2_noniso = D2_grid[:-1, 0]
d1_noniso = D1_grid[:-1, 0]
t_cop_noniso = t_cop_mat[:-1, 0]
t_crp_noniso = t_crp_mat[:-1, 0]

# Columns to plot: xx, CoP, and CrP transmission
columns = [
    ("xx", t_xx_mat, None),
    ("CoP", t_cop_mat, "r"),
    ("CrP", t_crp_mat, "b"),
]

# Plot magnitude/phase maps for each column, with the iso- and non-iso paths overlaid
for col, (sub, t_mat, iso_color) in enumerate(columns):
    ax[0, col].set(aspect="equal", title=rf"$|t_{{{sub}}}|$")
    ax[1, col].set(aspect="equal", title=rf"$\angle t_{{{sub}}}$")
    pc_abs = ax[0, col].pcolormesh(D2_grid, D1_grid, np.abs(t_mat), vmin=0, vmax=1)
    pc_phase = ax[1, col].pcolormesh(
        D2_grid, D1_grid, np.angle(t_mat), vmin=-np.pi, vmax=np.pi, cmap=plt.cm.twilight
    )
    if iso_color is not None:
        ax[0, col].plot(d2_iso, d1_iso, iso_color, lw=2)
        ax[1, col].plot(d2_iso, d1_iso, iso_color, lw=2)
        ax[0, col].plot(d2_noniso, d1_noniso, iso_color + "--", lw=2)
        ax[1, col].plot(d2_noniso, d1_noniso, iso_color + "--", lw=2)

# Shared colorbars for amplitude and phase
cax_abs = ax[0, -1].inset_axes([1.05, 0, 0.05, 1])
cax_phase = ax[1, -1].inset_axes([1.05, 0, 0.05, 1])
plt.colorbar(pc_abs, cax=cax_abs)
plt.colorbar(pc_phase, cax=cax_phase)

# Axis labels
for row in range(2):
    ax[row, 0].set(ylabel=r"$D_1$ (um)")
for col in range(3):
    ax[1, col].set(xlabel=r"$D_2$ (um)")

# Intensity and phase vs D1, with corresponding D2 shown on a secondary x-axis
fig, ax = plt.subplots(2, 1, sharex=True, tight_layout=True)
ax2 = [ax[0].twiny(), ax[1].twiny()]

xlim = (d1_iso.min(), d1_iso.max())
d1_xticks = np.arange(0.34, 0.53, 0.02)
d2_xticks = np.arange(0.21, 0.335, 0.02)

# Co- and cross-polarized transmitted intensity along the iso-curve
ax[0].plot(d1_iso, np.abs(t_cop_iso) ** 2, "r", label="Co-pol")
ax[0].plot(d1_iso, np.abs(t_crp_iso) ** 2, "b", label="Cross-pol")
ax[0].legend(frameon=False)
ax[0].set(xlim=xlim, xticks=d1_xticks, ylim=(0, 1), ylabel=r"Intensity, $|t|^2$")

# Corresponding transmission phase
ax[1].plot(d1_iso, np.angle(t_cop_iso), "r")
ax[1].plot(d1_iso, np.angle(t_crp_iso), "b")
ax[1].set(
    xlim=xlim,
    xticks=d1_xticks,
    xlabel=r"$D_1$ (um)",
    ylim=(-np.pi, np.pi),
    ylabel=r"Phase, $\angle t$ (rad)",
)

# Map D1 ticks to their corresponding D2 values on the top axis
d2_xticklabels = [f"{xt:.2f}" for xt in d2_xticks]
for a2, xticklabels in zip(ax2, [d2_xticklabels, []]):
    a2.tick_params(axis="x", colors="g")
    a2.spines["top"].set_color("g")
    a2.set(xlim=xlim, xticks=d2_to_d1(d2_xticks), xticklabels=xticklabels)
ax2[0].set_xlabel(r"$D_2$ (um)", color="g");

Meta-Unit Interpolation

Define an interpolation helper for a meta-unit library that maps the desired response \(S_3\) to geometrical parameters \((D_1, D_2)\).

S3_iso = (np.abs(t_crp_iso) ** 2 - np.abs(t_cop_iso) ** 2) / (
    np.abs(t_crp_iso) ** 2 + np.abs(t_cop_iso) ** 2
)


def S2D(s):
    d2_interp = np.interp(s, S3_iso[::-1], d2_iso[::-1])
    d1_interp = d2_to_d1(d2_interp)
    return d1_interp, d2_interp


# Compare interpolated D1/D2 curves to the simulated data
fig, ax = plt.subplots()
s = np.linspace(-1, 1, 101)
ax.plot(S3_iso, d1_iso, "kx", label=r"$D_1$, data")
ax.plot(s, S2D(s)[0], "m", label=r"$D_1$, interp")
ax.plot(S3_iso, d2_iso, "k.", label=r"$D_2$, data")
ax.plot(s, S2D(s)[1], "c", label=r"$D_2$, interp")
ax.legend(frameon=False)

ax.set(
    xlabel=r"$S_3/S_0$",
    xlim=(-1, 1),
    ylabel=r"$D$ (um)",
    ylim=(0.2, 0.53),
);

Define an interpolated non-isophase path in the same way, and separately interpolate the global phase \(\xi\) as a function of \(D_1\) so it can be compensated for in the reference design (optional).

# Compute S3 (LCP input) for the non-isophase path, and interpolate D1(S3)
S3_noniso = (np.abs(t_crp_noniso) ** 2 - np.abs(t_cop_noniso) ** 2) / (
    np.abs(t_crp_noniso) ** 2 + np.abs(t_cop_noniso) ** 2
)  # Assume LCP input

# Map a desired S3 to (D1, D2) for the non-isophase (reference) design


def S2D_noniso(s):
    d1_interp = np.interp(s, S3_noniso, d1_noniso)
    d2 = d2_noniso[0].item() * np.ones_like(s)
    return d1_interp, d2


# Compare interpolated D1 to the simulated data
fig, ax = plt.subplots(1, 2, sharey=True, tight_layout=True)
s = np.linspace(-1, 1, 101)
ax[0].plot(S3_noniso, d1_noniso, "kx", label=r"$D_1$, data")
ax[0].plot(s, S2D_noniso(s)[0], "m", label=r"$D_1$, interp")
ax[0].legend(frameon=False)

ax[0].set(
    xlabel=r"$S_3/S_0$",
    xlim=(-1, 1),
    ylabel=r"$D$ (um)",
    ylim=(0.2, 0.53),
    title=r"$S \rightarrow D_1$ interp",
)
# Global phase xi accumulated along the non-isophase path
t_xx_noniso = t_xx_mat[:-1, 0]
t_yy_noniso = t_yy_mat[:-1, 0]
xi_noniso = np.unwrap(np.angle(t_xx_noniso) + np.angle(t_yy_noniso)) / 2


# Interpolate the global phase as a function of D1, to compensate for it in the reference design
def D2Xi(d):
    return np.interp(d, d1_noniso, xi_noniso)


# Compare interpolated xi(D1) to the simulated data
ax[1].plot(xi_noniso, d1_noniso, "kx")
ax[1].plot(D2Xi(d1_noniso), d1_noniso, "m")
ax[1].set(xlabel=r"Global phase, $\xi$", title=r"$D_1 \rightarrow \xi$ interp");

Metasurface Simulation

Build a finite-aperture metasurface that converts an incident Gaussian beam into a target superposition of Laguerre-Gaussian modes, then simulate both the isophase design and a non-isophase reference for comparison.

Note: We use the metasurface aperture diameter \(D_\mathrm{MS} = 60\) instead of 200 µm in the paper to simplify the simulations.

# Laguerre-Gaussian beam profile, plus its Rayleigh range and divergence angle
def LaguerreGaussian(rho, phi, z, w0, l=0, p=0, k0=2 * np.pi / wl0):
    absl = np.abs(l)
    C_norm = np.sqrt(2 * math.factorial(p) / np.pi / math.factorial(p + absl))
    z_R = 0.5 * k0 * w0**2
    w_z = w0 * np.sqrt(1 + (z / z_R) ** 2)
    invR_z = z / (z**2 + z_R**2)
    psi_z = np.atan(z / z_R) * (2 * p + absl + 1)
    theta = 2 / (k0 * w0)
    u = (
        (C_norm / w_z)
        * (np.sqrt(2) * rho / w_z) ** absl
        * np.exp(-((rho / w_z) ** 2))
        * assoc_laguerre(2 * (rho / w_z) ** 2, p, absl)
        * np.exp(1j * (k0 * z + k0 * invR_z * rho**2 / 2 + l * phi - psi_z))
    )
    return u, (z_R, theta)


# Pillar center coordinates on the periodic grid
x_c, y_c = np.meshgrid(
    P * np.arange(-N, N + 1), P * np.arange(-N, N + 1), indexing="ij"
)
Rho = np.sqrt(x_c**2 + y_c**2)
Phi = np.atan2(y_c, x_c)

# Incident LG00 mode and target LG00/LG02 output modes
psi0, (z_R, _) = LaguerreGaussian(
    Rho,
    Phi,
    0,
    w0,
    l=0,
    p=0,
)
psi1, _ = LaguerreGaussian(
    Rho,
    Phi,
    0,
    w1,
    l=0,
    p=0,
)
psi2, (z_R2, _) = LaguerreGaussian(
    Rho,
    Phi,
    0,
    w2,
    l=2,
    p=0,
)

# Intensity profiles of the incident and target modes
incident = np.abs(psi0) ** 2
target1 = np.abs(psi1) ** 2
target2 = np.abs(psi2) ** 2

# Normalized modulation (mode-splitting ratio) the metasurface must apply
norm1 = np.abs(psi1) ** 2 / (np.abs(psi1) ** 2 + np.abs(psi2) ** 2)
norm2 = np.abs(psi2) ** 2 / (np.abs(psi1) ** 2 + np.abs(psi2) ** 2)

# Output intensity after applying the modulation to the incident beam
out1 = incident * norm1
out2 = incident * norm2

# Overlap integrals giving the coupling efficiency into each target mode
c1 = (np.sqrt(out1 * target1) * P**2).sum()
c2 = (np.sqrt(out2 * target2) * P**2).sum()
coupling_efficiency = np.abs(c1) ** 2 + np.abs(c2) ** 2

# Visualize incident, target, and modulated output profiles
fig, ax = plt.subplots(3, 1, sharex=True, tight_layout=True)
ax[0].plot(
    x_c[y_c == 0],
    incident[y_c == 0],
    "k",
    label=r"$\psi_\mathrm{inc}$, $\mathrm{LG}_{00}$, $w_0$",
)
ax[0].fill_between(
    x_c[y_c == 0],
    target1[y_c == 0],
    label=r"$\psi_\mathrm{tar1}$, $\mathrm{LG}_{00}$, $w_1$",
    alpha=0.5,
)
ax[0].fill_between(
    x_c[y_c == 0],
    target2[y_c == 0],
    label=r"$\psi_\mathrm{tar2}$, $\mathrm{LG}_{02}$, $w_2$",
    alpha=0.5,
)

ax[0].set(
    ylabel=r"$|\psi|^2$",
    ylim=(0, target1.max()),
    title="Input and target output profiles",
)
ax[0].legend(frameon=False)

ax[1].plot(x_c[y_c == 0], norm1[y_c == 0])
ax[1].plot(x_c[y_c == 0], norm2[y_c == 0])
ax[1].set(
    ylabel=r"$|t|^2$", ylim=(0, 1), title=r"Normalized target (modulation): $t_{1,2}$"
)

ax[2].plot(x_c[y_c == 0], out1[y_c == 0], "k--", label=r"$\psi_\mathrm{out,1}$")
ax[2].plot(x_c[y_c == 0], out2[y_c == 0], "k--", label=r"$\psi_\mathrm{out,2}$")
ax[2].fill_between(
    x_c[y_c == 0],
    target1[y_c == 0] * np.max(out1) / np.max(target1),
    label=r"$\psi_\mathrm{tar1}$",
    alpha=0.25,
)
ax[2].fill_between(
    x_c[y_c == 0],
    target2[y_c == 0] * np.max(out2) / np.max(target2),
    label=r"$\psi_\mathrm{tar2}$",
    alpha=0.25,
)
ax[2].set(
    ylabel=r"$|\psi|^2$",
    ylim=(0, out1.max()),
    title=r"Output profile: $\psi_\mathrm{out} = \psi_\mathrm{inc} * t_{1,2}$",
)

ax[2].legend(frameon=False)
ax[2].set(xlabel=r"$\rho$ (um)", xlim=(-D_MS / 2, D_MS / 2))
print(f"Coupling efficiency into output LG modes: {coupling_efficiency * 100:.2f}%")
Coupling efficiency into output LG modes: 98.91%

Assign each pillar’s \((D_1, D_2, \theta)\) from the interpolated \(S \to D\) maps to realize the target LG-mode superposition, then assemble the isophase and non-isophase pillar arrays into full metasurface simulations with far-field projection monitors at multiple propagation distances.

Since each pillar’s diameters depend only on radius and its rotation \(\theta\) scales with the local azimuthal angle, the metasurface layout exhibits 90° rotational symmetry. A single \(x\)-polarized GaussianBeam source with mirror symmetry enforced (symmetry=(-1,1,0)) is therefore enough: the \(y\)-polarized response needed to reconstruct the circular (LCP/RCP) output is obtained by rotating the simulated \(x\)-polarized far field by 90°, without running a second simulation.

# Build pillar geometries for the isophase and non-isophase designs from the interpolated S->D maps
def pillar_geometry(D1, D2, center, rotate):
    return make_pillar(D1, D2, center=center, rotate=rotate).geometry


def to_pillar_array(pillars):
    return td.Structure(
        geometry=td.GeometryGroup(geometries=pillars), medium=aSi, name="pillar array"
    )


pillars_sk = []
pillars_ref = []
for x, y, radius, phase, s in zip(
    x_c.reshape(-1),
    y_c.reshape(-1),
    Rho.reshape(-1),
    np.angle(psi2).reshape(-1),
    (norm2 - norm1).reshape(-1),
):
    if radius >= D_MS / 2:
        continue

    D1, D2 = S2D(s)
    D1_ref, D2_ref = S2D_noniso(s)
    pillars_sk.append(pillar_geometry(D1, D2, (x, y), phase / 2))
    pillars_ref.append(pillar_geometry(D1_ref, D2_ref, (x, y), phase / 2))

# Group pillars into single structures
pillar_array_sk = to_pillar_array(pillars_sk)
pillar_array_ref = to_pillar_array(pillars_ref)

# x-polarized Gaussian beam source matching the incident LG00 mode waist
gauss_xpol = td.GaussianBeam(
    center=(0, 0, -z_buffer1 / 4),
    size=(td.inf, td.inf, 0),
    source_time=td.GaussianPulse(freq0=f0, fwidth=fwidth),
    direction="+",
    pol_angle=0,
    waist_radius=w0,
    waist_distance=+z_buffer1 / 4,
)

# Far-field projection distances, in units of the Rayleigh range
proj_distances = z_R * np.array(
    [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.5, 2]
)

# Far-field monitor at each projection distance
mnts_far = []
for ii, dist in enumerate(proj_distances):
    r_mnt = D_MS * np.sqrt(1 + (dist / z_R) ** 2)
    N_mnt = int(r_mnt / (wl0 / 5))
    mnts_far.append(
        td.FieldProjectionCartesianMonitor(
            center=(0, 0, h + z_buffer2 / 2),
            size=(td.inf, td.inf, 0),
            freqs=[f0],
            normal_dir="+",
            custom_origin=(0, 0, 0),
            proj_axis=2,
            proj_distance=dist,
            x=wl0 / 5 * np.arange(-N_mnt, N_mnt + 1),
            y=wl0 / 5 * np.arange(-N_mnt, N_mnt + 1),
            far_field_approx=False,
            name=f"far_{ii}",
        )
    )

# Assemble the isophase metasurface simulation
ms_sim_size = (D_MS + 2 * P, D_MS + 2 * P, h + z_buffer1 + z_buffer2)
ms_sim_center_z = -z_buffer1 + ms_sim_size[2] / 2
sim_sk = td.Simulation(
    size=ms_sim_size,
    center=(0, 0, ms_sim_center_z),
    structures=[pillar_array_sk, substrate],
    sources=[gauss_xpol],
    monitors=mnts_far,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=20),
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
    run_time=td.RunTimeSpec(quality_factor=10),
    symmetry=(-1, 1, 0),
)

# Reference simulation using the non-isophase pillar array
sim_ref = sim_sk.updated_copy(structures=[pillar_array_ref, substrate])

Submit both metasurface simulations as a batch and download results.

# Submit both metasurface simulations as a batch job
batch_data = web.run_async(
    simulations={"skyrmion": sim_sk, "ref": sim_ref},
    folder_name="Skyrmion",
    path_dir="./data/metasurface",
    verbose=False,
)

Results

Postprocessing Field Data

Convert projected far fields into circular-polarization amplitudes, compute the normalized Stokes vector \(\mathbf{s}\), and integrate the skyrmion density to get \(N_\mathrm{sk}\), optionally plotting the LCP/RCP intensities and polarization texture.

def stokes_to_texture(comp1, comp2, lightness):
    """Encode a polarization state as an RGB texture: hue from atan2(comp2, comp1) via the HSV colormap, with lightness fading the hue to black (+1) or white (-1)."""
    hue = plt.cm.hsv((np.atan2(comp2, comp1) + np.pi) / (2 * np.pi))[..., :3]
    lightness = np.asarray(lightness).reshape(*hue.shape[:-1], 1)
    return np.where(
        lightness > 0,
        (1 - np.abs(lightness)) * hue,
        np.abs(lightness) + (1 - np.abs(lightness)) * hue,
    )


def postprocess(fielddata, r_int_factor=0.7, ax=None):
    # Extract far-field spherical components and coordinates
    Er = fielddata.Er.squeeze()
    Etheta = fielddata.Etheta.squeeze()
    Ephi = fielddata.Ephi.squeeze()
    x, y, z = Er.x, Er.y, Er.z
    _, theta, phi = td.Geometry.car_2_sph(x, y, z)
    rho = np.sqrt(x**2 + y**2)

    # Convert to Cartesian field components for the x-polarized source
    Ex_xpol, Ey_xpol, _ = td.Geometry.sph_2_car_field(Er, Etheta, Ephi, theta, phi)

    # Recover the y-polarized response by rotating the x-polarized far field 90 degrees
    Ex_ypol = -np.rot90(Ey_xpol, k=1)
    Ey_ypol = np.rot90(Ex_xpol, k=1)

    # Combine x- and y-polarized responses into a circular-polarization basis
    Ex = (Ex_xpol + 1j * Ex_ypol) / np.sqrt(2)
    Ey = (Ey_xpol + 1j * Ey_ypol) / np.sqrt(2)

    # LCP and RCP complex amplitudes
    psi_lcp = (Ex - 1j * Ey) / np.sqrt(2)
    psi_rcp = (Ex + 1j * Ey) / np.sqrt(2)

    # Normalized Stokes vector
    S0 = np.abs(psi_rcp) ** 2 + np.abs(psi_lcp) ** 2
    S1 = 2 * np.real(psi_lcp.conj() * psi_rcp)
    S2 = -2 * np.imag(psi_lcp.conj() * psi_rcp)
    S3 = np.abs(psi_rcp) ** 2 - np.abs(psi_lcp) ** 2
    S = np.stack([S1 / S0, S2 / S0, S3 / S0], axis=2)

    # Skyrmion density from the Stokes vector's spatial gradients
    dSdx = np.gradient(S, x, axis=0)
    dSdy = np.gradient(S, y, axis=1)

    density = np.sum(S * np.linalg.cross(dSdx, dSdy), axis=2)
    density = td.SpatialDataArray(density, coords={"x": x, "y": y})

    # Integrate the density within a radius scaled to the beam size, to get the skyrmion number
    r_int = D_MS / 2 * r_int_factor * np.sqrt(1 + (z / z_R2) ** 2)
    N_sk = (density * (rho < r_int)).integrate(coord=["x", "y"]) / (4 * np.pi)

    # Optional: plot LCP/RCP intensities and the Stokes texture
    if ax is not None:
        for a in ax:
            a.set(
                aspect="equal",
                xlim=(-r_int * 1.25, r_int * 1.25),
                ylim=(-r_int * 1.25, r_int * 1.25),
            )

        ax[0].pcolormesh(x, y, np.abs(psi_lcp) ** 2, cmap=plt.cm.inferno)
        ax[1].pcolormesh(x, y, np.abs(psi_rcp) ** 2, cmap=plt.cm.inferno)

        texture = stokes_to_texture(S1, S2, S3 / S0)
        ax[2].pcolormesh(x, y, texture)
        [a.add_patch(Circle((0, 0), r_int, color="w", fill=False, ls="--")) for a in ax]

    return (psi_lcp, psi_rcp, density), N_sk

Polarization Texture Visualizations

For each projection distance, plot the LCP/RCP intensities and polarization texture of the isophase design (left) against the non-isophase reference (right).

The dashed circle marks the integration domain radius r_int used below to compute the skyrmion number, chosen arbitrarily to cover roughly the extent of the \(\mathrm{LG}_{02}\) mode while excluding numerical noise from the near-to-far-field transform. It grows with propagation distance as \(\sqrt{1+(z/z_{R2})^2}\), tracking the beam’s divergence.

# Grid of subplots: one row per projection distance, isophase design vs. reference side by side
fig, ax = plt.subplots(
    len(proj_distances), 6, figsize=(6, len(proj_distances)), tight_layout=True
)

N_sk_list = []
N_sk_ref_list = []
# Postprocess and plot each projection distance for both designs
for ii, z in enumerate(proj_distances):
    fielddata_sk = batch_data["skyrmion"][f"far_{ii}"]
    fielddata_ref = batch_data["ref"][f"far_{ii}"]

    _, N_sk = postprocess(fielddata_sk, ax=ax[ii, :3])
    _, N_sk_ref = postprocess(fielddata_ref, ax=ax[ii, 3:6])
    [a.set(xticks=[], yticks=[]) for a in ax[ii]]

    N_sk_list.append(N_sk)
    N_sk_ref_list.append(N_sk_ref)
    ax[ii, 0].set_ylabel(rf"$z = ${z / z_R:.1f}$z_R$")

# Column titles
ax[0, 0].set_title(r"LCP")
ax[0, 1].set_title(r"RCP")
ax[0, 2].set_title(r"Texture")

ax[0, 3].set_title(r"LCP")
ax[0, 4].set_title(r"RCP")
ax[0, 5].set_title(r"Texture")

# Colormaps for intensities and polarization states
fig, ax = plt.subplots(figsize=(3, 1.5))
Theta_sphere, Phi_sphere = np.meshgrid(
    np.linspace(0, np.pi, 361), np.linspace(-np.pi, np.pi, 721), indexing="ij"
)

X_sphere = np.sin(Theta_sphere) * np.cos(Phi_sphere)
Y_sphere = np.sin(Theta_sphere) * np.sin(Phi_sphere)
Z_sphere = np.cos(Theta_sphere)

texture = stokes_to_texture(X_sphere, Y_sphere, Z_sphere)
ax.pcolormesh(Phi_sphere, Z_sphere, texture)
ax.set(
    xlabel=r"$\tan^{-1}(S_2/S_1)$",
    ylabel=r"$S_3/S_0$",
    xticks=(-np.pi, 0, np.pi),
    xticklabels=[r"$-\pi$", 0, r"$\pi$"],
    yticks=(-1, 0, 1),
    title="Intensity & Polarization",
)

axins = ax.inset_axes([0, 1.5, 1, 0.2])
cbar = plt.colorbar(
    plt.cm.ScalarMappable(cmap=plt.cm.inferno), cax=axins, orientation="horizontal"
)
cbar.ax.set_xlabel(r"$|\psi_{L,R}|^2$ (Arb.U.)");

Skyrmion Number Invariance

Plot the skyrmion number versus propagation distance to confirm it stays constant for the isophase design but drifts for the non-isophase reference.

# Skyrmion number vs. propagation distance, for both designs
fig, ax = plt.subplots()

ax.plot(proj_distances / z_R, N_sk_list, "ko-", label=r"Isophase design")
ax.plot(proj_distances / z_R, N_sk_ref_list, "go-", label=r"non-isophase design")
ax.set(
    xlabel=r"Projection distance, $z/z_R$", ylabel=r"Skyrmion number, $N_\mathrm{sk}$"
)

ax.legend(frameon=False);

The isophase design keeps the global phase \(\xi\) constant across the aperture, so its skyrmion number \(N_\mathrm{sk}\approx 1.9\) (theoretically, 2) stays invariant with propagation, while the non-isophase reference has a varying \(\xi\) that causes \(N_\mathrm{sk}\) to drift with distance. This also shows up in the polarization texture: the isophase design’s texture reaches full black (pure RCP, \(S_3/S_0=1\)) at a certain radius, while the non-isophase texture never reaches black, meaning it fails to fully sample the Poincaré sphere and cannot support a robust skyrmion.

Example library

More Examples

Read More
Dielectric metasurface absorber
Dielectric metasurface absorber
Gradient metasurface reflector
Gradient metasurface reflector
Graphene metamaterial absorber
Graphene metamaterial absorber

No matching results found.

Previous Tutorial Next Tutorial