Discontinuous Galerkin Deterministic Solvers for a Boltzmann-Poisson Model of Hot
Electron Transport by Averaged Empirical Pseudopotential Band Structures

1

José Morales-Escalante 2 , Irene M. Gamba3
ICES and Department of Mathematics, University of Texas, Austin, TX 78712

arXiv:1512.05403v2 [] 17 Jan 2018

Yingda Cheng4
Department of Mathematics, Michigan State University, East Lansing, MI 48824
Armando Majorana5
Department of Industrial Engineering, University of Catania, Catania, Italy
Chi-Wang Shu6
Division of Applied Mathematics, Brown University, Providence, RI 02912
and
James Chelikowsky

7

Department of Physics and ICES, University of Texas, Austin, TX 78712
February 2017

1

Support from the Institute of Computational Engineering and Sciences and the University of Texas Austin is
gratefully acknowledged.
2
E-mail: jmorales@ices.utexas.edu
3
E-mail: gamba@math.utexas.edu. JME and IMG are partially supported by NSF DMS 1109625, NSF CHE0934450, NSF-RNMS DMS-1107465 and the ICES Moncrief Grand Challenge Award.
4
E-mail: ycheng@math.msu.edu. Research supported by NSF grants DMS-1318186 and DMS-1453661.
5
E-mail:majorana@dii.unict.it. Research supported by J.T. Oden Visitors Program, ICES, The University of
Texas at Austin.
6
E-mail: shu@dam.brown.edu. Research supported by DOE grant DE-FG02-08ER25863 and NSF grant DMS1418750.
7
E-mail: jrc@ices.utexas.edu. JRC acknowledges support provided by the Scientific Discovery through Advanced
Computing (SciDAC) program funded by U.S. DOE, Office of Science, Advanced Scientific Computing Research and
Basic Energy Sciences under Award No. de-sc0008877.

1

Abstract
The purpose of this work is to incorporate numerically, in a discontinuous Galerkin (DG) solver
of a Boltzmann-Poisson model for hot electron transport, an electronic conduction band whose
values are obtained by the spherical averaging of the full band structure given by a local empirical
pseudopotential method (EPM) around a local minimum of the conduction band for silicon, as a
midpoint between a radial band model and an anisotropic full band, in order to provide a more
accurate physical description of the electron group velocity and conduction energy band structure
in a semiconductor. This gives a better quantitative description of the transport and collision
phenomena that fundamentally define the behaviour of the Boltzmann - Poisson model for electron
transport used in this work. The numerical values of the derivatives of this conduction energy
band, needed for the description of the electron group velocity, are obtained by means of a cubic
spline interpolation. The EPM-Boltzmann-Poisson transport with this spherically averaged EPM
calculated energy surface is numerically simulated and compared to the output of traditional analytic band models such as the parabolic and Kane bands, numerically implemented too, for the
case of 1D n+ − n − n+ silicon diodes with 400nm and 50nm channels. Quantitative differences
are observed in the kinetic moments related to the conduction energy band used, such as mean
velocity, average energy, and electric current (momentum).

Keywords: Deterministic numerical methods; discontinuous Galerkin schemes; Boltzmann-Poisson
systems; empirical pseudopotential method; statistical hot electron transport; semiconductor nano
scale devices

2

1

Introduction

The Boltzmann-Poisson (BP) system is a semi-classical model for electric charge transport in semiconductors. The BP system can be used to describe the hot electron transport in modern semiconductor devices at nano-scales. As stated in [1], this model describes the long range interactions over
charge carriers and the statistical evolution of its states that includes an account of the quantum
scattering events. The BP system treats charge carriers partly as classical particles by describing
them by means of a time-dependent probability density function f (t, x, k) over the phase space
(x, k), and using a Boltzmann equation to model the time evolution of the associated probability
density function in the phase space. The quantum nature of the carriers is considered in several
terms of the Boltzmann equation. The quantum crystal wave-vector k is used as the momentum
phase space variable in the model. The model for the local velocity of the charge carriers is their
respective group velocity v(k) = ~1 ∇k ε(k), related to the electronic energy band function ε(k) of
the considered semiconductor material. As usual, ~ is the Planck constant divided by 2π. The
collision integral operator models the quantum scattering mechanisms acting over the charge carriers. The flow of charge carriers is induced by the force over the electron charge −q, which is
assumed to be given by a mean electric field, F(t, x) = −q E(t, x). This effective electric field,
modeled by the Poisson Equation, takes into account long range interactions made of both internal
carrier self-consistent and external contributions, such as an applied potential (bias). Hence, timedependent solutions of the the BP system contain all the information on the transient of the carrier
distribution and the time evolution of the total electric field. A phenomenological derivation of the
BP model can be found in [1].
The semi-classical Boltzmann description of electron transport in semiconductors is, for a truly
3-D device, an equation in six dimensions plus time when the device is not in steady state. For a
1-D device model, under azimuthal symmetry assumptions, the dimensionality of the problem can
be reduced to 3 dimensions plus time. The heavy computational cost is the main reason why the
BP system had been traditionally solved numerically by means of Direct Simulation Monte Carlo
(DSMC) methods [3]. However, after the pioneer work [4], in recent years, deterministic solvers to
3

the BP system were proposed in [5, 6, 7, 8, 9, 10, 11, 12]. These methods provide accurate results
which, in general, agree well with those obtained from Monte Carlo (DSMC) simulations, often at
a fractional computational time. Moreover, these type of solvers can resolve transient details for
the electron probability density function f , which are difficult to compute with DSMC simulators.
The initial methods proposed in [7, 8, 9, 10] using weighted essentially non-oscillatory (WENO)
finite difference schemes to solve the Boltzmann-Poisson system, had the advantage that the scheme
is relatively simple to code and very stable even on coarse meshes for solutions containing sharp
gradient regions. However, a disadvantage of the WENO methods is that it requires smooth meshes
to achieve high order accuracy, hence it is not very flexible for adaptive meshes.
Motivated by the easy hp-adaptivity (h refers to the size of the element, and p refers to the polynomial degree of the space generated by the basis functions) and the simple communication pattern
of the discontinuous Galerkin (DG) methods for macroscopic (fluid level) models [13, 14, 15, 16], it
was proposed in [17, 18] to implement a DG solver to the full Boltzmann equation, that is capable
of capturing transients of the probability density function. In the previous work [17, 18], the first
DG solver for the BP system was proposed, and some numerical calculations were shown for one
and two-dimensional devices. In [19], the DG-LDG scheme for the Boltzmann-Poisson system was
carefully formulated, and extensive numerical studies were performed to validate the calculations.
Such scheme models electron transport along the conduction band for 1D diodes and 2D double
gate MOSFET devices with the energy band ε(k) = ε(|k|) given by the Kane band model (valid
close to a local minimum) in which the relation between the energy ε and the wavevector norm |k|
is given by the analytic formula, referred as the Kane band model:
ε (1 + αε) =

~2 |k|2
2m∗

(1)

where m∗ is the effective mass for the considered material, Silicon for the case of this paper, and α
is a non-parabolicity constant. This band model can be understood as a first order variation from
the parabolic band model, given by the particular case α = 0.
In all of the aforementioned deterministic solvers previous to [20], the energy-band function ε(k)
is given analytically, either by the parabolic band approximation or by the Kane non-parabolic
4

band model. The analytical band makes use of the explicit dependence of the carrier energy on
the quasimomentum, which significantly simplifies all expressions as well as implementation of
these techniques in the collision operator. However, some physical details of the band structure
are partly or totally ignored when using an analytic approximation, which hinders its application
to transport of hot carriers in high-field phenomena (the so called hot electron transport) where
the high anisotropy of the real band structure far from the conduction band minimum becomes
important. Full band models, on the other hand, are able to provide an accurate physical description
of the energy-band function, portraying this anisotropic band structure far from a conduction
band minimum. One of the most commonly used methods to compute full bands is the empirical
pseudopotential method (EPM). Such method gives a full band structure truncating the Fourier
series in the k-space [22] for a crystal lattice potential model given as the sum of potentials due
to individual atoms and associated electrons, with few parameters fitting empirical data such as
optical gaps, absorption rates, etc, to finally compute the energy eigenvalues of the Schrödinger
equation in Fourier space. A more detailed discussion of this method can be found in [22, 23]. While
full band models, as the ones given by EPM, have been widely used in DSMC simulators [3], their
inclusion in deterministic solvers for the transport Boltzmann Equation is more recent; on [24],
[25], full band models have also been combined with spherical harmonic expansion methods used
to solve the Boltzmann equation numerically. However, high order accuracy is not always achieved
by spherical harmonic expansion methods when energies vary strongly and only a few terms of the
expansion are usually employed [26]. In contrast, the simulations for the BP system developed in
our line of work, as in [17], [19], do not involve any asymptotics and so are very accurate for hot
electron transport regimes. A DG method for full conduction bands BP models was proposed in
[20], generalizing the solver that uses the Kane non-parabolic band and adapting it to treat the
full energy band case. A preliminary benchmark of numerical results shows that the the Dirac
delta functional in the energy variable can be applied in this case to reduce one dimension of the
collision integral, and so an accurate high-order simulation with comparable computational cost to
the analytic band cases is possible.

5

The work presented in this paper is focused on simulations for hot electron transport along a single
conduction band for Si computed by radially averaging an EPM full band structure. The band
obtained by this procedure represents a midpoint between a radial band model and a full band
anisotropic model, having then both the desired advantages of a band model with a dependence
on r and, at the same time having information of the anisotropic variation of the conduction
band in the k-space, by means of the numerical average performed over the angular domain of
the conduction band for a given k-sphere. The advantage of DG scheme in this framework is
shown in the accurate calculation of the Dirac delta functions in the collisional integrals based
on the weak formulation. The EPM-Boltzmann-Poisson transport with this spherically averaged
EPM calculated energy surface is numerically simulated and compared to the output of traditional
analytic band models such as the parabolic and Kane bands. Quantitative differences are observed
in the moments, demonstrating the significance of incorporating the physical band models in the
kinetic simulations.
The rest of the paper is organized as follows: in Sections 2 and 3, we present the BP model, and
the transformed equations under spherical coordinates in k. Section 4 contains the details of the
computations of the spherical average of a local EPM conduction band. The DG formulation is
presented in Section 5. The device configuration and numerical results are discussed in Sections 6
and 7. We conclude the paper in Section 8. Some technical details of the schemes are given in the
Appendix.

2

The Boltzmann-Poisson problem

We consider the probability density function (pdf ) for electrons along a single conduction band,
denoting it by f (t, x, k) We denote by Ωx the physical domain in the x-space, and similarly Ωk
as the domain in the k-space. Following [1], we recall the classical (strong) formulation of the
initial value problem for the BP system with boundary conditions, for a pdf of electrons on a single
conduction band:
Find f : R+ × Ωx × Ωk → R,

f (t, x, k) ≥ 0 and V (t, x) : R+ × Ωx → R, such that the Boltzmann

6

equation
∂f
1
q
+ ∇k ε(k) · ∇x f − E(t, x) · ∇k f = Q(f ) ,
∂t
~
~
with the linear collision operator Q(f ) describing the scattering over the electrons, where several
quantum mechanisms can be taken into account. In the low density approximation, the collisional
integral operator becomes linear in f , having the form:
Z
Q(f ) =



S(k0 , k)f (t, x, k0 ) − S(k, k0 )f (t, x, k) dk0

(2)

Ωk

where S(k, k0 ) is the scattering kernel, representing non-local interactions of electrons with a background density distribution.
In the case of silicon, for example, one of the most important collision mechanisms are electronphonon scatterings due to lattice vibrations of the crystal, which are modeled by acoustic (assumed
elastic) and optical (non-elastic) non-polar modes, the latter with a single frequency ωp , given by:
S(k, k0 ) = (nq + 1) K δ(ε(k0 ) − ε(k) + ~ωp )
+ nq K δ(ε(k0 ) − ε(k) − ~ωp ) + K0 δ(ε(k0 ) − ε(k)) ,

(3)

with K, K0 constants for silicon.
The symbol δ indicates the usual Dirac delta distribution corresponding to the well known Fermi’s
Golden Rule [2]. The constant nq is related to the phonon occupation factor:

−1


~ωp
−1
,
nq = exp
kB TL
where kB is the Boltzmann constant and TL = 300K is the constant lattice temperature. In the
Boltzmann equation, the energy band function ε(k) are often taken by a simple band model, e.g. the
Kane band model (1) or the parabolic band model. Since this term appears in both the transport
and collision part, it is evident that it plays an important role in the numerical simulation. We will
discuss more about the treatment of this term in later sections.
The self consistent field E(t, x) is solved from the Poisson equation
∇x · [r (x) ∇x V (t, x)] =

q
[ρ(t, x) − ND (x)] ,
0
7

E(t, x) = −∇x V (t, x)

where the parameter 0 is the dielectric constant in a vacuum, r (x) labels the relative dielectric
function which depends on the material. ρ(t, x), the electron charge density, is given by the integral
over the domain in the k-space Ωk :
Z
ρ(t, x) =

f (t, x, k) dk ,

(4)

Ωk

and ND (x) is the doping profile, representing an external fixed density of positive charge carriers.
The evolution is subject to the initial condition
f (0, x, k) = f0 (x, k) ∀ (x, k) ∈ Ωx × Ωk , t = 0,
and suitable boundary conditions for f on ∂Ωx × Ωk and Ωk × ∂Ωx , and for V on ∂Ωx are satisfied.
N
The boundary ∂Ωx is usually split for the Poisson Equation in Dirichlet ∂ΩD
x , Neumann ∂Ωx , and
I
N
Interface boundaries ∂ΩIx , such that ∂Ωx = ∂ΩD
x ∪ ∂Ωx ∪ ∂Ωx .

Examples of boundary conditions used for the Boltzmann Eq. include [1]:
• Charge neutrality [35], [10], [18], [19]:

fout (t, x, k) =

ND (x) fin (t, x, k)
ρin (t, x)

t ≥ 0, x ∈ ∂ΩD
x , k ∈ Ωk

(5)

This condition is usually employed at the device contacts (Dirichlet boundaries ∂ΩD
x ).
• Null x-flux:
n(x) · ∇x f (t, x, k) = 0

t ≥ 0, x ∈ ∂ΩN
x , k ∈ Ωk ,

(6)

where n(x) is the normal to the surface ∂ΩN
x at the point x. This condition is imposed on
the part of the physical domain with an insulating layer (Neumann boundaries ∂ΩN
x ).
• Vanishing boundary conditions in the k-space:
f (t, x, k) = 0

t ≥ 0, x ∈ Ωx , k ∈ ∂Ωk .

(7)

These conditions correspond to negligible densities for large energy values. We use these
vanishing conditions for the Boltzmann Equation in our work. We will just mention that, if
8

we had chosen Ωk as the first Brillouin zone, then periodic boundary conditions in the k-space
would be the correct physical conditions. However, it is difficult to apply these conditions on
the complex shape of the boundary of a truncated octahedron, which is the shape of the first
Brillouin zone for Silicon and Germanium crystals.
Boundary conditions related to the Poisson Equation could be:
• Applied potential (bias):
V (t, x) = V0 (t, x)

t ≥ 0, x ∈ ∂ΩD
x

(8)

This condition is imposed where we have device contacts (Dirichlet boundaries).
• Neumann boundary conditions for the electric potential:
n(x) · ∇x V (t, x) = 0

t ≥ 0 , x ∈ ∂ΩN
x

(9)

where n(x) is the normal to the surface ∂ΩN
x at the point x. This condition is imposed on
the part of the physical domain with an insulating layer, which is a Neumann boundary.
It is important to mention that the contact boundaries ∂ΩD
x for the Boltzmann and Poisson
equations must be the same.

3

Boltzmann equation in spherical coordinates for the k-vector

We show here the Boltzmann equation with the momentum k in spherical coordinates presented
in [20]. As opposed to the previous work in [19], the coordinate transformation based on the
Kane analytic band relation proposed in [5] can no longer be used for an energy band that does
not assume this analytic Kane band model and that takes into account anisotropy for ε(k). The
spherical coordinate system is used in k space instead of Cartesian coordinates because of the
higher resolution demands near the conduction band minimum (chosen as the origin k = 0), and
large cells in k-space are sufficient for describing the tail of the distribution function accurately.
The following change of variables into dimensionless quantities are introduced for a general problem:
√ ∗

p
t
x
2m kB TL √  p
t = , (~x, z) =
, k=
r µ, 1 − µ2 cos ϕ, 1 − µ2 sin ϕ
t∗
`∗
~
9

with r ≥ 0 , µ ∈ [−1, 1] , ϕ ∈ [−π, π].
1
ε(k)
kB TL
V (t∗ t, `∗ x, `∗ y, `∗ z)
,
Ψ(t, x, y, z) =
V∗
V∗
with cv =
and E∗ = 0.1 V∗ `−1
∗ .
`∗ E∗
ε(r, µ, ϕ) =

E = −cv ∇x V

where the spherical coordinate transformation maps the k-domain Ωk onto the set Ω of the (r, µ, ϕ)
space. Typical values for length, time and voltage are given by `∗ = 10−6 m, t∗ = 10−12 s and
V∗ = 1 Volt, respectively.
Thus, a new unknown “weighted” pdf function Φ is obtained by multiplying the pdf f by the
Jacobian of the spherical k-transformation:
√

r
f (t, x, y, z, r, µ, ϕ) ,
2

Φ(t, x, y, z, r, µ, ϕ) =

(10)

which can be interpreted as the probability density function of an electron being in the neighborhood
of the phase-space state (x, y, z, r, µ, ϕ) at time t.
Hence, writing the collisional integral in spherical coordinates and multiplying the Boltzmann equation by the Jacobian associated to the k-spherical transformation, yields the following Transformed
Boltzmann Equation (TBE) for the unknown Φ:

∂
∂
∂
∂
∂
∂
∂Φ
+
(a1 Φ) +
(a2 Φ) +
(a3 Φ) +
(a4 Φ) +
(a5 Φ) +
(a6 Φ) = C(Φ)
∂t
∂x
∂y
∂z
∂r
∂µ
∂ϕ

10

(11)

with the following transport terms ~a = (a1 , a2 , a3 , a4 , a5 , a6 )T :


√
∂ε 1 − µ2 ∂ε
a1 (·) = cD 2 r µ
,
+ √
∂r
r ∂µ
!
p
2 cos ϕ ∂ε
√ p
∂ε
1
−
µ
∂ε
µ
sin
ϕ
√
,
a2 (·) = cD 2 r 1 − µ2 cos ϕ
−
−√ p
∂r
∂µ
r
r 1 − µ2 ∂ϕ
!
p
2 sin ϕ ∂ε
√ p
µ
cos
ϕ
∂ε
1
−
µ
∂ε
√
−
+√ p
,
a3 (·) = cD 2 r 1 − µ2 sin ϕ
∂r
∂µ
r
r 1 − µ2 ∂ϕ
i
p
√ h
a4 (·) = −2 cE r µ Ex (t, x, y, z) + 1 − µ2 (cos ϕ Ey (t, x, y, z) + sin ϕ Ez (t, x, y, z)) ,
#
"
p
µ 1 − µ2
1 − µ2
√ Ex (t, x, y, z) −
√
(cos ϕ Ey (t, x, y, z) + sin ϕ Ez (t, x, y, z)) ,
a5 (·) = −cE
r
r
1
a6 (·) = −cE √ p
[ − sin ϕ Ey (t, x, y, z) + cos ϕ Ez (t, x, y, z)] ,
r 1 − µ2
and the linear collision operator
√ Z
r
S(r0 , µ0 , ϕ0 , r, µ, ϕ) Φ(t, x, y, z, r0 , µ0 , ϕ0 ) dr0 dµ0 dϕ0
C(Φ)(t, x, y, z, r, µ, ϕ) =
2 Ω
√
Z
r0 0 0 0
0 0
0
− Φ(t, x, y, z, r, µ, ϕ) S(r, µ, ϕ, r , µ , ϕ )
dr dµ dϕ ,
2
Ω

(12)

where the scattering kernel is
S(r, µ, ϕ, r0 , µ0 , ϕ0 ) = c+ δ(ε(r0 , µ0 , ϕ0 ) − ε(r, µ, ϕ) + αp )
+ c− δ(ε(r0 , µ0 , ϕ0 ) − ε(r, µ, ϕ) − αp ) + c0 δ(ε(r0 , µ0 , ϕ0 ) − ε(r, µ, ϕ)),
accounting for acoustic and optical electron-phonon interaction, the main scattering mechanisms
in silicon. The constants above are defined as
r

~ωp
t∗ qE∗
kB TL
, cE = √ ∗
, αp =
,
∗
2m
kB TL
2m kB TL
2m∗ t∗ p
(c+ , c− , c0 ) =
2 m∗ kB TL [(nq + 1)K, nq K, K0 ] .
~3
t∗
cD =
`∗

The details of the derivation of the TBE are collected in the appendix. The dimensionless Poisson
equation is
∂
∂x







∂Ψ
∂
∂Ψ
∂
∂Ψ
r
+
r
+
r
= cp [ρ(t, x, y, z, t) − ND (x, y, z)]
∂x
∂y
∂y
∂z
∂z

11

(13)

where
√
ND (x, y, z) =
Z
ρ(t, x, y, z) =

2 m∗ kB TL
~

√

−3
ND (`∗ x, `∗ y, `∗ z),

cp =

2 m∗ kB TL
~

3

`2∗ q
,
0 V∗

Φ(t, x, y, z, r0 , µ0 , ϕ0 ) dr0 dµ0 dϕ0 .

Ω

3.1

Geometrical interpretation of the force terms in the TBE

Although the terms (a1 , a2 , a3 ) related to the transport in the x-space due to the electron group
velocity in the TBE can be easily interpreted as just the gradient ∇k ε(k) expressed in spherical
coordinates, the terms (a4 , a5 , a6 ) related to the transport in the k-space due to the electric field
might be more obscure to understand. A simple expression for them can be identified.
a4 = −2 cE


p
√  p
√
r µ, 1 − µ2 cos ϕ, 1 − µ2 sin ϕ · E = −2 cE r êr · E,

p

1 − µ 2 p
1 − µ2
2
1 − µ , −µ cos ϕ, −µ sin ϕ · E = −cE √
a5 = −cE √
êµ · E,
r
r
1
1
a6 = −cE √ p
êϕ · E.
(0, − sin ϕ, cos ϕ) · E = −cE √ p
r 1 − µ2
r 1 − µ2

(14)

p

(15)
(16)

These transport terms express the acceleration field induced by E in spherical coordinates, as they
are related to the negatives of the directional cosines of E with respect to the unit vectors êr , êµ ,
êϕ . Hence, the TBE (11) is written in conservative, divergence form, as a flow in the k-space due to
the electric field decomposed in each of the orthogonal components of the spherical k-coordinates.
This can be easily derived from the expression for the divergence in general curvilinear coordinates,
applied to the particular case of spherical coordinates k(r, µ, ϕ). This calculation is included in the
Appendix.

4

Computation of the spherical average of a local EPM conduction
band for silicon

The motivation of this work is to incorporate numerically, in a DG solver of the BP system electronic
conduction bands whose values are obtained by the radial averaging of the full band structure given
by a local empirical pseudopotential method (EPM) around a local minimum of the conduction
band for silicon. By performing the radial averaging, it simplifies the discussion of the numerical
12

Figure 1: Local EPM conduction energy band-structure (ε) color plot in the k-space 1st octant
enclosing the Silicon Brillouin Zone. Conduction band local minimum: ko = (0.8562, 0, 0)(2π/a)
scheme significantly. This is done as a midpoint between a radial and an anisotropic full energy
band models, with the goal of providing a more accurate physical description of the electron group
velocity and of the scattering mechanisms by Fermi Golden Rule, and consequently improve the
transport and electron - phonon collision phenomena. The approximation of the electron group
velocity is obtained from the numerical values of the derivatives of the conduction band, which are
obtained by means of a cubic spline interpolation. The numerical values of the spherically averaged
EPM band and the derivatives are obtained as described below.
A local empirical pseudopotential method (EPM) code developed by Chelikowsky et al. [22] is
adapted to compute the conduction band structure of silicon in its Brillouin Zone in the k-space.
The local pseudopotentials are used in this EPM code to mimic a silicon semiconductor with crystal
diamond structure [23].
A color plot of the local EPM conduction band on the first octant of the k-space enclosing the
Brillouin Zone for silicon is shown in Fig. 1.
13

The calculated EPM band structure ε(kx , ky , kz ) = ε(k(r, µ, ϕ)) is then averaged over the k-spheres
rk around the local energy minimum point k0 = (0.8562, 0, 0)2π/a (where a is the lattice constant
for silicon) by means of a 10 point Gaussian quadrature on the angular space. Using the symmetry
of the silicon conduction band, the integration only needs to be performed in the (µ, ϕ) domain
[0, 1] × [0, π]
Z
ε̃(rk ) =

0

1Z π
0
Z

ε(rk , µ, ϕ) dµdϕ
10 X
10
X
≈
ωm ωn ε(rk , µm , ϕn ) .
Z
1
π
m=1
n=1
dµ dϕ

0

(17)

0

The values of the radius of these k-spheres are the grid points rk in the DG-BP simulations. In
this way we obtain a band model that has a dependence on r, and at the same time it uses the
information of the anisotropic energy band values in the angular k-domain via its numerical average.
As a midpoint between a radial band model and a full band anisotropic model, it has the desired
advantages of both. A cubic spline interpolation is then performed, using the numerical values of
the radial average ε̃(r) at the midpoints of the r-cells, and the derivative of this spline interpolation
is used to obtain a numerical approximation of the derivative dε̃/dr at these r-midpoints.
The spherical averages of the EPM conduction band ε̃(r) vs r ∝ |k − k0 |2 with the related spline
interpolation for Si are shown in Fig. 2 (in red). The parabolic (blue), which is a linear function
of r, and the Kane (green) analytic conduction band models for silicon are plotted as well.
It can be observed that there is a quantitative difference between the different energy band models.
The spherical average of the EPM band is below the Kane band model, which is below the Parabolic
band.

We show in Fig. 3 the relative l2 error norm of the spherical average EPM band with respect to
the local EPM data ε(r, µ, ϕ) as a function of r, given by the formula

[ε − ε̃]2
(rk ) ≈
hε2 i

P10

P10
ωm ωn [ε(rk , µm , ϕn ) − ε̃(rk )]2
m=1
P10 n=1
P10
2
m=1
n=1 ωm ωn [ε(rk , µm , ϕn )] .

14

Plot: Comparison of Energy (ε) vs r=|k−ko|2 for the EPMravg, Parabolic, and Kane Band Models
1
Parabolic
Kane, α=0.5/eV
Eravg

0.9
0.8
0.7

ε (eV)

0.6
0.5
0.4
0.3
0.2
0.1
0

0

0.005

0.01

0.015

0.02

0.025 0.03
r = |k−ko|2

0.035

0.04

0.045

0.05

Figure 2: Energy (E) vs r = |k − k0 |2 for the parabolic, Kane, & symmetrized EPM spherical
average band.

15

It can be observed in Fig. 3 that the relative l2 error increases with r, which indicates that far
away from the local minimum k0 the anisotropy of the conduction band becomes increasingly more
important.

5

DG formulation for the TBE and the Poisson equation

In this section, we will discuss the DG schemes based on the radially symmetric band. The scheme
is implemented based on a piecewise linear polynomial approximation.

5.1

Domain and Finite Element Space

Let’s consider a 2 D rectangular domain in the physical space and a rectangular domain Ωk in
momentum space. We use simple rectangular cells
h
i h
i
Ωijkmn = xi− 1 , xi+ 1 × yj− 1 , yj+ 1 × Kkmn
2

2

2

2

where
i h
i h
i
h
Kkmn = rk− 1 , rk+ 1 × µm− 1 , µm+ 1 × ϕn− 1 , ϕn+ 1
2

2

2

2

2

2

with
xi± 1 = xi ±
2

∆yj
∆xi
∆rk
, yj± 1 = yj ±
, rk± 1 = rk ±
··· ,
2
2
2
2
2

and i = 1, ..., Nx , j = 1, ..., Ny , k = 1, ...Nr , m = 1, ...Nµ , n = 1, ...., Nϕ .
The test functions ψ(x, y, r, µ, ϕ) belong to the linear function space

Vh1 = v : v|Ωijkmn ∈ P 1 (Ωijkmn ) ,
where P 1 (Ωijkmn ) is the set of polynomials of degree at most 1 on the cell Ωijkmn .
A set of piecewise linear basis functions for Vh1 in the open cell Ω̊ijkmn is given by

1, 2

(x − xi ) (y − yj ) (r − rk ) (µ − µm ) (ϕ − ϕn )
,2
,2
,2
,2
∆xi
∆yj
∆rk
∆µm
∆ϕn

16


(18)

Rel. l2−Norm Error of Radial avg. Conduction Band | a = 10.2610 A.U., ko = (0.8562,0,0).
0.18
0.16
0.14

Erelerr

0.12
0.1
0.08
0.06
0.04
0.02
0

0

0.005

0.01

0.015

0.02

0.025 0.03
r = |k−ko|2

0.035

0.04

0.045

0.05

Figure 3: Relative l2 -deviation in the angular space of EPM average from the EPM full band for
Si, as a function of r = |k − k0 |2

17

Hence, in the cell Ω̊ijkmn , we approximate our weighted-pdf Φ by a piecewise polynomial Φh of first
degree in Vh1 :
2(y − yj )
2(x − xi )
+ Yijkmn (t)
∆xi
∆yj
2(ϕ − ϕn )
2(r − rk )
2(µ − µm )
+ Pijkmn (t)
.
(19)
+ Rijkmn (t)
+ Mijkmn (t)
∆rk
∆µm
∆ϕn
h
i h
i
The charge density on xi− 1 , xi+ 1 × yj− 1 , yj+ 1 under this piecewise linear approximation is
Φh (t, x, y, r, µ, ϕ) = Tijkmn (t) + Xijkmn (t)

2

2

2

2


Nµ Nϕ 
Nr X
X
X
2(y − yj )
2(x − xi )
+ Yijkmn (t)
∆rk ∆µm ∆ϕn .
Tijkmn (t) + Xijkmn (t)
ρh (t, x, y) =
∆xi
∆yj
k=1 m=1 n=1

(20)
The problem is then reduced to find, by means of our numerical scheme, the unknowns:

Tijkmn (t), Xijkmn (t), Yijkmn (t), Rijkmn (t), Mijkmn (t), Pijkmn (t).

5.2

(21)

Discontinuous Galerkin Formulation for the TBE

The corresponding weak DG formulation and its corresponding approximation consists on finding
Φh ∈ Vh1 , such that for any test function vh ∈ Vh1 and a generic cell K of the decomposition of
Ωx × Ωk , solves:
Z
∂vh
∂vh
(a1 Φh ) dσ −
(a2 Φh ) dσ −
K
K ∂x
K ∂y
Z
Z
Z
∂vh
∂vh
∂vh
(a4 Φh ) dσ −
(a5 Φh ) dσ −
(a6 Φh ) dσ +
K ∂r
K ∂µ
K ∂ϕ
Z
Fx+ − Fx− + Fy+ − Fy− + Fr+ − Fr− + Fµ+ − Fµ− + Fϕ+ − Fϕ− =
C(Φh )vh dσ,
Z

∂Φh
vh dσ −
∂t

Z

(22)

K

where vh is a test function in Vh1 , dσ = dx dy dr dµ dϕ, and where the F ± ’s terms are boundary
integrals over four-dimensional boundary surfaces associated to each 5-dimensional volume element

18

Ωijkmn , that is:
Fx±

yj+ 1

Z

2

=

yj− 1
2

Fy± =

2

2

2

xi− 1
2

=
Z

=

xi+ 1
2

xi− 1
2

Z2

yj− 1
2

yj− 1
2

2

Z

Z2

Z2

rk+ 1
2

rk− 1

ϕn+ 1

2

Z

ϕn+ 1
2

ϕn− 1

a4 Φ̂h vh∓

2

dydrdµdϕ

yj± 1

dxdrdµdϕ

rk± 1

dxdydµdϕ

2

a5 Φ̂h vh∓

a6 Φ̂h vh∓

µm± 1

dxdydrdϕ

2

2

µm+ 1

µm− 1

2

a2 Φ̂h vh∓

2

2

xi± 1

2

ϕn− 1

rk+ 1

rk− 1

2

2

µm+ 1

Z

a1 Φ̂h vh∓

2

ϕn+ 1

ϕn− 1

2

2

yj+ 1

2

Z

2

2

Z

µm− 1

yj+ 1

Z

xi− 1

Fϕ±

2

ϕn+ 1

ϕn− 1

2

yj+ 1
2

xi+ 1

Z

Z

µm− 1

yj− 1

2

Fµ±

2

Z

2

2

2

xi+ 1

µm+ 1

µm+ 1

Z

rk− 1

2

Fr± =

Z

µm− 1

rk+ 1

Z

xi− 1

Z

2

rk− 1

xi+ 1

Z

rk+ 1

Z

ϕn± 1

dxdydrdµ

.

2

2

The values for vh∓ are the ones for the function vh on the interior of the considered cell. The upwind
numerical flux Φ̂h defines the value of Φh at the boundary. That means Φh might be discontinuous
at the boundary.
The collisional terms
Z
C(Φh )vh dσ
K

become a linear combination, with numerical constant coefficients, of the the unknowns (21) which
are precomputed and stored. The Poisson equation can be solved by either an integral formula,
projecting the solution to the electric field into the space Vh1 , for the 1D device case, or by means
of a LDG method, for higher dimensional cases. A Runge Kutta method is applied for the time
evolution of the time dependent coefficients (21) for the piecewise linear approximation Φh ∈ Vh1 .

5.3

Calculation of the Collision terms in the DG formulation

Here, we will show the details of the calculation of the collisional integrals since it is the most
demanding part of the simulation. For completeness, the calculation of the transport term is also
reported, and collected in the Appendix B.
h
i h
i h
i
Denote by Kkmn = rk− 1 , rk+ 1 × µm− 1 , µm+ 1 × ϕn− 1 , ϕn+ 1 the rectangular cells in the
2

2

2

2

2

2

spherical coordinates for k-space. Because the collisional operators only perform integrations in
p
k-space, it is convenient to write the basis functions in (18) as the product of two functions ηi,j
(x, y)

19

p
and ξk,m,n
(r, µ, ϕ), which are given in Ωijkmn = ΩI by



n
o
p
ηi,j
(~x)

p=0,1,..,5

=


n
o
p
ξk,m,n
(~r )

p=0,1,..,5

=

2(x − xi ) 2(y − yj )
1, 1, 1, 1 ,
,
∆xi
∆yj


,

(23)


2(r − rk ) 2(µ − µm ) 2(ϕ − ϕn )
1,
,
, 1, 1 ,
,
∆rk
∆µm
∆ϕn

(24)

where we define
~x = (x, y) ,

~r = (r, µ, ϕ) ,

~r 0 = (r0 , µ0 , ϕ0 ) ,

d ~r = dr dµ dϕ

I = (i, j, k, m, n)

1 if (~x, ~r ) ∈ Ω̊I
χI = χI (~x, ~r ) =
0 otherwise

(25)
(26)
(27)

and
WI0 (t) := TI (t) ,

WI1 (t) := RI (t) ,

WI2 (t) := MI (t) ,

WI4 (t) := XI (t) ,

WI5 (t) := YI (t) .

WI3 (t) := PI (t) ,

(28)
(29)

Then, in a piecewise continuous linear approximation of Φ, we have (almost everywhere), that


5
X
X
p
p
Φ(t, ~x, ~r ) =
χI (~x, ~r ) 
WIp (t) ηi,j
(~x) ξk,m,n
(~r ) 
(30)
p=0

I

Because the phonon collision scatterings only consider the Fermi Golden Rule [2] and the spherical
coordinates localize the negative part operator, there is a natural split of the collision operator in
gain and loss terms of probability density rates.
Gain Term of the collisional operator. The gain term, when using the piecewise linear
function (30), becomes
√ Z ∞Z 1 Z π
5
X
X
r
p
p
d~r 0 S(~r 0 , ~r )
χI (~x, ~r 0 )
WIp (t)ηi,j
(~x)ξk,m,n
(~r 0 ) ≈
2 0 −1 −π
p=0
I

5
X √r Z
X
p
p
S(~r 0 , ~r ) χI (~x, ~r 0 )
WIp (t) ηi,j
(~x) ξk,m,n
(~r 0 )d~r 0
2 Kkmn
p=0
I

5 √
XX
r
I

p=0

2

p
χij (~x) WIp (t) ηi,j
(~x)

Z

20

Kkmn

p
S(~r 0 , ~r ) ξk,m,n
(~r 0 ) d~r 0

=

with χI and WIp from (27 - 29)
(
χij (~x) =

1

h
i h
i
if (x, y) ∈ xi− 1 , xi+ 1 × yj− 1 , yj+ 1

0

otherwise

2

2

2

2

.

q
(~r) and an
In the weak formulation, the gain term is multiplied by the test function ηī,q j̄ (~x) ξk̄,
m̄,n̄

integral over the domain ΩI¯ , I¯ = (ī, j̄, k̄, m̄, n̄), with respect to (~x, ~r) is performed, obtaining:
5 √
XX
r

Z
ΩĪ

2

p=0

I

p
χij (~x) WIp (t) ηi,j
(~x)

Z

p
S(~r 0 , ~r ) ξk,m,n
(~r 0 ) d~r 0

Kkmn

q
ηī,q j̄ (~x) ξk̄,
(~r) d~x d~r
m̄,n̄

(31)

or
5
XX
I

WIp (t)

Z
Kk̄m̄n̄

p=0

√ Z
Z x 1Z y 1
j̄+
ī+
r
2
2
p
p
q
(~x) ηī,q j̄ (~x) d~x (32)
χij ηi,j
S(~r 0 , ~r )ξk,m,n
(~
r
)
d~
r
(~r 0 ) d~r 0 ξk̄,
m̄,n̄
2 Kkmn
x 1
y 1
ī−

2

j̄−

2

The integration with respect to x and y gives
Z x 1Z y 1
ī+ 2
j̄+ 2
p
χij (~x) ηi,j
(~x) ηī,q j̄ (~x) dx dy = δiī δj j̄ βpq ∆xī ∆yj̄
xī− 1

2

yj̄− 1
2

where the matrix βpq has the following terms:

1
 1

 1
(βpq ) = 
 1

 0
0

1
1
1
1
0
0

1
1
1
1
0
0

1
1
1
1
0
0

0
0
0
0

0
0
0
0
0

1
3





.




1
3

0

then equation (32) is reduced to:
5
X X

Wīpj̄kmn (t)βpq ∆xī ∆yj̄

k,m,n p=0


√ Z
r
q
p
0
0
0
(~r) d~r.
×
S(~r , ~r ) ξk,m,n (~r ) d~r ξk̄,
m̄,n̄
2 Kkmn
Kk̄m̄n̄
Z

(33)

Loss Term of the collisional operator. The weak formulation of the loss term of the collisional
operator gives
Z

+∞

Z
Φ(t, ~x, ~r)

ΩĪ

dr

Z

0

1

dµ

Φ(t, ~x, ~r) 

X Z

k,m,n

0

0

dϕ S(~r, ~r )
−π

r0 q
q
η (~x) ξk̄,
(~r) d~xd~r ≈
m̄,n̄
2 ī,j̄

(34)

√


ΩĪ

√

π

Z

−1

0

Z

0


r0
q
S(~r, ~r 0 )d~r 0  ηī,q j̄ (~x)ξk̄,
(~r) d~xd~r
m̄,n̄
2

Kkmn

Using the linear approximation of Φ given by (30), integral (34) becomes
Z
ΩĪ

5
X
p=0

"
X

χI (~x, ~r)

5
X

#
p
p
WIp (t) ηi,j
(~x) ξk,m,n
(~r) 

p=0

I

k,m,n


WI¯p (t)

Z
ΩĪ

X Z

p
ηī,p j̄ (~x)ξk̄,
(~r) 
m̄,n̄

X Z

k,m,n

Kkmn

√

r0

2

Kkmn

1√
2


q
r0 S(~r, ~r 0 ) d~r 0  × ηī,q j̄ (~x) ξk̄,
(~r) d~x d~r =
m̄,n̄


q
S(~r, ~r 0 )d~r 0  ηī,q j̄ (~x)ξk̄,
(~r) d~xd~r .
m̄,n̄

21

Therefore, equation (34) reduces to
5
X

WI¯p (t)βpq ∆xī ∆yj̄

p=0

k,m,n

√

"Z

X Z
Kk̄m̄n̄

Kkmn

#
r0
p
q
0
0
S(~r, ~r ) d~r ξk̄,
(~r) ξk̄,
(~r) d~r
m̄,n̄
m̄,n̄
2

(35)

In the case of an energy band function with radial dependance ε(r),
S(r, r0 ) = c0 δ(ε(r0 ) − ε(r)) + c+ δ(ε(r0 ) − ε(r) + αp ) + c− δ(ε(r0 ) − ε(r) − αp )
+1
X

=

cl δ(ε(r0 ) − ε(r) + lαp ) .

(36)

l=−1

The radial energy band function can be projected on the space of piecewise linear functions of r to
obtain

εh (r) =

Nr
X

χk [ε(rk ) + Ak (r − rk )] =

Nr
X

χk [ ε(rk ) + ∂r ε(rk ) (r − rk ) ] ,

(37)

k=1

k=1

and, after this projection, we can calculate the collision integrals involving a delta distribution with
the piecewise linear function in their argument. The computation of such collisional integrals are
given by
Z

Z

√

Z

r q
(~r) d~r
ξ
2 k̄,m̄,n̄
K
K
kmn
k̄
m̄n̄
p=0
k,m,n
#√
"Z
5
X
X Z
r0 0
p
p
q
0
−
WI¯ (t)βpq ∆xī ∆yj̄
S(r, r ) ξk̄,m̄,n̄ (~r) ξk̄,m̄,n̄ (~r) d~r
d~r
2
Kkmn
Kk̄m̄n̄

C(Φh )vh dσ =
K

5
X X

Wīpj̄kmn βpq ∆xī ∆yj̄

p=0

=

S(r

0

p
(~r 0 ) d~r 0
, r) ξk,m,n

k,m,n

X

5
X

Wīpj̄kmn (t)βpq ∆xī ∆yj̄ ×

(38)

k,m,n p=0

Kk̄m̄n̄

−

5
X

+1
X

"Z

Z

Kkmn l=−1

# √
p
cl δ(ε(r) − ε(r0 ) + lαp ) ξk,m,n
(~r 0 ) d~r 0

r q
ξ
(~r) d~r
2 k̄,m̄,n̄

WI¯p (t)βpq ∆xī ∆yj̄ ×

p=0

X Z
k,m,n Kkmn

"Z

+1
X

#√
cl δ(ε(r0 ) − ε(r) +

Kk̄m̄n̄ l=−1

In order to perform the integrations involving

√

p
q
lαp ) ξk̄,
(~r) ξk̄,
(~r)
m̄,n̄
m̄,n̄

d~r

r0 0
d~r .
2

r numerically by means of Gaussian quadrature,

the change of variables r = s2 is applied, so that the functions of s to be integrated are just
22

polynomials,
Z
C(Φh )vh dσ =
K
5
X X

Wīpj̄kmn (t)βpq ∆xī ∆yj̄

Z
Kk̄m̄n̄

k,m,n p=0
+1
X

"

−

#

Z

δ(ε(rk̄ ) + Ak̄ (s − rk̄ ) + lαp − ε(rk ) − Ak (r −

Kkmn

5
X

WI¯p (t)βpq ∆xī ∆yj̄

p=0
+1
X

"

=

p
rk )) ξk,m,n
(~r 0 ) d~r 0

dsdµdϕ

s0 0
2s ×
Kkmn 2

X Z
k,m,n

#

Z

02

δ(ε(s ) + lαp − ε(rk̄ ) − Ak̄ (r − rk̄ ))

cl
Kk̄m̄n̄

l=−1

0

2

cl

l=−1

s q
ξ
(s2 , µ, ϕ) 2s ×
2 k̄,m̄,n̄

5
X X

q
p
(~r)
(~r) ξk̄,
ξk̄,
m̄,n̄
m̄,n̄

d~r ds0 dµ0 dϕ0

Wīpj̄kmn (t)βpq ∆xī ∆yj̄ ×

(39)

k,m,n p=0

Z
Kk̄m̄n̄



+1
X



−

µm+ 1

Z

2

cl

µm− 1

l=−1
5
X


ε(rk̄ ) + Ak̄ (s2 − rk̄ ) + lαp − ε(rk ) + Ak rk
χk
×
Ak



Z ϕ 1
2
n+ 2
ε(rk̄ ) + Ak̄ (s − rk̄ ) + lαp − ε(rk ) + Ak rk 0 0
p
, µ , ϕ dµ0 dϕ0  s2 dsdµdϕ
ξk,m,n
Ak
ϕ 1


q
ξk̄,
(s2 , µ, ϕ)
m̄,n̄

n− 2

2

WI¯p (t)βpq ∆xī ∆yj̄

p=0



+1
X



Z
cl

2

µm̄− 1

l=−1

with r(s) =

µm̄+ 1
2

Z

ϕn̄+ 1

2

ϕn̄− 1

X Z
k,m,n Kkmn

0

0 02

dµ dϕ s χk̄



ε(rk ) + Ak (s02 − rk ) + lαp − ε(rk̄ ) + Ak̄ rk̄
Ak̄


q
p
(r(s), µ, ϕ) dµdϕ
(r(s), µ, ϕ) ξk̄,
ξk̄,
m̄,n̄
m̄,n̄


×

ds0 ,
r(s)

2

ε(rk )+Ak (s02 −rk )+lαp −ε(rk̄ )+Ak̄ rk̄
.
Ak̄

The integrals above involve only polynomials, which are numerically computed by Gaussian quadrature rules.

5.4

The algorithm for time evolution

Starting with given initial and boundary conditions, the algorithm advances from tn to tn+1 by the
following way:
1. Compute the density ρ.
2. Solve the Poisson equation and find the electric field E.

23

3. Compute the transport terms ai ’s.
4. Compute the collision part.
5. Solve the (large) system of ordinary differential equations for the coefficients of the linear
approximation of Φh (which are obtained from the DG formulation), by using a TVD Runge
- Kutta scheme.
6. Repeat the previous steps as needed.
Algorithm 1: DG-BP Algorithm
Data: Given Initial and Boundary Conditions
Result: Time Evolution of Probability Density Function
initialization;
while tl < Tmax do
Compute density ρ(x, tl );
Solve Poisson Eq. to find the electric field E(x, tl );
tl+1 = tl + ∆tl+1 ;
Compute the collision terms;
Compute the transport terms ai (x, k, tl+1 )’s;
Solve the large system of ODEs for the coefficients of Φh |tl+1 (obtained from the DG
formulation) by a TVD Runge-Kutta scheme. Repeat the previous steps if needed.
end

6

The n+ -n-n+ silicon diode

We consider the symmetric case of a 1D n+ -n-n+ diode, in which the conduction band energy
function is assumed to be of the form ε(|k|) = ε(r). This assumption preserves azimuthal symmetry
for the problem if the initial condition is independent of the azimuthal direction ϕ. Therefore, under
these assumptions the problem has azimuthal symmetry in k for all times t ≥ 0, so it suffices to
consider k = k(r, µ), reducing then the dimensionality of the problem to 1-D in x-space and 2-D in
k = k(r, µ), then the problem reduces to a 3-D plus time. Assuming E has null y and z components,
this symmetric case reduces the TBE to
∂Φ
∂
∂
∂
+
(a1 Φ) +
(a4 Φ) +
(a5 Φ) = C(Φ),
∂t
∂x
∂r
∂µ
where the terms a1 , a4 and a5 are now simplified.

24

(40)

The Poisson equation is reduced to
∂
∂x



∂Ψ
r
= cp [ρ(t, x) − ND (x)] .
∂x

(41)

For this case both the potential and electric field have analytic integral solutions, that are easily
computed numerically for the piecewise linear approximation of the density ρh . Then, such electric
field solution is projected in the Vh1 space of piecewise linear polynomials.

6.1

Device specifics

We consider first a diode of 1 µm length, with an n-channel of 400 nm length, doping of 5×1023 m−3
in the n+ region and 2 × 1021 m−3 in the n region. We also consider a 0.25µm diode with a 50 nm
channel with n+ -doping of 5 × 1024 m−3 , and n-doping of 1 × 1021 m−3 .

6.2

Numerical simulations

The space Vh1 of piecewise linear polynomials in (x, r, µ), with time dependent coefficients, is used
as both the trial and test space in our DG scheme. The input data of the numerical simulations is
• Computational domain: x ∈ [0, 1], r ∈ [0, rmax ], µ ∈ [−1, 1], where rmax is taken in the
numerical experiments such that Φ(t, x, r, µ) ≈ 0 for r ≈ rmax (for example, rmax = 36 for
Vbias = 0.5 Volts in the 400nm channel case)

• Initial condition: Φ(0, x, r, µ) = Πh


√
r −ε(r)
e
, where C constant is such that
CND (x)
2

ρ(x, 0) equals the doping ND (x) at t = 0. The initial condition is projected by Πh into Vh1 .
• Boundary conditions: Neutral charges at the endpoints x 1 = 0 and xNx + 1 = 1.
2

Φ(t, 0, r, µ) = ND (0)

2

Φ(t, x1 , r, µ)
Φ(t, xNx , r, µ)
and Φ(t, 1, r, µ) = ND (1)
.
ρ(t, x1 )
ρ(t, xNx )

Cut-off in the k-space Φ(t, x, rmax , µ) = 0.
• Applied potential - bias: Ψ(t, 0) = 0 and Ψ(t, 1) = V0 .
No boundary conditions are needed on r = 0, µ = ±1. Upwind fluxes in r and µ are analytically
zero at these boundaries, since they are related to points in k-space such as the origin and the
25

poles, which are transformed into boundaries when applying the spherical change of coordinates.
It is very simple to verify that a4 = 0 at r = 0, and a5 = 0 at µ = ±1.

7

Numerical results

The BP transport along the EPM spherical average energy band is numerically simulated by means
of our DG-BP solver, and compared to simulations where the values related to the analytical
Parabolic and Kane band models are implemented numerically. We compare simulations for two
n+ -n-n+ silicon diodes with different characteristics. The first one has a length of 1µm, an nchannel length of 400nm, n+ doping of 5 · 1023 m−3 , and n doping of 2 · 1021 m−3 . The other one
has a device length of 0.25µm, an n-channel length of 50nm, an n+ doping of 5 · 1024 m−3 , and n
doping of 1 · 1021 m−3 . We show simulations for a potential bias of V0 = 0.5 V . For the 400 nm
channel diode, the number of cells used in the simulations for each of the variables was: Nx = 120,
Nr = 80 and Nµ = 24. The interval size for r is taken as ∆r = 0.45, having then rmax = 36.
We use a mesh as in [19] which gives better resolution close to the first juncture at x = 0.3µm,
and which also has a finer refinement close to the pole in the direction of the electric field. It uses
∆x = 0.01 for the first 20 cells in x-space, ∆x = 0.005 for the next 40 cells, and ∆x = 0.01 for the
last 60 cells. Regarding µ, it uses 12 cells for µ ∈ [−1, 0.7], and 12 cells for µ ∈ [0.7, 1].
For the 50 nm channel diode, the number of cells used in the simulations for each of the variables
was: Nx = 64, Nr = 80 and Nµ = 20. The interval size for r is taken as ∆r = 0.8, having then
rmax = 64. As in [19], we use a mesh intended to give better resolution close to the junctures at
x = 0.1µm and x = 0.15µm, and which also has a finer refinement close to the pole in the direction
of the electric field. It uses ∆x = 0.01 for the first 9 cells in x-space, ∆x = 0.001 for the next 20 cells
close to the the first juncture at x = 0.1µm, ∆x = 0.005 for 6 cells at the center of the n-channel,
∆x = 0.001 for the next 20 cells close to the second juncture at x = 0.15µm, and ∆x = 0.01 for
the last 9 cells. Regarding µ, it uses 10 cells for µ ∈ [−1, 0.7], and 10 cells for µ ∈ [0.7, 1]. We let
the solver run until t = 5.0ps, a time when the simulations are close to a numerical stationary state.

26

Density(log-scale) vs. Position (x)

Density(log-scale) vs. Position (x)
1e+19
EPMravg
Kane
Parab

Density (1/cm3)

Density (1/cm3)

1e+17

1e+18

1e+16
EPMravg
Kane
Parab
1e+17
0

0.1

0.2

0.3

0.4
0.5
0.6
x (micro-meters)

0.7

0.8

0.9

1

0

0.05

0.1
0.15
x (micro-meters)

0.2

0.25

Figure 4: Density (ρ, in log-scale) vs. position (x) plots for different conduction band models:
parabolic, Kane, EPM average. t = 10.0ps. Left:400nm channel. Right: 50nm channel. 0.5 Volts
Bias.
1.2

2
EPMravg
Kane
Parab

EPMravg
Kane
Parab

1.8

1
1.6
1.4

<v> (107 cm/s)

<v> (107 cm/s)

0.8

0.6

1.2
1
0.8

0.4
0.6
0.4
0.2
0.2
0

0
0

0.1

0.2

0.3

0.4
0.5
0.6
x (micro-meters)

0.7

0.8

0.9

1

0

0.05

0.1
0.15
x (micro-meters)

0.2

0.25

Figure 5: Average velocity (v) vs. position (x) plots for different conduction band models:
parabolic, Kane, EPM average. t = 10.0ps. Left:400nm channel. Right: 50nm channel. 0.5
Volts Bias
We show plots of the average velocity, the average energy, the momentum (proportional to the
current), the electric field and potential, for both the 400nm channel and 50nm channel diodes.
There is a clear quantitative difference, particularly in kinetic moments such as average velocity,
average energy, and momentum (current), whose values depend on the energy band model used
in each case then. This should be expected since these kinetic moments are averages of quantities
related to ε(k) or its partial derivatives in k-space.

27

0.12

0.25
EPMravg
Kane
Parab

0.11

EPMravg
Kane
Parab
0.2

0.1

0.09

<e> (eV)

<e> (eV)

0.15
0.08

0.07

0.1
0.06

0.05

0.05

0.04

0.03

0
0

0.1

0.2

0.3

0.4
0.5
0.6
x (micro-meters)

0.7

0.8

0.9

1

0

0.05

0.1
0.15
x (micro-meters)

0.2

0.25

Figure 6: Average energy (ε) vs. position (x) for different conduction band models: parabolic,
Kane, EPM average. t = 10.0ps. Left:400nm channel. Right: 50nm channel. 0.5 Volts Bias.
6.5e+22

5e+24
EPMravg
Kane
Parab

EPMravg
Kane
Parab

4.5e+24

U ( 1/(s*cm2) )

U ( 1/(s*cm2) )

6e+22

5.5e+22

5e+22

4e+24

3.5e+24

4.5e+22

3e+24
0

0.1

0.2

0.3

0.4
0.5
0.6
x (micro-meters)

0.7

0.8

0.9

1

0

0.05

0.1
0.15
x (micro-meters)

0.2

0.25

Figure 7: Current (Momentum) vs. position (x) for different conduction band models: parabolic,
Kane, EPM average. t = 10.0ps. Left:400nm channel. Right: 50nm channel. 0.5 Volts Bias.
40

100
EPMravg
Kane
Parab

30

EPMravg
Kane
Parab
50

20
0

E (103 Volts/cm)

E (103 Volts/cm)

10

0

-10

-50

-100

-20
-150
-30
-200
-40

-50

-250
0

0.1

0.2

0.3

0.4
0.5
0.6
x (micro-meters)

0.7

0.8

0.9

1

0

0.05

0.1
0.15
x (micro-meters)

0.2

0.25

Figure 8: Electric field (E) vs. position (x) plots for different conduction band models: parabolic,
Kane, EPM average. t = 10.0ps. Left:400nm channel. Right: 50nm channel. 0.5 Volts Bias.
28

0.6

0.6
EPMravg
Kane
Parab

0.5

0.5

0.4

0.4

0.3

0.3

V (Volts)

V (Volts)

EPMravg
Kane
Parab

0.2

0.2

0.1

0.1

0

0

-0.1

-0.1
0

0.1

0.2

0.3

0.4
0.5
0.6
x (micro-meters)

0.7

0.8

0.9

1

0

0.05

0.1
0.15
x (micro-meters)

0.2

0.25

Figure 9: Electric Potential (V ) vs. position (x) plots for different conduction band models:
parabolic, Kane, EPM average. t = 10.0ps. Left:400nm channel. Right: 50nm channel. 0.5 Volts
Bias.

7.1

Current - Voltage Characteristic Curves

We also perform a study of the current - voltage characteristic curves for the considered pair of
devices with our 3 different band models. We show in Fig. 13 the plots of the Momentum (U)
vs. Bias (V ) for voltages of 0, 0.25, 0.5, 0.75, 1.0 V with the mesh described in the previous section
for the 400nm channel device, and in Fig. 15 the respective plots for the 50nm channel device.
We observe a clear quantitative difference between the currents according to the used band model.
The currents predicted using the parabolic band are bigger than the ones obtained using the Kane
band, and the EPM radial average predicts lower current values than the other two models. We
compute as well the IV curves with a refined mesh in the r variable subdividing the intervals of
the original mesh in half the size, and with a coarse mesh as well by joining pairs of subsequent
intervals, using in all cases the same domain. The current voltage characteristics predicted with
the coarser and finer r-meshes are quite close to the predictions with the original mesh used.
We compare as well in figures 14 and 16 the momentum predicted with a finer r-mesh for our set of
bias with the differences of these moments with the predicted ones by our original mesh (U − UR )
and a coarser mesh (UC − UR ). For the 400nm channel, we observe in Fig. 14 that the order of
magnitude of the current values for nonzero voltages is twice bigger than the difference between
those currents for the coarser and the finer mesh, and three times bigger than the difference between

29

2e+24
1.5e+24
1e+24
(1/scm^2) 5e+23
0
-5e+23
-1e+24
-1.5e+24

0

1

2

3

4

5

t (pico-seconds)

6

7

8

9

10 0

1
0.9
0.8
0.7
0.6
0.5
0.4
x (micro-meters)
0.3
0.2
0.1

Figure 10: Current (Momentum) vs. (t, x) for the 400nm channel diode, using the EPM average
band. From initial time to t = 10.0ps. 0.5 Volts bias. The initial condition is proportional to the
Maxwellian exp(−ε(r)) times the doping profile N (x). The initial oscillations are produced by the
initial state being far from the final steady state.

30

#10 -4

f (x o ; kx, ky)

2.5
2
1.5
1
0.5
0
6
10

4

ky

5
0

2
-5
0

-10

kx

Figure 11: PDF f (r, µ; t0 , x0 ) vs. (kx , ky ) coordinates (Azimuthal symmetry, with kz = 0) at the
point x0 = 0.3µm at t0 = 10.0ps, for the 1µm diode with a 400nm channel, using the EPM radial
average energy band. 0.5 Volts bias.
31

#10 -3

f (x o ; kx, ky)

1
0.8
0.6
0.4
0.2
0
6
10

4

ky

5
0

2
-5
0

-10

kx

Figure 12: PDF f (r, µ; t0 , x0 ) vs. (kx , ky ) coordinates (Azimuthal symmetry, with kz = 0) at the
point x0 = 0.7µm at t0 = 10.0ps, for the 1µm diode with a 400nm channel, using the EPM radial
average energy band. 0.5 Volts bias.
32

those currents for the original and final mesh. The differences are consistently lower for the EPM
average band when compared to the Kane and Parabolic bands.
For the 50nm channel, we observe in Fig. 16 that the order of magnitude of the current values
for nonzero voltages is twice bigger than the difference between those currents for the coarser and
the finer mesh, as for the difference between those currents for the original and final mesh. The
differences are consistently lower for the EPM average band when compared to the Kane and
Parabolic Band.
The set of plots presented for our devices point to several facts. In addition to the clear quantitative difference between the current voltage characteristics predicted according to which band
model is used, with the EPM radial average currents below the Kane predictions and these below
the Parabolic currents, the fact that the difference in the IV curves obtained with different meshes
is at least two orders of magnitude below the current values indicates that IV curves of comparable quantitative accuracy can be obtained with coarser meshes in the r-variable with a reduced
computational effort. We also notice that, for nonzero bias, the values of the differences U − UF ,
UC − UF between predicted currents is always positive, that is, the coarsening of the mesh comes
with a slight increase of the predicted current. Except for one point in Fig. 16 related to the Kane
band with a 1.0V bias, the differences UC − UF in currents with the coarser mesh are bigger than
the difference U − UF between the current predicted with the original mesh and the refined mesh.
Moreover, since we observe that both these differences UC − UF and U − UF , are lower for the EPM
average band model, the variations in the current prediction are smaller under the coarsening of
the mesh when using the EPM average as a band model.

33

Momentum (U) vs. Potential (V)

Momentum (U) vs. Potential (V)

Momentum (U) vs. Potential (V)

1e+23

1e+23

1e+23

8e+22

8e+22

8e+22

2

6e+22

4e+22

U ( 1/(s*cm2) )

1.2e+23

U ( 1/(s*cm2) )

1.2e+23

U ( 1/(s*cm ) )

1.2e+23

6e+22

4e+22

2e+22

4e+22

2e+22

2e+22

EPMravg
KANE0p5
PARAB

EPMravg
KANE0p5
PARAB

0

EPMravg
KANE0p5
PARAB

0
0

0.2

0.4

0.6

0.8

6e+22

1

0
0

0.2

V (Volts)

0.4

0.6

0.8

1

0

0.2

V (Volts)

0.4

0.6

0.8

1

V (Volts)

Figure 13: Momentum (U) at t = 10.0ps vs. Applied Bias (V ) plots for voltages of 0., 0.25, 0.5,
0.75, 1.0 V for the parabolic, Kane, and EPM average band models. Center: Original Mesh. Left:
Refined r-Mesh by a factor of 1/2. Right: Coarsened r-Mesh by a factor of 2. 400nm channel
device. There is a clear quantitative difference trend between the current values for the different
bands, with the parabolic band giving current values bigger than the Kane band, and the EPM
average band having lower current values than the other models.
Current (I) vs. Applied Potential (V)
1e+24

I ( 1/(s*cm2) )

1e+23

finerEPMravg
finerKANE0p5
finerPARAB
diffEPMravg
diffKane
diffParab
cdiffEPMravg
cdiffKane
cdiffParab

1e+22

1e+21

1e+20

1e+19
0

0.2

0.4

0.6

0.8

1

V (Volts)

Figure 14: Log plot of Momentum (U) at t = 10.0ps vs. Applied Bias (V ) for voltages of 0., 0.25,
0.5, 0.75, 1.0 V for the parabolic, Kane, and EPM average band models with a refined r-Mesh by
a factor of 1/2. Log plots of the difference of these values with the original mesh (U − UR ) and
with a coarser mesh (UC − UR ). 400nm channel device. Notice that the order of magnitude of the
current values for nonzero voltages is twice bigger than the difference between those currents for
the coarser and the finer mesh, and three times bigger than the difference between those currents
for the original and final mesh. The differences are consistently lower for the EPM average band
when compared to the Kane and Parabolic bands.
34

Momentum (U) vs. Potential (V)

Momentum (U) vs. Potential (V)

Momentum (U) vs. Potential (V)

9e+24

9e+24

1e+25

8e+24

8e+24

9e+24

7e+24

7e+24

6e+24

6e+24

5e+24

5e+24

8e+24
7e+24

3e+24

U ( 1/(s*cm2) )

U ( 1/(s*cm2) )

2

U ( 1/(s*cm ) )

6e+24

4e+24

4e+24
3e+24

5e+24
4e+24
3e+24

2e+24

2e+24

1e+24

1e+24

0

2e+24
1e+24

0

EPMravg
KANE0p5
PARAB

-1e+24

EPMravg
KANE0p5
PARAB

0

0.2

0.4

0.6

0.8

1

EPMravg
KANE0p5
PARAB

0

-1e+24

-1e+24
0

0.2

V (Volts)

0.4

0.6

0.8

1

0

0.2

V (Volts)

0.4

0.6

0.8

1

V (Volts)

Figure 15: Momentum (U) at t = 10.0ps vs. Applied Bias (V ) plots for voltages of 0., 0.25, 0.5,
0.75, 1.0 V for different conduction band models: parabolic, Kane, EPM average. Center: Original
Mesh. Left: Refined r-Mesh by a factor of 1/2. Right: Coarsened r-Mesh by a factor of 2. 50nm
channel device. There is a clear quantitative difference trend between the current values for the
different bands, with the parabolic band giving current values bigger than the Kane band, and the
EPM average band having lower current values than the other models.
Current (I) vs. Applied Potential (V)
1e+25

1e+24

I ( 1/(s*cm2) )

1e+23

1e+22

1e+21
finerEPMravg
finerKANE0p5
finerPARAB
diffEPMravg
diffKane
diffParab
cdiffEPMravg
cdiffKane
cdiffParab

1e+20

1e+19
0

0.2

0.4

0.6

0.8

1

V (Volts)

Figure 16: Log plot of Momentum (U) at t = 10.0ps vs. Applied Bias (V ) for voltages of 0., 0.25,
0.5, 0.75, 1.0 V for the parabolic, Kane, and EPM average band models with a refined r-Mesh by
a factor of 1/2. Log plots of the difference of these values with the original mesh (|U − UR |) and
with a coarser mesh (|UC − UR |) (these values are only negative for 0 V, therefore the absolute
value). 50nm channel device. Notice that the order of magnitude of the current values for nonzero
voltages is twice bigger than the difference between those currents for the coarser and the finer
mesh, as for the difference between those currents for the original and final mesh. The differences
are consistently lower for the EPM average band when compared to the Kane and Parabolic Band.
35

8

Conclusions

The implementation of a spherical average of an EPM full band structure as a conduction energy
band model in a DG solver for Boltzmann - Poisson represents a computational strategy that is
a midpoint between a radial and an anisotropic full band model. This difference in the values of
the energy band and its derivatives, introduced via the spherical average of EPM bandstructure
values over k-spheres and the spline interpolation of derivatives predicts a quantitative correction
in kinetic moments (averages) related to the energy band model, such as average velocity, energy,
and particularly the momentum (proportional to the current) given by our solver. It emphasizes
then the importance of an accurate physical modeling of the energy band structure and its partial
derivatives as these functions drive the mechanisms of collision (electron - phonon scattering) and
transport (via the electron group velocity) whose balance is the core modeling of electron transport
in semiconductor by the Boltzmann - Poisson system. This highlights the importance of band
models and their features (anisotropy, numerical approximation of its values and interpolation of
their derivatives, for example) in the BP numerical modeling of electron transport via DG schemes.
Future work will focus on the implementation of anisotropic EPM full bands in the DG solver for
BP, and on developing a positivity-preserving DG numerical scheme, along with error estimates for
the BP system under consideration.

9

Appendix

A

TBE in Divergence Form for k in Spherical Coordinates

The divergence in k in the standard spherical coordinates (|k|, θ, ϕ) for a vector field A(x, k, t) =
(A1 , A2 , A3 ) = A|k| ê|k| + Aθ êθ + Aϕ êϕ has the expression

∇k · A =

∂Aϕ
1 ∂(|k|2 A|k| )
1
∂ (Aθ sin θ)
1
+
+
2
|k|
∂|k|
|k| sin θ
∂θ
|k| sin θ ∂ϕ

(42)

The divergence of A, with a respective orthogonal decomposition: A = Ar êr +Aµ êµ +Aϕ êϕ , in terms
of the modified spherical coordinates (r, µ, ϕ) used in the TBE is obtained from (42) by taking into

36

p
√ dµ
dr
~
= √ ∗
= − 1 − µ2 ,
2 r,
d|k|
dθ
2m kB TL
following that êµ = −êθ , Aµ = −Aθ , and ê|k| = êr , A|k| = Ar . We have then


∂Aϕ
1 ∂(rAr ) dr
~
1
1
∂(Aθ sin θ) dµ
√
∇k · A =
+√ ∗
+√
r ∂r d|k|
∂µ
dθ
r sin θ
r sin θ ∂ϕ
2m kB TL
!#
"
!
p
~
1
2
∂
∂
1 − µ2
∂
p
.
= √ ∗
−
Aϕ
·√
(rAr ) +
Aθ +
∂µ
2
∂ϕ 2 1 − µ2
r ∂r
2m kB TL
account: |k|2 = (2m∗ kB TL /~2 )r, µ = cos θ. Therefore

We obtain then the divergence in the modified spherical coordinates used in this work:
!#
!
p
"

√
∂
r
1 − µ2
∂
1
~
∂
p
Aϕ
∇k · A = √ ∗
(rAr ) +
Aµ +
2
∂r
∂µ
2
∂ϕ 2 1 − µ2
2m kB TL

(43)

So, the k-transport term in the TBE (11) can be expressed in the divergence form (43) by using:
 √

 √
√ 
−q E(t, x)
r −q E(t, x)
r
r
· ∇k f =
∇k ·
f =
∇k · A
(44)
2
~
2
~
2
for the vector field
A=

−q E(t, x)
−q E(t, x)
Φ
f=
·√
~
~
r/2

(45)

The formula (43) mentioned above for A can be interpreted geometrically as a flow of electric field
in the orthogonal directions of the spherical coordinate geometry used, since by definition:
Ar = A · êr ,

B

Aµ = A · êµ ,

Aϕ = A · êϕ .

(46)

Details of calculation of the transport terms

The following notation for boundary terms will be needed. We denote by

Φ̂i± 1

= Φ̂(t, xi± 1 , y, r, µ, ϕ) ,

Φ̂j± 1 = Φ̂(t, x, yj± 1 , r, µ, ϕ)

Φ̂k± 1

= Φ̂(t, x, y, rk± 1 , µ, ϕ) ,

Φ̂m± 1 = Φ̂(t, x, y, r, µm± 1 , ϕ) ,

2

2

2

2

2

2

2

p
ηi±1,j
|i± 1

2

p
= ηi±1,j
(xi± 1 , y) ,

p
p
ηi,j
|i± 1 = ηi,j
(xi± 1 , y) ,

p
ηi,j±1
|j± 1

p
= ηi,j±1
(x, yj± 1 ) ,

p
p
ηi,j
|j± 1 = ηi,j
(x, yj± 1 ) ,

p
ξk,m,n
|k± 1

p
= ξk,m,n
(rk± 1 , µ, ϕ) ,

2

2

2

2

2

p
ξk±1,m,n
|k± 1
2

=

p
ξk,m,n±1
|n± 1

p
= ξk,m,n±1
(r, µ, ϕn± 1 ) ,

2

p
ξk±1,m,n
(rk± 1 , µ, ϕ) ,
2

Φ̂n± 1 = Φ̂(t, x, y, r, µ, ϕn± 1 )
2

2

(48)

2

2

2

(47)

2

p ∈ {0, 1, · · · , 5}.

2

p
p
ξk,m,n
|m± 1 = ξk,m,n
(r, µm± 1 , ϕ) ,
2

p
ξk,m±1,n
|m± 1
2

2

=

p
ξk,m±1,n
(r, µm± 1 , ϕ) ,
2

p ∈ {0, 1, · · · , 5}.

2

37

p
p
ξk,m,n
|n± 1 = ξk,m,n
(r, µ, ϕn± 1 )
2

2

(49)

We consider first the weak formulation of the transport terms in space, namely
∂
∂y

∂
∂x

(a1 Φ) and

(a2 Φ), related to the advection in x, where the first cartesian component a1 of the electron

group velocity is involved. Because their discretization forms are similar we only present the one
for

∂
∂x

(a1 Φ).
Z

(A1 )

=
=

∂
[a1 (r, µ, ϕ) Φ(t, ~x, ~r)] ψ(~x, ~r) d~x d~r =
∂x
ΩI
Z
Z
i
h
1
∂
a1 (~r) Φ̂i+ 12 ψi+ 12 − Φ̂i− 12 ψi− 21 d~xd~r −
a1 (~r)Φ(t, ~x, ~r) ψ(~x, ~r)d~xd~r .
∆xi ΩI
∂x
ΩI

Due to the upwind flux rule, we have to consider two cases depending on the sign of a1 . In the
sequel, the symbol ≈ will denote the approximation of given integral terms
If a1 (~r) > 0 in Kkmn , for q = 0, ..., 5, one obtains
(A1 ) ≈

Z
5
X
W p (t)
I

∆xi

p=0

−

∆xi

5
X

WIp (t)

p=0

=

5 Z
X
p=0

−

5
X

yj+ 1
2

yj− 1

ΩI

Z
ΩI

p
p
q
q
a1 (~r)ηi−1,j
|i− 21 ξk,m,n
(~r)ηi,j
|i− 12 ξk,m,n
(~r)d~xd~r

p
p
a1 (~r) ηi,j
(~x) ξk,m,n
(~r)

p
q
ηi,j
|i+ 21 ηi,j
|i+ 12 dy

Z
·

2

5 Z
X
p=0

−

p
p
q
q
a1 (~r)ηi,j
|i+ 12 ξk,m,n
(~r)ηi,j
(~r)d~xd~r
|i+ 12 ξk,m,n

Z
p
5
X
Wi−1
jkmn (t)
p=0

−

ΩI

yj+ 1
2

yj− 1

q
p
ηi−1,j
|i− 21 dy
|i− 21 ηi,j

Kkmn

 2 δq4
∆xi
p=0

Z

xi+ 1

2

xi− 1

2

Z

yj+ 1
2

p
q
a1 (~r)ξk,m,n
(~r)ξk,m,n
(~r)d~r · WIp (t)

Z

2



2 δq4 q
ξ
(~r) d~x d~r
∆xi k,m,n

Kkmn


p
ηi,j
(~x)

p
q
p
a1 (~r)ξk,m,n
(~r)ξk,m,n
(~r)d~r · Wi−1
jkmn (t)

Z

d~x

yj− 1

Kkmn

2

38

p
q
a1 (~r)ξk,m,n
(~r)ξk,m,n
(~r)



d~r · WIp (t) .

If a1 (~r) < 0 in Kkmn , for q = 0, ..., 5,
(A1 ) ≈

Z
p
5
X
Wi+1
jkmn (t)
∆xi

p=0

−

−

1
∆xi
5
X

5
X

WIp (t)

WIp (t)

Z
ΩI

5 Z y 1
X
j+

−

−

yj− 1

p
q
ηi,j
|i− 12 ηi,j
|i− 21 dy

2

5
X

p
p
(~x) ξk,m,n
(~r)
a1 (~r) ηi,j

2

5 Z y 1
X
j+
p=0

p
q
q
p
|i− 21 ξk,m,n
(~r)ηi,j
(~r) d~xd~r
|i− 12 ξk,m,n
a1 (~r)ηi,j

p
q
ηi+1,j
|i+ 12 ηi,j
|i+ 12 dy

2

p=0

ΩI

p=0

p=0

=

ΩI

Z

p
q
q
p
(~r)ηi,j
|i+ 12 ξk,m,n
(~r) d~xd~r
a1 (~r)ηi+1,j
|i+ 21 ξk,m,n

yj− 1

2

2 δq4
∆xi
p=0

Z

xi+ 1

2

xi− 1
2

Z

yj+ 1
2

2 δq4 q
ξ
(~r) d~x d~r
∆xi k,m,n

Z

p
q
p
a1 (~r)ξk,m,n
(~r)ξk,m,n
(~r)d~r Wi+1
jkmn (t)

Kkmn

Z
Kkmn

p
q
a1 (~r)ξk,m,n
(~r)ξk,m,n
(~r)d~r WIp (t)

p
ηi,j
(~x) dx dy

yj− 1

Z
Kkmn

2

p
q
a1 (~r) ξk,m,n
(~r) ξk,m,n
(~r) d~r WIp (t) .

We consider now the weak formulation for the transport terms in momentum space
∂
∂µ

(a5 Φ) +

∂
∂ϕ

∂
∂r

(a4 Φ) +

(a6 Φ) advected by the electric field. It can be noticed in Eq. 11 that all the terms

a4 , a5 , a6 are the sum of terms of the form a∗ (~r) E∗ (t, ~x), where E∗ (t, ~x) is a cartesian component
of the electric field.
The r-derivative including a4 can be split as a sum of terms as the following
Z

∂
[a∗ (~r) E∗ (t, ~x) Φ(t, ~x, ~r)] ψ(~x, ~r) d~x d~r
∂r
ΩI
Z
Z
1
1
=
a∗ |k+ 1 E∗ (t, ~x) Φ̂k+ 1 ψk+ 1 d~x d~r −
a∗ |k− 1 E∗ (t, ~x) Φ̂k− 1 ψk− 1 d~x d~r
2
2
2
2
2
2
∆rk ΩI
∆rk ΩI
Z
∂
−
a∗ (~r) E∗ (t, ~x) Φ(t, ~x, ~r) ψ(~x, ~r) d~x d~r .
∂r
ΩI

(A4∗ ) =

39

By the Upwind Flux Rule, if a∗ (~r) E∗ (t, ~x) > 0 in ∂r± ΩI , for q = 0, ..., 5 ,
(A4∗ ) ≈

−

Z
5
X
WIp
p
p
q
q
(~x)ξk,m,n
|k+ 12 ηi,j
(~x)ξk,m,n
|k+ 12 d~xd~r
a∗ |k+ 21 E∗ (t, ~x)ηi,j
∆r
k
ΩI
p=0
Z
p
5
X
Wijk−1mn
p=0
5
X

∆rk

p
q
p
q
(~x)ξk−1,m,n
|k− 21 ηi,j
(~x)ξk,m,n
|k− 12 d~xd~r
a∗ |k− 21 E∗ ηi,j

ΩI

Z

2 δq1 q
η (~x) d~x d~r
∆rk i,j
ΩI
p=0


5 Z x 1Z y 1
 Z µm+ 1Z ϕn+ 1
X
j+
i+
2
2
2
2
p
q

=
E∗ (t, ~x)ηi,j
(~x)ηi,j
(~x)dxdy  
a | 1 ξp
| 1 ξq
| 1 dµdϕ
 µ 1 ϕ 1 ∗ k+ 2 k,m,n k+ 2 k,m,n k+ 2
x
y
1
1
p=0
i−
j−
m−
n−
2
2
2
2




Z
Z µ 1Z ϕ 1

n+
m+
2 δq1
2
2
p
p
q
p
a∗ |k− 12 ξk−1,m,n
a∗ (~r) ξk,m,n
(~r) d~r WIp − 
|k− 12 ξk,m,n
|k− 21 dµ dϕ Wijk−1mn
.
−

∆rk Kkmn
µ
ϕ
1
1
−

WIp (t)

p
p
(~x) ξk,m,n
(~r)
a∗ (~r) E∗ (t, ~x) ηi,j

m−

2

n−

2

If a∗ (~r) E∗ (t, ~x) < 0 in ∂r± ΩI , q = 0, ..., 5 ,
(A4∗ ) ≈

−

−

Z
5
1 X p
p
p
q
q
Wijk+1mn (t)
a∗ |k+ 21 E∗ (t, ~x) ηi,j
(~x) ξk+1,m,n
|k+ 12 ηi,j
(~x) ξk,m,n
|k+ 12 d~x d~r
∆rk p=0
ΩI
Z
5
X
WIp
p
q
p
q
a∗ |k− 12 E∗ ηi,j
(~x)ξk,m,n
(~x)ξk,m,n
|k− 21 ηi,j
|k− 12 d~xd~r
∆r
k
Ω
I
p=0
5
X

WIp (t)

Z

p=0


Z
5
X

=
p=0

xi+ 1

2

xi− 1

ΩI

Z

p
p
a∗ (~r) E∗ (t, ~x) ηi,j
(~x) ξk,m,n
(~r)

yj+ 1
2

2 δq1 q
η (~x) d~x d~r
∆rk i,j


p
q
E∗ (t, ~x) ηi,j
(~x) ηi,j
(~x) dx dy 

yj− 1

2
2


Z
Z

ϕn+ 1
µm+ 1
2
2
p
q
p
× 
a∗ |k+ 12 ξk+1,m,n
|k+ 12 ξk,m,n
|k+ 12 dµ dϕ Wijk+1mn
(t)
 µ 1 ϕ 1
m−
n−
2
2



Z µ 1Z ϕ 1
Z
m+
n+
2 δq1
2
2
p
q
p
p

−
a∗ |k− 12 ξk,m,n |k− 12 ξk,m,n |k− 21 dµ dϕ +
a∗ (~r) ξk,m,n (~r) d~r WI (t) .
∆rk Kkmn
µ
ϕ
1
1
m−

2

n−

2

The weak form for the term related to

Z
(A5∗ ) =
ΩI

∂
∂µ

(a5 Φ) is

∂
[a∗ (~r) E∗ (t, ~x) Φ(t, ~x, ~r)] ψ(~x, ~r) d~x d~r .
∂µ

By the Upwind Flux rule, if a∗ (~r) E∗ (t, ~x) > 0 in ∂µ± ΩI , q = 0, ..., 5 ,

40

(A5∗ ) ≈


Z
5
X


xi+ 1

2

xi− 1

p=0

yj+ 1

Z

2


q
p
(~x) dx dy 
(~x) ηi,j
E∗ (t, ~x) ηi,j

yj− 1

2
 2

Z
 Z rk+ 1 Z ϕn+ 1
2 δq2
2
2
p
q
p
× 
a∗ |m+ 12 ξk,m,n
|m+ 21 ξk,m,n
|m+ 12 dr dϕ −
a∗ (~r) ξk,m,n
(~r) d~r WIp (t)
 r 1 ϕ 1
∆µm Kkmn
k−
n−
2
2



Z r 1Z ϕ 1

n+
k+
2
2
p
q
p
a∗ |m− 12 ξk,m−1,n
|m− 12 ξk,m,n
|m− 12 dr dϕ Wijkm−1n
(t) .
−

r 1
ϕ
1
k−

n−

2

2

If a∗ (~r) E∗ (t, ~x) < 0 in ∂µ± ΩI , q = 0, ..., 5 ,

Z
5
X

(A5∗ ) ≈
p=0

xi+ 1

2

Z

xi− 1

yj+ 1
2


q
p
(~x) dx dy 
(~x) ηi,j
E∗ (t, ~x) ηi,j

yj− 1

2
2


 Z rk+ 1 Z ϕn+ 1
2
2
p
× 
a | 1 ξp
| 1 ξq
| 1 dr dϕ Wijkm+1n
(t)
 r 1 ϕ 1 ∗ m+ 2 k,m+1,n m+ 2 k,m,n m+ 2
k−
n−
2
2



Z
Z r 1Z ϕ 1
n+
k+
2 δq2
2
2
p
p
q
p

a∗ |m− 12 ξk,m,n |m− 21 ξk,m,n |m− 21 dr dϕ +
a∗ (~r) ξk,m,n (~r) d~r WI (t) .
−
∆µm Kkmn
r 1
ϕ
1
k−

n−

2

2

Weak form for the term related to

Z
(A6∗ )
ΩI

∂
∂ϕ

(a6 Φ) is

∂
[a∗ (~r) E∗ (t, ~x) Φ(t, ~x, ~r)] ψ(~x, ~r) d~x d~r
∂ϕ

(50)

By the Upwind Flux Rule, if a∗ (~r) E∗ (t, ~x) > 0 in ∂ϕ± ΩI , q = 0, ..., 5 ,

Z
5
X

(A6∗ ) ≈
p=0

xi+ 1
2

xi− 1

Z

yj+ 1

2


p
q
E∗ (t, ~x) ηi,j
(~x) ηi,j
(~x) dx dy  ×

yj− 1

2
2


Z
 Z rk+ 1 Z µm+ 1
2 δq3
2
2
p
q
p

a∗ |n+ 21 ξk,m,n
|n+ 12 ξk,m,n
|n+ 12 drdµ −
a∗ (~r)ξk,m,n
(~r)d~r WIp (t)
 r 1 µ 1
∆ϕn Kkmn
k−
m−
2
2



Z r 1Z µ 1

k+
m+
2
2
p
q
p
−
a∗ |n− 21 ξk,m,n−1
|n− 12 ξk,m,n
|n− 12 dr dµ Wijkmn−1
(t) .

r 1
µ
1
k−

2

m−

2

If a∗ (~r) E∗ (t, ~x) < 0 in ∂ϕ± ΩI , q = 0, ..., 5 ,

41


Z
5
X

(A6∗ ) ≈
p=0

xi+ 1

2

xi− 1

Z

yj+ 1
2


q
p
(~x) dx dy  ×
(~x) ηi,j
E∗ (t, ~x) ηi,j

yj− 1

2
2


 Z rk+ 1 Z µm+ 1
2
2
p
q
p

a∗ |n+ 12 ξk,m,n+1
|n+ 12 ξk,m,n
|n+ 12 dr dµ Wijkmn+1
(t)
 r 1 µ 1
k−
m−
2
2



Z r 1Z µ 1
Z
k+
m+
2 δq3
2
2
p
q
p
−
a∗ |n− 21 ξk,m,n
|n− 21 ξk,m,n
|n− 21 dr dµ +
a∗ (~r) ξk,m,n
(~r) d~r WIp (t) .
∆ϕn Kkmn
rk− 1
µm− 1
2

2

References
[1] P. Markowich, C. Ringhofer and C. Schmeiser, Semiconductor Equations, Springer-Verlag, 1990.
[2] M. Lundstrom, Fundamentals of Carrier Transport, Cambridge University Press, 2000.
[3] C. Jacoboni and P. Lugli, The Monte Carlo Method for Semiconductor Device Simulation,
Spring-Verlag: Wien-New York, 1989.
[4] E. Fatemi and F. Odeh, finite difference solution of Boltzmann equation applied to electron
transport in semiconductor devices, Journal of Computational Physics, 108 (1993) 209-217.
[5] A. Majorana and R. Pidatella, A finite difference scheme solving the Boltzmann Poisson system
for semiconductor devices, Journal of Computational Physics, 174 (2001) 649-668.
[6] J.A. Carrillo, I.M. Gamba, A. Majorana and C.-W. Shu, A WENO-solver for 1D non-stationary
Boltzmann-Poisson system for semiconductor devices, Journal of Computational Electronics, 1
(2002) 365-375.
[7] J.A. Carrillo, I.M. Gamba, A. Majorana and C.-W. Shu, A direct solver for 2D non-stationary
Boltzmann-Poisson systems for semiconductor devices: a MESFET simulation by WENOBoltzmann schemes, Journal of Computational Electronics, 2 (2003) 375-380.
[8] J.A. Carrillo, I.M. Gamba, A. Majorana and C.-W. Shu, A WENO-solver for the transients
of Boltzmann-Poisson system for semiconductor devices. Performance and comparisons with
Monte Carlo methods, Journal of Computational Physics, 184 (2003) 498-525.
42

[9] M.J. Caceres, J.A. Carrillo, I.M. Gamba, A. Majorana and C.-W. Shu, Deterministic kinetic
solvers for charged particle transport in semiconductor devices, in Transport Phenomena and
Kinetic Theory Applications to Gases, Semiconductors, Photons, and Biological Systems. Series:
Modeling and Simulation in Science, Engineering and Technology. C. Cercignani and E. Gabetta
(Eds.), Birkhäuser (2006) 151-171.
[10] J.A. Carrillo, I.M. Gamba, A. Majorana and C.-W. Shu, 2D semiconductor device simulations
by WENO-Boltzmann schemes: efficiency, boundary conditions and comparison to Monte Carlo
methods, Journal of Computational Physics, 214 (2006) 55-80.
[11] M. Galler and A. Majorana, Deterministic and stochastic simulation of electron transport in
semiconductors, Bulletin of the Institute of Mathematics, Academia Sinica (New Series), 6th
MAFPD (Kyoto) special issue Vol. 2 (2) (2007) 349-365.
[12] M. Galler, Multigroup Equations for the Description of the Particle Transport in Semiconductors, Series on Advances in Mathematics for Applied Sciences - Vol. 70, World Scientific,
Singapore, 2005.
[13] Z. Chen, B. Cockburn, C. Gardner and J. Jerome, Quantum hydrodynamic simulation of
hysteresis in the resonant tunneling diode, Journal of Computational Physics, 274 (1995) 274280.
[14] Z. Chen, B. Cockburn, J. W. Jerome and C.-W. Shu, Mixed-RKDG finite element methods
for the 2-d hydrodynamic model for semiconductor device simulation, VLSI Design, 3 (1995)
145-158.
[15] Y.-X. Liu and C.-W. Shu, Local discontinuous Galerkin methods for moment models in device
simulations: formulation and one dimensional results, Journal of Computational Electronics, 3
(2004) 263-267.

43

[16] Y.-X. Liu and C.-W. Shu, Local discontinuous Galerkin methods for moment models in device
simulations: Performance assessment and two dimensional results, Applied Numerical Mathematics, 57 (2007) 629-645.
[17] Y. Cheng, I.M. Gamba, A. Majorana and C.-W. Shu, Discontinuous Galerkin Solver for the
Semiconductor Boltzmann Equation, SISPAD 07, T. Grasser and S. Selberherr, editors, Springer
(2007) 257-260.
[18] Y. Cheng, I. Gamba, A. Majorana and C.-W. Shu, Discontinuous Galerkin solver for
Boltzmann-Poisson transients, Journal of Computational Electronics, 7 (2008) 119-123.
[19] Y. Cheng, I. M. Gamba, A. Majorana and C.-W. Shu A discontinuous Galerkin solver for
Boltzmann-Poisson systems in nano-devices, Computer Methods in Applied Mechanics and
Engineering, 198 (2009) 3130-3150.
[20] Y. Cheng, I. M. Gamba, A. Majorana and C.W. Shu A discontinuous Galerkin solver for FullBand Boltzmann-Poisson Models, IWCE13 (13th International Workshop on Computational
Electronics) (2009).
[21] Y. Cheng, I. M. Gamba and J. Proft, Positivity-preserving discontinuous Galerkin schemes for
linear Vlasov-Boltzmann transport equations, Mathematics of Computation, 81 (2012) 153-190.
[22] J. Chelikowsky, and M. Cohen, Electronic structure of silicon, Physical Review B, 10 (1974)
5095-5107.
[23] M L. Cohen and J.R. Chelikowsky, Electronic Structure and Optical Properties of Semiconductors, Springer Series in Solid-State Sciences 75 (1989).
[24] S. Smirnov and C. Jungemann, A full band deterministic model for semiclassical carrier transport in semiconductors, Journal of Applied Physics, 99 (2006) 063707.
[25] C. Jungemann, A.T. Pham, B. Meinerzhagen, C. Ringhofer and M. Bollhöfer, Stable discretization of the Boltzmann equation based on spherical harmonics, box integration, and a maximum
entropy dissipation principle, Journal of Applied Physics, 100 (2006) 024502.
44

[26] A. Majorana, A comparison between bulk solutions to the Boltzmann equation and the spherical harmonic model for silicon devices, Progress in Industrial Mathematics at ECMI 2000 Mathematics in Industry, 1 (2002) 169-173.
[27] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite
element method for conservation laws II: general framework, Mathematics of Computation, 52
(1989) 411-435.
[28] B. Cockburn, S.-Y. Lin and C.-W. Shu, TVB Runge-Kutta local projection discontinuous
Galerkin finite element method for conservation laws III: one dimensional systems, Journal
of Computational Physics, 84 (1989) 90-113.
[29] B. Cockburn, S. Hou and C.-W. Shu, The Runge-Kutta local projection discontinuous Galerkin
finite element method for conservation laws IV: the multidimensional case, Mathematics of
Computation, 54 (1990) 545-581.
[30] B. Cockburn and C.-W. Shu, The Runge-Kutta local projection P1-discontinuous Galerkin
finite element method for scalar conservation laws, Mathematical Modelling and Numerical
Analysis, 25 (1991) 337-361.
[31] B. Cockburn and C.-W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation
laws V: multidimensional systems, Journal of Computational Physics, 141 (1998) 199-224.
[32] B. Cockburn and C.-W. Shu, Runge-Kutta discontinuous Galerkin methods for convectiondominated problems, Journal of Scientific Computing, 16 (2001) 173-261.
[33] B. Cockburn and C.-W. Shu, The local discontinuous Galerkin method for time-dependent
convection-diffusion systems, SIAM Journal on Numerical Analysis, 35 (1998) 2440-2463.
[34] D. Arnold, F. Brezzi, B. Cockburn, and L. Marini, Unified analysis of discontinuous Galerkin
methods for elliptic problems, SIAM Journal on Numerical Analysis, 39 (2002) 1749-1779.

45

[35] C. Cercignani, I.M. Gamba, C.L. Levermore,

A Drift-Collision Balance asymptotic for a

Boltzmann-Poisson System in Bounded Domains, SIAM J. Appl. Math., Vol 61, No. 6, (2001)
1932-1958.

46

