In this example we investigate the modulation of a silicon waveguide with an embedded doped silicon p++ - p - p++ junction inspired by the design presented in Chuyu Zhong, Hui Ma, Chunlei Sun, Maoliang Wei, Yuting Ye, Bo Tang, Peng Zhang, Ruonan Liu, Junying Li, Lan Li, and Hongtao Lin, "Fast thermo-optical modulators with doped-silicon heaters operating at 2 μm", Opt. Express 29, 23508-23516 (2021). The primary modulation effect is due to heating occurring inside the waveguide as a potential difference is applied across the p++ - p - p++ junction. For calculation of current flowing through the device we use our Charge solver. The non-isothermal analysis option allows the solver to also provide the temperature distribution from both Joule heating and generation/recombination heat sink/source.

The sketch below highlights the geometric features of the considered waveguide:

We begin by importing the required modules for this example.

Next, let's define some geometric properties of the waveguide as well as the doping. We define three doping regions within the Si waveguide with concentrations p++ - p - p++. The interface between the high doping concentration and low doping concentration is delimited by [-x_i, x_i] with $x=0$ being the center of the waveguide.

Let's now define some reference temperature and wavelength:

We define next dimensions for the cladding and buried oxide box:

Next we define some optical and thermal properties for Si and SiO2. Thermal properties can be extracted from www.efunda.com, www.periodic-table.org and www.azom.com

Heat-Charge simulation¶

We will begin by creating a HeatChargeSimulation object. To define it we will need to create at least one SemiconductorMedium. For this type of medium we can define the doping distributions that will be applied to all structures that have this as their medium.

Let's start, then, with creating the doping. There are two main ways of defining doping in Tidy3d: supplying a SpatialDataArray with the desired doping distribution or through doping boxes.

In what follows we'll use doping boxes. One thing to note is that these boxes are additive meaning that regions with boxes overlap will see the sum of the doping corresponding to each of the boxes.

Additionally, there are different types of doping boxes. Here we'll be using ConstantDoping and GaussianDoping. Gaussian doping allow for smooth transition of one doping level to another.

Next we'll use this doping to define the semiconductor properties.

You can see how materials are defined with the MultiPhysicsMedium class which allows to define physical group of properties, e.g., optical, thermal, and electrical properties.

The exception is air which is defined as Medium. This is because air structures will be removed from our Heat-Charge simulation.

Structures¶

We create now all the structures used in both the optic and Heat-Charge simulations.

For ease in defining boundary conditions we define here some auxiliary structures made up of a conductive medium defined with ChargeConductorMedium

Scene¶

We can now put those structures together into a scene and explore what it looks like

Let's also inspect the doping profiles and make sure it is what we were expecting. Note the smooth transition of doping between the p and p+ regions. These smooth transitions avoid sharp variations in electric potential hence avoiding artificially high electric field. This is relevant in non-isothermal simulations since some of the heat source terms are proportional to the electric field.

Boundary conditions¶

Since we're interested in the response of the system for different applied voltages we'll need to solve the Heat-Charge problem at each of these voltages.

First we deal with electric boundary conditions. For this case we apply VoltageBC to both conductors. This BC requires the definition of a source. We'll be using DCVoltageSource which can accept and array of voltages. When an array is provided a parameter scan will be run and the returned data will have the provided voltage values as a separate dimension.

Let's define forward bias values up to 5 V with a step of 0.5 V.

We next create thermal boundary conditions. We define two BCs: a temperature BC applied to the auxiliary Si structure; and convection boundary condition applied to the interface between the mediums air and SiO2.

Monitors¶

Since we're interested in obtaining the free carrier distribution we'll add a SteadyFreeCarrierMonitor to our Charge simulation. For visualization purposes we will also include a potential monitor, i.e., SteadyPotentialMonitor

Note that the below monitors have been defined with infinite size to cover the full domain. This will be later helpful when we interpolate carrier and temperature data onto the optics simulation domain.

Since we're also interested in the temperature field we need to define a temperature monitor:

Heat-Charge simulation object¶

When running Charge, we need to define the type of Charge simulation and some convergence settings. In the current case we set a relative tolerance of $1\cdot 10^{-5}$ and an absolute tolerance of $5\cdot 10^{10}$. The absolute tolerance may seem big though one should notice we have variables (electrons/holes) that take on values many orders of magnitude larger than the tolerance ($\approx 10^{20}$).

In the current case, we are going to run a non isothermal DC case which we can define the simulation with SteadyChargeDCAnalysis. In DC mode we can set the parameter convergence_dv which tells the solver to limit the size of the sweep, i.e., if we need to solve for a bias of 0 and 0.5 and we set convergence_dv=0.1, it will force the solver to go between 0 and 0.5 at intervals of 0.1.

The last element we need to define before we can run the heat simulation is the mesh which we build using Tidy3D function DistanceUnstructuredGrid(). We will give this function 5 arguments:

• dl_interfacedefines the grid size near the interface.
• dl_bulkdefines the grid size away from the interface.
• distance_interface defines the distance from the interface until which dl_interface will be enforced.
• distance_bulkdefines the distance from the interface after which dl_bulk will be enforced.
• non_refined_structures allows us to specify structures which will not be refined.
• mesh_refinements allows us to specify an array of refinement regions.

We'll use a spatial resolution of $0.005\mu m$ to refine the semiconductor region.

We now have all the required elements to define a Charge simulation object. Note that the simulation has 0 size in the $z$ direction. With this, we'll make sure that the simulation is 2D even if the structures themselves are not.

Next let's visualize our simulation domain. In the upper plot we visualize the domain with the thermal conductivity superimposed to the structures. It can be seen that the semiconductor region has higher (darker) thermal conductivity. The boundary conditions can also be seen: the red line indicates the convection BC whereas the yellow line at the bottom indicates the temperature BC.

The plot at the bottom shows the simulation structures of the Charge simulation along with the voltage BCs highlighted with a yellow line.

Now we can run the Charge simulation.

Solution visualization¶

We can begin by examining the electron, hole, and potential distributions at a couple of applied biases:

We can also explore the temperature field at two different applied biases:

We can further inspect how electrons and holes as well as the potential vary across the semiconductor at different biases by plotting the distribution along the coordinate $y=h\_side/2$:

As it can be seen the majority carrier (holes) is dominated by the doping profile wheres the minority carrier (electrons) clearly changes with the different applied voltages. The potential varies roughly linearly along $x$ as a result of the hole concentration following the doping concentration.

Similarly, we can check the temperature distribution along $y=h\_side/2$

Whenever we run a charge simulation, Charge will output some device characteristics. These are not linked to any monitor but instead they are related to the whole device. For the current DC isothermal simulation we have the entry steady_dc_current_voltage.

Perturbation mediums¶

We next create perturbation materials. These include the necessary $k$ and $n$ perturbations to account for the effect of both temperature and carrier concentration.

We can create a perturbed medium using PerturbationMedium() of Tidy3D. In this case, the permittivity of Si will change with charge and temperature; and the permittivity of SiO2 will only vary with temperature. It should be noted that a linear variation with temperature is considered, i.e., the temperature dependence of permittivity with temperature has been linearized at the given reference temperature.

For the Si charge-induced variation, we will use the empirical 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 which has been implemented in the class NedeljkovicSorefMashanovich.

Note also that since the permittivity of the material is a function of wave frequency, we require a reference wavelength. A reference temperature must also be given for the same reason.

The warning above is telling us that the perturbation model does not check whether the perturbed permittivity or conductivity are physical so it is up to us to check that (which we'll do below). Since this message is going to be repeated every time we create a perturbation we will temporarily set the logging level to "ERROR"

Combined Electric and Thermal Perturbations¶

In order to reuse the scene we created for the non isothermal Charge simulation, we need to update the mediums of the structures. This can be done by looping through the structures of the scene and updating the medium accordingly as below.

Note that since the perturbation models can't ensure that the resulting permittivity stays above 1 a warning is issued every time a structure is assigned a perturbed medium.

With the structures now having perturbation mediums we can next update the scene and visualize.

Before we can apply thermal and carrier perturbations we need to generate a reference optic simulation object:

Now that we have the reference optical simulation we can create a list of perturbed optical simulations which we can do with perturbed_mediums_copy().

We can now verify that all perturbations have the real part of the permittivity greater than 1:

Now that we have created the perturbed simulation we can examine the change in permittivity due to these perturbations:

Perturbation Mode Analysis¶

Instead of running a full wave simulation we will compute the modes of the waveguide. This is a smaller computation that will also highlight the influence of the temperature field over the refraction coefficient.

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

The next few lines create the mode simulation objects for each of the temperature and carrier fields and submits the computation.

Run all mode solver simulations in a batch.

And finally we can see the combined effect of both thermal and carrier density on the refraction index.

In the following we're computing the electric power as $P_{electric}=V \cdot i$. It is reminded here that since the Charge simulation was 2D the current has dimensions of $A/\mu m$.

Also note that we're setting the logging level back to "WARNING".

Based on the computed refractive index we can now compute both phase shift and loss over a waveguide length of 100 µm as a function of the input (heating) power. As it can be seen, we'd need to apply a heating power of about $\approx 40mW$ in order to get a phase shift of $\pi$.