Uniform grating coupler in BeamZ

Unknown · 2026-09-21

Uniform grating coupler in BeamZ

This notebook reproduces the workflow of Flexcompute's Tidy3D uniform grating coupler example using BeamZ geometry, sources, monitors, mode solving, and local execution.

BeamZ's tilted GaussianBeamSource is currently a 3D source. The fast design-loop model therefore uses a two-cell-wide 3D strip with PEC sidewalls, which is the native TE-equivalent of the reference notebook's y-periodic 2D model. The final model is fully 3D.

import gc
import os
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display

import beamz as bz
from beamz.analysis import mode_data_to_dataframe

um = bz.um

1. Initial design

We target 1.50–1.60 µm on an SOI stack. As in the reference notebook, silicon and silica are treated as nondispersive over this narrow band. The 260 nm silicon layer is etched by 160 nm, leaving a 100 nm slab. The uniform grating has 20 teeth and an 80% silicon fill fraction.

lambda0 = 1.55 * um
freq0 = bz.LIGHT_SPEED / lambda0
wavelengths = np.linspace(1.50, 1.60, 101) * um
freqs = bz.LIGHT_SPEED / wavelengths
fwidth = 0.5 * (np.max(freqs) - np.min(freqs))

n_si = 3.47
n_sio2 = 1.44
mat_air = bz.Material(permittivity=1.0)
mat_si = bz.Material(permittivity=n_si**2)
mat_sio2 = bz.Material(permittivity=n_sio2**2)

t_si = 0.26 * um
etch_depth = 0.16 * um
fill_fraction = 0.80
t_top_oxide = 0.68 * um
t_box = 2.0 * um
num_teeth = 20

theta = np.deg2rad(14.5)
mode_field_diameter = 10.8 * um
buffer = 0.60 * lambda0
source_gap = 0.50 * um
source_time = bz.GaussianPulse(freq0=freq0, fwidth=fwidth, offset=4.0)
mode_spec = bz.ModeSpec(num_modes=1, target_neff=n_si, polarization="te")

print(f"Silicon slab remaining after etch: {(t_si - etch_depth) / um:.2f} µm")
Silicon slab remaining after etch: 0.10 µm

Geometry helpers

The objects are added in painter's order: oxide and substrate first, the residual silicon slab second, and the teeth last. Each tooth restores the full 260 nm silicon thickness over its filled part of the period.

def box_from_bounds(lower, upper, material):
    lower = np.asarray(lower, dtype=float)
    upper = np.asarray(upper, dtype=float)
    return bz.Box(
        center=tuple(0.5 * (lower + upper)),
        size=tuple(upper - lower),
        material=material,
    )


def add_vertical_stack(design, bounds, *, z_offset=0.0):
    x_min, x_max, y_min, y_max, z_min, z_max = bounds
    design += box_from_bounds(
        (x_min, y_min, z_offset),
        (x_max, y_max, z_offset + t_top_oxide),
        mat_sio2,
    )
    design += box_from_bounds(
        (x_min, y_min, z_offset - t_box),
        (x_max, y_max, z_offset),
        mat_sio2,
    )
    design += box_from_bounds(
        (x_min, y_min, z_min),
        (x_max, y_max, z_offset - t_box),
        mat_si,
    )
    return design


def add_grating_teeth(
    design, period, y_min, y_max, *, x_offset=0.0, z_offset=0.0
):
    for index in range(num_teeth):
        x0 = x_offset + index * period
        x1 = x0 + fill_fraction * period
        design += box_from_bounds(
            (x0, y_min, z_offset + t_si - etch_depth),
            (x1, y_max, z_offset + t_si),
            mat_si,
        )
    return design

Parameterized design-loop simulation

The source is a 1 W Gaussian beam tilted by 14.5° toward −x and polarized along y. A mode monitor on the input slab measures the power coupled into the fundamental mode travelling toward −x. In BeamZ, angle_phi=np.pi selects that negative-x tangential wavevector for a -z beam.

def make_design_loop_sim(period, source_x):
    grating_length = num_teeth * period
    steps_per_wavelength = 25
    resolution = lambda0 / (steps_per_wavelength * n_si)
    # Keep at least two y cells: GaussianBeamSource is a 3D TF/SF source.
    strip_width = max(2.1 * resolution, 0.12 * um)

    reference_lower = np.array(
        (-buffer, -0.5 * strip_width, -t_box - buffer), dtype=float
    )
    reference_upper = np.array(
        (grating_length + buffer, 0.5 * strip_width, t_si + t_top_oxide + buffer),
        dtype=float,
    )
    # BeamZ domains are centered at the origin. Translate every geometry,
    # source, and monitor coordinate together; Tidy3D expresses this as
    # Simulation(center=sim_box.center, size=sim_box.size).
    shift = -0.5 * (reference_lower + reference_upper)
    lower = reference_lower + shift
    upper = reference_upper + shift
    x_min, y_min, z_min = lower
    x_max, y_max, z_max = upper
    x0, _, z0 = shift

    design = add_vertical_stack(
        bz.Design(background=mat_air),
        (x_min, x_max, y_min, y_max, z_min, z_max),
        z_offset=z0,
    )
    # Unetched input slab (reference x < 0).
    design += box_from_bounds(
        (x_min, y_min, z0), (x0, y_max, z0 + t_si), mat_si
    )
    # Residual 100 nm slab beneath the grating (reference x > 0).
    design += box_from_bounds(
        (x0, y_min, z0),
        (x_max, y_max, z0 + t_si - etch_depth),
        mat_si,
    )
    design = add_grating_teeth(
        design, period, y_min, y_max, x_offset=x0, z_offset=z0
    )

    pml_thickness = 12 * resolution
    source_half_width = min(
        mode_field_diameter,
        source_x - reference_lower[0] - pml_thickness,
        reference_upper[0] - pml_thickness - source_x,
    )
    source = bz.GaussianBeamSource(
        center=(source_x + x0, 0.0, t_top_oxide + source_gap + z0),
        size=(2 * source_half_width, strip_width, 0.0),
        source_time=source_time,
        direction="-z",
        angle_theta=theta,
        angle_phi=np.pi,
        pol_angle=np.pi / 2,
        waist_radius=0.5 * mode_field_diameter,
        background_index=1.0,
        wavelength=lambda0,
        power=1.0,
    )
    monitor = bz.ModeMonitor(
        center=(reference_lower[0] + pml_thickness + 2 * resolution + x0, 0.0, 0.5 * t_si + z0),
        size=(0.0, strip_width, 6 * t_si),
        freqs=freqs,
        mode_spec=mode_spec,
        name="mode",
    )
    simulation = bz.Simulation(
        domain=tuple(reference_upper - reference_lower),
        design=design,
        sources=(source,),
        monitors=(monitor,),
        boundaries=(bz.PML(formulation="cpml"),),
        grid_spec=bz.GridSpec.uniform(resolution),
        run_time=2.0e-12,
    )
    return simulation, {
        "source": source,
        "monitor": monitor,
        "bounds": (x_min, x_max, y_min, y_max, z_min, z_max),
        "shift": shift,
    }

sim0, preview = make_design_loop_sim(
    period=0.50 * um, source_x=(5.0) * um
)
print(f"Design-loop grid: {sim0.grid.shape} | time steps: {sim0.num_steps}")

# sim.plot() returns the standard xy/xz pair. The thin-strip xy view is not
# informative, so retain the reference notebook's single xz setup view.
fig, axes = sim0.plot(
    y=0.0,
    z=t_si - 0.5 * etch_depth,
    figsize=(10.5, 4.8),
    width_ratios=(0.001, 1.0),
    source_markers=True,
    monitor_markers=True,
    show=False,
)
axes = np.asarray(axes).flat
axes[0].set_visible(False)
axes[1].grid(False)
axes[1].set_title("Design-loop grating: xz cross-section at y=0")
plt.show()

Design-loop grid: (664, 7, 269) | time steps: 58711
Output

2. Estimate the grating period from slab modes

For first-order coupling,

$$ \Lambda = \frac{\lambda_0} {n_{\mathrm{eff}} - n_{\mathrm{clad}}\sin\theta_{\mathrm{clad}}}, \qquad n_{\mathrm{eff}} \approx f n_{\mathrm{unetched}} + (1-f)n_{\mathrm{etched}}. $$

Like ModeSimulation.from_simulation in the reference, the two BeamZ ModeSource.solve_modes calls reuse the exact simulation geometry and grid but move the solve plane to an unetched and an etched cross-section. These objects are solve requests only; they are not attached as FDTD sources.

shift_2d = np.asarray(preview["shift"])
unetched_probe = bz.ModeSource(
    center=(-0.10 * um + shift_2d[0], 0.0, 0.5 * t_si + shift_2d[2]),
    size=(0.0, sim0.design.height, 5 * t_si),
    source_time=source_time,
    direction="-",
    mode_spec=mode_spec,
)
etched_probe = unetched_probe.updated_copy(
    center=(
        num_teeth * 0.50 * um + 0.10 * um + shift_2d[0],
        0.0,
        0.5 * t_si + shift_2d[2],
    )
)

unetched_modes = unetched_probe.solve_modes(sim0, freqs=[freq0])
etched_modes = etched_probe.solve_modes(sim0, freqs=[freq0])
display(mode_data_to_dataframe(unetched_modes))
display(mode_data_to_dataframe(etched_modes))
wavelength n eff k eff loss (dB/cm) mode area
f mode_index
1.934145e+14 0 1.55 2.954398 0.0 0.0 0.039351
wavelength n eff k eff loss (dB/cm) mode area
f mode_index
1.934145e+14 0 1.55 2.188442 0.0 0.0 0.040381 
neff_unetched = float(np.real(np.asarray(unetched_modes.neffs)[0, 0]))
neff_etched = float(np.real(np.asarray(etched_modes.neffs)[0, 0]))

# Snell's law converts the 14.5° air angle to the angle inside the oxide.
theta_cladding = np.arcsin(np.sin(theta) / n_sio2)
weighted_neff = (
    fill_fraction * neff_unetched + (1 - fill_fraction) * neff_etched
)
estimated_period = lambda0 / (
    weighted_neff - n_sio2 * np.sin(theta_cladding)
)

print(f"n_eff unetched = {neff_unetched:.4f}")
print(f"n_eff etched   = {neff_etched:.4f}")
print(f"Estimated grating period = {estimated_period / bz.nm:.0f} nm")

n_eff unetched = 2.9544
n_eff etched   = 2.1884
Estimated grating period = 608 nm

3. Sweep period and fiber position

To keep local execution practical, the sweep samples 3 periods × 3 source positions, for 9 simulations total. This preserves the reference workflow's two design variables while bounding the runtime. BeamZ executes locally, so the dictionary of immutable simulations is followed by an explicit loop rather than a cloud batch.

periods = np.linspace(0.600, 0.630, 3) * um
source_positions = (np.linspace(3.0, 5.0, 3)) * um

simulations_2d = {
    f"p={period / um:.4f};source_x={source_x / um:.4f}": make_design_loop_sim(
        period, source_x
    )[0]
    for period in periods
    for source_x in source_positions
}
print(f"Prepared {len(simulations_2d)} immutable design-loop simulations.")

Prepared 9 immutable design-loop simulations.

sweep_results = {}
for index, (name, simulation) in enumerate(simulations_2d.items(), start=1):
    print(f"[{index:02d}/{len(simulations_2d):02d}] {name}")
    sweep_results[name] = simulation.run(
        progress=False,
        backend="cuda_streamed",
        performance=False,
    )
    # Each period realizes a slightly different domain shape. Evict compiled
    # plans as we go so a local sweep does not retain every device allocation.
    simulation.clear_compiled_cache()
    gc.collect()

[01/09] p=0.6000;source_x=3.0000

[02/09] p=0.6000;source_x=4.0000

[03/09] p=0.6000;source_x=5.0000

[04/09] p=0.6150;source_x=3.0000

[05/09] p=0.6150;source_x=4.0000

[06/09] p=0.6150;source_x=5.0000

[07/09] p=0.6300;source_x=3.0000

[08/09] p=0.6300;source_x=4.0000

[09/09] p=0.6300;source_x=5.0000

coupling_efficiency = np.empty((len(periods), len(source_positions)))
for i, period in enumerate(periods):
    for j, source_x in enumerate(source_positions):
        name = f"p={period / um:.4f};source_x={source_x / um:.4f}"
        mode_data = sweep_results[name].mode("mode")
        amplitude = (
            mode_data.amps.sel(direction="-", mode_index=0)
            .sel(f=freq0, method="nearest")
            .item()
        )
        coupling_efficiency[i, j] = np.abs(amplitude) ** 2

best_index = np.unravel_index(
    np.nanargmax(coupling_efficiency), coupling_efficiency.shape
)
best_period = periods[best_index[0]]
best_source_x = source_positions[best_index[1]]
best_efficiency = coupling_efficiency[best_index]

fig, ax = plt.subplots(figsize=(7.2, 4.6), constrained_layout=True)
mesh = ax.pcolormesh(
    source_positions / um,
    periods / bz.nm,
    100 * coupling_efficiency,
    shading="auto",
    cmap="RdBu_r",
)
colorbar = fig.colorbar(mesh, ax=ax)
colorbar.set_label("Coupling efficiency (%)")
ax.plot(best_source_x / um, best_period / bz.nm, "ko", ms=5, label="best sample")
ax.set(
    title="Fiber-to-waveguide coupling sweep",
    xlabel="Gaussian-beam center x (µm)",
    ylabel="Grating period (nm)",
)
ax.grid(False)
ax.legend(frameon=False)
plt.show()

print(f"Optimal period = {best_period / bz.nm:.1f} nm")
print(f"Optimal source x = {best_source_x / um:.2f} µm")
print(
    f"Optimal coupling = {100 * best_efficiency:.2f}% "
    f"({10 * np.log10(max(best_efficiency, 1e-15)):.2f} dB)"
)

# The 12-points/λ 3D grid is large. Release every design-loop result and
# compiled plan before rasterizing it so notebook execution stays within
# the host-memory budget. The compact sweep arrays and selected optimum remain.
sim0.clear_compiled_cache()
for simulation in simulations_2d.values():
    simulation.clear_compiled_cache()
del sweep_results, simulations_2d, mode_data
del unetched_modes, etched_modes, sim0, preview
plt.close('all')
gc.collect()
Output
Optimal period = 615.0 nm
Optimal source x = 4.00 µm
Optimal coupling = 26.47% (-5.77 dB)

21

4. Three-dimensional grating coupler

This is the same finite-width geometry as the Tidy3D reference: a 10.8 µm-wide grating connected to a 500 nm waveguide by a 50 µm linear taper. The mode monitor is halfway through the straight output buffer, exactly as in the reference.

Tidy3D's GridSpec.auto(min_steps_per_wvl=15) is a nonuniform material-aware mesh. BeamZ currently offers a uniform grid here; using 15 steps per silicon wavelength everywhere would create about 149 million cells and is not safe on a local 24 GB GPU. The BeamZ run therefore uses 8 uniform steps per silicon wavelength with subpixel material averaging. This changes numerical resolution, not device geometry.

waveguide_width = 0.50 * um
grating_width = mode_field_diameter
taper_length = 50.0 * um


def make_3d_sim(period, source_x):
    grating_length = num_teeth * period
    reference_lower = np.array(
        (-buffer - taper_length, -0.5 * grating_width - buffer, -t_box - buffer),
        dtype=float,
    )
    reference_upper = np.array(
        (grating_length + buffer, 0.5 * grating_width + buffer, t_si + t_top_oxide + buffer),
        dtype=float,
    )
    shift = -0.5 * (reference_lower + reference_upper)
    lower = reference_lower + shift
    upper = reference_upper + shift
    x_min, y_min, z_min = lower
    x_max, y_max, z_max = upper
    x0, _, z0 = shift  # reference coordinates x=0 and z=0 after translation
    bounds = (x_min, x_max, y_min, y_max, z_min, z_max)

    design = add_vertical_stack(
        bz.Design(background=mat_air), bounds, z_offset=z0
    )
    # Residual etched slab across the grating region, matching the reference.
    design += box_from_bounds(
        (x0, y_min, z0),
        (x_max, y_max, z0 + t_si - etch_depth),
        mat_si,
    )
    design = add_grating_teeth(
        design,
        period,
        -0.5 * grating_width,
        0.5 * grating_width,
        x_offset=x0,
        z_offset=z0,
    )

    # Exact six-vertex Tidy3D taper, including its straight output extension.
    taper_vertices = (
        (x0, 0.5 * grating_width),
        (x0, -0.5 * grating_width),
        (x0 - taper_length, -0.5 * waveguide_width),
        (x0 - taper_length - 2 * buffer, -0.5 * waveguide_width),
        (x0 - taper_length - 2 * buffer, 0.5 * waveguide_width),
        (x0 - taper_length, 0.5 * waveguide_width),
    )
    design += bz.Polygon(
        vertices=taper_vertices,
        depth=t_si,
        z=z0,
        material=mat_si,
    )

    steps_per_wavelength = 12
    resolution = lambda0 / (steps_per_wavelength * n_si)
    source = bz.GaussianBeamSource(
        center=(source_x + x0, 0.0, t_top_oxide + source_gap + z0),
        size=(2 * mode_field_diameter, 2 * mode_field_diameter, 0.0),
        source_time=source_time,
        direction="-z",
        angle_theta=theta,
        angle_phi=np.pi,
        pol_angle=np.pi / 2,
        waist_radius=0.5 * mode_field_diameter,
        background_index=1.0,
        wavelength=lambda0,
        power=1.0,
    )
    monitor = bz.ModeMonitor(
        center=(x0 - taper_length - 0.5 * buffer, 0.0, 0.5 * t_si + z0),
        size=(0.0, 4 * waveguide_width, 6 * t_si),
        freqs=freqs,
        mode_spec=mode_spec,
        name="mode",
    )
    domain = reference_upper - reference_lower
    field_xz = bz.FieldMonitor(
        center=(0.0, 0.0, 0.0),
        size=(domain[0], 0.0, domain[2]),
        freqs=[freq0],
        fields=("Ex", "Ey", "Ez"),
        name="field_xz",
    )
    simulation = bz.Simulation(
        domain=tuple(domain),
        design=design,
        sources=(source,),
        monitors=(monitor, field_xz),
        boundaries=(bz.PML(formulation="cpml"),),
        grid_spec=bz.GridSpec.uniform(
            resolution,
            max_total_cells=100_000_000,
        ),
        run_time=2.0e-12,
    )
    return simulation, {
        "source": source,
        "monitor": monitor,
        "field_xz": field_xz,
        "bounds": bounds,
        "shift": shift,
    }

sim_3d, setup_3d = make_3d_sim(best_period, best_source_x)
cell_count = int(np.prod(sim_3d.grid.shape))
print(f"3D grid: {sim_3d.grid.shape} = {cell_count / 1e6:.2f} M cells")
print(f"Time steps: {sim_3d.num_steps}")

# BeamZ uses x and y as the transverse axes of a z-normal source plane.
# For direction=-z and phi=pi, the tangential wavevector points toward -x.
beam_direction = np.array(
    (-np.sin(theta), 0.0, -np.cos(theta)), dtype=float
)
print(
    "Gaussian-beam propagation direction: "
    f"k_hat = ({beam_direction[0]:+.4f}, {beam_direction[1]:+.4f}, "
    f"{beam_direction[2]:+.4f})"
)

fig, ax = plt.subplots(figsize=(4.3, 4.3), constrained_layout=True)
ax.plot((0.0, 0.0), (0.0, -1.05), color="0.55", ls="--", lw=1.5)
ax.arrow(
    0.0,
    0.0,
    beam_direction[0],
    beam_direction[2],
    width=0.012,
    head_width=0.075,
    head_length=0.10,
    length_includes_head=True,
    color="#d81b60",
)
ax.text(-0.08, -0.42, f"{np.rad2deg(theta):.1f}°", ha="right", va="center")
ax.text(-0.27, -0.93, "toward −x", ha="right", va="top", color="#d81b60")
ax.set(
    title="Gaussian-beam propagation direction",
    xlabel="x (normalized)",
    ylabel="z (normalized)",
    xlim=(-0.42, 0.15),
    ylim=(-1.10, 0.10),
)
ax.set_aspect("equal")
ax.grid(False)
plt.show()

fig, axes = sim_3d.plot(
    z=setup_3d["shift"][2] + t_si - 0.5 * etch_depth,
    y=0.0,
    figsize=(13, 5.0),
    width_ratios=(1.35, 1.0),
    source_markers=True,
    monitor_markers=True,
    show=False,
)
for axis in np.asarray(axes).flat:
    axis.grid(False)
fig.suptitle("Finite-width 3D grating coupler", y=1.03)
plt.show()

3D grid: (1724, 341, 129) = 75.84 M cells
Time steps: 28182
Gaussian-beam propagation direction: k_hat = (-0.2504, +0.0000, -0.9681)
OutputOutput
results_3d = sim_3d.run(
    progress=True,
    backend="cuda_streamed",
    performance=True,
)
* Compiling simulation...
* Running simulation: [--------------------] 0% (0/28182 steps)
* Running simulation: [#-------------------] 5% (1410/28182 steps)
* Running simulation: [##------------------] 10% (2820/28182 steps)
* Running simulation: [###-----------------] 15% (4230/28182 steps)
* Running simulation: [####----------------] 20% (5640/28182 steps)
* Running simulation: [#####---------------] 25% (7050/28182 steps)
* Running simulation: [######--------------] 30% (8460/28182 steps)
* Running simulation: [#######-------------] 35% (9870/28182 steps)
* Running simulation: [########------------] 40% (11280/28182 steps)
* Running simulation: [#########-----------] 45% (12690/28182 steps)
* Running simulation: [##########----------] 50% (14100/28182 steps)
* Running simulation: [###########---------] 55% (15510/28182 steps)
* Running simulation: [############--------] 60% (16920/28182 steps)
* Running simulation: [#############-------] 65% (18330/28182 steps)
* Running simulation: [##############------] 70% (19740/28182 steps)
* Running simulation: [###############-----] 75% (21150/28182 steps)
* Running simulation: [################----] 80% (22560/28182 steps)
* Running simulation: [#################---] 85% (23970/28182 steps)
* Running simulation: [##################--] 90% (25380/28182 steps)
* Running simulation: [###################-] 95% (26790/28182 steps)
* Running simulation: [####################] 100% (28182/28182 steps)
Simulation runtime: 380.36 s
GCUPS: 5.619

5. Field propagation at 1.55 µm

The center-frequency DFT field on the symmetry plane provides a direct visual check of the source convention. The real $E_y$ close-up shows wavefronts descending toward −x onto the grating. The electric-field magnitude then shows the coupled field propagating left through the 50 µm taper. Material boundaries are overlaid by BeamZ's public field-plotting helper; grid lines remain disabled.

x0_3d = setup_3d["shift"][0]
grating_length_3d = num_teeth * best_period
grating_xlim = (
    (x0_3d - 2.0 * um) / um,
    (x0_3d + grating_length_3d + buffer) / um,
)
domain_x_um = sim_3d.domain[0] / um

fig, ax = results_3d.plot_field(
    monitor_name="field_xz",
    field_name="Ey",
    frequency=freq0,
    val="real",
    cmap="RdBu_r",
    xlim=grating_xlim,
    figsize=(12.5, 4.2),
    show=False,
)
ax.set_title("Incident and diffracted field: Re($E_y$), grating close-up")
ax.grid(False)
plt.show()

fig, ax = results_3d.plot_field(
    monitor_name="field_xz",
    field_name="E",
    frequency=freq0,
    val="abs",
    cmap="magma",
    xlim=(-0.5 * domain_x_um, 0.5 * domain_x_um),
    figsize=(12.5, 4.0),
    show=False,
)
ax.set_title("Electric-field magnitude through the full grating and taper")
ax.set_aspect("auto")
ax.grid(False)
plt.show()
OutputOutput

6. Coupling spectrum

The mode monitor reports source-normalized complex amplitudes. Because the Gaussian beam requests 1 W, the squared magnitude of the backward fundamental mode amplitude is the fiber-to-waveguide coupling efficiency.

mode_3d = results_3d.mode("mode")
amplitude_3d = np.asarray(mode_3d.amps.sel(direction="-", mode_index=0))
coupling_3d = np.abs(amplitude_3d) ** 2s

fig, ax = plt.subplots(figsize=(7.4, 4.5), constrained_layout=True)
ax.plot(
    wavelengths / um,
    10 * np.log10(np.maximum(coupling_3d, 1e-15)),
    color="#d81b60",
    lw=2.2,
)
ax.set(
    title="3D grating-coupler spectrum",
    xlabel="Wavelength (µm)",
    ylabel="Coupling efficiency (dB)",
    ylim=(-12, 0),
)
ax.grid(False)
plt.show()

center_index = int(np.argmin(np.abs(wavelengths - lambda0)))
print(
    f"Coupling at {wavelengths[center_index] / um:.3f} µm: "
    f"{100 * coupling_3d[center_index]:.2f}% "
    f"({10 * np.log10(max(coupling_3d[center_index], 1e-15)):.2f} dB)"
)
Output
Coupling at 1.550 µm: 38.64% (-4.13 dB)