How to model a multilevel blazed diffraction grating in Tidy3D FDTD

Multilevel blazed diffraction grating¶

In this example, we compute the grating efficiency of a multilevel diffraction grating whose design is inspired by M. Oliva, T. Harzendorf, D. Michaelis, U. D. Zeitner, and A. Tünnermann, "Multilevel blazed gratings in resonance domain: an alternative to the classical fabrication approach," Opt. Express 19, 14735-14745 (2011), DOI: 10.1364/OE.19.014735.

Tidy3D uses a near field to far field transformation specialized to periodic structures to compute the grating efficiency and its accuracy is verified through a comparison with the semi-analytical rigorous coupled wave analysis (RCWA) approach using the open-source library grcwa.

If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.

In [1]:

# basic python imports
import numpy as np
import matplotlib.pyplot as plt

# Tidy3D import
import tidy3d as td
from tidy3d import web

Normal incidence¶

We will first analyze the grating under normal incidence, as also studied in the paper. In this case, we can use periodic boundary conditions in both tangential directions.

Geometry setup¶

First, the structure and simulation geometry are defined. The structure includes a dielectric substrate with two dielectric patterned layers.

In [2]:

# Grating parameters (all lengths in um)
index = 1.46847
period = 1.866
width_layer1 = 0.519
width_layer2 = 1.202
height_layer1 = 0.333
height_layer2 = 0.281

# free space central wavelength
wavelength = 0.416

# Simulation domain geometry
space_above = wavelength * 3
height_substrate = wavelength * 3
space_below = wavelength * 3

# Define a buffer to make sure objects extend past the simulation boundary
buffer = 0.1

# Simulation domain along x and z
length_x = period
length_z = space_below + height_substrate + height_layer1 + height_layer2 + space_above

# Define the medium
grating_medium = td.Medium(permittivity=index**2)

# Create the substrate
substrate = td.Structure(
    geometry=td.Box(
        center=[0, 0, -length_z / 2 + height_substrate / 2 + space_below],
        size=[td.inf, td.inf, height_substrate],
    ),
    medium=grating_medium,
)

# Level 1 grating
center_L1 = [
    -buffer / 2 - length_x / 2 + width_layer1 / 2 + width_layer2 / 2,
    0,
    -length_z / 2 + space_below + height_substrate + height_layer2 / 2,
]
size_L1 = [width_layer1 + width_layer2 + buffer, td.inf, height_layer2]
grating_L1 = td.Structure(
    geometry=td.Box(center=center_L1, size=size_L1),
    medium=grating_medium,
)

# Level 2 grating
center_L2 = [
    -buffer / 2 - length_x / 2 + width_layer1 / 2,
    0,
    -length_z / 2 + space_below + height_substrate + height_layer2 + height_layer1 / 2,
]
size_L2 = [width_layer1 + buffer, td.inf, height_layer1]
grating_L2 = td.Structure(
    geometry=td.Box(center=center_L2, size=size_L2),
    medium=grating_medium,
)

# Collect all structures
structures = [substrate, grating_L1, grating_L2]

Source setup¶

Next, define the source plane wave impinging from above the grating at normal incidence.

In [3]:

# Central frequency in Hz
f0 = td.C_0 / wavelength

# Bandwidth
fwidth = f0 / 40.0

# Run time
run_time = 100 / fwidth

# Time dependence of source
gaussian = td.GaussianPulse(freq0=f0, fwidth=fwidth)

# Source
src_z = length_z / 2 - space_above / 2
angle_theta = np.pi / 10
source = td.PlaneWave(
    center=(0, 0, src_z),
    size=(td.inf, td.inf, 0),
    source_time=gaussian,
    direction="-",
    pol_angle=0,
    angle_theta=0,
    angle_phi=0,
)

Monitor setup¶

Here, we'll set up a field monitor to record and plot the frequency domain fields at a plane in the xz cross-section. We'll also set up two DiffractionMonitor, one for reflection, and the other for transmission.

In [4]:

# Fields
monitor_xz = td.FieldMonitor(
    center=[0, 0, 0], size=[td.inf, 0, td.inf], freqs=[f0], name="field_xz"
)

# The allowed orders will be computed automatically and returned as part of the results
monitor_r = td.DiffractionMonitor(
    center=[0.0, 0.0, length_z / 2 - wavelength],
    size=[td.inf, td.inf, 0],
    freqs=[f0],
    name="reflection",
    normal_dir="+",
)

monitor_t = td.DiffractionMonitor(
    center=[0.0, 0.0, -length_z / 2 + wavelength],
    size=[td.inf, td.inf, 0],
    freqs=[f0],
    name="transmission",
    normal_dir="-",
)

monitors = [monitor_xz, monitor_r, monitor_t]

Set up boundary conditions and initialize simualtion¶

For normal incidence, we can use periodic boundary conditions along the x and y directions. More generally, we need to use Bloch boundary conditions as will be illustrated below. We can also use Bloch boundaries with zero Bloch vector for normal incidence, but a simulation with Bloch boundaries uses complex fields and is twice more computationally expensive than a simulation with periodic boundaries, while the results for bloch_vec = 0 are equivalent.

Along z, a perfectly matched layer (PML) is applied to mimic an infinite domain. Because the diffraction grating introduces waves propagating at various angles, including steep angles with respect to the PML boundary, we use more layers than the default value to minimize spurious reflections at the PMLs.

In [5]:

# Simulation size
length_y = 0  # grating is translationally invariant along y
sim_size = (length_x, length_y, length_z)

# Resolution
min_steps_per_wvl = 60

# Boundaries
num_pml_layers = 40
boundary_spec = td.BoundarySpec(
    x=td.Boundary.periodic(),
    y=td.Boundary.periodic(),
    z=td.Boundary(
        minus=td.PML(num_layers=num_pml_layers), plus=td.PML(num_layers=num_pml_layers)
    ),
)

# Simulation
simulation = td.Simulation(
    size=sim_size,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl),
    structures=structures,
    sources=[source],
    monitors=monitors,
    run_time=run_time,
    boundary_spec=boundary_spec,
)

fig, ax = plt.subplots(1, 1, figsize=(5, 8))
simulation.plot(y=0, ax=ax)
plt.show()

Run simulation¶

In [6]:

sim_data = web.run(
    simulation,
    task_name="GratingEfficiency",
    path="data/GratingEfficiency.hdf5",
    verbose=True,
)

[16:27:36] Created task 'GratingEfficiency' with task_id                        webapi.py:139
           'fdve-b947ba35-c6ec-460f-a671-739525c62971v1'.

           View task using web UI at 'https://tidy3d.simulation.cloud/workbench webapi.py:141
           ?taskId=fdve-b947ba35-c6ec-460f-a671-739525c62971v1'.

Output()

[16:27:42] status = queued                                                      webapi.py:271

Output()

[16:27:49] status = preprocess                                                  webapi.py:265

[16:27:55] Maximum FlexCredit cost: 0.061. Use 'web.real_cost(task_id)' to get  webapi.py:288
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:292

           running solver                                                       webapi.py:302

Output()

[16:28:25] early shutoff detected, exiting.                                     webapi.py:316

           status = postprocess                                                 webapi.py:333

Output()

[16:28:39] status = success                                                     webapi.py:340

Output()

[16:28:51] loading SimulationData from data/GratingEfficiency.hdf5              webapi.py:512

[16:28:51] WARNING: updating Simulation from 2.1 to 2.2                             log.py:50

Diffraction data¶

Now we can extract the diffracted power from the output data structures and verfify that the sum across all reflection and transmission orders is close to 1. We can also access the diffraction angles and the complex power amplitudes for each order and polarization.

In [7]:

data_r = sim_data["reflection"]
data_t = sim_data["transmission"]
total_power = data_r.power.sum() + data_t.power.sum()
theta = data_t.angles[0].squeeze(drop=True).values

print(f"Total power: {total_power.values}")
print("Theta (degrees):", ", ".join(f"{t:.2f}" for t in theta * 180 / np.pi))
print(f"Amplitude data: \n{data_t.amps}")

Total power: <xarray.DataArray ()>
array(0.9980561)
Theta (degrees): 63.09, 41.98, 26.48, 12.88, 0.00, 12.88, 26.48, 41.98, 63.09
Amplitude data: 
<xarray.DataArray (orders_x: 9, orders_y: 1, f: 1, polarization: 2)>
array([[[[ 5.11645475e-18-1.13903179e-17j,
           3.12478968e-02-6.99604031e-02j]]],


       [[[-1.25584969e-17+7.76529007e-18j,
          -9.82699030e-02+6.07843732e-02j]]],


       [[[ 3.74715466e-17-2.01430094e-17j,
           3.04128878e-01-1.63471254e-01j]]],


       [[[ 3.59066936e-17+5.03400521e-17j,
           2.93110272e-01+4.10932509e-01j]]],


       [[[ 0.00000000e+00+0.00000000e+00j,
          -9.29520360e-03+5.36648245e-01j]]],


       [[[ 0.00000000e+00+0.00000000e+00j,
          -2.23863585e-01+2.80712710e-01j]]],


       [[[ 0.00000000e+00+0.00000000e+00j,
           1.24167995e-01-2.42265122e-01j]]],


       [[[ 0.00000000e+00+0.00000000e+00j,
           1.59876588e-01+9.83122650e-02j]]],


       [[[ 0.00000000e+00+0.00000000e+00j,
           1.26630822e-01+8.93053699e-02j]]]])
Coordinates:
  * orders_x      (orders_x) int64 -4 -3 -2 -1 0 1 2 3 4
  * orders_y      (orders_y) int64 0
  * f             (f) float64 7.207e+14
  * polarization  (polarization) <U1 's' 'p'

Reference results¶

To validate the accuracy of the results from Tidy3D, We will now compute the grating efficiency using the grcwa package. Be sure to install it in your Python environment first: pip install grcwa.

In [8]:

import grcwa

# Define lattice constants - size of the domain
# grcwa requires a finite non-zero size along each dimension
size_y = 1e-3
L1 = [sim_size[0], 0]
L2 = [0, size_y]

# Set truncation order
nG = 300

# Set up RCWA object
freq = f0 / td.C_0  # grcwa uses normalized coordinates where the speed of light is 1
obj = grcwa.obj(nG, L1, L2, freq, 0, 0, verbose=0)

# Set up the geometry (the layers are ordered top to bottom in grcwa)
num_patterned_layers = 3
thick_top = space_above
thick_layers = [height_layer1, height_layer2, height_substrate]
thick_bot = space_below

# Discretization points along x and y
num_x = 300
num_y = 100

# Permittivity info
eps_background = 1
eps_diel = index**2

# Add the layers to the grcwa model
obj.Add_LayerUniform(thick_top, eps_background)
for i in range(num_patterned_layers):
    obj.Add_LayerGrid(thick_layers[i], num_x, num_y)
obj.Add_LayerUniform(thick_bot, eps_background)

obj.Init_Setup(Gmethod=1)

if structures == []:
    eps_grid_substrate = np.ones((num_x, num_y)) * eps_background
    eps_grid_L1 = np.ones((num_x, num_y)) * eps_background
    eps_grid_L2 = np.ones((num_x, num_y)) * eps_background
elif len(structures) == 1:
    eps_grid_substrate = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L1 = np.ones((num_x, num_y)) * eps_background
    eps_grid_L2 = np.ones((num_x, num_y)) * eps_background
elif len(structures) == 2:
    eps_grid_substrate = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L1 = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L2 = np.ones((num_x, num_y)) * eps_background
else:

    eps_grid_substrate = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L1 = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L2 = np.ones((num_x, num_y)) * eps_diel

    # For each layer, we need to create a permittivity mask
    sim_center_rcwa = simulation.center
    sim_size_rcwa = [sim_size[0], size_y]

    # Create a grid of all possible coordinates
    x0 = np.linspace(
        sim_center_rcwa[0] - sim_size_rcwa[0] / 2,
        sim_center_rcwa[0] + sim_size_rcwa[0] / 2,
        num_x,
    )
    y0 = np.linspace(
        sim_center_rcwa[1] - sim_size_rcwa[1] / 2,
        sim_center_rcwa[1] + sim_size_rcwa[1] / 2,
        num_y,
    )
    x_sim, y_sim = np.meshgrid(x0, y0, indexing="ij")

    # Now mask out the coordinates that correspond to the dielectric regions
    center_L1 = grating_L1.geometry.center
    size_L1 = grating_L1.geometry.size

    center_L2 = grating_L2.geometry.center
    size_L2 = grating_L2.geometry.size

    def get_ind(x_grid, y_grid, diel_center, diel_size):
        """Get the anti-mask indices for a given dielectric slab."""
        ind = np.nonzero(
            (x_grid < diel_center[0] - diel_size[0] / 2)
            | (x_grid > diel_center[0] + diel_size[0] / 2)
            | (y_grid < diel_center[1] - diel_size[1] / 2)
            | (y_grid > diel_center[1] + diel_size[1] / 2)
        )
        return ind

    ind_L1 = get_ind(x_sim, x_sim, center_L1, size_L1)
    ind_L2 = get_ind(x_sim, x_sim, center_L2, size_L2)

    eps_grid_L1[ind_L1] = eps_background
    eps_grid_L2[ind_L2] = eps_background

# Combine the three layer masks
eps_grid = np.concatenate(
    (eps_grid_L2.flatten(), eps_grid_L1.flatten(), eps_grid_substrate.flatten())
)

# Apply these material masks to the model
obj.GridLayer_geteps(eps_grid)

# Set up the s-polarized plane wave source
planewave = {"p_amp": 1, "s_amp": 0, "p_phase": 0, "s_phase": 0}
obj.MakeExcitationPlanewave(
    planewave["p_amp"],
    planewave["p_phase"],
    planewave["s_amp"],
    planewave["s_phase"],
    order=0,
)

# Run grcwa to get the reflection and transmission efficiencies by order
R, T = obj.RT_Solve(normalize=1)
Ri, Ti = obj.RT_Solve(byorder=1)


def rcwa_order_index(orders_x, orders_y, obj, rcwa_data):
    """Helper function to extract data corresponding to particular order pairs."""
    ords = []
    out_data = []
    for order_y in orders_y:
        ords.append([])
        out_data.append([])
        for order_x in orders_x:
            order = [order_x, order_y]
            ords[-1].append(obj.G.tolist().index(order))
            out_data[-1].append(np.array(rcwa_data[ords[-1][-1]]))
    return ords, out_data


# Extract grcwa results at orders corresponding to those computed above by Tidy3D
r_ords, Ri_by_order = rcwa_order_index(
    data_r.orders_x,
    data_r.orders_y,
    obj,
    Ri,
)
t_ords, Ti_by_order = rcwa_order_index(
    data_t.orders_x,
    data_t.orders_y,
    obj,
    Ti,
)

Plot and compare diffracted power¶

Since this is essentially a 1D grating along x, we'll plot the power, which is normalized and corresponds to the grating efficiency, as a function of x for order 0 in y. The results are in excellent agreement with each other.

In [9]:

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

orders_x = data_r.orders_x
ax[0].plot(orders_x, data_r.power.sel(orders_y=0), "o-", label="Tidy3D")
ax[0].plot(orders_x, Ri_by_order[0], "s--", label="grcwa")
ax[0].set_title("Reflection")
ax[0].set_xlabel("Order along x")
ax[0].legend()

orders_x = data_t.orders_x
ax[1].plot(orders_x, data_t.power.sel(orders_y=0), "o-", label="Tidy3D")
ax[1].plot(orders_x, Ti_by_order[0], "s--", label="grcwa")
ax[1].set_title("Transmission")
ax[1].set_xlabel("Order along x")
ax[1].legend()

plt.show()

Plot and compare power amplitudes¶

We can also compare the transmitted complex power amplitude for each order and each polarization to those obtained via the grcwa package. The power amplitudes are complex and provide information about the phase difference among different orders.

Note that grcwa returns fields in Cartesian coordinates, while Tidy3D returns power amplitudes in the s and p polarization basis. Therefore, we will use some convenience methods to convert grcwa's fields to spherical coordinates before comparing the two solvers.

In [10]:

import xarray as xr

# get amplitudes from Tidy3D results
amps_sp = data_t.amps

# get amplitudes from grcwa

# we're sampling near-fields in the lower-most layer
layer = 4
# position above the simulation domain's bottom where we're sampling
z_offset = monitor_t.center[2] - (-sim_size[2] / 2)
amps_grcwa_xy, _ = obj.Solve_FieldFourier(layer, z_offset)

# Extract grcwa results at orders corresponding to those computed by Tidy3D
amps_grcwa_xy = [
    np.array(rcwa_order_index(data_t.orders_x, data_t.orders_y, obj, amps)[1])
    for amps in amps_grcwa_xy
]

# we need to swap the axes as below for the data to match Tidy3D data
amps_grcwa_xy = [np.swapaxes(amps, 0, 1) for amps in amps_grcwa_xy]

# to match the format of Tidy3D data, add a frequency dimension
amps_grcwa_xy = [amps[..., None] for amps in amps_grcwa_xy]

# convert to spherical coordinates
theta, phi = data_t.angles
amps_grcwa_tp = data_t.monitor.car_2_sph_field(
    amps_grcwa_xy[0], amps_grcwa_xy[1], amps_grcwa_xy[2], theta.values, phi.values
)

# make an xarray dataset for the rcwa amplitudes
coords = {}
coords["orders_x"] = np.atleast_1d(data_t.orders_x)
coords["orders_y"] = np.atleast_1d(data_t.orders_y)
coords["f"] = np.array(data_t.f)
coords["polarization"] = ["s", "p"]
amps_grcwa_sp = xr.DataArray(
    np.stack([amps_grcwa_tp[2], amps_grcwa_tp[1]], axis=3), coords=coords
)

# finally, we can compare the complex amplitudes for the y=0 order, as a function of orders along x
pol = "p"
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
orders_x = data_t.orders_x
data_tidy3d = data_t.amps.sel(polarization=pol).values[:, 0]
data_grcwa = amps_grcwa_sp.sel(polarization=pol).values[:, 0]

ax[0].plot(orders_x, np.real(data_tidy3d), "o-", label="Tidy3D")
ax[0].plot(orders_x, np.real(data_grcwa), "s--", label="grcwa")
ax[0].set_title(f"{pol}-pol, real part")
ax[0].set_ylabel("Transmission")
ax[0].set_xlabel("Order along x")
ax[0].legend()

ax[1].plot(orders_x, np.imag(data_tidy3d), "o-", label="Tidy3D")
ax[1].plot(orders_x, np.imag(data_grcwa), "s--", label="grcwa")
ax[1].set_title(f"{pol}-pol, imag part")
ax[1].set_ylabel("Transmission")
ax[1].set_xlabel("Order along x")
ax[1].legend()

plt.show()

Plot Field Distributions¶

Plot the frequency-domain electric field components at the center frequency along an xz cut of the domain. The y component is zero everywhere because the wave is entirely x-polarized.

In [11]:

data = sim_data["field_xz"]
fig, axs = plt.subplots(1, 4, tight_layout=True, figsize=(12, 5))
sim_data.simulation.plot(y=0, ax=axs[0])
sim_data.plot_field("field_xz", field_name="Ex", val="abs", ax=axs[1])
sim_data.plot_field("field_xz", field_name="Ey", val="abs", ax=axs[2])
sim_data.plot_field("field_xz", field_name="Ez", val="abs", ax=axs[3])
plt.show()

Off-normal incidence¶

To study a diffraction grating under an angled illumination, we need to use Bloch boundary conditions along the x and y directions that match the incoming wave. The Bloch vector can be automatically computed based on the plane wave parameters, as shown below.

If the illumination is such that the angle is non-zero along y, a slightly special treatment is required in that we cannot set the simulation size to zero in that direction anymore. Instead we should define the simulation length to be a small finite value, and set the mesh step in that direction to the same value.

In [12]:

# Angles
angle_theta = np.pi / 10
angle_phi = np.pi / 3

# Simulation size
length_y = 0.01  # needed when angle_phi is not zero
sim_size = (length_x, length_y, length_z)

# Angled source
src_z = length_z / 2 - space_above / 2
angle_theta = np.pi / 10
source = td.PlaneWave(
    center=(0, 0, src_z),
    size=(td.inf, td.inf, 0),
    source_time=gaussian,
    direction="-",
    pol_angle=0,
    angle_theta=angle_theta,
    angle_phi=angle_phi,
)

# Boundaries
bloch_x = td.Boundary.bloch_from_source(source=source, domain_size=sim_size[0], axis=0)
bloch_y = td.Boundary.bloch_from_source(source=source, domain_size=sim_size[1], axis=1)
boundary_spec = td.BoundarySpec(
    x=bloch_x,
    y=bloch_y,
    z=td.Boundary(
        minus=td.PML(num_layers=num_pml_layers), plus=td.PML(num_layers=num_pml_layers)
    ),
)

# Simulation
simulation = td.Simulation(
    size=sim_size,
    grid_spec=td.GridSpec(
        grid_x=td.AutoGrid(min_steps_per_wvl=60),
        grid_y=td.UniformGrid(dl=sim_size[1]),
        grid_z=td.AutoGrid(min_steps_per_wvl=60),
    ),
    structures=structures,
    sources=[source],
    monitors=monitors,
    run_time=run_time * 2,
    boundary_spec=boundary_spec,
)

fig, ax = plt.subplots(1, 1, figsize=(5, 8))
simulation.plot(y=0, ax=ax)
plt.show()

In [13]:

sim_data = web.run(
    simulation,
    task_name="GratingEfficiency",
    path="data/GratingEfficiency.hdf5",
    verbose=True,
)

[16:30:12] Created task 'GratingEfficiency' with task_id                        webapi.py:139
           'fdve-01b65d4a-0b85-4644-b85e-fd96e3b41f19v1'.

           View task using web UI at 'https://tidy3d.simulation.cloud/workbench webapi.py:141
           ?taskId=fdve-01b65d4a-0b85-4644-b85e-fd96e3b41f19v1'.

Output()

[16:30:14] status = queued                                                      webapi.py:271

Output()

[16:30:16] status = preprocess                                                  webapi.py:265

[16:30:22] Maximum FlexCredit cost: 0.649. Use 'web.real_cost(task_id)' to get  webapi.py:288
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:292

           running solver                                                       webapi.py:302

Output()

[16:32:49] early shutoff detected, exiting.                                     webapi.py:316

           status = postprocess                                                 webapi.py:333

Output()

[16:32:58] status = success                                                     webapi.py:340

Output()

[16:33:08] loading SimulationData from data/GratingEfficiency.hdf5              webapi.py:512

[16:33:08] WARNING: updating Simulation from 2.1 to 2.2                             log.py:50

Let's look at the same data as before.

In [18]:

data_r = sim_data["reflection"]
data_t = sim_data["transmission"]
total_power = data_r.power.sum() + data_t.power.sum()
theta = data_t.angles[0].squeeze(drop=True).values

print(f"Total power: {total_power.values}")
print("Theta (degrees):", ", ".join(f"{t:.2f}" for t in theta * 180 / np.pi))
print(f"Amplitude data: \n{data_t.amps}")

Total power: 1.000271117754282
Theta (degrees): 59.97, 41.10, 27.56, 18.00, 16.04, 23.30, 35.43, 51.66, 85.40
Amplitude data: 
<xarray.DataArray (orders_x: 9, orders_y: 1, f: 1, polarization: 2)>
array([[[[-0.05664602+0.05397184j, -0.08835355+0.0328627j ]]],


       [[[ 0.03238692+0.14692346j,  0.0747875 +0.24756416j]]],


       [[[-0.0848925 -0.16559587j, -0.25389771-0.37728042j]]],


       [[[ 0.01336718-0.00702349j, -0.44102546+0.21494824j]]],


       [[[ 0.08295825+0.25085525j, -0.08892292-0.31271915j]]],


       [[[ 0.14513272-0.20704379j, -0.05325255+0.069076j  ]]],


       [[[-0.01250289+0.08235093j, -0.02001391+0.01585265j]]],


       [[[-0.02789081-0.17130405j, -0.04793705+0.01526262j]]],


       [[[ 0.06931859-0.07539692j,  0.00913786-0.00718715j]]]])
Coordinates:
  * orders_x      (orders_x) int64 -3 -2 -1 0 1 2 3 4 5
  * orders_y      (orders_y) int64 0
  * f             (f) float64 7.207e+14
  * polarization  (polarization) <U1 's' 'p'

We still used a P-polarized source (pol_angle = 0 in the source defition), but now we also get nonzero polarization = "s" amplitudes in the data. This is because the input polarization is at an angle with respect to the x-axis, and is not preserved by the translational invariance along y. Similarly, now there is a nonzero Ey field component in the recorded near fields.

In [15]:

data = sim_data["field_xz"]
fig, axs = plt.subplots(1, 4, tight_layout=True, figsize=(12, 5))
sim_data.simulation.plot(y=0, ax=axs[0])
sim_data.plot_field("field_xz", field_name="Ex", val="abs", ax=axs[1])
sim_data.plot_field("field_xz", field_name="Ey", val="abs", ax=axs[2])
sim_data.plot_field("field_xz", field_name="Ez", val="abs", ax=axs[3])
plt.show()

Comparison to RCWA¶

We can also compare to RCWA results as before. Note that, because we are injecting backwards, we need to add an np.pi to angle_theta in the grcwa computation.

In [16]:

# Define lattice constants - size of the domain
# grcwa requires a finite non-zero size along each dimension
size_y = 1e-3
L1 = [sim_size[0], 0]
L2 = [0, size_y]

# Set truncation order
nG = 300

# Set up RCWA object
freq = f0 / td.C_0  # grcwa uses normalized coordinates where the speed of light is 1
obj = grcwa.obj(nG, L1, L2, freq, angle_theta + np.pi, angle_phi, verbose=0)

# Set up the geometry (the layers are ordered top to bottom in grcwa)
num_patterned_layers = 3
thick_top = space_above
thick_layers = [height_layer1, height_layer2, height_substrate]
thick_bot = space_below

# Discretization points along x and y
num_x = 300
num_y = 100

# Permittivity info
eps_background = 1
eps_diel = index**2

# Add the layers to the grcwa model
obj.Add_LayerUniform(thick_top, eps_background)
for i in range(num_patterned_layers):
    obj.Add_LayerGrid(thick_layers[i], num_x, num_y)
obj.Add_LayerUniform(thick_bot, eps_background)

obj.Init_Setup(Gmethod=1)

if structures == []:
    eps_grid_substrate = np.ones((num_x, num_y)) * eps_background
    eps_grid_L1 = np.ones((num_x, num_y)) * eps_background
    eps_grid_L2 = np.ones((num_x, num_y)) * eps_background
elif len(structures) == 1:
    eps_grid_substrate = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L1 = np.ones((num_x, num_y)) * eps_background
    eps_grid_L2 = np.ones((num_x, num_y)) * eps_background
elif len(structures) == 2:
    eps_grid_substrate = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L1 = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L2 = np.ones((num_x, num_y)) * eps_background
else:

    eps_grid_substrate = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L1 = np.ones((num_x, num_y)) * eps_diel
    eps_grid_L2 = np.ones((num_x, num_y)) * eps_diel

    # For each layer, we need to create a permittivity mask
    sim_center_rcwa = simulation.center
    sim_size_rcwa = [sim_size[0], size_y]

    # Create a grid of all possible coordinates
    x0 = np.linspace(
        sim_center_rcwa[0] - sim_size_rcwa[0] / 2,
        sim_center_rcwa[0] + sim_size_rcwa[0] / 2,
        num_x,
    )
    y0 = np.linspace(
        sim_center_rcwa[1] - sim_size_rcwa[1] / 2,
        sim_center_rcwa[1] + sim_size_rcwa[1] / 2,
        num_y,
    )
    x_sim, y_sim = np.meshgrid(x0, y0, indexing="ij")

    # Now mask out the coordinates that correspond to the dielectric regions
    center_L1 = grating_L1.geometry.center
    size_L1 = grating_L1.geometry.size

    center_L2 = grating_L2.geometry.center
    size_L2 = grating_L2.geometry.size

    def get_ind(x_grid, y_grid, diel_center, diel_size):
        """Get the anti-mask indices for a given dielectric slab."""
        ind = np.nonzero(
            (x_grid < diel_center[0] - diel_size[0] / 2)
            | (x_grid > diel_center[0] + diel_size[0] / 2)
            | (y_grid < diel_center[1] - diel_size[1] / 2)
            | (y_grid > diel_center[1] + diel_size[1] / 2)
        )
        return ind

    ind_L1 = get_ind(x_sim, x_sim, center_L1, size_L1)
    ind_L2 = get_ind(x_sim, x_sim, center_L2, size_L2)

    eps_grid_L1[ind_L1] = eps_background
    eps_grid_L2[ind_L2] = eps_background

# Combine the three layer masks
eps_grid = np.concatenate(
    (eps_grid_L2.flatten(), eps_grid_L1.flatten(), eps_grid_substrate.flatten())
)

# Apply these material masks to the model
obj.GridLayer_geteps(eps_grid)

# Set up the s-polarized plane wave source
planewave = {"p_amp": 1, "s_amp": 0, "p_phase": 0, "s_phase": 0}
obj.MakeExcitationPlanewave(
    planewave["p_amp"],
    planewave["p_phase"],
    planewave["s_amp"],
    planewave["s_phase"],
    order=0,
)

# Run grcwa to get the reflection and transmission efficiencies by order
R, T = obj.RT_Solve(normalize=1)
Ri, Ti = obj.RT_Solve(byorder=1)


def rcwa_order_index(orders_x, orders_y, obj, rcwa_data):
    """Helper function to extract data corresponding to particular order pairs."""
    ords = []
    out_data = []
    for order_y in orders_y:
        ords.append([])
        out_data.append([])
        for order_x in orders_x:
            order = [order_x, order_y]
            ords[-1].append(obj.G.tolist().index(order))
            out_data[-1].append(np.array(rcwa_data[ords[-1][-1]]))
    return ords, out_data


# Extract grcwa results at orders corresponding to those computed above by Tidy3D
r_ords, Ri_by_order = rcwa_order_index(
    data_r.orders_x,
    data_r.orders_y,
    obj,
    Ri,
)
t_ords, Ti_by_order = rcwa_order_index(
    data_t.orders_x,
    data_t.orders_y,
    obj,
    Ti,
)

In [17]:

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

orders_x = data_r.orders_x
ax[0].plot(orders_x, data_r.power.sel(orders_y=0), "o-", label="Tidy3D")
ax[0].plot(orders_x, Ri_by_order[0], "s--", label="grcwa")
ax[0].set_title("Reflection")
ax[0].set_xlabel("Order along x")
ax[0].legend()

orders_x = data_t.orders_x
ax[1].plot(orders_x, data_t.power.sel(orders_y=0), "o-", label="Tidy3D")
ax[1].plot(orders_x, Ti_by_order[0], "s--", label="grcwa")
ax[1].set_title("Transmission")
ax[1].set_xlabel("Order along x")
ax[1].legend()

plt.show()