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. Subsequent heat distribution and waveguide mode analysis is done with Tidy3D's heat and mode solvers.

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.

Charge simulation¶

We will begin by creating a Charge simulation. To define a Charge simulation 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. One of the methods we have to create doping profiles in Tidy3d is by defining doping boxes. We'll be using ConstantDoping. One thing to note is that these doping boxes are additive, i.e., if two donor doping boxes overlap the total concentration in the overlap region will be sum of these two overlapping doping boxes. To following further demonstrate this concept

We can use this doping in our semiconductor material

Charge structures¶

Now that we have our semiconductor material let's create some structures that will use it.

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

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

Charge boundary conditions¶

Since we're interested in the response of the system for different applied voltages we'll need to solve the charge problem at each of these voltages. In Charge this can readily be done since the VoltageBC can accept an array of voltages as source through DCVoltageSource. 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.

Charge 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 monitor has been defined in the plane.

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 an relative tolerance of $1\cdot 10^{-10}$ 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 an isothermal DC case which we can define the simulation with IsothermalSteadyChargeDCAnalysis. 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.

We'll use a spatial resolution of $0.005\mu m$.

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.

We can also plot here the simulation and some properties. In the properties plot, as well as the conductivity (which is invisible since the only conducting structures are the BC auxiliary ones) we can see the boundary conditions in blue and orange.

Solution visualization¶

We can explore the resulting hole density distribution as well as the resulting potential field at two different applied biases.

Whenever we run a charge simulation, Charge will output some device characteristics. These are not tight 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.

Heat and optics scene creation¶

Now that we have currents and carrier distributions for the waveguide we can start building the rest of the domain for heat and optics calculations. The waveguide is embedded in a silica cladding with the following dimensions:

Next we create the materials to be used. This is done through the functions Medium() and PerturbationMedium() of Tidy3D. The former is used for media with constant properties whereas we use the latter for materials with varying properties. In particular, 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. In this example we have encapsulated the material creation in the following function:

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

We will reuse the structures generated for the Charge simulation here so only the cladding and substrate are left to be created. Note that, since we need to use perturbation materials for the silicon structures, we will update the material for those structures.

Once we have the different structures we can create a scene and visualize it to make sure it's correct

Heat simulation¶

Let's now start creating our heat simulation

Boundary conditions¶

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.

Heat source¶

Next we create a heat source. The total input power to the system is the electric power, computed here as $P=V\cdot i$. Remember that the electric current computed previously comes from a 2D case and, therefore, it is current per unit length.

Additionally, the heat solver takes in volumetric sources applied over a defined volume. Since we will apply the heat source over the waveguide, we will divide the total power by the volume of the waveguide.

Temperature monitor¶

Next, we define a temperature monitor. We will be able to visualize the temperature field at the monitor.

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.

Heat simulation object¶

At this point we can finally create the heat simulation object as well as visualize it.

The red line indicates the convection BC whereas the yellow line depicts the temperature BC. The dotted area indicates where the heat source is being applied.

Note that we slightly expanded simulation in the y-direction to avoid structures' bounds to extend exactly to simulation edges and possibly causing an unexpected behavior.

The following runs the simulation.

Temperature field visualization¶

Once we have the previous heat simulation we can check the mesh and the resulting temperature field.

In the first plot we can check the thermal conductivity of the different materials. As it can be seen, silicon structures are much more conductive than the SiO2 cladding.

The last two pictures show the mesh and temperature field.

Once we have made sure that the mesh and heat simulation parameters are acceptable, we can run the simulations for each voltage in the voltages array. To do this we will first create a simulation object for each of the voltages:

Once the heat simulations have finished we can visualize a couple of them to make sure they're good:

Combined Electric and Thermal Perturbations¶

Once we have run the heat simulations, we can apply the temperature fields obtained along with the free carrier distributions obtained previously to our optical simulations.

As a first step we need to create a reference optics simulation:

Now that we have the reference optical simulation we can create a list of perturbed optical simulations. When passing charge and thermal fields to function perturbed_mediums_copy() one must make sure that the provided data fields are sampled on the same grid (structured or unstructured). Since in this example charge data and thermal data are coming from different grids, we first interpolate them onto a common grid, for which we select the Cartesian grid of the optic simulation.

Now that we have created the perturbed simulation we can examine the change 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.

Based on the computed refractive index we can now compute both phase shift and loss over a waveguide length of 100$\mu 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 3.5mW$ in order to get a phase shift of $\pi$.