Charge solver: a silicon electro-optic phase modulator

In this notebook, we simulate a silicon electro-optic phase modulator. The numerical results are compared with the experimental data provided in [1]. The modulator consists of a pn-junction with a high doping contact to allow for a less resistive electrical interface.

The sketch below illustrates the main dimensions of the modulator:

Outline

1. We start by setting up a Charge simulation in tidy3d. For this example, the isothermal Drift-Diffusion (DD) equations are solved. A basic summary of the equations solved and nomenclature used can be found our td.SemiconductorMedium documentation.

2. We show how the numerical results from this simulation can be used to determine the optical phase change.

References

[1] Chrostowski L, Hochberg M. Silicon Photonics Design: From Devices to Systems. Cambridge University Press; 2015.

[2] Schroeder, D., T. Ostermann, and O. Kalz. “Comparison of transport models far the simulation of degenerate semiconductors.” Semiconductor science and technology 9.4 (1994): 364.

1. Construct our Charge Simulation

import numpy as np
import tidy3d as td
from matplotlib import pyplot as plt

Specify our Modulator Dimensions

Based on the figure above, we can define the dimensions of our electro-optic phase modulator. These are defaulted to micro-meters.

Since 2D problems are not yet supported, we will extrude our geometry in the z dimension. This is perpendicular to the plane of the problem in z_size \(\mu m\).

# all units in um
w_core = 0.25
h_core = 0.22

w_clearance = 2.25
h_clearance = 0.09
w_side = 2.5
h_side = 0.22

w_contact = 1.2
h_contact = 1.

z_size = h_clearance/5

res = h_clearance/10

Create Multi-Physical Mediums

In this example, we want to simulate the optical and electronic behaviour of this device. This means we need to specify material models that can represent the behaviour in these two physical domains.

The td.MultiPhysicsMedium was introduced in order to compose existing tidy3d medium models together into a unified multi-domain representation. td.MultiPhysics medium will be fully interoperable with any tidy3d simulation definition or class that uses existing medium definitions, and should enable full flexibility in terms of creating unified material models progressively. In this example, we will focus on its usage for Charge simulations which it fully supports.

We still need to define a <Domain>Medium model for our Charge and Optical FDTD simulations. Let’s explore the different ways of doing this below.

Our silica (\(SiO_2\)) cladding behaves like an electronic insulator which can be defined with a td.ChargeInsulatorMedium. For the optical simulation, assuming we have an td.Medium frequency range valid to 0 Hz (electronic DC), we can reuse a the td.material_library optical medium models and create a td.MultiPhysicsMedium accordingly:

SiO2_optic = td.material_library['SiO2']['Palik_Lossless']
SiO2 = td.MultiPhysicsMedium(
    optical=SiO2_optic,
    charge=td.ChargeInsulatorMedium(permittivity=SiO2_optic.eps_model(frequency=0)),
    name="SiO2",
)

11:31:24 CET WARNING: frequency passed to 'Medium.eps_model()'is outside of     
             'Medium.frequency_range' = (59958491600000.0, 1998616386666666.8)  

/home/momchil/miniconda3/envs/flex/lib/python3.12/site-packages/pydantic/v1/validators.py:158: ComplexWarning: Casting complex values to real discards the imaginary part
  return float(v)

However, in this example, we will create our own silica multi-physics model. The above td.material_library['SiO2']['Palik_Lossless'] is not well defined for frequency=0. Hence, it is outside of the simulation validity region.

SiO2 = td.MultiPhysicsMedium(
    optical=td.Medium(permittivity=3.9),
    charge=td.ChargeInsulatorMedium(permittivity=3.9), # redefining permittivity
    name="SiO2",
)

And now we define the remaining auxiliary mediums used to define our boundary conditions:

aux = td.MultiPhysicsMedium(
    charge=td.ChargeConductorMedium(conductivity=1),
    name="aux"
)

air = td.MultiPhysicsMedium(
    heat=td.FluidSpec(),
    name="air"
)

A Drift-Diffusion Model defined in a td.SemiconductorMedium

We will begin defining the electronic properties of the intrinsic semiconductor: i.e., the td.SemiconductorMedium without any doping.

Within a SemiconductorMedium, the Drift-Diffusion (DD) equations will be solved. In this example, we will explore an isothermal (constant temperature) case. This simulation performed at $ T=300 $ K. The DD equations, their assumptions, and limitations for this type of problem are discussed below:

The electrostatic potential \(\psi\), the electron \(n\) and hole \(p\) concentrations are computed from the following coupled system of PDEs.

\[-\nabla \cdot \left( \varepsilon_0 \varepsilon_r \nabla \psi \right) = q \left( p - n + N_d^+ - N_a^- \right)\]

\[q \frac{\partial n}{\partial t} = \nabla \cdot \mathbf{J_n} - qR\]

\[q \frac{\partial p}{\partial t} = -\nabla \cdot \mathbf{J_p} - qR\]

The above system requires the definition of the generation-recombination rate \(R\) represented by our generation-recombination models and the flux functions (free carrier current density), $ \mathbf{J_n} $ and $ \mathbf{J_p} $. Their usual form is:

\[\mathbf{J_n} = q \mu_n \mathbf{F_{n}} + q D_n \nabla n\]

\[\mathbf{J_p} = q \mu_p \mathbf{F_{p}} - q D_p \nabla p\]

These depend on the carrier mobility for electrons \(\mu_n\) and holes \(\mu_p\) represented by our mobility models.

Where the effective field, defined in [2], is simplified to:

\[\mathbf{F_{n,p}} = \nabla \psi\]

This approximation does not consider the effects of bandgap narrowing and degeneracy on the effective electric field $ \mathbf{F_{n,p}} $, which is acceptable for non-degenerate semiconductors. Bandgap narrowing models which are important for doped semiconductor effects can also be defined.

Material properties are defined as class parameters or other classes:

Symbol	Parameter Name	Description
\[N_a\]	N_a	Ionized acceptors density
\[N_d\]	N_d	Ionized donors density
\[N_c\]	N_c	Effective density of states in the conduction band
\[N_v\]	N_v	Effective density of states in valence band
\[R\]	R	Generation-Recombination term
\[E_g\]	E_g	Bandgap Energy
\[\Delta E_g\]	delta_E_g	Bandgap Narrowing
\[\sigma\]	conductivity	Electrical conductivity
\[\varepsilon_r\]	permittivity	Relative permittivity
\[\varepsilon_0\]	tidy3d.constants.EPSILON_0	Free Space permittivity
\[q\]	tidy3d.constants.Q_e	Fundamental electron charge

Our material library already has a silicon td.MultiPhysicsMedium model valid for a 300K isothermal case which incorporates a td.SemiconductorMedium:

intrinsic_si = td.material_library['cSi'].variants['Si_MultiPhysics'].medium.charge

In case we want to build our own td.SemiconductorMedium, we can do so in the following manner (note that is identical to the one above):

# Create a semiconductor medium with mobility, generation-recombination, and bandgap narrowing models.
intrinsic_si = td.SemiconductorMedium(
    permittivity=11.7,
    N_c=2.86e19,
    N_v=3.1e19,
    E_g=1.11,
    mobility_n=td.CaugheyThomasMobility(
        mu_min=52.2,
        mu=1471.0,
        ref_N=9.68e16,
        exp_N=0.68,
        exp_1=-0.57,
        exp_2=-2.33,
        exp_3=2.4,
        exp_4=-0.146,
    ),
    mobility_p=td.CaugheyThomasMobility(
        mu_min=44.9,
        mu=470.5,
        ref_N=2.23e17,
        exp_N=0.719,
        exp_1=-0.57,
        exp_2=-2.33,
        exp_3=2.4,
        exp_4=-0.146,
    ),
    R=[
        td.ShockleyReedHallRecombination(tau_n=3.3e-6, tau_p=4e-6),
        td.RadiativeRecombination(r_const=1.6e-14),
        td.AugerRecombination(c_n=2.8e-31, c_p=9.9e-32),
    ],
    delta_E_g=td.SlotboomBandGapNarrowing(
        v1=6.92 * 1e-3,
        n2=1.3e17,
        c2=0.5,
        min_N=1e15,
    ),
    N_a=0,
    N_d=0,
)

Create our doping regions

In our intrinsic silicon model above, we set the ionized acceptors \(N_a\) and ionized donor \(N_d\) density to zero. Since we will be doping our semiconductors with acceptors and donors, we need to create these doping distributions.

In this examples, we’ll use Gaussian doping boxes using the td.GaussianDoping model. The Gaussian doping concentration \(N\) is defined in the following manner:

• \(N = N_{\text{max}}\) at locations more than width \(\mu m\) away from the sides of the box.

• \(N = N_{\text{ref}}\) at locations on the box sides.

• A Gaussian variation between \(N_{\text{max}}\) and \(N_{\text{ref}}\) at locations less than width µm away from the sides.

By definition, all sides of the box will have a concentration \(N_{\text{ref}}\) (except the side specified as the source), and the center of the box (width away from the box sides) will have a concentration \(N_{\text{max}}\).

\[N = \{N_{\text{max}}\} \exp \left[ - \ln \left( \frac{\{N_{\text{max}}\}}{\{N_{\text{ref}}\}} \right) \left( \frac{(x|y|z) - \{(x|y|z)_{\text{box}}\}}{\text{width}} \right)^2 \right]\]

# doping with boxes
acceptor_boxes = []
donor_boxes = []

acceptor_boxes.append(
    td.ConstantDoping.from_bounds(
        rmin=[-5, 0, -np.inf],
        rmax=[5, 0.22, np.inf],
        concentration=1e15
    )
)

# p implant
acceptor_boxes.append(
    td.GaussianDoping.from_bounds(
    rmin=[-6, -0.3, -np.inf],
    rmax=[-0.15, 0.098, np.inf],
    concentration=7e17,
    ref_con=1e6,
    width=0.1,
    source="ymax")
)

# n implant
donor_boxes.append(
    td.GaussianDoping.from_bounds(
        rmin=[0.15, -0.3, -np.inf],
        rmax=[6, 0.098, np.inf],
        concentration=5e17,
        ref_con=1e6,
        width=0.1,
        source="ymax"
    )
)

# p++
acceptor_boxes.append(
    td.GaussianDoping.from_bounds(
        rmin=[-6, -0.3, -np.inf],
        rmax=[-2, 0.22, np.inf],
        concentration=1e19,
        ref_con=1e6,
        width=0.1,
        source="ymax"
    )
)

# # n++
donor_boxes.append(
    td.GaussianDoping.from_bounds(
        rmin=[2, -0.3, -np.inf],
        rmax=[6, 0.22, np.inf],
        concentration=1e19,
        ref_con=1e6,
        width=0.1,
        source="ymax")
)


# p wg implant
acceptor_boxes.append(
    td.GaussianDoping.from_bounds(
        rmin=[-0.3, 0, -np.inf],
        rmax=[0.06, 0.255, np.inf],
        concentration=5e17,
        ref_con=1e6,
        width=0.12,
        source="xmin"
    )
)

# n wg implant
donor_boxes.append(
    td.GaussianDoping.from_bounds(
        rmin=[-0.06, 0.02, -np.inf],
        rmax=[0.25, 0.26, np.inf],
        concentration=7e17,
        ref_con=1e6,
        width=0.11,
        source="xmax")
)

Once we have our doping distributions, we can add it to our td.SemiconductorMedium and create our MultiPhysicsMedium to create our device structure with it. To include the doping to our semiconductor properties, we can implement this by creating a copy of it through the function updated_copy().

Si_2D_doping = td.MultiPhysicsMedium(
    charge=intrinsic_si.updated_copy(
        N_d=donor_boxes,
        N_a=acceptor_boxes,
    ),
    name="Si_doping",
)

Compose our device structures

We can now combine our medium definitions and geometry models to create our device structures.

# create objects
overlap_factor = 1.0001

oxide = td.Structure(
    geometry=td.Box(center=(0, h_core, 0), size=(10, 5, td.inf)),
    medium=SiO2,
    name="oxide"
)

core_p = td.Structure(
    geometry=td.Box(center=(-w_core/2, h_core/2, 0), size=(w_core, h_core, td.inf)),
    medium=Si_2D_doping,
    name="core_p"
)

core_n = td.Structure(
    geometry=td.Box(center=(w_core/2, h_core/2, 0), size=(w_core, h_core, td.inf)),
    medium=Si_2D_doping,
    name="core_n"
)

clearance_p = td.Structure(
    geometry=td.Box(center=(-w_core - w_clearance/2, h_clearance/2, 0), size=(w_clearance, h_clearance, td.inf)),
    medium=Si_2D_doping,
    name="clearance_p"
)

clearance_n = td.Structure(
    geometry=td.Box(center=(w_core + w_clearance/2, h_clearance/2, 0), size=(w_clearance, h_clearance, td.inf)),
    medium=Si_2D_doping,
    name="clearance_n"
)

side_p = td.Structure(
    geometry=td.Box(center=(-w_core - w_clearance - w_side/2, h_side/2, 0), size=(w_side, h_side, td.inf)),
    medium=Si_2D_doping,
    name="side_p"
)

side_n = td.Structure(
    geometry=td.Box(center=(w_core + w_clearance + w_side/2, h_side/2, 0), size=(w_side, h_side, td.inf)),
    medium=Si_2D_doping,
    name="side_n"
)

# create a couple structs to define the contacts
contact_p = td.Structure(
    geometry=td.Box(center=(-w_core - w_clearance - w_side + w_contact/2, h_side + h_contact/2, 0), size=(w_contact, h_contact, td.inf)),
    medium=aux,
    name="contact_p"
)

contact_n = td.Structure(
    geometry=td.Box(center=(w_core + w_clearance + w_side - w_contact/2, h_side + h_contact/2, 0), size=(w_contact, h_contact, td.inf)),
    medium=aux,
    name="contact_n"
)

A td.Scene helps us visualise that the structures are correctly positioned in this simulation.

# create a scene with the previous structures
all_structures = [oxide, core_p, core_n, clearance_n, clearance_p, side_p, side_n, contact_p, contact_n]

scene = td.Scene(
    medium=air,
    structures=all_structures,
)


_, ax = plt.subplots(2, 1, figsize=(6, 6))

scene.plot(z=0, ax=ax[0])
scene.plot(y=h_core/2, ax=ax[1])
plt.tight_layout()
plt.show()

And to make sure our doping concentration regions are correct, we can visualize them with the convencience function `plot_structures_property <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.Scene.html#tidy3d.Scene.plot_structures_property>`__.

In the plots below, we have used property="doping", but we can also use property="acceptors" or property="donors" to check individual doping concentration regions.

# plot doping from scene
_, ax = plt.subplots(1,2, figsize=(15, 4))
scene.plot_structures_property(z=0, property="doping", ax=ax[0], limits=[-1e18,1e18])
scene.plot_structures_property(z=0, property="doping", ax=ax[1], hlim=[-0.5,0.5], vlim=[-0.5,0.5], limits=[-5e17,5e17])

<Axes: title={'center': 'cross section at z=0.00 (μm)'}, xlabel='x (μm)', ylabel='y (μm)'>

Charge Boundary Conditions

In this example, we will apply a DC voltage bias across the modulator junction at different voltages. We can do this by using the VoltageBC defined at both the p and n regions of the simulation region. This boundary condition accepts a `td.DCVoltageSource <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.DCVoltageSource.html#tidy3d.DCVoltageSource>`__ model. This DCVoltageSource accepts a single float or an array of floats. In the latter case, a computation result for each of the voltages defined in the array will be calculated.

The modulator will operate in reverse bias. This means a zero voltage is applied to the p side and a scalar positive voltage is applied to the n side.

# create BCs
voltages = [-0.5, 0, 1, 2, 3, 4]

bc_v1 = td.HeatChargeBoundarySpec(
    condition=td.VoltageBC(source=td.DCVoltageSource(voltage=0)),
    placement=td.StructureBoundary(structure=contact_p.name),
)

bc_v2 = td.HeatChargeBoundarySpec(
    condition=td.VoltageBC(source=td.DCVoltageSource(voltage=voltages)),
    placement=td.StructureBoundary(structure=contact_n.name),
)

boundary_conditions = [bc_v1, bc_v2]

Create our Monitors

Monitors are used to determine exactly what data is to be extracted from the computation. In the case of a Charge simulation, we can extract electronic properties such as the potential \(\psi\), free carriers distribution \(\mathbf{J}\), and capacitance \(C\).

A td.SteadyCapacitanceMonitor is used to monitor small signal capacitance within the area/volume defined by the monitor.

The small signal-capacitance of electrons \(C_n\) and holes \(C_p\) is computed from the charge due to electrons \(Q_n\) and holes \(Q_p\) at an applied voltage \(V\) at a voltage difference \(\Delta V\) between two simulations.

\[C_{n,p} = \frac{Q_{n,p}(V + \Delta V) - Q_{n,p}(V)}{\Delta V}\]

Note that if only one solution is calculated (one single voltage data point is computed) the output for this monitor will be empty.

In order to extract the electrostatic potential \(\psi\) from our simulation, a td.SteadyPotentialMonitor results in a td.SteadyPotentialData.

A td.SteadyFreeCarrierMonitor let us record the steady-state of the free carrier concentration (\(\mathbf{J_n}\) electrons and \(\mathbf{J_p}\) holes) so that we can visualize results.

# capacitance monitors
capacitance_global_mnt = td.SteadyCapacitanceMonitor(
    center=(0, 0.14, 0), size=(td.inf, td.inf, 0), name="capacitance_global_mnt",
)

# charge monitor around the waveguide
capacitance_wvg_mnt = td.SteadyCapacitanceMonitor(
    center=(0, 0.14, 0), size=(0.6, 0.3, 0), name="capcitance_wvg_mnt",
)

# charge monitor around the waveguide
charge_3D_mnt = td.SteadyFreeCarrierMonitor(
    center=(0, 0.14, 0), size=(0.6, 0.3, td.inf), name="charge_3D_mnt", unstructured=True
)

charge_monitor_z0 = td.SteadyFreeCarrierMonitor(
    center=(0, 0.14, 0), size=(0.6, 0.3, 0), name="charge_z0", unstructured=True,
)

# voltage monitor around waveguide
voltage_monitor_z0 = td.SteadyPotentialMonitor(
    center=(0, 0.14, 0), size=(0.6, 0.3, 0), name="voltage_z0", unstructured=True,
)

# Will be used later for the mode simulations
charge_monitor_z0_big = td.SteadyFreeCarrierMonitor(
    center=(0, 0.14, 0), size=(td.inf, td.inf, 0), name="charge_z0_big", unstructured=True, conformal=True
)

Mesh and Convergence criteria

Before we can actually create the simulation object, we need to define some mesh parameters and convergence/tolerance criteria. One important note here is that in Charge simulations, we need to specify the type of simulation we’re interested in running. In this case, we want to run a steady state isothermal analysis which is specified with the class IsothermalSteadyChargeDCAnalysis. Note that this class takes an argument convergence_dv. To understand how this parameter affects convergence, let’s imagine a case where we want to solve for 0 V and 20 V bias. In this case the solver will first try to solve the 0V bias case and based off of this solution will go ahead and try to solve the 20V bias case. Since these are two very different biases, the solver may find it difficult to solve the second. By setting e.g.,convergence_dv=10 we’ll force the solver to go from 0V to 20V in 10V steps which eases convergence.

# mesh
mesh = td.DistanceUnstructuredGrid(dl_interface=res*1.14, dl_bulk=res*4, distance_interface=0.02*h_side,
    distance_bulk=1*h_side, relative_min_dl=0,
    uniform_grid_mediums=[Si_2D_doping.name],
    non_refined_structures=[oxide.name]
)

# devsim setting
convergence_settings = td.ChargeToleranceSpec(rel_tol=1e5, abs_tol=3e3, max_iters=400)

# currently we support Isothermal cases only. Temperature in K
analysis_type = td.IsothermalSteadyChargeDCAnalysis(temperature=300, convergence_dv=10, tolerance_settings=convergence_settings)

Create our HeatChargeSimulation

We’re now ready to create the actual simulation object.

# build heat simulation object
charge_sim = td.HeatChargeSimulation(
    sources=[],
    monitors=[capacitance_global_mnt, capacitance_wvg_mnt, charge_3D_mnt, charge_monitor_z0, voltage_monitor_z0, charge_monitor_z0_big],
    analysis_spec=analysis_type,
    center=(0,0,0),
    size=(10.5,10,z_size),
    structures=all_structures,
    medium=air,
    boundary_spec=boundary_conditions,
    grid_spec=mesh,
    symmetry=(0, 0, 0)
    )

# plot simulation
fig, ax = plt.subplots(1, 2, figsize=(10, 15))
charge_sim.plot(z=0, ax=ax[0])
charge_sim.plot_property(z=0, property="electric_conductivity", ax=ax[1])
plt.tight_layout()
plt.show()

Run Simulation

#%matplotlib inline
from tidy3d import web
charge_data=web.run(charge_sim, task_name="charge_junction", path="charge_junction.hdf5")

11:32:34 CET Created task 'charge_junction' with task_id
             'hec-a49326a6-6361-4c64-b7fe-b4c0fa3d8894' 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.

11:32:37 CET Maximum FlexCredit cost: 0.025. Minimum cost depends on task
             execution details. Use 'web.real_cost(task_id)' to get the billed
             FlexCredit cost after a simulation run.

11:32:38 CET 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.

11:32:54 CET status = preprocess

11:33:39 CET starting up solver

             running solver

11:45:11 CET status = postprocess

11:47:27 CET status = success

11:47:37 CET loading simulation from charge_junction.hdf5

11:47:38 CET WARNING: Warning messages were found in the solver log. For more   
             information, check 'SimulationData.log' or use                     
             'web.download_log(task_id)'.                                       

Post-process Charge simmulation

Let’s begin by visualizing the potential and electron fields at three different biases to make sure the solution looks reasonable. As it can be seen in the figure, as the (reverse) bias increases, the depletion region widens, as it is expected. The electric potential field also adapts to the applied biases and the amount of charge present in the waveguide area, which is also to be expected.

voltages = charge_data[charge_3D_mnt.name].holes.values.voltage.data

fig, ax = plt.subplots(2, 3)
for n, index in enumerate([0,2,5]):
    # let's read voltage first
    charge_data[voltage_monitor_z0.name].potential.sel(voltage=voltages[index]).plot(ax=ax[0][n], grid=False)

    # now let's plot some electrons
    np.log10(charge_data[charge_3D_mnt.name].electrons.sel(z=0, voltage=voltages[index])).plot(ax=ax[1][n], grid=False)
    ax[1][n].set_title("Bias {0:0.2f}V".format(voltages[index]))

plt.tight_layout()

/home/momchil/miniconda3/envs/flex/lib/python3.12/site-packages/xarray/core/computation.py:818: RuntimeWarning: divide by zero encountered in log10
  result_data = func(*input_data)

2. Experimental Data Validation

In the second section of this example, we will validate the simulation numerical results with the experimental results from [1].

Small Signal Capacitance

The td.HeatChargeSimulationData contains the extracted td.SteadyCapacitanceData specified by the monitor. The small signal capacitance is computed both with electrons and holes, and is averaged. As seen below, there is reasonable agreement between experimental data and simulation data (within the convergence tolerance).

# capacitance from monitor - waveguide area
CV_baehrjones = [
    [-0.4, -0.25, 0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4],
    [0.261, 0.248, 0.223, 0.208, 0.198, 0.190, 0.184, 0.175, 0.168, 0.157, 0.150]
]

mnt_v = np.array(charge_data[capacitance_global_mnt.name].electron_capacitance.coords["v"].data)
mnt_ce = np.array(charge_data[capacitance_global_mnt.name].electron_capacitance.data)
mnt_ch = np.array(charge_data[capacitance_global_mnt.name].hole_capacitance.data)

plt.plot(mnt_v, -0.5*(mnt_ce+mnt_ch) * 10 , 'k.-', label="Simulation Data")
plt.plot(CV_baehrjones[0], np.array(CV_baehrjones[1])*10, 'r-', label="Experimental Data")

plt.xlabel("Reverse bias (V)")
plt.ylabel("pF/cm")
plt.legend()
plt.grid()

Charge-to-Optical Coupling

In this section, we will explore how we can couple the Charge simulation data into an optical FDTD simulation. This can be implemented by defining a medium perturbation. This means that optical medium properties such as the refractive index change e.g. \(\Delta n \propto \varepsilon(N_e, N_h)\) and optical loss \(\alpha \propto \sigma(N_e, N_h)\) are proportional to changes to the permittivity \(\varepsilon\) and conductivity \(\sigma\) by the free electrons \(N_e\) and free holes \(N_h\) density.

We will use empiric relationships presented in M. Nedeljkovic, R. Soref and G. Z. Mashanovich, “Free-Carrier Electrorefraction and Electroabsorption Modulation Predictions for Silicon Over the 1–14- μm Infrared Wavelength Range,” IEEE Photonics Journal, vol. 3, no. 6, pp. 1171-1180, Dec. 2011, that state that changes in \(n\) and \(k\) of Si can be described by formulas

\[- \Delta n = \frac{dn}{dN_e}(\lambda) (\Delta N_e)^{a(\lambda)} + \frac{dn}{dN_h}(\lambda) (\Delta N_h)^{\beta(\lambda)}\]

\[\Delta \left( \frac{4 \pi k}{\lambda} \right) = \frac{dk}{dN_e}(\lambda) (\Delta N_e)^{\gamma(\lambda)} + \frac{dk}{dN_h}(\lambda) (\Delta N_h)^{\delta(\lambda)}\]

where \(\Delta N_e\) and \(\Delta N_h\) are electron and hole densities.

We start by defining the optical frequency range of interest:

wvl_um = 1.55
freq0 = td.C_0 / wvl_um

fwidth = freq0 / 5
freqs = np.linspace(freq0 - fwidth / 10, freq0 + fwidth / 10, 201)
wvls = td.C_0 / freqs

We extract the real (\(n\)) and imaginary (\(k\)) parts of the refractive index from our silicon optical medium from our material library.

si = td.material_library['cSi']['Palik_Lossless']
n_si, k_si = si.nk_model(frequency=td.C_0 / wvl_um)
si_non_perturb = td.Medium.from_nk(n=n_si, k=k_si, freq=freq0)

The charge parameters have the following values for wavelength of 1.55 \(\mu\)m:

\(\lambda\)	\(\frac{dn}{dN_e}\)	\(a\)	\(\frac{dn}{dN_h}\)	\(\beta\)	\(\frac{dk}{dN_e}\)	\(\gamma\)	\(\frac{dk}{dN_h}\)	\(\delta\)
\(1.55\)	\(5.40 \times 10^{-22}\)	\(1.011\)	\(1.53 \times 10^{-18}\)	\(0.838\)	\(8.88 \times 10^{-21}\)	\(1.167\)	\(5.84 \times 10^{-20}\)	\(1.109\)

ne_coeff = -5.4e-22 * 1.1
ne_pow = 1.011

nh_coeff = -1.53e-18 * 1.1
nh_pow = 0.838

k_factor = wvl_um * 1e-4 / 4 / np.pi  # factor for conversion from absorption coefficient into k

ke_coeff = k_factor * 8.88e-21 * 1.1
ke_pow = 1.167

kh_coeff = k_factor * 5.84e-20 * 1.1
kh_pow = 1.109

Given the nonlinear character of these dependencies we will incorporate them as sampled function on a rectangular grid formed by electron and hole density values. Specifically, we will sample the given \(n\) and \(k\) dependencies in the electron and hole density ranges up to \(10^{20}\) 1/cm\(^3\).

Ne_range = np.concatenate(([0], np.logspace(-6, 20, num=200)))
Nh_range = np.concatenate(([0], np.logspace(-6, 20, num=200)))

Ne_mesh, Nh_mesh = np.meshgrid(Ne_range, Nh_range, indexing='ij')

dn_mesh = ne_coeff * Ne_mesh ** ne_pow + nh_coeff * Nh_mesh ** nh_pow
dk_mesh = ke_coeff * Ne_mesh ** ke_pow + kh_coeff * Nh_mesh ** kh_pow

Perturbation Mediums Specifications

We need to create an optical simulation. This means we need to slightly modify our existing Scene used in our HeatChargeSimulation. In order to do this, we will define the charge to optical coupling through a PerturbationMedium. This medium will enable us to input the numerical results from our HeatChargeSimulationData and change the optical medium properties.

Now we convert sampled values of \(n\) and \(k\) into permittivity \(\varepsilon\) and conductivity \(\sigma\) values, and assemble a non-dispersive medium with perturbations (PerturbationMedium). Note that since the optical simulations don’t yet support MultiPhysicsMedium, we need to change the scene medium.

dn_data = td.ChargeDataArray(dn_mesh, coords=dict(n=Ne_range, p=Nh_range))
dk_data = td.ChargeDataArray(dk_mesh, coords=dict(n=Ne_range, p=Nh_range))

n_si_charge = td.CustomChargePerturbation(perturbation_values=dn_data)
k_si_charge = td.CustomChargePerturbation(perturbation_values=dk_data)

n_si_perturbation = td.ParameterPerturbation(
    charge=n_si_charge,
)

k_si_perturbation = td.ParameterPerturbation(
    charge=k_si_charge,
)

si_perturb = td.PerturbationMedium.from_unperturbed(
    medium=si_non_perturb,
    perturbation_spec=td.IndexPerturbation(
        delta_n=n_si_perturbation,
        delta_k=k_si_perturbation,
        freq=freq0,
    )
)

new_structs = []
for struct in charge_sim.structures:
    # NOTE: applying perturbation material to Si
    if struct.medium.name == Si_2D_doping.name:
        new_structs.append(struct.updated_copy(medium=si_perturb))

scene = td.Scene(
    medium=SiO2.optical, # currently td.Simulation cannot accpet a MultiphysicsMedium
    structures=new_structs,
)

Creating an Optical Mode Simulation

Let’s create an optical mode simulation with our updated scene.

span = 2*wvl_um

port_center = (0,h_core,-span/2)
port_size = (5, 5, 0)

buffer = 1*wvl_um
sim_size = (13+buffer, 10+buffer, span)

bc_spec = td.BoundarySpec(
    x = td.Boundary.pml(num_layers=20),
    y = td.Boundary.pml(num_layers=30),
    z = td.Boundary.periodic(),
)

sim = td.Simulation(
    center=(0,0,0),
    size=sim_size,
    medium=scene.medium,
    structures=scene.structures,
    run_time=6e-12,
    boundary_spec=bc_spec,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=60, wavelength=wvl_um),
)

_, ax = plt.subplots(1, 2, figsize=(15, 7))
sim.plot(y=sim.center[1], ax=ax[0])
sim.plot(z=sim.center[2], ax=ax[1])

plt.tight_layout()
plt.show()

Coupling our SteadyFreeCarrierData into our PerturbationMedium:

We already have configured our reference optical mode simulation with the default free carrier distribution. However, our td.HeatChargeSimulationData with from our td.SteadyFreeCarrierData can be used to perturb our silicon optical medium. Let’s see how we can do this below.

We can create a list of perturbed optical simulation by using the function perturbed_mediums_copy():

def apply_charge(charge_data):
    perturbed_sims = []
    for n, v in enumerate(charge_data[charge_monitor_z0_big.name].electrons.values.voltage.data):
        e_data = charge_data[charge_monitor_z0_big.name].electrons.sel(voltage=v)
        h_data = charge_data[charge_monitor_z0_big.name].holes.sel(voltage=v)
        perturbed_sims.append(
            sim.perturbed_mediums_copy(
                electron_density=e_data,
                hole_density=h_data,
            )
        )
    return perturbed_sims

perturbed_sims = apply_charge(charge_data)

Let’s check what the difference in permittivity looks like at different applied voltages:

_, ax = plt.subplots(1, 2, figsize=(15, 5))

sampling_region = td.Box(center=(0, h_core/2, 0), size=(6, 5, 1))
eps_undoped = sim.epsilon(box=sampling_region).sel(z=0, method="nearest")

for ax_ind, ind in enumerate([0, len(voltages)-1]):
    eps_doped = perturbed_sims[ind].epsilon(box=sampling_region).sel(z=0, method="nearest")
    eps_doped = eps_doped.interp(x=eps_undoped.x, y=eps_undoped.y)
    eps_diff = np.abs(np.real(eps_doped - eps_undoped))
    eps_diff.plot(x="x", ax=ax[ax_ind])

    ax[ax_ind].set_aspect("equal")
    ax[ax_ind].set_title(f"Bias: {voltages[ind]:1.1f} V")
    ax[ax_ind].set_xlabel("x (um)")
    ax[ax_ind].set_ylabel("y (um)")

plt.tight_layout()
plt.show()

Running Waveguide Mode Simulations

Instead of running a full FDTD simulation, we will compute the modes of the waveguide. This is a smaller computation that will also highlight the influence of free carriers fields over the refraction coefficient.

We will compute the waveguide modes on a plane centered around the core section of the device:

from tidy3d.plugins.mode import ModeSolver

mode_plane = td.Box(size=port_size)

# visualize
ax = sim.plot(z = 0)
mode_plane.plot(z=0, ax=ax, alpha=0.5)
plt.show()

mode_solvers = dict()
i = 0
for psim in perturbed_sims:
    ms = ModeSolver(
        simulation=psim,
        plane=mode_plane,
        freqs=np.linspace(freqs[0], freqs[-1], 11),
        mode_spec=td.ModeSpec(num_modes=1, precision="double"),
    )
    mode_solvers[i] = ms
    i +=1

# server mode computation
ms_batch_data = web.Batch(simulations=mode_solvers, reduce_simulation=True).run()

11:57:01 CET Started working on Batch containing 6 tasks.

11:57:15 CET Maximum FlexCredit cost: 0.123 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after the
             Batch has completed.

12:00:23 CET Batch complete.

Relative Phase Change \(\Delta \phi\)

Based on the computed effective refractive index \(n_{eff}\), we can now determine both phase shift \(\Delta \phi\) and \(\alpha\) optical loss per cm as a function of the applied bias voltage \(V_n\) using the following equations from [1]:

\[\Delta n_{\text{eff}} = \text{Re}(n_{\text{eff}}(V_n) - n_{\text{eff}}({V_0}))\]

\[\Delta \phi = \frac{2\pi \Delta n_{\text{eff}}}{\lambda_0}\]

\[\alpha_{\text{dB/cm}} = 10 \cdot 4 \pi \cdot \text{Im}(n_{\text{eff}}) / \lambda_0 \cdot 10^4 \cdot \log_{10}(e)\]

We compare this to simulation data from another popular simulation software.

n_eff_freq0 = [md.n_complex.sel(f=freq0, mode_index=0).values for md in ms_batch_data.values()]

ind_V0 = 1
delta_neff = np.real(n_eff_freq0 - n_eff_freq0[ind_V0])
rel_phase_change = 2*np.pi*delta_neff/wvl_um * 1e4
alpha_dB_cm = 10*4*np.pi*np.imag(n_eff_freq0)/wvl_um*1e4*np.log10(np.exp(1))

# other results
v_other = [-0.5, 0, 0.4, 0.7, 1.5, 2.4, 3.8]
pc_other=[-0.91, 0, 0.6, 1,   1.8, 2.7, 3.6]
loss_other = [4.9, 4.55, 4.34, 4.15, 3.77, 3.41, 2.95]

_, ax = plt.subplots(1, 2, figsize=(15, 4))
ax[0].plot(voltages, rel_phase_change, 'k.-', label="Simulation Data")
ax[0].plot(v_other, pc_other, 'r-', label="Other Software")
ax[0].set_xlabel("Reverse bias (V)")
ax[0].set_ylabel("Relative phase (rad/cm)")
ax[0].grid()
ax[0].legend()

ax[1].plot(voltages, alpha_dB_cm, 'k.-', label="Simulation Data")
ax[1].plot(v_other, loss_other, 'r-', label="Other Software")
ax[1].set_xlabel("Reverse bias (V)")
ax[1].set_ylabel("Loss (dB/cm)")
ax[1].grid()
ax[1].legend()
plt.tight_layout()
plt.show()

As it can be seen, we’d need to apply a bias of ~3V in order to get a phase shift of \(\pi\).