import numpy as np
import tidy3d as td
from matplotlib import pyplot as plt
from tidy3d import web
TIDY3D
LEARNING CENTER
Example Library
Example Library
/
Charge |
|
/
SISCAP modulator: from charge simulation to phase efficiency
SISCAP modulator: from charge simulation to phase efficiency¶
A SISCAP (Silicon-Insulator-Silicon CAPacitor) modulator is a MOS capacitor built inside a waveguide. A crystalline silicon layer, a thin gate oxide, and a poly-silicon layer overlap so that the gate oxide sits at the peak of the optical mode. When the capacitor is driven into accumulation, carrier sheets pile up on both sides of the gate oxide, exactly where the optical field is strongest. The overlap between the carrier sheets and the mode is therefore much larger than in a depletion pn phase shifter. The device in [1] reports a phase efficiency \(V_\pi L\) below 2 \(\text{V}\cdot\text{mm}\) at 1310 nm. That paper does not disclose its layer thicknesses or doping, so what follows is a SISCAP cross section of our own, simulated at 1.55 \(\mu\text{m}\), and the comparison with [1] stays qualitative.
This example follows the same two-stage workflow as the Charge solver example, applied to a SISCAP cross section:
A drift-diffusion simulation of the accumulation bias sweep, which gives the capacitance \(C(V)\) and the carrier distributions.
The carrier distributions are coupled into an optical mode solve through a perturbation medium, which gives the effective index change, the phase efficiency \(V_\pi L\), and the free carrier loss.
The charge solver ingredients (multi-physics media, drift-diffusion model, doping, boundary conditions, monitors) are described in detail in the Charge solver example, so the commentary here is short. The part that is specific to a MOS capacitor is the mesh: the gate oxide is 5 nm thick and the accumulation layers are about 1 nm thick, so the grid has to be roughly two orders of magnitude finer than for a pn junction.
The figure below shows the cross section simulated here: the two silicon layers overlapping across the gate oxide, each one with its own contact. The oxide cladding is not drawn and the contacts are sketched smaller than they are in the simulation.
References
[1] M. Webster, P. Gothoskar, V. Patel, D. Piede, S. Anderson, R. Tummidi, D. Adams, C. Appel, P. Metz, S. Sunder, B. Dama, and K. Shastri, “An efficient MOS-capacitor based silicon modulator and CMOS drivers for optical transmitters,” 11th International Conference on Group IV Photonics, 2014, DOI: 10.1109/Group4.2014.6961998.
[2] M. Nedeljkovic, R. Soref, and G. Z. Mashanovich, “Free-carrier electrorefraction and electroabsorption modulation predictions for silicon over the 1-14 \(\mu\text{m}\) infrared wavelength range,” IEEE Photonics Journal, vol. 3, no. 6, pp. 1171-1180, 2011, DOI: 10.1109/JPHOT.2011.2171930.
Charge Simulation
Cross Section
The waveguide is formed by the overlap of a p-type crystalline silicon layer and an n-type poly-silicon layer, separated by the gate oxide. Each layer extends laterally to its own contact, so the two layers are the two terminals of the capacitor. All lengths are in microns. As in the Charge solver example, the problem is solved in the \(xy\) plane, with \(y\) along the layer stack.
# all units in um
w_core = 0.5 # overlap width of the two silicon layers
h_csi = 0.11 # crystalline silicon (p-type) thickness
gap = 0.005 # gate oxide thickness
h_psi = 0.11 # poly-silicon (n-type) thickness
x_si = 3.0 # lateral extent of each layer, long enough to keep the metal off the mode
x_heavy = 1.5 # start of the heavily doped contact regions
w_contact = 1.0
h_contact = 0.5
y_gate = h_csi + gap / 2 # center of the gate oxide
h_stack = h_csi + gap + h_psi # top of the layer stack
n_body = 5e17 # doping of the two active layers
n_heavy = 1e19 # doping of the contact regionsMedia
Silicon is taken from the material library as a td.SemiconductorMedium valid at 300 K, the cladding is an insulator, and the contacts are auxiliary conductors used to place the boundary conditions. Every medium is wrapped in a td.MultiPhysicsMedium, so the same object carries the optical and the charge description.
# gate oxide and cladding: an insulator for the charge solve, the optical cladding later on
SiO2 = td.MultiPhysicsMedium(
optical=td.material_library["SiO2"]["Palik_LowLoss"],
charge=td.ChargeInsulatorMedium(permittivity=3.9),
name="SiO2",
)
# the contacts only carry the voltage boundary conditions, so any good conductor will do
aux = td.MultiPhysicsMedium(charge=td.ChargeConductorMedium(conductivity=1), name="aux")
# background medium of the simulation domain
air = td.MultiPhysicsMedium(heat=td.FluidSpec(), name="air")
# silicon drift-diffusion model at 300 K: mobility, recombination, bandgap narrowing
intrinsic_si = td.material_library["cSi"].variants["Si_MultiPhysics"].medium.chargeDoping
Both silicon layers are uniformly doped, so the doping is described with td.ConstantDoping boxes. The boxes are split at the center of the gate oxide, so acceptors are only added below it and donors only above it. Doping boxes add up, so the contact regions end up at n_body + n_heavy.
# the boxes are cut at the center of the gate oxide, so each layer only sees its own dopant
acceptors = [
# p-type crystalline silicon layer
td.ConstantDoping.from_bounds(
rmin=[-w_core / 2, -0.5, -td.inf], rmax=[x_si, y_gate, td.inf], concentration=n_body
),
# p++ contact region
td.ConstantDoping.from_bounds(
rmin=[x_heavy, -0.5, -td.inf], rmax=[x_si, y_gate, td.inf], concentration=n_heavy
),
]
donors = [
# n-type poly-silicon layer
td.ConstantDoping.from_bounds(
rmin=[-x_si, y_gate, -td.inf], rmax=[w_core / 2, 0.5, td.inf], concentration=n_body
),
# n++ contact region
td.ConstantDoping.from_bounds(
rmin=[-x_si, y_gate, -td.inf], rmax=[-x_heavy, 0.5, td.inf], concentration=n_heavy
),
]
si = td.MultiPhysicsMedium(
charge=intrinsic_si.updated_copy(N_a=acceptors, N_d=donors), name="Si"
)Structures
Each layer is a td.Structure: the two silicon layers overlap over w_core and each one runs out to a contact on its own side.
# one big oxide box: it fills the gate gap and the cladding in a single structure
oxide = td.Structure(
geometry=td.Box(center=(0, y_gate, 0), size=(12, 8, td.inf)), medium=SiO2, name="oxide"
)
si_bottom = td.Structure(
geometry=td.Box.from_bounds(rmin=(-w_core / 2, 0, -td.inf), rmax=(x_si, h_csi, td.inf)),
medium=si,
name="si_bottom",
)
si_top = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-x_si, h_csi + gap, -td.inf), rmax=(w_core / 2, h_stack, td.inf)
),
medium=si,
name="si_top",
)
# the contacts sit on top of each layer, touching but not overlapping it, and the voltage
# boundary condition is later placed on the contact boundary
contact_p = td.Structure(
geometry=td.Box.from_bounds(
rmin=(x_si - w_contact, h_csi, -td.inf), rmax=(x_si, h_csi + h_contact, td.inf)
),
medium=aux,
name="contact_p",
)
contact_n = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-x_si, h_stack, -td.inf), rmax=(-x_si + w_contact, h_stack + h_contact, td.inf)
),
medium=aux,
name="contact_n",
)
all_structures = [oxide, si_bottom, si_top, contact_p, contact_n]
scene = td.Scene(medium=air, structures=all_structures)
_, ax = plt.subplots(2, 1, figsize=(9, 3.4))
scene.plot(z=0, ax=ax[0])
ax[0].set_ylim(-0.5, 1.0)
ax[0].set_xlabel("")
scene.plot_structures_property(z=0, property="doping", ax=ax[1], limits=[-1e18, 1e18])
ax[1].set_xlim(-1.0, 1.0)
ax[1].set_ylim(-0.05, 0.28)
plt.tight_layout()
plt.show()
Boundary Conditions and Monitors
The gate, which is the n-type poly-silicon layer, is grounded, and the crystalline silicon contact is swept from 0 V to 2.4 V, both through a td.VoltageBC. This polarity attracts holes to the bottom side of the gate oxide and electrons to the top side, so both interfaces accumulate at the same time. The 0.2 V sweep step also sets the \(\Delta V\) used by the capacitance monitor.
Three more monitors record the solution itself: a td.SteadyPotentialMonitor and a td.SteadyFreeCarrierMonitor on a window around the gate, and a second free carrier monitor over the whole cross section, which the optical stage needs.
# one solve per entry; the 0.2 V spacing is also the dV the capacitance monitor differentiates
voltages = list(np.round(np.linspace(0, 2.4, 13), 3))
boundary_conditions = [
# swept terminal: the crystalline silicon layer
td.HeatChargeBoundarySpec(
condition=td.VoltageBC(source=td.DCVoltageSource(voltage=voltages)),
placement=td.StructureBoundary(structure=contact_p.name),
),
# grounded terminal: the poly-silicon gate
td.HeatChargeBoundarySpec(
condition=td.VoltageBC(source=td.DCVoltageSource(voltage=0)),
placement=td.StructureBoundary(structure=contact_n.name),
),
]
# small signal capacitance of the whole cross section, per unit length
capacitance_mnt = td.SteadyCapacitanceMonitor(
center=(0, y_gate, 0), size=(td.inf, td.inf, 0), name="capacitance"
)
# potential and carriers on a window around the gate, returned on the solver mesh
potential_mnt = td.SteadyPotentialMonitor(
center=(0, y_gate, 0), size=(1.0, 0.12, 0), name="potential_gate", unstructured=True
)
carrier_gate_mnt = td.SteadyFreeCarrierMonitor(
center=(0, y_gate, 0), size=(1.0, 0.12, 0), name="carriers_gate", unstructured=True
)
# full cross section, used later to perturb the optical medium
carrier_mnt = td.SteadyFreeCarrierMonitor(
center=(0, y_gate, 0), size=(td.inf, td.inf, 0), name="carriers", unstructured=True
)Mesh
This is the one part of the setup that differs substantially from a pn junction. Two length scales have to be resolved: the 5 nm gate oxide and the accumulation layers, which are one to two nanometers thick at these bias values. A td.GridRefinementRegion with dl_internal of 0.8 nm covers a thin band around the gate oxide, a second region covers the rest of the layer stack more coarsely, and the surrounding td.DistanceUnstructuredGrid relaxes to 150 nm out towards the contacts. Setting relative_min_dl = 0 is required; otherwise the element size is clamped relative to the size of the simulation domain and the sub-nanometer request is ignored.
# The band around the gate. Its height, 20 nm, is the 5 nm oxide plus roughly 8 nm of silicon on
# each side, which is a little more than one Debye length (about 6 nm at 5e17 cm^-3) and therefore
# covers the accumulation layers. Its width runs 75 nm past each edge of the overlap, where the
# fringing field bends around the corners of the two layers. The element size is a good fraction
# below the Debye length; 0.8 nm here also puts 6 elements across the oxide itself.
gate_refinement = td.GridRefinementRegion(
center=(0, y_gate, 0),
size=(w_core + 0.15, 0.02, 0),
dl_internal=0.0008,
transition_thickness=0.05, # grade back to the surrounding element size over 50 nm
)
# Outside that band the carriers only vary on the doping scale, so 10 nm is plenty. The region
# covers the whole 225 nm stack plus a small margin, and the slab on both sides of the overlap.
stack_refinement = td.GridRefinementRegion(
center=(0, y_gate, 0),
size=(w_core + 0.6, 0.28, 0),
dl_internal=0.01,
transition_thickness=0.2,
)
mesh = td.DistanceUnstructuredGrid(
dl_interface=0.01, # boundary layers form at every semiconductor-to-insulator interface
dl_bulk=0.15, # far from the junction nothing varies, so elements can be 15x larger
distance_interface=0.02,
distance_bulk=0.6,
relative_min_dl=0, # without this the element size is clamped relative to the domain size
sampling=800,
non_refined_structures=[oxide.name], # no need to mesh the cladding finely
mesh_refinements=[gate_refinement, stack_refinement],
)The HeatChargeSimulation Object
The td.HeatChargeSimulation collects the structures, the monitors, the boundary conditions and the mesh, and td.IsothermalSteadyChargeDCAnalysis selects the isothermal steady state DC analysis. Fermi statistics are turned on with fermi_dirac=True because the accumulated carrier densities exceed \(10^{20}\ \text{cm}^{-3}\), where the Boltzmann approximation overestimates the charge. Nothing else needs tuning. In particular convergence_dv caps the voltage step the solver takes internally between two solutions, so it only does something when consecutive bias points are far apart, which is not the case here: our sweep already advances in 0.2 V steps. Forcing 0.05 V internal steps together with a larger Newton iteration budget made this solve 2.4 times slower while moving the capacitance by less than 0.1 percent.
charge_sim = td.HeatChargeSimulation(
center=(0, y_gate, 0),
size=(2 * x_si + 1.0, 3.0, 0),
structures=all_structures,
medium=air,
monitors=[capacitance_mnt, potential_mnt, carrier_gate_mnt, carrier_mnt],
boundary_spec=boundary_conditions,
grid_spec=mesh,
# Fermi statistics matter here; everything else can stay at its default
analysis_spec=td.IsothermalSteadyChargeDCAnalysis(temperature=300, fermi_dirac=True),
)
ax = charge_sim.plot(z=0)
ax.set_ylim(-0.6, 1.0)
plt.show()
Mesh Inspection
Meshing runs as a separate task through td.VolumeMesher, so the mesh is worth checking before the solve: everything downstream depends on how well the gate region is resolved.
# meshing is a separate, cheap task, so the mesh can be inspected before paying for the solve
mesher = td.VolumeMesher(
simulation=charge_sim,
monitors=[td.VolumeMeshMonitor(size=(td.inf, td.inf, 0), name="mesh")],
)
mesher_job = web.Job(simulation=mesher, task_name="siscap_mesh")
mesher_data = mesher_job.run(path="siscap_mesh.hdf5")22:31:44 UTC Created task 'siscap_mesh' with resource_id
'vom-433ad2f9-996b-4eb7-b85e-5409341bf4fd' and task_type
'VOLUME_MESH'.
Tidy3D's VolumeMesher solver is currently in the beta stage. Cost
of VolumeMesher simulations is subject to change in the future.
↑ simulation.hdf5.gz ━━━━━━━━━━━━━━━━━━━━━━━━━ 100.0% • 2.2/2.2 kB • ? • 0:00:00
22:31:46 UTC Estimated FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to
get the billed FlexCredit cost after a simulation run.
22:31:50 UTC status = queued
To cancel the simulation, use 'web.abort(task_id)' or
'web.delete(task_id)' or abort/delete the task in the web UI.
Terminating the Python script will not stop the job running on the
cloud.
22:31:55 UTC starting up solver
running solver
22:32:05 UTC status = success
↓ simulation_data.hdf5.gz ━━━━━━━━━━━━━ 100.0% • 1.4/1.4 MB • 2.4 MB/s • 0:00:00
22:32:06 UTC Loading results from siscap_mesh.hdf5 _, ax = plt.subplots(1, 3, figsize=(15, 4))
mesher_data.plot_mesh("mesh", z=0, ax=ax[0], structures_fill=False)
ax[0].set_ylim(-0.6, 1.0)
mesher_data.plot_mesh("mesh", z=0, ax=ax[1], structures_fill=True)
ax[1].set_xlim(-0.5, 0.5)
ax[1].set_ylim(0.0, 0.24)
mesher_data.plot_mesh("mesh", z=0, ax=ax[2], structures_fill=True)
ax[2].set_xlim(-0.03, 0.03)
ax[2].set_ylim(0.104, 0.122)
ax[2].set_title("gate oxide")
plt.tight_layout()
plt.show()
Run the Charge Simulation
The meshing task is passed as a parent task so that it is not repeated by the solver run. Job.run is used rather than Job.step so that a second execution of the notebook restores the result from the local cache instead of failing on an already completed workflow.
job = web.Job(
simulation=charge_sim, task_name="siscap_charge", parent_tasks=[mesher_job.task_id]
)
job.estimate_cost()22:32:08 UTC Created task 'siscap_charge_solve' with resource_id
'hec-aa1eb82a-90a6-4dc8-b81e-0d044bcc5b1b' and task_type
'HEAT_CHARGE'.
Tidy3D's HeatCharge solver is currently in the beta stage. Cost of
HeatCharge simulations is subject to change in the future.
↑ simulation.hdf5.gz ━━━━━━━━━━━━━━━━━━━━━━━━━ 100.0% • 2.1/2.1 kB • ? • 0:00:00
22:32:10 UTC Estimated typical FlexCredit cost: 0.032. For charge simulations,
the billed cost depends on the number of solver iterations required
for convergence.
Maximum FlexCredit cost: 5.458. This assumes the charge solver
reaches its configured iteration limits for all applied biases. Use
'web.real_cost(task_id)' to get the billed FlexCredit cost after a
simulation run.
The FlexCredit estimate shown above is for the next workflow step
'solve' only. 5.457620872309933charge_data = job.run(path="siscap_charge.hdf5")22:32:11 UTC status = queued
To cancel the simulation, use 'web.abort(task_id)' or
'web.delete(task_id)' or abort/delete the task in the web UI.
Terminating the Python script will not stop the job running on the
cloud.
22:32:16 UTC status = preprocess
22:32:18 UTC status = queued
22:32:30 UTC status = preprocess
22:32:48 UTC starting up solver
running solver
22:34:46 UTC status = queued
22:35:33 UTC status = preprocess
22:35:56 UTC status = running
22:40:39 UTC status = success
↓ simulation_data.hdf5.gz ━━━━━━━━━━━━━ 100.0% • 8.9/8.9 MB • 9.2 MB/s • 0:00:00
22:40:44 UTC Loading results from siscap_charge.hdf5
WARNING: Warning messages were found in the solver log. For more
information, check 'SimulationData.log' or use
'web.download_log(task_id)'. The solver log reports that a couple of the thirteen bias points did not converge. That means the Newton iteration at those biases stopped on its iteration cap with the residual still above the requested tolerance, not that the solve failed. dc_convergence says which ones:
convergence = charge_data.device_characteristics.dc_convergence
for v, ok, iters in zip(
convergence.converged.coords["v"].data, convergence.converged.data, convergence.n_iters.data
):
print(f"{v:5.2f} V converged: {bool(ok)!s:5} Newton iterations: {int(iters)}") 0.00 V converged: True Newton iterations: 34
0.20 V converged: True Newton iterations: 13
0.40 V converged: False Newton iterations: 60
0.60 V converged: True Newton iterations: 56
0.80 V converged: True Newton iterations: 18
1.00 V converged: True Newton iterations: 18
1.20 V converged: True Newton iterations: 19
1.40 V converged: True Newton iterations: 23
1.60 V converged: True Newton iterations: 23
1.80 V converged: True Newton iterations: 33
2.00 V converged: True Newton iterations: 32
2.20 V converged: True Newton iterations: 28
2.40 V converged: True Newton iterations: 43The two that stall sit at 0.4 V and 0.6 V, on the steepest part of the C-V curve, where the device crosses from depletion into accumulation. That is the hardest region for a MOS capacitor: the surface potential swings fastest with applied bias there, and the carrier densities in the depleted sheets span twenty orders of magnitude, so the last decades of the residual are expensive to buy. Everywhere else the solver converges in 13 to 48 iterations.
It is worth checking whether that matters rather than assuming, and here it does not. The C-V curve and the carrier profiles are smooth through both points, the carrier fields are physical, and buying the missing residual with internal ramp steps and a larger iteration budget moves the capacitance by less than 0.1 percent, as mentioned above. A stalled point matters when the curve kinks or jumps at it, when the carrier densities collapse towards intrinsic values, or when the answer moves once you refine the solver settings. If any of that happens, refine the bias sweep or the mesh rather than trusting the point.
Potential and Carriers
At zero bias the built-in potential of the two oppositely doped layers already drops across the gate oxide. The applied bias adds to it, and at 2.4 V nearly the whole 2.4 V sits across the 5 nm oxide.
_, ax = plt.subplots(1, 2, figsize=(12, 2.6))
for a, v in zip(ax, [voltages[0], voltages[-1]]):
charge_data["potential_gate"].potential.sel(voltage=v).plot(ax=a, grid=False)
a.set_title(f"potential (V), bias {v} V")
plt.tight_layout()
plt.show()
A vertical cut through the middle of the capacitor shows the accumulation layers directly. At 0 V both interfaces are depleted. Around 0.8 V the bands are close to flat. Above that, holes build up in the crystalline silicon just below the oxide and electrons in the poly-silicon just above it, reaching a few times \(10^{19}\ \text{cm}^{-3}\) and decaying over a few nanometers. The interpolated values inside the oxide are not meaningful and are masked out.
# 30 nm of silicon on each side of the oxide, about five Debye lengths
y_cut = np.linspace(h_csi - 0.03, h_csi + gap + 0.03, 601)
in_oxide = (y_cut > h_csi) & (y_cut < h_csi + gap)
_, ax = plt.subplots(1, 2, figsize=(12, 4))
for v in [0.0, 0.8, 1.6, 2.4]:
# interpolate the unstructured solution onto the vertical line x = 0
holes = charge_data["carriers"].holes.sel(voltage=v).interp(x=[0], y=y_cut, z=[0])
electrons = charge_data["carriers"].electrons.sel(voltage=v).interp(x=[0], y=y_cut, z=[0])
for a, data in zip(ax, [holes, electrons]):
profile = np.array(data.values).ravel()
profile[in_oxide] = np.nan
a.semilogy(y_cut * 1e3, profile, label=f"{v} V")
for a, name in zip(ax, ["holes", "electrons"]):
a.axvspan(h_csi * 1e3, (h_csi + gap) * 1e3, color="0.85")
a.set_xlabel("y (nm)")
a.set_ylabel(f"{name} (cm$^{{-3}}$)")
a.set_title(f"{name}, vertical cut at x = 0, gate oxide shaded")
a.set_ylim(1e13, 1e20)
a.grid(which="both", alpha=0.3)
a.legend()
plt.tight_layout()
plt.show()
Capacitance
The td.SteadyCapacitanceMonitor returns the small-signal capacitance per unit length in fF/\(\mu\text{m}\), computed separately from the electron and the hole charge. The natural reference for a MOS capacitor is the geometric oxide capacitance \(C_{\text{ox}} = \varepsilon_{\text{ox}} w_{\text{core}} / t_{\text{ox}}\), which the device approaches once both interfaces are accumulated. The measured C-V of a SISCAP device has the same shape, see Fig. 1(c) of [1]: a shallow minimum where the two interfaces are depleted, then a steep rise into accumulation. That steep rise above about 1 V is also why [1] drives the device between 1.2 V and 2.2 V rather than between 0 V and 1 V.
capacitance = charge_data["capacitance"]
v_cap = np.array(capacitance.electron_capacitance.coords["v"].data)
c_e = np.array(capacitance.electron_capacitance.data)
c_h = np.array(capacitance.hole_capacitance.data)
# electrons and holes give the same capacitance, so averaging them is a free consistency check
c_fF_um = np.abs(0.5 * (c_e + c_h))
# td.EPSILON_0 is given in F/um, hence the 1e15 to get fF/um
c_ox = 3.9 * td.EPSILON_0 * w_core / gap * 1e15
plt.figure(figsize=(6, 4))
plt.plot(v_cap, c_fF_um, "k.-", label="simulation")
plt.axhline(c_ox, ls="--", c="tab:red", label=f"$C_{{ox}}$ = {c_ox:.2f} fF/$\\mu$m")
plt.xlabel("bias (V)")
plt.ylabel("capacitance (fF/$\\mu$m)")
plt.ylim(0, 1.15 * c_ox)
plt.legend()
plt.grid()
plt.show()
print(f"C at {v_cap[0]:.1f} V: {c_fF_um[0]:.3f} fF/um = {10 * c_fF_um[0]:.1f} pF/cm")
print(f"C at {v_cap[-1]:.1f} V: {c_fF_um[-1]:.3f} fF/um = {10 * c_fF_um[-1]:.1f} pF/cm")
print(f"C / C_ox at {v_cap[-1]:.1f} V: {c_fF_um[-1] / c_ox:.3f}")
C at 0.0 V: 0.993 fF/um = 9.9 pF/cm
C at 2.4 V: 3.384 fF/um = 33.8 pF/cm
C / C_ox at 2.4 V: 0.980Optical Response
The carrier distributions are now used to perturb the refractive index of silicon with the free-carrier dispersion model of [2], available as td.NedeljkovicSorefMashanovich. The perturbation, carried by a td.PerturbationMedium, is applied to the absolute carrier concentrations, so the loss at zero bias already contains the free-carrier absorption of the background doping.
wvl_um = 1.55
freq0 = td.C_0 / wvl_um
# start from lossless silicon so that all the loss we report comes from the free carriers
n_si = td.material_library["cSi"]["Palik_NoLoss"].nk_model(frequency=freq0)[0]
perturbation_model = td.NedeljkovicSorefMashanovich(ref_freq=freq0)
# a medium whose index and absorption follow whatever carrier density we hand it later
si_perturb = td.PerturbationMedium.from_unperturbed(
medium=td.Medium.from_nk(n=n_si, k=0, freq=freq0),
perturbation_spec=td.IndexPerturbation(
delta_n=perturbation_model.delta_n(),
delta_k=perturbation_model.delta_k(),
freq=freq0,
),
)Optical Simulation
The optical simulation reuses the two silicon layers with the perturbation medium, in an oxide background, at a wavelength of 1.55 \(\mu\text{m}\). The mode grid needs one refinement, a td.MeshOverrideStructure with a 0.1 nm step across a thin band around the gate oxide. Without it the 5 nm oxide is averaged away and the effective index is wrong in the second decimal. Refining the rest of the 110 nm layers beyond the automatic grid is not worth it: it shifts the effective index in the third decimal and leaves the index change with bias, which is what the modulator response depends on, essentially untouched.
# 0.1 nm across the gate band. A 1 nm step overestimates the index change by about 5 percent,
# while 0.1 nm and 0.025 nm agree to 0.1 percent, so 0.1 nm is where this converges.
gate_grid = td.MeshOverrideStructure(
geometry=td.Box(center=(0, y_gate, 0), size=(1.4, 0.04, td.inf)), dl=(0.01, 0.0001, 0.01)
)
# this simulation is only a container for the mode solver, so the run time and the boundaries are
# nominal; the grid and the media are what matter
opt_sim = td.Simulation(
center=(0, y_gate, 0),
size=(4.0, 3.0, 0.4),
medium=SiO2.optical,
structures=[
si_bottom.updated_copy(medium=si_perturb),
si_top.updated_copy(medium=si_perturb),
],
run_time=1e-12,
boundary_spec=td.BoundarySpec.all_sides(td.PML()),
grid_spec=td.GridSpec.auto(
wavelength=wvl_um,
min_steps_per_wvl=25,
override_structures=[gate_grid],
),
)
# the plane is a few wavelengths across, so the evanescent tails decay well inside it
mode_plane = td.Box(center=(0, y_gate, 0), size=(2.6, 2.0, 0))
# target_neff pins the search on the guided mode; double precision because the index change we
# are looking for is in the fourth decimal
mode_spec = td.ModeSpec(num_modes=1, precision="double", target_neff=n_si)
ax = opt_sim.plot(z=0)
mode_plane.plot(z=0, ax=ax, alpha=0.3)
ax.set_ylim(-0.6, 0.9)
plt.show()
The Optical Mode
The fundamental TE mode is centered on the overlap region and peaks inside the gate oxide. This is the property that makes the SISCAP structure efficient: the accumulation layers sit at the field maximum.
# perturbed_mediums_copy turns the carrier distribution into a custom optical medium
mode_sim_0 = td.ModeSimulation.from_simulation(
simulation=opt_sim.perturbed_mediums_copy(
electron_density=charge_data["carriers"].electrons.sel(voltage=0.0),
hole_density=charge_data["carriers"].holes.sel(voltage=0.0),
),
wavelength=wvl_um,
plane=mode_plane,
freqs=[freq0],
mode_spec=mode_spec,
)
mode_data_0 = mode_sim_0.run_local()field = np.abs(mode_data_0.modes.Ex.isel(mode_index=0, f=0))
_, ax = plt.subplots(1, 2, figsize=(11, 4))
field.plot(x="x", y="y", ax=ax[0], cmap="magma")
ax[0].set_xlim(-0.8, 0.8)
ax[0].set_ylim(y_gate - 0.4, y_gate + 0.4)
ax[0].set_title("$|E_x|$")
field.sel(x=0, method="nearest").plot(ax=ax[1])
ax[1].axvspan(h_csi, h_csi + gap, color="0.85")
ax[1].set_xlim(-0.05, 0.28)
ax[1].set_xlabel("y ($\\mu$m)")
ax[1].set_title("$|E_x|$ at x = 0, gate oxide shaded")
ax[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
print(f"n_eff = {float(mode_data_0.modes.n_eff.isel(mode_index=0, f=0)):.4f}")
n_eff = 2.5770Mode Solve at Every Bias
One td.ModeSimulation is built per bias point with perturbed_mediums_copy and solved with run_local: thirteen single-frequency mode solves run faster on this machine than the round trip to the cloud would take. Field output is switched off with fields=[], since only the complex effective index is needed.
mode_data = {}
for v in voltages:
# same waveguide at every bias, only the carrier distribution changes
perturbed_sim = opt_sim.perturbed_mediums_copy(
electron_density=charge_data["carriers"].electrons.sel(voltage=v),
hole_density=charge_data["carriers"].holes.sel(voltage=v),
)
mode_sim = td.ModeSimulation.from_simulation(
simulation=perturbed_sim,
wavelength=wvl_um,
plane=mode_plane,
freqs=[freq0],
mode_spec=mode_spec,
fields=[],
)
mode_data[f"{v:.1f}V"] = mode_sim.run_local()Phase Efficiency and Loss
From the complex effective index we get the phase shift, the phase efficiency, and the loss:
\[\Delta \phi = \frac{2 \pi \Delta n_{\text{eff}} L}{\lambda_0}, \qquad V_\pi L = \frac{\lambda_0}{2\, |dn_{\text{eff}} / dV|}, \qquad \alpha_{\text{dB/cm}} = \frac{40 \pi \, \mathrm{Im}(n_{\text{eff}})}{\lambda_0} \times 10^4 \log_{10} e\]
The accumulated carriers lower the index, so the phase shift is negative. Its magnitude is plotted below. The efficiency improves steeply as the capacitor enters accumulation, which is the behavior reported in Fig. 1(d) of [1], while the loss grows with the accumulated carrier density.
n_complex = np.array(
[complex(mode_data[f"{v:.1f}V"].modes.n_complex.sel(mode_index=0, f=freq0).values)
for v in voltages]
)
voltages_arr = np.array(voltages)
delta_neff = np.real(n_complex - n_complex[0])
phase_deg_mm = -360 * delta_neff / wvl_um * 1e3 # degrees per mm
# local efficiency, from the slope of n_eff against bias at each point
vpi_l = wvl_um * 1e-3 / (2 * np.abs(np.gradient(delta_neff, voltages_arr))) # V.mm
alpha_dB_cm = 10 * 4 * np.pi * np.imag(n_complex) / wvl_um * 1e4 * np.log10(np.e)
_, ax = plt.subplots(1, 3, figsize=(15, 4))
ax[0].plot(voltages_arr, phase_deg_mm, "k.-")
ax[0].set_ylabel("phase shift (deg/mm)")
ax[1].plot(voltages_arr, vpi_l, "k.-")
ax[1].set_ylabel("$V_\\pi L$ (V$\\cdot$mm)")
ax[1].set_ylim(0, min(12, 1.1 * np.max(vpi_l)))
ax[2].plot(voltages_arr, alpha_dB_cm / 10, "k.-")
ax[2].set_ylabel("loss (dB/mm)")
for a in ax:
a.set_xlabel("bias (V)")
a.grid()
plt.tight_layout()
plt.show()
Figures of Merit Over a 1 V Drive in Accumulation
The device in [1] is driven between 1.2 V and 2.2 V rather than between 0 V and 1 V, for the reason the C-V curve makes obvious. Averaged over that swing, this cross section gives the following. The loss is dominated by the accumulated carriers at the top of the swing.
i_lo = int(np.argmin(np.abs(voltages_arr - 1.2)))
i_hi = int(np.argmin(np.abs(voltages_arr - 2.2)))
# efficiency averaged over the swing, rather than the local slope plotted above
dneff_swing = delta_neff[i_hi] - delta_neff[i_lo]
dv_swing = voltages_arr[i_hi] - voltages_arr[i_lo]
vpi_l_swing = wvl_um * 1e-3 * dv_swing / (2 * abs(dneff_swing))
print(f"delta n_eff over 1.2 V to 2.2 V : {dneff_swing:.3e}")
print(f"phase shift : {abs(dneff_swing) / wvl_um * 360e3:.0f} deg/mm")
print(f"V_pi L : {vpi_l_swing:.2f} V.mm")
print(f"loss at 1.2 V : {alpha_dB_cm[i_lo] / 10:.2f} dB/mm")
print(f"loss at 2.2 V : {alpha_dB_cm[i_hi] / 10:.2f} dB/mm")
print(f"capacitance at 1.2 V : {c_fF_um[i_lo]:.2f} fF/um")
print(f"capacitance at 2.2 V : {c_fF_um[i_hi]:.2f} fF/um")delta n_eff over 1.2 V to 2.2 V : -3.591e-04
phase shift : 83 deg/mm
V_pi L : 2.16 V.mm
loss at 1.2 V : 1.69 dB/mm
loss at 2.2 V : 3.43 dB/mm
capacitance at 1.2 V : 2.76 fF/um
capacitance at 2.2 V : 3.37 fF/umThe two halves of the workflow give the quantities a link budget needs. The phase efficiency and the loss set the arm length and the insertion loss of a Mach-Zehnder built from this phase shifter, and the capacitance sets the drive power and, together with the series resistance of the two silicon layers, the RC bandwidth.
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: