Grating efficiency benchmark#

Run this notebook in your browser using Binder.

In this example, we compute the grating efficiency of a multilevel diffraction grating whose design is inspired by the work of Oliva et al. (2011).

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.

[1]:
# standard python imports
import numpy as np
import matplotlib.pyplot as plt

# tidy3D import
import tidy3d as td
from tidy3d import web
[15:00:37] INFO     Using client version: 1.7.0                                __init__.py:95

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.

[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

[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. Also, we’ll set up two DiffractionMonitors, one for reflection, and the other for transmission.

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

# We'll request all orders up to 20, but only those within the light cone
# will be returned in the final results.
orders = list(range(-20, 21))
monitor_r = td.DiffractionMonitor(
    center=[0., 0., length_z / 2 - wavelength],
    size=[td.inf, td.inf, 0],
    freqs=[f0],
    name='reflection',
    orders_x = orders,
    orders_y = orders,
    normal_dir='+',
)

monitor_t = td.DiffractionMonitor(
    center=[0., 0., -length_z / 2 + wavelength],
    size=[td.inf, td.inf, 0],
    freqs=[f0],
    name='transmission',
    orders_x = orders,
    orders_y = orders,
    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.

[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()
[15:00:38] INFO     Auto meshing using wavelength 0.4160 defined from        grid_spec.py:473
                    sources.                                                                 
<Figure size 360x576 with 1 Axes>
../_images/notebooks_GratingEfficiency_10_2.png

Run simulation#

[6]:
sim_data = web.run(simulation, task_name='GratingEfficiency', path='data/GratingEfficiency.hdf5')
           INFO     Using Tidy3D credentials from stored file.                     auth.py:70
[15:00:39] INFO     Authentication successful.                                     auth.py:30
[15:00:40] INFO     Uploaded task 'GratingEfficiency' with task_id              webapi.py:117
                    'bd1ce049-3903-47dd-b6a7-41bdb1738812'.                                  
[15:00:44] INFO     status = queued                                             webapi.py:261
[15:00:47] INFO     Maximum flex unit cost: 0.10                                webapi.py:252
[15:00:49] INFO     status = preprocess                                         webapi.py:273
[15:00:56] INFO     starting up solver                                          webapi.py:277
[15:01:06] INFO     running solver                                              webapi.py:283
[15:01:27] INFO     early shutoff detected, exiting.                            webapi.py:294
           INFO     status = postprocess                                        webapi.py:300
[15:01:39] INFO     status = success                                            webapi.py:306
[15:01:40] INFO     downloading file "output/monitor_data.hdf5" to              webapi.py:578
                    "data/GratingEfficiency.hdf5"                                            
[15:01:41] INFO     loading SimulationData from data/GratingEfficiency.hdf5     webapi.py:400

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.

[7]:
data_r = sim_data["reflection"]
data_t = sim_data["transmission"]
print(f"Total power: {sum(sum(data_r.power.values + data_t.power.values))}")
print(f"\nTheta angles (degrees): \n{data_t.angles[0] * 180 / np.pi}")
print(f"\nAmplitude data: \n{data_t.amps}")
           INFO     Auto meshing using wavelength 0.4160 defined from        grid_spec.py:473
                    sources.                                                                 
Total power: [0.99804998]

Theta angles (degrees):
<xarray.DataArray (orders_x: 9, orders_y: 1, f: 1)>
array([[[63.09360835]],

       [[41.97530865]],

       [[26.47924028]],

       [[12.88158202]],

       [[ 0.        ]],

       [[12.88158202]],

       [[26.47924028]],

       [[41.97530865]],

       [[63.09360835]]])
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

Amplitude data:
<xarray.DiffractionDataArray (orders_x: 9, orders_y: 1, polarization: 2, f: 1)>
array([[[[ 5.12910205e-18-1.13761677e-17j],
         [ 3.13276362e-02-6.98744631e-02j]]],


       [[[-1.25720226e-17+7.75219696e-18j],
         [-9.83758908e-02+6.06819934e-02j]]],


       [[[ 3.74934258e-17-2.01077500e-17j],
         [ 3.04306561e-01-1.63185034e-01j]]],


       [[[ 3.58631267e-17+5.03730697e-17j],
         [ 2.92754575e-01+4.11202056e-01j]]],


       [[[ 0.00000000e+00+0.00000000e+00j],
         [-9.76695105e-03+5.36655787e-01j]]],


       [[[ 0.00000000e+00+0.00000000e+00j],
         [-2.24140802e-01+2.80535888e-01j]]],


       [[[ 0.00000000e+00+0.00000000e+00j],
         [ 1.24422173e-01-2.42128602e-01j]]],


       [[[ 0.00000000e+00+0.00000000e+00j],
         [ 1.59764689e-01+9.85030251e-02j]]],


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

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.

[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.

[9]:
fig, ax = plt.subplots(1, 2, figsize=(10, 3))

orders_x = data_r.orders_x
ax[0].plot(orders_x, data_r.power.values[:,0], 'o-', label="Tidy3D")
ax[0].plot(orders_x, Ri_by_order[0], 'x--r', 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.values[:,0], 'o-', label="Tidy3D")
ax[1].plot(orders_x, Ti_by_order[0], 'x--r', label="grcwa")
ax[1].set_title("Transmission")
ax[1].set_xlabel("Order along x")
ax[1].legend()

plt.show()

<Figure size 720x216 with 2 Axes>
../_images/notebooks_GratingEfficiency_18_1.png

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.

[10]:
# 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
amps_grcwa_xy = [np.swapaxes(amps, 0, 1) for amps in amps_grcwa_xy]
# to match the format of Tidy3D data, add a fourth frequency dimension
amps_grcwa_xy = [amps[..., None] for amps in amps_grcwa_xy]

# to allow easily manipulating this data, we'll make a data array with grcwa results
# in the same format in which Tidy3D results are stored
coords = {}
coords["orders_x"] = np.atleast_1d(data_t.orders_x)
coords["orders_y"] = np.atleast_1d(data_t.orders_y)
coords["polarization"] = ["x", "y", "z"]
coords["f"] = np.array(data_t.frequencies)
amps_grcwa_car = td.DiffractionDataArray(
    np.stack(amps_grcwa_xy, axis=2), coords=coords
)

# since the data is now in Tidy3D-friendly format, we can use a convenience function to
# convert it to spherical coordinates
amps_grcwa_sp = data_t.car_2_sph(amps_grcwa_car)

# rename the polarization dimension to s and p from theta, phi
amps_grcwa_sp["polarization"] = ["p", "s"]

# 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), 'x--r', 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), 'x--r', 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()
<Figure size 864x288 with 2 Axes>
../_images/notebooks_GratingEfficiency_20_1.png

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.

[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_name="Ex", field_monitor_name="field_xz", ax=axs[1], val="abs")
sim_data.plot_field(field_name="Ey", field_monitor_name="field_xz", ax=axs[2], val="abs")
sim_data.plot_field(field_name="Ez", field_monitor_name="field_xz", ax=axs[3], val="abs")
plt.show()
<Figure size 864x360 with 7 Axes>
../_images/notebooks_GratingEfficiency_22_1.png

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.

[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()
[15:01:53] INFO     Auto meshing using wavelength 0.4160 defined from        grid_spec.py:473
                    sources.                                                                 
<Figure size 360x576 with 1 Axes>
../_images/notebooks_GratingEfficiency_24_2.png
[14]:
sim_data = web.run(simulation, task_name='GratingEfficiency', path='data/GratingEfficiency.hdf5')
[15:20:06] INFO     Uploaded task 'GratingEfficiency' with task_id              webapi.py:117
                    'ed46428d-d308-49aa-b5f8-eb071361723b'.                                  
[15:20:10] INFO     status = queued                                             webapi.py:261
[15:20:19] INFO     Maximum flex unit cost: 0.71                                webapi.py:252
           INFO     status = preprocess                                         webapi.py:273
[15:20:23] INFO     starting up solver                                          webapi.py:277
[15:20:38] INFO     running solver                                              webapi.py:283
[15:53:53] INFO     early shutoff detected, exiting.                            webapi.py:294
           INFO     status = postprocess                                        webapi.py:300
[15:54:36] INFO     status = success                                            webapi.py:306
           INFO     downloading file "output/monitor_data.hdf5" to              webapi.py:578
                    "data/GratingEfficiency.hdf5"                                            
[15:54:38] INFO     loading SimulationData from data/GratingEfficiency.hdf5     webapi.py:400
[15]:
data_r = sim_data["reflection"]
data_t = sim_data["transmission"]
print(f"Total power: {sum(sum(data_r.power.values + data_t.power.values))}")
print(f"\nTheta angles (degrees): \n{data_t.angles[0] * 180 / np.pi}")
print(f"\nAmplitude data: \n{data_t.amps}")
           INFO     Auto meshing using wavelength 0.4160 defined from        grid_spec.py:473
                    sources.                                                                 
Total power: [1.00028878]

Theta angles (degrees):
<xarray.DataArray (orders_x: 9, orders_y: 1, f: 1)>
array([[[59.96510993]],

       [[41.09624023]],

       [[27.56093082]],

       [[18.        ]],

       [[16.03511623]],

       [[23.30441275]],

       [[35.4338644 ]],

       [[51.65673818]],

       [[85.39539255]]])
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

Amplitude data:
<xarray.DiffractionDataArray (orders_x: 9, orders_y: 1, polarization: 2, f: 1)>
array([[[[-0.05679835+0.05389631j],
         [-0.08841113+0.03277967j]]],


       [[[ 0.03221876+0.14691651j],
         [ 0.07450146+0.24764791j]]],


       [[[-0.08472633-0.16561827j],
         [-0.25348129-0.37753052j]]],


       [[[ 0.01335681-0.0070062j ],
         [-0.44122592+0.21452508j]]],


       [[[ 0.08273648+0.25097149j],
         [-0.08859536-0.31285101j]]],


       [[[ 0.14537099-0.20693623j],
         [-0.05336622+0.06900587j]]],


       [[[-0.01253039+0.08236032j],
         [-0.02003819+0.01576802j]]],


       [[[-0.02751617-0.17130176j],
         [-0.047999  +0.01511375j]]],


       [[[ 0.06972028-0.07527757j],
         [ 0.00913691-0.00714857j]]]])
Coordinates:
  * orders_x      (orders_x) int64 -3 -2 -1 0 1 2 3 4 5
  * orders_y      (orders_y) int64 0
  * polarization  (polarization) <U1 's' 'p'
  * f             (f) float64 7.207e+14
Attributes:
    long_name:  diffraction amplitude

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.

[16]:
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_name="Ex", field_monitor_name="field_xz", ax=axs[1], val="abs")
sim_data.plot_field(field_name="Ey", field_monitor_name="field_xz", ax=axs[2], val="abs")
sim_data.plot_field(field_name="Ez", field_monitor_name="field_xz", ax=axs[3], val="abs")
plt.show()
<Figure size 864x360 with 7 Axes>
../_images/notebooks_GratingEfficiency_28_1.png

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.

[17]:
# 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,
)
[18]:
fig, ax = plt.subplots(1, 2, figsize=(10, 3))

orders_x = data_r.orders_x
ax[0].plot(orders_x, data_r.power.values[:,0], 'o-', label="Tidy3D")
ax[0].plot(orders_x, Ri_by_order[0], 'x--r', 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.values[:,0], 'o-', label="Tidy3D")
ax[1].plot(orders_x, Ti_by_order[0], 'x--r', label="grcwa")
ax[1].set_title("Transmission")
ax[1].set_xlabel("Order along x")
ax[1].legend()

plt.show()
<Figure size 720x216 with 2 Axes>
../_images/notebooks_GratingEfficiency_31_1.png