Parallel Waveguides Optical Force

This example computes the optical force between two parallel waveguides as a function of their separation distance, analyzing both symmetric and antisymmetric coupled modes.

Overview

Optical forces in coupled waveguides have applications in:

Optomechanics: Light-driven mechanical motion
MEMS/NEMS: Optical actuation of micro/nano devices
All-optical switching: Force-induced coupling changes
Sensing: Measuring nanoscale displacements

The force can be attractive or repulsive depending on the mode symmetry.

Simulation Parameters

Parameter	Default Value	Description
`resolution`	40	Pixels per μm
`Si`	n=3.45	Silicon waveguide index
`dpml`	1.0	PML thickness
`sx`	5	Cell size in x
`sy`	3	Cell size in y
`a`	1.0	Waveguide width/height
`s`	0.05-1.0	Separation range

Physical Setup

The simulation geometry:

Two parallel waveguides: Silicon (n=3.45) square cross-section
Variable separation: Sweep from 0.05a to 1.0a
Mode symmetries: Symmetric (even) and antisymmetric (odd) modes
Force calculation: Maxwell stress tensor integration

The optical force per unit length is normalized by the optical power: F/L = (force per length) × (c/P), giving a dimensionless quantity.

Python Code

"""
Parallel Waveguides Optical Force with OptixLog Integration

Computes the optical force between two parallel waveguides
as a function of separation distance.

Based on the Meep tutorial: parallel-wvgs-force.py
"""

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

api_key = os.getenv("OPTIX_API_KEY", "")
api_url = os.getenv("OPTIX_API_URL", "https://optixlog.com")
project_name = os.getenv("OPTIX_PROJECT", "MeepExamples")


def main():
    """Main simulation function for parallel waveguides force."""
    
    if not optixlog.is_master_process():
        return
    
    try:
        client = optixlog.init(
            api_key=api_key,
            api_url=api_url,
            project=project_name,
            run_name="parallel_waveguides_force",
            config={
                "simulation_type": "optical_force",
            },
            create_project_if_not_exists=True
        )
        
        # Simulation parameters
        resolution = 40
        Si = mp.Medium(index=3.45)
        dpml = 1.0
        sx = 5
        sy = 3
        a = 1.0
        k_point = mp.Vector3(z=0.5)
        
        cell = mp.Vector3(sx + 2 * dpml, sy + 2 * dpml, 0)
        pml_layers = [mp.PML(dpml)]
        
        client.log(
            step=0,
            resolution=resolution,
            silicon_index=3.45,
            waveguide_width=a,
        )
        
        def parallel_waveguide(s, xodd, sim_index):
            """Simulate parallel waveguides with separation s."""
            geometry = [
                mp.Block(
                    center=mp.Vector3(-0.5 * (s + a)),
                    size=mp.Vector3(a, a, mp.inf),
                    material=Si,
                ),
                mp.Block(
                    center=mp.Vector3(0.5 * (s + a)),
                    size=mp.Vector3(a, a, mp.inf),
                    material=Si
                ),
            ]
            
            symmetries = [
                mp.Mirror(mp.X, phase=-1 if xodd else 1),
                mp.Mirror(mp.Y, phase=-1)
            ]
            
            sim = mp.Simulation(
                resolution=resolution,
                cell_size=cell,
                geometry=geometry,
                boundary_layers=pml_layers,
                symmetries=symmetries,
                k_point=k_point,
            )
            
            sim.init_sim()
            
            EigenmodeData = sim.get_eigenmode(
                0.22,
                mp.Z,
                mp.Volume(center=mp.Vector3(), size=mp.Vector3(sx, sy)),
                2 if xodd else 1,
                k_point,
                match_frequency=False,
                parity=mp.ODD_Y,
            )
            
            fcen = EigenmodeData.freq
            sim.reset_meep()
            
            eig_sources = [
                mp.EigenModeSource(
                    src=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
                    size=mp.Vector3(sx, sy),
                    center=mp.Vector3(),
                    eig_band=2 if xodd else 1,
                    eig_kpoint=k_point,
                    eig_match_freq=False,
                    eig_parity=mp.ODD_Y,
                )
            ]
            
            sim.change_sources(eig_sources)
            
            flux_reg = mp.FluxRegion(
                direction=mp.Z, center=mp.Vector3(), size=mp.Vector3(sx, sy)
            )
            wvg_flux = sim.add_flux(fcen, 0, 1, flux_reg)
            
            force_reg1 = mp.ForceRegion(
                mp.Vector3(0.49 * s), direction=mp.X, weight=1, size=mp.Vector3(y=sy)
            )
            force_reg2 = mp.ForceRegion(
                mp.Vector3(0.5 * s + 1.01 * a), direction=mp.X, weight=-1, size=mp.Vector3(y=sy)
            )
            wvg_force = sim.add_force(fcen, 0, 1, force_reg1, force_reg2)
            
            sim.run(until_after_sources=1500)
            
            flux = mp.get_fluxes(wvg_flux)[0]
            force = mp.get_forces(wvg_force)[0]
            
            sim.reset_meep()
            return flux, force, fcen
        
        # Separation range
        s_values = np.arange(0.05, 1.05, 0.05)
        forces_odd = np.zeros(len(s_values))
        forces_even = np.zeros(len(s_values))
        fluxes_odd = np.zeros(len(s_values))
        fluxes_even = np.zeros(len(s_values))
        
        for k, s in enumerate(s_values):
            fluxes_odd[k], forces_odd[k], _ = parallel_waveguide(s, True, k*2)
            fluxes_even[k], forces_even[k], _ = parallel_waveguide(s, False, k*2+1)
        
        norm_forces_odd = -forces_odd / fluxes_odd
        norm_forces_even = -forces_even / fluxes_even
        
        client.log(
            step=1,
            num_separations=len(s_values),
            max_antisymmetric_force=float(np.max(norm_forces_odd)),
            max_symmetric_force=float(np.max(norm_forces_even))
        )
        
    except Exception as e:
        print(f"Simulation Error: {e}")


if __name__ == "__main__":
    main()

How to Run

# Set your OptixLog API key
export OPTIX_API_KEY="your_api_key_here"

# Run the force calculation
python parallel-wvgs-force.py

Results and Analysis

Force vs Separation

The simulation reveals:

Antisymmetric mode: Repulsive force (positive)
Symmetric mode: Attractive force (negative)
Both forces decay exponentially with separation

Mode Splitting

As waveguides approach:

Mode frequencies split (symmetric vs antisymmetric)
Coupling strength increases
Force magnitude increases

Key Findings

Mode	Force Sign	Physical Origin
Symmetric	Attractive	Fields add in gap
Antisymmetric	Repulsive	Fields cancel in gap

OptixLog Metrics

separation: Waveguide gap distance
mode_type: "symmetric" or "antisymmetric"
frequency: Mode frequency
flux: Optical power
normalized_force: F/L × (ac/P)

The force scales as ~exp(-s/δ) where δ is the evanescent decay length, roughly λ/(2π√(n²-1)).

Physical Insights

Maxwell Stress Tensor

The optical force is computed by integrating the Maxwell stress tensor:

T_ij = ε₀(E_i E_j - ½δ_ij E²) + μ₀(H_i H_j - ½δ_ij H²)
Force = ∮ T·n̂ dA

Gradient Force

For coupled waveguides, the dominant force is:

Proportional to field intensity gradient
Directed toward higher field regions
Sign depends on mode symmetry

Mode Decomposition - Mode analysis
Coupler - Directional coupler design
Straight Waveguide - Single waveguide modes

Parallel Waveguides Force