# Numerical analysis of a finite volume scheme for charge transport in
perovskite solar cells

Dilara Abdel Weierstrass Institute (WIAS), Mohrenstr. 39, 10117 Berlin,
Germany Claire Chainais-Hillairet Univ. Lille, CNRS, Inria, UMR 8524 -
Laboratoire Paul Painlevé, F-59000 Lille, France Patricio
Farrell11footnotemark: 1 Maxime Herda Inria, Univ. Lille, CNRS, UMR 8524 -
Laboratoire Paul Painlevé, F-59000 Lille, France

###### Abstract

In this paper, we consider a drift-diffusion charge transport model for
perovskite solar cells, where electrons and holes may diffuse linearly
(Boltzmann approximation) or nonlinearly (e.g. due to Fermi-Dirac statistics).
To incorporate volume exclusion effects, we rely on the Fermi-Dirac integral
of order $-1$ when modeling moving anionic vacancies within the perovskite
layer which is sandwiched between electron and hole transport layers. After
non-dimensionalization, we first prove a continuous entropy-dissipation
inequality for the model. Then, we formulate a corresponding two-point flux
finite volume scheme on Voronoi meshes and show an analogous discrete entropy-
dissipation inequality. This inequality helps us to show the existence of a
discrete solution of the nonlinear discrete system with the help of a
corollary of Brouwer’s fixed point theorem and the minimization of a convex
functional. Finally, we verify our theoretically proven properties
numerically, simulate a realistic device setup and show exponential decay in
time with respect to the $L^{2}$ error as well as a physically and
analytically meaningful relative entropy.

###### Contents

  1. 1 Introduction
  2. 2 Charge transport model for perovskite solar cells
    1. 2.1 Boundary conditions
    2. 2.2 Statistics functions
    3. 2.3 Thermodynamic free energy
    4. 2.4 Non-dimensionalization of the model
  3. 3 Continuous entropy-dissipation inequality
    1. 3.1 Entropy functions
    2. 3.2 Proof of the entropy-dissipation inequality
  4. 4 Discrete version of charge transport model
    1. 4.1 Definition of discretization mesh
    2. 4.2 Finite volume discretization
    3. 4.3 Discrete entropy-dissipation inequality
  5. 5 Existence of a discrete solution
    1. 5.1 Existence of electric potential
    2. 5.2 A priori estimates
    3. 5.3 Existence of quasi Fermi potentials
  6. 6 Numerical experiments
    1. 6.1 Verifying the properties of the scheme
    2. 6.2 PSC simulation setup
  7. 7 Conclusion and outlook
  8. A Estimates on statistics and entropy functions
    1. A.1 Proof of Lemma 3.1
    2. A.2 Proof of Lemma 3.2
    3. A.3 Boltzmann and Fermi-Dirac entropy functions
  9. B A technical result

##  1 Introduction

In recent years, Perovskite Solar Cells (PSCs) have become one of the fastest
growing photovoltaic technologies [30, 38]. Two advantages of PSCs stand out:
On the one hand, certain architectures have significantly lower production
costs than conventional solar cells. On the other hand, silicon-perovskite
tandem cells have become more efficient than classical single-junction silicon
solar cells. However, the commercialization of PSCs is still in its early
stages and several challenges need to be overcome, the most prominent being
the relatively fast degradation of these devices. Apart from solar cells,
perovskite materials also show promise for use in LEDs, photodetectors and
memristors.

There is well-established mathematical literature concerning drift-diffusion
mathematical models to describe charge transport in classical or organic
semiconductors and similar physical systems (see for instance and non-
exhaustively [20, 21, 32, 33, 35, 39]). In the perovskite material, a major
difference is that apart from electrons and holes, ion migration plays a
fundamental role. Therefore, to correctly reflect the physical behavior,
perovskite models must contain additional freely moving ion species. In PSCs,
the various charge carrier species live in different parts of the domain,
evolve on different time scales and obey different diffusion laws. Moreover,
one must take into account the photogeneration effect. Because of the
differences and new features compared to classical semiconductors, there is a
need for new mathematical and numerical modeling and analysis for these
devices. Several models have been proposed recently in the literature – with
the exception of [3], nearly all of them are in one dimension [10, 15].

In this paper, we consider a three-layer drift-diffusion perovskite charge
transport model for three different charge carriers. Within the entire device
electrons and holes may diffuse linearly (Boltzmann approximation) or
nonlinearly. The nonlinear diffusion may be governed, for example, by Fermi-
Dirac statistics [22] but our model allows even more general statistical
relationships between densities and potentials. In the middle perovskite
region, we include anion vacancies. It has been demonstrated that the ion
migration despite being considerably slow cannot be neglected [10]. To
correctly incorporate volume exclusion effects, we rely on Fermi-Dirac
statistics of order $-1$. All three drift-diffusion equations are coupled
self-consistently to a nonlinear Poisson equation.

Given that perovskite models are relatively new, the arising Partial
Differential Equations (PDE) model has not been studied yet mathematically. In
comparison, classical drift-diffusion models have been studied in detail [20,
21, 32, 33, 35]. For drift-diffusion models, an essential a priori estimate is
based on the evolution of the physical free energy or a related functional. It
allows to study the well-posedness of the equations as well as the asymptotic
behavior of its solution [8, 20, 35]. The techniques relying on a well-chosen
physically relevant Lyapunov functional have been used for many systems of
dissipative PDEs and are usually referred to as entropy methods (see [26] and
references therein), which is why in the following we will use the term
_entropy_ instead of _free energy_.

The design of numerical schemes for drift-diffusion models is also an mature
but still very active field of research (see for instance [9, 11, 25, 28, 29,
34, 35, 36, 40]). In order to ensure the quality of the numerical simulation
and the stability of the numerical method, efforts have been made towards the
design of structure preserving schemes [6, 12, 24, 31, 34]. Their aim is to
preserve physical features of the original model such as the decay of free-
energy or non-negativity of solutions. Because of the stiffness in drift-
diffusion models arising from small parameters such as the Debye length, fully
implicit in time numerical methods are usually preferred. They yield
robustness of the scheme as well as asymptotic preserving properties [7, 12].
Finally, the treatment of nonlinear diffusion to handle general statistics
function has also been investigated [5, 19, 25].

The main goal of this paper is the analysis of an implicit in time two-point
flux approximation (TPFA) finite volume scheme for the perovskite model. The
model and the scheme originate from [3]. The scheme relies on the the excess
chemical potential flux scheme which appears to be used for the first time in
[40] and was later numerically analyzed in [11, 23] and compared in [1, 28].

The main tool for our analysis is an entropy-dissipation inequality for the
perovskite model. After non-dimensionalizing the perovskite model, we
establish such an inequality in the continuous setting in Theorem 3.3. Then,
we adapt the arguments to show in Theorem 4.4 that the discrete counterpart of
this inequality also holds for a solution of the implicit finite volume
scheme. This _a priori_ estimate on the scheme allows us to prove the
existence of a discrete solution at each time step in Theorem 5.5. The proof
relies on a corollary of Brouwer’s fixed point theorem for the quasi Fermi
potentials, coupled with the minimization of a convex functional for the
electric potential.

We illustrate and complement our theoretical results with numerical
experiments. We investigate the convergence in space and witness second order
accuracy, as expected. By introducing a relative free energy with respect to
the steady solution, we illustrate the long time behavior of transient
solutions and their convergence towards steady state solutions exponentially
fast in time. Finally, we investigate the large time behavior of the
perovskite model at a constant applied voltage which physically corresponds to
investigating the influence of preconditioning a PSC before current-voltage
measurements.

The remainder of the paper is organized as follows: In Section 2, we introduce
the original and the non-dimensionalized perovskite model for general
statistics functions and the underlying free energy. In Section 3, we then
prove a continuous entropy-dissipation inequality for the perovskite model. In
Section 4, we present the corresponding finite volume scheme, for which we
prove a discrete entropy-dissipation inequality. This inequality will allow us
to deduce the existence of a discrete solution in Section 5. Finally, in
Section 6, we corroborate and complement our theoretical observations
numerically before we conclude and discuss future research in Section 7.

##  2 Charge transport model for perovskite solar cells

Let $\mathbf{\Omega}\subseteq\mathbb{R}^{d},d\leq 3$, be an open, connected
and bounded spatial domain, which is partitioned into three pairwise disjoint,
open subdomains
$\mathbf{\Omega}=\left(\cup_{k}\overline{\mathbf{\Omega}}_{k}\right)^{\circ},$
$k\in\\{\text{intr},\text{HTL},\text{ETL}\\}$. Here,
$\mathbf{\Omega}_{\text{intr}}$ refers to the intrinsic perovskite region and
$\mathbf{\Omega}_{\text{HTL}},\mathbf{\Omega}_{\text{ETL}}$ to the doped
electron and hole transport layers, respectively. We denote the interface
between the transport layers and the perovskite layer by
$\mathbf{\Sigma}_{\text{HTL}}=\partial\mathbf{\Omega}_{\text{HTL}}\cap\partial\mathbf{\Omega}_{\text{intr}}$
and
$\mathbf{\Sigma}_{\text{ETL}}=\partial\mathbf{\Omega}_{\text{ETL}}\cap\partial\mathbf{\Omega}_{\text{intr}}$.
The boundaries of the hole and electron transport layer do not intersect,
$\partial\mathbf{\Omega}_{\text{HTL}}\cap\partial\mathbf{\Omega}_{\text{ETL}}=\emptyset$,
see Figure 2.1 for a potential device geometry. The unknowns of the charge
transport model are given by the electric potential $\psi(\mathbf{x},t)$ and
the quasi Fermi potentials (frequently called electrochemical potentials) of
moving charge carriers, denoted by $\varphi_{\alpha}(\mathbf{x},t)$,
$\alpha\in\\{n,p,a\\}$. Here, the quantities $n$, $p$ and $a$ refer to the
electrons, holes and anion vacancies. Unlike quasi Fermi potentials for
electrons and holes $\varphi_{n},\varphi_{p}$, which are defined for any
$\mathbf{x}\in\mathbf{\Omega}$, the quasi Fermi potential of anion vacancies
$\varphi_{a}$ is defined only in the intrinsic domain
$\mathbf{\Omega}_{\text{intr}}$. The density of a charge carrier is denoted by
$n_{\alpha}$, $\alpha\in\\{n,p,a\\}$. We examine the model for the charge
transport in PSCs formulated in [3], where for $t\geq 0$ the mass balances are
given by

$\displaystyle z_{n}q\partial_{t}n_{n}+\nabla\cdot\mathbf{j}_{n}$
$\displaystyle=z_{n}q\Bigl{(}G(\mathbf{x})-R(n_{n},n_{p})\Bigr{)},$
$\displaystyle{\bf x}\in\mathbf{\Omega},\ t\geq 0,$ (2.1a) $\displaystyle
z_{p}q\partial_{t}n_{p}+\nabla\cdot\mathbf{j}_{p}$
$\displaystyle=z_{p}q\Bigl{(}G(\mathbf{x})-R(n_{n},n_{p})\Bigr{)},$
$\displaystyle{\bf x}\in\mathbf{\Omega},\ t\geq 0,$ (2.1b) $\displaystyle
z_{a}q\partial_{t}n_{a}+\nabla\cdot\mathbf{j}_{a}$ $\displaystyle=0,$
$\displaystyle{\bf x}\in\mathbf{\Omega}_{\text{intr}},\ t\geq 0,$ (2.1c)

self-consistently coupled with the nonlinear and region-wise defined Poisson
equation

$-\nabla\cdot(\varepsilon_{s}\nabla\psi)=\left\\{{\begin{aligned}
&q\Bigl{(}z_{n}n_{n}+z_{p}n_{p}+C({\bf x})\Bigr{)},&&{\bf
x}\in{\mathbf{\Omega}}_{\text{HTL}}\cup{\mathbf{\Omega}}_{\text{ETL}},\ t\geq
0,\\\ &q\Bigl{(}z_{n}n_{n}+z_{p}n_{p}+z_{a}n_{a}+C({\bf x})\Bigr{)},&&{\bf
x}\in{\mathbf{\Omega}}_{\text{intr}},\ t\geq 0.\end{aligned}}\right.$ (2.2)

The charge numbers of the three moving charge carriers are given by $z_{n}$,
$z_{p}$ and $z_{a}$, the elementary charge is denoted by $q$ and
$\varepsilon_{s}$ refers to the region-dependent dielectric permittivity.
Throughout this paper, we assume the standard charge numbers $z_{n}=-1$ and
$z_{p}=1$ for electrons and holes. For the anionic vacancies, we will simply
assume $z_{a}>0$ and later for the numerical experiments we set $z_{a}=1$.
Note that all of the following computations can be extended to include
movement of cation vacancies, i.e. movement of negative ion vacancy charge
numbers. However, this case seems to be of less physical interest for PSCs.
The charge carrier densities are linked to the set of unknowns
$(\psi,\varphi_{n},\varphi_{p},\varphi_{a})$ by the state equations [18]

$n_{\alpha}=N_{\alpha}\mathcal{F}_{\alpha}\Bigl{(}\eta_{\alpha}(\varphi_{\alpha},\psi)\Bigr{)},\quad\eta_{\alpha}=z_{\alpha}\frac{q(\varphi_{\alpha}-\psi)+E_{\alpha}}{k_{B}T},\quad\alpha\in\\{n,p,a\\},$
(2.3)

where $\mathcal{F}_{\alpha}$ is called statistics function which will be
discussed in Section 2.2. The quantities $N_{n},N_{p}$ denote the effective
conduction and valence density of states, whereas $N_{a}$ is given by the
maximal vacancy concentration. The argument $\eta_{\alpha}$ of the statistics
function is called chemical potential. It depends on the band-edge energies
$E_{\alpha}$, a constant temperature $T$ and the Boltzmann constant $k_{B}$.
Moreover, the electric currents $\mathbf{j}_{\alpha}$ for each species are
given by

$\mathbf{j}_{\alpha}=-z_{\alpha}^{2}\mu_{\alpha}n_{\alpha}\nabla\varphi_{\alpha},\quad\alpha\in\\{n,p,a\\},$
(2.4)

with the carrier mobility $\mu_{\alpha}$. Concerning the right-hand side of
the continuity equations (2.1a), (2.1b) and Poisson equation (2.2) we assume
that the doping profile $C$ is bounded, i.e. $C\in
L^{\infty}(\mathbf{\Omega})$ and that the photogeneration rate satisfies
$0\leq G\in L^{\infty}(\mathbf{\Omega})$. In other words the carrier dependent
doping profile and the photogeneration rate are constant in time. It is common
to assume a Beer-Lambert generation profile, describing an exponential decay
in the $z$ direction, see Figure 2.1,

$G(\mathbf{x})=F_{\text{ph}}\alpha_{g}\exp(-\alpha_{g}z),\quad\mathbf{x}=(x,y,z)^{T},$

where $F_{\text{ph}}$ denotes the incident photon flux and $\alpha_{g}$ the
material absorption coefficient. Lastly, the recombination rate $R$ is of the
form [18]

$\displaystyle
R(n_{n},n_{p})=r(n_{n},n_{p})n_{n}n_{p}\Bigl{(}1-\exp\left(\varphi_{n}-\varphi_{p}\right)\Bigr{)},\quad
r(n_{n},n_{p})=\sum_{r}r_{r}(n_{n},n_{p}),$

where $r_{r}(n_{n},n_{p})\geq 0$ is given by the sum of all present
recombination processes, which, for PSCs, are radiative and trap-assisted
Shockley-Read-Hall recombination

$\displaystyle r_{\text{rad}}=r_{0}\quad\text{and}\quad
r_{\text{SRH}}=\frac{1}{\tau_{p}(n_{n}+n_{n,\tau})+\tau_{p}(n_{p}+n_{p,\tau})},$

where $r_{0}$ is a constant rate coefficient, $\tau_{n},\tau_{p}$ are the
carrier life times and $n_{n,\tau},n_{p,\tau}$ the reference carrier
densities. Furthermore, we supply the system (2.1) with initial conditions for
$t=0$

$\displaystyle\varphi_{n}(\mathbf{x},0)=\varphi_{n}^{0}(\mathbf{x}),\quad$
$\displaystyle\quad\varphi_{p}(\mathbf{x},0)=\varphi_{p}^{0}(\mathbf{x})\quad$
$\displaystyle\text{for\;}\mathbf{x}\in\boldsymbol{\Omega},\hphantom{xx}$
(2.5a) $\displaystyle\varphi_{a}(\mathbf{x},0)$
$\displaystyle=\varphi_{a}^{0}(\mathbf{x}),\quad$
$\displaystyle\text{for\;}\mathbf{x}\in\boldsymbol{\Omega}_{\text{intr}},$
(2.5b)

where we assume $\varphi_{n}^{0},\varphi_{p}^{0}\in
L^{\infty}(\mathbf{\Omega})$ and $\varphi_{a}^{0}\in
L^{\infty}(\mathbf{\Omega}_{\text{intr}})$. Correspondingly, we define initial
densities
$n_{\alpha}^{0}(\mathbf{x})=N_{\alpha}\mathcal{F}_{\alpha}(\eta_{\alpha}(\varphi_{\alpha}^{0},\psi(\mathbf{x},0)))$.

###  2.1 Boundary conditions

The outer boundary of $\mathbf{\Omega}$ is decomposed into two ohmic contacts
modeled by Dirichlet conditions $\mathbf{\Gamma}^{D}$ and an isolated
interface $\mathbf{\Gamma}^{N}$, where we impose no flux Neumann boundary
conditions. We assume
$\mathbf{\Gamma}^{D}\cap\overline{\mathbf{\Omega}}_{\text{intr}}=\emptyset$,
i.e. the ohmic contacts are solely located at the outer boundary of transport
layers. More precisely, let the Dirichlet values $\psi^{D},\varphi^{D}\in
W^{1,\infty}(\mathbf{\Omega})$ be given. Then, the outer boundary conditions
are modeled via

$\displaystyle\psi(\mathbf{x},t)=\psi^{D}(\mathbf{x}),\;\;\varphi_{n}(\mathbf{x},t)=\varphi_{p}(\mathbf{x},t)$
$\displaystyle=\varphi^{D}(\mathbf{x}),\quad~{}$
$\displaystyle\mathbf{x}\in\mathbf{\Gamma}^{D},\ t\geq 0,$ (2.6a)
$\displaystyle~{}~{}\nabla\psi(\mathbf{x},t)\cdot\boldsymbol{\nu}(\mathbf{x})=\mathbf{j}_{n}(\mathbf{x},t)\cdot\boldsymbol{\nu}(\mathbf{x})=\mathbf{j}_{p}(\mathbf{x},t)\cdot\boldsymbol{\nu}(\mathbf{x})$
$\displaystyle=0,\quad~{}$ $\displaystyle\mathbf{x}\in\mathbf{\Gamma}^{N},\
t\geq 0,$ (2.6b)

where $\boldsymbol{\nu}$ is the outward pointing unit normal to
$\mathbf{\Gamma}^{N}$. Note that the same Dirichlet value $\varphi^{D}$ is
imposed on both quasi Fermi potentials. Concerning the anion vacancies, we
impose no flux Neumann boundary conditions on the whole intrinsic boundary,
namely

$\mathbf{j}_{a}(\mathbf{x},t)\cdot\boldsymbol{\nu}_{\text{intr}}(\mathbf{x})=0,\quad\mathbf{x}\in\partial\mathbf{\Omega}_{\text{intr}},\
t\geq 0,$ (2.7)

where $\boldsymbol{\nu}_{\text{intr}}$ is the outward pointing unit normal to
$\partial\mathbf{\Omega}_{\text{intr}}$.

Figure 2.1: A two-dimensional three-layer device with the relevant potentials
stated per subdomain.

The traces of the potentials $\psi$, $\varphi_{n}$ and $\varphi_{p}$ coincide
on both sides of the internal boundaries $\mathbf{\Sigma}_{\text{HTL}}$ and
$\mathbf{\Sigma}_{\text{ETL}}$. Moreover, the corresponding fluxes are also
continuous across internal boundaries. More precisely, for $t\geq 0$

$\displaystyle\Bigl{(}\varepsilon_{s}\nabla\psi(\mathbf{x},t;k)-\varepsilon_{s}\nabla\psi(\mathbf{x},t;\text{intr})\Bigr{)}\cdot\boldsymbol{\nu}_{\text{intr}}(\mathbf{x})$
$\displaystyle=0,\quad\mathbf{x}\in\mathbf{\Sigma}_{k},\
k\in\\{\text{HTL},\text{ETL}\\},$ (2.8a)
$\displaystyle~{}\Bigl{(}\mathbf{j}_{n}(\mathbf{x},t;k)-\mathbf{j}_{n}(\mathbf{x},t;\text{intr})\Bigr{)}\cdot\boldsymbol{\nu}_{\text{intr}}(\mathbf{x})$
$\displaystyle=0,\quad\mathbf{x}\in\mathbf{\Sigma}_{k},\
k\in\\{\text{HTL},\text{ETL}\\},$ (2.8b)
$\displaystyle~{}\Bigl{(}\mathbf{j}_{p}(\mathbf{x},t;k)-\mathbf{j}_{p}(\mathbf{x},t;\text{intr})\Bigr{)}\cdot\boldsymbol{\nu}_{\text{intr}}(\mathbf{x})$
$\displaystyle=0,\quad\mathbf{x}\in\mathbf{\Sigma}_{k},\
k\in\\{\text{HTL},\text{ETL}\\}.$ (2.8c)

Here, we use the notation that for any function $f$ the expression
$f(\mathbf{x},t;k)$ denotes the trace of $f$, restricted onto
$\mathbf{\Omega}_{k}$, $k\in\\{\text{HTL},\text{ETL},\text{intr}\\}$,
evaluated at the respective interface between transport and perovskite layer.

######  Remark 2.1 (Conservation of mass for anion vacancies).

Observe that by integrating (2.1c) over $\mathbf{\Omega}_{\text{intr}}$ and
using the Neumann boundary conditions (2.7) the total mass of anionic
vacancies is conserved, namely

$\int_{\mathbf{\Omega}_{\text{intr}}}n_{a}(\mathbf{x},t)\,d\mathbf{x}=\int_{\mathbf{\Omega}_{\text{intr}}}n_{a}^{0}(\mathbf{x})\,d\mathbf{x},\quad\text{for
all}\;t\geq 0.$

An equivalent condition does not hold for electron and hole densities due to
the boundary conditions and the recombination/generation terms.

###  2.2 Statistics functions

Lastly, we need to discuss the choice of statistics functions
$\mathcal{F}_{\alpha}$ in (2.3) depending on the charge carrier species
$\alpha\in\\{n,p,a\\}$. Our results will hold under general abstract
assumptions on these functions. However, we also provide specific examples
below.

#### Electric Charge Carriers $\boldsymbol{\alpha\in\\{n,p\\}}$.

For electrons and holes $\alpha\in\\{n,p\\}$ we assume

$\left\\{\begin{aligned}
&\mathcal{F}_{n},\mathcal{F}_{p}:\mathbb{R}\rightarrow(0,\infty)\text{ are
}C^{1}\text{- diffeomorphisms};\\\\[5.0pt]
&0<\mathcal{F}_{\alpha}^{\prime}(\eta)\leq\mathcal{F}_{\alpha}(\eta)\leq\exp(\eta),\quad\eta\in\mathbb{R},\
\alpha\in\\{n,p\\}.\end{aligned}\right.$ (H1)

There are two important statistics functions, satisfying assumption (H1),
which are commonly used for modeling electric charge transport in PSCs. The
first one is given by the Fermi-Dirac integral of order $\mathit{1/2}$ defined
as

$F_{1/2}(\eta)=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}\frac{\xi^{1/2}}{\exp(\xi-\eta)+1}\>\mathrm{d}\xi,$
(2.9)

which is fundamental in the simulation of inorganic three-dimensional
semiconductors [18, 37]. The function $F_{1/2}$ behaves like $\eta^{3/2}$ when
the chemical potential $\eta$ tends to $+\infty$, namely in the large density
limit. In the low density limit, when the chemical potential $\eta$ tends to
$-\infty$, it behaves like the Boltzmann statistics function

$F_{B}(\eta)=\exp(\eta),$ (2.10)

which is another important statistics functions for electrons and holes.
Observe that the choice $\mathcal{F}_{n}=\mathcal{F}_{p}=F_{1/2}$ leads to
nonlinear diffusion in the electric currents (2.4), whereas
$\mathcal{F}_{n}=\mathcal{F}_{p}=F_{B}$ yields linear diffusion.

#### Ionic Charge Carriers $\boldsymbol{\alpha=a}$.

We assume that the statistics function for ionic charge carriers satisfies the
following assumption

$\left\\{\begin{aligned} &\mathcal{F}_{a}:\mathbb{R}\rightarrow(0,1)\text{ is
a }C^{1}\text{- diffeomorphism};\\\\[5.0pt]
&0<\mathcal{F}_{a}^{\prime}(\eta)\leq\mathcal{F}_{a}(\eta)\leq\exp(\eta),\quad\eta\in\mathbb{R}.\end{aligned}\right.$
(H2)

Observe that the boundedness of the image of $\mathcal{F}_{a}$ reflects the
boundedness of the anion vacancy density. Such a choice necessarily leads to
nonlinear diffusion for the densities in (2.1a). For anion vacancies in PSCs
the Fermi-Dirac integral of order $\mathit{-1}$ is chosen to reflect the
limitation of ion concentration by the lattice sites available in the crystal
[3]. This particular statistics function reads

$\displaystyle F_{-1}(\eta)=\frac{1}{\exp(-\eta)+1}.$ (2.11)

Note that for both assumptions (H1) and (H2), the positivity of the statistics
functions reflects the positivity of the number densities of charge carriers.

###  2.3 Thermodynamic free energy

The thermodynamic free energy for the discussed model is given by the sum of
different energy contributions. Following [4, 29], on the one hand, the
contribution of electrons and holes can be derived from a quasi-free Fermi
gas. On the other hand, the electric contribution to the total energy is given
by the electrostatic field energy. Lastly, assuming an ideal lattice gas [3]
we can derive a consistent energy contribution of anion vacancies which
extends the electric free energy formulation in [29]. Hence, in total the free
energy functional for the PSC model reads

$\displaystyle\mathbb{E}_{f}(t)=$
$\displaystyle\frac{1}{2}\int_{\mathbf{\Omega}}\varepsilon_{s}|\nabla\psi|^{2}\,d\mathbf{x}~{}+\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}\left[k_{B}TN_{\alpha}\Phi_{\alpha}\left(\frac{n_{\alpha}}{N_{\alpha}}\right)-z_{\alpha}E_{\alpha}n_{\alpha}\right]\,d\mathbf{x}$
$\displaystyle+\int_{\mathbf{\Omega}_{\text{intr}}}\left[k_{B}TN_{a}\Phi_{a}\left(\frac{n_{a}}{N_{a}}\right)-z_{a}E_{a}n_{a}\right]\,d\mathbf{x},$

where $\Phi_{\alpha}$ is an antiderivative of $\mathcal{F}_{\alpha}^{-1}$ and
where for the sake of simplicity we neglected external interaction effects of
the electric potential. Thus, for non-degenerate semiconductors, i.e. for
$\mathcal{F}_{n}=\mathcal{F}_{p}=\exp$ and $\mathcal{F}_{a}$ chosen as Fermi-
Dirac integral of order $-1$ the free energy simplifies to

$\begin{split}\mathbb{E}_{f}(t)=&\frac{1}{2}\int_{\mathbf{\Omega}}\varepsilon_{s}|\nabla\psi|^{2}\,d\mathbf{x}~{}+\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}\left[k_{B}Tn_{\alpha}\left(\log\left(\frac{n_{\alpha}}{N_{\alpha}}\right)-1\right)-z_{\alpha}E_{\alpha}n_{\alpha}\right]\,d\mathbf{x}\\\
~{}&+\int_{\mathbf{\Omega}_{\text{intr}}}\left[k_{B}T\left(n_{\alpha}\log\left(\frac{n_{\alpha}}{N_{\alpha}}\right)+\left(N_{\alpha}-n_{\alpha}\right)\log\left(1-\frac{n_{\alpha}}{N_{\alpha}}\right)\right)-z_{a}E_{a}n_{a}\right]\,d\mathbf{x}.~{}\end{split}$
(2.12)

###  2.4 Non-dimensionalization of the model

In this subsection, we derive the relevant non-dimensional parameters of the
model, following [15] and [32, Section 2.4]. Starting from the charge
transport model in (2.1)-(2.4), we rewrite the equations in terms of the
scaled variables given as the ratio of the unscaled physical quantity and the
scaling factors defined in Table 1 (where $U_{T}=k_{B}T/q$ shall denote the
thermal voltage).

In the following we make several simplifications in order to simplify the
presentation and the forthcoming computations. More precisely, we assume from
now on and until the end of Section 5 that the mobilities $\mu_{n}$ and
$\mu_{p}$, the dielelectric permittivity $\varepsilon_{s}$, and the effective
conduction and valence density of states $N_{n}$ and $N_{p}$ are constant in
the domain $\mathbf{\Omega}$. Moreover, we assume that $\mu_{n}=\mu_{p}$ and
$N_{n}=N_{p}$. In practice, the previous quantities vary in each subdomain.
Finally, the band-edge energy $E_{\alpha}$ is assumed to be null for all
moving charge carriers $\alpha$. In Section 6.2, we perform numerical
simulations with heterogeneous parameters and non-zero band-edge energies. All
the analysis of Section 3 and Section 4 can be adapted without the previous
simplifications. However, apart from creating notational overhead, the key
ideas remain the same.

Symbol | Meaning | Scaling factor | Order of magnitude
---|---|---|---
$\>\>\mathbf{x}$ | space variable | $l$ |  $10^{-5}\,$cm
$\psi$, $\psi^{D}$, $\varphi_{\alpha}$, $\varphi^{D}$ | electric and quasi Fermi potentials | $U_{T}$ | $10^{-2}\,\text{V}$
$n_{n}$, $n_{p}$ | densities of electrons and holes | $\tilde{N}$ |  $10^{18}\,$cm-3
$n_{a}$ | density of anion vacancies | $\tilde{N}_{a}$ |  $10^{21}\,$cm-3
$C$ | doping profile | $\tilde{N}$ |  $10^{18}\,$cm-3
$\mu_{n},\mu_{p}$ | electron and hole mobility | $\tilde{\mu}$ |  $10^{0}\,$cm2V-1s-1
$\mu_{a}$ | anion vacancy mobility | $\tilde{\mu}_{a}$ |  $10^{-12}\,$cm2V-1s-1
$t$ | time variable | $\displaystyle\frac{l^{2}}{\tilde{\mu}_{a}U_{T}}$ |  $10^{4}\,$s
$\mathbf{j}_{n},\mathbf{j}_{p}$ | current density for electrons and holes | $\displaystyle\frac{qU_{T}\tilde{N}\tilde{\mu}}{l}$ |  $10^{2}\,$Acm-2
$\mathbf{j}_{a}$ | current density for anion vacancies | $\displaystyle\frac{qU_{T}\tilde{N}_{a}\tilde{\mu}_{a}}{l}$ |  $10^{-7}\,$Acm-2
$R$ | recombination rate | $\displaystyle\frac{\tilde{\mu}U_{T}\tilde{N}}{l^{2}}$ |  $10^{26}\,$cm-3s-1
$G$ | photogeneration rate | $F_{\text{ph}}\alpha_{g}$ |  $10^{22}\,$cm-3s-1
Table 1: Scaling factors for a PSC device based on the default parameters of
[13] at $T=298K$.

In Table 1, the time scale is chosen to be that of the anion vacancies. By
replacing $\tilde{\mu}_{a}$ with $\tilde{\mu}$ in the scaling factor of the
time variable, one could write the dimensionless version adapted to the
electrons and holes time scale. We assume that the scaling factor $\tilde{N}$
is exactly equal to $N_{n}=N_{p}$ and that $\tilde{N}_{a}=N_{a}$. Similarly,
since we assumed, for simplicity, that the mobilities are constant in the
domain, we take $\tilde{\mu}_{a}=\mu_{a}$ and $\tilde{\mu}=\mu_{n}=\mu_{p}$.
By denoting the scaled quantities with the same symbol as the corresponding
unscaled quantities the dimensionless version of the model (2.1)-(2.4) reads

$\displaystyle\nu\,z_{n}\partial_{t}n_{n}+\nabla\cdot\mathbf{j}_{n}$
$\displaystyle=z_{n}\Bigl{(}\gamma\,G(\mathbf{x})-R(n_{n},n_{p})\Bigr{)},$
$\displaystyle{\bf x}\in\mathbf{\Omega},\ t\geq 0,$ (2.13a)
$\displaystyle\nu\,z_{p}\partial_{t}n_{p}+\nabla\cdot\mathbf{j}_{p}$
$\displaystyle=z_{p}\Bigl{(}\gamma\,G(\mathbf{x})-R(n_{n},n_{p})\Bigr{)},$
$\displaystyle{\bf x}\in\mathbf{\Omega},\ t\geq 0,$ (2.13b) $\displaystyle
z_{a}\partial_{t}n_{a}+\nabla\cdot\mathbf{j}_{a}$ $\displaystyle=0,$
$\displaystyle{\bf x}\in\mathbf{\Omega}_{\text{intr}},\ t\geq 0,$ (2.13c)

coupled to the Poisson equation

$-\lambda^{2}\Delta\psi=\left\\{{\begin{aligned}
&\delta(z_{n}n_{n}+z_{p}n_{p}+C({\bf x})),&&{\bf
x}\in{\mathbf{\Omega}}_{\text{HTL}}\cup{\mathbf{\Omega}}_{\text{ETL}},\ t\geq
0,\\\ &z_{a}n_{a}+\delta(z_{n}n_{n}+z_{p}n_{p}+C({\bf x})),&&{\bf
x}\in{\mathbf{\Omega}}_{\text{intr}},\ t\geq 0.\end{aligned}}\right.$ (2.14)

The state equation can be rewritten as

$n_{\alpha}=\mathcal{F}_{\alpha}\Bigl{(}z_{\alpha}(\varphi_{\alpha}-\psi)\Bigr{)},\quad\alpha\in\\{n,p,a\\},$
(2.15)

and we have the following expressions for the charge carrier currents

$\displaystyle\mathbf{j}_{\alpha}=-z_{\alpha}^{2}n_{\alpha}\nabla\varphi_{\alpha},\quad\alpha\in\\{n,p,a\\}.$
(2.16)

There are four dimensionless parameters, the rescaled Debye length, which is
taken with respect to the anion vacancies

$\lambda=\sqrt{\frac{\varepsilon_{s}U_{T}}{l^{2}q\tilde{N}_{a}}},$ (2.17)

the relative mobility of anion vacancies with respect to the mobility of
electrons and holes

$\nu=\frac{\tilde{\mu}_{a}}{\tilde{\mu}},$ (2.18)

the relative concentration of electric carriers with respect to the anion
vacancy concentration

$\delta=\frac{\tilde{N}}{\tilde{N}_{a}},$ (2.19)

and the rescaled photogeneration rate

$\gamma=\frac{F_{\text{ph}}\alpha_{g}l^{2}}{\tilde{\mu}U_{T}\tilde{N}}.$
(2.20)

The parameter $\nu$ can also be interpreted as the ratio between the electric
and ionic carrier time scale. In a typical device all of these parameters are
small. More precisely, $\nu\approx 10^{-12}$, $\lambda\approx 10^{-2}$,
$\delta\approx 10^{-3}$ and $\gamma\approx 10^{-5}$. In particular, the
parameters $\nu$ and $\lambda^{2}$ generate important stiffness in the model,
which motivates the use of a robust implicit-in-time numerical scheme (see
Section 4).

##  3 Continuous entropy-dissipation inequality

For drift-diffusion systems in semiconductor modeling, the natural _a priori_
estimate [20, 21] is based on the evolution of a global quantity which has the
physical meaning of a free energy. In the following, we call this quantity
_(relative) entropy_ , in the sense of _entropy method_ for PDEs rather than
in the physical sense.

###  3.1 Entropy functions

For $\alpha\in\\{n,p,a\\}$, we define the relative entropy function
$\Phi_{\alpha}$, associated with the statistics $\mathcal{F}_{\alpha}$, to be
an anti-derivative of the inverse statistics function namely

$\Phi_{\alpha}^{\prime}(x)=\mathcal{F}_{\alpha}^{-1}(x),\quad x\geq 0.$ (3.1)

Observe that (H1) and (H2) imply that the statistics function is strictly
increasing and therefore $\Phi_{\alpha}$ is strictly convex. Of course
equation (3.1) does not define $\Phi_{\alpha}$ uniquely, but the value of the
constant is not crucial for $\alpha=n,p$ in what follows because we will
introduce relative entropies. The constant may be taken in general to ensure
that $\Phi_{\alpha}$ is non-negative and vanishes at only one point, which is
indeed necessary for $\alpha=a$.

We also define the relative entropy $H_{\alpha}$ by

$\displaystyle
H_{\alpha}(x,y)=\Phi_{\alpha}(x)-\Phi_{\alpha}(y)-\Phi_{\alpha}^{\prime}(y)(x-y),\quad
x\geq 0.$ (3.2)

Observe that $H_{\alpha}$ is non-negative due to the convexity of
$\Phi_{\alpha}$.

#### Examples.

Let us give two examples for the typical statistics functions of the electric
and ionic charge carriers. In the case of the Boltzmann statistics, one has

$\mathcal{F}_{n}(\eta)=\mathcal{F}_{p}(\eta)=e^{\eta},\qquad\Phi_{n}(x)=\Phi_{p}(x)=x\log(x)-x+1.$

In the case of the Fermi-Dirac integral of order $-1$ for $\mathcal{F}_{a}$,
one has

$\mathcal{F}_{a}(\eta)=\frac{1}{\exp(-\eta)+1},\qquad\Phi_{a}(x)=x\log(x)+(1-x)\log(1-x)+\log(2).$

Note that both examples for the mathematical entropy functions coincide with
the respective physical free energy contributions in (2.12).

#### Properties of the entropy and relative entropy function.

Let us state some useful results for the entropy functions. The proofs may be
found in Appendix A.

######  Lemma 3.1.

One has the following bounds on the entropy functions (3.1) and (3.2).

  * (i)

Let $\mathcal{F}_{\alpha}$ be a statistics function satisfying (H1) and
$H_{\alpha}$ be the associated relative entropy function. Then, for any
$\varepsilon>0$ and $y_{0}\geq 0$, there exists a constant
$c_{y_{0},\varepsilon}>0$ such that

$\displaystyle x\leq c_{y_{0},\varepsilon}+\varepsilon
H_{\alpha}(x,y),\quad\text{for all}\;x\geq 0,\ y\in[0,y_{0}].$

  * (ii)

Let $\mathcal{F}_{a}$ be a statistics function satisfying (H2) and $\Phi_{a}$
be the associated entropy function. Then, for any $\varepsilon>0$, there
exists a constant $c_{\varepsilon}>0$ such that

$\displaystyle x\leq c_{\varepsilon}+\varepsilon\Phi_{a}(x),\quad\text{for
all}\;x\geq 0.$ (3.3)

Under a last assumption on the statistics functions for electrons and holes

$\lim_{x\rightarrow+\infty}\frac{H_{\alpha}(x,y_{0})}{\mathcal{F}_{\alpha}^{-1}(x)}=+\infty,\quad\text{for}\;y_{0}\geq
0\quad\text{and}\quad\alpha=n,p,$ (H3)

we have the following result.

######  Lemma 3.2.

Let $\mathcal{F}_{\alpha}$ with $\alpha=n,p$ be a statistics function
satisfying (H1) and (H3). Then, for any $\varepsilon>0$ and $y_{0}\geq 0$,
there exists a constant $c_{y_{0},\varepsilon}>0$ such that

$\displaystyle\max(\mathcal{F}_{\alpha}^{-1}(x),0)\leq
c_{y_{0},\varepsilon}+\varepsilon H_{\alpha}(x,y),\quad\text{for all}\;x\geq
0,\ y\in[0,y_{0}].$

We will also show in Appendix A that the Boltzmann statistics and the Fermi-
Dirac statistics of order $1/2$ both satisfy (H1), (H3), while the Fermi-Dirac
statistics of order $-1$ satisfies (H2).

###  3.2 Proof of the entropy-dissipation inequality

The thermodynamic free energy introduced in Subsection 2.3 is of physical
interest. Now, however, we would like to prove an entropy-dissipation
inequality whose discrete counterpart will allow us to prove the existence of
a discrete solution and the stability of the scheme. For this reason, we
introduce a variation of this functional which from now we will refer to as
total relative entropy in agreement with the mathematical literature. Adapting
the functional of [7, 27] to our system, the total relative entropy with
respect to the Dirichlet boundary values $\psi^{D},\varphi^{D}$ is given by

$\mathbb{E}(t)=\frac{\lambda^{2}}{2}\int_{\mathbf{\Omega}}|\nabla(\psi-\psi^{D})|^{2}\,d\mathbf{x}+\int_{\mathbf{\Omega}_{\text{intr}}}\Phi_{a}(n_{a})\,d\mathbf{x}+\delta\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}H_{\alpha}(n_{\alpha},n_{\alpha}^{D})\,d\mathbf{x},$
(3.4)

where the entropy functions $\Phi_{a}$, $H_{n}$ and $H_{p}$ are given by
(3.1), (3.2) and $n_{\alpha}^{D}$ can be calculated by inserting
$\varphi^{D},\psi^{D}$ into the state equation (2.3). Note that due to our
specific choice for $\Phi_{a}$ the middle term is non-negative as well which
implies that the entropy is non-negative. Taking into account the fact that
$z_{\alpha}^{2}=1$ for $\alpha\in\\{n,p\\}$, the associated non-negative
dissipation $\mathbb{D}$ is defined as

$\displaystyle\mathbb{D}(t)=\frac{\delta}{\nu}\int_{\mathbf{\Omega}}R(n_{n},n_{p})\left(\varphi_{p}-\varphi_{n}\right)\,d\mathbf{x}+\frac{z_{a}^{2}}{2}\int_{\mathbf{\Omega_{\text{intr}}}}n_{a}|\nabla\varphi_{a}|^{2}\,d\mathbf{x}+\frac{\delta}{2\nu}\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}n_{\alpha}|\nabla\varphi_{\alpha}|^{2}\,d\mathbf{x}.$
(3.5)

######  Theorem 3.3.

(Continuous entropy-dissipation inequality) Consider a smooth solution to the
model (2.13)–(2.16), with initial conditions (2.5) and boundary conditions
(2.6), (2.7), (2.8). Then, for any $\varepsilon>0$, there is a constant
$c_{\varepsilon,\mathbf{\Omega}}>0$

$\displaystyle\frac{\text{d}}{\text{d}t}\mathbb{E}(t)+\mathbb{D}(t)\leq
c_{\varepsilon,\mathbf{\Omega}}+\varepsilon\mathbb{E}(t),$ (3.6)

where the entropy is defined in (3.4) and the dissipation of entropy in (3.5).
The constant $c_{\varepsilon,\mathbf{\Omega}}$ depends only on $\varepsilon$,
$\mathbf{\Omega}$, on the boundary data and the photogeneration term via the
norms $\|G\|_{L^{\infty}},\|\varphi^{D}\|_{W^{1,\infty}}$ and
$\|\psi^{D}\|_{W^{1,\infty}}$, on $z_{a}$ and on the dimensionless parameters
$\delta$, $\gamma$ and $\nu$.

######  Remark 3.4 (Thermodynamic equilibrium).

If the boundary data is at thermodynamic equilibrium, _i.e._
$\nabla\varphi^{D}=\nabla\psi^{D}=\mathbf{0}$ and without external generation
of electric carriers, _i.e._ $G=0$, then the entropy-dissipation inequality
simplifies to

$\frac{\text{d}}{\text{d}t}\mathbb{E}(t)+\mathbb{D}(t)\leq 0.$

Indeed, while we do not specify precisely the dependencies of the constant on
the data it is clear from the proof of Theorem 3.3 that the right hand-side of
(3.6) vanishes in this setting. In this case the entropy decays in time and
the solution is expected to converge exponentially fast towards the
thermodynamic equilibrium
$(\varphi_{n}^{\text{eq}},\varphi_{p}^{\text{eq}},\varphi_{a}^{\text{eq}},\psi^{\text{eq}})$.
This thermodynamic equilibrium is such that the quasi Fermi potentials for
electrons and holes are constant on $\mathbf{\Omega}$

$\varphi_{p}^{\text{eq}}=\varphi_{n}^{\text{eq}}=\varphi^{D}$

and $\varphi_{a}^{\text{eq}}$ is constant on
$\mathbf{\Omega}_{\textrm{intr}}$, determined by the conservation of mass for
anion vacancies

$\int_{\mathbf{\Omega}_{\text{intr}}}n_{a}(\varphi_{a}^{\text{eq}},\psi^{\text{eq}})\,d\mathbf{x}=\int_{\mathbf{\Omega}_{\text{intr}}}n_{a}(\mathbf{x},0)\,d\mathbf{x},$

where the electric potential $\psi^{\text{eq}}$ satisfies the following
nonlinear Poisson equation

$-\lambda^{2}\Delta\psi^{\text{eq}}=\left\\{{\begin{aligned}
&\delta\left(n_{p}(\varphi_{p}^{\text{eq}},\psi^{\text{eq}})-n_{n}(\varphi_{n}^{\text{eq}},\psi^{\text{eq}})+C({\bf
x})\right),&&{\bf
x}\in{\mathbf{\Omega}}_{\text{HTL}}\cup{\mathbf{\Omega}}_{\text{ETL}},\\\
~{}&z_{a}n_{a}(\varphi_{a}^{\text{eq}},\psi^{\text{eq}})+\delta\left(n_{p}(\varphi_{p}^{\text{eq}},\psi^{\text{eq}})-n_{n}(\varphi_{n}^{\text{eq}},\psi^{\text{eq}})+C({\bf
x})\right),&&{\bf x}\in{\mathbf{\Omega}}_{\text{intr}}.\
\end{aligned}}\right.$

The system is supplemented with the Dirichlet and Neumann boundary conditions
(2.6a) and (2.6b) for the electric potential. The proof of this asymptotic
behavior is beyond the scope of the present paper but could be investigated
following the lines of the seminal work of Gajewski [20].

###### Proof of Theorem 3.3.

First, let us take the derivative of (3.4) with respect to time

$\displaystyle\frac{\text{d}}{\text{d}t}\mathbb{E}(t)$
$\displaystyle=\lambda^{2}\int_{\mathbf{\Omega}}(\partial_{t}\nabla\psi)\cdot\nabla(\psi-\psi^{D})\,d\mathbf{x}~{}+\int_{\mathbf{\Omega}_{\text{intr}}}\mathcal{F}_{a}^{-1}(n_{a})\partial_{t}n_{a}\,d\mathbf{x}$
(3.7)
$\displaystyle+\delta\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}\left(\mathcal{F}_{\alpha}^{-1}(n_{\alpha})-\mathcal{F}_{\alpha}^{-1}(n_{\alpha}^{D})\right)\partial_{t}n_{\alpha}\,d\mathbf{x}.$

By integrating the first term by parts and using the Poisson equation (2.2)
one obtains

$\lambda^{2}\int_{\mathbf{\Omega}}(\partial_{t}\nabla\psi)\cdot\nabla(\psi-\psi^{D})\,d\mathbf{x}\
=\
\delta\int_{\mathbf{\Omega}}(z_{n}\partial_{t}n_{n}+z_{p}\partial_{t}n_{p})(\psi-\psi^{D})\,d\mathbf{x}+\int_{\mathbf{\Omega}_{\text{intr}}}z_{a}\partial_{t}n_{a}(\psi-\psi^{D})\,d\mathbf{x},$

where all the boundary terms cancel thanks to the boundary conditions (2.6)
and (2.8a). Plugging this back into (3.7) and using the state equation (2.15),
we have

$\displaystyle\frac{\text{d}}{\text{d}t}\mathbb{E}(t)=\int_{\mathbf{\Omega}_{\text{intr}}}z_{a}\left(\varphi_{a}-\psi^{D}\right)\partial_{t}n_{a}\,d\mathbf{x}\\!\\!~{}+\delta\\!\\!\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}z_{\alpha}\left(\varphi_{\alpha}-\varphi^{D}\right)\partial_{t}n_{\alpha}\,d\mathbf{x}.$

Next, we insert the balance equations (2.13) and the definition of the current
densities (2.16)

$\displaystyle\frac{\text{d}}{\text{d}t}\mathbb{E}(t)=$
$\displaystyle-\int_{\mathbf{\Omega}_{\text{intr}}}\nabla\cdot\mathbf{j}_{a}\left(\varphi_{a}-\psi^{D}\right)\,d\mathbf{x}-\frac{\delta}{\nu}\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}\nabla\cdot\mathbf{j}_{\alpha}\left(\varphi_{\alpha}-\varphi^{D}\right)\,d\mathbf{x}$
$\displaystyle+\frac{\delta}{\nu}\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}z_{\alpha}\left(\gamma
G-R\right)\left(\varphi_{\alpha}-\varphi^{D}\right)\,d\mathbf{x}$
$\displaystyle~{}=$
$\displaystyle-\int_{\mathbf{\Omega}_{\text{intr}}}z_{a}^{2}n_{a}\nabla\varphi_{a}\cdot\nabla\left(\varphi_{a}-\psi^{D}\right)\,d\mathbf{x}-\frac{\delta}{\nu}\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}n_{\alpha}\nabla\varphi_{\alpha}\cdot\nabla\left(\varphi_{\alpha}-\varphi^{D}\right)\,d\mathbf{x}$
$\displaystyle+\frac{\delta}{\nu}\int_{\mathbf{\Omega}}\left(\gamma
G-R\right)(\varphi_{p}-\varphi_{n})\,d\mathbf{x},$

where we used $z_{n}=-1=-z_{p}$ and integrated by parts with boundary terms
vanishing again thanks to (2.6b), (2.7), (2.8b) and (2.8c). By expanding the
first terms and using Young’s inequality we get

$\frac{\text{d}}{\text{d}t}\mathbb{E}(t)+\mathbb{D}(t)\leq\frac{\delta\gamma}{\nu}\int_{\mathbf{\Omega}}G\left(\varphi_{p}-\varphi_{n}\right)\,d\mathbf{x}+\frac{z_{a}^{2}}{2}\int_{\mathbf{\Omega}_{\text{intr}}}n_{a}|\nabla\psi^{D}|^{2}\,d\mathbf{x}+\frac{\delta}{2\nu}\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}n_{\alpha}|\nabla\varphi^{D}|^{2}\,d\mathbf{x}.$
(3.8)

It remains to bound the terms of the right-hand side. For the first term on
the right hand side of (3.8) we use the state equation (2.15) and Lemma 3.2 to
find for some $\varepsilon>0$

$\displaystyle\frac{\delta\gamma}{\nu}\int_{\mathbf{\Omega}}G\left(\varphi_{p}-\varphi_{n}\right)\,d\mathbf{x}$
$\displaystyle=\frac{\delta\gamma}{\nu}\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}G\mathcal{F}^{-1}_{\alpha}(n_{\alpha})\,d\mathbf{x}~{}\leq\frac{\delta\gamma}{\nu}||G||_{\infty}\\!\\!\\!\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}\max(\mathcal{F}^{-1}_{\alpha}(n_{\alpha}),0)\,d\mathbf{x}$
$\displaystyle\leq\frac{\delta\gamma}{\nu}||G||_{\infty}\sum_{\alpha\in\\{n,p\\}}\left(c_{y_{\alpha}^{D},\varepsilon}|\mathbf{\Omega}|+\varepsilon\int_{\mathbf{\Omega}}H_{\alpha}\left(n_{\alpha},n_{\alpha}^{D}\right)\,d\mathbf{x}\right),$

where $H_{\alpha}$ is defined in (3.2) and $c_{y_{\alpha}^{D},\varepsilon}$ is
the corresponding constant introduced in Lemma 3.2, where for any species
$\alpha\in\\{n,p\\}$ we introduce

$y_{\alpha}^{D}={\mathcal{F}}_{\alpha}(\|\varphi^{D}\|_{\infty}+\|\psi^{D}\|_{\infty}).$

With help of Lemma 3.1 the second remainder term of (3.8) is estimated by

$\displaystyle\frac{z_{a}^{2}}{2}\int_{\mathbf{\Omega}_{\text{intr}}}n_{a}|\nabla\psi^{D}|^{2}\,d\mathbf{x}+\frac{\delta}{2\nu}\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}n_{\alpha}|\nabla\varphi^{D}|^{2}\,d\mathbf{x}$
$\displaystyle~{}\leq$
$\displaystyle~{}\frac{z_{a}^{2}}{2}||\nabla\psi^{D}||_{\infty}^{2}\int_{\mathbf{\Omega}_{\text{intr}}}n_{a}\,d\mathbf{x}+\frac{\delta}{2\nu}||\nabla\varphi^{D}||_{\infty}^{2}\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}n_{\alpha}\,d\mathbf{x}$
$\displaystyle~{}\leq$
$\displaystyle~{}\max\left\\{\frac{z_{a}^{2}}{2}||\nabla\psi^{D}||_{\infty}^{2},\frac{\delta}{2\nu}||\nabla\varphi^{D}||_{\infty}^{2}\right\\}\left((c_{y_{n}^{D},\varepsilon}+c_{y_{p}^{D},\varepsilon}+c_{\varepsilon})|\mathbf{\Omega}|+3\varepsilon\mathbb{E}\right)$

since the first term in (3.4) is non-negative. Plugging these estimates back
into (3.8) proves the entropy-dissipation estimate (up to a redefinition of
$\varepsilon$). ∎

Using Grönwall’s lemma, an immediate consequence of Theorem 3.3 is that, as
functions of time, the entropy $t\mapsto\mathbb{E}(t)$ and the dissipation
$t\mapsto\mathbb{D}(t)$ are respectively locally bounded and locally
integrable. More precisely, one has the following result.

######  Corollary 3.5.

For any $\varepsilon>0$, one has

$\mathbb{E}(t)+\int_{0}^{t}\mathbb{D}(s)\mathrm{d}s\leq e^{\varepsilon
t}\mathbb{E}(0)+\frac{c_{\varepsilon,\mathbf{\Omega}}}{\varepsilon}(e^{\varepsilon
t}-1)\,,\quad t\geq 0.$

##  4 Discrete version of charge transport model

In this section, we introduce our numerical scheme for (2.13)-(2.16). It is a
finite volume scheme with a two-point flux approximation of the fluxes and a
backward Euler scheme in time. As in the continuous setting, we will show that
an entropy-dissipation relation also holds at the discrete level, ensuring
stability and preservation of the physical structure of the model.

###  4.1 Definition of discretization mesh

First, we introduce the time discretization and the spatial mesh of the domain
$\mathbf{\Omega}$. The mesh, given by the triplet
$\left(\mathcal{T},\mathcal{E},\\{\mathbf{x}_{K}\\}_{K\in\mathcal{T}}\right)$,
will be assumed to be admissible in the sense of [17]. Let $\mathcal{T}$
denote a family of non-empty, convex, open and polygonal control volumes
$K\in\mathcal{T}$, whose Lebesgue measure is denoted by $m_{K}$. For
$K,L\in\mathcal{T}$ with $K\neq L$ we assume that the intersection is empty.
Also, we infer that the union of the closure of all control volumes is equal
to the closure of the domain, i.e.

$\mathbf{\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu}=\bigcup_{K\in\mathcal{T}}\mkern
1.5mu\overline{\mkern-1.5muK\mkern-1.5mu}\mkern 1.5mu.$

The subset of cells contained in the intrisic domain is denoted by
$\mathcal{T}_{\text{intr}}\subset\mathcal{T}$. It is assumed that the closure
of the control volumes in the intrinsic domain form a partition of
$\mathbf{\Omega}_{\text{intr}}$, namely

$\mathbf{\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern
1.5mu_{\text{intr}}}=\bigcup_{K\in\mathcal{T}_{\text{intr}}}\mkern
1.5mu\overline{\mkern-1.5muK\mkern-1.5mu}\mkern 1.5mu.$

Further, we call $\mathcal{E}$ a family of faces, where $\sigma\in\mathcal{E}$
is a closed subset of $\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{\Omega}\mkern-1.5mu}\mkern 1.5mu$ contained
in a hyperplane of $\mathbb{R}^{d}$. Each $\sigma$ has a strictly positive
$(d-1)$-dimensional measure, denoted by $m_{\sigma}$. We use the abbreviation
$\sigma=K|L=\partial K\cap\partial L$ for the intersection between two
distinct control volumes which is either empty or reduces to a face contained
in $\mathcal{E}$. Also, for any $K\in\mathcal{T}$ we assume that there exists
a subset $\mathcal{E}_{K}$ of $\mathcal{E}$ such that the boundary of a
control volume can be described by $\partial
K=\bigcup_{\sigma\in\mathcal{E}_{K}}\sigma$ and, consequently, it follows that
$\mathcal{E}=\bigcup_{K\in\mathcal{T}}\mathcal{E}_{K}$. The set of faces
contained in the intrinsic domain are denoted by

$\mathcal{E}_{\text{intr}}=\\{\sigma\in\mathcal{E}\ \text{s.t.}\
\sigma\subset\mathbf{\overline{\Omega}}_{\text{intr}}\\}.$

Now, we distinguish the faces that are on the boundary of $\mathbf{\Omega}$ by
the notations

$\mathcal{E}^{D}=\\{\sigma\in\mathcal{E}\ \text{s.t.}\
\sigma\subset\mathbf{\Gamma}^{D}\\},\quad\mathcal{E}^{N}=\\{\sigma\in\mathcal{E}\
\text{s.t.}\ \sigma\subset\mathbf{\Gamma}^{N}\\}.$

These sets form partitions of $\mathbf{\Gamma}^{D}$ and $\mathbf{\Gamma}^{N}$,
respectively. We also introduce the set of interior faces in the whole and the
intrinsic domain, respectively

$\mathcal{E}^{\text{int}}=\\{\sigma\in\mathcal{E}\ \text{s.t.}\
\sigma\not\subset\partial\mathbf{\Omega}\\},\quad\mathcal{E}_{\text{intr}}^{\text{int}}=\\{\sigma\in\mathcal{E}_{\text{intr}}\
\text{s.t.}\ \sigma\not\subset\partial\mathbf{\Omega}_{\text{intr}}\\}.$

To each control volume $K\in\mathcal{T}$ we assign a cell center
$\mathbf{x}_{K}\in K$ and we assume that the family of cell centers
$(x_{K})_{K\in{\mathcal{T}}}$ satisfies the orthogonality condition: If $K$
and $L$ share a face $\sigma=K|L$, then the vector

$\overline{\mathbf{x}_{K}\mathbf{x}_{L}}\text{ is orthogonal to }\sigma=K|L.$

For each edge $\sigma\in{\mathcal{E}}$, we define $d_{\sigma}$ as the
Euclidean distance between $\mathbf{x}_{K}$ and $\mathbf{x}_{L}$, if
$\sigma=K|L$ or between $\mathbf{x}_{K}$ and the affine hyperplane spanned by
$\sigma$, if $\sigma\subset\partial\mathbf{\Omega}$. Lastly, we introduce the
transmissibility of the edge $\sigma$:

$\tau_{\sigma}=\frac{m_{\sigma}}{d_{\sigma}}.$

The notations are illustrated in Figure 4.1.

(a)

(b)

Figure 4.1: Neighboring control volumes (a) in the interior of the domain and
(b) near outer boundaries $\mathbf{\Gamma}^{D}$ and $\mathbf{\Gamma}^{N}$.

We assume that the mesh is regular in the following sense. There is a constant
$\xi>0$, which does not depend on the size of the mesh
$h_{\mathcal{T}}=\max_{K}(\mathrm{diam}(K))$ such that

$\forall K\in\mathcal{T},\ \forall\sigma\in\mathcal{E}_{K},\
\left\\{\begin{array}[]{l}d_{\sigma}\geq\xi\,\mathrm{diam}(K)\,,\\\
m_{K}\geq\xi\,\sum_{\sigma\in\mathcal{E}_{K}}m_{\sigma}d_{\sigma}\,.\end{array}\right.$

The regularity assumptions have to be understood as an asymptotic property as
$h_{\mathcal{T}}\to 0$ which are always satisifed on a given mesh due to the
finite number of cells. We remark that Voronoi meshes satisfy all the
assumptions stated in this section.

For the time discretization we decompose the interval $[0,t_{F}]$, for a given
end time $t_{F}>0$ into a finite and increasing number of time steps
$0=t^{1}<\ldots<t^{M}=t_{F}$ with a step-size $\tau^{m}=t^{m}-t^{m-1}$ at time
step $m=2,\dots,M$. We finally introduce $\Delta
t=\max_{m=2,\dots,M}\tau^{m}$.

###  4.2 Finite volume discretization

Now, we introduce the finite volume discretization for (2.13)–(2.16). In what
follows, the quantity $u_{K}^{m}$ represents an approximation of the mean
value of $u(\mathbf{x},t)$ on the cell $K$ at time $t^{m}$, where $u$ is one
of the potentials $\varphi_{n},\varphi_{p}$, $\psi$. In this case, we define
${\mathbf{u}}^{m}=(u_{K}^{m})_{K\in{\mathcal{T}}}$. For $\varphi_{a}$ the
approximation is only given for $K\in\mathcal{T}_{\text{intr}}$, so that we
define
${{\boldsymbol{\varphi}}}_{a}^{m}=(\varphi_{a,K}^{m})_{K\in{\mathcal{T}}_{\textrm{intr}}}$.
The discretizations of the doping profile $C$, the photogeneration rate $G$
and the boundary data $\varphi^{D}$, $\psi^{D}$ are given by

$\chi_{K}=\frac{1}{m_{K}}\int_{K}\chi(\mathbf{x})d\mathbf{x},\quad
K\in\mathcal{T},\ \chi=C,\,G,\,\psi^{D}\text{ or }\varphi^{D},$

and

$\chi_{\sigma}=\frac{1}{m_{\sigma}}\int_{\sigma}\chi(\gamma)d\gamma,\quad\sigma\in\mathcal{E}^{D},\
\chi=\psi^{D},\,\varphi^{D}.$

We discretize in the same way the initial conditions $\varphi_{n}^{0}$,
$\varphi_{p}^{0}$, $\varphi_{a}^{0}$, which lead to the corresponding vectors
${{\boldsymbol{\varphi}}}_{n}^{0}$, ${{\boldsymbol{\varphi}}}_{p}^{0}$ and
${{\boldsymbol{\varphi}}}_{a}^{0}$. The finite volume scheme is formulated as
follows. First, the discrete mass balance equations for the three charge
carriers are given by

$\displaystyle\nu
z_{n}m_{K}\frac{n_{n,K}^{m}-n_{n,K}^{m-1}}{\tau^{m}}+\sum_{\sigma\in\mathcal{E}_{K}}J_{n,K,\sigma}^{m}$
$\displaystyle=z_{n}m_{K}\Bigl{(}\gamma
G_{K}-R(n_{n,K}^{m},n_{p,K}^{m})\Bigr{)},$ $\displaystyle K\in\mathcal{T},\
m\in\mathbb{N},$ (4.1a) $\displaystyle\nu
z_{p}m_{K}\frac{n_{p,K}^{m}-n_{p,K}^{m-1}}{\tau^{m}}+\sum_{\sigma\in\mathcal{E}_{K}}J_{p,K,\sigma}^{m}$
$\displaystyle=z_{p}m_{K}\Bigl{(}\gamma
G_{K}-R(n_{n,K}^{m},n_{p,K}^{m})\Bigr{)},$ $\displaystyle K\in\mathcal{T},\
m\in\mathbb{N},$ (4.1b) $\displaystyle
z_{a}m_{K}\frac{n_{a,K}^{m}-n_{a,K}^{m-1}}{\tau^{m}}+\sum_{\sigma\in\mathcal{E}_{K}}J_{a,K,\sigma}^{m}$
$\displaystyle=0,$ $\displaystyle K\in\mathcal{T}_{\text{intr}},\
m\in\mathbb{N}.$ (4.1c)

They are coupled via the discrete Poisson equation

$-\lambda^{2}\sum_{\sigma\in\mathcal{E}_{K}}\tau_{\sigma}D_{K,\sigma}\boldsymbol{\psi}^{m}=\left\\{{\begin{aligned}
&\delta
m_{K}(z_{n}n_{n,K}^{m}+z_{p}n_{p,K}^{m}+C_{K}),&&K\in\mathcal{T}\setminus\mathcal{T}_{\text{intr}},\
m\in\mathbb{N},\\\ &m_{K}z_{a}n_{a,K}^{m}+\delta
m_{K}(z_{n}n_{n,K}^{m}+z_{p}n_{p,K}^{m}+C_{K}),&&K\in\mathcal{T}_{\text{intr}},\
m\in\mathbb{N}.\end{aligned}}\right.$ (4.2)

In the previous equation, the notation $D_{K,\sigma}$ denotes the finite
difference operator acting on vectors of unknowns
${\boldsymbol{u}}=(u_{K})_{K}$ and is given by

$D_{K,\sigma}{\boldsymbol{u}}=\left\\{\begin{array}[]{ll}u_{L}-u_{K},&\text{if
}\sigma=K|L,\\\ u_{\sigma}^{D}-u_{K},&\text{if }\sigma\in\mathcal{E}^{D},\\\
0,&\text{otherwise.}\end{array}\right.$ (4.3)

The discrete densities are given by the state equation (2.15) inside the
domain and at the Dirichlet boundary, namely

$\displaystyle{\boldsymbol{n}}_{\alpha}^{m}$
$\displaystyle=\mathcal{F}_{\alpha}\left(z_{\alpha}({\boldsymbol{\varphi}}_{\alpha}^{m}-{\boldsymbol{\psi}}^{m})\right),$
$\displaystyle\alpha\in\\{n,p,a\\},\ m\in\mathbb{N},$ (4.4a)
$\displaystyle{\boldsymbol{n}}_{\alpha}^{D}$
$\displaystyle=\mathcal{F}_{\alpha}\left(z_{\alpha}({\boldsymbol{\varphi}}_{\alpha}^{D}-{\boldsymbol{\psi}}^{D})\right),$
$\displaystyle\alpha\in\\{n,p\\},$ (4.4b) $\displaystyle
n_{\alpha,\sigma}^{D}$
$\displaystyle=\mathcal{F}_{\alpha}\left(z_{\alpha}(\varphi_{\sigma}^{D}-\psi^{D}_{\sigma})\right),$
$\displaystyle\alpha\in\\{n,p\\},\ \sigma\in\mathcal{E}^{D},$ (4.4c)

where the statistics function is applied pointwise to the input vector. Let us
remark that the discrete values of the boundary densities, defined by (4.4b)
and (4.4c), are bounded but that the upper bound may differ from
$\|n_{\alpha}^{D}\|_{\infty}$. Indeed, for $\alpha\in\\{n,p\\}$, we have

$\max\Bigl{(}\max_{K\in{\mathcal{T}}}n_{\alpha,K}^{D},\max_{\sigma\in{\mathcal{E}}^{D}}n_{\alpha,\sigma}^{D}\Bigr{)}\leq{\mathcal{F}}_{\alpha}(\|\varphi^{D}\|_{\infty}+\|\psi^{D}\|_{\infty})=y_{\alpha}^{D}.$
(4.5)

We use the excess chemical potential scheme (frequently called Sedan scheme)
as TPFA scheme for $J_{\alpha,K,\sigma}^{m}$. The earliest reference, we could
find for this thermodynamically consistent flux discretization scheme is [40].
The definition of the numerical flux is based on the following reformulation
of the currents:

${\mathbf{j}}_{\alpha}=-z_{\alpha}\left(\nabla
n_{\alpha}+n_{\alpha}\nabla(z_{\alpha}\varphi_{\alpha}-\log
n_{\alpha})\right)$

and on the approximation of convection-diffusion fluxes by Scharfetter-Gummel
numerical fluxes, see [11, 36]. For the electrons and holes, it reads

$\forall\alpha\in\\{n,p\\},\quad
J_{\alpha,K,\sigma}^{m}=\left\\{\begin{array}[]{ll}-z_{\alpha}\tau_{\sigma}\Bigl{(}B\left(-Q_{\alpha,K,\sigma}^{m}\right)n_{\alpha,L}^{m}-B\left(Q_{\alpha,K,\sigma}^{m}\right)n_{\alpha,K}^{m}\Bigr{)},&\text{if
}\sigma=K|L,\\\
-z_{\alpha}\tau_{\sigma}\Bigl{(}B\left(-Q_{\alpha,K,\sigma}^{m}\right)n_{\alpha,\sigma}^{D}-B\left(Q_{\alpha,K,\sigma}^{m}\right)n_{\alpha,K}^{m}\Bigr{)},&\text{if
}\sigma\in\mathcal{E}^{D},\\\ 0,&\text{otherwise}.\end{array}\right.$ (4.6)

For the anion vacancies it is given by

$J_{a,K,\sigma}^{m}=\left\\{\begin{array}[]{ll}-z_{a}\tau_{\sigma}\Bigl{(}B\left(-Q_{a,K,\sigma}^{m}\right)n_{a,L}^{m}-B\left(Q_{a,K,\sigma}^{m}\right)n_{a,K}^{m}\Bigr{)},&\text{if
}\sigma=K|L\in\mathcal{E}^{\text{int}}_{\text{intr}},\\\
0,&\text{otherwise}.\end{array}\right.$ (4.7)

The quantity $Q_{\alpha,K,\sigma}^{m}$ is defined as

$Q_{\alpha,K,\sigma}^{m}=D_{K,\sigma}\left(z_{\alpha}{\boldsymbol{\varphi}}_{\alpha}^{m}-\log{\boldsymbol{n}}_{\alpha}^{m}\right)$
(4.8)

with $K\in\mathcal{T}$,
$\sigma\in\mathcal{E}_{K}\cap(\mathcal{E}^{\text{int}}\cup\mathcal{E}^{D})$
for electrons and holes ($\alpha\in\\{n,p\\}$) and
$K\in\mathcal{T}_{\text{intr}}$,
$\sigma\in\mathcal{E}_{K}\cap\mathcal{E}^{\text{int}}_{\text{intr}}$ in the
case of anion vacancies ($\alpha=a$). In the previous formula, the logarithm
is applied componentwise. Lastly, the function $B$ denotes the Bernoulli
function

$B(x)=\frac{x}{\exp(x)-1},\;\text{for}\;x\neq 0\quad\text{and}\quad B(0)=1.$
(4.9)

Note that the fluxes are locally conservative in the sense that for
$\sigma=K|L$

$0=J_{\alpha,K,\sigma}^{m}+J_{\alpha,L,\sigma}^{m}=Q_{\alpha,K,\sigma}^{m}+Q_{\alpha,L,\sigma}^{m}=D_{K,\sigma}{\boldsymbol{\psi}}^{m}+D_{L,\sigma}{\boldsymbol{\psi}}^{m}=D_{K,\sigma}{\boldsymbol{\varphi}}_{\alpha}^{m}+D_{L,\sigma}{\boldsymbol{\varphi}}_{\alpha}^{m}.$
(4.10)

Since the fluxes $J_{\alpha,K,\sigma}^{m}$ and $J_{\alpha,L,\sigma}^{m}$ agree
up to sign for any interior edge, we introduce the notation

$D_{\sigma}{\boldsymbol{u}}=|D_{K,\sigma}{\boldsymbol{u}}|\textrm{ for }\
\sigma\in\mathcal{E}_{K},\textrm{ for
}{\boldsymbol{u}}={\boldsymbol{\psi}}^{m},{\boldsymbol{\varphi}}_{\alpha}^{m},\
\alpha\in\\{n,p,a\\}.$

######  Remark 4.1 (Boundary conditions).

Observe that all the boundary conditions have been considered in the
definition of the scheme. The external boundary conditions (2.6) for the
electric potential are handled in the definition of (4.3). For the quasi Fermi
potentials of electrons and holes external boundary conditions are included in
the definition of (4.4c) and (4.6) as well as (4.3) and (4.8). The Neumann
boundary conditions for anion vacancies (2.7) are included in the definition
of (4.7). Finally, observe that the continuity of fluxes of electrons, holes
and electric potential through the interfaces $\mathbf{\Sigma}_{\text{ETL}}$
and $\mathbf{\Sigma}_{\text{HTL}}$ is automatically ensured thanks to (4.10).

######  Remark 4.2 (Mobilities, permittivity and band-edge energies).

As explained at the beginning of Section 2.4, we made several simplifications
concerning the mobilities, permittivity and band-edge energy in order to
lighten the presentation. The generalization of the present scheme to take
into account non-constant mobilities $\mu_{n},\mu_{p}$ and permittivity
$\varepsilon_{s}$ amounts to introducing a consistent prefactor depending on
the edge $\sigma$ in formula (4.6) and in the sum of the left hand side of
(4.2) respectively. To take into account non-zero band-edge energies
$E_{\alpha}$ one needs to add the corresponding term to the quasi Fermi
potential $\boldsymbol{\varphi}_{\alpha}$ in (4.8).

At first glance, it might not be obvious why the fluxes (4.6) and (4.7) are
discrete versions of (2.16). It turns out that one can define

$\overline{n}^{m}_{\alpha,\sigma}:=\left\\{\begin{array}[]{ll}\displaystyle\frac{B(-Q_{\alpha,K,\sigma}^{m})n_{\alpha,L}^{m}-B(Q_{\alpha,K,\sigma}^{m})n_{\alpha,K}^{m}}{z_{\alpha}(\varphi_{\alpha,L}^{m}-\varphi_{\alpha,K}^{m})},&\sigma=K|L,\
\alpha\in\\{n,p\\},\\\\[10.00002pt]
\displaystyle\frac{B(-Q_{\alpha,K,\sigma}^{m})n_{\alpha,\sigma}^{D}-B(Q_{\alpha,K,\sigma}^{m})n_{\alpha,K}^{m}}{z_{\alpha}(\varphi_{\sigma}^{D}-\varphi_{\alpha,K}^{m})},&\sigma\in\mathcal{E}^{D}\cap\mathcal{E}_{K},\,\
\alpha\in\\{n,p\\},\end{array}\right.$ (4.11)

and

$\overline{n}^{m}_{a,\sigma}:=\displaystyle\frac{B(-Q_{a,K,\sigma}^{m})n_{a,L}^{m}-B(Q_{a,K,\sigma}^{m})n_{a,K}^{m}}{z_{\alpha}(\varphi_{a,L}^{m}-\varphi_{a,K}^{m})},\
\ \sigma=K|L\in{\mathcal{E}}_{\textrm{intr}}^{\textrm{int}},$ (4.12)

so that the fluxes can be rewritten to

${J}_{\alpha,K,\sigma}^{m}=-\tau_{\sigma}z_{\alpha}^{2}\overline{n}^{m}_{\alpha,\sigma}\,D_{K,\sigma}\boldsymbol{\varphi}_{\alpha}^{m},\,\quad\text{for
all}\;\alpha\in\\{n,p,a\\}.$ (4.13)

Observe that $\overline{n}^{m}_{\alpha,\sigma}$ is well-defined in the sense
that thanks to (4.10) it depends only on the edge (and not nodal values) as
well as the fact that a boundary edge has only one associated control volume.
The reformulation of the fluxes (4.13) now is closer to (2.16) but the analogy
would not be complete, if $\overline{n}^{m}_{\alpha,\sigma}$ is not consistent
with the density at the interface $\sigma$. This is actually the case as the
following lemma shows. It is adapted from [11, Lemma 3.1].

######  Lemma 4.3.

The interface value $\overline{n}^{m}_{\alpha,\sigma}$ defined by (4.11) is a
convex combination of $n^{m}_{\alpha,K}$ and $n^{m}_{\alpha,L}$ (_resp._
$n^{D}_{\alpha,\sigma}$), if $\sigma=K|L$ (_resp._
$\sigma\in\mathcal{E}^{D}$). In particular it is framed between the minimum
and maximum of the two values. The same result holds for
$\overline{n}^{m}_{a,\sigma}$ defined by (4.12) for
$\sigma=K|L\in{\mathcal{E}}_{\textrm{intr}}^{\textrm{int}}$.

###### Proof.

It suffices to observe that for $\sigma=K|L$ (the boundary case can be readily
adapted),

$\overline{n}^{m}_{\alpha,\sigma}=\frac{B(y)-B(x)}{x-y}n_{\alpha,L}^{m}+\frac{B(-x)-B(-y)}{x-y}n_{\alpha,K}^{m},$

with $x=D_{K,\sigma}\log{\boldsymbol{n}}_{\alpha}^{m}$ and
$y=-Q_{\alpha,K,\sigma}^{m}$. To see this, use the expression of the Bernoulli
function $B$ to get that the coefficients are non-negative and sum to $1$. We
refer to [11] for additional details concerning this computation. ∎

###  4.3 Discrete entropy-dissipation inequality

In the following, we derive a discrete counterpart of (3.6) for the discrete
relative entropy ($m\in\mathbb{N}$)

$\mathbb{E}_{\mathcal{T}}^{m}=\,\frac{\lambda^{2}}{2}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}\left(D_{\sigma}(\boldsymbol{\psi}^{m}-\boldsymbol{\psi}^{D})\right)^{2}+\sum_{K\in\mathcal{T}_{\text{intr}}}m_{K}\Phi_{a}(n_{a,K}^{m})\\\
+\delta\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}m_{K}H_{\alpha}(n_{\alpha,K}^{m},n_{\alpha,K}^{D}).$
(4.14)

We recall that the entropy functions $\Phi_{a}$ and $H_{\alpha}$ are defined
in (3.1) and (3.2). The corresponding discrete non-negative dissipation
$\mathbb{D}_{\mathcal{T}}^{m}$ for $m\in\mathbb{N}$ is given by

$\mathbb{D}_{\mathcal{T}}^{m}=\frac{z_{a}^{2}}{2}\sum_{\sigma\in\mathcal{E}_{\text{intr}}^{\text{int}}}\tau_{\sigma}\overline{n}^{m}_{a,\sigma}(D_{\sigma}\boldsymbol{\varphi}_{a}^{m})^{2}+\frac{\delta}{2\nu}\sum_{\alpha\in\\{n,p\\}}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}\overline{n}^{m}_{\alpha,\sigma}(D_{\sigma}\boldsymbol{\varphi}_{\alpha}^{m})^{2}\\\
+\frac{\delta}{\nu}\sum_{K\in\mathcal{T}}m_{K}R(n^{m}_{n,K},n^{m}_{p,K})\left(\varphi_{p,K}^{m}-\varphi_{n,K}^{m}\right).$
(4.15)

######  Theorem 4.4.

(Discrete entropy-dissipation inequality) For any solution to the finite
volume scheme (4.1)–(4.8) one has the following entropy-dissipation
inequality: For any $\varepsilon>0$, there is a constant
$c_{\varepsilon,\mathbf{\Omega},\xi}>0$ such that for any $m\in\mathbb{N}$,
one has

$\frac{\mathbb{E}_{\mathcal{T}}^{m}-\mathbb{E}_{\mathcal{T}}^{m-1}}{\tau^{m}}+\mathbb{D}_{\mathcal{T}}^{m}\leq
c_{\varepsilon,\mathbf{\Omega},\xi}+\varepsilon\mathbb{E}_{\mathcal{T}}^{m}.$
(4.16)

The constant $c_{\varepsilon,\mathbf{\Omega},\xi}$ depends solely on
$\varepsilon$, the measure of $\mathbf{\Omega}$, the mesh regularity $\xi>0$,
the boundary data and the photogeneration term via the norms
$\|G\|_{L^{\infty}},\|\varphi^{D}\|_{W^{1,\infty}}$ and
$\|\psi^{D}\|_{W^{1,\infty}}$, as well as on $z_{a}$ and the dimensionless
parameters $\delta$, $\gamma$ and $\nu$. If $G=0$ and
$\nabla\varphi^{D}=\nabla\psi^{D}=\mathbf{0}$, then the right hand-side of
(4.16) vanishes.

###### Proof.

Let us start by considering the difference of the entropies at time $t^{m}$
and $t^{m-1}$, that is

$\displaystyle\mathbb{E}_{\mathcal{T}}^{m}-\mathbb{E}_{\mathcal{T}}^{m-1}=$
$\displaystyle\,\frac{\lambda^{2}}{2}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}\left(\left(D_{\sigma}(\boldsymbol{\psi}^{m}-\boldsymbol{\psi}^{D})\right)^{2}-\left(D_{\sigma}(\boldsymbol{\psi}^{m-1}-\boldsymbol{\psi}^{D})\right)^{2}\right)$
$\displaystyle+\sum_{K\in\mathcal{T}_{\text{intr}}}m_{K}\left(\Phi_{a}(n_{a,K}^{m})-\Phi_{a}(n_{a,K}^{m-1})\right)$
$\displaystyle+\delta\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}m_{K}\left(\Phi_{\alpha}(n_{\alpha,K}^{m})-\Phi_{\alpha}(n_{\alpha,K}^{m-1})-\Phi_{\alpha}^{\prime}\left(n_{\alpha,K}^{D}\right)(n_{\alpha,K}^{m}-n_{\alpha,K}^{m-1})\right).$

Using a convexity inequality in every sum one finds

$\displaystyle\mathbb{E}_{\mathcal{T}}^{m}-\mathbb{E}_{\mathcal{T}}^{m-1}\leq$
$\displaystyle\,\lambda^{2}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}D_{K,\sigma}\left(\boldsymbol{\psi}^{m}-\boldsymbol{\psi}^{D}\right)D_{K,\sigma}\left(\mathbf{\psi}^{m}-\mathbf{\psi}^{m-1}\right)$
$\displaystyle+\sum_{K\in\mathcal{T}_{\text{intr}}}m_{K}\mathcal{F}^{-1}_{a}(n_{a,K}^{m})\Bigl{(}n_{a,K}^{m}-n_{a,K}^{m-1}\Bigr{)}$
$\displaystyle+\delta\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}m_{K}\left(\mathcal{F}^{-1}_{\alpha}(n_{\alpha,K}^{m})-\mathcal{F}^{-1}_{\alpha}\left(n_{\alpha,K}^{D}\right)\right)\left(n_{\alpha,K}^{m}-n_{\alpha,K}^{m-1}\right).$

In order to compute the first sum, we use a discrete integration by parts
(consisting in reordering sums by using the conservativity relations on
fluxes) and the discrete Poisson equation (4.2) to get

$\displaystyle\lambda^{2}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}D_{K,\sigma}\left(\boldsymbol{\psi}^{m}-\boldsymbol{\psi}^{D}\right)D_{K,\sigma}\left(\boldsymbol{\psi}^{m}-\boldsymbol{\psi}^{m-1}\right)$
$\displaystyle=$
$\displaystyle-\lambda^{2}\sum_{K\in\mathcal{T}}\sum_{\sigma\in\mathcal{E}_{K}}\tau_{\sigma}D_{K,\sigma}\left(\boldsymbol{\psi}^{m}-\boldsymbol{\psi}^{m-1}\right)\left(\psi_{K}^{m}-\psi_{K}^{D}\right)$
$\displaystyle=$ $\displaystyle\
\delta\\!\\!\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}m_{K}z_{\alpha}\left(n_{\alpha,K}^{m}-n_{\alpha,K}^{m-1}\right)\left(\psi_{K}^{m}-\psi_{K}^{D}\right)~{}+\sum_{K\in\mathcal{T}_{\text{intr}}}m_{K}z_{a}\left(n_{a,K}^{m}-n_{a,K}^{m-1}\right)\left(\psi_{K}^{m}-\psi_{K}^{D}\right).$

Plugging this relation back into the initial estimate and using the relation
in (4.4), we obtain

$\displaystyle\mathbb{E}_{\mathcal{T}}^{m}-\mathbb{E}_{\mathcal{T}}^{m-1}\leq$
$\displaystyle\
\delta\\!\\!\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}m_{K}z_{\alpha}\left(\varphi_{\alpha,K}^{m}-\varphi_{K}^{D}\right)\left(n_{\alpha,K}^{m}-n_{\alpha,K}^{m-1}\right)$
$\displaystyle+\sum_{K\in\mathcal{T}_{\text{intr}}}m_{K}z_{a}\left(\varphi_{a,K}^{m}-\psi_{K}^{D}\right)\left(n_{a,K}^{m}-n_{a,K}^{m-1}\right).$

Now, divide by the time step size $\tau^{m}$ and insert the mass balances in
(4.1)

$\frac{\mathbb{E}_{\mathcal{T}}^{m}-\mathbb{E}_{\mathcal{T}}^{m-1}}{\tau^{m}}\leq-\frac{\delta}{\nu}\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}\sum_{\sigma\in\mathcal{E}_{K}}J_{\alpha,K,\sigma}^{m}\left(\varphi_{\alpha,K}^{m}-\varphi_{K}^{D}\right)-\sum_{K\in\mathcal{T}_{\text{intr}}}\sum_{\sigma\in\mathcal{E}_{K}}J_{a,K,\sigma}^{m}\left(\varphi_{a,K}^{m}-\psi_{K}^{D}\right)\\\
+\frac{\delta}{\nu}\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}z_{\alpha}m_{K}\Bigl{(}\gamma
G_{K}-R(n_{n,K}^{m},n_{p,K}^{m})\Bigr{)}\left(\varphi_{\alpha,K}^{m}-\varphi_{K}^{D}\right).$

Next, we insert the formulas for the fluxes (4.13) with $z_{\alpha}^{2}=1$,
$\alpha=n,p$, and perform a discrete integration by parts to deduce

$\frac{\mathbb{E}_{\mathcal{T}}^{m}-\mathbb{E}_{\mathcal{T}}^{m-1}}{\tau^{m}}\leq-\frac{\delta}{\nu}\sum_{\alpha\in\\{n,p\\}}\sum_{\sigma\in\mathcal{E}^{\text{int}}\cup\mathcal{E}^{D}}\tau_{\sigma}\overline{n}^{m}_{\alpha,\sigma}\,D_{K,\sigma}\boldsymbol{\varphi}_{\alpha}^{m}D_{K,\sigma}\left(\boldsymbol{\varphi}_{\alpha}^{m}-\boldsymbol{\varphi}^{D}\right)\\\
~{}-z_{a}^{2}\sum_{\sigma\in\mathcal{E}^{\text{int}}_{\text{intr}}}\tau_{\sigma}\overline{n}^{m}_{a,\sigma}\,D_{K,\sigma}\boldsymbol{\varphi}_{a}^{m}D_{K,\sigma}\left(\boldsymbol{\varphi}_{a}^{m}-\boldsymbol{\psi}^{D}\right)\\\
~{}-\frac{\delta}{\nu}\sum_{K\in\mathcal{T}}m_{K}R(n_{n,K}^{m},n_{p,K}^{m})\left(\varphi_{p,K}^{m}-\varphi_{n,K}^{m}\right)+\frac{\delta\gamma}{\nu}\sum_{K\in\mathcal{T}}m_{K}G_{K}\left(\varphi_{p,K}^{m}-\varphi_{n,K}^{m}\right).$

After using the inequality $-a(a-b)\leq-(a^{2}-b^{2})/2$ in the first two
sums, we obtain

$\begin{split}\frac{\mathbb{E}_{\mathcal{T}}^{m}-\mathbb{E}_{\mathcal{T}}^{m-1}}{\tau^{m}}+\mathbb{D}_{\mathcal{T}}^{m}&\leq~{}\frac{\delta}{2\nu}\sum_{\alpha\in\\{n,p\\}}\sum_{\sigma\in\mathcal{E}^{\text{int}}\cup\mathcal{E}^{D}}\tau_{\sigma}\overline{n}^{m}_{\alpha,\sigma}\,(D_{\sigma}\boldsymbol{\varphi}^{D})^{2}+\frac{z_{a}^{2}}{2}\sum_{\sigma\in\mathcal{E}^{\text{int}}_{\text{intr}}}\tau_{\sigma}\overline{n}^{m}_{a,\sigma}\,(D_{\sigma}\boldsymbol{\psi}^{D})^{2}\\\
~{}&\quad\;+\frac{\delta\gamma}{\nu}\sum_{K\in\mathcal{T}}m_{K}G_{K}\left(\varphi_{p,K}^{m}-\varphi_{n,K}^{m}\right).\end{split}$
(4.17)

Observe that at this stage it is obvious that, if $G=0$ and
$\nabla\varphi^{D}=\nabla\psi^{D}=\mathbf{0}$, then the entropy-dissipation
inequality of the theorem holds for a vanishing right hand-side. In the
general case, it remains to estimate the different remainder terms in the
right-hand-side of (4.17).

For the first and second remainder terms in (4.17) we need the following
intermediate result. Let $\sigma=K|L\in\mathcal{E}^{\text{int}}$, then

$\displaystyle\frac{D_{\sigma}\boldsymbol{\varphi}^{D}}{d_{\sigma}}$
$\displaystyle\leq\frac{1}{d_{\sigma}m_{K}m_{L}}\int_{K}\int_{L}|\varphi^{D}(\mathbf{x})-\varphi^{D}(\mathbf{y})|\,d\mathbf{x}\,d\mathbf{y}$
$\displaystyle\leq\frac{\mathrm{diam}(K)+\mathrm{diam}(L)}{d_{\sigma}}\|\nabla\varphi^{D}\|_{L^{\infty}}\leq\frac{2}{\xi}\|\nabla\varphi^{D}\|_{L^{\infty}}\,.$

The last inequality holds also, if $\sigma\in\mathcal{E}^{D}$ or by replacing
$\varphi^{D}$ with $\psi^{D}$. Let us now go back to the first remainder term
of (4.17). We set

$S_{1}=\frac{\delta}{2\nu}\sum_{\alpha\in\\{n,p\\}}\sum_{\sigma\in\mathcal{E}^{\text{int}}\cup\mathcal{E}^{D}}\tau_{\sigma}\overline{n}^{m}_{\alpha,\sigma}\,(D_{\sigma}\boldsymbol{\varphi}^{D})^{2}.$

Since $\overline{n}^{m}_{\alpha,\sigma}$ is the convex combination of two non-
negative unknowns (see Lemma 4.3) it is bounded from above by the sum of these
unknowns. It yields

$\displaystyle S_{1}$
$\displaystyle\leq\frac{2\delta}{\nu\xi^{2}}\|\nabla\varphi^{D}\|_{L^{\infty}}^{2}\sum_{\alpha\in\\{n,p\\}}\sum_{\sigma\in\mathcal{E}^{\text{int}}\cup\mathcal{E}^{D}}m_{\sigma}\,d_{\sigma}\overline{n}^{m}_{\alpha,\sigma}$
$\displaystyle\leq\frac{2\delta}{\nu\xi^{2}}\|\nabla\varphi^{D}\|_{L^{\infty}}^{2}\sum_{\alpha\in\\{n,p\\}}\left(2\sum_{K\in\mathcal{T}}n_{\alpha,K}^{m}\sum_{\sigma\in\mathcal{E}_{K}}m_{\sigma}\,d_{\sigma}+\sum_{\sigma\in\mathcal{E}^{D}}m_{\sigma}\,d_{\sigma}y_{\alpha}^{D}\right)$
$\displaystyle\leq\frac{2\delta}{\nu\xi^{2}}\|\nabla\varphi^{D}\|_{L^{\infty}}^{2}\sum_{\alpha\in\\{n,p\\}}\left(\frac{2}{\xi}\sum_{K\in\mathcal{T}}m_{K}n_{\alpha,K}^{m}+|\mathbf{\Omega}|y_{\alpha}^{D}\right)$
$\displaystyle\leq\frac{2\delta}{\nu\xi^{2}}\|\nabla\varphi^{D}\|_{L^{\infty}}^{2}\left(\frac{2}{\xi}\varepsilon\mathbb{E}_{\mathcal{T}}^{m}+\sum_{\alpha\in\\{n,p\\}}\left(\frac{2}{\xi}|\mathbf{\Omega}|c_{y_{\alpha}^{D},\varepsilon}+|\mathbf{\Omega}|y_{\alpha}^{D}\right)\right)\,,$

where $y_{\alpha}^{D}$ has been defined in (4.5) and
$c_{y_{\alpha}^{D},\varepsilon}>0$ is the constant of the inequality (i) in
Lemma 3.1. Similarly, by using (ii) in Lemma 3.1 we obtain that the second
remainder term satisfies

$\displaystyle
S_{2}=\frac{z_{a}^{2}}{2}\sum_{\sigma\in\mathcal{E}^{\text{int}}_{\text{intr}}}\tau_{\sigma}\overline{n}^{m}_{a,\sigma}\,(D_{\sigma}\boldsymbol{\psi}^{D})^{2}$
$\displaystyle\leq\frac{4z_{a}^{2}}{\xi^{3}}\|\nabla\psi^{D}\|_{L^{\infty}}^{2}(\varepsilon\mathbb{E}_{\mathcal{T}}^{m}+c_{\varepsilon}|\mathbf{\Omega}|)\,.$

For the last remainder term coming from the photogeneration we set

$S_{3}=\frac{\delta\gamma}{\nu}\sum_{K\in\mathcal{T}}m_{K}G_{K}\left(\varphi_{p,K}^{m}-\varphi_{n,K}^{m}\right)=\frac{\delta\gamma}{\nu}\sum_{K\in\mathcal{T}}m_{K}G_{K}\left(\mathcal{F}^{-1}_{n}(n_{n,K}^{m})+\mathcal{F}^{-1}_{p}(n_{p,K}^{m})\right),$

and we use the state equation (4.4a) and Lemma 3.2 to estimate

$\displaystyle S_{3}$
$\displaystyle\leq\frac{\delta\gamma}{\nu}\left(\max_{K\in\mathcal{T}}G_{K}\right)\sum_{K\in\mathcal{T}}m_{K}\Bigl{(}\max(\mathcal{F}^{-1}_{n}(n_{n,K}^{m}),0)+\max(\mathcal{F}^{-1}_{p}(n_{p,K}^{m}),0)\Bigr{)}$
$\displaystyle\leq\frac{\delta\gamma}{\nu}||G||_{L^{\infty}}\\!\\!\\!\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}m_{K}\left(c_{y_{\alpha}^{D},\varepsilon}+\varepsilon
H_{\alpha}(n_{\alpha,K}^{m},n_{\alpha,K}^{D})\right)$
$\displaystyle\leq\frac{\delta\gamma}{\nu}||G||_{L^{\infty}}\left(|\mathbf{\Omega}|c_{y_{\alpha}^{D},\varepsilon}+\varepsilon\mathbb{E}_{\mathcal{T}}^{m}\right)\,.$

For the last step it is important to remember that each term in the definition
of the entropy is non-negative. Therefore, if we combine everything back into
(4.17) we find

$\frac{\mathbb{E}_{\mathcal{T}}^{m}-\mathbb{E}_{\mathcal{T}}^{m-1}}{\tau^{m}}+\mathbb{D}_{\mathcal{T}}^{m}~{}\leq~{}c_{1,\xi}\Bigl{(}\varepsilon\mathbb{E}_{\mathcal{T}}^{m}+c_{2,\varepsilon,\mathbf{\Omega}}\Bigr{)},$

for some constants $c_{1,\xi},c_{2,\varepsilon,\mathbf{\Omega}}>0$ depending
on all the aforementioned quantities. Since $c_{1,\xi}$ does not depend on
$\varepsilon$, this is equivalent to the desired inequality (4.16) up to a
redefinition of $\varepsilon$. ∎

From the discrete entropy-dissipation inequality (4.16), we can deduce some
bounds on the entropy $\mathbb{E}_{\mathcal{T}}^{m}$ and on the cumulated
dissipation $\sum_{k=1}^{m}\tau^{k}\mathbb{D}_{\mathcal{T}}^{k}$ for any $m>0$
thanks to a discrete Grönwall’s Lemma. Corollary 4.5 states the discrete
counterpart of Corollary 3.5.

######  Corollary 4.5.

Provided that $\varepsilon<(\Delta t)^{-1}$, one has for any $m\geq 1$ that

$\mathbb{E}_{\mathcal{T}}^{m}+\sum_{k=1}^{m}\tau^{k}\mathbb{D}_{\mathcal{T}}^{k}\leq(\mathbb{E}_{\mathcal{T}}^{0}+c_{\varepsilon,\mathbf{\Omega},\xi}t^{m})(1-\varepsilon\Delta
t)^{-m}.$ (4.18)

###### Proof.

For $j\in\mathbb{N}$, we define

$w^{j}=\mathbb{E}_{\mathcal{T}}^{j}\prod_{k=1}^{j}(1-\varepsilon\tau^{k}).$

Using (4.16), we obtain that

$w^{j}-w^{j-1}+\tau^{j}\mathbb{D}_{\mathcal{T}}^{j}\prod_{k=1}^{j-1}(1-\varepsilon\tau^{k})\leq
c_{\varepsilon,\mathbf{\Omega},\xi}\tau^{j}\prod_{k=1}^{j-1}(1-\varepsilon\tau^{k})$

and summing over $m$, we get

$w^{m}-w^{0}+\sum_{j=1}^{m}\tau^{j}\mathbb{D}_{\mathcal{T}}^{j}\prod_{k=1}^{j-1}(1-\varepsilon\tau^{k})\leq
c_{\varepsilon,\mathbf{\Omega},\xi}\sum_{j=1}^{m}\tau^{j}\prod_{k=1}^{j-1}(1-\varepsilon\tau^{k}).$

We can now multiply the last inequality by
$\prod_{k=1}^{m}(1-\varepsilon\tau^{k})^{-1}$ in order to come back to the
discrete entropy. It yields

$\mathbb{E}_{\mathcal{T}}^{m}+\sum_{j=1}^{m}\tau^{j}\mathbb{D}_{\mathcal{T}}^{j}\prod_{k=j}^{m}(1-\varepsilon\tau^{k})^{-1}\leq\prod_{k=1}^{m}(1-\varepsilon\tau^{k})^{-1}\mathbb{E}_{\mathcal{T}}^{0}+c_{\varepsilon,\Omega,\xi}\sum_{j=1}^{m}\tau^{j}\prod_{k=j}^{m}(1-\varepsilon\tau^{k})^{-1}.$

But, since $\varepsilon\Delta t<1$, we have

$1\leq\prod_{k=j}^{m}(1-\varepsilon\tau^{k})^{-1}\leq\prod_{k=1}^{m}(1-\varepsilon\tau^{k})^{-1}\leq(1-\varepsilon\Delta
t)^{-m},$

which yields (4.18) as $\sum_{j=1}^{m}\tau^{j}=t^{m}.$ ∎

##  5 Existence of a discrete solution

In this section, we will now establish the existence of a solution to the
finite volume scheme (4.1)–(4.8), which consists of a nonlinear system of
equations at each time step. Knowing the solution at step $m-1$, we want to
establish the existence of a solution at time step $m$. We may consider that
the unknowns of the nonlinear system of equations are the quasi Fermi
potentials and the electrostatic potential, as the densities of electrons,
holes and anion vacancies are defined as functions of these potentials through
(4.4). The proof consists of three main parts: we start in Section 5.1 showing
the existence and uniqueness of a discrete electric potential for given quasi
Fermi potentials associated to the Poisson equation and continue in Section
5.2 with proving some a priori estimates on the quasi Fermi and electrostatic
potentials, obtained as consequences of the bounds on the entropy and the
dissipation. Then, in Section 5.3 the existence of quasi Fermi potentials is
shown which finalizes the proof. For this, forgetting the superscript $m$, we
denote by $\mathbf{X}$ the vector containing the unknown quasi Fermi
potentials which is defined by

$\mathbf{X}=\Bigl{(}(\varphi_{n,K}-\varphi_{n,K}^{D})_{K\in\mathcal{T}},(\varphi_{p,K}-\varphi_{p,K}^{D})_{K\in\mathcal{T}},(\varphi_{a,K}-\psi_{K}^{D})_{K\in\mathcal{T}_{\text{intr}}}\Bigr{)}.$
(5.1)

###  5.1 Existence of electric potential

The aim of the first lemma is to show the existence of a unique
$\boldsymbol{\psi}=\left(\psi_{K}\right)_{K\in\mathcal{T}}$ dependent on
$\mathbf{X}$.

######  Lemma 5.1.

Let $\mathbf{X}$ denote the vector containing the unknown quasi Fermi
potentials as defined in (5.1). Then, there exists a unique solution
$\boldsymbol{\psi}(\mathbf{X})$ to the discrete nonlinear Poisson equation
(4.2). Further, the mapping $\mathbf{X}\mapsto\boldsymbol{\psi}(\mathbf{X})$
is continuous.

###### Proof.

Let us define the discrete functional

$\displaystyle\mathcal{J}(\boldsymbol{\Psi})=\,$
$\displaystyle\frac{\lambda^{2}}{2}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}|D_{\sigma}\boldsymbol{\Psi}|^{2}+\delta\\!\\!\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}m_{K}\mathcal{G}_{\alpha}(z_{\alpha}(\varphi_{\alpha,K}-\Psi_{K}))+\\!\\!\\!\sum_{K\in\mathcal{T}_{\text{intr}}}m_{K}\mathcal{G}_{a}(z_{a}(\varphi_{a,K}-\Psi_{K}))$
$\displaystyle-\delta\sum_{K\in\mathcal{T}}m_{K}C_{K}\Psi_{K},$

where $\mathcal{G}_{\alpha}$ denotes the primitive of $\mathcal{F}_{\alpha}$
which vanishes at $-\infty$. We can compute $\nabla\mathcal{J}$, where the
$K$-th component is given by

$\frac{\partial\mathcal{J}}{\partial\Psi_{K}}=\left\\{\begin{aligned}
&-\lambda^{2}\sum_{\sigma\in\mathcal{E}_{K}}\tau_{\sigma}D_{K,\sigma}\boldsymbol{\Psi}-\delta
m_{K}\Bigl{(}C_{K}+\\!\\!\sum_{\alpha\in\\{n,p\\}}z_{\alpha}\mathcal{F}_{\alpha}(z_{\alpha}(\varphi_{\alpha,K}-\Psi_{K}))\Bigr{)},&&K\in\mathcal{T}\setminus\mathcal{T}_{\text{intr}},\\\\[4.30554pt]
~{}&-\lambda^{2}\sum_{\sigma\in\mathcal{E}_{K}}\tau_{\sigma}D_{K,\sigma}\boldsymbol{\Psi}-\delta
m_{K}\Bigl{(}C_{K}+\\!\\!\sum_{\alpha\in\\{n,p\\}}z_{\alpha}\mathcal{F}_{\alpha}(z_{\alpha}(\varphi_{\alpha,K}-\Psi_{K}))\Bigr{)}\\\
&-m_{K}z_{a}\mathcal{F}_{a}(z_{a}(\varphi_{a,K}-\psi_{K})),&&K\in\mathcal{T}_{\text{intr}}.\end{aligned}\right.$

We conclude that a solution $\boldsymbol{\psi}$ to the discrete Poisson
equation (4.2) satisfies $\nabla\mathcal{J}(\boldsymbol{\psi})=\mathbf{0}$.
The existence of a global minimum of $\mathcal{J}$ is guaranteed through its
continuity and coercivity. The coercivity follows from the coercivity of
$\boldsymbol{\Psi}\mapsto\frac{\lambda^{2}}{2}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}|D_{\sigma}\boldsymbol{\Psi}|^{2}-\delta\sum_{K\in\mathcal{T}}m_{K}C_{K}\Psi_{K}$
(by a discrete Poincaré inequality [17]) and the boundedness from below of the
two other contributions. The strict convexity of $\mathcal{J}$, due to being a
sum of a strictly convex and convex functions, gives the uniqueness of this
global minimum $\boldsymbol{\psi}$. Lastly, the continuity of
$\mathbf{X}\mapsto\boldsymbol{\psi}(\mathbf{X})$ follows from the implicit
function theorem applied to $\nabla\mathcal{J}$ since the Hessian of
$\mathcal{J}$ with respect to $\boldsymbol{\Psi}$ is strictly row diagonally
dominant. ∎

As a consequence of Lemma 5.1 we can interpret in the following the electric
potential as a continuous map
$\boldsymbol{\psi}=\boldsymbol{\psi}(\mathbf{X})$.

###  5.2 A priori estimates

The discrete entropy defined by (4.14) can be denoted by
${\mathbb{E}}_{\mathcal{T}}(\mathbf{X})$ and its associated dissipation
defined by (4.15) can be denoted by ${\mathbb{D}}_{\mathcal{T}}(\mathbf{X})$.
In this dissipation, we may distinguish the contributions of electrons, holes,
anion vacancies and of the recombination-generation terms. Therefore, we
introduce the following notations

$\displaystyle{\mathbb{D}}_{\mathcal{T},a}(\mathbf{X})$
$\displaystyle=\frac{z_{a}^{2}}{2}\sum_{\sigma\in\mathcal{E}_{\rm intr}^{\rm
int}}\tau_{\sigma}\overline{n}_{a,\sigma}(D_{\sigma}\boldsymbol{\varphi}_{a})^{2},$
(5.2) $\displaystyle{\mathbb{D}}_{\mathcal{T},\alpha}(\mathbf{X})$
$\displaystyle=\frac{\delta}{2\nu}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}\overline{n}_{\alpha,\sigma}(D_{\sigma}\boldsymbol{\varphi}_{\alpha})^{2},\quad\text{for
all}\;\alpha\in\\{n,p\\}.$ (5.3)

In this section, the letter $R$ refers to a positive number, not to the
recombination term.

######  Lemma 5.2.

Assume that there exists $M_{E}>0$, such that
${\mathbb{E}}_{\mathcal{T}}(\mathbf{X})\leq M_{E}$. Then, there exists some
$R>0$ depending on $M_{E}$, $\lambda$ and on the mesh $\mathcal{T}$, such that

$-R\leq\psi_{K}-\psi_{K}^{D}\leq R,\quad\forall K\in\mathcal{T}.$ (5.4)

###### Proof.

As the entropy contributions of anion vacancies and of the relative entropies
for electrons and holes are non-negative, the bound on
${\mathbb{E}}_{\mathcal{T}}(\mathbf{X})$ directly implies a bound on the
electric energy,

$\frac{\lambda^{2}}{2}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}(D_{\sigma}(\boldsymbol{\psi}-\boldsymbol{\psi}^{D}))^{2}\leq
M_{E}.$

For each edge $\sigma\in\mathcal{E}^{D}$, such that
$\sigma\in\mathcal{E}_{K}$, we deduce a bound on $|\psi_{K}-\psi_{K}^{D}|$
depending on $M_{E}$, $\lambda$ and on the mesh. The same bound applies to
$|(\psi_{K}-\psi_{K}^{D})-(\psi_{L}-\psi_{L}^{D})|$ for any interior edge
$\sigma=K|L$. From these bounds, and by using the connectedness of the mesh
and the finite number of control volumes one can inductively get a uniform
finite bound for all $(\psi_{K}-\psi_{K}^{D})_{K\in\mathcal{T}}$. ∎

######  Lemma 5.3.

Assume that there exists $M_{D}>0$ such that
${\mathbb{D}}_{\mathcal{T},a}(\mathbf{X})\leq M_{D}$ and that there also
exists ${\bar{n}}\in(0,1)$, such that

$\displaystyle\frac{1}{|\mathbf{\Omega}|}\sum_{K\in\mathcal{T}_{\textrm{intr}}}m_{K}n_{a,K}={\bar{n}}.$
(5.5)

Then, there exists some $R>0$ depending on $M_{D}$, ${\bar{n}}$ and
$\mathcal{T}$, such that

$-R\leq\varphi_{a,K}\leq R,\quad\forall K\in\mathcal{T}_{\textrm{intr}}.$
(5.6)

Let us first note that due to hypothesis (H2) on ${\mathcal{F}}_{a}$, the
result stated in Lemma 5.3 is equivalent to the fact that there exists an
$\varepsilon\in(0,1)$ satisfying

$\varepsilon\leq n_{a,K}\leq 1-\varepsilon,\quad\forall
K\in\mathcal{T}_{\textrm{intr}}.$

This result is a direct consequence of [11, Lemma 3.2]. Its proof follows the
main lines of the proof of Lemma 3.7 in [11] and is left to the reader.
Lastly, we prove bounds on the quasi Fermi potentials of electric charge
carriers.

######  Lemma 5.4.

Let $\alpha\in\\{n,p\\}$. Assume that there exists $M_{E}>0$, such that
${\mathbb{E}}_{\mathcal{T}}(\mathbf{X})\leq M_{E}$ and $M_{D}>0$ such that
${\mathbb{D}}_{\mathcal{T},\alpha}(X)\leq M_{D}$. Then, there exists some
$R>0$ depending on $M_{E}$, $M_{D}$, $\mathbf{\Omega}$, $\mathcal{T}$,
$\psi^{D}$, $\varphi^{D}$, such that

$-R\leq\varphi_{\alpha,K},\leq R,\quad\forall K\in\mathcal{T}.$ (5.7)

###### Proof.

In order to prove Lemma 5.4, we will still stay close to the proof of Lemma
3.7 in [11]. It needs an adaptation of Lemma 3.2 in [11] due to the different
hypotheses on the statistics function ${\mathcal{F}}_{\alpha}$ and the
different kind of boundary conditions.

Let us first rewrite ${\mathbb{D}}_{\mathcal{T},\alpha}(\mathbf{X})$ by using
the reformulation of the fluxes (4.13) based on the definition (4.11) of
${\overline{n}}_{\alpha,\sigma}$

${\mathbb{D}}_{\mathcal{T},\alpha}(\mathbf{X})=-\frac{\delta}{2\nu}\sum_{\sigma\in\mathcal{E}}J_{\alpha,K,\sigma}D_{K,\sigma}{\boldsymbol{\varphi}}_{\alpha},~{}$

where the flux discretization is defined through (4.6) and (4.8). Introducing
the function ${\mathcal{K}}_{\alpha}:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$
defined by

${\mathcal{K}}_{\alpha}(x,a)=\log({\mathcal{F}}_{\alpha}(x-a))-x,\quad\forall(x,a)\in\mathbb{R}\times\mathbb{R},$

we note that

$Q_{\alpha,K,\sigma}={\mathcal{K}}_{\alpha}(z_{\alpha}\varphi_{\alpha,K},z_{\alpha}\psi_{K})-{\mathcal{K}}_{\alpha}(z_{\alpha}\varphi_{\alpha,K,\sigma},z_{\alpha}\psi_{K,\sigma}),$

where $\varphi_{\alpha,K,\sigma}$ and $\psi_{K,\sigma}$ stand for
$\varphi_{\alpha,L},\psi_{L}$, if $\sigma=K|L\in\mathcal{E}^{\textrm{int}}$
and for $\varphi_{\sigma}^{D},\psi_{\sigma}^{D}$, if
$\sigma\in\mathcal{E}^{D}$. Thus,
${\mathbb{D}}_{\mathcal{T},\alpha}(\mathbf{X})$ can be rewritten as

${\mathbb{D}}_{\mathcal{T},\alpha}(\mathbf{X})=\frac{\delta}{2\nu}\sum_{\sigma\in\mathcal{E}}\tau_{\sigma}{\mathcal{D}}_{\alpha}\bigl{(}z_{\alpha}\varphi_{\alpha,K},z_{\alpha}\varphi_{\alpha,K,\sigma},z_{\alpha}\psi_{K},z_{\alpha}\psi_{K,\sigma}\bigr{)},$

with ${\mathcal{D}}_{\alpha}:\mathbb{R}^{4}\to\mathbb{R}$ defined by

${\mathcal{D}}_{\alpha}(x,y,a,b)=(x-y)\Bigl{[}B\Bigl{(}{\mathcal{K}}_{\alpha}(x,a)-{\mathcal{K}}_{\alpha}(y,b)\Bigr{)}{\mathcal{F}}_{\alpha}(x-a)-B\Bigl{(}{\mathcal{K}}_{\alpha}(y,b)-{\mathcal{K}}_{\alpha}(x,a)\Bigr{)}{\mathcal{F}}_{\alpha}(y-b)\Bigr{]}.$

Following the strategy of proof of Lemma 3.7 in [11], we introduce
$\Upsilon_{\overline{\Phi},\overline{\Psi}}:\mathbb{R}\to\mathbb{R}$ defined
by

$\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)=\inf\left\\{{\mathcal{D}}_{\alpha}(x,y,a,b);\
-\overline{\Phi}\leq y\leq\overline{\Phi},-\overline{\Psi}\leq
a,b\leq\overline{\Psi}\right\\}$

and we establish (see Appendix B, Lemma B.1) that

$\lim_{x\to-\infty}\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)=+\infty\quad\text{and}\quad\lim_{x\to+\infty}\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)=+\infty.$
(5.8)

Then, we use that the discrete values of the electrostatic potential are
bounded thanks to Lemma 5.2 and that the Dirichlet boundary conditions ensure
that there exists at least one
$\varphi_{\alpha,K,\sigma}=\varphi_{\sigma}^{D}$ which is bounded. Lastly, the
bound on ${\mathbb{D}}_{\mathcal{T},\alpha}(\mathbf{X})$ implies that the
value $\varphi_{\alpha,K}$ is also bounded thanks to (5.8), where this
property propagates from cell to cell, such that it holds on the whole domain.
∎

###  5.3 Existence of quasi Fermi potentials

Finally, we can formulate and prove the existence of discrete solutions in
Theorem 5.5.

######  Theorem 5.5.

For all $m\geq 1$, the finite volume scheme (4.1)–(4.8) for the perovskite
model has at least one solution
$(\boldsymbol{\varphi}_{n}^{m},\boldsymbol{\varphi}_{p}^{m},\boldsymbol{\varphi}_{a}^{m},\boldsymbol{\psi}^{m})\in\mathbb{R}^{\theta}$
with
$\theta=3\textrm{Card}(\mathcal{T})+\textrm{Card}(\mathcal{T}_{\textrm{intr}})$.
Moreover, this solution satisfies the following $L^{\infty}$ bounds. There
exists $R>0$ depending on the data and on the mesh such that

$-R\leq\boldsymbol{\varphi}_{n}^{m},\boldsymbol{\varphi}_{p}^{m},\boldsymbol{\varphi}_{a}^{m},\boldsymbol{\psi}^{m}\leq
R,\quad\text{for all}\;m\geq 1,$

holds component-wise.

The discrete mass balances in (4.1) at step $m$ constitute a nonlinear system
of equations. More precisely, we can introduce a continuous vector field
$\mathbf{P}_{m}:\mathbb{R}^{\theta_{\mathbf{X}}}\to\mathbb{R}^{\theta_{\mathbf{X}}}$
with
$\theta_{\mathbf{X}}=2\textrm{Card}(\mathcal{T})+\textrm{Card}(\mathcal{T}_{\textrm{intr}})$,
such that $\mathbf{P}_{m}(\mathbf{X}^{m})=\mathbf{0}$ is equivalent to (4.1),
where $\mathbf{X}^{m}$ is defined by (5.1), noting that we have omitted the
superscript $m$ there. In order to formulate the electron and hole components
of $\mathbf{P}_{m}(\mathbf{X}^{m})$, we put every term of the equations (4.1a)
and (4.1b) on the left-hand side and rescale by a factor $\delta\tau^{m}/\nu$.
The anion related components are given by (4.1c) rescaled by $\tau^{m}$. In
order to prove Theorem 5.5, we apply a corollary of Brouwer’s fixed point
theorem [16, Section 9.1] to a regularized version of $\mathbf{P}_{m}$, and
then take limits in the regularization parameter. The fixed point lemma reads
as follows.

######  Lemma 5.6.

Let $N\in\mathbb{N}$ and $\mathbf{P}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ be a
continuous vector field. Assume that there exists $R>0$, such that
$\mathbf{P}(\mathbf{X})\cdot\mathbf{X}\geq 0$, if $\|\mathbf{X}\|=R$. Then,
there exists $\mathbf{X}^{\ast}\in\mathbb{R}^{N}$ such that
$\mathbf{P}(\mathbf{X}^{\ast})=\mathbf{0}$ and $\|\mathbf{X}\|\leq R$.

###### Proof of Theorem 5.5.

First, we prove the existence of quasi Fermi potentials $\mathbf{X}$, where
for the sake of readability, we omit the superscript $m$. We recall that Lemma
5.1 guarantees the existence of a continuous and uniquely determined map
$\mathbf{X}\mapsto\boldsymbol{\psi}(\mathbf{X})$ solving the nonlinear Poisson
equation (4.2) for any given quasi Fermi potentials $\mathbf{X}$. Thus,
$\mathbf{P}_{m}$ is well-defined and continuous. The scalar product
$\mathbf{P}_{m}(\mathbf{X})\cdot\mathbf{X}$ is given by

$\displaystyle\mathbf{P}_{m}(\mathbf{X})\cdot\mathbf{X}=$
$\displaystyle~{}\sum_{\alpha\in\\{n,p\\}}\sum_{K\in\mathcal{T}}\left(\delta
z_{\alpha}m_{K}(n_{\alpha,K}-n_{\alpha,K}^{m-1})(\varphi_{\alpha,K}-\varphi_{K}^{D})~{}+\frac{\delta\tau^{m}}{\nu}\sum_{\sigma\in\mathcal{E}_{K}}J_{\alpha,K,\sigma}(\varphi_{\alpha,K}-\varphi_{K}^{D})\right)$
$\displaystyle+\frac{\delta\tau^{m}}{\nu}\sum_{K\in\mathcal{T}}m_{K}R(n_{n,K},n_{p,K})(\varphi_{p,K}-\varphi_{n,K})-\frac{\delta\gamma\tau^{m}}{\nu}\sum_{K\in\mathcal{T}}m_{K}G_{K}(\varphi_{p,K}-\varphi_{n,K})$
$\displaystyle+\sum_{K\in\mathcal{T}}z_{a}m_{K}(n_{a,K}-n_{a,K}^{m-1})(\varphi_{a,K}-\psi_{K}^{D})+\tau^{m}\sum_{K\in\mathcal{T}}\sum_{\sigma\in\mathcal{E}_{K}}J_{a,K,\sigma}(\varphi_{a,K}-\psi_{K}^{D}).$

We have established the following inequality within the proof of Theorem 4.4
for $\varepsilon>0$

$\mathbf{P}_{m}(\mathbf{X})\cdot\mathbf{X}\geq(1-\varepsilon\tau^{m}){\mathbb{E}}_{\mathcal{T}}(\mathbf{X})-{\mathbb{E}}_{\mathcal{T}}(\mathbf{X}^{m-1})+\tau^{m}{\mathbb{D}}_{\mathcal{T}}(\mathbf{X})-\tau^{m}c_{\varepsilon,\mathbf{\Omega},\xi},$

where $\mathbf{X}^{m-1}$ denotes the known solution at the previous time step
$m-1$. For suitable $\varepsilon$, there exists $M>0$, such that

$\mathbf{P}_{m}(\mathbf{X})\cdot\mathbf{X}\geq\frac{1}{2}{\mathbb{E}}_{\mathcal{T}}(\mathbf{X})+\tau^{m}{\mathbb{D}}_{\mathcal{T}}(\mathbf{X})-M.$

Our goal is to use Lemma 5.6 to show the existence of a solution at time step
$t^{m}$. Instead of showing now the non-negativity of the scalar product
$\mathbf{P}_{m}(\mathbf{X})\cdot\mathbf{X}$, we introduce a parameter-
dependent regularization of $\mathbf{P}_{m}$ which satisfies the assumptions
of Lemma 5.6.

For a given $\mu>0$, we define
$\mathbf{P}_{m}^{\mu}(\mathbf{X})=\mathbf{P}_{m}(\mathbf{X})+\mu\mathbf{X}$,
which satisfies

$\mathbf{P}_{m}^{\mu}(\mathbf{X})\cdot\mathbf{X}=\mathbf{P}_{m}(\mathbf{X})\cdot\mathbf{X}+\mu\|\mathbf{X}\|^{2}\geq\mu\|\mathbf{X}\|^{2}-M\>\geq
0,\quad\text{for}\>\|\mathbf{X}\|\geq\sqrt{M/\mu}.$

Then, Lemma 5.6 shows the existence of $\mathbf{X}^{m,\mu}\in
B(\mathbf{0},\sqrt{M/\mu})$, such that
$\mathbf{P}_{m}^{\mu}(\mathbf{X}^{m,\mu})=\mathbf{0}$. Next, we need to show
$\mathbf{X}^{m,\mu}$ is actually uniformly bounded in $\mu$. Let us check the
hypotheses of Lemmas 5.2, 5.3 and 5.4. We take the scalar product of
$\mathbf{P}_{m}^{\mu}(\mathbf{X}^{m,\mu})$ with the vector
$\mathbf{V}=({\mathbf{0}}_{\mathcal{T}},{\mathbf{0}}_{\mathcal{T}},{\mathbf{1}}_{\mathcal{T}_{\rm
intr}})$. Since the sum over all fluxes in the intrinsic region vanishes, we
obtain

$\sum_{K\in\mathcal{T}_{\rm
intr}}z_{a}m_{K}n_{a,K}^{m,\mu}-\sum_{K\in\mathcal{T}_{\rm
intr}}z_{a}m_{K}n_{a,K}^{m-1}+\mu\mathbf{X}^{m,\mu}\cdot\mathbf{V}=0,$

and therefore, after rescaling with the measure of $\mathbf{\Omega}$, we have

$\left|\frac{1}{|\mathbf{\Omega}|}\sum_{K\in\mathcal{T}_{\rm
intr}}z_{a}m_{K}n_{a,K}^{m,\mu}-\frac{1}{|\mathbf{\Omega}|}\sum_{K\in\mathcal{T}_{\rm
intr}}z_{a}m_{K}n_{a,K}^{m-1}\right|\leq\frac{\mu}{|\mathbf{\Omega}|}\|\mathbf{X}^{m,\mu}\|\|\mathbf{V}\|\leq\frac{\sqrt{M\mu}}{|\mathbf{\Omega}|}\|\mathbf{V}\|.$
(5.9)

But since the solution at the previous time step exists and hence is bounded,
there exists $\varepsilon^{(m-1)}\in(0,1)$, such that
$\frac{1}{|\mathbf{\Omega}|}\sum_{K\in\mathcal{T}_{\rm
intr}}m_{K}n_{a,K}^{m-1}\in(\varepsilon^{(m-1)},1-\varepsilon^{(m-1)})$. Thus,
we deduce from (5.9) that, for $\mu$ sufficiently small, $\mathbf{X}^{m,\mu}$
satisfies for $\varepsilon^{(m)}=\varepsilon^{(m-1)}/2$

$\frac{1}{|\mathbf{\Omega}|}\sum_{K\in\mathcal{T}_{\rm
intr}}m_{K}n_{a,K}^{m,\mu}\in(\varepsilon^{(m)},1-\varepsilon^{(m)}).$

Moreover, as
$\mathbf{P}_{m}^{\mu}(\mathbf{X}^{m,\mu})\cdot\mathbf{X}^{m,\mu}=0\geq\frac{1}{2}{\mathbb{E}}_{\mathcal{T}}(\mathbf{X}^{m,\mu})+\tau^{m}{\mathbb{D}}_{\mathcal{T}}(\mathbf{X}^{m,\mu})-M$,
we see that ${\mathbb{E}}_{\mathcal{T}}(\mathbf{X}^{m,\mu})$ and
${\mathbb{D}}_{\mathcal{T}}(\mathbf{X}^{m,\mu})$ are uniformly bounded in
$\mu$ by $M$. Hence, we can apply Lemmas 5.2, 5.3, 5.4 to deduce that
$\|\mathbf{X}^{m,\mu}\|$ is bounded uniformly in $\mu$ (for $\mu$ sufficiently
small). Finally, we can extract a subsequence, which converges to a limit
denoted by $\mathbf{X}^{m}$ as $\mu$ tends to $0$. This limit satisfies
$\mathbf{P}_{m}^{0}(\mathbf{X}^{m})=\mathbf{P}_{m}(\mathbf{X}^{m})=\mathbf{0}$.
Thus, we have found quasi Fermi potentials which solve the discrete system
(4.1). It remains to show the existence of a uniquely determined
$\boldsymbol{\psi}(\mathbf{X}^{m})$ which solves (4.2). However, this follows
from Lemma 5.1, which ends the proof of Theorem 5.5.

∎

##  6 Numerical experiments

The numerical examples were performed with ChargeTransport.jl, a Julia package
for the simulation of charge transport in semiconductors [2]. In a first step
the aim is to verify properties of the finite volume scheme (4.1)-(4.8) such
as a special case of the entropy-dissipation inequality in Theorem 4.4 as well
as the spatial convergence rate. In a second step, the charge transport model
(2.13)-(2.20) is simulated for a physical meaningful set of parameters. In all
simulation setups we are interested in the large time behavior of the model.
For this reason, we introduce an entropy with respect to the steady state

$\mathbb{E}_{\infty}(t)=\frac{\lambda^{2}}{2}\int_{\mathbf{\Omega}}|\nabla(\psi-\psi^{\infty})|^{2}\,d\mathbf{x}+\int_{\mathbf{\Omega}_{\text{intr}}}H_{a}(n_{a},n^{\infty}_{a})\,d\mathbf{x}+\delta\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}H_{\alpha}(n_{\alpha},n_{\alpha}^{\infty})\,d\mathbf{x},$
(6.1)

where $H_{\alpha}$ is defined in (3.2). The non-negative functional
$E_{\infty}$ can be seen as a measure of the distance between a solution at
time $t$ and the steady state of the model which vanishes, if and only if the
solution at time $t$ and the steady state coincide almost everywhere.
Furthermore, from an analytical point of view $E_{\infty}$ may help to prove
the convergence of the discrete solution to the discrete steady state [12].

###  6.1 Verifying the properties of the scheme

Within this section we assume a one-dimensional domain $\Omega=(0,6)$ and set
$\Omega_{\text{HTL}}=(0,2),\Omega_{\text{intr}}=(2,4),\Omega_{\text{ETL}}=(4,6)$.
We choose $513$ nodes per subdomain, resulting in a total number of $1537$
nodes with a grid spacing $h\approx$3.9\text{\times}{10}^{-3}$$. The time
domain is given by $[0,80]$ which we discretize with a time step of $\Delta
t=$1.0\text{\times}{10}^{-1}$$. We set the rescaled Debye length to
$\lambda=1$, the relative mobility of anion vacancies to $\nu=1$, the relative
concentration to $\delta=1$, and the rescaled photogeneration rate to
$\gamma=1$.

#### Thermal Equilibrium boundary conditions

Let us first study the implications of the assumptions in Remark 3.4. To this
end, we assume a constant doping $C=0.1$ and no generation and recombination,
i.e. $G=R=0$. The Dirichlet functions (2.6a) are chosen as constant functions
$\varphi^{D}=0.5$ and $\psi^{D}=\text{arcsinh}(C/2)+0.5$. The sinusoidal
initial conditions for electrons, holes and the electric potential as well as
the constant initial condition for anion vacancies along with the steady state
solutions are depicted in Figure 6.1 on the left panel. On the right panel we
show the steady state densities.

Figure 6.1: Steady state solutions
$(\varphi_{n}^{\infty},\varphi_{p}^{\infty},\varphi_{a}^{\infty},\psi^{\infty})$
with the respective initial conditions as dotted lines (left) and the
associated steady state densities of charge carriers (right) calculated via
(2.3).

Since for these specific choices, we have
$\mathbf{0}=\nabla\varphi^{D}=\nabla\psi^{D}$ and $G=0$ the discrete entropy-
dissipation inequality in Theorem 4.4 indicates that the relative entropy with
respect to the Dirichlet boundary values (4.14) does not increase in time.
This result can be numerically verified, see Figure 6.2. Due to a non-constant
electric potential $\psi^{\infty}$ we observe that the relative entropy (4.14)
(in blue on the left panel) levels off after an initial decrease. Furthermore,
the relative entropy with respect to the steady state (6.1) (in green on the
left panel of Figure 6.2) as well as the quadratic $L^{2}$ errors between the
steady state and a solution at time $t$ (right panel) decay exponentially with
a similar slope, reaching machine precision at a similar time.

Figure 6.2: Left: Time evolution of the relative entropy with respect to the
Dirichlet boundary functions (3.4) as well as the relative entropy with
respect to the steady state (6.1). Right: Time evolution of the quadratic
$L^{2}$ errors between the computed and the steady state solutions.

#### Non-constant boundary values

Next, we adjust the doping and the boundary values. Let us assume that the
doping $C$ is a piecewise constant function given by $0.5$ in
$\Omega_{\text{ETL}}$ and by $-0.5$ in
$\Omega_{\text{HTL}}\cup\Omega_{\text{intr}}$. The boundary values are set to
$\varphi|_{x=0}^{D}=1,\varphi|_{x=6}^{D}=0$,
$\quad\psi|_{x=0}^{D}=\text{arcsinh}(-0.5/2)+1,\psi|_{x=6}^{D}=\text{arcsinh}(0.5/2)$.
We choose quadratic initial conditions for $\varphi_{n},\varphi_{p},\psi$ and
a constant initial condition for $\varphi_{a}$. The initial conditions are
additionally to the steady state solutions depicted in Figure 6.3 as dotted
lines.

Figure 6.3: Steady state potentials with the corresponding initial conditions
(left) and the associated steady state densities of charge carriers (right)
calculated via the state equation (2.3).

Again, the relative entropy with respect to the steady state and the quadratic
$L^{2}$ errors decay exponentially and reach machine precision with a similar
slope, see Figure 6.4 in the left and middle panel.

Finally, we complete this section with an investigation of the spatial
convergence behavior. Suppose $n_{*}\in\mathbb{N}$ is given, then $n=2\cdot
2^{n_{*}-1}+1$ nodes are chosen in each of the three subdomain, i.e. in total
we have $n_{\text{tot}}=3\cdot 2\cdot 2^{n_{*}-1}+1$ nodes. We calculate a
reference solution on a grid with $n_{*}=9$ corresponding to $1537$ nodes with
a grid spacing $h_{*}\approx$3.9\text{\times}{10}^{-3}$$. The $L^{2}$ errors
between the solution $u_{n_{*}}$, for $n_{*}=2,\ldots,8$, and the reference
solution projected onto the coarser mesh evaluated at the final time are shown
in Figure 6.4, right panel. Since for the final time $t_{F}=80$ the system is
already within machine precision of the steady state, the error shown is
purely due to the spatial discretization. We observe second order experimental
convergence.

Figure 6.4: Time evolution of the relative entropy with respect to steady
state (6.1) for non-constant boundary values (left) and of the quadratic
$L^{2}$ errors between steady state and solutions at time $t$ (middle). On the
right, the $L^{2}$ error with respect to the grid spacing $h$ is shown.

###  6.2 PSC simulation setup

In the final simulation setup, we choose the rescaling factors and non-
dimensionalized parameters in such a way that the resulting solutions
correspond to a realistic PSC device. All parameters are chosen in agreement
with Section 2.4 and with non-zero band-edge energies. Apart from the
additional non-scaled parameters
$N_{a}=$1.0\text{\times}{10}^{21}$\,\text{cm}^{-3}$ and
$E_{a}=-4.45\;\text{eV}$, we used the parameter set provided in [14]. The mesh
is given by $385$ nodes with a uniform grid spacing in each layer, namely
$h_{\text{ETL}}\approx$7.8\text{\times}{10}^{-8}$$cm,
$h_{\text{HTL}}\approx$1.6\text{\times}{10}^{-7}$$cm in the transport layers
and $h_{\text{intr}}\approx$3.1\text{\times}{10}^{-7}$$cm in the perovskite
layer. The uniform time mesh is built with a step size of $\Delta
t=0.5\,\text{s}$ and the final time is given by $t_{F}=220\,\text{s}$.

Figure 6.5:  Steady state potentials with the respective initial conditions
(left) and the associated densities of charge carriers (right) for a PSC
three-layer device at an applied voltage $0.9\,\text{V}$.

Usually a PSC device is held for several seconds at a constant voltage,
ensuring that ionic charges equilibrate. This procedure is often called
preconditioning protocol [13]. Afterwards, scan protocols with a time-
dependent applied voltage – incorporated via time-dependent Dirichlet boundary
conditions – are performed to study the device physics. Thus, the steady
potentials and their respective densities depicted in Figure 6.5 can be
regarded as the solutions after a successful preconditioning scan. Within the
presented configuration the applied voltage is chosen such that the steady
state electric potential $\psi^{\infty}$ is constant which can be observed
well in Figure 6.5. The depicted initial conditions correspond to a solution
of the charge transport model with a non-constant vacancy concentration. As
before, we consider the large time behavior of the quadratic $L^{2}$ errors
and the relative entropy with respect to the steady state. Taking the
thermodynamic free energy (2.12) into account, we can reformulate the
dimensional relative entropy with respect to the steady state (6.1)

$\mathbb{E}_{\infty}(t)=\frac{1}{2}\int_{\mathbf{\Omega}}\varepsilon_{s}|\nabla(\psi-\psi^{\infty})|^{2}\,d\mathbf{x}+\int_{\mathbf{\Omega}_{\text{intr}}}H_{a}(n_{a},n^{\infty}_{a})\,d\mathbf{x}+\sum_{\alpha\in\\{n,p\\}}\int_{\mathbf{\Omega}}H_{\alpha}(n_{\alpha},n_{\alpha}^{\infty})\,d\mathbf{x},$
(6.2)

where
$H_{\alpha}(x,y)=\Phi_{\alpha}(x)-\Phi_{\alpha}(y)-\Phi_{\alpha}^{\prime}(y)(x-y)$,
as defined in (3.2), but with

$\displaystyle\Phi_{n}(x)=k_{B}Tx$
$\displaystyle\left(\log\left(\frac{x}{N_{n}}\right)-1\right)-z_{n}E_{n}x,\quad\Phi_{p}(x)=k_{B}Tx\left(\log\left(\frac{x}{N_{p}}\right)-1\right)-z_{p}E_{p}x,\quad$
(6.3) $\displaystyle~{}\Phi_{a}(x)$
$\displaystyle=k_{B}T\left(x\log\left(\frac{x}{N_{a}}\right)+\left(N_{a}-x\right)\log\left(1-\frac{x}{N_{a}}\right)\right)-z_{a}E_{a}x,$
(6.4)

i.e. we extend the contributions of the relative entropy with respect to the
steady state such that they are consistent with the thermodynamic free energy
(2.12). As before, Figure 6.6 indicates an exponential decay towards zero of
the relative entropy (6.2) as well as of the quadratic $L^{2}$ errors with
respect to time. In contrast to the observations made in previous section, the
relative entropy with respect to the steady state (6.2) vanishes faster than
the quadratic $L^{2}$ errors. This may be explained by the additional terms in
(6.3) due to non-zero band-edge energies which influence the convergence
behavior. Still, we see a similar convergence rate of the two introduced
measures for the deviation of a solution at time $t$ from the steady state.

Figure 6.6: Time evolution of the relative entropy with respect to the steady
state (6.2) (left) and of the quadratic $L^{2}$ errors between the computed
and the steady state solutions (right).

##  7 Conclusion and outlook

For a charge transport model for perovskite solar cells, we discussed and
proved a continuous entropy-dissipation inequality. We allowed general
statistics function for the electric charge carriers. Moreover, we proved an
analogous entropy-dissipation inequality for a finite volume scheme based on
the excess chemical potential flux. The entropy-dissipation inequality helped
us to prove the existence of a discrete solution at every time step.
Furthermore, for a model in thermodynamic equilibrium we proved the decay of
the continuous and discrete relative entropy with respect to the boundary
conditions and numerically verified this result. A spatial convergence of
order $2$ was also numerically shown. In the last experiment, we studied the
numerical convergence towards the steady state for a setup which can be
physically interpreted as preconditioning a PSC device before applying a
measurement protocol. Especially the relative entropy with respect to the
steady state for non-zero band-edge energies decays exponentially towards
zero. Studying the model behavior with non-zero and even irregular band-edge
energies is from a mathematical and physical point of view of interest in the
future. Also time-dependent Dirichlet functions coinciding with physically
realistic measurement techniques for perovskite solar cells can be a topic of
future research. Finally, deriving reduced models by ignoring small
dimensionless parameters could also be investigated.

Acknowledgments This work was partially supported by the Leibniz competition
as well as the CEMPI (ANR-11-LABX-0007) and the scientific department of the
French Embassy in Germany.

##  Appendix A Estimates on statistics and entropy functions

In this section we deal with the estimates concerning statistics and entropy
functions. First, we provide the proof of Lemma 3.1 and Lemma 3.2 which are
stated under abstract assumptions (H1), (H2) and/or (H3) on the statistics.
Then, we show that the examples of Fermi-Dirac and Boltzmann statistics (2.9),
(2.10) and (2.11) satisfy these assumptions.

###  A.1 Proof of Lemma 3.1

Let us show point (i) of the Lemma. Let $\mathcal{F}_{\alpha}$ with
$\alpha=n,p$ be a statistics function satisfying (H1) and $H_{\alpha}$ be the
associated relative entropy function. Let $\varepsilon>0$ and $y_{0}\geq 0$.
For $x\geq 0$ and $y\in[0,y_{0}]$ one has

$\displaystyle x\leq\sup_{x\in\mathbb{R}}\left(x-\varepsilon
H_{\alpha}(x,y)\right)+\varepsilon H_{\alpha}(x,y).$

The first term on the right-hand side is the Legendre transform of
$x\mapsto\varepsilon H_{\alpha}(x,y)$ evaluated at $1$. It is exactly given by
$\bar{x}-\varepsilon H_{\alpha}(\bar{x},y)$ with
$\bar{x}=\mathcal{F}_{\alpha}\left(\frac{1}{\varepsilon}+\mathcal{F}_{\alpha}^{-1}(y)\right)$.
In turn, one has

$\displaystyle x\leq c_{y_{0},\varepsilon}+\varepsilon
H_{\alpha}(x,y),\quad\text{for all}\;x\geq 0,\ y\in[0,y_{0}],$

for
$c_{y_{0},\varepsilon}=\mathcal{F}_{\alpha}\left(\frac{1}{\varepsilon}+\mathcal{F}_{\alpha}^{-1}(y_{0})\right)$,
since $\mathcal{F}_{\alpha}$ is increasing.

For point (ii), where $\mathcal{F}_{a}$ is a statistics function satisfying
(H2) and $\Phi_{a}$ is the associated entropy function, an analogous
calculation will prove the estimate with
$c_{\varepsilon}=\mathcal{F}_{a}(\frac{1}{\varepsilon})$

$\displaystyle x\leq c_{\varepsilon}+\varepsilon\Phi_{a}(x),\quad\text{for
all}\;x\geq 0.$

###  A.2 Proof of Lemma 3.2

Now, let us consider $\mathcal{F}_{\alpha}$ with $\alpha=n,p$ as a statistics
function satisfying (H1) and (H3). Let $\varepsilon>0$ and $y_{0}\geq 0$. On
the one hand, because of (H3) there exists $x_{0}=x_{0}(y_{0},\varepsilon)\geq
0$ such that

$\displaystyle 0\leq\mathcal{F}_{\alpha}^{-1}(x)\leq\varepsilon
H_{\alpha}(x,y_{0}),\quad\text{for all}\;x\geq x_{0}.$

Further, the calculation
$\partial_{y}H_{\alpha}(x,y)=-\Phi_{\alpha}^{\prime\prime}(y)(x-y)$ reveals
that $y\mapsto H_{\alpha}(x,y)$ is non-increasing for all $y\in[0,y_{0}],y\leq
x$, due to the convexity of $\Phi_{\alpha}$. Hence,

$\displaystyle 0\leq\mathcal{F}_{\alpha}^{-1}(x)\leq\varepsilon
H_{\alpha}(x,y)\leq\varepsilon H_{\alpha}(x,y_{0}),\quad\text{for
all}\;y\in[0,y_{0}],\;x\geq\max\\{x_{0},y_{0}\\}=:\bar{x}.$

On the other hand, for $0\leq x\leq\bar{x}$

$\displaystyle\mathcal{F}_{\alpha}^{-1}(x)$
$\displaystyle\leq\mathcal{F}_{\alpha}^{-1}(\bar{x}).$

Hence, in total the claim is proven with
$c_{y_{0},\varepsilon}=\mathcal{F}_{\alpha}^{-1}(\bar{x})$

$\displaystyle\max(\mathcal{F}_{\alpha}^{-1}(x),0)\leq
c_{y_{0},\varepsilon}+\varepsilon H_{\alpha}(x,y),\quad\text{for all}\;x\geq
0,\quad y\in[0,y_{0}].$

###  A.3 Boltzmann and Fermi-Dirac entropy functions

We now relate the statistics (2.9), (2.10) to the hypotheses (H1) and (H3),
and the statistics (2.11) to the hypothesis (H2). We only give a proof in the
case of the Fermi-Dirac statistics of order $1/2$ since the results are
essentially trivial for the other statistics.

######  Lemma A.1 (Boltzmann).

Assume that the statistics for electrons and holes is the Boltzmann statistics
(2.10), namely

$\mathcal{F}_{\alpha}(\eta)=F_{B}(\eta)=e^{\eta}\,,\quad\eta\in\mathbb{R}\,,\
\alpha\in\\{n,p\\}.$

Then, $\mathcal{F}_{\alpha}$ satisfies (H1) and (H3).

######  Lemma A.2 (Fermi-Dirac of order $-1$).

Assume that the statistics for anion vacancies is the Fermi-Dirac integral of
order $-1$ in (2.11), namely

$\mathcal{F}_{a}(\eta)=F_{-1}(\eta)=\frac{1}{\exp(-\eta)+1}\,,\quad\eta\in\mathbb{R}\,.$

Then, $\mathcal{F}_{a}$ satisfies (H2).

######  Lemma A.3 (Fermi-Dirac of order $1/2$).

Assume that the statistics for electrons and holes is the Fermi-Dirac integral
of order $1/2$ in (2.9), namely

$\mathcal{F}_{\alpha}(\eta)=F_{1/2}(\eta)=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}\frac{\xi^{1/2}}{\exp(\xi-\eta)+1}\>\mathrm{d}\xi\,,\quad\eta\in\mathbb{R}\,,\
\alpha\in\\{n,p\\}.$

Then, $\mathcal{F}_{\alpha}$ satisfies (H1) and (H3).

###### Proof.

First, observe that $\mathcal{F}_{\alpha}$ is smooth and strictly increasing
with limits $0$ and $+\infty$, when $\eta\to-\infty$ and $\eta\to+\infty$
respectively. Then,

$F_{1/2}(\eta)\exp(-\eta)=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}\frac{\xi^{1/2}}{\exp(\xi)+\exp(\eta)}\>\mathrm{d}\xi\leq\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}\xi^{1/2}\exp(-\xi)\>\mathrm{d}\xi=1.$

Moreover,

$F_{1/2}^{\prime}(\eta)=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}\frac{\xi^{1/2}\exp(\xi-\eta)}{(\exp(\xi-\eta)+1)^{2}}\>\mathrm{d}\xi\leq
F_{1/2}(\eta).$

This proves (H1). Now let us focus on the behavior at infinity of $F_{1/2}$.
We claim the existence of constants $c_{1},c_{2}>0$, such that

$\displaystyle c_{1}\eta^{3/2}\leq F_{1/2}(\eta)\leq
c_{2}\eta^{3/2},\quad\text{for}\;\eta\geq 1.$ (A.1)

With (A.1) we can conclude $F_{1/2}^{-1}(s)=\mathcal{O}(s^{2/3})$ for
$s\to\infty$. Therefore, the associated entropy function behaves like
$\mathcal{O}(s^{5/3})$ and (H3) readily follows. To see that (A.1) is indeed
satisfied, let us consider (2.9) on the two intervals $[0,\eta]$ and
$[\eta,\infty)$, where the respective integrals are denoted by $I_{1},I_{2}$
and substitute $z=\xi-\eta$. This yields
$F_{1/2}(\eta)=\frac{2}{\sqrt{\pi}}(I_{1}+I_{2})$ with

$I_{1}=\int_{0}^{\eta}\frac{\xi^{1/2}}{\exp\left(\xi-\eta\right)+1}\,d\xi\quad\mbox{
and }\quad I_{2}=\int_{0}^{\infty}\frac{(z+\eta)^{1/2}}{\exp z+1}\,dz.$

We bound $I_{1}$ and $I_{2}$ separately. On the one hand, since
$0\leq\exp\left(\xi-\eta\right)\leq 1$ for $\xi\leq\eta$, we obtain

$\displaystyle\frac{1}{3}\eta^{3/2}\leq I_{1}\leq\frac{2}{3}\eta^{3/2}.$

On the other hand, we split $I_{2}$ into an integral over $[0,\eta]$ and one
over $[\eta,+\infty)$ and bound each term

$I_{2}=\int_{0}^{\eta}\frac{(z+\eta)^{1/2}}{\exp
z+1}\,dz+\int_{\eta}^{\infty}\frac{(z+\eta)^{1/2}}{\exp z+1}\,dz\leq
2^{1/2}\eta^{1/2}\int_{0}^{\eta}\frac{1}{\exp
z+1}\,dz+2^{1/2}\int_{\eta}^{\infty}\frac{z^{1/2}}{\exp z+1}\,dz.$

But,

$\int_{0}^{\eta}\frac{1}{\exp z+1}\,dz\leq\int_{0}^{\eta}\exp(-z)dz\leq
1\quad\mbox{ and }\quad\int_{\eta}^{\infty}\frac{z^{1/2}}{\exp
z+1}\,dz\leq\int_{0}^{\infty}\frac{z^{1/2}}{\exp z}\,dz=\Gamma(3/2),$

where $\Gamma$ is the Euler’s Gamma function, satisfying
$\Gamma(3/2)=\frac{\sqrt{\pi}}{2}$. Hence, assuming $\eta\geq 1$, we receive
for $I_{2}$

$\displaystyle 0\leq
I_{2}\leq\sqrt{2}\Bigl{(}\eta^{1/2}+\frac{\sqrt{\pi}}{2}\Bigr{)}\leq\sqrt{2}\eta^{3/2}\Bigl{(}1+\frac{\sqrt{\pi}}{2}\Bigr{)}.$

And the claim in (A.1) is shown with the constants
$c_{1}=\frac{2}{3\sqrt{\pi}}$ and
$c_{2}=\frac{2}{\sqrt{\pi}}\left(\frac{2}{3}+\sqrt{2}\Bigl{(}1+\frac{\sqrt{\pi}}{2}\Bigr{)}\right)$.
∎

##  Appendix B A technical result

In this section, we establish a technical result, stated in Lemma B.1, which
is crucial for the proof of bounds satisfied by the quasi Fermi potentials of
electrons and holes, see Lemma 5.4.

######  Lemma B.1.

Assume that the statistics function $\mathcal{F}_{\alpha}$ satisfies the
hypothesis (H1). Let us define the functions
${\mathcal{K}}_{\alpha}:(x,a)\in\mathbb{R}^{2}\mapsto{\mathcal{K}}_{\alpha}(x,a)\in\mathbb{R}$
and
${\mathcal{D}}_{\alpha}:(x,y,a,b)\in\mathbb{R}^{4}\mapsto{\mathcal{D}}_{\alpha}(x,y,a,b)\in\mathbb{R}$
by

$\displaystyle{\mathcal{K}}_{\alpha}(x,a)=\log({\mathcal{F}}_{\alpha}(x-a))-x,$
$\displaystyle{\mathcal{D}}_{\alpha}(x,y,a,b)=(x-y)\Bigl{[}B\Bigl{(}{\mathcal{K}}_{\alpha}(x,a)-{\mathcal{K}}_{\alpha}(y,b)\Bigr{)}{\mathcal{F}}_{\alpha}(x-a)-B\Bigl{(}{\mathcal{K}}_{\alpha}(y,b)-{\mathcal{K}}_{\alpha}(x,a)\Bigr{)}{\mathcal{F}}_{\alpha}(y-b)\Bigr{]}.$

Then, for all $\overline{\Phi},\overline{\Psi}\in\mathbb{R}$, the function
$\Upsilon_{\overline{\Phi},\overline{\Psi}}:\mathbb{R}\to\mathbb{R}$ defined
by

$\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)=\inf\left\\{{\mathcal{D}}_{\alpha}(x,y,a,b);\
-\overline{\Phi}\leq y\leq\overline{\Phi},-\overline{\Psi}\leq
a,b\leq\overline{\Psi}\right\\}$

verifies

$\displaystyle\lim_{x\to-\infty}\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)=+\infty\quad\text{and}\quad\lim_{x\to+\infty}\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)=+\infty.$

###### Proof.

Let us first remark that the function ${\mathcal{K}}_{\alpha}$ is non-
increasing with respect to its both variables $x$ and $a$ and that the
Bernoulli function $B$ is also non-increasing on $\mathbb{R}$. We assume that
$\overline{\Phi}$ and $\overline{\Psi}$ are given. The regularity of the
functions ensure that there exist positive constants
$\lambda,\underline{\mu},\overline{\mu}$ such that, for
$y\in[-\overline{\Phi},\overline{\Phi}]$ and
$a,b\in[-\overline{\Psi},\overline{\Psi}]$, we have

$\displaystyle-\lambda\leq{\mathcal{K}}_{\alpha}(y,b)\leq\lambda~{}\quad\text{and}\quad~{}\underline{\mu}\leq{\mathcal{F}}_{\alpha}(y-b)\leq\overline{\mu},~{}\quad\text{for
all}\;-\overline{\Phi}\leq y\leq\overline{\Phi},\ -\overline{\Psi}\leq
b\leq\overline{\Psi}.$

This implies the following inequalities, for $x\in\mathbb{R}$,
$y\in[-\overline{\Phi},\overline{\Phi}]$ and
$a,b\in[-\overline{\Psi},\overline{\Psi}]$,

$\displaystyle\hphantom{-}B\Bigl{(}{\mathcal{K}}_{\alpha}(x,-\overline{\Psi})+\lambda\Bigr{)}\leq$
$\displaystyle\\!\\!\\!\hphantom{-}B\Bigl{(}{\mathcal{K}}_{\alpha}(x,a)-{\mathcal{K}}_{\alpha}(y,b)\Bigr{)}\\!\\!\\!$
$\displaystyle\leq\hphantom{-}B\Bigl{(}{\mathcal{K}}_{\alpha}(x,\overline{\Psi})-\lambda\Bigr{)},$
$\displaystyle~{}-B\Bigl{(}-\lambda-{\mathcal{K}}_{\alpha}(x,-\overline{\Psi})\Bigr{)}\leq$
$\displaystyle\\!\\!\\!-B\Bigl{(}{\mathcal{K}}_{\alpha}(y,b)-{\mathcal{K}}_{\alpha}(x,a)\Bigr{)}\\!\\!\\!$
$\displaystyle\leq-B\Bigl{(}\lambda-{\mathcal{K}}_{\alpha}(x,\overline{\Psi})\Bigr{)},$

yielding

$B\Bigl{(}{\mathcal{K}}_{\alpha}(x,-\overline{\Psi})+\lambda\Bigr{)}{\mathcal{F}}_{\alpha}(x-\overline{\Psi})-B\Bigl{(}-\lambda-{\mathcal{K}}_{\alpha}(x,-\overline{\Psi})\Bigr{)}{\overline{\mu}}~{}\leq~{}\frac{{\mathcal{D}}_{\alpha}(x,y,a,b)}{x-y}~{}\leq\\\
~{}B\Bigl{(}{\mathcal{K}}_{\alpha}(x,\overline{\Psi})-\lambda\Bigr{)}{\mathcal{F}}_{\alpha}(x+\overline{\Psi})-B\Bigl{(}\lambda-{\mathcal{K}}_{\alpha}(x,\overline{\Psi})\Bigr{)}{\underline{\mu}}.$
(B.1)

Let us first consider that $x\leq-\overline{\Phi}$. Then, we deduce from
(B.1), that

$\frac{{\mathcal{D}}_{\alpha}(x,y,a,b)}{x-y}\leq
B\Bigl{(}{\mathcal{K}}_{\alpha}(-\overline{\Phi},\overline{\Psi})-\lambda\Bigr{)}{\mathcal{F}}_{\alpha}(x+\overline{\Psi})-B\Bigl{(}\lambda-{\mathcal{K}}_{\alpha}(-\overline{\Phi},\overline{\Psi})\Bigr{)}{\underline{\mu}}$

But, due to (H1),
$\displaystyle\lim_{x\to-\infty}{\mathcal{F}}_{\alpha}(x+\overline{\Psi})=0$,
which implies that, for $-x$ large enough, the right-hand-side of the last
inequality is negative. Therefore, for such an $x$ with $x\leq y$ we have

${\mathcal{D}}_{\alpha}(x,y,a,b)\geq(x-y)\Bigl{[}B\Bigl{(}{\mathcal{K}}_{\alpha}(-\overline{\Phi},\overline{\Psi})-\lambda\Bigr{)}{\mathcal{F}}_{\alpha}(x+\overline{\Psi})-B\Bigl{(}\lambda-{\mathcal{K}}_{\alpha}(-\overline{\Phi},\overline{\Psi})\Bigr{)}{\underline{\mu}}\Bigr{]}$

and, taking the infimum in $y\in[-\overline{\Phi},\overline{\Phi}]$, we obtain

$\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)\geq(x+\overline{\Phi})\Bigl{[}B\Bigl{(}{\mathcal{K}}_{\alpha}(-\overline{\Phi},\overline{\Psi})-\lambda\Bigr{)}{\mathcal{F}}_{\alpha}(x+\overline{\Psi})-B\Bigl{(}\lambda-{\mathcal{K}}_{\alpha}(-\overline{\Phi},\overline{\Psi})\Bigr{)}{\underline{\mu}}\Bigr{]}.$

As the first product in the right-hand-side tends to 0, while the second one
tends to $+\infty$, we deduce that

$\lim_{x\to-\infty}\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)=+\infty.$

We may now consider that $x\geq\overline{\Phi}$. From (B.1), we deduce that

$\frac{{\mathcal{D}}_{\alpha}(x,y,a,b)}{x-y}\geq\Bigl{(}B\bigl{(}{\mathcal{K}}_{\alpha}(\overline{\Phi},-\overline{\Psi})+\lambda\bigr{)}{\mathcal{F}}_{\alpha}(x-\overline{\Psi})-B\bigl{(}-\lambda-{\mathcal{K}}_{\alpha}(\overline{\Phi},-\overline{\Psi})\bigr{)}{\overline{\mu}}\Bigr{)}.$

For $x$ sufficiently large, the right-hand-side of the last inequality is
positive and

$\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)\geq(x-\overline{\Phi})\Bigl{(}B\bigl{(}{\mathcal{K}}_{\alpha}(\overline{\Phi},-\overline{\Psi})+\lambda\bigr{)}{\mathcal{F}}_{\alpha}(x-\overline{\Psi})-B\bigl{(}-\lambda-{\mathcal{K}}_{\alpha}(\overline{\Phi},-\overline{\Psi})\bigr{)}{\overline{\mu}}\Bigr{)}.$

Therefore, we get

$\lim_{x\to+\infty}\Upsilon_{\overline{\Phi},\overline{\Psi}}(x)=+\infty.$

∎

## References

  * [1] D. Abdel, P. Farrell, and J. Fuhrmann.  Assessing the quality of the excess chemical potential flux scheme for degenerate semiconductor device simulation.  Optical and Quantum Electronics, 53(163), 2021.
  * [2] D. Abdel, P. Farrell, and J. Fuhrmann.  ChargeTransport.jl: Simulating charge transport in semiconductors.  https://github.com/PatricioFarrell/ChargeTransport.jl, 2022.
  * [3] D. Abdel, P. Vágner, J. Fuhrmann, and P. Farrell.  Modelling charge transport in perovskite solar cells: Potential-based and limiting ion depletion.  Electrochimica Acta, 390:138696, 2021.
  * [4] G. Albinus, H. Gajewski, and R. Hünlich.  Thermodynamic design of energy models of semiconductor devices.  Nonlinearity, 15(2):367–383, 2002.
  * [5] M. Bessemoulin-Chatard.  A finite volume scheme for convection–diffusion equations with nonlinear diffusion derived from the scharfetter–gummel scheme.  Numerische Mathematik, 121(4):637–670, 2012.
  * [6] M. Bessemoulin-Chatard and C. Chainais-Hillairet.  Exponential decay of a finite volume scheme to the thermal equilibrium for drift–diffusion systems.  Journal of Numerical Mathematics, 25(3):147–168, 2017.
  * [7] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and M.-H. Vignal.  Study of a finite volume scheme for the drift-diffusion system. Asymptotic behavior in the quasi-neutral limit.  SIAM J. Numer. Anal., 52(4):1666–1691, 2014.
  * [8] P. Biler and J. Dolbeault.  Long time behavior of solutions to nernst-planck and debye-hückel drift-diffusion systems.  Annales Henri Poincaré, 1:461–472, 2000.
  * [9] F. Brezzi, L. Marini, S. Micheletti, P. Pietra, R. Sacco, and S. Wang.  Discretization of semiconductor device problems (i).  Handbook of numerical analysis, 13:317–441, 2005.
  * [10] P. Calado, A. Telford, D. Bryant, X. Li, J. Nelson, B. O’Regan, and P. R. F. Barnes.  Evidence for ion migration in hybrid perovskite solar cells with minimal hysteresis.  Nature Communications, 7, 2016.
  * [11] C. Cancès, C. Chainais-Hillairet, J. Fuhrmann, and B. Gaudeul.  A numerical-analysis-focused comparison of several finite volume schemes for a unipolar degenerate drift-diffusion model.  IMA Journal of Numerical Analysis, 41(1):271–314, 07 2020.
  * [12] C. Chainais-Hillairet and M. Herda.  Large-time behaviour of a family of finite volume schemes for boundary-driven convection–diffusion equations.  IMA Journal of Numerical Analysis, 40(4):2473–2504, 2019.
  * [13] N. E. Courtier.  Modelling ion migration and charge carrier transport in planar perovskite solar cells.  PhD thesis, University of Southampton, 2019.
  * [14] N. E. Courtier, J. M. Cave, A. B. Walker, G. Richardson, and J. M. Foster.  Ionmonger: a free and fast planar perovskite solar cell simulator with coupled ion vacancy and charge carrier dynamics.  Journal of Computational Electronics, 18:1435–1449, 2019.
  * [15] N. E. Courtier, G. Richardson, and J. M. Foster.  A fast and robust numerical scheme for solving models of charge carrier transport and ion vacancy motion in perovskite solar cells.  Applied Mathematical Modelling, 2018.
  * [16] L. C. Evans.  Partial Differential Equations: Second Edition, volume 19 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, R.I., 2010.
  * [17] R. Eymard, T. Gallouët, and R. Herbin.  Finite volume methods.  In Handbook of numerical analysis, Vol. VII, pages 713–1020. North-Holland, Amsterdam, 2000.
  * [18] P. Farrell, D. H. Doan, M. Kantner, J. Fuhrmann, T. Koprucki, and N. Rotundo.  Drift-diffusion models.  In Optoelectronic Device Modeling and Simulation: Fundamentals, Materials, Nanostructures, LEDs, and Amplifiers, pages 733–771. CRC Press Taylor & Francis Group, 2017.
  * [19] P. Farrell, M. Patriarca, J. Fuhrmann, and T. Koprucki.  Comparison of thermodynamically consistent charge carrier flux discretizations for fermi–dirac and gauss–fermi statistics.  Optical and Quantum Electronics, 50(2):1–10, 2018.
  * [20] H. Gajewski.  On existence, uniqueness and asymptotic behavior of solutions of the basic equations for carrier transport in semiconductors.  Z. Angew. Math. Mech., 65:101–108, 1985.
  * [21] H. Gajewski and K. Gröger.  On the basic equations for carrier transport in semiconductors.  J. Math. Anal. Appl., 113:12–35, 1986.
  * [22] H. Gajewski and K. Gröger.  Semiconductor equations for variable mobilities based on boltzmann statistics or fermi-dirac statistics.  Mathematische Nachrichten, 140(1):7–36, 1989.
  * [23] B. Gaudeul and J. Fuhrmann.  Entropy and convergence analysis for two finite volume schemes for a Nernst-Planck-Poisson system with ion volume constraints.  Numerische Mathematik, 151:99–149, 2022.
  * [24] A. Glitzky.  Uniform exponential decay of the free energy for voronoi finite volume discretized reaction-diffusion systems.  Mathematische Nachrichten, 284(17-18):2159–2174, 2011.
  * [25] A. Jüngel.  Numerical approximation of a drift-diffusion model for semiconductors with nonlinear diffusion.  ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 75(10):783–799, 1995.
  * [26] A. Jüngel.  Entropy methods for diffusive partial differential equations, volume 804.  Springer, 2016.
  * [27] A. Jüngel and Y.-J. Peng.  A hierarchy of hydrodynamic models for plasmas. Quasi-neutral limits in the drift-diffusion equations.  Asymptotic Anal., 28(1):49–73, 2001.
  * [28] M. Kantner.  Generalized Scharfetter–Gummel schemes for electro-thermal transport in degenerate semiconductors using the Kelvin formula for the Seebeck coefficient.  Journal of Computational Physics, 402:109091, 2020.
  * [29] M. Kantner and T. Koprucki.  Non-isothermal Scharfetter–Gummel Scheme for Electro-Thermal Transport Simulation in Degenerate Semiconductors.  In Finite Volumes for Complex Applications IX - Methods, Theoretical Aspects, Examples, pages 173–182. Springer International Publishing, 2020.
  * [30] J. Y. Kim, J.-W. Lee, H. S. Jung, H. Shin, and N.-G. Park.  High-efficiency perovskite solar cells.  Chemical Reviews, 120(15):7867–7918, 2020.
  * [31] C. Liu, C. Wang, S. Wise, X. Yue, and S. Zhou.  A positivity-preserving, energy stable and convergent numerical scheme for the poisson-nernst-planck system.  Mathematics of Computation, 90(331):2071–2106, 2021.
  * [32] P. A. Markowich.  The stationary semiconductor device equations.  Springer Science & Business Media, 1985.
  * [33] P. A. Markowich, C. A. Ringhofer, and C. Schmeiser.  Semiconductor equations.  Springer Science & Business Media, 2012.
  * [34] J. Moatti.  A structure preserving hybrid finite volume scheme for semi-conductor models with magnetic field on general meshes.  arXiv preprint arXiv:2207.02567, 2022.
  * [35] M. S. Mock.  Analysis of mathematical models of semiconductor devices, volume 3.  Boole Press, 1983.
  * [36] D. Scharfetter and H. Gummel.  Large-signal analysis of a silicon read diode oscillator.  IEEE Transactions on electron devices, 16(1):64–77, 1969.
  * [37] S. M. Sze and K. K. Ng.  Physics of Semiconductor Devices.  Wiley, 2006.
  * [38] N. Tessler and Y. Vaynzof.  Insights from device modeling of perovskite solar cells.  ACS Energy Letters, 5(4):1260–1270, 04 2020.
  * [39] W. Van Roosbroeck.  Theory of the flow of electrons and holes in germanium and other semiconductors.  The Bell System Technical Journal, 29(4):560–607, 1950.
  * [40] Z. Yu and R. Dutton.  SEDAN III – A one-dimensional device simulator.  www-tcad.stanford.edu/tcad/programs/sedan3.html, 1988.