Drift-Diffusion Simulation of High-Speed Optoelectronic Devices

Abstract

In this paper, we address the problem of research and development of the advanced optoelectronic devices designed for on-chip optical interconnections in integrated circuits. The development of the models, techniques, and applied software for the numerical simulation of carrier transport and accumulation in high-speed A^IIIB^V (A and B refer to group III and V semiconductors, respectively) optoelectronic devices is the purpose of the paper. We propose the model based on the standard drift-diffusion equations, rate equation for photons in an injection laser, and complex analytical models of carrier mobility, generation, and recombination. To solve the basic equations of the model, we developed the explicit and implicit techniques of drift-diffusion numerical simulation and applied software. These aids are suitable for the stationary and time-domain simulation of injection lasers and photodetectors with various electrophysical, constructive, and technological parameters at different control actions. We applied the model for the simulation of the lasers with functionally integrated amplitude and frequency modulators and uni-travelling-carrier photodetectors. According to the results of non-stationary simulation, it is reasonable to optimize the parameters of the lasers-modulators and develop new construction methods aimed at the improvement of photodetectors’ response time.

1. Introduction

Nowadays, on-chip metal interconnections are very close to the physical limit of their application in state-of-the-art integrated circuits (ICs). It is caused by the appreciable degradation of response time, energy and technological efficiency, channel capacity, noise immunity, reliability, and other characteristics of integrated metal conductors at dimensions within several tens of nanometers and less, which are typical for modern integrated electronics [1]. The sustainable progress in IC production technology requires the immediate solution of the interconnection problem. It can be performed by means of the development of new interconnection types or optimization and modification of traditional interconnecting techniques. Both approaches have advantages and disadvantages and evolve concurrently.

This paper is focused on the research and development of novel high-performance on-chip interconnections for next-generation ICs. Super- and nanoconductors, radio frequency and optical commutation lines, carbon nanotubes, and graphene are the promising types of advanced integrated interconnections [2,3,4,5]. However, optical interconnecting is characterized by some important advantages over its counterparts and can solve the problem for various integrated devices in the immediate future [6,7,8].

The optoelectronic approach deals with the constructive and technological integration of optical interconnections with silicon-based electronic IC elements. Typical on-chip optical interconnection consists of a source of optical radiation, a high-speed modulator, an integrated waveguide, and a photodetector. Optoelectronics has made significant progress in the improvement of these devices, but a set of theoretical and practical problems remains unsolved. As of now, on-chip optoelectronic systems are suitable for highly specialized devices with sufficiently large scales or high-levels of interconnection in ICs. For example, optical interconnections can be used for the inter-core commutation in multi-core ultra-large-scale ICs [9].

Silicon photonics seem to be one of the advanced research directions in the field of IC optical interconnection [10,11,12]. This approach provides some important benefits, such as small losses, adaptability to streamlined manufacture, and high efficiency of waveguides. The main problem of silicon photonics is the development of on-chip optical radiation sources. This problem remains unsolved and causes the relevance of the design of hybrid A^IIIB^V-silicon systems [13,14,15] (A and B refer to group III and V semiconductors, respectively).

In this paper we consider the modelling of dynamic characteristics of A^IIIB^V high-speed optoelectronic devices for on-chip optical interconnections in ICs. The development of the appropriate model, simulation technique, and software is the purpose of the research.

The semiclassical and quantum-mechanical methods are the essential approaches to the simulation of semiconductor devices in modern computational electronics [16]. We propose to apply the drift-diffusion approximation of the semiclassical approach for the simulation of lasers and photodetectors. In spite of several limitations, it is utilized for the research of semiconductor devices and allows us to obtain adequate simulation results [17].

2. Models

In this paper we propose a numerical model designed for the simulation of charge carrier transport in high-speed A^IIIB^V optoelectronic devices within the drift-diffusion approximation of the semiclassical approach [16]. The model includes the following basic differential equations:

• The drift-diffusion equations for the electron and hole current densities $\overrightarrow{j_{\mathrm{n}}}$ and $\overrightarrow{j_{\mathrm{p}}}$ [18]:

  $$\overrightarrow{j_{\mathrm{n}}}=-q{\mathsf{\mu }}_{\mathrm{n}}\left[n\cdot \nabla \left(\mathsf{\phi }+V_{\mathrm{n}}\right)-{\mathsf{\phi }}_{\mathrm{T}}\cdot \nabla n\right];$$

  (1)

  $$\overrightarrow{j_{\mathrm{p}}}=-q{\mathsf{\mu }}_{\mathrm{p}}\left[p\cdot \nabla \left(\mathsf{\phi }-V_{\mathrm{p}}\right)+{\mathsf{\phi }}_{\mathrm{T}}\cdot \nabla p\right],$$

  (2)

  where $q$ is the electron charge; ${\mathsf{\mu }}_{\mathrm{n}}$, ${\mathsf{\mu }}_{\mathrm{p}}$ are the electron and hole mobilities; $n$, $p$ are the electron and hole densities; $\mathsf{\phi }$ is the electrostatic potential; $V_{\mathrm{n}}$, $V_{\mathrm{p}}$ are the heterostructure potentials in the conduction and valence bands of a semiconductor; ${\mathsf{\phi }}_{\mathrm{T}}$ is the temperature potential;
• The continuity equations for electrons and holes [18]:

  $$\frac{\partial n}{\partial t}=\frac{1}{q}\cdot \nabla \overrightarrow{j_{\mathrm{n}}}+G-R;$$

  (3)

  $$\frac{\partial p}{\partial t}=-\frac{1}{q}\cdot \nabla \overrightarrow{j_{\mathrm{p}}}+G-R,$$

  (4)

  where $t$ is time; $G$, $R$ are the generation and recombination rates of electron-hole pairs;
• The Poisson’s equation for electrostatic field [18]:

  $$\nabla \left(\mathsf{\epsilon }\cdot \nabla \mathsf{\phi }\right)=\frac{q}{{\mathsf{\epsilon }}_0}\cdot \left(n-p-N_{\mathrm{d}}+N_{\mathrm{a}}\right),$$

  (5)

  where $\mathsf{\epsilon }$ is the dielectric permittivity of a semiconductor; ${\mathsf{\epsilon }}_0$ is the permittivity of vacuum; $N_{\mathrm{d}}$, $N_{\mathrm{a}}$ are the densities of ionized donors and acceptors;
• The rate equation for photons in an injection laser [19]:

  $$\frac{\partial n_{\mathrm{ph}}}{\partial t}=-\frac{n_{\mathrm{ph}}}{{\mathsf{\tau }}_{\mathrm{ph}}}+\frac{\mathsf{\beta }\cdot {\left(np-n_i^2\right)}^{1/2}}{{\mathsf{\tau }}_s}+v_g\cdot g\cdot n_{\mathrm{ph}},$$

  (6)

  where $n_{\mathrm{ph}}$ is the photon density; ${\mathsf{\tau }}_{\mathrm{ph}}$ is the photon lifetime in the resonant cavity of a laser; $\mathsf{\beta }$ is the fraction of spontaneous emission in the laser mode; $n_{\mathrm{i}}$ is the intrinsic carrier concentration; ${\mathsf{\tau }}_s$ is the time constant of spontaneous recombination; $v_g$ is the photon velocity in the active region of a laser; $g$ is the optical gain coefficient.

If one considers a photodetector, it is enough to solve Equations (1)–(5), which represent the standard drift-diffusion formulation. The simulation of an injection laser requires the solution of additional Equation (6), describing the dynamics of photons leaving from the resonant cavity and photon generation caused by the spontaneous and stimulated radiative recombination of charge carriers.

To calculate heterostructure potentials $V_{\mathrm{n}}$ and $V_{\mathrm{p}}$, the following formulas are used:

$$V_{\mathrm{n}}=-{\mathsf{\phi }}_{\mathrm{C}}+{\mathsf{\phi }}_{\mathrm{C},\mathrm{ref}};V_{\mathrm{p}}={\mathsf{\phi }}_{\mathrm{V}}-{\mathsf{\phi }}_{\mathrm{V},\mathrm{ref}},$$

(7)

where ${\mathsf{\phi }}_{\mathrm{C}}$, ${\mathsf{\phi }}_{\mathrm{V}}$ are the potentials of conduction band bottom and valence band top; ${\mathsf{\phi }}_{\mathrm{C},\mathrm{ref}}$, ${\mathsf{\phi }}_{\mathrm{V},\mathrm{ref}}$ are the reference potentials of the corresponding bands. It is reasonable to choose the minimum of ${\mathsf{\phi }}_{\mathrm{C}}$ and maximum of ${\mathsf{\phi }}_{\mathrm{V}}$ in the energy band diagram of a heterostructure as the reference potentials for Equation (7). The data required for the band structure calculations are represented in [20].

The optical gain coefficient $g$ in Equation (6) depends on charge carrier and photon concentrations $n$ and $n_{\mathrm{ph}}$ in the active region of the laser. Traditionally, it is defined by the following analytical model [19]:

$$g=\{\begin{array}{l}{g}_0{\left(1+{\mathsf{\epsilon }}_{\mathrm{E}}\cdot \mathsf{\Gamma }\cdot n_{\mathrm{ph}}\right)}^{-1/2}\ln \left[\frac{A{\left(np-n_{\mathrm{i}}^2\right)}^{1/2}+B\left(np-n_{\mathrm{i}}^2\right)+C{\left(np-n_{\mathrm{i}}^2\right)}^{3/2}}{A{n}_0+B{n}_0^2+C{n}_0^3}\right],{\mathsf{\phi }}_{\mathrm{n}}-{\mathsf{\phi }}_{\mathrm{p}}\ge {\mathsf{\phi }}_{\mathrm{C}}-{\mathsf{\phi }}_{\mathrm{V}};\hfill \\ 0,{\mathsf{\phi }}_{\mathrm{n}}-{\mathsf{\phi }}_{\mathrm{p}}<{\mathsf{\phi }}_{\mathrm{C}}-{\mathsf{\phi }}_{\mathrm{V}},\hfill \end{array}$$

(8)

where $g_0$ is the proportionality factor; ${\mathsf{\epsilon }}_{\mathrm{E}}$ is the resulting factor of optical gain reduction; $\mathsf{\Gamma }$ is the factor of optical confinement; $A$, $B$, $C$ are the coefficients of the trap, radiative, and Auger recombination; $n_0$ is the threshold carrier concentration; ${\mathsf{\phi }}_{\mathrm{n}}$, ${\mathsf{\phi }}_{\mathrm{p}}$ are the quasi-Fermi levels for electrons and holes. In Equation (8), as well as in Equation (6), the generalized carrier concentration $N=n=p$ is replaced by the expression ${\left(np-n_{\mathrm{i}}^2\right)}^{1/2}$, which takes into account the electron, hole, and photon densities.

The following Dirichlet boundary conditions for Equations (1)–(6) are applied at ohmic and Schottky contacts in the case of fixed voltage mode [21]:

$$n=\left\{\frac{N_{\mathrm{d}}-N_{\mathrm{a}}}{2}+{\left[{\left(\frac{N_{\mathrm{d}}-N_{\mathrm{a}}}{2}\right)}^2+n_{\mathrm{i}}^2\right]}^{1/2}\right\}\cdot \exp \left(-\frac{{\mathsf{\phi }}_{\mathrm{Sch}}}{{\mathsf{\phi }}_{\mathrm{T}}}\right);$$

(9)

$$p=\left\{-\frac{N_{\mathrm{d}}-N_{\mathrm{a}}}{2}+{\left[{\left(\frac{N_{\mathrm{d}}-N_{\mathrm{a}}}{2}\right)}^2+n_{\mathrm{i}}^2\right]}^{1/2}\right\}\cdot \exp \left(\frac{{\mathsf{\phi }}_{\mathrm{Sch}}}{{\mathsf{\phi }}_{\mathrm{T}}}\right);$$

(10)

$$\mathsf{\phi }={\mathsf{\phi }}_{\mathrm{T}}\cdot \ln \left(\frac{n}{n_{\mathrm{i}}}\right)+U(t)-{\mathsf{\phi }}_{\mathrm{Sch}}=-{\mathsf{\phi }}_{\mathrm{T}}\cdot \ln \left(\frac{p}{n_{\mathrm{i}}}\right)+U(t)-{\mathsf{\phi }}_{\mathrm{Sch}};$$

(11)

$$n_{\mathrm{ph}}=0,$$

(12)

where ${\mathsf{\phi }}_{\mathrm{Sch}}$ is the Schottky barrier height (at ohmic contacts ${\mathsf{\phi }}_{\mathrm{Sch}}=0$); $U(t)$ is the bias voltage applied to a contact at the moment $t$. Equations (9)–(11) were derived under the assumptions of thermodynamic equilibrium and infinite recombination velocity at a contact.

If the fixed current mode is considered, the boundary conditions in Equations (9), (10), and (12) remain unchanged, and Equation (11) is replaced by the following formula [22]:

$$\frac{\partial \mathsf{\phi }}{\partial \xi }=\left[\frac{j(t)}{q}+{\mathsf{\mu }}_{\mathrm{n}}\left({\mathsf{\phi }}_{\mathrm{T}}\cdot \frac{\partial n}{\partial \xi }-n\frac{\partial V_{\mathrm{n}}}{\partial \xi }\right)-{\mathsf{\mu }}_{\mathrm{p}}\left({\mathsf{\phi }}_{\mathrm{T}}\cdot \frac{\partial p}{\partial \xi }-p\frac{\partial V_{\mathrm{p}}}{\partial \xi }\right)\right]\cdot {\left({\mathsf{\mu }}_{\mathrm{n}}n+{\mathsf{\mu }}_{\mathrm{p}}p\right)}^{-1},$$

(13)

where $\xi$ is the normal to the surface of a contact; $j(t)$ is the density of current flowing through a contact at the moment $t$.

At the contact-free boundaries of a device, the following Neumann conditions are set [22]:

$$\frac{\partial n}{\partial \xi }=0;\frac{\partial p}{\partial \xi }=0;\frac{\partial \mathsf{\phi }}{\partial \xi }=0;\frac{\partial n_{\mathrm{ph}}}{\partial \xi }=0.$$

(14)

To obtain the proper initial conditions for the non-stationary drift-diffusion model, the stationary formulation of Equations (1)–(6) should be solved numerically.

In this research, we apply the following analytical models for the calculation of charge carrier mobilities ${\mathsf{\mu }}_{\mathrm{n}}$ and ${\mathsf{\mu }}_{\mathrm{p}}$:

• The Caughey–Thomas model, which takes into account the dependence of electron and hole mobilities on the lattice temperature $T$ and impurity concentrations $N_{\mathrm{d}}$ and $N_{\mathrm{a}}$ [20]:

  $${\mathsf{\mu }}_{c,\mathrm{LF}}={\mathsf{\mu }}_{c,\min }+\frac{{\mathsf{\mu }}_{c,\min ,300}\cdot {\left(T/300\mathrm{K}\right)}^{\gamma }-{\mathsf{\mu }}_{c,\min }}{1+{\left[\left(N_{\mathrm{d}}+N_{\mathrm{a}}\right)/N_{c,\mathrm{ref}}\right]}^{{\mathsf{\alpha }}_c}},$$

  (15)

  where $c$ is the type of conductivity ($n$ or $p$); ${\mathsf{\mu }}_{c,\min }$, ${\mathsf{\mu }}_{c,\min ,300}$, $\gamma$, $N_{c,\mathrm{ref}}$, ${\mathsf{\alpha }}_c$ are the material-dependent fitting parameters;
• The model of high-field hole mobility, which describes the effect of hole drift velocity saturation [20]:

  $${\mathsf{\mu }}_{\mathrm{p}}={\mathsf{\mu }}_{\mathrm{p},\mathrm{LF}}\cdot {\left[1+{\left(\frac{{\mathsf{\mu }}_{\mathrm{p},\mathrm{LF}}\cdot E_{\mathrm{p}}}{v_{\mathrm{p},\mathrm{sat}}}\right)}^{{\mathsf{\beta }}_{\mathrm{p}}}\right]}^{-1/{\mathsf{\beta }}_{\mathrm{p}}};$$

  (16)

  $$E_{\mathrm{p}}=\left|\nabla {\mathsf{\phi }}_{\mathrm{p}}\right|=\left|\nabla \left(\mathsf{\phi }-V_{\mathrm{p}}\right)+\frac{1}{p}\cdot \nabla \left({\mathsf{\phi }}_{\mathrm{T}}\cdot p\right)\right|,$$

  (17)

  where ${\mathsf{\mu }}_{\mathrm{p},\mathrm{LF}}$ is the Caughey–Thomas mobility of holes; $E_{\mathrm{p}}$ is the gradient of quasi-Fermi level for holes; $v_{\mathrm{p},\mathrm{sat}}$ is the saturation velocity of holes; ${\mathsf{\beta }}_{\mathrm{p}}$ is the material-dependent fitting parameter;
• The model of high-field electron mobility, which addresses the overshoot and saturation of field-velocity characteristics for electrons in A^IIIB^V semiconductors [20]:

  $${\mathsf{\mu }}_{\mathrm{n}}=\left({\mathsf{\mu }}_{\mathrm{n},\mathrm{LF}}+v_{\mathrm{n},\mathrm{sat}}\cdot \frac{E_{\mathrm{n}}^{{\mathsf{\beta }}_{\mathrm{n}}-1}}{E_0^{{\mathsf{\beta }}_{\mathrm{n}}}}\right)\cdot {\left(1+\frac{E_n^{{\mathsf{\beta }}_{\mathrm{n}}}}{E_0^{{\mathsf{\beta }}_{\mathrm{n}}}}\right)}^{-1};$$

  (18)

  $$E_{\mathrm{n}}=\left|\nabla {\mathsf{\phi }}_{\mathrm{n}}\right|=\left|\nabla \left(\mathsf{\phi }+V_{\mathrm{n}}\right)-\frac{1}{n}\cdot \nabla \left({\mathsf{\phi }}_{\mathrm{T}}\cdot n\right)\right|,$$

  (19)

  where ${\mathsf{\mu }}_{\mathrm{n},\mathrm{LF}}$ is the Caughey–Thomas mobility of electrons; $E_{\mathrm{n}}$ is the gradient of quasi-Fermi level for electrons; $v_{\mathrm{n},\mathrm{sat}}$ is the saturation velocity of electrons in the upper valley; ${\mathsf{\beta }}_{\mathrm{n}}$ is the material-dependent fitting parameter; $E_0$ is the critical field intensity, at which the electron drift velocity reaches the maximum.

It is worth noting that the mobility models in Equations (15)–(19) are phenomenological. It means that the models are suitable for the adequate estimation of high-field effects, but their accurate simulation requires the experimentally measured mobility or application of the hydrodynamic approach.

We take into account the following processes of carrier recombination in high-speed A^IIIB^V optoelectronic devices:

• The stimulated radiative recombination in a laser [19]:

  $$R_{\mathrm{st}}=v_{\mathrm{g}}\cdot g\cdot n_{\mathrm{ph}};$$

  (20)
• The spontaneous radiative recombination in a laser [19]:

  $$R_{\mathrm{sp}}=\frac{{\left(np-n_{\mathrm{i}}^2\right)}^{1/2}}{{\mathsf{\tau }}_{\mathrm{S}}};$$

  (21)
• The Auger recombination [20]:

  $$R_{\mathrm{aug}}=\left(C_{\mathrm{n},\mathrm{aug}}\cdot n+C_{\mathrm{p},\mathrm{aug}}\cdot p\right)\cdot \left(np-n_{\mathrm{i}}^2\right),$$

  (22)

  where $C_{\mathrm{n},\mathrm{aug}}$, $C_{\mathrm{p},\mathrm{aug}}$ are the material-dependent Auger coefficients for electrons and holes;
• The direct recombination in materials with a direct bandgap [20]:

  $$R_{\mathrm{dir}}=C_{\mathrm{dir}}\cdot \left(np-n_{\mathrm{i}}^2\right),$$

  (23)

  where $C_{\mathrm{dir}}$ is the material-dependent coefficient of the direct recombination.

We neglect the Shockley–Read–Hall recombination because its inertness is more durable than the time intervals being considered [20]. Furthermore, in direct bandgap materials, that process dominates only under conditions of very low carrier density [23].

The following analytical model is usually utilized for the calculation of the optical generation rate $G_{\mathrm{opt}}$ in resonant-cavity-enhanced (RCE) photodetectors [21]:

$$G_{\mathrm{opt}}(t)=\frac{P_{\mathrm{in}}(t)\cdot Q}{V_{\mathrm{res}}\cdot E_{\mathrm{ph}}};$$

(24)

$$Q=\left(1-R_1\right)\cdot \left[1-\exp \left(-\mathsf{\alpha }L\right)\right]\cdot \frac{1+R_2\cdot \exp \left(-\mathsf{\alpha }L\right)}{{\left[1-\sqrt{R_1{R}_2}\cdot \exp \left(-\mathsf{\alpha }L\right)\right]}^2},$$

(25)

where $P_{\mathrm{in}}(t)$ is the power of incident optical radiation at the moment $t$; $Q$ is the quantum efficiency of a device; $V_{\mathrm{res}}$ is the volume of the absorption region of a photodetector; $E_{\mathrm{ph}}$ is the photon energy equal to the bandgap energy of the absorption region; $R_1$, $R_2$ are the reflection coefficients of the semi- and totally-reflecting mirrors, which form a resonant cavity; $\mathsf{\alpha }$ is the absorption coefficient of a semiconductor; $L$ is the length of a resonator. Equation (25) is valid at resonant wavelengths $\lambda =2L_{\mathrm{res}}/m$, where $m$ is a positive integer and gives the peak value of quantum efficiency. The model in Equations (24) and (25) assumes the instant propagation of light throughout a resonant cavity and does not take into account the coordinate variation of the optical generation rate. In two- and three-dimensional cases, the last feature leads to the overestimated results of photoresponse calculation.

In this paper we propose an improved model for the calculation of the optical generation rate in RCE photodetectors.

If a monochromatic light propagates through a semiconductor material, the absorption is described by the Beer–Lambert–Bouguer law [24]:

$$I(y)=I_0\exp \left(-\mathsf{\alpha }y\right),$$

(26)

where $y$ is the thickness of the absorbing layer; $I_0$ is the intensity of the incoming light beam; $I$ is the intensity of an outgoing light beam; $\mathsf{\alpha }$ is the absorption coefficient of a semiconductor at the wavelength of incident light.

Here we consider the waveguide p-i-n RCE photodetector shown in Figure 1. The resonant cavity is formed by semi- and totally-reflecting mirrors (positions 2 and 4). The laser beam (position 1) with time-dependent power $P_{\mathrm{in}}(t)$ falls orthogonally to the surface of semi-reflecting mirror (position 2). The $y$ coordinate axis is directed along the incident laser beam and has the reference point near the semi-reflecting mirror. Outside of the resonant cavity, the photon density $n_{\mathrm{in}}(t)$ is given by the following equation:

$$n_{\mathrm{in}}(t)=\frac{P_{\mathrm{in}}(t)}{c{E}_{\mathrm{ph}}HW},$$

(27)

where $c$ is the speed of light in vacuum; $E_{\mathrm{ph}}$ is the photon energy equal to the bandgap energy of the absorption region; $W$, $H$ are the width and height of the resonant cavity or absorption region, as shown in Figure 1. After the passing of light through the semi-reflecting mirror, the photon density $n_{\mathrm{ph}}(y,t)$ at the distance $y$ and the moment $t$ is defined by the following formula resulting from Equation (26):

$$n_{\mathrm{ph}}(y,t)=P_{\mathrm{in}}\left(t-\frac{y}{c/n_{\mathrm{r}}}\right)\cdot \frac{1-R_1}{c{E}_{\mathrm{ph}}HW}\cdot \exp \left(-\mathsf{\alpha }y\right),$$

(28)

where $n_{\mathrm{r}}$ is the refractive index of the absorption region; $\frac{y}{c/n_{\mathrm{r}}}$ is the instant, at which photons overcome the $y$ distance in the resonant cavity. The reduction of photon density caused by the absorption and $M$-fold reflection of light from the resonator mirrors is described by the following equation:

$$n_{\mathrm{ph}}(y,t)={\displaystyle \sum _{m=1}^M\left[P_{\mathrm{in}}\left(t-\frac{2L\lfloor m/2\rfloor +{\left(-1\right)}^{m+1}y}{c/n_{\mathrm{r}}}\right)\cdot \frac{1-R_1}{c{E}_{\mathrm{ph}}HW}\cdot R_1^{\lfloor \frac{m-1}{2}\rfloor }{R}_2^{\lfloor \frac{m}{2}\rfloor }\cdot \exp \left(-2\lfloor \frac{m}{2}\rfloor \mathsf{\alpha }L+{\left(-1\right)}^m\mathsf{\alpha }y\right)\right]},$$

(29)

where $M=\lceil \frac{c\cdot t}{n_{\mathrm{r}}\cdot L}\rceil$ is the multiplicity of photon beam reflection; $\lfloor \rfloor$, $\lceil \rceil$ are the floor and ceiling functions. In the case of one-dimensional simulation along the $x$ axis, the total optical power $P_{\mathrm{abs}}$ absorbed in the whole resonant cavity and optical generation rate $G_{\mathrm{opt}}$ are given by the following formulas:

$$P_{\mathrm{abs}}(t)=\frac{E_{\mathrm{ph}}HW}{\Delta t}\cdot \left[1-\exp \left(-\mathsf{\alpha }\frac{c}{n_{\mathrm{r}}}\Delta t\right)\right]\cdot {\displaystyle \underset{0}{\overset{L}{\int }}{n}_{\mathrm{ph}}(y,t)dy};$$

(30)

$$G_{\mathrm{opt}}(t)=\frac{P_{\mathrm{abs}}(t)}{L\cdot W\cdot H\cdot E_{\mathrm{ph}}},$$

(31)

where $\Delta t$ is the time step. Otherwise, ${\displaystyle \underset{0}{\overset{L}{\int }}{n}_{\mathrm{ph}}(y,t)dy}$ in Equation (30) should be replaced by $n_{\mathrm{ph}}(y,t)\cdot \Delta y$, where $\Delta y$ is the step of coordinate grid along the $y$ axis, and $L\cdot W\cdot H$ in Equation (31) should be replaced by $\Delta y\cdot W\cdot \Delta x$, where $\Delta x$ is the step of coordinate grid along the $x$ axis.

The model in Equations (29)–(31) is suitable for one- and two-dimensional non-stationary simulation of high-speed RCE photodetectors and takes into account the dynamics of photon propagation in the resonant cavity. The model does not consider wave effects, but the introduction of special fitting parameters in Equation (30) allows us to fit the model to the experimental data.

3. Simulation Techniques

The implementation of the model in Equations (1)–(6), with the boundary conditions in Equations (9)–(14) and the additional parameter models in Equations (7), (8), (15)–(31), requires the application of numerical methods.

To provide a stable and convergent numerical solution of a system of differential equations, an efficient discretization scheme and a high-performance algorithm are needed. These application-oriented components form a technique of numerical simulation. The development of a simulation technique implies the application of various mathematical methods and approaches, the design of applied software, and the conduction of computational experiments.

Nowadays, a lot of problems in the field of computational electronics remains unsolved. In our opinion, the universal approach to the efficient implementation of the drift-diffusion models has not been developed yet. The non-stationary simulation of heterostructural devices is characterized by the following challenges:

• The instability of discretization schemes and algorithms;
• Significant time and resource consumption;
• The inadequacy of simulation results caused by calculation errors.

Sometimes, it is quite difficult to deal with these challenges due to the features of the drift-diffusion model, the variety of parameter values, and the lack of appropriate mathematical investigations.

In this paper, we propose two techniques designed for the numerical simulation of high-speed optoelectronic devices within the drift-diffusion approximation. Both techniques are based on the finite difference method [25], in which derivatives are approximated by the finite difference expressions. The first technique applies the explicit, upwind, and Gummel’s methods and is suitable for one- and two-dimensional non-stationary simulation of lasers and photodetectors. The second one is based on the implicit, Newton’s, and Slotboom’s methods and is suitable for two-dimensional stationary simulation of lasers and photodetectors and one-dimensional non-stationary simulation of photosensitive devices. The distinctive features of the proposed simulation techniques are discussed below.

3.1. Explicit Technique

The explicit simulation technique applies the combined variable base of $\left\{n,p,\mathsf{\phi },F_{\mathrm{n}},F_{\mathrm{p}}\right\}$, where $F_{\mathrm{n}}=\exp \left(-{\mathsf{\phi }}_{\mathrm{n}}\right)$ and $F_{\mathrm{p}}=\exp \left({\mathsf{\phi }}_{\mathrm{p}}\right)$. It is necessary for the reduction of computational error caused by the difference between the magnitude orders of carrier densities and electrostatic potential. The transformation equations for electron and hole concentrations have the following form [26]:

$$n=n_{\mathrm{i}}\cdot F_{\mathrm{n}}\cdot \exp \left(\frac{\mathsf{\phi }+V_{\mathrm{n}}}{{\mathsf{\phi }}_{\mathrm{T}}}\right);$$

(32)

$$p=n_{\mathrm{i}}\cdot F_{\mathrm{p}}\cdot \exp \left(\frac{-\mathsf{\phi }+V_{\mathrm{p}}}{{\mathsf{\phi }}_{\mathrm{T}}}\right).$$

(33)

Current density Equations (1) and (2) are considered in terms of the Slotboom’s variable base of $\left\{F_{\mathrm{n}},F_{\mathrm{p}},\mathsf{\phi }\right\}$ [27], and Equations (3)–(5) are considered in terms of the standard variable base of $\left\{n,p,\mathsf{\phi }\right\}$. After the substitution of Equations (1) and (2) in Equations (3) and (4) and the scaling of all equations by the standard factors [17,18], the drift-diffusion formulation in terms of $\left\{n,p,\mathsf{\phi },F_{\mathrm{n}},F_{\mathrm{p}}\right\}$ is given by

$$\frac{\partial n}{\partial t}=\nabla \left[{\mathsf{\mu }}_{\mathrm{n}}\cdot \exp \left(\mathsf{\phi }+V_{\mathrm{n}}\right)\cdot \nabla F_{\mathrm{n}}\right]+G-R;$$

(34)

$$\frac{\partial p}{\partial t}=\nabla \left[{\mathsf{\mu }}_{\mathrm{p}}\cdot \exp \left(-\mathsf{\phi }+V_{\mathrm{p}}\right)\cdot \nabla F_{\mathrm{p}}\right]+G-R;$$

(35)

$$\nabla \left(\mathsf{\epsilon }\cdot \nabla \mathsf{\phi }\right)=n-p-N_{\mathrm{d}}+N_{\mathrm{a}}.$$

(36)

To discretize Equations (34) and (35), we use the explicit difference scheme [28] for the time derivatives, and the first-order upwind scheme [29] for the coordinate derivatives. Equation (36) is discretized as usual. In the two-dimensional case, the resulting finite difference equations have the following form:

$$\begin{array}{c}{n}_{i,j}^{k+1}=n_{i,j}^k+\frac{2\Delta t^k}{\Delta x_i+\Delta x_{i-1}}\{[{\mathsf{\mu }}_{\mathrm{n},i+1,j}^k\cdot \exp \left({\mathsf{\phi }}_{i+1,j}^k+V_{\mathrm{n},i+1,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{n},i+1,j}^k-F_{\mathrm{n},i,j}^k\right)+{\mathsf{\mu }}_{\mathrm{n},i,j}^k\exp \left({\mathsf{\phi }}_{i,j}^k+V_{\mathrm{n},i,j}\right)\cdot \\ \mathsf{\theta }\left(F_{\mathrm{n},i,j}^k-F_{\mathrm{n},i+1,j}^k\right)]\frac{F_{\mathrm{n},i+1,j}^k-F_{\mathrm{n},i,j}^k}{\Delta x_i}-[{\mathsf{\mu }}_{\mathrm{n},i,j}^k\exp \left({\mathsf{\phi }}_{i,j}^k+V_{\mathrm{n},i,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{n},i,j}^k-F_{\mathrm{n},i-1,j}^k\right)+{\mathsf{\mu }}_{\mathrm{n},i-1,j}^k\cdot \\ \exp \left({\mathsf{\phi }}_{i-1,j}^k+V_{\mathrm{n},i-1,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{n},i-1,j}^k-F_{\mathrm{n},i,j}^k\right)]\frac{F_{\mathrm{n},i,j}^k-F_{\mathrm{n},i-1,j}^k}{\Delta x_{i-1}}\}+\frac{2\Delta t^k}{\Delta y_j+\Delta y_{j-1}}\{[{\mathsf{\mu }}_{\mathrm{n},i,j+1}^k\exp \left({\mathsf{\phi }}_{i,j+1}^k+V_{\mathrm{n},i,j+1}\right)\cdot \\ \mathsf{\theta }\left(F_{\mathrm{n},i,j+1}^k-F_{\mathrm{n},i,j}^k\right)+{\mathsf{\mu }}_{\mathrm{n},i,j}^k\cdot \exp \left({\mathsf{\phi }}_{i,j}^k+V_{\mathrm{n},i,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{n},i,j}^k-F_{\mathrm{n},i,j+1}^k\right)]\cdot \frac{F_{\mathrm{n},i,j+1}^k-F_{\mathrm{n},i,j}^k}{\Delta y_j}-\\ [{\mathsf{\mu }}_{\mathrm{n},i,j}^k\cdot \exp \left({\mathsf{\phi }}_{i,j}^k+V_{\mathrm{n},i,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{n},i,j}^k-F_{\mathrm{n},i,j-1}^k\right)+{\mathsf{\mu }}_{\mathrm{n},i,j-1}^k\cdot \exp \left({\mathsf{\phi }}_{i,j-1}^k+V_{\mathrm{n},i,j-1}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{n},i,j-1}^k-F_{\mathrm{n},i,j}^k\right)]\cdot \\ \frac{F_{\mathrm{n},i,j}^k-F_{\mathrm{n},i,j-1}^k}{\Delta y_{j-1}}\}+\left(G_{i,j}^k-R_{i,j}^k\right)\cdot \Delta t^k;\end{array}$$

(37)

$$\begin{array}{c}{p}_{i,j}^{k+1}=p_{i,j}^k+\frac{2\Delta t^k}{\Delta x_i+\Delta x_{i-1}}\{[{\mathsf{\mu }}_{\mathrm{p},i+1,j}^k\cdot \exp \left(-{\mathsf{\phi }}_{i+1,j}^k+V_{\mathrm{p},i+1,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{p},i+1,j}^k-F_{\mathrm{p},i,j}^k\right)+{\mathsf{\mu }}_{\mathrm{p},i,j}^k\exp \left(-{\mathsf{\phi }}_{i,j}^k+V_{\mathrm{p},i,j}\right)\cdot \\ \mathsf{\theta }\left(F_{\mathrm{p},i,j}^k-F_{\mathrm{p},i+1,j}^k\right)]\frac{F_{\mathrm{p},i+1,j}^k-F_{\mathrm{p},i,j}^k}{\Delta x_i}-[{\mathsf{\mu }}_{\mathrm{p},i,j}^k\exp \left(-{\mathsf{\phi }}_{i,j}^k+V_{\mathrm{p},i,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{p},i,j}^k-F_{\mathrm{p},i-1,j}^k\right)+{\mathsf{\mu }}_{\mathrm{p},i-1,j}^k\cdot \\ \exp \left(-{\mathsf{\phi }}_{i-1,j}^k+V_{\mathrm{p},i-1,j}\right)\mathsf{\theta }\left(F_{\mathrm{p},i-1,j}^k-F_{\mathrm{p},i,j}^k\right)]\frac{F_{\mathrm{p},i,j}^k-F_{\mathrm{p},i-1,j}^k}{\Delta x_{i-1}}\}+\frac{2\Delta t^k}{\Delta y_j+\Delta y_{j-1}}\{[{\mathsf{\mu }}_{\mathrm{p},i,j+1}^k\exp \left(-{\mathsf{\phi }}_{i,j+1}^k+V_{\mathrm{p},i,j+1}\right)\cdot \\ \mathsf{\theta }\left(F_{\mathrm{p},i,j+1}^k-F_{\mathrm{p},i,j}^k\right)+{\mathsf{\mu }}_{\mathrm{p},i,j}^k\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^k+V_{\mathrm{p},i,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{p},i,j}^k-F_{\mathrm{p},i,j+1}^k\right)]\cdot \frac{F_{\mathrm{p},i,j+1}^k-F_{\mathrm{p},i,j}^k}{\Delta y_j}-\\ [{\mathsf{\mu }}_{\mathrm{p},i,j}^k\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^k+V_{\mathrm{p},i,j}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{p},i,j}^k-F_{\mathrm{p},i,j-1}^k\right)+{\mathsf{\mu }}_{\mathrm{p},i,j-1}^k\cdot \exp \left(-{\mathsf{\phi }}_{i,j-1}^k+V_{\mathrm{p},i,j-1}\right)\cdot \mathsf{\theta }\left(F_{\mathrm{p},i,j-1}^k-F_{\mathrm{p},i,j}^k\right)]\cdot \\ \frac{F_{\mathrm{p},i,j}^k-F_{\mathrm{p},i,j-1}^k}{\Delta y_{j-1}}\}+\left(G_{i,j}^k-R_{i,j}^k\right)\cdot \Delta t^k;\end{array}$$

(38)

$$\begin{array}{c}\frac{1}{\Delta x_i+\Delta x_{i-1}}\cdot \left[\left({\mathsf{\epsilon }}_{i+1,j}+{\mathsf{\epsilon }}_{i,j}\right)\cdot \frac{{\mathsf{\phi }}_{i+1,j}^{k+1}-{\mathsf{\phi }}_{i,j}^{k+1}}{\Delta x_i}-\left({\mathsf{\epsilon }}_{i,j}+{\mathsf{\epsilon }}_{i-1,j}\right)\cdot \frac{{\mathsf{\phi }}_{i,j}^{k+1}-{\mathsf{\phi }}_{i-1,j}^{k+1}}{\Delta x_{i-1}}\right]+\\ \frac{1}{\Delta y_j+\Delta y_{j-1}}\cdot \left[\left({\mathsf{\epsilon }}_{i,j+1}+{\mathsf{\epsilon }}_{i,j}\right)\cdot \frac{{\mathsf{\phi }}_{i,j+1}^{k+1}-{\mathsf{\phi }}_{i,j}^{k+1}}{\Delta y_j}-\left({\mathsf{\epsilon }}_{i,j}+{\mathsf{\epsilon }}_{i,j-1}\right)\cdot \frac{{\mathsf{\phi }}_{i,j}^{k+1}-{\mathsf{\phi }}_{i,j-1}^{k+1}}{\Delta y_{j-1}}\right]=n_{i,j}^{k+1}-p_{i,j}^{k+1}-N_{\mathrm{d},i,j}+N_{\mathrm{a},i,j},\end{array}$$

(39)

where $i$, $j$ are the indexes of coordinate grid points on the y and x axes; $k$ is the index of time grid points; $\Delta x_i=x_{i+1}-x_i$ is the i-th step of coordinate grid along the x axis; $\Delta y_j=y_{j+1}-y_j$ is the j-th step of coordinate grid along the y axis; $\Delta t^k$ is the k-th step of time grid.

To solve the system formed by finite-difference Equations (37)–(39), we use the modified Gummel’s iterative method [30]. At each time step we calculate carrier concentrations according to Equations (37) and (38), and then we solve the system of linear algebraic equations (SLAE) in order to compute the distribution of electrostatic potential. Carrier mobilities, and generation and recombination rates for the next instant are estimated after the solution of the Poisson’s equation at the actual moment of time. In paper [21] we showed that the application of Gummel’s method together with the explicit and upwind difference schemes allows for the reduction in consumption of computational resources, but requires a shorter time step than other approaches.

If the explicit technique is utilized for the simulation of an injection laser, Equation (6) is added to the basic drift-diffusion formulation. After the scaling and explicit discretization, the rate equation for photons has the following form:

$$n_{\mathrm{ph},i,j}^{k+1}=n_{\mathrm{ph},i,j}^k+\Delta t^k\left(-\frac{n_{\mathrm{ph},i,j}^k}{{\mathsf{\tau }}_{\mathrm{ph}}}+\mathsf{\beta }\cdot \exp \left(\frac{V_{\mathrm{n},i,j}+V_{\mathrm{p},i,j}}{2}\right)\sqrt{F_{\mathrm{n},i,j}^k{F}_{\mathrm{p},i,j}^k-1}+v_{\mathrm{g}}{g}_{i,j}^k{n}_{\mathrm{ph},i,j}^k\right).$$

(40)

The calculation of photon density is performed after the solution of the Poisson’s equation at each time step.

3.2. Implicit Technique

The implicit technique of the drift-diffusion numerical simulation applies the Slotboom‘s variable base of $\left\{F_{\mathrm{n}},F_{\mathrm{p}},\mathsf{\phi }\right\}$, which allows us to reduce the computational error and calculate the adequate coordinate distributions of current density. The technique requires the setting of boundary conditions for $F_{\mathrm{n}}$ and $F_{\mathrm{p}}$, which are given by the following equations:

• At ohmic and Schottky contacts:

  $$F_{\mathrm{n}}=\exp \left(\frac{-U(t)+{\mathsf{\phi }}_{\mathrm{Sch}}}{{\mathsf{\phi }}_{\mathrm{T}}}\right);$$

  (41)

  $$F_{\mathrm{p}}=\exp \left(\frac{U(t)-{\mathsf{\phi }}_{\mathrm{Sch}}}{{\mathsf{\phi }}_{\mathrm{T}}}\right);$$

  (42)
• At contact-free boundaries:

  $$\frac{\partial F_{\mathrm{n}}}{\partial \xi }=0;\frac{\partial F_{\mathrm{p}}}{\partial \xi }=0.$$

  (43)

The discretization scheme is based on the implicit method [28] and central differences. In the two-dimensional case, the normalized finite difference drift-diffusion equations have the following form:

$$\begin{array}{c}\frac{\Delta y_j+\Delta y_{j-1}}{\Delta x_i}\cdot \left[{\mathsf{\mu }}_{\mathrm{n},i+1,j}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i+1,j}^{k+1}+V_{\mathrm{n},i+1,j}\right)+{\mathsf{\mu }}_{\mathrm{n},i,j}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{n},i,j}\right)\right]\cdot \left(F_{\mathrm{n},i+1,j}^{k+1}-F_{\mathrm{n},i,j}^{k+1}\right)-\\ \frac{\Delta y_j+\Delta y_{j-1}}{\Delta x_{i-1}}\cdot \left[{\mathsf{\mu }}_{\mathrm{n},i,j}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{n},i,j}\right)+{\mathsf{\mu }}_{\mathrm{n},i-1,j}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i-1,j}^{k+1}+V_{\mathrm{n},i-1,j}\right)\right]\cdot \left(F_{\mathrm{n},i,j}^{k+1}-F_{\mathrm{n},i-1,j}^{k+1}\right)+\\ \frac{\Delta x_i+\Delta x_{i-1}}{\Delta y_j}\cdot \left[{\mathsf{\mu }}_{\mathrm{n},i,j+1}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i,j+1}^{k+1}+V_{\mathrm{n},i,j+1}\right)+{\mathsf{\mu }}_{\mathrm{n},i,j}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{n},i,j}\right)\right]\cdot \left(F_{\mathrm{n},i,j+1}^{k+1}-F_{\mathrm{n},i,j}^{k+1}\right)-\\ \frac{\Delta x_i+\Delta x_{i-1}}{\Delta y_{j-1}}\cdot \left[{\mathsf{\mu }}_{\mathrm{n},i,j}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{n},i,j}\right)+{\mathsf{\mu }}_{\mathrm{n},i,j-1}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i,j-1}^{k+1}+V_{\mathrm{n},i,j-1}\right)\right]\cdot \left(F_{\mathrm{n},i,j}^{k+1}-F_{\mathrm{n},i,j-1}^{k+1}\right)-\\ \frac{\left(\Delta x_i+\Delta x_{i-1}\right)\left(\Delta y_j+\Delta y_{j-1}\right)}{\Delta t^k}\cdot \left[F_{\mathrm{n},i,j}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{n},i,j}\right)-F_{\mathrm{n},i,j}^k\cdot \exp \left({\mathsf{\phi }}_{i,j}^k+V_{\mathrm{n},i,j}\right)\right]+\\ \left(G_{i,j}^{k+1}-R_{i,j}^{k+1}\right)\left(\Delta x_i+\Delta x_{i-1}\right)\left(\Delta y_j+\Delta y_{j-1}\right)=0;\end{array}$$

(44)

$$\begin{array}{c}\frac{\Delta y_j+\Delta y_{j-1}}{\Delta x_i}\cdot \left[{\mathsf{\mu }}_{\mathrm{p},i+1,j}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i+1,j}^{k+1}+V_{\mathrm{p},i+1,j}\right)+{\mathsf{\mu }}_{\mathrm{p},i,j}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{p},i,j}\right)\right]\cdot \left(F_{\mathrm{p},i+1,j}^{k+1}-F_{\mathrm{p},i,j}^{k+1}\right)-\\ \frac{\Delta y_j+\Delta y_{j-1}}{\Delta x_{i-1}}\cdot \left[{\mathsf{\mu }}_{\mathrm{p},i,j}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{p},i,j}\right)+{\mathsf{\mu }}_{\mathrm{p},i-1,j}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i-1,j}^{k+1}+V_{\mathrm{p},i-1,j}\right)\right]\cdot \left(F_{\mathrm{p},i,j}^{k+1}-F_{\mathrm{p},i-1,j}^{k+1}\right)+\\ \frac{\Delta x_i+\Delta x_{i-1}}{\Delta y_j}\cdot \left[{\mathsf{\mu }}_{\mathrm{p},i,j+1}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i,j+1}^{k+1}+V_{\mathrm{p},i,j+1}\right)+{\mathsf{\mu }}_{\mathrm{p},i,j}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{p},i,j}\right)\right]\cdot \left(F_{\mathrm{p},i,j+1}^{k+1}-F_{\mathrm{p},i,j}^{k+1}\right)-\\ \frac{\Delta x_i+\Delta x_{i-1}}{\Delta y_{j-1}}\cdot \left[{\mathsf{\mu }}_{\mathrm{p},i,j}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{p},i,j}\right)+{\mathsf{\mu }}_{\mathrm{p},i,j-1}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i,j-1}^{k+1}+V_{\mathrm{p},i,j-1}\right)\right]\cdot \left(F_{\mathrm{p},i,j}^{k+1}-F_{\mathrm{p},i,j-1}^{k+1}\right)-\\ \frac{\left(\Delta x_i+\Delta x_{i-1}\right)\left(\Delta y_j+\Delta y_{j-1}\right)}{\Delta t^k}\cdot \left[F_{\mathrm{p},i,j}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{p},i,j}\right)-F_{\mathrm{p},i,j}^k\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^k+V_{\mathrm{p},i,j}\right)\right]+\\ \left(G_{i,j}^{k+1}-R_{i,j}^{k+1}\right)\left(\Delta x_i+\Delta x_{i-1}\right)\left(\Delta y_j+\Delta y_{j-1}\right)=0;\end{array}$$

(45)

$$\begin{array}{c}\frac{\Delta y_j+\Delta y_{j-1}}{\Delta x_i}\cdot \left({\mathsf{\epsilon }}_{i+1,j}+{\mathsf{\epsilon }}_{i,j}\right)\cdot \left({\mathsf{\phi }}_{i+1,j}^{k+1}-{\mathsf{\phi }}_{i,j}^{k+1}\right)-\frac{\Delta y_j+\Delta y_{j-1}}{\Delta x_{i-1}}\cdot \left({\mathsf{\epsilon }}_{i,j}+{\mathsf{\epsilon }}_{i-1,j}\right)\cdot \left({\mathsf{\phi }}_{i,j}^{k+1}-{\mathsf{\phi }}_{i-1,j}^{k+1}\right)+\\ \frac{\Delta x_i+\Delta x_{i-1}}{\Delta y_j}\cdot \left({\mathsf{\epsilon }}_{i,j+1}+{\mathsf{\epsilon }}_{i,j}\right)\cdot \left({\mathsf{\phi }}_{i,j+1}^{k+1}-{\mathsf{\phi }}_{i,j}^{k+1}\right)-\frac{\Delta x_i+\Delta x_{i-1}}{\Delta y_{j-1}}\cdot \left({\mathsf{\epsilon }}_{i,j}+{\mathsf{\epsilon }}_{i,j-1}\right)\cdot \left({\mathsf{\phi }}_{i,j}^{k+1}-{\mathsf{\phi }}_{i,j-1}^{k+1}\right)+\\ \left(\Delta x_i+\Delta x_{i-1}\right)\left(\Delta y_j+\Delta y_{j-1}\right)\left[F_{\mathrm{n},i,j}^{k+1}\cdot \exp \left({\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{n},i,j}\right)-F_{\mathrm{p},i,j}^{k+1}\cdot \exp \left(-{\mathsf{\phi }}_{i,j}^{k+1}+V_{\mathrm{p},i,j}\right)-N_{\mathrm{d},i,j}+N_{\mathrm{a},i,j}\right]=0.\end{array}$$

(46)

To solve the system formed by finite difference Equations (44)–(46), we apply the standard Newton’s method [25]. At each instant of time, the continuity equations together with the Poisson’s one are solved iteratively. After the solution, carrier mobilities, and generation and recombination rates are recalculated for the next iteration step. The aforementioned approach is characterized by quadratic convergence (about several tens of iterations are necessary) but requires a lot of computational resources for the processing of large matrixes. That is why we do not recommend it for non-stationary two-dimensional simulation.

4. Simulation Results and Discussion

We developed the specialized software package for the numerical simulation of carrier transport and accumulation in high-speed on-chip optoelectronic devices based on A^IIIB^V heterostructures. The package implements the drift-diffusion model and finite difference simulation techniques discussed above. We wrote the program in the Octave programming language and executed it using the GNU Octave software. The proposed modelling aids allow the simulation of photodetectors and injection lasers with various electrophysical, constructive, and technological parameters operating at different control actions.

Figure 2a shows the cross-section of the injection laser with a functionally integrated modulator based on a double A^IIIB^V nanoheterostructure.

The laser-modulator applies the principle of the controlled spatial relocation of charge carrier density peaks within quantum regions of valence and conduction bands for the generation of amplitude- or frequency-modulated optical signals [31]. Unlike conventional laser diodes, the laser-modulator consists of heavily-doped p⁺ and n⁺ regions, as shown in positions 4 and 5 in Figure 2, with the power contacts (position 8), the amplitude or frequency modulator (position 6) with the control contacts (position 9), and the Schottky (positions 5 and 10) and p-n (positions 3 and 4) junctions. The aforementioned design provides the alteration of intensity or wavelength of laser radiation by means of the change in the transverse field applied to the control contacts (position 9). During the operation of the device, the pumping current remains unchanged.

In this paper, we considered two types of the laser-modulator. The first one includes the amplitude modulator with the energy band diagram represented in Figure 2b. The second one has the frequency modulator shown in Figure 2c. The operation principles and structures of the aforementioned types of the laser-modulator are thoroughly addressed in papers [31,32,33,34]. Since the pumping current does not change during the operation, the peak modulation frequency of the devices is not limited by slow transients in the power circuit. The performance is determined by the inertness of controlled relocation of carrier density maximums between the quantum wells of heterostructures.

Figure 3a,b show the results of drift-diffusion two-dimensional numerical simulation of transients in the laser with a functionally integrated amplitude modulator. The active region of the device is formed by the spatially displaced quantum wells in conduction, as shown in positions 11 and 12 in Figure 2b, and valence (positions 12 and 13) bands. The configuration of the wells is such that they intersect each other at the thin narrow-bandgap region (position 12). In the operation mode being considered, the change of control voltage at constant pumping current leads to the spatial relocation of electron and hole density peaks within the quantum wells. During that process, the total number of charge carriers in the wells remains nearly unchanged. If the maximums of electron and hole densities are shifted to the overlap region (position 12) by the control electric field, the raising of the laser radiation intensity occurs. On the contrary, the relocation of carrier density peaks to the space division regions (positions 11 and 13), resulting in the reduction of intensity of generated optical radiation. The curves in Figure 3a demonstrate the dynamics of electron and hole concentrations in the region of quantum well overlap (position 12).

In Figure 3a,b, the series of 5-ps pulses is applied to the control contacts at the instant of 60 ps. Before that, the supply transient occurred due to the application of the fixed pumping current to the power contacts. The variation of control voltage leads to the relocation of electron and hole density maximums within the potential wells of the heterostructure, as shown in Figure 3a. This process results in the change of photon linear density represented in Figure 3b. According to the results, the duration of electron and hole relocations are very short and do not exceed 0.1 ps. However, the response of photon density is 50 times slower than the carrier one. This feature is caused by the durable lifetime of photons in the resonant cavity of the laser-modulator (${\mathsf{\tau }}_{\mathrm{ph}}\approx 3\mathrm{ps}$). The shortening of ${\mathsf{\tau }}_{\mathrm{ph}}$ is very problematic because it leads to the growth of the threshold pumping current density and degradation of some other parameters.

Figure 3c,d show the results of time-domain drift-diffusion two-dimensional simulation of the laser with a functionally integrated frequency modulator. The active region of the device consists of two quantum wells in the valence band, as shown in positions 14 and 16 in Figure 2c, and one quantum well in the conduction band (position 15). The quantum wells in the valence band have different depths depending on b and c coefficients of GaSb_bAs_1-b and GaSb_cAs_1-c solid solutions. Regions 14, 15, and 16, as shown in Figure 2c, form two heterojunctions with staggered gaps (type II) and different bandgap energies at boundaries 14-15 (boundary between regions 14 and 15) and 15-16. The wavelength of generated laser radiation is specified by the direction of the control electric field. At one direction, the peaks of electron and hole densities are relocated to the 14-15 boundary, and the laser-modulator generates the radiation with ${\lambda }_1$ wavelength. At the opposite direction, the maximums are overlapped at the 15-16 boundary, and the device generates the radiation with ${\lambda }_2$ wavelength. Analogously to the laser with amplitude modulation, the total number of charge carriers in the quantum wells remains nearly unchanged during the device operation.

In Figure 3c,d, the interval between the reference time and the instant of 93 ps corresponds to the transient process in the power-supply circuit of the laser-modulator. The pumping current density of 3∙10⁵ A/cm² and zero control voltage remain unchanged during the supply transient. Operation mode of the device starts at 93 ps and lasts until 137 ps. During this period, we apply the rectangular voltage pulses with the amplitude of 1 V to the control contacts. The duration of each pulse is 10 ps. The pumping current density is not changed. The time of the carrier density peak relocation in the quantum wells of the laser-modulator takes less than 0.1 ps, as shown in Figure 3c. This process is controlled by the transverse electric field. Therefore, the relocation is a rapid process which does not limit the modulation frequency of laser radiation. The simulation results demonstrate that the aforementioned characteristic is affected by the photon lifetime in the resonant cavity of the laser-modulator. In the case being considered, the photon lifetime is the same as for the device with the amplitude modulation. Consequently, the duration of the modulating pulse of control voltage should be more than 3 ps in order to provide the adequate modulation coefficient.

Figure 4 shows the results of the application of the model in Equations (29)–(31) for the simulation of transients in the GaAs p-i-n RCE photodiode, as shown in Figure 1. The device is illuminated by 0.2-ps rectangular optical pulses at the moments of 0 and 0.4 ps. The maximum power of the laser radiation outside the resonator is about 1 mW. The reflection coefficients of the semi- and totally-reflecting mirrors are 0.5 and 0.99. The length, height, and width of the resonant cavity are 4000, 200, and 1000 nm.

Figure 4a–c demonstrate the dynamics of photon propagation through the resonant cavity of the photodiode, which is taken into account by the model in Equations (29)–(31). In Figure 4a,b, the first pass of photons through the resonant cavity is shown. Figure 4c represents the near-stationary state after the end of resonator transient. In this state, the further change of photon concentration is insignificant. In Figure 4d, the time dependence of the optical power absorbing by the whole resonant cavity is shown. According to the curve, the duration of the transient caused by the finite velocity of photons in the resonator is about 0.1 ps. If the input optical power grows gradually due to the inertness of the laser-modulator, the resonant cavity transient is more durable and significantly affects the performance of the photodetector.

Figure 5 represents the results of non-stationary drift-diffusion simulation of the uni-travelling-carrier (UTC) photodiode [35] based on InP/In_0.53Ga_0.47As heterojunctions. UTC photodiodes utilize only electrons as active carriers and, in general, provide better performance (about 1 ps and less according to [36]) than conventional p-i-n photodetectors. That is why UTC diodes are promising devices for the operation as parts of on-chip optical interconnections together with the high-speed lasers-modulators. In Figure 5, the lengths of the p-type absorption and i-type collection layers are 150 and 200 nm. The device is illuminated by 5-ps rectangular laser pulse at the initial time moment. The maximum power of the laser radiation outside the resonant cavity is about 0.3 mW. The photodetector operates in the fixed voltage mode, and the absolute value of the reverse bias voltage equals 0.5 V. The length of the resonant cavity is 1000 nm. The light propagates orthogonally to the direction of the built-in electric field (analogously to the p-i-n photodetector shown in Figure 1).

Figure 5a–c show the energy band diagram of the UTC photodiode being considered and spatial distributions of electric field intensity and of carrier concentrations in its structure. The results are computed at the moment of 4.8 ps after the start of illumination. The graphs demonstrate the effect of electric field screening, which occurs in the absorption region of the photodetector at quite high powers of absorbing laser radiation. If the device is illuminated during several picoseconds and the carrier leaving is sufficiently slow, the significant quantities of photogenerated electrons and holes are accumulated in the absorption region, as shown in Figure 5c. The presence of excess charge carriers causes the damping of the built-in electric field in the device according to Figure 5b. The screening effect leads to the growth of carrier mobility and to the reduction of drift carrier transport within the absorption region.

Figure 5d demonstrate the photoresponse of the aforementioned UTC photodetector to the illumination by the rectangular laser pulse. The front edge of the photocurrent pulse can be divided into two intervals. The first one is characterized by a steep slope and lasts from 0 till approximately 1 ps. The second interval has a shallow slope and lasts from 1 ps till the end of illumination at 5 ps. During the second interval, carrier transport in the UTC photodiode is affected by the electric field screening. The back edge of the photocurrent pulse is slower than the front one due to the aforementioned screening effect. Several picoseconds are necessary for the clean-out of excess photogenerated carriers and for the recovery of the built-in electric field.

5. Conclusions

In this paper we proposed the complex model designed for one- and two-dimensional simulation of the transients and stationary states in A^IIIB^V high-speed optoelectronic devices. The model is based on the drift-diffusion approximation of the semiclassical approach to the description of carrier transport and accumulation in semiconductor devices. It includes the basic drift-diffusion formulation, rate equations for photons in an injection laser, Dirichlet and Neumann boundary conditions for the fixed voltage and current modes, and additional analytical models for carrier mobilities, and generation and recombination rates.

To implement the model, we developed two techniques of computer-aided numerical simulation. The first one exploits the explicit, first-order upwind, and Gummel’s methods. It is characterized by resource efficiency, but requires a very short time step for the stable solution of the equations. We applied the explicit technique for two-dimensional non-stationary simulation of the lasers-modulators and photodetectors. The second technique is based on the implicit difference scheme and Newton’s method. It provides the quadratic convergence of a numerical solution, but takes a lot of computational resources due to the simultaneous solution of three differential equations. We used the implicit technique for the computation of initial conditions and one-dimensional time-domain simulation of the high-speed photodetectors.

We developed a specialized software package for the implementation of the proposed models and numerical simulation techniques. It allows for the simulation of injection lasers and photodetectors with various electrophysical, constructive, and technological parameters operating at different control actions.

The proposed simulation aids can be used for the development of new optoelectronic devices, estimation of their expected characteristics, and optimization of basic parameters.

We performed the numerical drift-diffusion simulation of the lasers with functionally integrated optically modulators and UTC photodetectors. According to the obtained results, it is reasonable to continue the research of the lasers-modulators and optimize their parameters, especially the lifetime of photons in the resonant cavity. Besides that, the application of the lasers-modulators as parts of on-chip optical interconnections in ICs requires the development of the advanced methods aimed at the reduction of response time of high-speed photosensitive devices.