Drift–Diffusion Theory: From Boltzmann Transport to Energy Balance

This page presents a first-principles derivation of the drift–diffusion equations used in semiconductor device simulation, starting from the Boltzmann Transport Equation in the relaxation-time approximation. By taking successive moments of the Boltzmann equation, we show how the electron and hole continuity equations, drift–diffusion current relations, and the energy-transport (hot-carrier) equations arise in a single, self-consistent framework. Particular emphasis is placed on heterojunctions, band-edge driving, density-of-states effects, and thermal transport, leading to the solver-ready equations implemented in OghmaNano.

1. Introduction

OghmaNano’s electrical engine is a 1D/2D/3D drift–diffusion framework whose defining feature is its support for dynamic (non-equilibrium) trap states. Rather than assuming traps instantaneously equilibrate with the free carriers, OghmaNano can evolve trap occupancies explicitly in both energy and space, which is essential for correctly modelling disordered semiconductors and transient measurements (e.g. ToF, CELIV) as well as steady-state operation.

There are many legitimate routes to drift–diffusion. For example, one can derive it from the Boltzmann Transport Equation via moment expansions and controlled closures, from irreversible thermodynamics (Onsager / entropy-production arguments), or from a random-walk / Fokker–Planck viewpoint that connects microscopic hopping to macroscopic drift and diffusion. These approaches often produce equations that look similar at first glance, but they are not equally reliable when you move beyond the simplest “homogeneous semiconductor” case.

In particular, once you introduce heterojunctions, spatially varying effective density of states, position-dependent effective mass, band-edge gradients, or non-trivial statistics, it is no longer safe to bolt on extra terms to a drift–diffusion current and hope the result remains self-consistent. The same is true when you begin to include thermal driving or couple electrical transport to heat generation. What you want is a derivation framework that naturally generates the correct extra terms and provides a clear path to extensions, rather than a collection of ad hoc modifications.

For that reason, the derivation presented here starts from the Boltzmann Transport Equation. Taking moments of the BTE yields, in a systematic way, the continuity equations and the drift–diffusion constitutive relations as a low-inertia limit, and it also shows how the same framework “lifts” to the next level: the energy-balance (energy-transport) equations. Even if you do not solve the full hydrodynamic model, the energy framework is valuable because it identifies consistent electrical heat-source terms and clarifies when carrier heating must be accounted for. The derivation below forms the foundation of OghmaNano's electrical model.

2. Boltzmann Transport Equation (RTA)

At the most fundamental level used in device modelling, charge transport is described in terms of a distribution function \(f(\mathbf{r},\mathbf{k},t)\). This function represents the probability of finding a charge carrier at position \(\mathbf{r}\), with crystal momentum \(\hbar\mathbf{k}\), at time \(t\). All macroscopic electrical quantities — carrier density, current density, energy density — can be obtained by taking suitable moments of this distribution.

The Boltzmann Transport Equation (BTE) is the equation of motion for this distribution function. It is a conservation law in phase space: it accounts for how carriers move through real space, how their momentum changes under applied forces, and how scattering processes redistribute carriers in momentum space. Starting from the BTE therefore provides a single, unified framework from which drift–diffusion, energy transport, and related models can all be derived in a consistent way.

In its full form the collision (scattering) term of the BTE is complicated and material-specific. For practical device modelling it is common to use the relaxation-time approximation (RTA), in which scattering is assumed to drive the distribution toward a local quasi-equilibrium form \(f^0\) over a characteristic time \(\tau\). This approximation preserves the essential physics of momentum and energy relaxation while keeping the equations tractable.

With this approximation, the semiclassical Boltzmann transport equation can be written as

\[ \frac{\partial f}{\partial t} + \mathbf{v}\cdot \nabla_{\mathbf{r}} f + \frac{\mathbf{F}}{\hbar}\cdot \nabla_{\mathbf{k}} f = -\frac{f - f^0}{\tau}. \]

The three terms on the left-hand side have a clear physical meaning: the first describes the explicit time evolution of the distribution, the second describes transport through real space with carrier velocity \(\mathbf{v}\), and the third describes acceleration in momentum space due to an applied force \(\mathbf{F}\) (for example \(\mathbf{F}=-q\mathbf{E}\) in an electric field). The right-hand side represents scattering, which relaxes the system toward \(f^0\).

Drift–diffusion theory does not attempt to solve this equation directly. Instead, it proceeds by taking moments of the BTE — integrals over momentum space — to obtain evolution equations for physically meaningful quantities such as carrier density, current density, and energy density. The following sections show how this procedure leads naturally to the familiar drift–diffusion equations and, at the next level, to energy-balance (hot-carrier) models.

3. Taking moments of the Boltzmann Transport Equation

The Boltzmann Transport Equation describes carrier dynamics in terms of a distribution function \(f(\mathbf{r},\mathbf{k},t)\), which gives the occupation of electronic states as a function of position, momentum, and time. In equilibrium this distribution reduces to the familiar Fermi–Dirac function, while under bias it evolves in response to electric fields, band-edge variations, and scattering processes. To obtain equations that are useful at the device scale, we do not attempt to solve directly for \(f(\mathbf{r},\mathbf{k},t)\). Instead, we derive evolution equations for macroscopic quantities by taking moments of the entire Boltzmann equation, systematically connecting the microscopic carrier statistics to drift–diffusion, energy transport, and related models.

Formally, a moment equation is obtained by multiplying the full BTE by a weight \(A(\mathbf{k})\) and integrating over all momentum space:

\[ \int A(\mathbf{k}) \left( \frac{\partial f}{\partial t} + \mathbf{v}\cdot\nabla_{\mathbf{r}} f + \frac{\mathbf{F}}{\hbar}\cdot\nabla_{\mathbf{k}} f \right) \mathrm{d}^3k = \int A(\mathbf{k})\left(-\frac{f-f^0}{\tau}\right)\mathrm{d}^3k. \]

Each choice of the weighting function \(A(\mathbf{k})\) generates a balance equation: an evolution equation expressing how a macroscopic quantity is driven and relaxed. Setting \(A=1\) yields the particle balance equation, from which the familiar electron and hole continuity equations are obtained. Taking the first moment, with \(A=\hbar\mathbf{k}\) (or equivalently \(m^\ast\mathbf{v}\)), produces a momentum balance equation that describes how carrier momentum responds to forces and scattering; the standard drift–diffusion current is recovered by taking the overdamped limit of this equation, in which momentum relaxes rapidly and inertial terms are neglected. Taking the next moment, with \(A=W(\mathbf{k})\), gives an energy balance equation, which describes the transport and relaxation of carrier energy and provides the minimal extension needed to model hot carriers, thermal driving, and electrical heat generation. In this way, continuity, drift–diffusion, and energy-transport models emerge as successive levels of approximation within a single, self-consistent framework.

Standard drift–diffusion retains only the particle balance equation and a simplified constitutive relation for the current obtained from the momentum balance equation. In doing so, it does not conserve carrier energy.

4. Zeroth moment: particle balance (continuity equation)

The first and most important moment of the Boltzmann Transport Equation is obtained by setting the weighting function to unity, \(A(\mathbf{k}) = 1\). This corresponds to counting carriers: integrating the full Boltzmann equation over all momenta gives a balance law for the total carrier density at each point in space. For this reason the zeroth moment is called the particle balance (or particle continuity) equation.

4.1 Derivation (zeroth moment / particle balance)

Start from the relaxation-time form of the Boltzmann equation:

\[ \frac{\partial f}{\partial t} + \mathbf{v}\cdot \nabla_{\mathbf{r}} f + \frac{\mathbf{F}}{\hbar}\cdot \nabla_{\mathbf{k}} f = -\frac{f - f^0}{\tau}. \]

Multiply by \(A(\mathbf{k})=1\) and integrate over all \(\mathbf{k}\):

\[ \int \left( \frac{\partial f}{\partial t} + \mathbf{v}\cdot \nabla_{\mathbf{r}} f + \frac{\mathbf{F}}{\hbar}\cdot \nabla_{\mathbf{k}} f \right)\mathrm{d}^3k = -\int \frac{f - f^0}{\tau}\,\mathrm{d}^3k. \]

The first term defines the carrier density \[ n(\mathbf{r},t) = \int f(\mathbf{r},\mathbf{k},t)\,\mathrm{d}^3k, \] so that \[ \int \frac{\partial f}{\partial t}\,\mathrm{d}^3k = \frac{\partial n}{\partial t}. \]

The second term becomes a divergence in real space:

\[ \int \mathbf{v}\cdot \nabla_{\mathbf{r}} f \,\mathrm{d}^3k = \nabla_{\mathbf{r}}\cdot \int \mathbf{v}\, f \,\mathrm{d}^3k \equiv \nabla\cdot(n\mathbf{u}), \]

where the average carrier velocity is \[ \mathbf{u}(\mathbf{r},t)=\frac{1}{n}\int \mathbf{v}(\mathbf{k})\, f(\mathbf{r},\mathbf{k},t)\,\mathrm{d}^3k. \]

The third term (force term) vanishes under standard assumptions that the distribution decays rapidly as \(|\mathbf{k}|\rightarrow\infty\), so the corresponding \(\mathbf{k}\)-space surface integral is zero:

\[ \int \nabla_{\mathbf{k}}\cdot(\cdots)\,\mathrm{d}^3k \approx 0. \]

Collecting terms, the zeroth-moment equation becomes

\[ \frac{\partial n}{\partial t} + \nabla\cdot(n\mathbf{u}) = -\int \frac{f - f^0}{\tau}\,\mathrm{d}^3k. \]

The right-hand side represents processes that create or remove free carriers when viewed at the drift–diffusion level: generation, recombination, and (in materials with traps) exchange with trap states. We therefore write it compactly as \(G - R\) (with trap capture/emission included in the effective \(R\) term or written explicitly in the trapping model).

4.2 Continuity equation in drift–diffusion form

Introducing the electron current density \[ \mathbf{J}_n = -q\,n\,\mathbf{u} \] the particle balance equation becomes the familiar electron continuity equation:

\[ \frac{\partial n}{\partial t} = \frac{1}{q}\nabla\cdot \mathbf{J}_n + G - R. \]

The corresponding hole continuity equation is

\[ \frac{\partial p}{\partial t} = -\frac{1}{q}\nabla\cdot \mathbf{J}_p + G - R. \]

5. First moment: momentum balance → drift–diffusion current

The drift–diffusion current equations are obtained from the first moment of the Boltzmann Transport Equation. This moment corresponds to a momentum balance, describing how carrier momentum is driven by forces such as electric fields and band-edge gradients, and reduced by scattering. Rather than solving the full momentum balance equation, standard drift–diffusion is obtained by expressing the carrier velocity directly in terms of local driving forces, without retaining an explicit evolution equation for momentum or energy. This approximation is appropriate whenever higher-order momentum and energy transport effects are not required.

5.1 From the first moment of the BTE to the drift–diffusion current

We start from the Boltzmann Transport Equation in the relaxation-time approximation, written with a momentum-relaxation time \(\tau_p\):

\[ \frac{\partial f}{\partial t} + \mathbf{v}\cdot \nabla_{\mathbf{r}} f + \frac{\mathbf{F}}{\hbar}\cdot \nabla_{\mathbf{k}} f = -\frac{f - f^0}{\tau_p}. \]

To obtain the momentum balance, multiply the entire equation by the momentum-like weight \(A(\mathbf{k}) = m^\ast \mathbf{v}\) (equivalently \(\hbar\mathbf{k}\) for a parabolic band), and integrate over momentum space:

\[ \int m^\ast\mathbf{v} \left( \frac{\partial f}{\partial t} + \mathbf{v}\cdot \nabla_{\mathbf{r}} f + \frac{\mathbf{F}}{\hbar}\cdot \nabla_{\mathbf{k}} f \right)\mathrm{d}^3k = -\int m^\ast\mathbf{v}\,\frac{f-f^0}{\tau_p}\,\mathrm{d}^3k. \]

5.2 Identifying the velocity moment, current, and pressure

The first-moment integrals introduce the quantities that appear in device-scale transport equations. We begin by defining the carrier density and the mean carrier velocity as

\[ n = \int f\,\mathrm{d}^3k, \qquad \mathbf{u} = \frac{1}{n}\int \mathbf{v} f\,\mathrm{d}^3k. \]

The velocity \(\mathbf{u}\) defined above is the mean (or average) carrier velocity. It represents the net drift of the carrier population at a given point in space. Individual carriers typically move much faster than \(\mathbf{u}\), executing random thermal motion superimposed on this slow collective drift. The mean velocity is therefore not the typical speed of a carrier, but the small residual velocity that remains after averaging over all microscopic motion.

This distinction is crucial. Random thermal motion does not contribute to electrical current because it averages to zero, whereas the mean drift velocity \(\mathbf{u}\) captures the imbalance in carrier motion caused by electric fields, concentration gradients, and temperature gradients. It is this mean velocity that determines the macroscopic current density.

With these definitions, the particle flux is \(n\mathbf{u} = \int \mathbf{v} f\,\mathrm{d}^3k\), and the electron current density is

\[ \mathbf{J}_n = -q\,n\,\mathbf{u}. \]

Substituting these definitions back into the first-moment equation allows the momentum-space integrals to be rewritten entirely in terms of \(n\), \(\mathbf{u}\), and higher-order velocity moments. In particular, the collision term simplifies because the reference distribution \(f^0\) represents a local equilibrium state with zero mean drift velocity, so that

\[ \int m^\ast \mathbf{v}\, f^0 \,\mathrm{d}^3k = 0. \]

As a result, the collision term contributes only a momentum-relaxation term proportional to the mean velocity \(\mathbf{u}\).

The remaining spatial-transport term in the first-moment equation contains integrals of the form \(\int m^\ast \mathbf{v}\mathbf{v}\, f\,\mathrm{d}^3k\), which represent the flux of momentum through space. To make the physical content of this term explicit, we decompose the carrier velocity into a mean part and a fluctuation:

\[ \mathbf{v} = \mathbf{u} + (\mathbf{v}-\mathbf{u}). \]

Substituting this decomposition into the second velocity moment and using \(\int (\mathbf{v}-\mathbf{u}) f\,\mathrm{d}^3k = 0\) splits the momentum flux into a convective contribution and a fluctuation term. The fluctuation contribution is identified as the carrier pressure (or stress) tensor,

\[ \mathbf{P} = m^\ast \int (\mathbf{v}-\mathbf{u})(\mathbf{v}-\mathbf{u})\, f \,\mathrm{d}^3k. \]

Physically, \(\mathbf{P}\) represents the transport of momentum associated with velocity fluctuations about the mean drift velocity. Its divergence gives rise to diffusion and thermal driving terms in the reduced transport equations.

Collecting all contributions, the first-moment equation can now be written as

\[ \frac{\partial}{\partial t}(n m^\ast \mathbf{u}) + \nabla \cdot \mathbf{P} + n\,\mathbf{F} = -\frac{n m^\ast \mathbf{u}}{\tau_p}. \]

This equation is the general momentum balance obtained directly from the first moment of the Boltzmann Transport Equation. No drift–diffusion assumptions have yet been made; these enter in the next step, where the momentum balance is reduced to a local force balance.

5.3 Reducing the momentum balance to a force balance

In drift–diffusion modelling, we are not interested in explicitly resolving the time evolution or spatial transport of momentum itself. Instead, we assume that carrier momentum adjusts locally to the applied forces, so that any transient or non-local momentum transport effects can be neglected. Under this assumption, the momentum balance simplifies by dropping the explicit momentum-transport and inertia terms, while keeping the force and relaxation terms in place.

With this approximation, the momentum balance reduces to

\[ n\,\mathbf{F} - \nabla \cdot \mathbf{P} = \frac{n m^\ast \mathbf{u}}{\tau_p}. \]

Written in this way, the equation should be read as a local force balance: the forces driving carrier motion are balanced by momentum loss to scattering.

The remaining term \(\nabla \cdot \mathbf{P}\) represents how momentum is redistributed by random carrier motion. The tensor \(\mathbf{P}\) itself encodes correlations between velocity components, and in general allows for anisotropic momentum transport. In many semiconductor devices, the carrier distribution is close to isotropic in velocity space. In this case, there is no preferred direction for the velocity fluctuations, and the pressure tensor reduces to a scalar pressure multiplying the identity tensor:

\[ \mathbf{P} \approx p\,\mathbf{I}. \]

Physically, this means that random thermal motion contributes equally in all directions. When this form is substituted into the momentum balance, the divergence of the pressure tensor simplifies to the gradient of a scalar pressure:

\[ \nabla \cdot \mathbf{P} = \nabla \cdot (p\,\mathbf{I}) = \nabla p. \]

This term is the origin of diffusion and thermal driving in the drift–diffusion equations. In the next section, this force-balance equation is solved for the mean carrier velocity \(\mathbf{u}\), leading directly to the drift–diffusion current.

5.4 From the force balance to the drift–diffusion current

We now return to the reduced local momentum balance derived in the previous section. For electrons, it is convenient to express the driving force in terms of the conduction-band edge \(E_c(\mathbf{r})\), which incorporates both electrostatic potential and material offsets. We therefore write the force on electrons as

\[ \mathbf{F} = -\nabla E_c, \]

where \(E_c(\mathbf{r})\) may be written (up to a reference) as \(E_c = \chi - q\phi\), so that spatial variations in either electron affinity or electrostatic potential contribute naturally to carrier driving. Substituting this into the reduced momentum balance gives

\[ -\,n\,\nabla E_c - \nabla p = \frac{n m^\ast}{\tau_p}\,\mathbf{u}. \]

This equation expresses a local balance between band-edge driving, pressure-driven transport, and momentum loss to scattering. Solving for the mean carrier velocity yields

\[ \mathbf{u} = -\frac{\tau_p}{m^\ast} \left( \nabla E_c + \frac{1}{n}\nabla p \right). \]

Using the definition of the electron current density, \(\mathbf{J}_n = -q n \mathbf{u}\), we obtain

\[ \mathbf{J}_n = q\,\frac{q\tau_p}{m^\ast}\,n\,\nabla E_c + q\,\frac{\tau_p}{m^\ast}\,\nabla p. \]

Introducing the mobility \(\mu_n \equiv q\tau_p/m^\ast\), this becomes

\[ \mathbf{J}_n = q\mu_n n\,\nabla E_c + \mu_n\,\nabla p. \]

At this point, the current is expressed in terms of the pressure gradient \(\nabla p\), and no assumption has yet been made about carrier statistics. To proceed, we express the pressure and density as moments of the underlying distribution function.

The carrier density and pressure may be written as

\[ n = \int g(E)\, f(E)\,\mathrm{d}E, \qquad p = \frac{2}{3} \int (E - E_c)\, g(E)\, f(E)\,\mathrm{d}E, \]

where \(g(E)\) is the density of states and \(f(E)\) is the Fermi–Dirac distribution. For a parabolic band, these integrals reduce to standard Fermi–Dirac integrals,

\[ n = N_c\, F_{1/2}(\eta), \qquad p = n\, k_B T\, \frac{F_{3/2}(\eta)}{F_{1/2}(\eta)}, \]

with reduced chemical potential \(\eta = (E_{Fn}-E_c)/(k_B T)\), and \(F_j(\eta)\) the complete Fermi–Dirac integral of order \(j\).

Taking the gradient of the pressure gives

\[ \nabla p = k_B T\, \frac{F_{3/2}(\eta)}{F_{1/2}(\eta)}\,\nabla n \;+\; n\,k_B T\, \nabla\!\left( \frac{F_{3/2}(\eta)}{F_{1/2}(\eta)} \right) \;+\; n\,k_B \frac{F_{3/2}(\eta)}{F_{1/2}(\eta)}\,\nabla T. \]

Substituting this expression into the current shows that diffusion and thermal driving are governed by a generalized Einstein relation,

\[ D_n = \frac{\mu_n k_B T}{q}\, \frac{F_{3/2}(\eta)}{F_{1/2}(\eta)}. \]

This expression is exact for a parabolic band and remains valid from the non-degenerate to the strongly degenerate regime. In the non-degenerate limit (\(\eta \ll -1\)), \(F_{3/2}(\eta)/F_{1/2}(\eta) \rightarrow 1\), and the generalized Einstein relation reduces to the familiar form

\[ D_n = \frac{\mu_n k_B T}{q}. \]

The first term represents drift driven by spatial variations of the conduction-band edge, including both electrostatic fields and heterojunction band offsets. The second term, involving the pressure gradient, is the origin of diffusion and thermal driving.

In the non-degenerate case, the carrier pressure is \(p = n k_B T\), so that

\[ \nabla p = k_B T\,\nabla n + n k_B \nabla T. \]

Substituting this expression into the current gives

\[ \mathbf{J}_n = q\mu_n n\,\nabla E_c + \mu_n k_B T\,\nabla n + \mu_n n k_B \nabla T. \]

Using the Einstein relation for non-degenerate carriers, \(D_n = \mu_n k_B T / q\), the current can be written in the familiar drift–diffusion form with explicit thermal driving:

\[ \mathbf{J}_n = \underbrace{q\mu_n n\,\nabla E_c}_{\text{band-edge drift}} \;+\; \underbrace{q D_n \nabla n}_{\text{diffusion}} \;+\; \underbrace{q D_n \frac{n}{T}\nabla T}_{\text{thermal driving}}. \]

In heterostructures, the effective density of states \(N_c(\mathbf{r},T)\) may vary spatially due to changes in bandstructure parameters such as effective mass. A self-consistent reduction then produces an additional material-driving term proportional to \(\nabla\ln N_c\).

Including this contribution, the drift–diffusion current used in device simulators can be written as

\[ \mathbf{J}_n = \underbrace{q\mu_n n\,\nabla E_c}_{\text{band-edge drift}} \;+\; \underbrace{q D_n \nabla n}_{\text{diffusion}} \;+\; \underbrace{q D_n \frac{n}{T}\nabla T}_{\text{thermal driving}} \;-\; \underbrace{q D_n n\,\nabla\ln N_c}_{\text{DOS / effective-mass driving}}. \]

Written in this form, all driving mechanisms arise naturally from the same momentum-balance framework. Electrostatic fields, heterojunction band offsets, temperature gradients, and spatial variations in material parameters all enter on exactly the same footing.

An entirely analogous derivation applies to holes. Repeating the same steps using the valence-band edge \(E_v(\mathbf{r})\) and hole statistics yields the hole current density

\[ \mathbf{J}_p = \underbrace{q\mu_p p\,\nabla E_v}_{\text{band-edge drift}} \;-\; \underbrace{q D_p \nabla p}_{\text{diffusion}} \;-\; \underbrace{q D_p \frac{p}{T}\nabla T}_{\text{thermal driving}} \;+\; \underbrace{q D_p p\,\nabla\ln N_v}_{\text{DOS / effective-mass driving}}, \]

where \(N_v(\mathbf{r},T)\) is the effective density of states of the valence band. The opposite signs of the diffusion and thermal-driving terms reflect the positive charge of holes.

• Thermal driving (thermodiffusion / Seebeck-like term): spatial variations in carrier temperature modify the local quasi-equilibrium distribution and generate a net carrier flux proportional to \(\nabla T\), even in the absence of an electric field.
• Density-of-states / effective-mass driving: in heterostructures or spatially varying materials, gradients in bandstructure parameters (such as effective mass) lead to gradients in \(N_c\) and \(N_v\). These terms are essential at material interfaces and ensure that carrier transport remains consistent with the underlying statistics and bandstructure.

7. Energy transport / hot-carrier extension

In the drift–diffusion model derived earlier, carrier transport is described by particle continuity and a momentum balance reduced to an algebraic current law. Carrier energy is assumed to relax rapidly to the lattice, so no explicit equation is solved for it.

When this assumption is relaxed, the next equation in the Boltzmann moment hierarchy must be retained: the energy balance equation. This equation governs how carriers gain energy from electric fields and band-edge gradients, transport that energy through the device, and lose it to the lattice.

7.1 Second moment of the Boltzmann equation

The energy balance is obtained by multiplying the Boltzmann Transport Equation by the single-particle energy \(W(\mathbf{k})\) and integrating over momentum space. For electrons, the energy of a state is written as

\[ W(\mathbf{k},\mathbf{r}) = E_c(\mathbf{r}) + \varepsilon(\mathbf{k}), \]

where \(E_c(\mathbf{r})\) is the conduction-band edge and \(\varepsilon(\mathbf{k})\) is the kinetic energy relative to the band edge (for a parabolic band, \(\varepsilon=\hbar^2 k^2/2m^\ast\)).

Defining the mean carrier energy per particle as

\[ \bar W = \frac{1}{n}\int W f\,\mathrm{d}^3k, \]

the corresponding energy density is \(n\bar W\). This quantity plays the same role for energy that \(n\) plays for particle number.

7.2 Energy flux and relation to drift–diffusion currents

The spatial-transport term in the second moment contains integrals of the form \(\int W\,\mathbf{v} f\,\mathrm{d}^3k\), which motivates the definition of the energy flux

\[ \mathbf{q} \equiv \int W\,\mathbf{v} f\,\mathrm{d}^3k. \]

This quantity is directly analogous to the particle flux \(n\mathbf{u} = \int \mathbf{v} f\,\mathrm{d}^3k\), and therefore to the current density \(\mathbf{J}_n = -q n\mathbf{u}\). Energy transport is thus naturally coupled to carrier transport.

7.3 Work done by the electric field and band-edge gradients

The force term in the Boltzmann equation produces a distinct contribution in the second moment. Using integration by parts in momentum space and the identity \(\nabla_{\mathbf{k}} W = \hbar\mathbf{v}\), one finds

\[ \int W\,\frac{\mathbf{F}}{\hbar}\cdot\nabla_{\mathbf{k}} f\,\mathrm{d}^3k = -\,\mathbf{F}\cdot\int \mathbf{v} f\,\mathrm{d}^3k = -\,n\,\mathbf{u}\cdot\mathbf{F}. \]

Writing the force on electrons as \(\mathbf{F}=-\nabla E_c\), this term becomes

\[ n\,\mathbf{u}\cdot\nabla E_c = -\frac{1}{q}\,\mathbf{J}_n\cdot\nabla E_c. \]

This is the electrical work done on the carrier population. In the special case where \(E_c=-q\phi\), it reduces to the familiar Joule heating term \(\mathbf{J}_n\cdot\mathbf{E}\).

7.4 Energy balance equation and connection to drift–diffusion

Collecting all contributions, the second-moment equation can be written as

\[ \frac{\partial}{\partial t}(n\bar W) + \nabla\cdot\mathbf{q} - \frac{1}{q}\,\mathbf{J}_n\cdot\nabla E_c = -\frac{n(\bar W-\bar W_0)}{\tau_W}. \]

This equation is directly coupled to the drift–diffusion model through the current density \(\mathbf{J}_n\). No additional ad hoc heating terms are required: electrical power dissipation emerges naturally from the same band-edge driving that appears in the current law.

7.5 Relation to the standard drift–diffusion limit

To see how this connects to familiar drift–diffusion models, consider the non-degenerate, near-equilibrium case. The kinetic contribution to the mean energy is

\[ \bar W - E_c = \frac{3}{2}k_B T_e, \]

so that

\[ n\bar W = nE_c + \frac{3}{2}n k_B T_e. \]

In steady-state drift–diffusion with rapid energy relaxation, \(T_e \approx T_L\) and \(\bar W \approx \bar W_0\), the time derivative \(\partial_t(n\bar W)\) vanishes and the energy equation collapses to a balance between electrical heating and energy loss to the lattice.

Retaining the energy balance equation lifts this restriction. The carrier temperature \(T_e\) becomes a dynamic variable, allowing the model to capture field heating, velocity overshoot, and non-equilibrium transport while remaining fully compatible with the drift–diffusion framework derived earlier.

Putting everything together, the energy balance equation implemented in a drift–diffusion–based solver can be written in the following form:

\[ \frac{\partial}{\partial t} \!\left( \frac{3}{2} n k_B T_e \right) \;+\; \nabla\!\cdot\! \left( \frac{3}{2} k_B T_e\,\frac{\mathbf{J}_n}{q} - \kappa \nabla T_e \right) \;-\; \frac{1}{q}\,\mathbf{J}_n\cdot\nabla E_c \;=\; -\,\frac{3}{2}\frac{n k_B (T_e - T_L)}{\tau_W}. \]

Here \(T_e\) is the carrier (electron) temperature, \(T_L\) is the lattice temperature, \(\kappa\) is the carrier thermal conductivity, and \(\tau_W\) is the energy-relaxation time. This equation is solved self-consistently with Poisson’s equation, the carrier continuity equations, and the drift–diffusion current relations.

8. Final coupled model solved in OghmaNano

OghmaNano solves a coupled system of partial differential equations. The exact set of equations depends on the selected physical model (drift–diffusion only, or drift–diffusion with energy transport). The table below summarises the equations and when they are used.