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¶

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$.

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:

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.

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

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:

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 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] $$

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().

Compose our device structures¶

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

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

And to make sure our doping concentration regions are correct, we can visualize them with the convenience function 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.

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 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 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.

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.

Create our HeatChargeSimulation¶

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

Run Simulation¶

Post-process Charge simulation¶

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.

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).

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 start by defining the optical frequency range of interest:

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

Perturbation Mediums Specifications¶

To compute the pertubationWe 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.

For convenience this model has been implemented as NedeljkovicSorefMashanovich

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.

Creating an Optical Mode Simulation¶

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

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 our td.SteadyFreeCarrierData can be used to perturb our silicon optical medium. Let's see how we can do this below. Note that we will be using a subset of the applied voltages to reduce the computational resources

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

Let's check what the difference in permittivity looks like at different applied voltages, compared to zero voltage.

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:

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.

As it can be seen, we'd need to apply a bias of ~3.3V in order to get a phase shift of $\pi$.

It is worth noting here that the refractive index and optical loss perturbations from the model in 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 change significantly at the specified wavelength of 1.55 $\mu m$ and results will also change accordingly. For instance, if one computed the perturbations based on a reference wavelength of 1.575$\mu m$ the following results are obtained: