Simudo: a device model for intermediate band materials

TL;DR

Simudo is a novel free device modeling tool using finite element methods to simulate intermediate band materials, offering higher convergence rates and extensibility for semiconductor device analysis.

Contribution

It introduces the first freely available FEM-based device model capable of simulating intermediate band materials with robust convergence.

Findings

Simudo converges faster and more accurately than Synopsys Sentaurus on benchmark tests.

It effectively models complex semiconductor phenomena with high numerical stability.

The tool is designed for extensibility to include various physical processes.

Abstract

We describe Simudo, a free Poisson/drift-diffusion steady state device model for semiconductor and intermediate band materials, including self-consistent optical absorption and generation. Simudo is the first freely available device model that can treat intermediate band materials. Simudo uses the finite element method (FEM) to solve the coupled nonlinear partial differential equations in two dimensions, which is different from the standard choice of the finite volume method in essentially all commercial semiconductor device models. We present the continuous equations that Simudo solves, show the FEM formulations we have developed, and demonstrate how they allow robust convergence with double-precision floating point arithmetic. With a benchmark semiconductor pn-junction device, we show that Simudo has a higher rate of convergence than Synopsys Sentaurus, converging to high accuracy with a considerably smaller mesh. Simudo includes many semiconductor phenomena and parameters and is designed for extensibility by the user to include many physical processes.

Tables

Table 1. Table 1: Comparison of selected device models

	IB optics with photofilling	IB nonradiative processes	IB transport	Junctions	2D
Sentaurus	N	Y	limited	Y	Y
PC1D Clugston97	N	Y	N	Y	N
SCAPS Burgelman00	Y	Y	N	Y	N
TiberCAD Maur08	N	Y	Y	Y	Y
Martí Marti02	N	N	Y	N	N
Strandberg Strandberg11	Y	N	Y	N	N
Tobias Tobias11 , Yoshida Yoshida12a	Y	N	Y	Y	N
Simudo	Y	Y	Y	Y	Y

Table 2. Table 2: Common symbols used in this manuscript.

Symbol	Definition
$u_k$	Carrier density in band $k$
${ℰ}_k$	Band edge energy of band $k$; central energy of IB $k$
$s_k$	Sign of carriers in band $k$ (+ for VB, - for CB)
$w_k$	Quasi-Fermi level of carriers in band $k$
${\overrightarrow{j}}_k$	Current density in band $k$
${\mu }_k$	Mobility of carriers in band $k$
${𝒟}_k$	Diffusion coefficient of carriers in band $k$
$S_k$	Surface recombination velocity of carriers in band $k$
$N_k$	Effective density of states for nondegenerate band $k$
	Integrated density of states for intermediate band $k$
$f_k$	Filling fraction $f_k=u_k/N_k$ of IB $k$
$f_{k,\pm }$	$\left(1-f_k\right)$ for positive sign and $f_k$ for negative sign
$f_{k,0}$	Charge neutral filling fraction of IB $k$
$g_k$	Net generation in band $k$ due to all generation and recombination processes
$T$	Temperature
$k_B$	Boltzmann constant
$q$	Elementary charge
$\varphi ,\overrightarrow{E},\rho$	Electrostatic potential, electric field, and charge density
${\alpha }_{fi,\lambda }$	Optical absorption coefficient from band $i$ to $f$ at vacuum wavelength $\lambda$
${\sigma }_{fi}^{\text{opt}}$	Optical cross section from band $i$ to $f$
${\Phi }_{\lambda ,\hat{s}}$	Photon spectral flux density at wavelength lambda in direction $\hat{s}$
${\Phi }_{[{\lambda }_1,{\lambda }_2],\hat{s}}$	Photon flux in direction $\hat{s}$ from ${\lambda }_1$ to ${\lambda }_2$, i.e., ${\int }_{{\lambda }_1}^{{\lambda }_2}{\Phi }_{\lambda ,\hat{s}}d\lambda$
$\hat{n}$	Surface normal vector
$\hat{s}$	Direction of light propagation

Table 3. Table 3: Parameters modeled on the device from Marti2002a , for Figs.  8 – 10 .

Value	Definition
${ℰ}_C=1.67\mathrm{eV}$	Conduction band edge energy
${ℰ}_I=1.10\mathrm{eV}$	Intermediate band energy
${ℰ}_V=0\mathrm{eV}$	Valence band edge energy
$N_C=N_V=5\times {10}^{18}{\mathrm{cm}}^{-3}$	CB and VB effective density of states
$N_I={10}^{17}{\mathrm{cm}}^{-3}$	IB density of states
${\mu }_C={\mu }_V=2000{\mathrm{cm}}^2/V/s$	CB and VB mobility
${\mu }_I=0.001-300{\mathrm{cm}}^2/V/s$	IB mobility
${\alpha }_{CV}={10}^4{\mathrm{cm}}^{-1}$	Absorption coefficient for CV process
${\sigma }_{CI}^{\text{opt}}=\left(2-10\right)\times {10}^{-13}{\mathrm{cm}}^2$	Optical cross section for CI process
${\sigma }_{IV}^{\text{opt}}=2\times {10}^{-13}{\mathrm{cm}}^2$	Optical cross section for IV process
$\epsilon =13{\epsilon }_0$	Dielectric constant
$T_s=6000K$	Sun temperature
$T_c=300K$	Cell temperature
$X=1000$	Solar concentration factor
$w_{\text{IB}}=1.3\mu m$	IB region length
$f_{I,0}=1/2$	Charge-neutral IB filling fraction

Equations

$${\vec{j}}_k$$

$$\frac{\partial u_k}{\partial t}$$

$$u_k$$

$$\vec{\nabla }{u}_k$$

$$=-s_k\frac{u_k}{k_{B}T}\vec{\nabla }(w_k+q\varphi ),$$

$${\vec{j}}_k$$

$$u_I$$

$$\vec{\nabla }{\Phi }_{\lambda ,\hat{s}}\cdot \hat{s}$$

$${\alpha }_{\lambda }=\sum _{i,f}{\alpha }_{fi,\lambda },$$

$$g_k^{\text{opt}}$$

$${\alpha }_{CI,\lambda }$$

$${\alpha }_{IV,\lambda }$$

$${\alpha }_{kI,\lambda }={\sigma }_{kI,\lambda }^{\text{opt}}{u}_{I,s_k},$$

$$r^{\text{SRH}}$$

$$r_{Ik}^{\text{SR}}$$

$$r_{Ik}^{\text{rad}}=\left[e^{-s_k(w_k-w_I)/k_{B}T}-1\right]u_{I,s_k}{\mathcal{I}}_{Ik},$$

$${\mathcal{I}}_{Ik}$$

$$u_{I,s_k}=\frac{u_{I,-s_k}{u}_k}{u_{k,1}}{e}^{s_k(w_k-w_I)/k_{B}T},$$

$$r_{Ik}^{\text{rad}}=\left[1-e^{s_k(w_k-w_I)/k_{B}T}\right]\frac{u_{I,-s_k}{u}_k}{u_{k,1}}{\mathcal{I}}_{Ik}.$$

$${\vec{j}}_k\cdot \hat{n}{|}_{\Gamma }=S_k(u_k-u_{k0}){|}_{\Gamma },$$

$$\hat{s}\cdot \vec{\nabla }(\hat{s}\cdot \vec{\nabla }{\Phi }_{\lambda ,\hat{s}})+\hat{s}\cdot \vec{\nabla }({\alpha }_{\lambda }{\Phi }_{\lambda ,\hat{s}})$$

$${\Phi }_{\lambda ,\hat{s}}$$

$$\hat{s}\cdot \vec{\nabla }{\Phi }_{\lambda ,\hat{s}}+{\alpha }_{\lambda }{\Phi }_{\lambda ,\hat{s}}$$

$${\Phi }_{[{\lambda }_1,{\lambda }_2],\hat{s}}={\int }_{{\lambda }_1}^{{\lambda }_2}{\Phi }_{\lambda ,\hat{s}}\mathrm{d}\lambda ,$$

$$\hat{s}\cdot \vec{\nabla }(\hat{s}\cdot \vec{\nabla }{\Phi }_{[{\lambda }_1,{\lambda }_2],\hat{s}})+\hat{s}\cdot \vec{\nabla }({\alpha }_{\lambda }{\Phi }_{[{\lambda }_1,{\lambda }_2],\hat{s}})=0.$$

$$\vec{\nabla }\cdot (\epsilon \vec{\nabla }\varphi )$$

$$\vec{\nabla }\varphi$$

$$\vec{\nabla }\cdot (\epsilon \vec{E})$$

$$\rho$$

$$\underset{\text{natural BC}}{\underbrace{{\oint }_{\Gamma }\vec{\psi }\cdot \hat{n}{\varphi }_{\text{BC}}}}-{\int }_{\Omega }(\nabla \cdot \vec{\psi })\varphi +{\int }_{\Omega }\vec{\psi }\cdot \vec{E}$$

$$\displaystyle\int_{\Omega}v\,\big{(}\nabla\cdot(\varepsilon\vec{E})\big{)}-\int_{\Omega}v\,\rho$$

$$0$$

Full text

1 Introduction

Device models are essential components of the development of semiconductor devices, from transistors to solar cells to lasers. Standard semiconductor device models, such as Synopsys Sentaurus, treat materials with 0, 1, or 2 bands (i.e., dielectrics, metals, and semiconductors, respectively) along with an electrostatic potential. At a given location in a given material, each band has its own carrier concentration, with particle motion given by diffusion and electric-field-induced drift. Since the electric field itself depends on particle motion, the resulting Poisson/drift-diffusion (PDD) equations are nonlinear and require numerical solution in the general case Bank1983 ; Fichtner1983a ; Markowich1986 ; Piprek18 ; Schenk1998 .

A new class of materials, called intermediate band (IB) materials, has been developed over the last 20 years with the goal of improving solar cell efficiency and producing effective infrared photodetectors Luque97 ; Okada15 ; Mailoa14 ; Berencen17 ; Wang18 . These IB materials are like semiconductors except they have an extra band of allowed electronic energy levels above the valence band (VB) and below the conduction band (CB), as shown in Fig. 1. Such a band structure permits optical absorption from VB to IB and from IB to CB, which is the key to the increased solar cell efficiency Luque97 . It is also possible to consider multiple IBs, though such materials have not yet been realized in practice Brown02 .

Where IB devices have been made, they have not generally been highly efficient, which is believed to be largely due to fast nonradiative recombination processes Okada15 ; Marti06 ; Wang09 ; Lopez11 ; Mailoa14 ; Berencen17 ; Wang18 ; Sullivan15 . It has not been possible, however, to perform standard device modeling to optimize these devices, to determine the ideal layer thicknesses, doping levels, etc., since standard semiconductor device models do not allow the possibility of treating a third band. Therefore, we do not know what efficiencies existing IB materials could permit, if they were optimized. Interpreting experiments on IB materials and designing the best devices require device modeling capabilities.

In order to describe the basic physics of IB devices, one must be able to describe

Optical processes between CB, VB, and IB, with rates dependent on IB filling fraction $f$, 2. 2.

Nonradiative processes between CB, VB, and IB, with rates dependent on $f$, 3. 3.

Carrier transport within the IB, 4. 4.

Junctions with standard semiconductors.

There is a large array of standard numerical semiconductor device models based on the coupled Poisson and carrier-continuity equations, including general purpose ones, such as Synopsys Sentaurus and Silvaco, as well as more specialized models such as Crosslight, which includes modelling of quantum well physics and a coupled treatment of carrier-density dependent optics for lasers, and TiberCAD Maur08 . Nextnano++ also includes features specific to quantum structures and can solve an 8-band k.p model self-consistently with the Poisson and drift-diffusion calculations Birner2007 . There are also more focused ones, such as PC1D Clugston97 ; Haug16 , AFORS-HETVarache15 , SCAPS Burgelman00 , and Solcore Alonso-Alvarez2018 , which are 1D models focused on solar cells. Many of these models allow treatment of deep-lying states inside the semiconductor band gap, primarily as Shockley-Read-Hall (SRH) trapping and recombination centers Shockley52 . Sentaurus and PC1D, for example, do not permit optical generation from the deep-lying states. SCAPS does permit both thermal and optical processes, but does not consider transport of carriers inside the defect band. We do not attempt a full characterization of all the available device models, but Table 1 shows which of these requirements are met by these device models.

There have been a number of device models developed specifically for IB materials, mostly for solar cells, all in steady state. These include traditional Luque97 ; Cuadra04a ; Levy08 ; Hu10 and Boltzmann-approximation Strandberg17 ; Strandberg17a detailed balance models, semianalytic models in the drift Lin09 and diffusion Marti02 ; Navruz08 ; Krich14 limits, and PDD models Yoshida10 ; Strandberg11 ; Tobias11 ; Yoshida12a . The semianalytic models are specific to either the drift or diffusive limits, while the PDD models allow treatment of IB regions that are neither fully depleted nor fully quasi-neutral. A comparison of the features of the PDD models is also included in Table 1. To our knowledge, none has been released as open-source software.

Here we introduce Simudo, a free and open source steady state PDD solver with self-consistent optics for arbitrary numbers of bands. Simudo uses the finite element method (FEM) to solve the coupled Poisson, drift-diffusion, and Beer-Lambert optical propagation equations self-consistently, when necessary including changing $f$ according to local generation and recombination, with associated changes in the optical absorption coefficient. Simudo has built-in radiative recombination, Shockley-Read trapping, and SRH recombination models in the non-degenerate limit and is straightforward to extend to include other models of generation or recombination. All of the band parameters, from energies to mobilities to cross sections, can vary in space or as functions of other parameters.

Simudo has a number of innovations in its formulation of the problem, described below, and allows high-accuracy simulation of benchmark semiconductor problems while working with 64-bit arithmetic, making it useful both for standard semiconductor simulations as well as for IB devices. It is written in the Python programming language, using the FEniCS platform to solve the FEM problem AlnaesBlechta2015a . It exposes an easy-to-use API for defining problems and extracting results. It is designed for two-dimensional systems and is available for download at https://github.com/simudo/simudo.

Semiconductor device modeling is typically performed using the finite volume method (FVM), which ensures local charge conservation at each cell in the domain Eymard00 . Finite element methods are less common in semiconductor device models, though there are a number of examples He91 ; Maur08 ; Nachaoui99 ; Bochev15 . FEM is widely used for related advection-diffusion problems in computational fluid dynamics (CFD) studies Cockburn00 . Commercial packages for CFD, as for semiconductor device modeling, are generally based on the FVM method. TiberCAD, a commercial device modeling package with many novel features, uses FEM with continuous basis functionsMaur08 . Finite element methods including discontinuous local basis functions, called discontinuous Galerkin (DG) methods, also permit local charge conservation Cockburn00 , and they have recently begun to be applied to semiconductor device problems Kumar16 ; Kumar17 . FEM methods simplify consideration of complicated simulation domains and in theory allow higher-order convergence of solutions, but performance of such methods can only be determined with testing. We use such a DG-FEM method here to produce a general purpose steady state PDD solver capable of treating IB systems, and we show that Simudo realizes the higher-order convergence with mesh size, converging much more rapidly than Synopsys Sentaurus as the mesh spacing is reduced. As shown in Sec. 3.5, for a reference pn-diode, Simudo demonstrates quartic self-convergence with mesh density while FVM-based Synopsys Sentaurus demonstrates only quadratic convergence. In the reference problem, Simudo achieves 5-6 digits of convergence with 193 mesh points while Sentaurus requires more than 3000. Simudo provides both a flexible framework for the study of IB devices and also a freely available example of a DG-FEM semiconductor device model.

In Section 2, we define the coupled partial differential equations (PDEs) Simudo solves. Section 3 describes the heart of Simudo, giving in detail the conversion of the equations of Section 2 to the weak forms solved using FEM. This section describes the choices for dynamical variables, the weak forms used for FEM, and how these choices enable Simudo to achieve accuracy despite the problems of finite precision arithmetic. This section concludes with a comparison to Synopsys Sentaurus on a benchmark pn-diode, showing the high quality of Simudo’s results. Section 4.1 gives examples of setting up a simple problem using the API, including examples of its convenient topology definitions and Section 4.2 demonstrates the extensibility of Simudo to include new physical processes (in this case, Auger recombination). Section 4.3 demonstrates the use of Simudo to analyze a system originally studied in Ref. Marti02 , showing that its model works better than had been anticipated in the case with equal subgap optical absorption cross sections, but that unequal subgap absorption cross sections produce more complicated phenomena that require IB transport to describe properly.

2 Statement of problem

In this section, we describe the mathematical model of the steady state PDD and optical problems we use in Simudo. Carriers both drift in response to electric fields and diffuse. Carriers are generated optically and recombine using a variety of models. The local carrier concentration determines both the electric field and the optical absorption coefficients, so the transport, Poisson, and optical propagation equations are all coupled. Symbols used in this manuscript are summarized in Table 2.

2.1 Carrier transport and generation

We consider a CB, a VB, and some number of IBs under the assumption that the carrier population in each band is in local quasi-equilibrium with a temperature $T$ and quasi-Fermi level $w_{k}$, where $k$ can be one of $\{C,V,I\}$ for the CB, VB, and IB, respectively. In the case of multiple IBs, $k$ can take values $I_{1},I_{2},\dots$, indexing the various IBs, but we simplify the following discussion to consider the case of just one IB, indexed as $I$.

In the most common approximation of semiconductor device modeling, the carrier dynamics in each band can be described by the drift-diffusion equation and the continuity equation. Letting $u_{k}$ represent the carrier concentration in band $k$, $u_{V}$ and $u_{C}$ are the hole and electron concentrations, respectively, which we use interchangeably with their standard symbols, $p$ and $n$. We let $s_{k}=\pm 1$ give the charge of the carriers in band $k$, $+1$ for the VB and -1 for the CB. Then

[TABLE]

where $\vec{j}_{k}$ is the current density of carriers in band $k$, $\mu_{k}$ is the carrier mobility, $\mathcal{D}_{k}$ is the carrier diffusion constant, $\vec{E}$ is the electric field, $q$ is the elementary charge, and $g_{k}$ contains all the generation, trapping, and recombination processes (see Section 2.2). For non-degenerate bands in which $w_{k}$ is sufficiently far from the band edge $\mathcal{E}_{k}-q\phi$, we can write

[TABLE]

where $N_{k}$ is the effective density of states of band $k$, $\phi$ is the electrostatic potential, and $k_{B}$ is Boltzmann’s constant.

Then, assuming $\mathcal{E}_{k}$ is spatially constant,

[TABLE]

For such nondegenerate bands, the Einstein relation gives $\mu_{k}=q\mathcal{D}_{k}/k_{B}T$, from which Eq. 1a gives Fichtner1983a

[TABLE]

which we use instead of Eq. 1a. Equation 5 also applies to the case of degenerate bands, as shown in Poupaud1991a , even though the Einstein relation requires a modification. Moreover, Eq. 5 applies in the case of spatially-varying band structure (e.g., spatially-varying $N_{c},\mathcal{E}_{c}$)Marshak1984a , so it is considerably more general than this derivation.

Since an intermediate band is often partially filled, we cannot model it using the non-degenerate approximation of Eq. 2. We write $D_{I}(E)$ for the density of states of the IB, such that $N_{I}=\int\mathrm{d}E\>D_{I}(E)$ is the total density of IB states. If the IB has quasi-Fermi level $w_{I}$, the electron concentration is

[TABLE]

If the bandwidth of the IB is narrow relative to $k_{B}T$, we can approximate the IB density of states as a Dirac delta $D_{I}(E)=N_{I}\delta(E-\mathcal{E}_{I})$, and so

[TABLE]

where $f_{I}$ is the filling fraction of the IB, and can be written as $f_{I}=f(\mathcal{E}_{I}-w_{I}-q\phi)$ where $f(E)$ is the Fermi function. We work in this limit for the remainder of this manuscript. Extending beyond this sharp-IB case is not difficult but requires more cumbersome notation.

2.2 Carrier generation and recombination

Each band’s continuity equation (Eq. 1b) has a generation term $g_{k}$. This term is the sum of contributions from all generation and recombination processes to the band, which depend on which physical models are included in the simulations. We now specify the details of optical generation $g_{k}^{\text{opt}}$ and a variety of recombination processes $r_{k}$, each of which enters either as a negative or positive contribution to $g_{k}$, as required for the process.

2.2.1 Optical carrier generation

Modeling optical carrier generation requires modeling the changing light intensity through the device. We use a simple Beer-Lambert model for optical propagation and absorption

[TABLE]

where $\Phi_{\lambda,\hat{s}}$ is the photon spectral flux at vacuum wavelength $\lambda$ and direction of propagation $\hat{s}$ and $\alpha_{\lambda}$ is the total absorption coefficient, which can be written as

[TABLE]

where $\alpha_{fi,\lambda}$ is the absorption coefficient for the optical process at wavelength $\lambda$ that moves a carrier from band $i$ to band $f$. In the usual semiconductor case, $\alpha_{VC,\lambda}=0$ and $\alpha_{CV,\lambda}$ is finite for $\lambda$ corresponding to energies larger than the band gap. Free-carrier absorption is included in $\alpha_{ii,\lambda}.$ The carrier generation rate in band $k$ due to optical processes is then

[TABLE]

Further details of the optical propagation model are described in Section 2.3.

In nondegenerate bands, there are always enough carriers to excite in or out of a band. That is, the valence band always has electrons available, and the conduction band has empty states available to be filled, so the absorption coefficient $\alpha_{CV,\lambda}$ is insensitive to the free carrier density in the bands. In an IB, however, the VB$\rightarrow$IB process requires empty states in the IB while the IB$\rightarrow$CB process requires filled states in the IB. To capture this phenomenon, we write

[TABLE]

where $\sigma_{fi,\lambda}^{\text{opt}}$ is the optical capture cross section from band $i$ to $f$ at wavelength $\lambda$ of a single intermediate state. We can combine these equations into a single expression,

[TABLE]

where $u_{I,-}$ is just $u_{I}$ and $u_{I,+}=N_{I}-u_{I}$ is the number of holes in band $I$, and $\alpha_{kI}$ is understood to be $\alpha_{IV}$ when $k=V$. Since $\alpha_{fi,\lambda}$ depends on the carrier concentrations, and the carrier concentrations depend on $\alpha_{fi,\lambda}$ (through the generation rate $g_{k}^{\textnormal{opt}}$), the transport and the optical models feed into each other, so they must be solved in a self-consistent manner.

2.2.2 Recombination and trapping

Simudo offers several built-in radiative and nonradiative recombination and trapping mechanisms using the non-degenerate limit for the CB and VB, each including an equivalent thermal generation. An example is the SRH recombination model with a single trap level at energy $\mathcal{E}_{I}$ Shockley52 , in which two trapping processes (of an electron and a hole) produce a recombination event, with recombination rate

[TABLE]

where $\tau_{p},\tau_{n}$ are the carrier lifetimes and $p_{1},n_{1}$ are the carrier concentrations of holes and electrons, respectively, if their quasi-Fermi levels were equal to $\mathcal{E}_{I}$. This $r^{\text{SRH}}$ appears as a negative contribution to $g_{k}$ for both CB and VB.

We can model traps as intermediate bands with $w_{I}$ tracked explicitly, in which case we implement standard Shockley-Read trapping Shockley52 ,

[TABLE]

where $f_{I,-s_{k}}$ is the IB filling fraction of carriers with charge $-s_{k}$, and $\tau_{k}$ is the Shockley-Read lifetime for band $k$, as in Eq. 13 Shockley52 . Note that $r_{IC}^{\text{SR}}$ makes a negative contribution to $g_{C}$ and a positive contribution to $g_{I}$, while $r_{IV}^{\text{SR }}$ makes a negative contribution to both $g_{V}$ and $g_{I}$.

Simudo also implements radiative trapping from band $k=C,V$ to $I$. When we use Boltzmann statistics rather than Bose statistics for the emitted photons, which is valid when $|w_{k}-w_{I}|$ remains at least a few $k_{B}T$ below $|\mathcal{E}_{k}-\mathcal{E}_{I}|$, as in Ref. Strandberg11 , then the radiative trapping can be written

[TABLE]

where

[TABLE]

where $n_{r}$ is the index of refraction, and $u_{k,1}$ is either $n_{1}$ or $p_{1}$ for $k=C,V$, respectively. Note that Ref. Strandberg11 includes only the recombination term, and we add the corresponding thermal generation term, which is the -1 in Eq. 15. We can re-express Eq. 15 in a similar form to the nonradiative terms by using the relation

[TABLE]

which follows from Eqs. 2, 7. Then

[TABLE]

As with Eq. 2, Eqs. 14-16 are valid in the non-degenerate limit where $w_{k}$ does not approach $\mathcal{E}_{k}$ but full degenerate statistics are used for the IB. Extensions to the degenerate limit can be added, if desired. Simudo also treats standard radiative recombination between conduction and valence bands Nelson03 .

We also treat surface recombination at external surfaces $\Gamma$ of the device, which imposes a boundary condition

[TABLE]

where $S_{k}$ is the surface recombination velocity of carriers in band $k$ at boundary $\Gamma,$ $\hat{n}$ is the normal to $\Gamma$, and $u_{k0}$ is the carrier concentration at equilibrium McIntosh14 . The current release of Simudo supports only $S_{k}=0$ or $\infty$, which impose $\vec{j}_{k}\cdot\hat{n}|_{\Gamma}=0$ or $(u_{k}-u_{k0})|_{\Gamma}=0$, respectively.

2.3 Optical equations

For each wavelength, we need to solve the optical propagation according to Eq. 8. For stability of the numerical solution, it is convenient to use a second-order equation so that we can apply boundary conditions on both the inlet and outlet boundaries Zhao2013a . We take the derivative of Eq. 8 with respect to the direction of propagation,

[TABLE]

With no reflection from the back, the boundary conditions are then

[TABLE]

where $\Phi_{\lambda,\hat{s}}^{0}$ is a spectral photon flux at the inlet boundary.

In the case where $\alpha_{\lambda}$ is constant for $\lambda$ in an interval $[\lambda_{1},\lambda_{2}]$, the optical flux at all wavelengths in that range obeys Eq. 19 and can thus be treated together. We can write

[TABLE]

where $\Phi_{[\lambda_{1},\lambda_{2}],\hat{s}}$ is a photon flux (where $\Phi_{\lambda,\hat{s}}$ is a spectral photon flux). In this case, we have

[TABLE]

Simudo uses this form, which allows simple treatment of piecewise constant absorption coefficients with a small number of optical fields $\Phi_{[\lambda_{1},\lambda_{2}],\hat{s}}$. When optical fields with only one propagation direction $\hat{s}$ are considered, we write the spectral flux density $\Phi_{\lambda}$ and the flux density $\Phi_{[\lambda_{1},\lambda_{2}]}$.

2.4 Poisson’s equation

In electrostatics, Poisson’s equation relates $\phi$, the charge density $\rho$, and the permittivity $\varepsilon$,

[TABLE]

It can also be split into two equations

[TABLE]

where $\vec{E}$ is the electric field.

The charge density $\rho$ is the sum of the static charge and the mobile charge in each band. In an IB material,

[TABLE]

where $N_{A}$, $N_{D}$ are the shallow acceptor and donor doping concentrations, respectively, the mobile charge in the IB is $qN_{I}(f_{I}-f_{I,0})$, with $f_{I,0}$ the IB filling fraction of the bulk IB material at $T=0\mathrm{\>K}$. For a donor-type IB $f_{I,0}=1$, and for an acceptor-type IB $f_{I,0}=0$. Note that in writing the shallow dopant terms $N_{D}$ and $N_{A}$, we are assuming complete ionization of these impurities.

3 Numerical method

Simudo uses the finite element method (FEM) to solve the coupled Poisson/drift-diffusion and optics problems, given by Eqs. 1b,5,19, and 22. The FEM method divides the simulation domain $\Omega$ into cells, which are generally triangles in 2D, and enforces a weak form of the desired PDE’s with a set of test functions defined on those cells, with boundary conditions applied on the domain boundary $\Gamma$. The method is well-described in many reference texts Johnson1987a ; Gockenbach2006a ; Donea2003a . In this section, we detail the weak forms used for these coupled equations and the solution method for the resulting nonlinear system. We benchmark Simudo against the industry standard Synopsys Sentaurus commercial simulator on a standard semiconductor problem to show the quality of our results.

3.1 Solution method

The PDD problem is a coupled nonlinear system of PDEs, which we solve iteratively using Newton’s method as implemented in the FEniCS package. The solution procedure is outlined in Fig. 2. The goal is to find a solution $y=(\phi,\vec{E},w_{C},w_{I},w_{V},\vec{j}_{C},\vec{j}_{I},\vec{j}_{V})$ that satisfies Eqs. 1b,5,22 and associated $\Phi_{\lambda}$ that obeys Eq. 19. The optical problem is solved alongside the PDD problem in a self-consistent manner. That is, the PDD subproblem produces the absorption coefficient $\alpha(\lambda)$ (which, for processes involving the IB, depends on the filling fraction). The optical subproblem is then solved using these absorption coefficients, yielding a new photon flux $\Phi_{\lambda}$, which is fed back into the PDD where it enters in the optical carrier generation process, and the cycle iterates until a self-consistent solution is found.

The convergence of Newton’s method depends on the quality of the initial guess. Steps 1 and 2 in Fig. 2 are the pre-solver, which is used once to make the initial guess for the main Newton solver, illustrated in step 3. The full procedure is:

For each point in space, calculate the equilibrium Fermi level $E_{F}$ at that point assuming local charge neutrality and $\phi=0$.

Physically, this step finds the $E_{F}$ in each location that makes it charge neutral, before any charge is allowed to flow. 2. 2.

Determine the built-in potential $\phi_{\mathrm{bi}}$ of the equilibrium system, using $\phi_{\text{0}}=E_{F}/q$ as the initial guess. That is, solve only Eq. 22 for $\phi$ while keeping all $w_{k}=0$ and thus all $\vec{j}_{k}=0$.

Physically, this step allows charge to move, forming depletion regions as the carriers move to achieve a zero-current configuration that satisfies Poisson’s equation. The carrier density inside bulk-like regions of space changes little from the bulk equilibrium value in step 1, making $\phi_{\textnormal{0}}$ an excellent guess in large regions of space. 3. 3.

Main solver loop Adaptively ramp up light intensity and/or bias, starting with thermal equilibrium (dark, no bias). Each solution requires a loop of Newton iterations. Within each Newton iteration, do the following:

1. (a)

**Optical calculation **For each optical field $\lambda$, solve for the photon flux $\Phi_{\lambda}$ given the latest value of $\alpha_{\lambda}$. Note that Eq. 19 is linear when $\alpha_{\lambda}$ is fixed. 2. (b)

**PDD Newton step **Perform one Newton step of the PDD problem.

1. i.

Solve for $\delta y=\delta(\phi,\vec{E},w_{C},w_{I},w_{V},\vec{j}_{C},\vec{j}_{I},\vec{j}_{V})$. Use the value of $\Phi_{\lambda}$ (and thus optical carrier generation) computed in the previous step. 2. ii.

Update $y\leftarrow y+\delta y$.

As an option, logarithmic damping can be applied to $\delta y$ to prevent Newton’s method from diverging, e.g., $y\leftarrow y+\textnormal{LogDamping}(\delta y)$ where $\textnormal{LogDamping}(z)=\mathop{\mathrm{sgn}}\nolimits(z)\log(1+c|z|)/c$ for $c=1.72$ or other user-defined value Gaury2019a .

We now describe the weak forms that we use for each of Eqs. 1b,5, 19, and 21. There is much flexibility in the choice of particular weak forms, all of which can be equivalent to the same strong form. In Section 3.3.2 we illustrate the use of a partitioned offset representation for $w_{k}$, which allows internal currents to be calculated accurately with double-precision arithmetic.

3.2 Poisson equation

Here we introduce the formulation we use to implement Eq. 22. We use a mixed method to solve for both $\phi$ and $\vec{E}$ explicitly Brezzi1985a ; Roberts1991a . The potential $\phi$ is represented as a superposition of discontinuous Galerkin (DG) basis functions of order $d_{\textnormal{poisson}}-1$ (cell-wise discontinuous polynomials), and $\vec{E}$ is represented using Brezzi-Douglas-Marini basis functions of order $d_{\textnormal{poisson}}$ (cellwise discontinuous polynomials with continuous normal component on cell boundaries) Brezzi1985a . We use the BDM space for all vector quantities, including $\vec{E}$ and $\vec{j}_{k}$. The BDM space is H(div) conforming, meaning the divergence is accurately calculated and fluxes between cells are preserved, which makes it a natural choice for conserved or almost conserved vector quantities.111In the BDM space, the normal fluxes are shared by adjacent elements. The flux exiting the perimeter of a collection of cells exactly equals the sum of fluxes out of each of the cells, with exact arithmetic. While $\vec{j}_{k}$ is not a conserved quantity, due to generation and recombination that occur inside of cells, the BDM space ensures that $\vec{j}_{k}$ is accurately preserved when passing between cells. A method using CG or DG functions for $\vec{j}$ would be susceptible to numerical errors associated with non-conservation of currents between cells, and we show in Sec. 3.3.2 that Simudo conserves current well in a pn diode. In the results below, $d_{\textnormal{poisson}}=2$.

We multiply Eq. 22a by test function $\vec{\psi}\in\mathrm{BDM}(d_{\textnormal{poisson}})$ and Eq. 22b by test function $v\in\mathrm{DG}(d_{\textnormal{poisson}}-1)$, then integrate each spatially, giving the weak forms

[TABLE]

which must hold for every test function $\vec{\psi}$ and $v$, where $\Omega$ is the full domain and $\Gamma$ is the boundary. Note that Eq. 24 includes an integration by parts. In this case, the electric field BC is an essential BC, imposed by reducing the set of test functions to those that satisfy the BC, while the potential BC is a natural BC.

3.3 Transport equations

The drift-diffusion equations are often numerically challenging to solve in semiconductors. In carrier density-based formulations, poor resolution of the gradients of $u_{k}$ makes linear interpolation of current density unstable, which the Scharfetter-Gummel box method corrects for FVM methods Bank1983 . Additionally, catastrophic cancellation can occur in Eq. 1a, e.g., for the majority carrier in a quasi-neutral region of a semiconductor, when the drift and diffusion contributions are nearly equal in magnitude. The current is given by the difference and can be hard to evaluate with finite precision arithmetic. We address these issues by using a quasi-Fermi-level-based representation for carrier density Cummings09 ; Nachaoui99 . Calculating $\vec{\nabla}w_{k}$ in finite precision for Eq. 5 can also be challenging when $w_{k}$ is very flat, and in Sec. 3.3.2 we introduce a partitioned offset representation for $w_{k}$ to allow accurate determination of $\vec{\nabla}w_{k}$ with essentially no extra computational cost. We use a mixed FEM method that solves explicitly for both $w_{k}$ and the current density $\vec{j}_{k}$. As described in Sec. 3.2, the BDM space of basis functions enforces local current conservation in the solutions, which also enables local current densities to be well determined. Without the mixed method, local current conservation is enforced indirectly, and we were not able to obtain well-converged results for local currents.

3.3.1 Quasi-Fermi level formulation

The quasi-Fermi level $w_{k}$ is represented as a superposition of DG basis functions of order $d_{\textnormal{transport}}-1$, and the current density $\vec{j}_{k}$ is represented using BDM basis functions of order $d_{\textnormal{transport}}$. Section 3.2 contains a discussion of these functions’ properties and of the mixed method. In the results below, $d_{\textnormal{transport}}=2$.

We derive weak forms of Eq. 1b and Eq. 5, multiplying Eq. 1b by test function $v\in\mathrm{DG}(d_{\textnormal{transport}}-1)$, taking the dot product of Eq. 5 with the test function $\vec{\psi}\in\mathrm{BDM}(d_{\textnormal{transport}})$, and integrating each equation spatially, giving

[TABLE]

where the second equation was obtained by a further integration by parts.

3.3.2 Quasi-Fermi level offset partitioning

As written, Eq. 26b still suffers from a form of catastrophic cancellation in its last term, which corresponds to the gradient term in Eq. 5. Since $\int_{\Omega}\vec{\nabla}\cdot\vec{\psi}=0$ for $\vec{\psi}\in\text{BDM}(d_{\text{transport}})$, the last term is nonzero only if $w_{k}$ varies within the domain where $\psi$ is nonzero. $w_{k}$ can be extremely flat, for example in quasineutral regions, which makes this integral hard to calculate with finite arithmetic precision. This difficulty is more apparent in Eq. 5: if $\vec{j}_{k}/\mu_{k}u_{k}=\vec{\nabla}w_{k}$ is small,222relative to $|w_{k}|/(\textnormal{mesh size})$ a representation of $w_{k}$ that stores its value on the nodes of the mesh cannot resolve such small changes in $w_{k}$ across space.

We circumvent this issue by using an offset representation for $w_{k}$. The idea is to give each cell in the domain its own (spatially constant) base quasi-Fermi level $w_{k0}$ relative to which the new dynamical variable $\delta w_{k}$ is expressed. That is, $w_{k}=w_{k0}+\delta w_{k}$ where $\delta w_{k}$ is the quantity we actually solve for instead of $w_{k}$. Before every Newton iteration step, the $w_{k0}$ of each cell is initialized to the cell average of $w_{k}$ from the previous iteration. This representation allows small spatial changes of $\delta w_{k}$ to be accurately represented, enabling accurate determination of the current.

The last remaining question is how to adjoin regions with different base $w_{k0}$ values. We connect them by adding a surface integral jump term to Eq. 26b, resulting in

[TABLE]

where $[w_{k0}]$ is the jump operator, which takes the difference between the values of a discontinuous expression on either side of a facet. The rest of this section is dedicated to deriving that term and comparing the result to a formulation without the offset representation.

We substitute $w_{k}=w_{k0}+\delta w_{k}$ into Eq. 26b, and we obtain

[TABLE]

Our goal now is to rewrite the $(\star)$ term. Since $w_{k0}$ is constant on each cell $K$, $\vec{\nabla}w_{k0}=\vec{0}$ within each cell. We integrate by parts using

[TABLE]

yielding

[TABLE]

Summing over all cells $K$,

[TABLE]

which, plugged into Eq. 28, yields Eq. 27.

We perform a test of Simudo, which uses the partitioned offset representation, against the identical model without the offset representation. We consider a standard silicon pn-junction diode with symmetric doping of $10^{18}$ cm*-3* and SRH lifetimes of 1 ns and 1 $\mu$s in the p- and n-type regions, respectively. Each region has a length of 250 nm, for total device length of 500 nm. Although the problem is one-dimensional, we consider a 2D region with a height of $1\mathrm{\>\mu m}$. At each contact, the majority carrier has an infinite surface recombination velocity while the minority carrier has zero surface recombination. We use a mesh with 769 points in the x-direction, which is tightest near the contacts and junction and expands out geometrically toward the middle of the quasi-neutral regions. The mesh has 2 points in y-direction, and further details of the mesh are given in Sec. 3.5.

The offset representation allows internal current densities to be resolved accurately throughout the device. Figure 3(a) shows the electron and hole currents under 0.16 V bias. Without the offset representation, the majority currents are poorly resolved, due to the inability to resolve $\vec{\nabla}w_{k}$ with double-precision arithmetic. The majority currents in the no-offset model become worse as the mesh density increases (not shown), as expected for approximations of $\vec{\nabla}w_{k}$. Figure 3(b) shows the total charge current density at the contact, and the results with and without the offset representation agree to 5 digits. We conclude that the offset representation allows robust extraction of internal current densities but does not seem to be important for the overall current density of the test device. The offset representation imposes essentially no extra computational cost on Simudo while enabling robust determination of internal current densities.

3.4 Optics

The optical problem is solved by self-consistently iterating through the optical flux variables $\Phi_{[\lambda_{i,\textnormal{min}},\lambda_{i,\textnormal{max}}],\hat{s}_{i}}$ and independently solving Eq. 19 for each one. For convenience, we write $\Phi=\Phi_{[\lambda_{i,\textnormal{min}},\lambda_{i,\textnormal{max}}],\hat{s}_{i}}$, $\hat{s}=\hat{s}_{i}$, and $\alpha=\alpha_{[\lambda_{i,\textnormal{min}},\lambda_{i,\textnormal{max}}]}$ for the remainder of this section. We represent $\Phi$ using CG basis functions of order $d_{\textnormal{optical}}=2$.

We now derive the weak form used in Simudo to solve each optical propagation problem. We follow closely the derivation in Zhao2013a of the modified second order radiative transfer equation (MSORTE) method, without the scattering matrix. Integrating Eq. 19 with a test function $v\in\mathrm{CG}(d_{\textnormal{optical}})$ gives

[TABLE]

where $\vec{v}=v\,\hat{s}$ and $\zeta=\hat{s}\cdot\vec{\nabla}\Phi$. Using Eq. 29, we obtain

[TABLE]

Inserting the outlet boundary condition Eq. 20b into the first term, we obtain the final weak form

[TABLE]

The inlet boundary condition Eq. 20a is applied directly on $\Phi$ as an essential boundary condition.

3.5 Sentaurus benchmark comparison

To validate Simudo, we benchmark it against the industry standard Synopsys Sentaurus device simulator. Since Sentaurus does not support intermediate band materials, the benchmark is limited to standard semiconductors. Our test problem is the same silicon pn-junction as considered in Section 3.3.2, with the overall $J(V)$ as shown in Fig. 3(b). See below for a discussion of the differences in implementation of the Ohmic condition between Simudo and Sentaurus.

We study the convergence of the results as the mesh is refined, using a number of mesh points in the x-direction ranging from 45 to 12289, with the same meshes used for both Sentaurus and Simudo. The mesh spacing is nonuniform in the x-direction since the carrier and current densities vary most rapidly near the contacts and the junction. The meshes are generated by splitting the structure into two regions, one extending from the p-contact to the junction and one from the n-contact to the junction. In each region, a mesh spacing $d_{0}$ is applied to the cell adjacent to the contact and the cell adjacent to the junction. The mesh spacing increases geometrically toward the center of each region with a growth factor of 1.2. To generate the finest mesh sizes, these cells are further subdivided into 4, 16 or 64 equal parts. There are 2 points in the y-direction. The computational cost generally increases with the number of degrees of freedom rather than with the number of cells in the mesh. Simudo has more degrees of freedom associated with each cell than Sentaurus, due to its higher-order basis functions. For this mesh, Simudo has 36 degrees of freedom per triangle while Sentaurus has 9 per triangle, with 2 triangles per mesh point. In this mesh, each triangle has at least one edge on the boundary of the device, which increases the number of degrees of freedom per triangle compared to a mesh where most triangles share sides; this effect is similar for both methods.

Figure 4(a) shows the results of the study, where we plot the relative error in $J(V)$ for each simulation, with the reference current $J_{\text{ref}}$ taken from Sentaturus with the densest mesh. With 193 mesh points, Simudo converges approximately as well as Sentaurus with 3073 mesh points. Above 769 mesh points, the Simudo results show no further improvement in error, indicating either that Simudo and Sentaurus converge to results that differ at the $10^{-5}$ level or that Sentaurus with 12289 mesh points is only converged to the $10^{-5}$ level, while Simudo may have converged more precisely. The simulations were performed on different machines, so we do not report timing data.

The figure makes clear that Simudo converges much more rapidly with mesh size than Sentaurus, which demonstrates the higher order convergence that FEM is supposed to provide over FVM. Figure 4(b) shows the scaling of the errors with the number of degrees of freedom in each simulation, at 0.4 V. The solid blue line shows that Sentaurus’ self-convergence scales like $N^{-2}$ with the number of degrees of freedom. The dashed green line shows Simudo’s self-convergence, with $J_{\text{ref}}$ taken from the Simudo simulation with 12289 mesh points. It shows that Simudo’s self-convergence scales like $N^{-4}$ with the number of degrees of freedom. For all but the smallest meshes, Simudo’s convergence is superior to Sentaurus’ at the same number of degrees of freedom. Taken together, these figures show that with 193 mesh points, Simudo’s result is as good as Sentaurus’ with 64 times as many mesh points, which is equivalent to 16 times as many degrees of freedom.

Note that the boundary conditions at the contacts are not precisely the same for the Sentaurus and Simudo simulations. Both are intended to simulate Ohmic contacts for the majority carrier and surface recombination velocities of 0 for the minority carrier. The Simudo simulations are performed with surface recombination velocity $S=\infty$ and 0 for the majority and minority carriers, respectively, imposing equilibrium carrier concentration at the boundary for majority carriers and setting $\vec{j}_{k}\cdot\hat{n}|_{\Gamma}=0$ for minority carriers, as described in Sec. 2.2.2. The Sentaurus simulations are performed with $S=0$ for the minority carriers, in agreement with Simudo, and the default “Ohmic contact” boundary condition for the majority carriers, which imposes charge neutrality and equilibrium carrier concentration at the contact. Under small and reverse bias, these two sets of boundary conditions should be equivalent, but under large forward bias, the default Sentaurus boundary condition is expected to give incorrect results due to its imposition of charge neutrality Synopsys15 . Sentaurus provides a “Modified Ohmic” boundary condition, which should be closer to the Simudo boundary condition, but we were unable to attain convergence using it. As a result, at larger biases the Simudo and Sentaurus results diverge from each other, and we do not include them in Fig. 4. For biases larger than 1 V, the diode is in high injection, and the Boltzmann approximation used in this calculation is not accurate, regardless.

4 Examples and results

In this section, we give examples of using Simudo. Section 4.1 shows how to set up a simple 1-dimensional pn-junction device and demonstrates the helpful tools that Simudo provides for defining regions and boundaries. Section 4.2 shows the extensibility of Simudo by illustrating the code required to add a new Auger recombination process. Section 4.3 illustrates the use of Simudo to study a system first considered in Ref. Marti02 .

4.1 PN junction and topology definitions

We include in the supplementary material the code listing equilibrium.py describing a simple pn-junction device in Simudo. This example constructs the device and implements steps 1-4 of the pre-solver shown in Fig. 2. Here, we discuss some of the pieces of that code and illustrate the useful topology construction operations built in to Simudo.

In the 1-dimensional pn-junction example, the object ls contains information about the layers, including their sizes, positions, and mesh. The object pdd sets up the Poisson/drift-diffusion solver and has information about the bands in each material, including recombination processes and boundary conditions. In this example, there are only two bands (VB, CB); for a problem including an IB, pdd would have a third band, too.

Simudo is designed for 2-dimensional simulations, and it has sophisticated tools to define the arrangement of materials, dopings, contacts, meshing regions, or other user-defined spatial properties. In many FEM solvers, interfaces must be tracked manually, including their orientation, to ensure that integrals over those interfaces are added together properly. Simudo introduces a set of topology tools that instead allow users to define the regions and interfaces in which they are interested, and Simudo takes care of all the bookkeeping. The user defines regions as desired (e.g., emitter, base, defective-region), which can then be given properties, whether they be doping levels, recombination parameters, or other desired properties. These regions are initially defined abstractly, without having any coordinates in the device, using CellRegions and FacetRegions, and are later connected to geometry and materials by the mesh generator.

Full details are given in the documentation accompanying Simudo, but we give a further illustration of these methods in Figs. 5 and 6. That example illustrates the creation of arbitrary CellRegion objects, including unions and intersections, and edges that connect them. When R is a CellRegions container, accessing a nonexistent attribute (such as R.domain) causes its creation. The user can define new CellRegion objects by applying Boolean operations on previous ones and new FacetRegion objects by using the boundary method. For example, consider the region R.region1. Then

R.region1.boundary(R.region2) creates a signed boundary from region1 to region2, as illustrated in Fig. 5. All of these custom regions are kept as symbolic expressions and evaluated by Simudo only when needed (e.g., when asked to apply a boundary condition or when asked to compute a volume or surface integral). This layer of abstraction allows the user not to worry about the details of mesh markers, entity indices, and facet orientations LoggMardalEtAl2012a , and is described more fully in the documentation that accompanies Simudo.

The examples in Figs. 5-6 illustrate another useful concept. The mesh generation interprets the external region as being outside the simulation domain, allowing convenient definitions for boundary conditions and current flow. The FacetRegions are used in the pn-junction example shown in the supplementary material to define the boundary conditions, which – in step 2 – are conductive at the left and right contacts and nonconductive at the top and bottom surfaces. That example also shows how the mesh can be refined by adding extra mesh points near the contacts.

4.2 Extensibility: Adding Auger recombination

The initial release of Simudo contains radiative and Shockley-Read trapping and recombination processes in the non-degenerate limits for VB, CB. The user can easily add modified physics to their problems, which we demonstrate here with an example of adding an Auger recombination process to Simudo, with the form

[TABLE]

where $C_{n},C_{p}$ are the Auger coefficients, and $p_{0}$ and $n_{0}$ are the hole and electron concentrations at thermal equilibrium, respectively Nelson03 . The code is listed in Fig. 7. The function get_generation_user(band) adds a negative local generation rate in the CB and VB and returns 0 for all other bands. This recombination process moves particles between two bands, the src_band and the dst_band. In this case, where the electrons and holes have opposite charge, the Auger process destroys both particles simultaneously; if both carrier types involved in the process had the same charge (e.g., for a CB-to-IB trapping process), the process would represent a particle-number-conserving transfer (rather than a recombination) from the src_band to the dst_band, with the appropriate sign for the recombination process determined by the get_band_generation_sign method (inherited from TwoBandEOPMixin). This method’s sign convention is that the dst_band always gains carriers through the generation process, while the src_band gains or loses as required by conservation of charge.

4.3 P[IB]N junction

In Ref. Marti2002a , the authors consider a quantum-dot-based IB solar cell with a p-n-IB-p-n structure. They present a drift-diffusion model for the IB region only, with the carrier density and current density boundary conditions obtained from a depletion approximation and law of the junction. This model assumes that transport is diffusion-dominated in the IB region, and drift can therefore be neglected. This early device model gave important insights into the behavior of IB devices.

In testing the self-consistency of the model, the authors estimate the IB mobility required to remain in the diffusion-dominated regime, finding that an IB mobility greater than $62\mathrm{\>cm^{2}/V/s}$ is required to make their model consistent. This claim raises an immediate question: does something interesting happen when the IB mobility goes below that threshold? Since Simudo is a full drift-diffusion device model, we can directly answer that question.

We model a similar device with a simpler p-IB-n structure. This device has the same band and absorption parameters as the device described in Marti2002a , summarized in Table 3. The incident light is a blackbody spectrum at $6000\mathrm{\>K}$ with a solar concentration factor $X=1000$. The device has equal subgap optical cross sections and nearly current matched (within 10%) incident photon fluxes for the two subgap transitions. In this case, the local IV and CI generation rates are nearly identical throughout the IB device. The code to set up this problem is included in the supplementary material, in marti2002.py.

We simulate this device with $\mu_{I}$ ranging from $0.001$ to $300\mathrm{\>cm^{2}/V/s}$, with resulting J-V curves shown as solid lines in Fig. 8a, which is tightly zoomed and still shows only minor effects of this over-$10^{5}$ change in $\mu_{I}$. In fact, the IB and drifts currents contribute negligibly to the transport inside the IB region and the current remains diffusion-dominated throughout, as shown in Fig. 9a.

The device behavior is approximately independent of $\mu_{I}$ because the IV and CI generation rates are roughly equal at each point in the device, so IB transport is not required to enable sub-gap current matching, but this behavior is not generic for all IB devices. We illustrate this effect by increasing $\sigma_{CI}^{\textnormal{opt}}$ by a factor of five (while keeping $\sigma_{IV}^{\textnormal{opt}}$ unchanged). In this case, the device is still globally approximately current-matched, but the CI absorption process occurs preferentially at the top of the device while the IV generation occurs deeper in the IB region. The device thus relies on IB transport for the CI and IV generation rates to balance over the full device. These effects are shown in the dashed curves of Fig. 8, which show a stronger dependence on $\mu_{I}$ than in the matched case. When IB mobility is low, the excess CI generation in the front of the device instead causes local CI trapping, with equivalent local IV trapping toward the back of the device, reducing overall current. In the high-mobility case, the overall current is slightly larger with the mismatched absorptions, due to the increased optical depth. Figure 8b shows how $J_{\text{sc}}$ varies with $\mu_{I}$ in both of these cases, where the greater dependence on $\mu_{I}$ in the mismatched case is apparent. Figure 9b shows that in the matched-$\sigma^{\text{opt}}$ case, $j_{I}$ is never particularly large, while it grows to be five times larger in the mismatched case, showing the role of IB currents in internally balancing the optical absorptions.

In the low-mobility limit, where $j_{I}$ is always small, when local CI and IV current generations are imbalanced, the filling fraction $f$ of the IB must shift to equalize generation and recombination at each point Strandberg09 . This effect is visible in Fig. 10, where at low mobility, the mismatched-$\sigma^{\text{opt }}$ case has photodepletion at the front side and photofilling at the back side, consistent with excess CI generation at the front and excess IV generation at the back. Both the matched-$\sigma^{\text{opt }}$ and the high-mobility mismatched-$\sigma^{\textnormal{opt}}$ cases maintain an approximately uniform IB filling fraction.

These examples together show the utility of Simudo to explore the performance of IB devices and resolve an assertion made in earlier device models without the benefit of a coupled PDD/optics solver.

5 Conclusion

The availability of a device model for intermediate band materials should enable both understanding of this new class of materials and optimization of IB devices. Simudo’s use of the FEM and its methods for overcoming catastrophic cancellation may also prove useful in standard semiconductor device simulation. Simudo has been validated against Synopsys Sentaurus for standard semiconductor devices and shown to converge more rapidly with mesh size. This self-consistent solution of the Poisson/drift-diffusion and optical propagation equations provides a platform for studying a wide range of optoelectronic materials and devices, including solar cells and photodetectors, with tools to enable extensibility to arbitrary generation and recombination models, thermal effects, and more. The near-term roadmap for Simudo includes explicit heterojunction support and nonlocal tunneling, which will be available with future releases at github.com/simudo/simudo. We hope that the free and open source nature of this software will enable further development of IB materials and device simulation more broadly.

Acknowledgments

We acknowledge funding from US Army Research Laboratory (W911NF-16-2-0167), the Natural Sciences and Engineering Research Council of Canada TOP-SET training program, and computing resources from Compute Canada. We thank Emily Zinnia Zhang for alpha testing Simudo, contributing the first code implementing trapping processes, and valuable conversations.