Skip to content

Temperature and carrier maps

A heater warms a waveguide; a PN junction injects or removes free carriers. Both change the refractive index. This guide shows how to put that index change into a PhotonHub simulation from a map you already have: a temperature map from a heat simulation, a carrier-density map from a device simulation, or a formula you write yourself.

PhotonHub does not solve heat or carrier transport. The map is an input.

CPU solver only

Media built from a map run on the local CPU solver. A GPU run that carries one is refused before it starts. Run with device="cpu", the default for ph.run_local.

Describe the map

ph.SpatialMap holds a scalar on a rectilinear grid, in microns, in the frame you draw your structures in:

import numpy as np
import photonhub as ph

x = np.linspace(-1.0, 1.0, 41)        # um
z = np.linspace(-0.5, 0.8, 27)        # um
X, Z = np.meshgrid(x, z, indexing="ij")
T = 300.0 + 60.0 * np.exp(-(X / 0.8) ** 2) * np.exp(-((Z - 0.6) / 0.5) ** 2)  # kelvin

heat = ph.SpatialMap(x_um=x, y_um=(0.0,), z_um=z, values=T[:, None, :])

values has the shape (len(x_um), len(y_um), len(z_um)). Each coordinate list is strictly increasing and may be non-uniform. An axis given a single coordinate means the map does not vary along it. Above, the map is a cross-section of a straight heater that runs along y. Beyond the outermost node on an axis, the map keeps the edge value.

To sample a formula, use ph.SpatialMap.from_function(fn, x_um, y_um, z_um), where fn(x, y, z) takes arrays in microns. To read a file exported by another tool, load it into arrays with NumPy and pass the coordinates and values.

Make the medium

ph.perturbed_medium turns a base medium and one or more maps into a medium whose permittivity, and loss for carriers, is given at the map's nodes:

si = ph.Medium(permittivity=3.48 ** 2)

hot_si = ph.perturbed_medium(si, temperature=heat, dn_dT=1.8e-4,
                             reference_temperature_k=300.0)

core = ph.Structure(
    geometry=ph.Box(center_um=(0.0, 0.0, 0.11), size_um=(0.5, 40.0, 0.22)),
    medium=hot_si)

The geometry decides where the medium is; the map decides what it is there. Any geometry works: box, polygon, cylinder or sphere. Give each material that the heat reaches its own perturbed medium, with its own dn_dT, on the same map.

For carriers, pass excess electron and hole densities in cm^-3:

doped_si = ph.perturbed_medium(si, electrons_cm3=ne_map, holes_cm3=nh_map,
                               wavelength_um=1.55)

wavelength_um is required with carrier maps: give the wavelength your simulation runs at (see the models below).

Temperature and carriers can be given together; their index changes add. The maps must then share one grid. The base may also be a material library entry, such as ph.materials.get("Si"), which is evaluated at wavelength_um (1.55 um when omitted).

The models

Thermo-optic. First order in temperature, at each node:

n = n0 + dn_dT * (T - T_ref)

n0 is the base medium's index and dn_dT is yours. For orientation, silicon is about 1.8e-4 per kelvin near 300 K and 1550 nm, and silica about 1e-5.

Carrier plasma dispersion in silicon at 1550 nm. Two published fits are built in, selected with carrier_model=. Densities are in cm^-3 and the absorption change in 1/cm.

carrier_model index change absorption change
"nedeljkovic2011_1550" (default) -5.4e-22 Ne^1.011 - 1.53e-18 Nh^0.838 8.88e-21 Ne^1.167 + 5.84e-20 Nh^1.109
"soref_bennett1987_1550" -8.8e-22 Ne - 8.5e-18 Nh^0.8 8.5e-18 Ne + 6.0e-18 Nh
  • M. Nedeljkovic, R. Soref, G. Z. Mashanovich, "Free-carrier electrorefraction and electroabsorption modulation predictions for silicon over the 1-14 um infrared wavelength range", IEEE Photonics Journal 3(6), 1171-1180 (2011), doi:10.1109/JPHOT.2011.2171930. The model is the paper's eqs. (4) and (5) (p. 1175) with the coefficients of its Table 1, row 1.55 um (p. 1178). The fit was made over 3.2e17 to 4e19 cm^-3 electrons and 5e17 to 1e20 cm^-3 holes; outside that span it is an extrapolation.
  • R. A. Soref, B. R. Bennett, "Electrooptical effects in silicon", IEEE J. Quantum Electron. 23(1), 123-129 (1987), doi:10.1109/JQE.1987.1073206. The coefficients are the commonly quoted 1.55 um expressions for this paper. They have not yet been checked against the paper itself. Computed from them, 1e18 holes per cm^3 give an index change of -2.135e-3.

ph.carrier_plasma(electrons_cm3, holes_cm3, model) returns the two changes on their own, and ph.perturbation.CARRIER_MODELS lists the coefficients.

The absorption is carried as a conductivity at each node, sigma = eps0 * c * n * alpha, which reproduces alpha at every frequency for small loss. The fits themselves hold near 1550 nm. With carrier maps, perturbed_medium needs wavelength_um and refuses one outside 1.50 to 1.60 um. Only that argument is checked, not the band your simulation runs: a medium built for 1.55 um and run at 1.31 um uses the 1550 nm values without a warning. Inside the allowed range, the 1550 nm values can be off by up to about 8% in the index change and up to about 14% in the absorption change at 1.50 and 1.60 um. These figures interpolate between the 1.3, 1.55 and 2.0 um rows of Table 1 in Nedeljkovic 2011; the paper has no rows at 1.50 or 1.60 um.

Put it in a simulation

sim = ph.Simulation(
    structures=[core],
    background=ph.Background(permittivity=1.444 ** 2),
    ports=[ph.Port("in", (0.0, -15.0, 0.11), "y", 0.5, thickness_um=0.22),
           ph.Port("out", (0.0, 15.0, 0.11), "y", 0.5, thickness_um=0.22)],
    source="in", wlens_um=(1.5, 1.55, 1.6),
    domain=ph.Domain(), mesh=ph.Mesh(cells_per_wlen=20))
data = ph.run_local(sim)

The map moves with the structures when the domain is fitted, so you draw and import in the same frame.

Ports may sit inside a heated region. The port mode is solved on the perturbed cross-section at the port plane, so a port in a uniformly heated guide launches and reads the heated mode. This is checked end to end (see below).

Inside a doped region, the port mode is the mode of the index-shifted core. Port modes are lossless, so the carrier loss is not part of it. That mode is checked, but a transmission read at a port inside a doped region has not yet been checked against a reference run. Keep ports outside doped regions when the number matters.

A phase shift is a difference of two runs. Run the heater off and on, on the same mesh, and compare the transmitted phase or the transmission. Most of the mesh error is common to the two runs and cancels.

Limits

  • CPU only. GPU runs refuse these media before they start.
  • Subpixel averaging is off. A simulation with these media turns subpixel averaging off, so material interfaces follow the mesh in steps. Absolute effective indices carry that error; a difference between two runs on one mesh removes most of it. Refine the mesh and compare before you rely on a value.
  • Interpolation. The solver interpolates the permittivity and loss linearly between map nodes, not the temperature or density. On a map finer than its own variation the difference is negligible.
  • First order. The thermo-optic model is linear in temperature, and the carrier models are fits for silicon near 1550 nm only, over the density span of the fit.
  • The base medium must be non-dispersive, isotropic and not PEC.
  • No heat or carrier solver. Compute the map elsewhere, or write it as a formula.
  • Reloading a result. A result saved by 0.1.6 with these media has no client state, so ph.RunResult(path) reads it without its declared ports and in the solver's frame. Runs on this SDK reload in full.

How it was checked

These checks are part of the solver's validation suite. The first three use a plane wave at normal incidence on a 0.5 um slab, uniform mesh, subpixel averaging off, non-dispersive media and 12 PML layers.

  • A uniformly heated slab (n = 2.0, 200 K, dn/dT = 2e-4): the transmitted phase change against the closed-form slab. The error falls from 3.8e-4 to 4.4e-5 to 1.7e-5 rad as the mesh goes from 25 to 8.3 to 5 nm, on a phase change of about 0.08 rad.
  • A linear temperature ramp on a non-uniform map not aligned to the mesh, against a finely sliced layered-medium reference. The reference uses the same rule the solver applies between map nodes (permittivity linear in position), not the physical ramp (index linear in position). Compared with the physical ramp, that rule alone moves the reflectance by 8.7e-5 and the phase by 1.6e-4 rad. Against the reference, the reflectance error is 4.7e-3, 5.6e-4, 2.2e-4 and 1.2e-4, and the phase-change error 4.8e-4, 6.2e-5, 2.4e-5 and 1.4e-5 rad, at 25, 8.3, 5 and 3.6 nm.
  • A silicon slab with 4e19 cm^-3 electrons and 7e19 cm^-3 holes, inside the default fit's density span: transmittance against the closed form of the complex index. The error is 2.91e-2, 3.27e-3, 1.21e-3 and 6.3e-4 at 25, 8.3, 5 and 3.6 nm. It tracks the mesh error of the same slab without carriers (3.23e-2, 3.64e-3, 1.33e-3, 7.0e-4). The absorption alone, the ratio of the two slabs on one mesh, is off by 2.45e-3, 2.59e-4, 1.13e-4 and 5.1e-5. This checks how the solver carries the loss map, not the carrier model itself.
  • A 500 x 220 nm silicon strip heated by 200 K, ports inside the heated guide, at 40 and 30 nm meshes. The port effective index equals that of an independent solve of a strip whose index is uniformly shifted. A second mode solver in the SDK, with a node-collocated discretization, agrees on the index change within 0.5% at its 10 and 6.7 nm meshes (up to 0.9% off at 20 nm). Both solvers read the same painted cross-section, so this checks the mode solvers, not the painting. The heated strip's transmission matches the unheated strip's to 8.5e-5 (40 nm) and 5.8e-5 (30 nm), and its reflection stays at the unheated strip's level. Neither strip reads a transmission of exactly 1: that offset is the launch and readout floor of these meshes, and it grows as the mesh is refined, so the check is the match between the two strips.