Skip to content

Tutorial for computing the emission from a cleaved stack using a 2D simulation #3169

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
stevengj merged 3 commits into from
Mar 20, 2026
Merged

Tutorial for computing the emission from a cleaved stack using a 2D simulation #3169

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0b6eb53
tutorial for computing the emission from a cleaved stack using a 2D s…
oskooi Mar 13, 2026
ed2cd24
update axis label of plot of flux vs. kz and provide additional detai…
oskooi Mar 14, 2026
62ef37d
Apply suggestion from @stevengj
stevengj Mar 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 23 additions & 0 deletions doc/docs/Python_Tutorials/Custom_Source.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -652,3 +652,26 @@ As a final note, computing the [extraction efficiency of an LED](Local_Density_o
delta_length
)
```

Dipole Emission of a Cleaved Multilayer Stack
---------------------------------------------

The emission from a cleaved facet of a multilayer stack (relevant for [edge-emitting semiconductor lasers](https://www.optica-opn.org/home/articles/volume_30/february_2019/features/semiconductor_lasers_for_3-d_sensing/)) can be modeled using a 2D simulation via a Fourier-series expansion of line-current sources. A schematic of the device geometry and simulation layout is shown below. This example involves computing the extraction efficiency of a dipole (wavelength of 1 μm) within an active region ($n$ = 2.3) surrounded by cladding layers ($n$ = 1.5) and placed on a substrate ($n$ = 3.2).

![](../images/edge_emitter_layout_2D.png#center)

There are three important things to note regarding the setup of this calculation:

- A 2D simulation (in $xy$) with out-of-plane wavevector ($k_z$) requires specifying the parameter [`kz_2d`](../2d_Cell_Special_kz.md) of the `Simulation` constructor to `"real/imag"` and using the default value for `force_complex_fields` of `False`. Alternatively, `kz_2d` can be set to "complex" (default value) or "3d" but this requires specifying `force_complex_fields` to `True` (to ensure that $k_z=0$ is handled consistently with $k_z \ne 0$, since the latter necessarily use complex fields).
- The Fourier-series expansion requires only *non-negative* wavevectors $k_z$ since the flux is the same for $\pm k_z$.
- The criteria for truncating the Fourier series is determined by when the *total* power (computed using the local density of states) has fallen below a given threshold (which needs to be small enough to ensure the results are converged).

A plot of the total and extracted power versus $k_z$ (in units of $\omega$) is shown below for an $E_x$ dipole positioned 0.32 μm from the edge facet. When $k_z > \omega$ (below the light line of air), the extracted power is essentially zero and coupling into guided substrate modes dominates. This is why the total power reaches a local minimum at $k_z = \omega$ and starts increasing, reaches a maximum, and eventually becomes negligible when $k_z$ is sufficiently large. The extraction efficiency computed from this data is 24.12%. This example involved 76 2D simulations.

![](../images/edge_emitter_flux_vs_kz.png#center)
Comment thread
oskooi marked this conversation as resolved.

The results of the 2D calculation can be validated using a 3D simulation. The table below shows the extraction efficiency computed using 2D and 3D simulations for three different dipole test cases at a resolution of 25 pixels/μm. The 2D and 3D results are in reasonable agreement.

![](../images/edge_emitter_2D_vs_3D_table.png#center)

The simulation script is [examples/edge_emitter_2D.py](https://github.com/NanoComp/meep/blob/master/python/examples/edge_emitter_2D.py). The script used for validating the results using a 3D simulation is [examples/edge_emitter_3D.py](https://github.com/NanoComp/meep/blob/master/python/examples/edge_emitter_3D.py).
Binary file added doc/docs/images/edge_emitter_2D_vs_3D_table.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added doc/docs/images/edge_emitter_flux_vs_kz.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added doc/docs/images/edge_emitter_layout_2D.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
234 changes: 234 additions & 0 deletions python/examples/edge_emitter_2D.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,234 @@
"""Computes the extraction efficiency of a cleaved stack using a 2D sim."""

import argparse
from typing import Tuple

import meep as mp
import matplotlib.pyplot as plt
import numpy as np


RESOLUTION_UM = 25
WAVELENGTH_UM = 1.0
CLADDING_UM = 0.5
ACTIVE_UM = 0.3
SUBSTRATE_UM = 1.0
PML_UM = 0.5
AIR_UM = 1.0
SIDE_UM = 3.0
N_SUBSTRATE = 3.2
N_CLADDING = 1.5
N_ACTIVE = 2.3
FIELD_DECAY_PERIOD = 25.0
FIELD_DECAY_TOL = 1e-6
DELTA_KZ = 0.02
DIPOLE_FLUX_TOL = 0.02
DEBUG_OUTPUT = 0


def line_current_in_cleaved_stack(
dipole_pol: str, dipole_pos_um: float, kz: float
) -> Tuple[float, float]:
"""Computes the flux of a line current in a cleaved multilayer stack.

Args:
dipole_pol: polarization of the dipole ("x", "y", or "z").
dipole_pos_um: position of the dipole relative to the edge of the
active region.
kz: z-component of the Bloch-periodic wavevector.

Returns:
The flux emitted by a line current in the active region and into air as
a 2-tuple.
"""
if dipole_pol == "x":
src_cmpt = mp.Ex
elif dipole_pol == "y":
src_cmpt = mp.Ey
elif dipole_pol == "z":
src_cmpt = mp.Ez
else:
raise ValueError("dipole_pol must be x, y, or z.")

cell_x_um = PML_UM + SIDE_UM + AIR_UM + PML_UM
cell_y_um = (
PML_UM + SUBSTRATE_UM + CLADDING_UM + ACTIVE_UM + CLADDING_UM + AIR_UM + PML_UM
)
cell_size = mp.Vector3(cell_x_um, cell_y_um, 0)

boundary_layers = [mp.PML(thickness=PML_UM)]

frequency = 1 / WAVELENGTH_UM

src_pos = mp.Vector3(
-0.5 * cell_x_um + PML_UM + SIDE_UM - dipole_pos_um,
-0.5 * cell_y_um + PML_UM + SUBSTRATE_UM + CLADDING_UM + 0.5 * ACTIVE_UM,
0,
)
sources = [
mp.Source(
src=mp.GaussianSource(frequency, fwidth=0.1 * frequency),
center=src_pos,
component=src_cmpt,
)
]

geometry = [
mp.Block(
material=mp.Medium(index=N_SUBSTRATE),
center=mp.Vector3(
-0.5 * cell_x_um + 0.5 * (PML_UM + SIDE_UM),
-0.5 * cell_y_um + 0.5 * (PML_UM + SUBSTRATE_UM),
0,
),
size=mp.Vector3(PML_UM + SIDE_UM, PML_UM + SUBSTRATE_UM, mp.inf),
),
mp.Block(
material=mp.Medium(index=N_CLADDING),
center=mp.Vector3(
-0.5 * cell_x_um + 0.5 * (PML_UM + SIDE_UM),
-0.5 * cell_y_um + PML_UM + SUBSTRATE_UM + 0.5 * CLADDING_UM,
0,
),
size=mp.Vector3(PML_UM + SIDE_UM, CLADDING_UM, mp.inf),
),
mp.Block(
material=mp.Medium(index=N_ACTIVE),
center=mp.Vector3(
-0.5 * cell_x_um + 0.5 * (PML_UM + SIDE_UM),
-0.5 * cell_y_um
+ PML_UM
+ SUBSTRATE_UM
+ CLADDING_UM
+ 0.5 * ACTIVE_UM,
0,
),
size=mp.Vector3(PML_UM + SIDE_UM, ACTIVE_UM, mp.inf),
),
mp.Block(
material=mp.Medium(index=N_CLADDING),
center=mp.Vector3(
-0.5 * cell_x_um + 0.5 * (PML_UM + SIDE_UM),
-0.5 * cell_y_um
+ PML_UM
+ SUBSTRATE_UM
+ CLADDING_UM
+ ACTIVE_UM
+ 0.5 * CLADDING_UM,
0,
),
size=mp.Vector3(PML_UM + SIDE_UM, CLADDING_UM, mp.inf),
),
]

sim = mp.Simulation(
resolution=RESOLUTION_UM,
cell_size=cell_size,
k_point=mp.Vector3(0, 0, kz),
kz_2d="real/imag",
force_complex_fields=False,
sources=sources,
geometry=geometry,
boundary_layers=boundary_layers,
)

flux_mon = sim.add_flux(
frequency,
0,
1,
mp.FluxRegion(
center=mp.Vector3(
-0.5 * cell_x_um + PML_UM + 0.5 * (cell_x_um - 2 * PML_UM),
0.5 * cell_y_um - PML_UM,
0,
),
size=mp.Vector3(cell_x_um - 2 * PML_UM, 0, 0),
),
mp.FluxRegion(
center=mp.Vector3(
0.5 * cell_x_um - PML_UM,
0.5 * cell_y_um - PML_UM - 0.5 * (cell_y_um - 2 * PML_UM),
0,
),
size=mp.Vector3(0, cell_y_um - 2 * PML_UM, 0),
),
mp.FluxRegion(
center=mp.Vector3(
0.5 * cell_x_um - PML_UM - +0.5 * AIR_UM, -0.5 * cell_y_um + PML_UM, 0
),
size=mp.Vector3(AIR_UM, 0, 0),
weight=-1.0,
),
)

if DEBUG_OUTPUT:
if mp.am_master():
fig, ax = plt.subplots()
sim.plot2D(ax=ax)
fig.savefig("edge_emitter_layout.png", dpi=150, bbox_inches="tight")

sim.run(
mp.dft_ldos(frequency, 0, 1),
until_after_sources=mp.stop_when_fields_decayed(
FIELD_DECAY_PERIOD,
src_cmpt,
src_pos,
FIELD_DECAY_TOL,
),
)

pixel_area = (1 / RESOLUTION_UM) ** 2
dipole_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * pixel_area

air_flux = mp.get_fluxes(flux_mon)[0]

return dipole_flux, air_flux


if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument(
"dipole_pol",
type=str,
choices=["x", "y", "z"],
help="polarization of the electric dipole (x, y, z)",
)
parser.add_argument(
"dipole_pos_um",
type=float,
help="position of the dipole relative to the edge of the cleaved facet",
)
args = parser.parse_args()

num_kz = 0
max_dipole_flux = 0
dipole_flux = []
air_flux = []

# Brillouin-zone integration over kz.
while True:
kz = num_kz * DELTA_KZ
dipole_flux_kz, air_flux_kz = line_current_in_cleaved_stack(
args.dipole_pol, args.dipole_pos_um, kz
)
print(f"kz:, {num_kz}, {kz:.2f}, {dipole_flux_kz:.2f}, " f"{air_flux_kz:.2f}")

dipole_flux.append(dipole_flux_kz)
air_flux.append(air_flux_kz)
num_kz += 1

if dipole_flux_kz > max_dipole_flux:
max_dipole_flux = dipole_flux_kz
elif (dipole_flux_kz / max_dipole_flux) < DIPOLE_FLUX_TOL:
break

dipole_flux = np.array(dipole_flux)
air_flux = np.array(air_flux)
extraction_efficiency = (air_flux[0] + 2 * np.sum(air_flux[1:])) / (
dipole_flux[0] + 2 * np.sum(dipole_flux[1:])
)

print(
f"extraction_efficiency:, {args.dipole_pol}, {args.dipole_pos_um:.2f},"
f" {100 * extraction_efficiency:.2f}%"
)
Loading
Loading