Far-field projection limitations with non-confined fields

Far-field projections are a great way to reduce computational cost by using the equivalence principle to calculate far fields based on near-field data. The general workflow is demonstrated in this notebook.

However, as discussed in our technical article, there are two main limitations of the far-field projection: the near fields must either be calculated on a closed surface with a uniform index, or the fields must decay at the borders of a 2D monitor.

This can be a challenging limitation for scatterers placed on a substrate, where a closed monitor will cross regions with different indices, and the fields are likely not decaying at the boundaries.

In this notebook, we will explore the limitations of this scenario and discuss caveats and workarounds.

import numpy as np
import pandas as pd
import tidy3d as td
from matplotlib import pyplot as plt
from tidy3d import web

Simulation Setup

We will analyze a simple case of a silicon box scatterer on a silicon oxide substrate.

First, we will create a function to build the base simulation for our analysis, using the simulation plane size as a parameter. To compute the far-field projection, we will place several 2D monitors of the FieldProjectionCartesianMonitor class above the scatterer, using different window_size values. This parameter defines the size of the transition region of the windowing function, ensuring that the recorded near fields decay to zero near the monitor edges.

We will also add one 3D FieldProjectionCartesianMonitor surrounding the scatterer, with the z-plane aligned to the substrate interface so it can be excluded from the far-field projection calculation.

For reference, we will add a FieldMonitor 10 µm away from the scatterer, and another one at the same position as the 2D far-field projection monitors.

Since the field projection distance is small, the far_field_approx argument is set to False.

freq0 = td.C_0 / 1.55
fwidth = 0.1 * freq0

def getSim(size):
    source = td.PlaneWave(
        name="source",
        center=[0, 0, -0.5],
        size=[1000, 1000, 0],
        source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
        direction="+",
    )

    field_monitor = td.FieldMonitor(
        center=[0, 0, 10],
        size=[size, size, 0],
        name="field_monitor",
        freqs=[freq0],
    )

    near_field_monitor = td.FieldMonitor(
        center=(0.0, 0.0, 2.1),
        size=(size, size, 0),
        name="near_field_monitor",
        freqs=[freq0],
        colocate=False,
    )

    far_field_monitors_2d = []
    for window_size in [0, 0.2, 0.4, 0.6, 0.8]:
        far_field_monitor_2d = td.FieldProjectionCartesianMonitor(
            center=(0.0, 0.0, 2.1),
            size=(size, size, 0),
            name=f"far_field_monitor_{window_size:.2f}",
            freqs=[freq0],
            normal_dir="+",
            far_field_approx=False,
            proj_axis=2,
            proj_distance=10,
            x=np.linspace(-field_monitor.size[0] / 2, field_monitor.size[0] / 2, 101),
            y=np.linspace(-field_monitor.size[1] / 2, field_monitor.size[1] / 2, 101),
            window_size=(window_size, window_size),
            custom_origin=(0, 0, 0),
        )
        far_field_monitors_2d.append(far_field_monitor_2d)

    far_field_monitor_3d = td.FieldProjectionCartesianMonitor(
        center=(0, 0, 1.5),
        size=(0.9 * size, 0.9 * size, 2.9),
        name="far_field_monitor_3D",
        freqs=[freq0],
        exclude_surfaces=["z-"],
        far_field_approx=False,
        proj_axis=2,
        proj_distance=10,
        x=np.linspace(-field_monitor.size[0] / 2, field_monitor.size[0] / 2, 101),
        y=np.linspace(-field_monitor.size[1] / 2, field_monitor.size[1] / 2, 101),
        custom_origin=(0, 0, 0),
    )

    SiO2 = td.Medium(
        name="medium_1",
        permittivity=1.45**2,
    )
    substrate = td.Structure(
        geometry=td.Box(center=[0, 0, -500], size=[2000, 2000, 1000]), name="substrate", medium=SiO2
    )

    Si = td.Medium(
        name="Si",
        permittivity=3.5**2,
    )
    scatterer = td.Structure(
        geometry=td.Box(center=[0, 0, 0.5], size=[1, 1, 1]),
        name="scatterer",
        background_medium=td.MultiPhysicsMedium(),
        medium=Si,
    )

    sim = td.Simulation(
        center=[0, 0, 5],
        size=[size, size, 12],
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=25),
        run_time=2e-12,
        sources=[source],
        monitors=[field_monitor, near_field_monitor, far_field_monitor_3d] + far_field_monitors_2d,
        structures=[substrate, scatterer],
    )

    return sim


sim_base = getSim(10)

Now we can take a look at the setup:

sim_base.plot(y=0)
plt.show()

sim_base.plot_3d()

sim_data = web.run(sim_base, task_name="far_field_1", folder_name="FFP")

04:15:38 UTC Created task 'far_field_1' with resource_id
             'fdve-c595a9ff-1437-4773-8c02-d8dbe5e22142' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-c595a9ff-143
             7-4773-8c02-d8dbe5e22142'.

             Task folder: 'default'.

04:15:40 UTC Estimated FlexCredit cost: 0.432. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the
             billed FlexCredit cost after a simulation run.

             status = success

04:15:43 UTC Loading results from simulation_data.hdf5

Now we can visualize side by side the field projection using the 2D monitor, with the window_size argument set to (0, 0), i.e., using the raw near-field data without filtering, and the far-field projection from the 3D monitor.

For reference, we will also plot the field monitor with the expected field profile, as well as the field monitor showing the raw near fields.

fig, ax = plt.subplots(2, 2, figsize=(8, 8))

sim_data["far_field_monitor_0.00"].power.plot(cmap="magma", ax=ax[0][0])
sim_data["far_field_monitor_3D"].power.plot(cmap="magma", ax=ax[0][1])
sim_data.plot_field("field_monitor", "E", "abs", ax=ax[1][0])
sim_data.plot_field("near_field_monitor", "E", "abs", ax=ax[1][1])

plt.show()

As we can see, the near fields do not decay at the boundaries, and consequently, the field projections do not capture the correct field profile.

Next, we can visualize the far-field projection for different values of window_size.

fig, ax = plt.subplots(2, 2, figsize=(12, 10))

sim_data["far_field_monitor_0.20"].power.plot(cmap="magma", ax=ax[0][0])
sim_data["far_field_monitor_0.40"].power.plot(cmap="magma", ax=ax[0][1])
sim_data["far_field_monitor_0.60"].power.plot(cmap="magma", ax=ax[1][0])
sim_data["far_field_monitor_0.80"].power.plot(cmap="magma", ax=ax[1][1])

ax[0][0].set_title("window_size = 0.20")
ax[0][1].set_title("window_size = 0.40")
ax[1][0].set_title("window_size = 0.60")
ax[1][1].set_title("window_size = 0.80")

for axes in ax.flatten():
    axes.set_aspect("equal")

plt.show()

Varying the Simulation Plane Size

Now, we will perform the same analysis, but increasing the simulation plane size.

sims = {}
for size in [20, 30, 40, 50]:
    sims[str(size)] = getSim(size)

batch = web.Batch(simulations=sims, folder_name="TestesFFPSize")

batch_data = batch.run(path_dir="FFP_Size")

04:15:49 UTC Started working on Batch containing 4 tasks.

04:15:53 UTC Maximum FlexCredit cost: 15.808 for the whole batch.

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

04:15:55 UTC Batch complete.

As we can see, for larger plane sizes, it is possible to obtain a correct far-field projection. For higher values of the window_size argument, the projection becomes closer to the true fields, at the cost of a smaller region of interest.

Therefore, in cases where a 3D field projection monitor is not feasible, and the fields do not sufficiently decay at the edges of a 2D projection monitor, a compromise between the simulation plane size and the usable region of interest is necessary.

fig, Ax = plt.subplots(4, 6, figsize=(24, 17))
import matplotlib as mpl

for i, size in enumerate(sims.keys()):
    sim_data2 = batch_data[str(size)]

    sim_data2["far_field_monitor_0.00"].power.plot(
        cmap="magma", ax=Ax[i][0], add_colorbar=False, vmin=0
    )
    sim_data2["far_field_monitor_0.20"].power.plot(
        cmap="magma", ax=Ax[i][1], add_colorbar=False, vmin=0
    )
    sim_data2["far_field_monitor_0.40"].power.plot(
        cmap="magma", ax=Ax[i][2], add_colorbar=False, vmin=0
    )
    sim_data2["far_field_monitor_0.60"].power.plot(
        cmap="magma", ax=Ax[i][3], add_colorbar=False, vmin=0
    )

    sim_data2["far_field_monitor_3D"].power.plot(
        cmap="magma", ax=Ax[i][4], add_colorbar=False, vmin=0
    )

    E = (
        np.abs(sim_data2["field_monitor"].Ex) ** 2
        + np.abs(sim_data2["field_monitor"].Ey) ** 2
        + np.abs(sim_data2["field_monitor"].Ez) ** 2
    )
    E.plot(cmap="magma", ax=Ax[i][5], add_colorbar=False, vmin=0)

    Ax[i][0].set_title(f"size = {size}")
    Ax[i][1].set_title("window_size = 0.20")
    Ax[i][2].set_title("window_size = 0.40")
    Ax[i][3].set_title("window_size = 0.60")
    Ax[i][4].set_title("3D projection monitor")
    Ax[i][5].set_title("field_monitor")

    Ax[i][1].axis("off")
    Ax[i][2].axis("off")
    Ax[i][3].axis("off")
    Ax[i][4].axis("off")
    Ax[i][5].axis("off")
for axes in Ax.flatten():
    axes.set_aspect("equal")


plt.show()