arXiv:1512.09210v3 [math.NA] 26 Feb 2018

Galerkin Methods for Boltzmann-Poisson
transport with reflection conditions on rough
boundaries
José A. Morales Escalante1 , Irene M. Gamba2
1

2

TU Wien - Institute for Analysis and Scientific Computing
The University of Texas at Austin - Institute for Computational
Engineering and Sciences & Department of Mathematics
Abstract
We consider in this paper the mathematical and numerical modelling
of reflective boundary conditions (BC) associated to Boltzmann - Poisson systems, including diffusive reflection in addition to specularity, in
the context of electron transport in semiconductor device modelling at
nano scales, and their implementation in Discontinuous Galerkin (DG)
schemes. We study these BC on the physical boundaries of the device
and develop a numerical approximation to model an insulating boundary
condition, or equivalently, a pointwise zero flux mathematical condition
for the electron transport equation. Such condition balances the incident
and reflective momentum flux at the microscopic level, pointwise at the
boundary, in the case of a more general mixed reflection with momentum
dependant specularity probability p(~k). We compare the computational
prediction of physical observables given by the numerical implementation
of these different reflection conditions in our DG scheme for BP models,
and observe that the diffusive condition influences the kinetic moments
over the whole domain in position space.
Keywords: Galerkin; Boltzmann-Poisson; boundary; reflection; diffusive; specular.

1

Introduction

The dynamics of electronic transport in modern semiconductor devices can be
described by the semiclassical Boltzmann-Poisson (BP) model
X
qi ~
∂fi
1
+ ∇~k εi · ∇~x fi − E
· ∇~k fi =
Qi,j ,
∂t
~
~
j
∇~x · ( ∇~x V ) =

X

qi ρi − N (~x),

~ = −∇~x V,
E

(1.1)

(1.2)

i

where fi (~x, ~k, t) is the probability density function (pdf) over phase space (~x, ~k)
of a carrier in the i-th energy band in position ~x, with crystal momentum ~~k at
1

time t. The collision operators Qi,j (fi , fj ) model i-th and j-th carrier recom~ x, t) is the electric
binations, collisions with phonons or generation effects. E(~
~
field, V (~x, t) is the electric potential, εi (k) is the i-th energy band surface, the
i-th charge density ρi (t, ~x) is the k-average of fi , −qi is the electric charge of
the i-th carrier, N (~x) is the doping profile, and  is the electric permittivity of
the material.
The BP model for electron transport on a single conduction energy band for
electrons has the form
∂f
1
q~
+ ∇~k ε(~k) · ∇~x f − E(~
x, t) · ∇~k f = Q(f ),
∂t
~
~

(1.3)

~ = −∇~x V,
E

(1.4)

∇~x · ( ∇~x V ) = q [ρ(~x, t) − N (~x)] ,

ε(~k), and the electron
with the quantum mechanical electron group velocity
R
~
~
density ρ(~x, t) = Ω~ f (~x, k, t) dk. The collision integral operator Q(f ) describes
k
the scattering over the electrons, where several mechanisms of quantum nature
can be taken into account. In the low density regime, the collisional integral
operator can be approximated as linear in f , having the form
Z h
i
Q(f ) =
S(~k 0 , ~k)f (t, ~x, ~k 0 ) − S(~k, ~k 0 )f (t, ~x, ~k) d~k 0 ,
(1.5)
1
k
~ ∇~

Ω~k

where S(~k, ~k 0 ) is the scattering kernel, representing non-local interactions of
electrons with a background density distribution. For example, in the case
of silicon, one of the most important collision mechanisms are electron-phonon
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, ~k 0 )

=

(nq + 1) K δ(ε(~k 0 ) − ε(~k) + ~ωp )
+ nq K δ(ε(~k 0 ) − ε(~k) − ~ωp ) + K0 δ(ε(~k 0 ) − ε(~k)) ,

(1.6)

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



−1
~ωp
nq = exp
−1
,
KB TL
where KB is the Boltzmann constant and TL = 300K is the lattice temperature.
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. 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 [14]. However, after the
pioneer work [15], in recent years, deterministic solvers to the BP system were
proposed in [16, 17, 18, 19, 20, 21, 22]. These methods provide accurate results
which, in general, agree well with those obtained from Monte Carlo (DSMC)
2

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 [18, 19, 20, 21] using weighted essentially nonoscillatory (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 and the simple communication pattern of
the discontinuous Galerkin (DG) methods for macroscopic (fluid level) models
[23, 24, 25, 26], it was proposed in [27, 28] 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 [27, 28], the first DG solver for (1.1)-(1.2) was proposed,
and some numerical calculations were shown for one and two-dimensional devices. In [29], 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 an analytic Kane energy band model.
A DG method for full conduction bands BP models was proposed in [30], following the lines of the schemes in [27, 28, 29], 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 direct evaluation of the Dirac delta function can be avoided, and so an accurate high-order
simulation with comparable computational cost to the analytic band cases is
possible. It would be more difficult or even unpractical to produce the full band
computation with other transport scheme. It is worth to notice that a highorder positivity-preserving DG scheme for linear Vlasov-Boltzmann transport
equations, under the action of quadratically confined electrostatic potentials, independent of the electron distribution, has been developed in [31]. The authors
there show that these DG schemes conserve mass and preserve the positivity
of the solution without sacrificing accuracy. In addition, the standard semidiscrete schemes were studied showing stability and error estimates.
The type of DG method discussed in this paper, as was done in [29], belongs to
a class of finite element methods originally devised to solve hyperbolic conservation laws containing only first order spatial derivatives, e.g. [32, 33, 34, 35, 36].
Using a piecewise polynomial space for both the test and trial functions in the
spatial variables, and coupled with explicit and nonlinearly stable high order
Runge-Kutta time discretization, the DG method is a conservative scheme that
has the advantage of flexibility for arbitrarily unstructured meshes, with a compact stencil, and with the ability to easily accommodate arbitrary hp-adaptivity.
For more details about DG scheme for convection dominated problems, we refer
3

to the review paper [37], later generalized to the Local DG (LDG) method to
solve the convection diffusion equations [38] and elliptic equations [39].
Regarding Boundary Conditions (BC), there are several kinds of BC for BP
semiconductor models. They vary according to the considered device and physical situation. We list below several examples of BC that could arise in the case
of electron transport along a single conduction band.

Charge neutrality boundary conditions, given by [4]
fout (t, ~x, ~k)

Γ

= ND (~x)

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

,

Γ

subset of ∂Ω~x ,

(1.7)

Γ

where Ω~x is the position domain. This BC is imposed in source and drain
boundaries, where electric currents enter or exit the device, to achieve neutral
charges there, as ρout (~x, t) − ND (~x) = 0.

Reflective BC happen in insulating boundaries, usually defined by a Neumann
boundary ΓN , of 2D and 3D devices. In general, reflective BC can be formulated
as the values of the pdf at the inflow boundary being dependent on the outflow
boundary values

f (~x, ~k, t)|ΓN − = FR f |ΓN + ,
(1.8)

FR f |ΓN + denoting that the reflection boundary condition is a function of the
outflow boundary values of the probability density function, where the Neumann
Inflow Boundary is defined as
Γ−
x, ~k) | ~x ∈ ΓN , ~k ∈ Ω~k , ~v (~k) · η(~x) < 0},
N = {(~

(1.9)

1
∇~ ε(~k) ,
(1.10)
~ k
with Ω~k the momentum domain, η(~x) outward unit normal, and the Neumann
Outflow Boundary is defined as
~v (~k) =

Γ+
x, ~k) | ~x ∈ ΓN , ~k ∈ Ω~k , ~v (~k) · η(~x) > 0} .
N = {(~

(1.11)

Specular Reflection BC over the Neumann Inflow Boundary is given by
f |− (~x, ~k, t) = FS (f |+ ) = f |+ (~x, ~k 0 , t)
(~x, ~k 0 ) ∈ Γ+
N,

k~0

for

(~x, ~k) ∈ Γ−
N,

t > 0,

s.t. ~v (~k 0 ) = ~v (~k) − 2 η(~x) · ~v (~k) η(~x) .

(1.12)
(1.13)

Diffusive reflection is a known condition from kinetic theory, in which the distribution function at the Inflow boundary is proportional to a Maxwellian [1],
[2] with T = TW = TW (~x) the temperature at the wall
~
f |− (~x, ~k, t) = FD (f |+ ) = C σ {f |+ } (~x, t) e−ε(k)/KB T , (~x, ~k) ∈ Γ−
N ,
Z
σ {f |+ } (~x, t) =
~v (~k) · η(~x)f |+ (~x, ~k, t)dk ,
~
v (~
k)·η>0

4

(1.14)
(1.15)

Z
C = C {η(~x)} =

|~v · η| e

−ε(~
k)/KB TL

d~k

−1
.

~
v ·η<0

Mixed reflection BC models the effect of a physical surface on electron transport
in metals and semiconductors, giving the reflected pdf representing the electrons
as a linear convex combination of specular and diffuse components, as in the
formula
f |− (~x, ~k, t)

= FM (f |+ ) = p FS (f |+ ) + (1 − p) FD (f |+ )
= p f |+ (~x, ~k 0 , t) + (1 − p) C 0 σ 0 {f |+ } (~x, t) e

(1.16)
ε(~
k)
BT

−K

,

(~x, ~k) ∈ Γ−
N .

p is sometimes called specularity parameter. It can either be constant or a function, dependant of the momentum. For example, the work by Soffer [5] studies a
statistical model for the reflection from a rough surface in electrical conduction.
It derives a specularity parameter p(~k) which depends on the momentum, given
by
2
2
2
p(~k) = e−4lr |k| cos Θ ,
(1.17)
where lr is the rms height of the rough interface, and Θ is th angle between the
incident electron and the interface surface normal.
Reflection BC is a widely studied topic in the context of the kinetic theory of
gases modelled by Boltzmann Equations. However, in the context of kinetic
models for electron transport in semiconductors, there is less extensive previous
work related to the study of the effect of reflection boundary conditions such
as diffusive, specular, or mixed reflection. An example of the list of references
where reflection BC are studied for Boltzmann equations in the context of kinetic theory of gases would include the works of Cercignani [3] and Sone [1],
where the specular, diffusive, and mixed reflection BC are formulated for the
Boltzmann Eq. for gases. V. D. Borman, S. Yu. Krylov, A. V. Chayanov [9]
study the nonequilibrium phenomena at a gas-solid interface. The recent paper of Brull, Charrier, Mieussens [10] studies the gas-surface interaction at a
nano-scale and the boundary conditions for the associated Boltzmann equation.
The recent work of Struchtrup [11] studies as well the Maxwell boundary condition and velocity dependent accommodation coefficients in the context of gases
mentioned. It considers the convex combination of specular reflection, isotropic
scattering, and diffusive reflection, incorporating velocity dependent coefficients
into a Maxwell-type reflection kernel. It develops a modification of Maxwell’s
BC, extending the Maxwell model by allowing it to incorporate velocity dependent accomodation coefficients into the microscopic description and satisfying
conditions of reciprocity and unitary probability normalization.
Regarding reflectivity in the context of Boltzmann models of electron transport,
Fuchs [6] proposed a boundary condition for the probability density function of
free electrons incident in the material surface, which is a convex combination of
specular & diffuse reflection with a constant specularity parameter p. Greene
([7], [8]) studied conditions for the Fuchs BC in which the specularity parameter
p(~k) is dependant on the angle of the momentum ~k, deriving a boundary condition for electron distributions at crystal surfaces valid for metal, semimetal, &
5

semiconductor surfaces, and showing that Fuchs’ reflectivity parameter differs
from the kinetic specularity parameter in physical significance and in magnitude.
It considers the unperturbed electron states of a crystal with an ideal perfectly
specular surface as standing wave states, and the diffusive reflection killing partially the incoming wave function. Soffer [5] studies a statistical model for the
electrical conduction, and derives under certain assumptions, such as a rough
surface random model with a Gaussian probability of height above or below
a horizontal plane, analytical formulas for a momentum dependant specularity
parameter p(~k) = exp(−4lr2 |k|2 cos2 Θ) associated to this physical phenomena,
abovementioned in (1.17). As mentioned before, lr is the rms height of rough
interface, and Θ is th angle between the incident electron and the interface surface normal.
The reference book of Markowich, Ringhofer, & Schmeiser [12] for semiconductor equations discusses the mathematical definition of boundaries according to
the physical phenomena, and defines accordingly the kind of BC to be imposed
at those boundaries: Dirichlet, Neumann, Inflow and Outflow boundaries. A
work of particular importance for us is the one by Cercignani, Gamba, and
Levermore [4]. They study high field approximations to a Boltzmann-Poisson
system and boundary conditions in a semiconductor. The BP system for electrons in a semiconductor in the case of high fields and small devices is considered. Boundary conditions are proposed at the kinetic level that yield charge
neutrality at ohmic contacts, which are Dirichlet boundaries, and at insulating Neumann boundaries. Both BC, either the one yielding charge neutrality at
Dirichlet boundaries, or the one rendering zero flux of electrons at the boundary,
assume that the pdf is proportional to a ground state associated to an asymptotic expansion of a dimensionless Boltzmann-Poisson system. Then they study
closures of moment equations and BC for both the pdf and for the moment
closures. The paper [40] also comments on the study of boundary conditions for
kinetic and macroscopic approximations for the Boltzmann - Poisson system in
bounded domains. Jüngel mentions in his semiconductors book [2] the different
kinds of reflection BC common on the kinetic theory of gases, specular, diffusive,
and mixed reflection but no further study of diffusive and mixed reflection BC
in the context of semiconductors is pursued.
We intend to present in this work a mathematical, numerical, and computational
study of the effect of diffusive, specular, and mixed reflection BC in BoltzmannPoisson models of electron transport in semiconductors, solved by means of
Discontinuous Galerkin FEM solvers. We study the mathematical formulation
of these reflection BC in the context of BP models for semiconductors, and
derive equivalent numerical formulations of the diffusive and mixed reflection BC
with non-constant p(~k), such that an equivalent numerical zero flux condition is
satisfied pointwise at the insulating Neumann boundaries at the numerical level.
We present numerical simulations for a 2D silicon diode and a 2D double gated
MOSFET, comparing the effects of specular, diffusive, and mixed reflection
boundary conditions in the physical observable quantities obtained from the
simulations.

6

2

BP system with ~k coordinate transformation
assuming a Kane Energy Band

The Kane Energy Band Model is a dispersion relation between the conduction
energy band ε (measured from a local minimum) and the norm of the electron
wave vector |k|, given by the analytical function (α is a constant parameter, m∗
is the electron reduced mass for Si, and ~ is Planck’s constant)
ε(1 + αε) =

~2 |k|2
2m∗

.

(2.18)

For our preliminary numerical studies we will use a Boltzmann-Poisson model
as in [29] , in which the conduction energy band is assumed to be given by a
Kane model. We use the following dimensionalized variables, with the related
characteristic parameters
t = t/t∗ , (x, y) = ~x/`∗ , `∗ = 10−6 m, t∗ = 10−12 s, V∗ = 1V .
A transformed Boltzmann transport equation is used as in [29] as well, where the
coordinates used to describe ~k are: µ, the cosine of the polar angle, the azimuthal
angle ϕ, and the dimensionless Kane Energy w = ε/KB T , which is assumed as
the conduction energy band. We will assume that the wall temperature is equal
to the lattice temperature, so TW = T = TL , and αK = αKB T . So ~k(w, µ, ϕ),
where
√ ∗
 p

p
p
~k = 2m KB TL w(1 + αK w) µ, 1 − µ2 cos ϕ, 1 − µ2 sin ϕ . (2.19)
~
A new unknown function Φ is used in the transformed Boltzmann Eq. [29] ,
which is proportional to the Jacobian of the transformation and to the density
of states (up to a constant factor)
Φ(t, x, y, w, µ, ϕ) = s(w)f (t, ~x, ~k) ,
where
s(w) =

p
w(1 + αK w)(1 + 2αK w) .

(2.20)

The transformed Boltzmann transport equation for Φ used in [29] is
∂Φ
∂
∂
∂
∂
∂
+
(g1 Φ) +
(g2 Φ) +
(g3 Φ) +
(g4 Φ) +
(g5 Φ) = C(Φ). (2.21)
∂t
∂x
∂y
∂w
∂µ
∂ϕ
The vector (g1 , g2 ) represent the 2D cartesian components of the electron velocity ~1 ∇~k ε(~k), in the coordinate system (w, µ, ϕ). The triplet (g3 , g4 , g5 )
represent the transport in the phase space of the new momentum coordinates
(w, µ, ϕ) due to the self consistent electric field
~ x, y) = (Ex (t, x, y), Ey (t, x, y), 0) ,
E(t,

7

with
p

w(1 + αK w)
µ,
1 + 2αK w
p
w(1 + αK w) p
g2 (·) = cx
1 − µ2 cos ϕ ,
1 + 2αK w
p
i
p
2 w(1 + αK w) h
µ Ex (t, x, y) + 1 − µ2 cos ϕ Ey (t, x, y) ,
g3 (·) = − ck
1 + 2αK w
p
2 w(1 + αK w)
~ x, y) ,
= − ck
êw · E(t,
1 + 2αK w
p
hp
i
1 − µ2
1 − µ2 Ex (t, x, y) − µ cos ϕ Ey (t, x, y) ,
g4 (·) = − ck p
w(1 + αK w)
p
1 − µ2
~ x, y) ,
= − ck p
êµ · E(t,
w(1 + αK w)
− sin ϕ
p
Ey (t, x, y)
g5 (·) = −ck p
w(1 + αK w) 1 − µ2
1
~ x, y) ,
p
= −ck p
êϕ · E(t,
w(1 + αK w) 1 − µ2
r
t∗ 2 KB TL
t∗ qE∗
cx =
and ck = √ ∗
,
`∗
m∗
2m KB TL
and êw , êµ , êϕ the orthonormal vector basis in our momentum coordinate space.
The right hand side of (2.21) is the collision operator (having applied the Dirac
Delta’s due to electron-phonon scattering, which depend on the energy differences between transitions)
 Z π
Z 1
C(Φ)(t, x, y, w, µ, ϕ) = s(w) c0
dϕ0
dµ0 Φ(t, x, y, w, µ0 , ϕ0 )
g1 (·)

=

cx

−1

0

Z
+
0

π

dϕ0

Z

1


dµ0 [c+ Φ(t, x, y, w + γ, µ0 , ϕ0 ) + c− Φ(t, x, y, w − γ, µ0 , ϕ0 )]

−1

− Φ(t, x, y, w, µ, ϕ) 2π [c0 s(w) + c+ s(w − γ) + c− s(w + γ)] ,
with the dimensionless parameters
2m∗ t∗ p
2 m∗ KB TL (K0 , (nq + 1)K, nq K) ,
(c0 , c+ , c− ) =
~3
The electron density is
Z
n(t∗ t, `∗ x, `∗ y) =

R3

√
f (t∗ t, `∗ x, `∗ y, k) dk =

where
Z
ρ(t, x, y) =

+∞

Z

1

dw
0

Z
dµ

−1

2 m∗ KB TL
~

γ=

~ωp
.
KB TL

3
ρ(t, x, y) ,

π

dϕ Φ(t, x, y, w, µ, ϕ) .

(2.22)

0

Hence, the dimensionless Poisson equation is




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

(2.23)

with
√
ND (x, y) =

3

2 m∗ KB TL
~

√

−3
ND (`∗ x, `∗ y) and cp =

2 m∗ KB TL
~

3

`2∗ q
.
0

Discontinuous Galerkin Method for Transformed
Boltzmann - Poisson System and Implementation of Boundary Conditions

The domain of the devices to be considered can be represented by means of
a rectangular grid in both position and momentum space. This rectangular
grid, bidimensional in position space and tridimensional in momentum space, is
defined as
Ωijkmn = Xij × Kkmn ,
i h
i
h
Xij = xi− 21 , xi+ 12 × yj− 12 , yj+ 21 ,
h
i h
i h
i
Kkmn = wk− 12 , wk+ 12 × µm− 12 , µm+ 12 × ϕn− 21 , ϕn+ 21 ,
where i = 1, . . . Nx , j = 1, . . . Ny , k = 1, . . . Nw , m = 1, . . . Nµ , n = 1, . . . Nϕ ,
xi± 21 = xi ±

∆xi
,
2

yj± 12 = yj ±

∆yj
,
2

∆wk
∆µm
∆ϕn
, µm± 21 = µm ±
, ϕn± 21 = ϕn ±
.
2
2
2
The finite dimensional space used to approximate the functions is the space
of piecewise continuous polynomials which are piecewise linear in (x, y) and
piecewise constant in (w, µ, ϕ) ,
wk± 12 = wk ±

Vh = {v : v|Ωijkmn ∈ Q1,0 (Ωijkmn ) = P 1 (Xij ) ⊗ P 0 (Kkmn )},

(3.24)

1,0
with the
h ) of tensor iproduct polynomials, linear over the element
h set Q (Ωiijkmn

Xij = xi− 21 , xi+ 12 × yj− 12 , yj+ 21 , and constant over the element
h
i h
i h
i
Kkmn = wk− 12 , wk+ 12 × µm− 21 , µm+ 12 × ϕn− 12 , ϕn+ 21 .

The function Φh will denote the piecewise polynomial approximation of Φ over
elements ΩI ,

X 
(x − xi )
(y − yj )
Φh =
χI TI (t) + XI (t)
+ YI (t)
, I = (i, j, k, m, n).
∆xi /2
∆yj /2
I

The density ρh (t, x, y) on the cell [xi− 21 , xi+ 12 ] × [yj− 21 , yj+ 21 ] is, under this

9

approximation,
ρh


Nµ Nϕ 
Nw X
X
X
(x − xi )
(y − yj )
Tijkmn + Xijkmn
+ Yijkmn
∆wk ∆µm ∆ϕn
∆xi /2
∆yj /2
m=1 n=1

=

k=1

Nµ Nϕ
Nw X
X
X

=

Tijkmn ∆wk ∆µm ∆ϕn

k=1 m=1 n=1


+ 

Nµ Nϕ
Nw X
X
X



(x − xi )
∆xi /2

Xijkmn ∆wk ∆µm ∆ϕn 

k=1 m=1 n=1


+ 

Nµ Nϕ
Nw X
X
X


Yijkmn ∆wk ∆µm ∆ϕn 

k=1 m=1 n=1

3.1

(y − yj )
.
∆yj /2

DG Formulation for Transformed Boltzmann Equation

The Discontinuous Galerkin formulation for the Boltzmann equation (2.21) is
as follows. Find Φh ∈ Vh , s.t.
Z
Z
Z
(Φh )t vh dΩ −
g1 Φh (vh )x dΩ −
g2 Φh (vh )y dΩ
Ωijkmn
Fx+ − Fx−

+

+

Fy+

−

Fy−

Ωijkmn
+ Fw+ −

Fw−

+

Fµ+

−

Fµ−

Ωijkmn
+ Fϕ+ −

Fϕ−

Z
=

C(Φh ) vh dΩ.

(3.25)

Ωijkmn

for any test function vh ∈ Vh . In (3.25), the boundary integrals are given by
Z y 1Z w 1Z µ 1Z ϕ 1
k+
m+
n+
j+
2
2
2
2
∓
gd
Fx± =
1 Φ vh (xi± 21 , y, w, µ, ϕ)dy dw dµ dϕ,
yj− 1

wk− 1

2

Fy± =

xi+ 1

Z

2

xi− 1

Z

wk+ 1

Z

2

2

=

Z

xi− 1

xi+ 1

Z

2

xi− 1

Fϕ± =

Z

2

2

2

yj+ 1

2

Z

2

2

Z

Z

µm+ 1
2

Z

wk+ 1
2

∓
gd
2 Φ vh (x, yj± 21 , w, µ, ϕ)dx dw dµ dϕ,

∓
gd
3 Φ vh (x, y, wk± 12 , µ, ϕ)dx dy dµ dϕ,

2

ϕn+ 1

Z

∓
gd
4 Φ vh (x, y, w, µm± 21 , ϕ)dx dy dw dϕ,

2

2

2

2

2

ϕn− 1

wk+ 1

wk− 1

ϕn+ 1

ϕn− 1

2

Z

2

2

µm− 1

Z

ϕn+ 1

ϕn− 1

2

yj+ 1

yj− 1

2

wk− 1

2

xi+ 1

xi− 1

yj+ 1

yj− 1

2

µm+ 1

2

2

Z

2

µm− 1

yj− 1

2

Fµ± =

Z

2

xi+ 1

ϕn− 1

2

wk− 1

2

Fw±

µm− 1

2

Z

µm+ 1
2

µm− 1

∓
gd
5 Φ vh (x, y, w, µ, ϕn± 21 )dx dy dw dµ,

2

where the upwind numerical fluxes gd
s Φ, s = 1, ..., 5 are defined as
10




g1 − |g1 |
g1 + |g1 |
−
Φh |xi±1/2 +
Φh |+
xi±1/2 ,
2
2




g2 + |g2 |
g2 − |g2 |
Φ h |−
Φh | +
yj±1/2 +
yj±1/2 ,
2
2




g3 + |g3 |
g3 − |g3 |
Φ h |−
Φh |+
wk±1/2 +
wk±1/2 ,
2
2




g4 + |g4 |
g4 − |g4 |
Φ h |−
Φh | +
µm±1/2 +
µm±1/2 ,
2
2




g5 + |g5 |
g5 − |g5 |
Φ h |−
Φh |+
ϕn±1/2 +
ϕn±1/2 . (3.26)
2
2


3.2

gd
1 Φ|xi±1/2

=

gd
2 Φ|yj±1/2

=

gd
3 Φ|wk±1/2

=

gd
4 Φ|µm±1/2

=

gd
5 Φ|ϕn±1/2

=

Poisson Equation - Local Discontinuous Galerkin (LDG)
Method

The Poisson equation (2.23) is solved by the LDG method as in [29] .
By means of this scheme we find a solution Ψh , qh , sh ∈ Wh1 , where (q, s) =
(∂x Ψ, ∂y Ψ) and Wh1 = {v : v|Xij ∈ P 1 (Xij )}, P 1 (Xij ) the set of linear polynomials on Xij . It involves rewriting the equation into the form

∂Ψ
∂Ψ

 q=
,
s=
,
∂x
∂y
∂
∂


(r q) +
(r s) = R(t, x, y) ,
∂x
∂y

(3.27)

where R(t, x, y) = cp [ρ(t, x, y) − ND (x, y)] is a known function that can be
computed at each time step once Φ is solved from (3.25), and the coefficient
r depends on x, y. The
h Poisson system
i h is only onithe (x, y) domain. Hence,

we use the grid Iij = xi− 21 , xi+ 12 × yj− 12 , yj+ 12 , with i = 1, . . . , Nx , j =
1, . . . , Ny + My , where j = Ny + 1, . . . , Ny + My denotes the oxide-silicon region,
and the grid in j = 1, . . . , Ny is consistent with the five-dimensional rectangular
grid for the Boltzmann equation in the silicon region. The approximation space
is defined as
Wh` = {v : v|Iij ∈ P ` (Iij )}.
(3.28)

Here P ` (Iij ) denotes the set of all polynomials of degree at most ` on Iij . The

11

LDG scheme for (3.27) is: to find qh , sh , Ψh ∈ Vh` , such that
Z
Z y 1h
i
j+
2
Ψ̂h vh+ (xi− 12 , y) − Ψ̂h vh− (xi+ 21 , y) dy,
0 =
[qh vh + Ψh (vh )x ] dxdy +
Iij

yj− 1
2

Z
0

=

Z

xi+ 1
2

[sh wh + Ψh (wh )y ] dxdy +
Iij

xi− 1

Z

yj+ 1
2

r qh (ph )x dxdy +
Ii,j

yj− 1

−
c
r q h ph (xi+ 12 , y)dy

Z

Ii,j

2

xi− 1

2

−

yj− 1

+
c
r q h ph (xi− 21 , y)dy

2

xi+ 1

r sh (ph )y dxdy +

yj+ 1

Z

2

Z
−

i
Ψ̃h wh+ (x, yj− 21 ) − Ψ̃h wh− (x, yj+ 12 ) dx,

2

Z
−

h

−
f
r sh ph (x, yj+ 12 )dx −

2

Z

xi+ 1
2

xi− 1

+
f
r sh ph (x, yj− 12 )dx

2

Z
=

R(t, x, y)ph dxdy ,

(3.29)

Ii,j

hold true for any vh , wh , ph ∈ Wh` . In the above formulation, we choose the flux
+
as follows, in the x-direction, we use Ψ̂h = Ψ−
r q h = r qh − [Ψh ]. In the
h , c
−
+
y-direction, we use Ψ̃h = Ψh , f
r sh = r sh − [Ψh ]. On some part of the domain
boundary, the above flux needs to be changed to accommodate various boundary
conditions. For example, in the case of a double gate MOSFET device, for the
boundary condition of the Poisson equation, Ψ = ΨS at source, Ψ = ΨD at drain
and Ψ = ΨG at gate. For the rest of the boundary regions, we have homogeneous
Neumann boundary conditions, i.e., ∂Ψ
∂n = 0. The relative dielectric constant
in the oxide-silicon region is r = 3.9, in the silicon region is r = 11.7. Near
the drain then, we are given Dirichlet boundary condition, so we need to flip
the flux in x−direction: let Ψ̂h (xi+ 12 , y) = Ψ+
r q h (xi+ 21 , y) =
h (xi+ 21 , y) and c
−
r qh (xi+ 21 , y) − [Ψh ](xi+ 21 , y), if the point (xi+ 21 , y) is at the drain. For the
gate, we need to flip the flux in y−direction: let Ψ̃h (x, yj+ 12 ) = Ψ+
h (x, yj+ 12 )
−
and f
r sh (x, yj+ 21 ) = r sh (x, yj+ 12 ) − [Ψh ](x, yj+ 12 ), if the point (x, yj+ 21 ) is at
the gate. For the bottom, we need to use the Neumann condition, and flip the
−
flux in y-direction, i.e., Ψ̃h = Ψ+
r sh = r sh . This scheme described above
h , f
∂Ψ
will enforce the continuity of Ψ and r ∂n across the interface of silicon and
oxide-silicon interface. The solution of (3.29) gives us approximations to both
the potential Ψh and the electric field (Ex )h = −cv qh , (Ey )h = −cv sh .

3.3

RK-DG Algorithm for BP, from tn to tn+1

The following RK-DG algorithm for BP is a dynamic extension of the Gummel
iteration map. We write below the steps to evolve from time tn to time tn+1 .
1. Compute the electron density ρh (x, y, t).
2. Solve Poisson Eq. for the given ρh (x, y, t) by Local DG, obtaining the
potential Ψh and the electric field Eh = −(qh , sh ). Compute then the
respective transport terms gs , s = 1, ..., 5.
3. Solve by DG the advection and collision part of the Boltzmann Equation.
A Method of Lines (an ODE system) for the time dependent coefficients
of Φh (degrees of freedom) is obtained.
12

4. Evolve ODE system by Runge-Kutta from tn to tn+1 . (If partial time
step necessary, repeat Step 1 to 3 as needed).

4

Boundary Conditions Implementation for 2D~x, 3D-~k devices at x, w, µ, ϕ Boundaries

We will consider in this work 2D devices in position space, which need a 3D
momentum description for kinetic equations modeling semiconductors. For example, a common device of interest is a 2D double gate MOSFET. A schematic
plot of it is given in Figure 4.1. The shadowed region denotes the oxide-silicon
region, whereas the rest is the silicon region. Potential bias are applied at the
source, drain, and gates. The problem is symmetric about the x-axis.
Another possible 2D problem is the case of a bi-dimensional bulk silicon diode,
for which the doping is constant all over the physical domain, and which would
have just an applied potential (bias) between the source x = 0 and the drain
x = Lx (no gates), with insulating reflecting boundaries at y = 0 and y = Ly .

Figure 4.1: Schematic representation of a 2D double gate MOSFET device.
From 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, v198 (2009), p. 3143.
We consider in the following sections the different kinds of boundary conditions
for 2D devices and their numerical implementation, either at ~x-boundaries or
at w-boundaries.
~

4.1

Poisson Equation Boundary Condition

The BC for Poisson Eq. are imposed over the (x, y)-domain.
For example, for the case of a 2D Double gated MOSFET, Dirichlet BC would
be imposed to the potential Ψ, as we have three different applied potentials
biases, Ψ = 0.5235 Volts at the source x = 0, Ψ = 1.5235 Volts at the drain
x = Lx , Ψ = 1.06 Volts at the gates. Homogeneous Neumann BC would be
imposed for the rest of the boundaries, that is, ∂n̂ Ψ = 0.
For the case of a 2D bulk silicon diode, we impose Dirichlet BC for the difference
of potential Ψ between source and drain, Ψ = 0.5235 Volts at the source x = 0,
13

Ψ = 1.5235 Volts at the drain x = Lx . For the boundaries y = 0, Ly we impose
Homogeneous Neumann BC too, that is, ∂y Ψ|y0 = 0, y0 = 0, Ly .

4.2

Charge Neutrality BC

As in [29], at the source and drain contacts, we implement the charge neutrality
boundary condition (1.7). Ghost cells for i = 0 and i = Nx + 1 at the respective
boundaries are used, implementing numerically the boundary conditions
Φ(i = 0) = Φ(i = 1)

ND (i = 1)
,
ρ(i = 1)

Φ(i = Nx + 1) = Φ(i = Nx )

4.3

ND (i = Nx )
.
ρ(i = Nx )

Cut - Off BC

In the (w, µ, ϕ)-space, we only need to apply a cut-off Boundary Condition. At
w = wmax , Φh is made machine zero,
Φh (x, y, w, µ, ϕ, t)|w=wmax = 0.

(4.30)

No other boundary condition is necessary for w-boundaries,
~
since analytically
we have that
• at w = 0, g3 = 0,
• at µ = ±1, g4 = 0,
• at ϕ = 0, π, g5 = 0,
so, at such regions, the numerical flux always vanishes.

5

Reflection BC on BP

Reflection Boundary Conditions can be expressed in the form
f (~x, ~k, t)|ΓN − = FR (f |ΓN + ),

(5.31)

such that the following pointwise zero flux condition is satisfied at reflecting
boundaries, so
Z
0 = η(~x) · J(~x, t) = η(~x) ·
~v (~k) f (~x, ~k, t) d~k ,
(5.32)
Ω~k

Z

Z
~
~
~
η(~x) · ~v (~k) f (~x, ~k, t)|ΓN − d~k ,
0 =
η(~x) · ~v (k) f (~x, k, t)|ΓN + dk +
η·~
v <0
η·~
v >0
Z
Z
~
0 =
~v · η f |ΓN + dk +
~v · η FR (f |ΓN + ) d~k ,
~
v ·η>0

~
v ·η<0

as in Cercignani, Gamba, Levermore[4], where the given BC at Neumann boundary regions at the kinetic level is such that the particle flow vanishes.
14

For simplicity we write ~v = ~v (~k) = ∇~k ε(~k)/~. We will study three kinds of
reflective boundary conditions: specular, diffusive, and mixed reflection. The
last one is a convex combination of the previous two, but the convexity parameter can be either constant or momentum dependant, p(~k). We go over the
mathematics and numerics related to these conditions below.

5.1

Specular Reflection

It is clear that, at the analytical level, the specular reflection BC (1.12) satisfies
the zero flux condition pointwise at reflecting boundaries, since
Z
Z
~
~
~
d~k = 0.
|η(~x) · ~v (~k)| f (~x, ~k 0 , t)
dk −
|η(~x) · ~v (k)| f (~x, k, t)
Γ +
N

η·~
v >0

Γ

−η·~
v <0

N+

Specular reflection BC in our transformed Boltzmann Eq. for the new coordinate system is mathematically formulated in our problem as
(x, y, w, µ, ϕ) ∈ Γ−
N.

Φ|− (x, y, w, µ, ϕ, t) = Φ|+ (x, y, w, µ, π − ϕ, t),

(5.33)

To impose numerically specular reflection BC at y = 0, Ly in the DG method,
we follow the procedure of [29]. We relate the inflow values of the pdf, associated
to the outer ghost cells, to the outflow values of the pdf, which are associated
to the interior cells adjacent to the boundary, as given below by
Φh |− (x, y1/2 , w, µ, ϕ, t) = Φh |+ (x, y1/2 , w, µ, π − ϕ, t),

y1/2 = 0,

Φh |− (x, yNy + 12 , w, µ, ϕ, t) = Φh |+ (x, yNy + 12 , w, µ, π − ϕ, t),

(5.34)

yNy + 12 = Ly .

In the case of the boundary y1/2 = 0, assuming ∆y0 = ∆y1 , ∆ϕn0 = ∆ϕn , with
n0 = Nϕ − n + 1, if (x, y1/2 − y, w, µ, ϕ) ∈ Ωi0kmn then (x, y1/2 + y, w, µ, π − ϕ) ∈
Ωi1kmn0 . The values of Φh |±
y1/2 at the related inner and outer boundary cells
Ωi0kmn (j = 0) and Ωi1kmn0 (j = 1) must be equal at the boundary y1/2 = 0.
Indeed
+
Φh |−
Ωi0kmn (x, y1/2 , w, µ, ϕ, t) = Φh |Ωi1kmn0 (x, y1/2 , w, µ, π − ϕ, t) =⇒

(y1/2 − y0 )
(x − xi )
+ Yi0kmn
=
∆xi /2
∆y0 /2
(y1/2 − y1 )
(x − xi )
+ Xi1kmn0
+ Yi1kmn0
.
∆xi /2
∆y1 /2

Ti0kmn + Xi0kmn
Ti1kmn0

Therefore, from the equality above we find the relation between the coefficients
of Φh at inner and outer adjacent boundary cells, given by
Ti0kmn = Ti1kmn0 , Xi0kmn = Xi1kmn0 , Yi0kmn = −Yi1kmn0 .

(5.35)

Following an analogous procedure for the boundary yNy +1/2 , we have
+
Φh |−
Ωi,Ny +1,kmn (x, yNy + 21 , w, µ, ϕ, t) = Φh |Ωi,N

y ,kmn0

(x, yNy + 12 , w, µ, π − ϕ, t) .

Then
(yNy + 12 − yNy +1 )
(x − xi )
+ Yi,Ny +1,kmn
=
∆xi /2
∆yNy +1 /2
(yNy + 12 − yNy )
(x − xi )
+ Xi,Ny ,kmn0
+ Yi,Ny ,kmn0
,
(5.36)
∆xi /2
∆yNy /2

Ti,Ny +1,kmn + Xi,Ny +1,kmn
Ti,Ny ,kmn0

15

and hence
Ti,Ny +1,kmn = Ti,Ny ,kmn0 , Xi,Ny +1,kmn = Xi,Ny ,kmn0 , Yi,Ny +1,kmn = −Yi,Ny ,kmn0 .

5.2

Diffusive Reflection

The diffusive reflection BC can be formulated as
~
f (~x, ~k, t)|− = FD (f |+ ) = C σ {f |+ } (~x, t) e−ε(k)/KB TL ,

(~x, ~k) ∈ Γ−
N , (5.37)

where σ {f |+ } (~x, t) = σ(~x, t) and C = C{η(~x)} are the function and parameter
such that the zero flux condition is satisfied at each of the points of the Neumann
Boundary, so
Z
Z
i
h
~
0 =
~v · η f |ΓN + d~k +
~v · η Cσ(~x, t)e−ε(k)/KB TL d~k ,
~
v ·η>0
~
v ·η<0
Z
Z
~
0 =
~v · η f |ΓN + d~k − σ(~x, t) · C
|~v · η| e−ε(k)/KB TL d~k .
~
v ·η>0

~
v ·η<0

It follows then that
Z
σ {f |+ } (~x, t) =

~
v (~
k)·η>0

Z

~v · η f |ΓN + (~x, ~k, t) d~k ,

~
|~v · η| e−ε(k)/KB TL d~k

C {η(~x)} =

(5.38)

−1
,

(5.39)

~
v ·η<0
~

e−ε(k)/KB TL
f (~x, ~k, t)|− =
R

R

~
v ·η<0

~
v (~
k)·η>0

~v · η f |ΓN + (~x, ~k, t) d~k

|~v · η| e−ε(~k)/KB TL d~k

.

(5.40)

The diffusive reflection BC, formulated in terms of the unknown function Φ of
the transformed Boltzmann Equation 2.21, is expressed as
Φ|− (x, y, w, µ, ϕ, t) = FD (Φ|+ ) = C σ {Φ|+ } (x, y, t) e−w s(w) ,
Z
σ(x, y, t) =
η · (g1 , g2 )(w, µ, ϕ) Φ|+ dwdµdϕ ,

(5.41)
(5.42)

(g1 ,g2 )·η>0

!−1

Z
|(g1 , g2 ) · η| e

C(η) =

−w

s(w) dwdµdϕ

.

(5.43)

(g1 ,g2 )·η<0

We have, over the portion of the boundary considered, that η = (0, −1, 0) for
y = 0 and η = (0, 1, 0) for y = Ly . Therefore
R
e−w s(w) −g2 >0 |g2 | Φ|+ dwdµdϕ
Φ|− (x, yb , w, t) = R
, yb = 0 ,
(5.44)
|g | e−w s(w) dwdµdϕ
−g2 <0 2
e−w s(w)
Φ|− (x, yb , w, t) = R

R

|g | Φ|+ dwdµdϕ
+g2 >0 2
−w
|g | e s(w) dwdµdϕ
+g2 <0 2

16

,

yb = Ly .

(5.45)

5.2.1

Numerical Formulation of Diffusive BC for DG

For the DG numerical method, we have to project the boundary conditions to
be imposed in the space Vh . Our goal is to have at the numerical level an
equivalent pointwise zero flux condition at the reflection boundary regions.
We formulate then the diffusive BC for the DG method as
Φh |− (x, yb , w, µ, ϕ, t)

=
=

Πh {FD (Φh |+ )}

Πh C σh {Φh |+ } (x, yb , t) e−w s(w) ,

yb = 0, Ly ,

where Πh is the projector of functions into the finite element space Vh , σh ∈ Vh
is a function in our piecewise polynomial space for (x, y) and C is a parameter
such that the zero flux condition is satisfied numerically, so
Z
Z
0 =
~g · η Φh |+ dw
~+
~g · η Φh |− dw
~
~
g ·η>0
~
g ·η<0
Z
Z
=
~g · η Φh |+ dw
~+
~g · η Πh {FD (Φh |+ )} dw
~
(5.46)
~
g ·η>0
~
g ·η<0
Z
Z

=
~g · η Φh |+ dw
~+
~g · η Πh C σh {Φh |+ } (x, yb , t) e−w s(w) dw
~.
~
g ·η>0

~
g ·η<0

In the space Vh of piecewise continuous polynomials which are tensor products
of polynomials of degree p in ~x and of degree q in w,
~ it holds that
Πh {f1 (~x)f2 (w)}
~ = Πh {f1 (~x)} Πh {f2 (w)}
~ ,
p,q

Vh = {v : v|Ωijkmn ∈ Q

p

(5.47)
q

(Ωijkmn ) = P (Xij ) ⊗ P (Kkmn )}.

Therefore, for our particular case we have


Πh C σh (x, yb , t) e−w s(w) = C σh (x, yb , t) Πh e−w s(w) ,

(5.48)

so for the numerical zero flux condition pointwise we have that
Z
Z

0 =
~g · η Φh |+ dw
~+
~g · η C σh {Φh |+ } (x, yb , t) Πh e−w s(w) dw
~
~
g ·η>0
~
g ·η<0
Z
Z

0 =
~g · η Φh |+ dw
~ − σh {Φh |+ } (x, yb , t) C
|~g · η| Πh e−w s(w) dw
~.
~
g ·η>0

~
g ·η<0

We observe then that we can obtain a numerical equivalent of the pointwise zero
flux condition if we define
Z
−1
 −w
C {η} = C {±ŷ} =
|~g · η| Πh e s(w) dw
~
, η = ±ŷ .
±g2 =~
g ·η<0
Z
σh {Φh |+ } (x, yb , t) =
~g · η Φh |+ dw
~ = σ {Φh |+ } (x, yb , t) , yb = 0, Ly .
±ŷ·~
g >0

In our particular case, in which we have chosen our function space as piecewise linear in (x, y) and piecewise constant in (w, µ, ϕ), the projection of the
Maxwellian is a piecewise constant approximation representing its average value
over each momentum cell , that is
R wk+ −w
R
−w
X
X
e s(w)dw
 −w
e
s(w)dwdµdϕ
w
=
χkmn k−
.
Πh e s(w) =
χkmn kmn
∆wk ∆µm ∆ϕn
∆wk
k,m,n

k,m,n

17

Therefore, for the particular space we have chosen, we have that
Z
σh {Φh |+ } (x, yb , t) =
±g2 Φh |+ dw
~ = σ {Φh |+ } (x, yb , t) ,

(5.49)

±g2 >0

(η = −ŷ) , or yb = Ly = yNy +1/2 (η = +ŷ) ,
Z wk+1/2
Z
1
e−w s(w) dw
C −1 =
|g2 | dw dµ dϕ, η = ±ŷ ,
∆wk wk−1/2
k,m,n
k,m,n

Φh |− (x, yb , w, µ, ϕ, t) = C σh {Φh |+ } (x, yb , t) Πh e−w s(w) , yb = 0, Ly ,
R −w
R
P±g2 <0
e
s(w) dw
k
|g
|
Φ
|
d
w
~
χ
2
h
+
kmn
k,m,n
∆wk
±g2 >0
R
Φh |− (x, yb , w, µ, ϕ, t) =
.
P±g2 <0 R
e−w s(w) dw
|g | dw dµ dϕ k ∆wk
k,m,n
kmn 2
yb = 0 = y1/2
±g
2 <0
X

By the upper index ±g2 < 0 in a sum we mean to say that the sum is taken
over the values of k, m, n for which ±g2 = ~g · η < 0. We notice that the
polynomial approximation σh is equal to the analytical function σ operating on
the polynomial approximation Φh |+ . However, the constant C needed in order
to achieve the zero flux condition numerically is not equal to the value of this
parameter in the analytical solution. In this case C is an approximation of the
analytical value using a piecewise constant approximation of the Maxwellian
(its average over cells).
The approximate operator σh {Φh |+ } (x, y, t) gives a piecewise linear polynomial
dependant on (x, y) with time dependent coefficients. We have that
Z
Φh |+ ∈ Vh =⇒ σh {Φh |+ } (x, y, t) =
|g2 | Φh |+ dwdµdϕ ∈ Vh ,
± cos ϕ>0

where Φh |+ is such that, at the boundary y = yb of the cell Ωijkmn , it is given
by
Φh |+ (t, x, y, w, µ, ϕ) = Tijkmn (t) + Xijkmn (t)

2(y − yj )
2(x − xi )
+ Yijkmn (t)
.
∆xi
∆yj

We define I = ijkmn, so in ΩI = Xij × Kkmn . Then,
σh (x, y, t) = σI0 (t) + σIx (t)

(x − xi )
(y − yj )
+ σIy (t)
.
∆xi /2
∆yj /2

(5.50)

We summarize the main results of these calculations for σh and Φh |− , by showing
just the ones related to y = Ly (the case y = 0 is analogous). At the boundary
y = Ly , the inner cells associated to outflow have j = Ny , adjacent to the
boundary, whereas the ghost cells related to inflow have the index j = Ny + 1.
We compute the integral σh as
p
Z
w(1 + αK w) p
σh {Φh |+ } (x, y, t) =
1 − µ2 cos ϕ Φh |+ dwdµdϕ
1 + 2αK w
cos ϕ≥0
N
p
n≤ 2p Z
X
w(1 + αK w) p
1 − µ2 cos ϕ Φh |+ dwdµdϕ.
=
1 + 2αK w
Kkmn
k,m,n

18

Therefore, we have, with I = (i, j, k, m, n), j = Ny below, that
n≤

σI0

Nϕ

X2

=

TiNy kmn
wk−1/2

k,m,n

n≤

σIx

σIy

p

XiNy kmn
wk−1/2

w(1 + αK w)
dw
1 + 2αK w

Z

µm+1/2

Z
p
1 − µ2 dµ

µm−1/2

ϕn+1/2

cos ϕdϕ,

ϕn−1/2

Nϕ

X2

=

wk+1/2

Z

k,m,n

n≤

p
Z ϕn+1/2
Z µm+1/2 p
w(1 + αK w)
2
1 − µ dµ
cos ϕdϕ,
dw
1 + 2αK w
ϕn−1/2
µm−1/2

Nϕ

X2

=

wk+1/2

Z

Z

wk+1/2

p

YiNy kmn
wk−1/2

k,m,n

w(1 + αK w)
dw
1 + 2αK w

Z

µm+1/2

Z
p
1 − µ2 dµ

µm−1/2

ϕn+1/2

cos ϕdϕ.

ϕn−1/2

Once the coefficients of σh have been computed, we use them to obtain the
polynomial approximation Φh |− , with j = Ny + 1, from (5.49)
Nϕ

Φh |−
y=Ly

=

X2
X n≥
i

k,m,n


 R −w
e s(w)dw
y
0
x (x − xi )
χiNy kmn C σI + σI
+ σI · 1 k
.
∆xi /2
∆wk

We have at the same time, by definition, that
n≥

Φh |−
y=Ly =

Nϕ

X2



(x − xi )
χi,Ny +1,kmn Ti,Ny +1,k,m,n + Xi,Ny +1,k,m,n
− 1 · Yi,Ny +1,k,m,n .
∆xi /2

ikmn

Therefore, the coefficients for Φh |−
y=Ly are
R
Ti,Ny +1,kmn (t) =

0
CσiN
(t) k
y kmn

e−w s(w)dw
,
∆wk

e−w s(w)dw
,
∆wk
R −w
e s(w)dw
y
k
Yi,Ny +1,kmn (t) = −1 · CσiN
,
(t)
y kmn
∆wk

(5.51)

R

Xi,Ny +1,kmn (t) =

x
CσiN
(t) k
y kmn

(5.52)

(5.53)

keeping in mind that our parameter C is given by the formula
n≥

C

−1

=

Np

X2

kmn

5.3

R
k

e−w s(w)dw
∆wk

Z p
Z p
Z
w(1 + αK w)
2
dw
1 − µ dµ cos ϕdϕ.
1 + 2αK w
k
m
n

Mixed Reflection

The mixed reflection condition is a convex combination of the specular and
diffusive reflections:
~
f (~x, ~k, t)|− = pf |+ (~x, ~k 0 , t) + (1 − p)Cσ {f |+ } (~x, t)e−ε(k)/KB T ,

19

(~x, ~k) ∈ Γ−
N ,

p is the Specularity Parameter, 0 ≤ p ≤ 1. p can be either constant or p = p(~k),
a function of the wave vector momentum.
For p constant, it can be shown easily that the previous formulas obtained for the
specular and diffusive BC, in particular the previous formulas for σ C(x), works
also in this case to obtain a zero flux condition at the Neumann boundaries:


Z
Z
−ε(~
k)
0
K B TL
~
~
σ(~x, t) d~k
η·J =
~v · ηf |+ dk +
~v · η pf (~x, k , t)|+ + (1 − p)Ce
~
v ·η>0
~
v ·η<0
Z
Z
Z
−ε(~
k)
0
~
~
=
~v · ηf |+ dk + p
~v · ηf |+ dk + (1 − p) σC
~v · ηe KB TL d~k
~
v ·η>0

~
v ·η<0

~
v ·η<0

= σ(~x, t) − pσ(~x, t) + (1 − p) σ(~x, t) (−1) = 0 .
However, for p(~k) a function of the crystal momentum the same choice of σ(~x, t)
and C(x) as in the diffusive case does not necessarily guarantee that the zero flux
condition will be satisfied at Neumann boundaries. Therefore, a new condition
for C in order to satisfy this condition must be derived. We derive it below.
The general mixed reflection BC can be formulated as
~
f (~x, ~k, t)|− = p(~k)f |+ (~x, ~k 0 , t) + (1−p(~k)) C 0 σ 0 {f |+ } (~x, t) e−ε(k)/KB T ,

(~x, ~k) ∈ Γ−
N

where σ 0 {f |+ } (~x, t) and C 0 are the function and parameter such that the pointwise zero flux condition is satisfied at the Neumann boundaries
0

= η(~x) · J(~x, t)
Z
Z
=
~v · η f |+ d~k +
~
v ·η>0

Since
Z
0=



−ε(~
k)
~v · η p(~k)f (~x, ~k 0 , t)|+ + (1 − p(~k))C 0 e KB TL σ 0 (~x, t) d~k.

~
v ·η<0

~v ·ηf |+ d~k+

Z

~
v ·η>0

~v ·η p(~k)f 0 |+ d~k − σ 0 (~x, t) C 0

~
v ·η<0

we conclude then that
Z
0
σ {f |+ } (~x, t) =

−ε

(1−p(~k))|~v ·η| e KB TL d~k ,

~
v ·η<0

~v · η f |+ d~k −

~
v ·η>0

C 0 {η(~x)} =

Z

Z

|~v · η| p(~k) f (~x, ~k 0 , t)|+ d~k , (5.54)

~
v ·η<0

Z

−ε

(1 − p(~k))|~v · η| e KB TL d~k

−1
.

(5.55)

~
v ·η<0

The general mixed reflection BC then has the specific form
f (~x, ~k, t)|−

= p(~k) f |+ (~x, ~k 0 , t)
+

with

(1 − p(~k))e

(~x, ~k) ∈ Γ−
N ,

ε(~
k)
−K T
B

R
~
v ·η>0


R
~v · ηf |+ d~k − ~v·η<0 |~v · η|p(~k)f (~x, ~k 0 , t)|+ d~k
,
−ε(~
k)
R
~k))|~v · η| e KB TL d~k
(1
−
p(
~
v ·η<0

(~x, ~k 0 ) ∈ Γ+
v (~k 0 ) = ~v (~k) − 2(~v (~k) · η)η .
N s.t. ~

Notice that the product C 0 σ 0 (~x, t) has the form

20


R
~k −
~k) f (~x, ~k 0 , t)|+ d~k
~
v
·
η
f
|
d
|~
v
·
η|
p(
+
~
v ·η>0
~
v ·η<0
−ε(~
k)
R
~k))|~v · η| e KB TL d~k
(1
−
p(
~
v ·η<0

R
C 0 σ 0 (~x, t) =

(5.56)

which for the case of p constant, it reduces to the original function σ(~x, t) and
parameter C {η(~x)}.
If

p

C 0 σ 0 (~x, t)

=
=

=

=

ct,
R


R
~v · η f |+ d~k − p ~v·η<0 |~v · η| f (~x, ~k 0 , t)|+ d~k
−ε(~
k)
R
(1 − p)|~v · η| e KB TL d~k
~
v ·η<0
R
(1 − p) ~v·η>0 ~v · η f |+ d~k
−ε(~
k)
R
(1 − p) ~v·η<0 |~v · η| e KB TL d~k
R
~v · η f |+ d~k
~
v ·η>0
−ε(~
k)
R
|~v · η| e KB TL d~k
~
v ·η<0
~
v ·η>0

= C σ (~x, t) .
However, for the non-constant case p(~k) the new function and parameter σ 0 (~x, t),
C 0 (η) need to be used instead, as the previous σ(~x, t), C(η) will not satisfy the
zero flux condition in general for p(~k), since
Z
Z
Z
−ε
0 =
~v · ηf |+ d~k +
~v · ηp(~k)f 0 |+ d~k − σ 0 C 0
(1 − p(~k))|~v · η|e KB TL d~k
~
v ·η>0

~
v ·η<0

~
v ·η<0

R
~v · η f |+ d~k + ~v·η<0 ~v · η p(~k) f (~x, ~k 0 , t)|+ d~k
~
v ·η>0
−ε
R
(1 − p(~k))|~v · η| e KB TL d~k
~
v ·η<0
R
~v · η f |+ d~k
~
v ·η>0
= Cσ(~x, t) in general for p(~k).
−ε(~
k)
R
KB TL
~
|~v · η| e
dk
~
v ·η<0
R

C 0 σ0

=

6=

A more general possible case of mixed reflection BC would have a specularity
parameter p(~x, ~k, t) dependent on position, momentum, and time. The related
reflective BC would then be


f |− (~x, ~k, t) = p(~x, ~k, t)f |+ (~x, ~k 0 , t) + 1 − p(~x, ~k, t) C ∗ (~x, t)σ ∗ (~x, t)M (~x, ~k)
(~x, ~k) ∈ ΓN − , and

(~x, ~k 0 ) ∈ ΓN + ,

(5.57)

where M (~x, ~k) is the equilibrium probability distribution (not necessarily a
Maxwellian) according to which the electrons diffusively reflect on the physical boundary. σ ∗ (~x, t) and C ∗ (~x, t) are the functions such that the zero flux

21

condition is satisfied pointwise at insulating boundaries
Z
Z
Z
0 = η(~x) ·
~v (~k)f d~k =
η(~x) · ~v (~k)f |+ d~k +
~
v ·η>0

Ω~k

Z

η · ~v f |+ d~k +

=

Z

~
v ·η>0

η(~x) · ~v (~k)f |− d~k

~
v ·η<0

h
i
η · ~v p(~x, ~k, t)f 0 |+ + (1 − p) C ∗ (~x, t)σ ∗ (~x, t)M (~x, ~k) d~k

~
v ·η<0

Z

η · ~v f |+ d~k +
~
v ·η>0
Z
− σ ∗ (~x, t) C ∗ (~x, t)

Z

η · ~v p(~x, ~k, t)f |+ (~x, ~k 0 , t)d~k

=

~
v ·η<0



|η · ~v | 1 − p(~x, ~k, t) M (~x, ~k)d~k .

~
v ·η<0

Therefore we conclude for this reflection case that
Z
Z
∗
~
σ {f |+ } (~x, t) =
|η · ~v |f |+ dk −
|η · ~v | p(~x, ~k, t)f |+ (~x, ~k 0 , t)d~k ,
~
v ·η>0

~
v ·η<0

(5.58)
C ∗ (~x, t) =

Z

−1

,
|η · ~v | 1 − p(~x, ~k, t) M (~x, ~k)d~k


(5.59)

~
v ·η<0

and then the full BC formula for the p(~x, ~k, t) reflection case is
f |− (~x, ~k, t) = p(~x, ~k, t)f |+ (~x, ~k 0 , t) +
i
hR
R


|η · ~v |f |+ d~k − ~v·η<0 |η · ~v | p(~x, ~k, t)f |+ (~x, ~k 0 , t)d~k
~
v ·η>0


1 − p(~x, ~k, t) M (~x, ~k)
R
|η · ~v | 1 − p(~x, ~k, t) M (~x, ~k)d~k .
~
v ·η<0
Remark: p(~x, ~k, t) can be any iid random variable in (~x, ~k, t).
5.3.1

Numerical Implementation

The numerical implementation of the general mixed reflection with specularity
parameter p(~k) is done in such a way that a numerical equivalent of the pointwise
zero flux condition is achieved.
The general mixed reflection boundary condition in our DG numerical scheme
is


Φh |− = Πh FM Φh |+
(5.60)

0
0 0
−w
= Πh p(w)Φ
~ h |+ (~x, w
~ , t) + (1 − p(w))C
~
σh {Φh |+ } (~x, t) e s(w) .
We will be using the notation
w
~ = (w, µ, ϕ),

dw
~ = dw dµ dϕ ,

w
~ 0 = (w, µ, π − ϕ).

(5.61)

The specific form of C 0 and σ 0 will be deduced from the numerical analogous
of the mixed reflection boundary condition. We want to satisfy numerically the

22

zero flux condition
Z
0 = η(~x) ·
~v (w)
~ Φh dw
~
(5.62)
Ωw
~
Z
Z
=
~v (w)
~ · η Φh |+ dw
~ +
~v (w)
~ · η Φh |− dw
~
~
v ·η>0
~
v ·η<0
Z
Z

0 0
=
~v · η Φh |+ dw
~ +
~v · η Πh p(w)Φ
~ 0h |+ + (1 − p(w))C
~
σh (~x, t)e−w s(w) dw
~
~
v ·η>0
~
v ·η<0
Z
Z
=
~v · η Φh |+ dw
~ −
|~v · η| Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} dw
~
~
v ·η>0
~
v ·η<0
Z

0 0
+
~v · η Πh (1 − p(w))C
~
σh (~x, t)e−w s(w) dw.
~
(5.63)
~
v ·η<0

In the space Vh of piecewise continuous polynomials which are tensor products
of polynomials of degree p in ~x and of degree q in w,
~ it holds that
Πh {f1 (~x)f2 (w)}
~ = Πh {f1 (~x)} Πh {f2 (w)}
~ ,
p,q

Vh = {v : v|Ωijkmn ∈ Q

p

(5.64)
q

(Ωijkmn ) = P (Xij ) ⊗ P (Kkmn )}.

Therefore, we have for our particular case that


R
−w
X
(1
−
p(
w))e
~
s(w)d
w
~

K
0 0

Πh (1 − p(w))C
~
σh (~x, t)e−w s(w) = C 0 σh0 (~x, t) 
χkmn kmn R
dw
~
Kkmn
k,m,n

Using this, our numerical pointwise zero flux condition is
Z

0

=

+
=

+
=
−
=
−

Z
~v · η Φh |+ dw
~ −
|~v · η| Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} dw
~
~
v ·η>0
~
v ·η<0


R
Z
−w
X
(1
−
p(
w))e
~
s(w)d
w
~
 dw
~v · η C 0 σh0 (~x, t) 
χkmn Kkmn R
~
d
w
~
~
v ·η<0
K
kmn
k,m,n
Z
Z
~v · η Φh |+ dw
~ −
|~v · η| Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} dw
~
~
v ·η>0
~
v ·η<0


R
Z
−w
X
(1
−
p(
w))e
~
s(w)d
w
~
K
 dw
C 0 σh0 (~x, t)
~v · η 
χkmn kmn R
~
d
w
~
~
v ·η<0
K
kmn
k,m,n
Z
Z
~v · η Φh |+ dw
~ −
|~v · η| Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} dw
~
~
v ·η>0
~
v ·η<0
R
Z
X
(1 − p(w))e
~ −w s(w)dw
~
0 0
C σh (~x, t)
χkmn
|~v · η| dw
~ Kkmn R
d
w
~
~
v ·η<0
Kkmn
k,m,n
Z
Z
~v · η Φh |+ dw
~ −
|~v · η| Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} dw
~
~
v ·η>0
~
v ·η<0
R
Z
X
(1 − p(w))e
~ −w s(w)dw
~
σh0 (~x, t) C 0
|~v · η| dw
~ Kkmn R
.
dw
~
Kkmn
Kkmn
k,m,n, ~
v ·η<0

23

We conclude then that we can achieve a numerical equivalent of the pointwise
zero flux condition by defining
Z
Z
σh0 {Φh |+ } (~x, t) =
~v · η Φh |+ dw
~−
|~v · η| Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} dw,
~
~
v ·η>0

0

(C {η})

−1

X

=

k,m,n, ~
v ·η<0

~
v ·η<0

R

Z

Kkmn

|~v · η| dw
~

(1 − p(w))e
~ −w s(w)dw
~
∆wk ∆µm ∆ϕn

Kkmn

. (5.65)

Therefore, the inflow BC in our DG numerical method is given by the expression
Φh |−

=
+

Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)}

Z
0
Πh (1 − p(w))C
~

~v · η Φh |+ dw
~

~
v ·η>0

Z

0

−



−w

|~v · η| Πh {p(w)Φ
~ h |+ (~x, w
~ , t)} dw
~ e


s(w) .

~
v ·η<0

The particular form of the coefficients defining the piecewise polynomial approximation Φh |− for the general mixed reflection BC is presented below for the
boundary y = Ly , since the calculations for the case of the boundary y = 0 are
analogous.
For the boundary yNy +1/2 = Ly , η · ~v ∝ +ŷ · ~g = g2 ∝ cos ϕ , which defines
the sign of g2 . Outflow cells have the index j = Ny . They are cells inside the
domain adjacent to the boundary. Inflow cells have the index j = Ny + 1. They
are ghost cells adjacent to the boundary. We have in our case that
Z
Z
σh0 =
g2 Φh |+ dw
~ −
|g2 | Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} dw
~
cos ϕ>0
n≤

=

cos ϕ<0

Nϕ
2

k,m,n

Nϕ
X2 Z

n>

X Z

g2 Φh |+ dw
~ −

Kkmn

|g2 | Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} dw
~.

Kkmn

k,m,n

If I = (i, Ny +1, k, m, n) (inflow), I 0 = (i, Ny , k, m, n0 ), n0 = Nϕ0 −n+1 (outflow),
the projection integrand is given by
n>Nϕ /2
0

Πh {p(w)Φ
~ h |+ (~x, w
~ , t)} =

X
I

R
χI



p(w)d
~ w
~
(x − xi )
TI 0 + XI 0
+ YI 0 (+1) .
∆xi /2
dw
~
kmn

kmn
R

The coefficients of σh0 are given below. We have now that I = (i, Ny , k, m, n),

24

I 0 = (i, Ny , k, m, n0 ),
n≤Np /2
0
σ 0 i,Ny

Z p
Z
Z p
w(1 + αK w)
dw
1 − µ2 dµ cos ϕdϕ
TI (t)
1 + 2αK w
m
n
k

X

=

n0 = Nϕ0 − n + 1 , so from the previous two formulas then

k,m,n
n>Np /2

−

x
σ 0 i,Ny

n≤Np /2

Z p

X

XI (t)
k

k,m,n
n>Np /2

XI 0 (t)
k

n≤Np /2
y

X

=

Z p
YI (t)
k

k,m,n
n>Np /2

X

−

Z p

Z p

X
k,m,n

σ 0 i,Ny

w(1 + αK w)
dw
1 + 2αK w

T (t)
k

−

R
Z p
Z
p(w)d
~ w
~
2
R
,
1 − µ dµ | cos ϕ|dϕ kmn
d
w
~
m
n
kmn

I0

k,m,n

=

w(1 + αK w)
dw
1 + 2αK w

Z p

X

YI 0 (t)
k

Z
cos ϕdϕ
n

Z p

1−

R

Z

µ2 dµ

| cos ϕ|dϕ

m

Z p

w(1 + αK w)
dw
1 + 2αK w

1 − µ2 dµ

m

w(1 + αK w)
dw
1 + 2αK w

w(1 + αK w)
dw
1 + 2αK w

Z p

k,m,n

(5.66)

n

1 − µ2 dµ

m

p(w)d
~ w
~
,
d
w
~
kmn

kmn
R

Z
cos ϕdϕ
n

R
Z p
Z
p(w)d
~ w
~
2
R
1 − µ dµ | cos ϕ|dϕ kmn
.
d
w
~
m
n
kmn

Since on one hand we have
−

Φh |Ly

=


0 0
Πh {p(w)Φ
~ h |+ (~x, w
~ 0 , t)} + Πh (1 − p(w))C
~
σh {Φh |+ } (~x, t) e−w s(w)
n>

Nϕ

R


p dw
~
(x − xi )
=
χi,Ny +1,kmn Rkmn
Ti,Ny ,k,m,n0 + Xi,Ny ,k,m,n0
+ Yi,Ny ,k,m,n0
∆xi /2
dw
~
kmn
ikmn
R
n>Nϕ /2
X
(1 − p(w))e
~ −w s(w)dw
~
R
+
χi,Ny +1,k,m,n kmn
×
dw
~
kmn
i,k,m,n


(x − xi )
y
0
x
+ σ 0 i,Ny (+1) ,
× C 0 σ 0 i,Ny + σ 0 i,Ny
∆xi /2
X2

and on the other hand
n>
−
Φh |yN +1/2
y

=

Nϕ

X2

i,k,m,n



(x − xi )
χi,Ny +1,k,m,n Ti,Ny +1,k,m,n + Xi,Ny +1,k,m,n
− Yi,Ny +1,k,m,n ,
∆xi /2

we conclude that the coefficients for Φh |− are
R
R
(1 − p(w))e
~ −w s(w)dw
~
p(w)d
~ w
~
0
0
0
kmn
R
+ C σ i,Ny kmn
,
Ti,Ny +1,k,m,n = TI 0 R
dw
~
dw
~
kmn
kmn
R
R
p(w)d
~ w
~
(1 − p(w))e
~ −w s(w)dw
~
0 0x
kmn
R
R
0
Xi,Ny +1,k,m,n = XI
+ C σ i,Ny kmn
,
d
w
~
d
w
~
kmn
kmn
R
R


p(w)d
~ w
~
(1 − p(w))e
~ −w s(w)dw
~
y
R
R
Yi,Ny +1,k,m,n = − YI 0 kmn
+ C 0 σ 0 i,Ny kmn
,
dw
~
dw
~
kmn
kmn
I0

=

(i, Ny , k, m, n0 ),

I = (i, Ny , k, m, n) .
25

(5.67)

6
6.1

Numerical Results
2D bulk silicon

We present results of numerical simulations for the case of n 2D bulk silicon
diode with an applied bias between the boundaries x = 0, Lx , and reflection
BC at the boundaries y = 0, Ly (Figs. 6.2). The required dimensionality in
momentum space is a 3D ~k(w, µ, ϕ). The specifics of our simulations are:
Initial Condition: Φ(w)|t=0 = Πh {N e−w s(w)}. Final Time: 1.0ps
Boundary Conditions (BC):
~k-space: Cut-off - at w = wmax , Φ is machine zero.
Only needed BC in (w, µ, ϕ): transport normal to the boundary analitically zero
at ’singular points’ boundaries:
At w = 0, g3 = 0. At µ = ±1, g4 = 0. At ϕ = 0, π, g5 = 0.
~x-space: Charge Neutrality at boundaries x = 0, x = 0.15µm.
Bias - Potential: V |x=0 = 0.5235 V, V |x=0.15µm = 1.5235 V.
Neumann BC for Potential at y = 0, Ly = 12nm: ∂y V |y=0, Ly = 0.
Reflection BC at y = 0, y = 12nm: Specular, Diffusive, Mixed Reflection with
constant specularity p = 0.5, and Mixed Reflection using a momentum dependent specularity p(~k) = exp(−4η 2 |k|2 sin2 ϕ), the nondimensional roughness rms
height coefficient being η = 0.5.
We observe an influence of the Diffusive and Mixed Reflection in macroscopic
observables. It is particularly noticeable in the kinetic moments. For example, the charge density slightly increases with diffusivity close to the reflecting
boundaries, and, due to mass conservation, alters the density profile over the
domain. Momentum & mean velocity increase with diffusive reflection over the
domain, while the energy is decreased by diffusive reflection over the domain.
There is a negligible difference in the electric field x component below its orders
of magnitude for the different reflection cases.

6.2

2D double gated MOSFET

We present as well the results of numerical simulations for the case of a 2D
double gated MOSFET device (Figs. 6.3). On one hand, the BC for the Poisson
Eq. for this device would be the Dirichlet BC Ψ = 0.5235 Volts at the source
x = 0, Ψ = 1.5235 Volts at the drain x = Lx , and Ψ = 1.06 Volts at the gates.
On the other hand, Homogeneous Neumann BC ∂n̂ Ψ = 0 are imposed at the
rest of the boundaries. Specular reflection is applied at the boundary y = 0
because the solution is symmetric with respect to y = 0 for our 2D double gate
MOSFET (Fig. 4.1). At the boundary y = Ly we apply specular, diffusive,
and mixed reflection BC, both with constant p = 0.5, and with a momentum
dependent p(~k) = exp(−4η 2 |k|2 sin2 ϕ) with roughness coefficient η = 0.5. We
use again the initial condition: Φ(w)|t=0 = Πh {ND (x, y)e−w s(w)}, running the
simulations up to the physical time of 1.0ps. We use again as well a cut-off BC
in the boundary of the momentum domain, so Φ is machine zero at w = wmax ,
and we apply charge neutrality BC at x = 0, x = 0.15µm.
We observe a quantitative difference in the kinetic moments and other observables between the different cases of reflective BC, with the physical quantities

26

being of the same order of magnitude. The electron density increases close to
the gates with diffusive reflection, and close to the center of the device, given
by the boundary y = 0, the density profile is greater for specular reflection.
The energy moment clearly decreases with diffusive reflection over the physical
domain. The momentum x-component for specular reflection is less than for
the other reflective cases. There is a difference in the profile of the electric field
x-component between the specular reflection and the other cases that include
diffusivity, increasing it with diffusive reflection close to the drain. The electric
field y-component increases with diffusive reflection close to the boundary y = 0
representing the center of the device. The electric potential is greater for the
cases including diffusive reflection than for the perfectly specular case.

6.3

Electrons reentering the 2D domain with reflective BC
in y and periodic BC in x: comparison of bulk silicon
with collisionless plasma

We consider in this case almost the same physical situation and parameters
for the previous section on the 2D bulk silicon, except that instead of using
the charge neutrality conditions we apply periodic boundary conditions in the
x-boundaries, simulating then that the electrons reenter the material on the
opposite x-boundary after the outflow exits the domain. We compare these
simulations with ones in which no collisions are considered, corresponding the
latter to the case of a collisionless plasma with reflective BC in y and periodic
BC in x. For both cases, bulk silicon with electron-phonon collisions and the
collisionless electron gas, we still apply an external potential such that V = 0
at x = 0 and V = 1Volt at x = Lx . This can be understood in the framework of periodic BC in x as a periodic sawtooth wave with period equal to the
length of the x-domain. We do this comparison in order to study the effect
of the reflective boundary conditions in y, with and without the influence of
the collisions over electrons, and we let the electrons re-enter the domain under
periodic boundary conditions in x, eliminating then the charge neutrality conditions in x and any possible effect due to the latter. Since due to the periodic
BC in x the electrons re-enter the domain after they exit it in outflow, the effect
of boundary conditions is exclusively related to the reflection in the transport
domain in the y-boundaries. For example, in Figs. 6.6 we present the plots of
Relative Mass vs Time (ps) for Specular, Diffusive, Mixed with constant and
momentum dependent specularity for different sets of simulations. The top figure is related to simulations for bulk silicon with charge neutrality conditions
on the non-reflecting boundaries, the middle figure is associated to simulations
for bulk silicon with periodic boundary conditions on the non-reflecting boundaries, and the bottom figure is related to the simulations for collisionless electron
transport with periodic boundary conditions on the non-reflecting boundaries.
The last two sets of simulations mentioned conserve the mass during all the
time, and these sets isolate the effect of reflection boundary conditions by using
periodic boundary conditions instead of charge neutrality conditions. The first
set associated to charge neutrality conditions in adition to reflection boundary conditions, however, have a slight increase in the relative mass of less than
0.5%. This slight increase then is associated only to the inclusion of charge neutrality conditions and a possible accumulation of numerical error due solely to it.

27

We notice in our comparison then the following effects of the collision operator
in comparison with the collisionless plasma case. As expected, the main effect
of collisions is to decrease the magnitude of the average energy, average velocity
and momentum (therefore the current) of electrons over the domain (Fig. 6.4).
The effect of collisions on the distribution of the electron density profile over the
domain is negligible. Regarding the isolated effects of the reflection boundary
conditions in the kinetic moments and other physical observables of interest by
considering the collisionless plasma with periodic BC in x, we notice, as earlier
in the section for bulk silicon, the slight increase of the density profile close to
the reflecting boundaries when adding diffusivity in the boundary conditions,
and by conservation of mass, a decrease of the density profile over the center
of the domain. The mean energy decreases over the position domain with the
inclusion of diffusive reflection BC, as well as the x components (which are the
dominant) of the momentum and velocity (Fig. 6.5). It is important to notice
this expected effect of the isolated reflection BC in the collisionless plasma
case, since for the case that includes electron-phonon collisions combined with
adding diffusive reflection BC gives actually an increase in the x components
of the momentum and velocity compared to the purely specular reflection case
(Fig. 6.4). The collisionless plasma with periodic BC in x and reflection BC
in y isolates the effect of the latter then and shows the expected behaviour of
a decrease in the mean energy, velocity and momentum x-compoments when
adding diffusivity in the reflection boundary conditions.

7

Conclusions

We have considered the mathematical and numerical modeling of Reflective
Boundary Conditions in 2D devices and their implementation in DG-BP schemes.
We have studied the specular, diffusive and mixed reflection BC on the boundaries of the position domain of the device. We developed a numerical equivalent
of the zero flux condition at the position domain boundaries for the case of a
more general mixed reflection with a momentum dependant specularity parameter p(~k). We compared the influence of these different reflection cases in the
computational prediction of moments after implementing numerical BC equivalent to the respective reflective BC, each one satisfying a mathematical zero
flux condition at insulating boundaries. There are effects due to the inclusion of
diffusive reflection boundary conditions over the moments and physical observables of the probability density function, whose influence is not only restricted
to the boundaries but actually to the whole domain. Particularly noticeable
effects of the inclusion of diffusivity in kinetic moments are the increase of the
density close to the reflecting boundary, the decrease of the mean energy over
the domain and, in the case when electron-phonon collisions for silicon are included, the increase of the x-components of the mean velocity and momentum
over the domain, whereas for the collisionless case, for which only the effects
of the reflection boundary conditions are considered (such as when electrons
are allowed to reenter the material via periodic boundary conditions in x), a
decrease in those x-components of mean velocity and momentum is observed,
as expected when adding diffusivity to the reflection boundary conditions.

28

To summarize, specular boundary conditions have the physical meaning of a reflection with a perfectly smooth surface with no roughness. Diffusive boundary
conditions are the opposite case, with the physical meaning of a rough surface
that completely diffuses the momentum. That means, for a fixed ~x at a time
t, since f |− (~x, ~k, t) = Cσ(~x, t) exp(−ε(~k)), the probability density function is
higher for momentum vectors with a lower energy band value, therefore the momentum with the highest f |− (~x, ~k, t) value, for that given ~x at time t, occurs at
~k = ~0, where the origin of the momentum space has been chosen as the position
of the local energy band valley. Hence, our physical interpretation of the diffusive reflection condition is that it diffuses the momentum giving a higher probability for lower magnitude momentum values, with highest probability density
value at ~k = ~0. The mixed reflection condition is a convex combination of both,
meaning that the reflection process is partially specular and partially diffuses
the momentum, with the probability of specularity potentially depending on
the momentum variable. Estimating which boundary condition is more physical, we believe mixed reflection case may be the most suitable, as no surface is
in practice perfectly specular, and the diffusive reflection is the case in the other
end of the spectrum that minimizes the total reflected momentum. Regarding
on how to understand the quantitative differences between the output for the
considered boundary conditions, we propose, in Section 6.3, a numerical study
of the different reflection conditions for collisionless electrons with periodic conditions on the other two boundaries. This case provides the best understanding
of the boundary condition role in the simulation as it isolates effects of the reflection conditions, since there is no dissipation from collision mechanisms, and
the electrons reenter the domain after exiting a periodic boundary. In fact, we
show in Figures 6.5 that the diffusivity in the boundary condition, as expected,
lowers the momentum in the main direction of transport of the electrons, which
is x̂ (the transport in the ŷ direction is negligible), and it also lowers the energy
average. Both of these quantitative differences are expected from the reflection
of the electrons with a rough boundary. Regarding the quantitative difference
between the density output, we observe that particles tend to stay closer to the
rough boundaries when increasing the degree of diffusivity in the boundary condition, since the diffusivity decreases the total reflected momentum as it is more
probable to have a reflected momentum with lower magnitudes. Therefore the
density profile increases for more diffusive conditions as it tends to accumulate
more particles in the boundary by lowering their momentum after reflection.
Future research will consider, for example, the inclusion of surface roughness
scattering mechanisms in the collision operator for our diffusive reflection problem in silicon devices. It will be related as well to the inclusion of diffusive
reflection BC with a DG-BP-EPM full energy band. More importantly, another
line of work of our interest for future research will be the more general case of a
p(~x, ~k) specular probability dependant on momentum and position as well, considering in addition to its mathematical aspects the related numerical issues and
the respective computational modelling, intending to use experimental values of
p(~x, ~k) as input for the simulations.

29

Acknowledgment
The authors’ research was partially supported by NSF grants NSF CHE-0934450,
NSF-RNMS DMS-1107465 and DMS 143064, and the ICES Moncrief Grand
Challenge Award. The computational work was partially performed by means
of TACC resources under project A-ti4. Support from the Institute of Computational Engineering and Sciences and the University of Texas Austin is gratefully
acknowledged.

References
[1] Y. Sone, Molecular Gas Dynamics: Theory, Techniques, and Applications,
Birkhauser (2007).
[2] A. Jüngel, Transport Equations for Semiconductors, Springer Verlag (2009).
[3] C. Cercignani, The Boltzmann Equation and Its Applications, SpringerVerlag, Appl. Math. Sc. 67, (1988).
[4] C. Cercignani, I. M. Gamba, C.D. Levermore, High Field Approximations
to a Boltzmann - Poisson System and Boundary Conditions in a Semiconductor, Appl. Math. Lett. 10, 4, 111-117 (1997).
[5] S. Soffer, Statistical Model for the size effect in Electrical Conduction, Journal of Applied Physics 38 1710 (1967).
[6] K. Fuchs, Proc. Cambridge Phil. Soc. 34 100 (1938).
[7] R. F. Greene, Boundary Conditions for Electron Distributions at Crystal
Surfaces, Physical Review 141 687 (1966).
[8] R. F. Greene, R. W. O’Donnell, Scattering of Conduction Electrons by Localized Surface Charges, Physical Review 147 599 (1966).
[9] V. D. Borman, S. Yu. Krylov, A. V. Chayanov, Theory of nonequilibrium
phenomena at a gas-solid interface, Sov. Phys. JETP 67 (10), 1988.
[10] Brull, Charrier, Mieussens, Gas-surface interaction and boundary conditions for the Boltzmann equation, Kinetic & Related Models (2014).
[11] Struchtrup, H. Maxwell boundary condition and velocity dependent accommodation coefficient, Phys. Fluids 25, 112001 (2013).
[12] P. Markowich, C. Ringhofer and C. Schmeiser, Semiconductor Equations,
Springer-Verlag, 1990.
[13] M. Lundstrom, Fundamentals of Carrier Transport, Cambridge University
Press, 2000.
[14] C. Jacoboni and P. Lugli, The Monte Carlo Method for Semiconductor
Device Simulation, Spring-Verlag: Wien-New York, 1989.
[15] 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.
30

[16] 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.
[17] 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.
[18] 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 WENO-Boltzmann schemes, Journal of Computational Electronics, 2 (2003) 375-380.
[19] 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.
[20] 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.
[21] 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.
[22] 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.
[23] 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) 274-280.
[24] 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.
[25] 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.
[26] 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.
[27] 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.
31

[28] 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.
[29] 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.
[30] Y. Cheng, I. M. Gamba, A. Majorana and C.W. Shu A discontinuous
Galerkin solver for Full-Band Boltzmann-Poisson Models, IWCE13 (13th
International Workshop on Computational Electronics) (2009).
[31] 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.
[32] 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.
[33] 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.
[34] 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.
[35] B. Cockburn and C.-W. Shu, The Runge-Kutta local projection P1discontinuous Galerkin finite element method for scalar conservation laws,
Mathematical Modelling and Numerical Analysis, 25 (1991) 337-361.
[36] 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.
[37] B. Cockburn and C.-W. Shu, Runge-Kutta discontinuous Galerkin methods for convection-dominated problems, Journal of Scientific Computing, 16
(2001) 173-261.
[38] 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.
[39] 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.
[40] 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.

32

Density (1/m3) vs. Position (x,y)

Energy Moment <e> vs. Position (x,y)

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

1.02e+21

0.3

1.015e+21

0.25

1.01e+21

0.2

1.005e+21

(eV)
0.15

1e+21

0.1

9.95e+20
9.9e+20

0.05
0.012

0

0.02

0.01
0.008
0.006
y (microns)
0.004

0.04

0.06
0.08
x (microns)

0.1

0.12

0

0.02

0.04

0.06
0.08
x (microns)

0.002
0.14

0.160

Momentum X-component (Ux: Current) vs. Position (x,y)

0.012
0.01
0.008
0.006
y (microns)
0.004
0.1

0.12

0.002
0.14

0.160

Momentum Y-component (Uy: Current) vs. Position (x,y)

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

0.0000014

4e-09
3e-09
2e-09
1e-09
0
-1e-09
-2e-09
-3e-09
-4e-09

0.0000013
0.0000012
0.0000011
0.0000010
0.0000009
0.0000008

0

0.02

0.012
0.01
0.008
0.006
y (microns)
0.004

0.04

0.06
0.08
x (microns)

0.1

0.12

0

0.02

0.04

0.06
0.08
x (microns)

0.002
0.14

0.160

0.012
0.01
0.008
0.006
y (microns)
0.004

Electric Field X-component (Ex) vs. Position (x,y)

0.1

0.12

0.002
0.14

0.160

Electric Field Y-component (Ey) vs. Position (x,y)

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

-66.661
-66.6615
-66.662
-66.6625
-66.663
-66.6635
-66.664
-66.6645
-66.665
-66.6655
-66.666

0.0004
0.0003
0.0002
0.0001
0
-0.0001
-0.0002
-0.0003

0

0.02

0.014
0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

0

0.02

0.04

0.06
0.08
x (microns)

0.160

0.014
0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.1

0.12

0.14

0.160

Electric Potential (V) vs. Position (x,y)
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

1.6
1.5
1.4
1.3
1.2
V (V) 1.1
1
0.9
0.8
0.7
0.6
0.5

0

0.02

0.014
0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

0.160

Mean energy e(eV ), Momentum
Figure 6.2:
Density ρ (m−3 ),
28 cm−2
Ux , Uy (10
),
Electric
Field
Components
Ex and Ey , and Potential
s
V (V olts) vs Position (x, y) in (µm) plot for Specular, Diffusive, Mixed p = 0.5
& Mixed p(~k) = exp(−4η 2 |k|2 sin2 ϕ), η = 0.5 Reflection for 2D bulk silicon.

33

Density (1/m3) vs. Position (x,y)

Energy Moment <e> vs. Position (x,y)

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

6.5e+24

0.4
0.35
0.3
0.25
(eV)
0.2
0.15
0.1
0.05

6e+24
5.5e+24
5e+24
4.5e+24
4e+24
3.5e+24

0

0.02

0.012
0.01
0.008
0.006
y (microns)
0.004

0.04

0.06
0.08
x (microns)

0.1

0.12

0

0.02

0.04

0.06
0.08
x (microns)

0.002
0.14

0.160

Momentum X-component (Ux: Current) vs. Position (x,y)

0.012
0.01
0.008
0.006
y (microns)
0.004
0.1

0.12

0.002
0.14

0.160

Momentum Y-component (Uy: Current) vs. Position (x,y)

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

0.0055

0.0002
0.00015
0.0001
5e-05
0
-5e-05
-0.0001
-0.00015
-0.0002
-0.00025

0.005
0.0045
0.004
0.0035
0.003
0.0025

0

0.02

0.012
0.01
0.008
0.006
y (microns)
0.004

0.04

0.06
0.08
x (microns)

0.1

0.12

0

0.02

0.04

0.06
0.08
x (microns)

0.002
0.14

0.160

Electric Field X-component (Ex, 105 V/m) vs. Position (x,y)

0.012
0.01
0.008
0.006
y (microns)
0.004
0.1

0.12

0.002
0.14

0.160

Electric Field Y-component (Ey, 105 V/m) vs. Position (x,y)

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

-20
-30
-40
-50
-60
-70
-80
-90
-100
-110
-120
-130

50
0
-50
-100
-150
-200
-250

0

0.02

0.014
0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

0

0.02

0.04

0.06
0.08
x (microns)

0.160

0.014
0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.1

0.12

0.14

0.160

Electric Potential (V) vs. Position (x,y)
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

V

1.6
1.5
1.4
1.3
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5

0

0.02

0.014
0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

0.160

Figure 6.3:
Density ρ (m−3 ),
Mean energy e(eV ), Momentum
28 cm−2
Ux , Uy (10
),
Electric
Field
Components
Ex and Ey , and Potential
s
V (V olts) vs Position (x, y) in (µm) plot for Specular, Diffusive, Mixed p = 0.5
& Mixed p(~k) = exp(−4η 2 |k|2 sin2 ϕ), η = 0.5 Reflection for a 2D double gated
MOSFET.

34

Density (1/m3) vs. Position (x,y)

Energy Moment <e> vs. Position (x,y)
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

1.02e+21

0.3

1.015e+21

0.25

1.01e+21

(eV)

0.2

1.005e+21
0.15

1e+21

0.1

9.95e+20
9.9e+20

0.05
0.012

0

0.02

0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

0

0.02

0.04

0.06
0.08
x (microns)

0.160

Momentum X-component (Ux: Current) vs. Position (x,y)

0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.1

0.12

0.14

0.160

Momentum Y-component (Uy: Current) vs. Position (x,y)

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

1.4e-06

2e-13
1.5e-13
1e-13
5e-14
0
-5e-14
-1e-13
-1.5e-13
-2e-13

1.3e-06
1.2e-06
1.1e-06
1e-06
9e-07
8e-07

0

0.02

0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

0

0.02

0.04

0.06
0.08
x (microns)

0.160

0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

Velocity X-component <Vx> vs. Position (x,y)

0.1

0.12

0.14

0.160

Electric Potential (V) vs. Position (x,y)
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

0.14
0.13
0.12

V (V)

0.11
0.1
0.09
0.08

0

0.02

0.012
0.01
0.008
0.006
y
(microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

1.6
1.5
1.4
1.3
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5

0

0.02

0.06
0.08
x (microns)

0.160

0.014
0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.1

0.12

0.14

0.160

Figure 6.4:
Density ρ (m−3 ),
Mean energy e(eV ), Momentum
28 cm−2
Ux , Uy (10
s ), Average Velocity Component Vx , and Potential V (V olts)
vs Position (x, y) in (µm) plot for Specular, Diffusive, Mixed p = 0.5 & Mixed
p(~k) = exp(−4η 2 |k|2 sin2 ϕ), η = 0.5 Reflection for electrons in 2D bulk silicon
with reflective BC in y and periodic BC in x.

35

Density (1/m3) vs. Position (x,y)

Energy Moment <e> vs. Position (x,y)
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

1.03e+21
1.025e+21
1.02e+21
1.015e+21
1.01e+21
1.005e+21
1e+21
9.95e+20
9.9e+20

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

(eV)

0

0.02

0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

0.14

0

0.02

0.04

0.06
0.08
x (microns)

0.160

Momentum X-component (Ux: Current) vs. Position (x,y)

0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.1

0.12

0.14

0.160

Momentum Y-component (Uy: Current) vs. Position (x,y)

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

6.5e-06
6e-06
5.5e-06
5e-06
4.5e-06
4e-06
3.5e-06
3e-06
2.5e-06
2e-06
1.5e-06
1e-06

8e-11
6e-11
4e-11
2e-11
0
-2e-11
-4e-11
-6e-11
-8e-11

0

0.02

0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

0

0.02

0.04

0.06
0.08
x (microns)

0.160

0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

Velocity X-component <Vx> vs. Position (x,y)

0.1

0.12

0.14

0.160

Electric Potential (V) vs. Position (x,y)
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

0.65
0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1

Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

V (V)

0

0.02

0.012
0.01
0.008
0.006
y
(microns)
0.004
0.002

0.04

0.06
0.08
x (microns)

0.1

0.12

0.14

1.6
1.5
1.4
1.3
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5

0

0.02

0.06
0.08
x (microns)

0.160

0.014
0.012
0.01
0.008
0.006
y (microns)
0.004
0.002

0.04

0.1

0.12

0.14

0.160

Figure 6.5:
Density ρ (m−3 ),
Mean energy e(eV ), Momentum
28 cm−2
Ux , Uy (10
s ), Average Velocity Component Vx , and Potential V (V olts)
vs Position (x, y) in (µm) plot for Specular, Diffusive, Mixed p = 0.5 & Mixed
p(~k) = exp(−4η 2 |k|2 sin2 ϕ), η = 0.5 Reflection for 2D collisionless electrons
with reflective BC in y and periodic BC in x.

36

Relative Mass at each time tl
1.01
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

M(t)/M(0)

1.005

1

0.995

0.99
0

0.2

0.4

0.6
t (ps)

0.8

1

1.2

Relative Mass at each time tl
1.01
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

M(t)/M(0)

1.005

1

0.995

0.99
0

0.2

0.4

0.6
t (ps)

0.8

1

1.2

Relative Mass at each time tl
1.01
Specular
Diffusive
Mixed Ct p=0.5
p(k) Soffer eta=0.5

M(t)/M(0)

1.005

1

0.995

0.99
0

0.2

0.4

0.6
t (ps)

0.8

1

1.2

Figure 6.6: Relative Mass vs Time (ps) plot for Specular, Diffusive, Mixed
p = 0.5 & Mixed p(~k) = exp(−4η 2 |k|2 sin2 ϕ), η = 0.5 Reflection. The figure on
top is related to the simulations for bulk silicon with charge neutrality conditions on the non-reflecting boundaries. The figure in the middle is associated
to the simulations for bulk silicon with periodic boundary conditions on the
non-reflecting boundaries. The bottom figure is related to the simulations for
collisionless electron transport with periodic boundary conditions on the nonreflecting boundaries. The figures show the conservation of mass when isolating
the effect of reflection boundary conditions in the simulations, as the slight increase in the relative mass of less than 0.5% is associated only to simulations
that also include charge neutrality conditions and possibly an accumulation of
numerical error due solely to it.
37

