Calculating optical forces involves integrating the Maxwell stress tensor (MST) over a closed surface surrounding the region where the force is to be evaluated.

In the frequency domain, the time-averaged Maxwell stress tensor is defined as

$$ \langle T_{ij} \rangle = \frac{1}{2}\operatorname{Re}\!\left[\epsilon E_i E_j^* + \mu H_i H_j^* - \frac{1}{2}\left(\epsilon |\mathbf{E}|^2 + \mu |\mathbf{H}|^2\right)\delta_{ij}\right] $$

where $\mathbf{E}$ and $\mathbf{H}$ are the complex electric and magnetic fields, $\epsilon$ and $\mu$ are the permittivity and permeability of the medium, and $\delta_{ij}$ is the Kronecker delta.

The total time-averaged force acting on the enclosed object is obtained by integrating the MST over a closed surface $S$,

$$ \mathbf{F} = \oint_S \langle \mathbf{T} \rangle \cdot \hat{\mathbf{n}} \, dS, $$

where $\hat{\mathbf{n}}$ is the outward-pointing surface normal.

In Tidy3D, one can define eight FluxMonitors) surrounding the region of interest to form this closed surface and evaluate the surface integrals numerically.

In this notebook, we demonstrate the complete workflow:

1) Define auxiliary functions to generate the surface monitors around the region of interest, compute the MST, and evaluate the resulting force vector.

2) Validate the force calculation using a plane wave incident on a dielectric interface.

3) Compute the optical forces associated with the symmetric and antisymmetric modes of coupled waveguides, reproducing the results presented in
Michelle L. Povinelli, Marko Lončar, Mihai Ibanescu, Elizabeth J. Smythe, Steven G. Johnson, Federico Capasso, and John D. Joannopoulos, “Evanescent-wave bonding between optical waveguides”, Optics Letters 30(22), 2005. DOI: https://doi.org/10.1364/OL.30.003042

[1]:

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

Auxiliary Functions ¶

Creating Surface Monitors¶

Now, we will define a function to create eight FieldMonitor objects around the region of interest.

[2]:

def surface_mon(center, size, freqs):
    """Create paired field and permittivity monitors on the six faces of a box."""

    mons = []

    size = np.array(size, dtype=float)
    center = np.array(center, dtype=float)

    # Map axis labels to array indices for easy access.
    axis = {"x": 0, "y": 1, "z": 2}

    # Unit vectors along each Cartesian axis, used to translate monitor centers.
    normal = {
        "x": np.array([1, 0, 0]),
        "y": np.array([0, 1, 0]),
        "z": np.array([0, 0, 1]),
    }

    # Mask describing which dimensions should span each face monitor.
    mask = {
        "x": np.array([0, 1, 1]),
        "y": np.array([1, 0, 1]),
        "z": np.array([1, 1, 0]),
    }

    for c in ["x", "y", "z"]:
        i = axis[c]

        for d in [1, -1]:
            # Position the monitor on the +/- face along axis c.
            center_mon = tuple(center + d * normal[c] * size[i] / 2)
            size_mon = tuple(size * mask[c])

            name = f"{c}{'+' if d == 1 else '-'}"

            mons.append(
                td.FieldMonitor(
                    center=center_mon,
                    size=size_mon,
                    name=name,
                    freqs=freqs,
                    colocate=False,
                )
            )

            # Include matching permittivity monitor so material data is available
            # when evaluating the Maxwell stress tensor.
            mons.append(
                td.PermittivityMonitor(
                    center=center_mon,
                    size=size_mon,
                    name=f"e_{name}",
                    freqs=freqs,
                )
            )

    return mons

Maxwell Stress Tensor Calculation¶

Now, we will define a function to compute the Maxwell Stress Tensor from a Tidy3D SimulationData object, using the monitors defined by the surface_mon function.

[3]:

def stressTensor(sim_data, monitor_name, eps_mon):
    """Compute the Maxwell stress tensor on the plane defined by `monitor_name`, following the naming convention of the `surface_mon` function."""

    if monitor_name[0] == "x":
        i, j = ("y", "z")

        Ex = sim_data[monitor_name].Ex[:, 1:-1, 1:-1, :]
        Ey = sim_data[monitor_name].Ey[:, 1:-1, 1:-1, :]
        Ez = sim_data[monitor_name].Ez[:, 1:-1, 1:-1, :]

        Hx = sim_data[monitor_name].Hx[:, 1:-1, 1:-1, :]
        Hy = sim_data[monitor_name].Hy[:, 1:-1, 1:-1, :]
        Hz = sim_data[monitor_name].Hz[:, 1:-1, 1:-1, :]

        # Interpolate permittivity onto the field grid covering this face.
        eps_xx = sim_data[eps_mon].eps_xx.interp(y=Ex.coords[i], z=Ex.coords[j])
        eps_yy = sim_data[eps_mon].eps_yy.interp(y=Ey.coords[i], z=Ey.coords[j])
        eps_zz = sim_data[eps_mon].eps_zz.interp(y=Ez.coords[i], z=Ez.coords[j])

    elif monitor_name[0] == "y":
        i, j = ("x", "z")

        Ex = sim_data[monitor_name].Ex[1:-1, :, 1:-1, :]
        Ey = sim_data[monitor_name].Ey[1:-1, :, 1:-1, :]
        Ez = sim_data[monitor_name].Ez[1:-1, :, 1:-1, :]

        Hx = sim_data[monitor_name].Hx[1:-1, :, 1:-1, :]
        Hy = sim_data[monitor_name].Hy[1:-1, :, 1:-1, :]
        Hz = sim_data[monitor_name].Hz[1:-1, :, 1:-1, :]

        eps_xx = sim_data[eps_mon].eps_xx.interp(x=Ex.coords[i], z=Ex.coords[j])
        eps_yy = sim_data[eps_mon].eps_yy.interp(x=Ey.coords[i], z=Ey.coords[j])
        eps_zz = sim_data[eps_mon].eps_zz.interp(x=Ez.coords[i], z=Ez.coords[j])

    elif monitor_name[0] == "z":
        i, j = ("x", "y")

        Ex = sim_data[monitor_name].Ex[1:-1, 1:-1, :, :]
        Ey = sim_data[monitor_name].Ey[1:-1, 1:-1, :, :]
        Ez = sim_data[monitor_name].Ez[1:-1, 1:-1, :, :]

        Hx = sim_data[monitor_name].Hx[1:-1, 1:-1, :, :]
        Hy = sim_data[monitor_name].Hy[1:-1, 1:-1, :, :]
        Hz = sim_data[monitor_name].Hz[1:-1, 1:-1, :, :]

        eps_xx = sim_data[eps_mon].eps_xx.interp(x=Ex.coords[i], y=Ex.coords[j])
        eps_yy = sim_data[eps_mon].eps_yy.interp(x=Ey.coords[i], y=Ey.coords[j])
        eps_zz = sim_data[eps_mon].eps_zz.interp(x=Ez.coords[i], y=Ez.coords[j])

    coords = [Ex.x.squeeze(), Ex.y.squeeze(), Ex.z.squeeze()]

    # Convert xarray data to numpy arrays.
    Ey = Ey.values
    Ez = Ez.values
    Hx = Hx.values
    Hy = Hy.values
    Hz = Hz.values
    eps_xx = eps_xx.values
    eps_yy = eps_yy.values
    eps_zz = eps_zz.values

    # Displacement (D) and magnetic flux (B) densities.
    Dx = Ex * eps_xx * td.EPSILON_0
    Dy = Ey * eps_yy * td.EPSILON_0
    Dz = Ez * eps_zz * td.EPSILON_0

    Bx = Hx * td.MU_0
    By = Hy * td.MU_0
    Bz = Hz * td.MU_0

    # Scalar energy-density terms appearing in the stress tensor formula.
    E_mod = Ex * np.conj(Dx) + Ey * np.conj(Dy) + Ez * np.conj(Dz)
    H_mod = Hx * np.conj(Bx) + Hy * np.conj(By) + Hz * np.conj(Bz)

    # Maxwell stress tensor components prior to taking the real (time-averaged) part.
    Txx = (Ex * np.conj(Dx) - 0.5 * E_mod) + (Hx * np.conj(Bx) - 0.5 * H_mod)
    Txy = (Ex * np.conj(Dy)) + (Hx * np.conj(By))
    Txz = (Ex * np.conj(Dz)) + (Hx * np.conj(Bz))

    Tyy = (Ey * np.conj(Dy) - 0.5 * E_mod) + (Hy * np.conj(By) - 0.5 * H_mod)
    Tyx = (Ey * np.conj(Dx)) + (Hy * np.conj(Bx))
    Tyz = (Ey * np.conj(Dz)) + (Hy * np.conj(Bz))

    Tzz = (Ez * np.conj(Dz) - 0.5 * E_mod) + (Hz * np.conj(Bz) - 0.5 * H_mod)
    Tzx = (Ez * np.conj(Dx)) + (Hz * np.conj(Bx))
    Tzy = (Ez * np.conj(Dy)) + (Hz * np.conj(By))

    T = np.array([[Txx, Txy, Txz], [Tyx, Tyy, Tyz], [Tzx, Tzy, Tzz]])

    # Time-average via the real part and include the 1/2 prefactor.
    return 0.5 * T.real, coords

Force Calculation¶

Finally, we will define a function to integrate all monitors from a SimulationData object and compute the force vector.

[4]:

def integrateStressTensor(sim_data):
    """Integrate the Maxwell stress tensor over a closed box of monitors."""

    Fx_x = Fy_x = Fz_x = 0.0
    Fx_y = Fy_y = Fz_y = 0.0
    Fx_z = Fy_z = Fz_z = 0.0

    for c in ["x", "y", "z"]:
        for d in [1, -1]:
            name = c + ("+" if d == 1 else "-")
            T, coords = stressTensor(sim_data, name, "e_" + name)

            if c == "x":
                # Extract tensor components T_{ix} on the x± faces.
                Txx = T[0, 0].squeeze()
                Tyx = T[1, 0].squeeze()
                Tzx = T[2, 0].squeeze()

                iy = np.argsort(coords[1])
                iz = np.argsort(coords[2])

                y_sorted = np.asarray(coords[1])[iy]
                z_sorted = np.asarray(coords[2])[iz]

                Txx = Txx[iy][:, iz]
                Tyx = Tyx[iy][:, iz]
                Tzx = Tzx[iy][:, iz]

                Fx_x += d * np.trapz(np.trapz(Txx, y_sorted, axis=0), z_sorted, axis=0)
                Fy_x += d * np.trapz(np.trapz(Tyx, y_sorted, axis=0), z_sorted, axis=0)
                Fz_x += d * np.trapz(np.trapz(Tzx, y_sorted, axis=0), z_sorted, axis=0)

            elif c == "y":
                # Integrate T_{iy} over the y± faces.
                Txy = T[0, 1].squeeze()
                Tyy = T[1, 1].squeeze()
                Tzy = T[2, 1].squeeze()

                ix = np.argsort(coords[0])
                iz = np.argsort(coords[2])

                x_sorted = np.asarray(coords[0])[ix]
                z_sorted = np.asarray(coords[2])[iz]

                Txy = Txy[ix][:, iz]
                Tyy = Tyy[ix][:, iz]
                Tzy = Tzy[ix][:, iz]

                Fx_y += d * np.trapz(np.trapz(Txy, x_sorted, axis=0), z_sorted, axis=0)
                Fy_y += d * np.trapz(np.trapz(Tyy, x_sorted, axis=0), z_sorted, axis=0)
                Fz_y += d * np.trapz(np.trapz(Tzy, x_sorted, axis=0), z_sorted, axis=0)

            elif c == "z":
                # Integrate T_{iz} over the z± faces.
                Txz = T[0, 2].squeeze()
                Tyz = T[1, 2].squeeze()
                Tzz = T[2, 2].squeeze()

                ix = np.argsort(coords[0])
                iy = np.argsort(coords[1])

                x_sorted = np.asarray(coords[0])[ix]
                y_sorted = np.asarray(coords[1])[iy]

                Txz = Txz[ix][:, iy]
                Tyz = Tyz[ix][:, iy]
                Tzz = Tzz[ix][:, iy]

                Fx_z += d * np.trapz(np.trapz(Txz, x_sorted, axis=0), y_sorted, axis=0)
                Fy_z += d * np.trapz(np.trapz(Tyz, x_sorted, axis=0), y_sorted, axis=0)
                Fz_z += d * np.trapz(np.trapz(Tzz, x_sorted, axis=0), y_sorted, axis=0)

    # Sum contributions from all faces and apply the outward normal convention.
    Fx = Fx_x + Fx_y + Fx_z
    Fy = Fy_x + Fy_y + Fy_z
    Fz = Fz_x + Fz_y + Fz_z

    return Fx, Fy, Fz

Force at a Planar Interface With Oblique Incidence ¶

Consider an incident plane wave in medium 1 (index $n_1$) impinging on a planar interface with medium 2 (index $n_2$), with incidence angle $\theta_i$ measured from the surface normal $\hat{\mathbf{n}}$ pointing from medium 1 to medium 2. With no absorption, the net time-averaged normal force comes from the change in normal momentum flux of the incident, reflected, and transmitted waves.

If $I$ is the incident intensity in medium 1, and $R,T$ are the power reflectance/transmittance (so $R+T=1$), then the normal force on an illuminated area $A$ can be written as

$$ F_{normal} \;=\; \frac{I\,A}{c_0}\Big[\cos(\theta_i)\,(1+R)\;-\;n_2\,\cos(\theta_t)\,T\Big]. $$

To implement this in Tidy3D, we will create a function that generates the simulation object as a function of the permittivity of medium 2 and the incident angle. We will define a simple simulation with PlaneWave excitation and Bloch boundary conditions. To create the monitors, we will use the predefined

[5]:

wl = 1
freq0 = td.C_0 / 1
fwidth = 0.1 * freq0


def getSim(theta, permittivity=3.5**2, direction="z", min_steps_per_wvl=20):
    """Build a planar-scattering simulation for a slab and surrounding monitors."""

    plane_size = 1.3

    medium = td.Medium(permittivity=permittivity)

    # Frequencies to sample (single point for steady-state response).
    freqs = [freq0]

    # Convert incidence and polarization angles to radians.
    theta_rad = np.deg2rad(theta)

    # Select injection plane and Bloch axes depending on propagation direction.
    if direction == "z":
        center = (0, 0, 1)
        size = [td.inf, td.inf, 0]
        axis_1 = 0
        axis_2 = 1

    elif direction == "y":
        center = (0, 1, 0)
        size = [td.inf, 0, td.inf]
        axis_1 = 0
        axis_2 = 2

    elif direction == "x":
        center = (1, 0, 0)
        size = [0, td.inf, td.inf]
        axis_1 = 1
        axis_2 = 2

    # Define the incident plane wave illumination.
    planewave_0 = td.PlaneWave(
        size=size,
        center=center,
        source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
        direction="-",
        angle_theta=theta_rad,
    )

    # Bloch boundaries.
    bloch_1 = td.Boundary.bloch_from_source(
        source=planewave_0, domain_size=plane_size, axis=axis_1, medium=td.Medium(permittivity=1)
    )

    bloch_2 = td.Boundary.bloch_from_source(
        source=planewave_0, domain_size=plane_size, axis=axis_2, medium=td.Medium(permittivity=1)
    )

    # Create closed surface monitors used later for force calculations.
    surface_mons = surface_mon(center=(0, 0, 0), size=(1, 1, 1), freqs=freqs)

    if direction == "x":
        bspec = td.BoundarySpec(y=bloch_1, z=bloch_2, x=td.Boundary.absorber(num_layers=80))
        center_box = [-55.5, 500, 500]
        size_box = [111, 1222, 1222]
        sim_size = [4, plane_size, plane_size]
    elif direction == "y":
        bspec = td.BoundarySpec(x=bloch_1, z=bloch_2, y=td.Boundary.absorber(num_layers=80))
        center_box = [500, -55.5, 500]
        size_box = [1222, 111, 1222]
        sim_size = [plane_size, 4, plane_size]
    elif direction == "z":
        bspec = td.BoundarySpec(x=bloch_1, y=bloch_2, z=td.Boundary.absorber(num_layers=80))
        center_box = [500, 500, -55.5]
        size_box = [1222, 1222, 111]
        sim_size = [plane_size, plane_size, 4]

    # Slab structure occupying the high-index region.
    box_0 = td.Structure(
        geometry=td.Box(center=center_box, size=size_box), name="box_0", medium=medium
    )

    # Assemble the full simulation with sources, monitors, and structure.
    sim = td.Simulation(
        size=sim_size,
        boundary_spec=bspec,
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl),
        run_time=5e-12,
        sources=[planewave_0],
        monitors=surface_mons,
        structures=[box_0],
    )

    return sim


sim = getSim(45)
_ = sim.plot(y=0)

plt.show()

No description has been provided for this image

Now we will vary the incidence angle and compare the results with the analytical solution. To test all three directions, we will create three batches, corresponding to plane-wave incidence along the x, y, and z directions.

[6]:

sims_x = {}
sims_y = {}
sims_z = {}
thetas = [0, 10, 20, 30, 40, 50]
for theta in thetas:
    sims_x[f"{theta}"] = getSim(theta, direction="x")
    sims_y[f"{theta}"] = getSim(theta, direction="y")
    sims_z[f"{theta}"] = getSim(theta, direction="z")

batch_x = web.Batch(simulations=sims_x, folder_name="OpticalForces_x")
batch_data_thetas_x = batch_x.run(path_dir="batchThetas_x")

batch_y = web.Batch(simulations=sims_y, folder_name="OpticalForces_y")
batch_data_thetas_y = batch_y.run(path_dir="batchThetas_y")

batch_z = web.Batch(simulations=sims_z, folder_name="OpticalForces_z")
batch_data_thetas_z = batch_z.run(path_dir="batchThetas_z")

Output()

15:16:27 -03 Started working on Batch containing 6 tasks.

15:16:37 -03 Maximum FlexCredit cost: 1.961 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.

Output()

15:16:52 -03 Batch complete.

Output()

15:17:00 -03 Started working on Batch containing 6 tasks.

15:17:09 -03 Maximum FlexCredit cost: 1.961 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.

Output()

15:17:25 -03 Batch complete.

Output()

15:17:32 -03 Started working on Batch containing 6 tasks.

15:17:42 -03 Maximum FlexCredit cost: 1.961 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.

Output()

15:17:58 -03 Batch complete.

[7]:

axes = ["x", "y", "z"]
components = ["Fx", "Fy", "Fz"]

# Pre-allocate containers for the force components keyed by face axis.
forces = {axis: {comp: [] for comp in components} for axis in axes}

for theta_key in batch_data_thetas_x.keys():
    print(theta_key)

    # Evaluate the stress tensor on each pair of monitors.
    evaluations = {
        "x": integrateStressTensor(batch_data_thetas_x[theta_key]),
        "y": integrateStressTensor(batch_data_thetas_y[theta_key]),
        "z": integrateStressTensor(batch_data_thetas_z[theta_key]),
    }

    # Append the results to the appropriate component lists.
    for axis, (fx_val, fy_val, fz_val) in evaluations.items():
        for comp, value in zip(components, (fx_val, fy_val, fz_val)):
            forces[axis][comp].append(np.array(value))

# Unpack back into the original variable names for downstream cells.
FxT_x, FyT_x, FzT_x = forces["x"]["Fx"], forces["x"]["Fy"], forces["x"]["Fz"]
FxT_y, FyT_y, FzT_y = forces["y"]["Fx"], forces["y"]["Fy"], forces["y"]["Fz"]
FxT_z, FyT_z, FzT_z = forces["z"]["Fx"], forces["z"]["Fy"], forces["z"]["Fz"]

15:18:02 -03 Loading simulation from                                            
             batchThetas_x/fdve-99e2b64c-257a-44d5-b3c5-2372b1f68314.hdf5

15:18:03 -03 Loading simulation from                                            
             batchThetas_y/fdve-4c073505-6e95-4b5b-875e-fc5148f3be13.hdf5

15:18:04 -03 Loading simulation from                                            
             batchThetas_z/fdve-0126fde7-5d90-43c4-b788-cdf6d429840c.hdf5

15:18:05 -03 Loading simulation from                                            
             batchThetas_x/fdve-e6bc84cd-be55-4c7d-9357-61ad393681c8.hdf5

15:18:06 -03 Loading simulation from                                            
             batchThetas_y/fdve-e67995e5-a9e5-4231-8854-b7833e9d63f1.hdf5

15:18:07 -03 Loading simulation from                                            
             batchThetas_z/fdve-cc605a19-cf3f-4ca5-8155-613ccf4e0f12.hdf5

15:18:09 -03 Loading simulation from                                            
             batchThetas_x/fdve-4fe84d05-0779-466e-a925-e800f6c152f4.hdf5

15:18:10 -03 Loading simulation from                                            
             batchThetas_y/fdve-e1d1d5da-2bc6-4609-b02b-4b80640c114f.hdf5

15:18:11 -03 Loading simulation from                                            
             batchThetas_z/fdve-d5fc6a02-73e8-4d64-a906-fbf5b50cd3b9.hdf5

15:18:12 -03 Loading simulation from                                            
             batchThetas_x/fdve-de429a97-17d4-4019-a319-1f0b0a6d85c2.hdf5

15:18:13 -03 Loading simulation from                                            
             batchThetas_y/fdve-e4acba72-c067-49af-8116-9d81364eac99.hdf5

15:18:14 -03 Loading simulation from                                            
             batchThetas_z/fdve-fe6f5d3f-7f2c-4f88-adb2-7b0505df29d4.hdf5

15:18:15 -03 Loading simulation from                                            
             batchThetas_x/fdve-3153c32d-9779-4c34-9e3a-de8a0ca505e8.hdf5

15:18:16 -03 Loading simulation from                                            
             batchThetas_y/fdve-68719354-7b94-46cf-ae26-320807537525.hdf5

15:18:18 -03 Loading simulation from                                            
             batchThetas_z/fdve-eb1d756f-1af8-49be-b0c6-792deab76fe4.hdf5

15:18:19 -03 Loading simulation from                                            
             batchThetas_x/fdve-741932f0-18e7-42a9-aa9a-71a683edc5ca.hdf5

15:18:20 -03 Loading simulation from                                            
             batchThetas_y/fdve-7fb06ad8-d125-4aa4-b515-52d6d52ac840.hdf5

15:18:21 -03 Loading simulation from                                            
             batchThetas_z/fdve-1e217314-839f-41a0-9595-3e1f46f9580a.hdf5

Now we will define the analytical function. Note that the PlaneWave source is normalized to inject 1 W of power. Therefore, the total power crossing the monitor area must be corrected by the simulation plane size. More information on source normalization can be found here.

[8]:

def interface_force_eps(eps, I=1, A=1, theta_i=0):
    import numpy as np

    theta_i = np.deg2rad(theta_i)
    n1 = 1.0
    n2 = np.sqrt(eps)

    theta_t = np.arcsin(np.sin(theta_i) / n2)

    # Fresnel power coefficients (TE and TM)
    R = (
        (n2 * np.cos(theta_i) - n1 * np.cos(theta_t))
        / (n2 * np.cos(theta_i) + n1 * np.cos(theta_t))
    ) ** 2

    T = 1 - R  # lossless

    Fz = (I * A / td.C_0) * (np.cos(theta_i) * (1 + R) - n2 * np.cos(theta_t) * T)

    return np.abs(Fz)


# Adjusting the incident intensity to match the simulation results
I0 = (1 / 1.3) ** 2
Fz_a_theta = interface_force_eps(eps=3.5**2, theta_i=np.array(thetas), I=I0)

Now we can visualize the results. The force normal to the incident plane follows the analytical expression for all three directions, and the other components are close to zero, as expected.

[9]:

fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharey=False)
axes[0].plot(thetas, np.array([np.sum(i) for i in FxT_x]), label="Fx")
axes[0].plot(thetas, np.array([np.sum(i) for i in FyT_x]), label="Fy")
axes[0].plot(thetas, np.array([np.sum(i) for i in FzT_x]), label="Fz")
axes[0].plot(thetas, np.array(Fz_a_theta), "o", label="Analytic")
axes[0].set_title("x")
axes[0].legend()

axes[1].plot(thetas, np.array([np.sum(i) for i in FxT_y]), label="Fx")
axes[1].plot(thetas, np.array([np.sum(i) for i in FyT_y]), label="Fy")
axes[1].plot(thetas, np.array([np.sum(i) for i in FzT_y]), label="Fz")
axes[1].plot(thetas, np.array(Fz_a_theta), "o", label="Analytic")
axes[1].set_title("y")
axes[1].legend()

axes[2].plot(thetas, np.array([np.sum(i) for i in FxT_z]), label="Fx")
axes[2].plot(thetas, np.array([np.sum(i) for i in FyT_z]), label="Fy")
axes[2].plot(thetas, np.array([np.sum(i) for i in FzT_z]), label="Fz")
axes[2].plot(thetas, np.array(Fz_a_theta), "o", label="Analytic")
axes[2].set_title("z")
axes[2].legend()


for ax in axes:
    ax.set_xlabel("Angle (deg)")
    ax.set_ylabel("Force (N)")


plt.show()

Force in Coupled Waveguides ¶

Next, we calculate the attraction and repulsion forces in a system of two coupled waveguides, as described in the reference paper.

The system consists of two square waveguides with size $a = 0.31$ µm and varying separation $d$. Two types of modes can be supported: symmetric modes (same $E_z$ field sign in both waveguides) and antisymmetric modes (opposite signs).

We will define a function to create the setup based on the separation $d$, add a surface monitor around one of the waveguides, and calculate the net force in the direction perpendicular to the waveguides.

[10]:

from typing import List, Optional


def build_double_waveguide_sim(
    a: float = 0.310,
    d: float = 0.25,
    wavelength: float = 1.55,
    wg_medium: Optional[td.Medium] = td.Medium(permittivity=3.45**2),
    x_span: float = 6.0,
    buffer_y: float = 2,
    buffer_z: float = 2,
    source_offset: float = 0.2,
    mode_index=1,
) -> td.Simulation:
    """Create a Tidy3D simulation containing two identical square waveguides"""

    if a <= 0.0:
        raise ValueError("Waveguide size 'a' must be positive.")
    if d < 0.0:
        raise ValueError("Waveguide separation 'd' must be non-negative.")
    if source_offset >= x_span / 2:
        raise ValueError("Source offset must be smaller than half of the simulation span.")

    # Distance from simulation centre to each waveguide centre along y.
    half_center_sep = 0.5 * (a + d)
    waveguide_centers: List[float] = [-half_center_sep, half_center_sep]

    # Simulation extents along each dimension (µm).
    y_half_span = a + d / 2 + buffer_y
    z_half_span = a / 2 + buffer_z
    sim_size = (x_span, 2 * y_half_span, 2 * z_half_span)

    structures = [
        td.Structure(
            geometry=td.Box(center=(0.0, yc, 0.0), size=(td.inf, a, a)),
            medium=wg_medium,
        )
        for yc in waveguide_centers
    ]

    freq0 = td.C_0 / wavelength
    pulse = td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10)
    source_x = -0.5 * x_span + source_offset

    sources = [
        td.ModeSource(
            center=(source_x, 0, 0.0),
            size=(0.0, td.inf, td.inf),
            source_time=pulse,
            direction="+",
            mode_spec=td.ModeSpec(num_modes=4),
            mode_index=mode_index,
        )
    ]

    sim = td.Simulation(
        size=sim_size,
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=50),
        structures=structures,
        sources=sources,
        run_time=5e-12,
        monitors=surface_mon(
            center=(0, (d / 2 + a / 2 + 0.5), 0), size=(1, a + d / 2 + 1, 1), freqs=[freq0]
        ),
        symmetry=(0, 0, -1),
    )

    ms = ModeSolver(
        simulation=sim,
        plane=sources[0].bounding_box,
        freqs=[sources[0].source_time.freq0],
        mode_spec=sources[0].mode_spec,
    )

    return sim, ms


sim, ms = build_double_waveguide_sim(d=0.1)

First, we will run a mode solver simulation to identify the effective indices of the symmetric and antisymmetric modes.

[11]:

mode_data = web.run(ms, "mode")

15:18:22 -03 Created task 'mode' with resource_id                               
             'mo-813a64c5-c8e7-44ee-820b-72b8bd09eade' and task_type            
             'MODE_SOLVER'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=mo-813a64c5-c8e7-
             44ee-820b-72b8bd09eade'.

             Task folder: 'default'.

Output()

15:18:25 -03 Estimated FlexCredit cost: 0.006. Minimum cost depends on task     
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

15:18:28 -03 status = success

Output()

15:18:32 -03 Loading simulation from simulation_data.hdf5

[12]:

fig, axes = plt.subplots(1, 2, figsize=(12, 4), constrained_layout=True)

for ax, (mode_index, title) in zip(
    axes,
    [
        (0, "E_z (symmetric, mode 0)"),
        (1, "E_z (antisymmetric, mode 1)"),
    ],
):
    ms.plot_field("Ez", "real", mode_index=mode_index, ax=ax)
    ax.set_title(title)

# 0 → symmetric fundamental mode
# 1 → antisymmetric higher-order mode

plt.show()

Next, we will sweep the values of $d$ and create two batch simulations for each symmetry.

[13]:

d_values = np.linspace(0.03, 0.31, 21)
sims_1 = {i: build_double_waveguide_sim(d=i, mode_index=0)[0] for i in d_values}
batch_1 = web.Batch(simulations=sims_1, folder_name="double_waveguide_1")
batch_data_1 = batch_1.run(path_dir="double_waveguide_1")

sims_2 = {i: build_double_waveguide_sim(d=i, mode_index=1)[0] for i in d_values}
batch_2 = web.Batch(simulations=sims_2, folder_name="double_waveguide_2")
batch_data_2 = batch_2.run(path_dir="double_waveguide_2")

Output()

15:18:46 -03 Started working on Batch containing 21 tasks.

15:19:19 -03 Maximum FlexCredit cost: 30.935 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.

Output()

15:20:12 -03 Batch complete.

Output()

15:20:38 -03 Started working on Batch containing 21 tasks.

15:21:12 -03 Maximum FlexCredit cost: 30.935 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.

Output()

15:22:05 -03 Batch complete.

Finally, we calculate the forces for each simulation and plot the y component.

As we can see, the results closely match Fig. 1(b) of the reference paper.

[14]:

F_1 = [integrateStressTensor(batch_data_1[i]) for i in batch_data_1.keys()]
F_2 = [integrateStressTensor(batch_data_2[i]) for i in batch_data_2.keys()]

15:22:17 -03 Loading simulation from                                            
             double_waveguide_1/fdve-429e6f98-0488-4bad-8734-6b2077937952.hdf5

15:22:19 -03 Loading simulation from                                            
             double_waveguide_1/fdve-dad279bb-f717-42cd-a7be-bb1f9ef67f24.hdf5

15:22:21 -03 Loading simulation from                                            
             double_waveguide_1/fdve-7c3f1732-90ba-47f5-b581-42bd4f551efd.hdf5

15:22:23 -03 Loading simulation from                                            
             double_waveguide_1/fdve-8e4321b1-be10-49c1-b573-a0d18e4df520.hdf5

15:22:25 -03 Loading simulation from                                            
             double_waveguide_1/fdve-26ca36e4-859a-4717-aa74-bcaa5b933dbc.hdf5

15:22:27 -03 Loading simulation from                                            
             double_waveguide_1/fdve-78cfadb3-0d5e-43b9-ad73-8abe662781f2.hdf5

15:22:29 -03 Loading simulation from                                            
             double_waveguide_1/fdve-e9f01ddf-b99e-48ea-a92e-a19072296fae.hdf5

15:22:31 -03 Loading simulation from                                            
             double_waveguide_1/fdve-a8aa28bf-4b73-495b-901f-f809a65befc1.hdf5

15:22:33 -03 Loading simulation from                                            
             double_waveguide_1/fdve-6fcfe4aa-f4a5-4352-ba2c-4642b4f3e1e6.hdf5

15:22:35 -03 Loading simulation from                                            
             double_waveguide_1/fdve-c18c14e6-b37f-4b4e-9f2b-e79bea743d76.hdf5

15:22:37 -03 Loading simulation from                                            
             double_waveguide_1/fdve-27e00ab1-badd-4936-93cd-03320bddb275.hdf5

15:22:39 -03 Loading simulation from                                            
             double_waveguide_1/fdve-16c4e51e-eb70-4f72-9806-d2d427401c95.hdf5

15:22:41 -03 Loading simulation from                                            
             double_waveguide_1/fdve-77cd041f-32c9-495d-8542-fbed16fd4558.hdf5

15:22:42 -03 Loading simulation from                                            
             double_waveguide_1/fdve-2b00e780-a1b3-4b2c-b2e3-9e6366df914c.hdf5

15:22:44 -03 Loading simulation from                                            
             double_waveguide_1/fdve-2e5ebaf8-ebd5-4cb0-986a-43c84aca5a5b.hdf5

15:22:46 -03 Loading simulation from                                            
             double_waveguide_1/fdve-dca03e55-7c04-4105-a9b4-a6c914534d97.hdf5

15:22:48 -03 Loading simulation from                                            
             double_waveguide_1/fdve-6adf16d7-d260-4e43-a9af-c937b6bd393b.hdf5

15:22:50 -03 Loading simulation from                                            
             double_waveguide_1/fdve-a5277208-7c3e-47ca-9ef8-d57047be2635.hdf5

15:22:52 -03 Loading simulation from                                            
             double_waveguide_1/fdve-345414a6-f1d3-4ace-a109-2b89db77c1bb.hdf5

15:22:54 -03 Loading simulation from                                            
             double_waveguide_1/fdve-f31c4c82-3db2-43a7-bdd1-099718d90e14.hdf5

15:22:56 -03 Loading simulation from                                            
             double_waveguide_1/fdve-302b98b2-daf6-432d-9445-ae18bb60f6ba.hdf5

15:22:58 -03 Loading simulation from                                            
             double_waveguide_2/fdve-c440d48a-5b8d-4f07-a4ab-5e82838efbf0.hdf5

15:23:00 -03 Loading simulation from                                            
             double_waveguide_2/fdve-a2630227-eba8-48a7-abce-10c9b00c6e4c.hdf5

15:23:02 -03 Loading simulation from                                            
             double_waveguide_2/fdve-8d7c75e1-a14e-4ad7-8db6-82170f32ea74.hdf5

15:23:04 -03 Loading simulation from                                            
             double_waveguide_2/fdve-1e993cf9-d8a1-4caf-bc0f-1c48ced4c3dc.hdf5

15:23:06 -03 Loading simulation from                                            
             double_waveguide_2/fdve-1cab0798-d5fa-47fd-83c2-9d4891c1042a.hdf5

15:23:08 -03 Loading simulation from                                            
             double_waveguide_2/fdve-1cae816f-4edd-48a7-b6f0-2f0894d49af7.hdf5

15:23:10 -03 Loading simulation from                                            
             double_waveguide_2/fdve-8eea836a-c3b7-4d31-8727-a7177d682210.hdf5

15:23:12 -03 Loading simulation from                                            
             double_waveguide_2/fdve-e00e118b-23bc-4e71-b56a-4e5e2f9f4650.hdf5

15:23:14 -03 Loading simulation from                                            
             double_waveguide_2/fdve-9909c949-8fa6-4000-a5e5-d7ea94efd6cf.hdf5

15:23:16 -03 Loading simulation from                                            
             double_waveguide_2/fdve-08b1355b-0870-4ea0-bcc4-b49228ac05ef.hdf5

15:23:18 -03 Loading simulation from                                            
             double_waveguide_2/fdve-2bd6b712-43b1-4a10-91e0-9a77d1bc684d.hdf5

15:23:20 -03 Loading simulation from                                            
             double_waveguide_2/fdve-89882346-844f-49fd-9f08-bc621cd9a462.hdf5

15:23:22 -03 Loading simulation from                                            
             double_waveguide_2/fdve-11f993ce-866e-45d4-9c5d-44d284567ade.hdf5

15:23:24 -03 Loading simulation from                                            
             double_waveguide_2/fdve-5e7875af-f5c5-4456-adaa-dede06ee059f.hdf5

15:23:26 -03 Loading simulation from                                            
             double_waveguide_2/fdve-98de4e3d-10fc-4336-86ad-a2062c908eb3.hdf5

15:23:28 -03 Loading simulation from                                            
             double_waveguide_2/fdve-37eb5f7b-4151-4de6-ad14-b4efe5c1059d.hdf5

15:23:30 -03 Loading simulation from                                            
             double_waveguide_2/fdve-1dcf955a-6b14-4500-8269-2fdf01d1abba.hdf5

15:23:31 -03 Loading simulation from                                            
             double_waveguide_2/fdve-db0c3add-71fd-481d-b2eb-10e3e542bd5c.hdf5

15:23:33 -03 Loading simulation from                                            
             double_waveguide_2/fdve-18de35cb-c1c0-4630-89ea-809d14c36d45.hdf5

15:23:35 -03 Loading simulation from                                            
             double_waveguide_2/fdve-68aa8314-8865-493a-8faa-c1d51e3a2ac9.hdf5

15:23:37 -03 Loading simulation from                                            
             double_waveguide_2/fdve-4fd09745-b0e4-4d95-a31b-5850af122d8b.hdf5

[15]:

fig, ax = plt.subplots()
factor = 0.31 * td.C_0 / 1  # Source is normalized to 1W
ax.plot(d_values / 0.31, np.array([i[1] for i in F_2]) * factor, label="Antisymmetric")
ax.plot(d_values / 0.31, np.array([i[1] for i in F_1]) * factor, label="Symmetric")
ax.set_xlabel("d/a")
ax.set_ylabel("Force (N/L)(ac/P)")

ax.legend()
plt.show()