which yields (1.50).

<figure><img src="image_15.png" /></figure>

We finally state the main result of Section 1.3.3, about the existence of (positive) solutions to the nonlinear scheme (1.45). The proof of this result is the subject of the next subsection.

**Theorem 1** (Existence of positive solutions). Let $u^{in} \in L^2(\Omega)$ be a non-negative function such that $\int_{\Omega} u^{in} = M > 0$. There exists at least one positive solution $(u_D^n \in V_D)_{n \ge 1}$ to the nonlinear scheme (1.45). Moreover, there exists $\varepsilon > 0$, depending on $\Lambda, \phi, u^{in}, M, \Omega, d, \Delta t$, and $\mathcal{D}$ such that

$$ \forall n \ge 1, \quad u_K^n \ge \varepsilon \quad \forall K \in \mathcal{M} \quad \text{and} \quad u_\sigma^n \ge \varepsilon \quad \forall \sigma \in \mathcal{E}. \qquad (1.51) $$

The uniform-in-time positivity result (1.51) on discrete solutions is the equivalent in the HFV context of [56, Lemma 3.7] and [52, Lemma 3.5] obtained, respectively, in the VAG and DDFV contexts.

### Existence of discrete solutions

The existence of discrete solutions to the nonlinear scheme (1.45) is proved in two steps. First, we introduce a regularised scheme, for which we prove the existence of solutions by a fixed-point argument, inspired from the proof of existence in [25]. Then, we prove that sequences of regularised solutions satisfy uniform a priori bounds, which allows us to pass to the limit in the regularisation parameter. Notice that our proof of existence uses the same estimates, but follows a quite different path than the ones in the VAG [56] and DDFV [52] contexts, in which the proof is based on the topological degree, together with a monotonicity argument. Henceforth, we reason in the $w_D$ variable, and we recall that $u_D = u_D^\infty \times \exp(w_D)$ according to (1.47). The advantage of doing so is that we can seek for solutions $w_D$ in the whole space $V_D$, with bijective correspondence with solutions $u_D$ that are automatically positive. Recalling the definition (1.49) of the discrete relative entropy and dissipation, and using (1.44) combined with (1.40), we let, for all $w_D \in V_D$,

$$ \mathbb{E}(w_M) = \sum_{K \in \mathcal{M}} |K| u_K^\infty \Phi_1(e^{w_K}), \quad \mathbb{D}(w_D) = \sum_{K \in \mathcal{M}} r_K (u_K^\infty \times \exp(w_K)) \delta_K w_K \cdot A_K \delta_K w_K, \quad (1.52) $$

in such a way that $\mathbb{E}^n = \mathbb{E}(w_M(u_M^n))$ and $\mathbb{D}^n = \mathbb{D}(w_D(u_D^n))$ for all $n \ge 1$. Using the fact that $\Phi_1(0) = 1$, we extend the definition of $\mathbb{E}(w_M)$ to the case where some $w_K$'s are equal to $-\infty$.

Before proceeding with the proof of Theorem 1, we state two preliminary lemmas. The first one, that can be found, e.g., in [119, Section 9.1], is a corollary of Brouwer's fixed-point theorem. This result is instrumental to show the existence of solutions to the regularised scheme.

**Lemma 1.** Let $N \in \mathbb{N}^*$, and let $P: \mathbb{R}^N \to \mathbb{R}^N$ be a continuous vector field. Assume that there is $r > 0$ such that

$$ P(x) \cdot x \ge 0 \quad \text{if } |x| = r. $$

Then, there exists a point $x_0 \in \mathbb{R}^N$ such that $P(x_0) = 0$ and $|x_0| \le r$.

The second lemma, whose proof is postponed until Appendix 1.C.1, establishes sufficient boundedness conditions on the discrete mass and (relative) dissipation so that a priori bounds hold for vectors of discrete unknowns. This result is instrumental to show that sequences of regularised solutions satisfy (uniform) a priori bounds.
**Lemma 2.** Let $\underline{w}_D \in \underline{V}_D$, and assume that there exist $C_\# > 0$, and $M_\# \ge M_b > 0$ such that

$$M_b \le \sum_{K \in \mathcal{M}} |K| u_K^\infty e^{w_K} \le M_\# \quad \text{and} \quad \mathbb{D}(\underline{w}_D) \le C_\#. \qquad (1.53)$$

Then, there exists $C > 0$, depending on $\Lambda$, $u_b^\infty$, $u_\#^\infty$, $M_b$, $M_\#$, $C_\#$, $\Omega$, $d$, and $\mathcal{D}$ such that

$$|w_K| \le C \quad \forall K \in \mathcal{M} \quad \text{and} \quad |w_\sigma| \le C \quad \forall \sigma \in \mathcal{E}.$$

We can now proceed with the proof of Theorem 1. Let us first define the following inner product and corresponding norm on the space $\underline{V}_D$: for all $\underline{z}_D, \underline{v}_D \in \underline{V}_D$,

$$\langle \underline{z}_D, \underline{v}_D \rangle = \sum_{K \in \mathcal{M}} z_K v_K + \sum_{\sigma \in \mathcal{E}} z_\sigma v_\sigma \quad \text{and} \quad \| \underline{v}_D \| = \sqrt{\langle \underline{v}_D, \underline{v}_D \rangle}.$$

Letting $N = |\mathcal{M}| + |\mathcal{E}|$, and identifying $\underline{V}_D$ to $\mathbb{R}^N$, the inner product $\langle \cdot, \cdot \rangle$ is nothing but the standard inner product on $\mathbb{R}^N$. For all $K \in \mathcal{M}$, and all $\sigma \in \mathcal{E}_K$, we let, for $\underline{u}_K \in \underline{V}_K$ positive,

$$\mathcal{F}_{K,\sigma}^{\text{nl}}(\underline{u}_K) = r_K(\underline{u}_K) \sum_{\sigma' \in \mathcal{E}_K} A_K^{\sigma\sigma'} \left( \log\left(\frac{u_K}{u_K^\infty}\right) - \log\left(\frac{u_{\sigma'}}{u_{\sigma'}^\infty}\right) \right), \qquad (1.54)$$

where the $A_K^{\sigma\sigma'}$ are defined by (1.10). Combining (1.39) and (1.47) with (1.11) and (1.9), there holds that $T_K(\underline{u}_K, \underline{w}_K, \underline{v}_K) = \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{K,\sigma}^{\text{nl}}(\underline{u}_K)(v_K - v_\sigma)$ for all $\underline{v}_K \in \underline{V}_K$. In what follows, we let $n \in \mathbb{N}^*$ and $u_M^{n-1} \ge 0$ be given. We assume that $M^{n-1} = \int_{\Omega} u_M^{n-1} > 0$ and that $u_M^{n-1} > 0$ if $n > 1$. We also assume that $\mathbb{E}(w_M^{n-1}) > 0$. If $\mathbb{E}(w_M^{n-1}) = 0$ (which is equivalent to $u_M^{n-1} = u_M^\infty$), then necessarily, by (1.50), $\underline{u}_D^n = \underline{u}_D^\infty$ uniquely solves (1.45a). Letting, for any $\underline{u}_D \in \underline{V}_D$ positive, $\underline{G}_D^n(\underline{u}_D)$ be the element of $\underline{V}_D$ such that

$$\mathcal{G}_K^n(\underline{u}_D) = |K| \frac{u_K - u_K^{n-1}}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{K,\sigma}^{\text{nl}}(\underline{u}_K) \quad \forall K \in \mathcal{M}, \qquad (1.55a)$$

$$\mathcal{G}_\sigma^n(\underline{u}_D) = -\left(\mathcal{F}_{K,\sigma}^{\text{nl}}(\underline{u}_K) + \mathcal{F}_{L,\sigma}^{\text{nl}}(\underline{u}_L)\right) \quad \forall \sigma = K \mid L \in \mathcal{E}_{\text{int}}, \qquad (1.55b)$$

$$\mathcal{G}_\sigma^n(\underline{u}_D) = -\mathcal{F}_{K,\sigma}^{\text{nl}}(\underline{u}_K) \quad \forall \sigma \in \mathcal{E}_{\text{ext}} \text{ with } \mathcal{M}_\sigma = \{K\}, \qquad (1.55c)$$

we infer that, for all $\underline{v}_D \in \underline{V}_D$,

$$\frac{1}{\Delta t} (u_M - u_M^{n-1}, v_M)_\Omega + T_D(\underline{u}_D, \underline{w}_D(\underline{u}_D), \underline{v}_D) = \langle \underline{G}_D^n(\underline{u}_D), \underline{v}_D \rangle. \qquad (1.56)$$

Hence, a positive vector $\underline{u}_D^n \in \underline{V}_D$ is a solution to the nonlinear equation (1.45a) if and only if $\underline{G}_D^n(\underline{u}_D^n) = \underline{0}_D$. With this observation in hand, we now detail the two steps of the proof.

**Step 1:** Using the relation $\underline{u}_D = \underline{u}_D^\infty \times \exp(\underline{w}_D)$, we define the vector field $\underline{P}_D^{n,\mu} : \underline{V}_D \to \underline{V}_D$ such that, for all $\underline{w}_D \in \underline{V}_D$,

$$\underline{P}_D^{n,\mu}(\underline{w}_D) = \underline{G}_D^n(\underline{u}_D^\infty \times \exp(\underline{w}_D)) + \mu \underline{w}_D, \qquad (1.57)$$

with $\underline{G}_D^n$ defined by (1.55) and $\mu \ge 0$. Notice that, unlike $\underline{G}_D^n$, the vector field $\underline{P}_D^{n,\mu}$ is continuous on the whole space $\underline{V}_D$ for any $\mu \ge 0$. If $\underline{w}_D^n \in \underline{V}_D$ satisfies $\underline{P}_D^{n,0}(\underline{w}_D^n) = \underline{0}_D$, then letting $\underline{u}_D^n =$
$\underline{u}_D^\infty \times \exp(\underline{w}_D^n)$, we have $\underline{G}_D^n(\underline{u}_D^n) = \underline{0}_D$, therefore $\underline{u}_D^n$ is a (positive) solution to (1.45a). For $\mu > 0$, the problem of finding $\underline{w}_D^{n,\mu} \in \underline{V}_D$ such that $\underline{P}_D^{n,\mu}(\underline{w}_D^{n,\mu}) = \underline{0}_D$ can thus be seen as a regularisation of the original problem. By (1.57) and (1.56), for all $\underline{w}_D \in \underline{V}_D$, we have

$$ \langle \underline{P}_D^{n,\mu}(\underline{w}_D), \underline{w}_D \rangle = \sum_{K \in \mathcal{M}} \frac{|K|}{\Delta t} (u_K^\infty e^{w_K} - u_K^{n-1}) w_K + \sum_{K \in \mathcal{M}} r_K (\underline{u}_K^\infty \times \exp(\underline{w}_K)) \delta_K \underline{w}_K \cdot \mathbb{A}_K \delta_K \underline{w}_K + \mu \|\underline{w}_D\|^2. \quad (1.58) $$

By (1.52), we recognise in the second term of the right-hand side the quantity $\mathbb{D}(\underline{w}_D) \ge 0$. As far as the first term is concerned, for $n > 1$, by positivity of the $(u_K^{n-1})_{K \in \mathcal{M}}$, there exist real numbers $(w_K^{n-1})_{K \in \mathcal{M}}$ such that $u_K^{n-1} = u_K^\infty e^{w_K^{n-1}}$ for all $K \in \mathcal{M}$, and since $\Phi_1$ is convex,

$$ \begin{aligned} \sum_{K \in \mathcal{M}} \frac{|K|}{\Delta t} (u_K^\infty e^{w_K} - u_K^{n-1}) w_K &= \sum_{K \in \mathcal{M}} \frac{|K|}{\Delta t} u_K^\infty w_K (e^{w_K} - e^{w_K^{n-1}}) \\ &\ge \sum_{K \in \mathcal{M}} \frac{|K|}{\Delta t} u_K^\infty (\Phi_1(e^{w_K}) - \Phi_1(e^{w_K^{n-1}})) = \frac{\mathbb{E}(w_\mathcal{M}) - \mathbb{E}(w_\mathcal{M}^{n-1})}{\Delta t}, \end{aligned} \quad (1.59) $$

where we have used the definition (1.52) of $\mathbb{E}(w_\mathcal{M})$. For $n = 1$, now, it may happen that $u_K^0$ be zero for some $K \in \mathcal{M}$, and then $w_K^0$ such that $u_K^0 = u_K^\infty e^{w_K^0}$ cannot be defined. However, letting in that case $w_K^0 = -\infty$, the inequality above still holds since $\Phi_1(0) = 1$ and $\Phi_1(e^s) - 1 \le s e^s$ for all $s \in \mathbb{R}$. By non-negativity of $\mathbb{E}(w_\mathcal{M})$, we finally infer from (1.58) and (1.59) that

$$ \langle \underline{P}_D^{n,\mu}(\underline{w}_D), \underline{w}_D \rangle \ge \mu \|\underline{w}_D\|^2 - \frac{\mathbb{E}(w_\mathcal{M}^{n-1})}{\Delta t}, $$

so that, for $\mu > 0$, there holds $\langle \underline{P}_D^{n,\mu}(\underline{w}_D), \underline{w}_D \rangle \ge 0$ if $\|\underline{w}_D\| = \sqrt{\frac{\mathbb{E}(w_\mathcal{M}^{n-1})}{\mu \Delta t}} > 0$. By Lemma 1, we then conclude about the existence of solutions to the regularised scheme. There exists $\underline{w}_D^{n,\mu} \in \underline{V}_D$ such that

$$ \underline{P}_D^{n,\mu}(\underline{w}_D^{n,\mu}) = \underline{0}_D \quad \text{and} \quad \|\underline{w}_D^{n,\mu}\| \le \sqrt{\frac{\mathbb{E}(w_\mathcal{M}^{n-1})}{\mu \Delta t}}. \quad (1.60) $$

**Step 2:** Since $\langle \underline{P}_D^{n,\mu}(\underline{w}_D^{n,\mu}), \underline{w}_D^{n,\mu} \rangle = 0$, by (1.58) and (1.59), we have

$$ \frac{\mathbb{E}(w_\mathcal{M}^{n,\mu})}{\Delta t} + \mathbb{D}(\underline{w}_D^{n,\mu}) + \mu \|\underline{w}_D^{n,\mu}\|^2 \le \frac{\mathbb{E}(w_\mathcal{M}^{n-1})}{\Delta t}. $$

The three terms on the left-hand side being non-negative, we infer that

$$ \mathbb{D}(\underline{w}_D^{n,\mu}) \le C_\# \quad (1.61) $$

with $C_\# = \frac{\mathbb{E}(w_\mathcal{M}^{n-1})}{\Delta t} > 0$. Moreover, since $\langle \underline{P}_D^{n,\mu}(\underline{w}_D^{n,\mu}), \underline{1}_D \rangle = 0$, by (1.57) and (1.56), we have

$$ \sum_{K \in \mathcal{M}} |K| u_K^\infty e^{w_K^{n,\mu}} - M^{n-1} = -\mu \Delta t \langle \underline{w}_D^{n,\mu}, \underline{1}_D \rangle. $$
Applying a Cauchy–Schwarz inequality, and recalling the bound (1.60), we obtain

$$ \left| \sum_{K \in \mathcal{M}} |K| u_K^{\infty} e^{w_K^{n, \mu}} - M^{n-1} \right| \le \mu \Delta t \|\underline{1}_D\| \|\underline{w}_D^{n, \mu}\| \le \sqrt{\mu} \sqrt{N \Delta t \mathbb{E}(w_M^{n-1})}, $$

so that, letting $\mu_0 = \frac{(M^{n-1})^2}{4N\Delta t \mathbb{E}(w_M^{n-1})} > 0$, the following holds for all $0 < \mu \le \mu_0$:

$$ 0 < \frac{M^{n-1}}{2} = M_b \le \sum_{K \in \mathcal{M}} |K| u_K^{\infty} e^{w_K^{n, \mu}} \le M_{\#} = \frac{3M^{n-1}}{2}. \quad (1.62) $$

By (1.62) and (1.61), we infer that $\underline{w}_D^{n, \mu}$ satisfies (1.53) for $\mu$ sufficiently small with constants that are uniform in $\mu$, so that by Lemma 2 the family $(\underline{w}_D^{n, \mu})_{0<\mu\le\mu_0}$ is bounded uniformly in $\mu$. As a consequence, by compactness, there is $\underline{w}_D^n \in \underline{V}_D$ such that, up to extraction, $\underline{w}_D^{n, \mu}$ converges towards $\underline{w}_D^n$ when $\mu$ tends to zero. Since $\underline{P}_D^{n, \mu}$ converges to $\underline{P}_D^{n, 0}$ as $\mu$ tends to zero, we finally infer that $\underline{P}_D^{n, 0}(\underline{w}_D^n) = \underline{0}_D$ (also, $\sum_{K \in \mathcal{M}} |K| u_K^{\infty} e^{w_K^n} = M^{n-1}$).

**Conclusion:** Letting $\underline{u}_D^n = \underline{u}_D^\infty \times \exp(\underline{w}_D^n)$, we have $\underline{G}_D^n(\underline{u}_D^n) = \underline{0}_D$, therefore $\underline{u}_D^n$ is a (positive) solution to (1.45a). By Propositions 4 and 5, and since $\underline{D}(\underline{w}_D^n) = \underline{D}^n$, $\underline{E}^n$ is non-negative, and $\underline{E}^n$ is non-increasing in $n$ according to (1.50), we deduce that

$$ \sum_{K \in \mathcal{M}} |K| u_K^{\infty} e^{w_K^n} = M \quad \text{and} \quad \underline{D}(\underline{w}_D^n) \le \frac{\underline{E}^0}{\Delta t}. $$

By Lemma 2 (recall also that $u_b^\infty, u_\#^\infty$ only depend on $M, \phi$, and $\Omega$), there exists $C > 0$, depending on $\Lambda, \phi, u^{in}, M, \Omega, d, \Delta t$, and $\mathcal{D}$, but not on $n \in \mathbb{N}^*$, such that $-C\underline{1}_D \le \underline{w}_D^n \le C\underline{1}_D$. By (1.48), we finally infer that $\underline{u}_D^n \ge \varepsilon\underline{1}_D$, with $\varepsilon = u_b^\infty e^{-C} > 0$ still independent of $n \in \mathbb{N}^*$. This concludes the proof of Theorem 1.

## 1.4 Long-time behaviour

In this section, we analyse the long-time behaviour of the three HFV schemes we have introduced in Section 1.3, thereby proving the main results of this paper.

**Remark 8 (Linear schemes and nonhomogeneous data).** In order to stay consistent with Section 1.3, we here below state our asymptotic results of Theorems 2 and 3, which respectively concern the (linear) standard and exponential fitting schemes, for discrete problems that feature homogeneous data (i.e., $g^D = 0$ when $|\Gamma^D| > 0$, or $M = 0$ otherwise). Nonetheless, Theorems 2 and 3 remain valid in the general case of nonhomogeneous data (we refer to Remark 2 for the straightforward adaptation of the schemes to this situation). Indeed, the proofs of the latter results only hinge on the fact that the difference between the discrete transient and steady-state solutions belongs to the homogeneous space $\underline{V}_{D,0}$ or $\underline{V}_{D,0}^\omega$, which is always true. This remark does not apply, however, to the nonlinear scheme.

### 1.4.1 Asymptotic behaviour of the standard HFV scheme

We recall that $u$ is the solution to Problem (1.1), and that $u^\infty$ is the corresponding steady-state, solution to Problem (1.2), and we consider the following definition of the relative entropy
and dissipation:

$$ \mathbb{E}(t) = \frac{1}{2} \|u(t) - u^{\infty}\|_{L^2(\Omega)}^2, \quad \mathbb{D}(t) = \int_{\Omega} (\Lambda \nabla(u(t) - u^{\infty}) - (u(t) - u^{\infty})V^{\phi}) \cdot \nabla(u(t) - u^{\infty}). $$

It can be easily verified that the following entropy/dissipation relation holds at the continuous level:

$$ \frac{d}{dt} \mathbb{E}(t) + \mathbb{D}(t) = 0. $$

It is assumed that $V^{\phi}$ is such that $\mathbb{D}(t) \ge C \|\nabla(u(t) - u^{\infty})\|_{L^2(\Omega; \mathbb{R}^d)}^2$ for some $C > 0$, so that $\mathbb{D}(t)$ indeed defines a dissipation.

At the discrete level, recalling that $(\underline{u}_D^n \in \underline{V}_{D,0})_{n \ge 1}$ is the solution to Problem (1.26), and that $\underline{u}_D^\infty \in \underline{V}_{D,0}$ is the corresponding steady-state, solution to Problem (1.20), we consider the following equivalents of the relative entropy and dissipation: for all $n \in \mathbb{N}^*$,

$$ \mathbb{E}^n = \frac{1}{2} \|u_M^n - u_M^\infty\|_{L^2(\Omega)}^2, \quad \mathbb{D}^n = a_D(\underline{u}_D^n - \underline{u}_D^\infty, \underline{u}_D^n - \underline{u}_D^\infty), $$

where the discrete bilinear form $a_D$ is defined by (1.15) and (1.19). The definition of the relative entropy is seamlessly extended to the case $n = 0$. Our main result on the standard HFV scheme is the following.

**Theorem 2 (Asymptotic stability).** Assume that the advection field $V^{\phi}$ satisfies the conditions (1.21) of Proposition 2, with constant $\beta$. Then, the following discrete entropy/dissipation relation holds true:

$$ \forall n \in \mathbb{N}, \quad \frac{\mathbb{E}^{n+1} - \mathbb{E}^n}{\Delta t} + \mathbb{D}^{n+1} \le 0. \qquad (1.63) $$

Furthermore, the discrete entropy decays exponentially fast in time: there is $v = \frac{2\kappa}{C_p^2} > 0$, where $\kappa$ is the constant of (1.22) (only depending on $\Lambda$, $\beta$, $\Omega$, $d$, $\Gamma^D$, and $\theta_D$), and $C_P$ is either equal to $C_{PW}$ if $|\Gamma^D| = 0$ or to $C_{P,\Gamma^D}$ otherwise (where $C_{PW}$, $C_{P,\Gamma^D}$ are the Poincaré constants of Proposition 7), such that

$$ \forall n \in \mathbb{N}, \quad \mathbb{E}^{n+1} \le (1 + \nu \Delta t)^{-1} \mathbb{E}^n. \qquad (1.64) $$

Consequently, the discrete solution converges exponentially fast in time towards its associated discrete steady-state: for all $n \in \mathbb{N}^*$,

$$ \|u_M^n - u_M^\infty\|_{L^2(\Omega)} \le (1 + \nu \Delta t)^{-\frac{n}{2}} \|u_M^0 - u_M^\infty\|_{L^2(\Omega)}. \qquad (1.65) $$

*Proof.* Let $n \in \mathbb{N}$. One has

$$ \mathbb{E}^{n+1} - \mathbb{E}^n = \sum_{K \in \mathcal{M}} \frac{|K|}{2} \left( (u_K^{n+1} - u_K^\infty)^2 - (u_K^n - u_K^\infty)^2 \right). $$

Since $x \mapsto x^2$ is convex, for all $x, y \in \mathbb{R}$, we have $y^2 - x^2 \le 2y(y - x)$, therefore

$$ \mathbb{E}^{n+1} - \mathbb{E}^n \le \sum_{K \in \mathcal{M}} |K|(u_K^{n+1} - u_K^n)(u_K^{n+1} - u_K^\infty) = (u_M^{n+1} - u_M^n, u_M^{n+1} - u_M^\infty)_\Omega. \qquad (1.66) $$
Now, testing (1.26a) with $\underline{v}_D^{n+1} = \underline{u}_D^{n+1} - \underline{u}_D^\infty \in \underline{V}_{D,0}$ yields

$$(u_M^{n+1} - u_M^n, v_M^{n+1})_{\Omega} = -\Delta t a_D (\underline{u}_D^{n+1}, \underline{v}_D^{n+1}) + \Delta t ((f, v_M^{n+1})_{\Omega} + (g^N, v_{\mathcal{E}}^{n+1})_{\Gamma^N}).$$

By definition (1.20) of the discrete steady-state $\underline{u}_D^\infty$, we also have

$$(f, v_M^{n+1})_{\Omega} + (g^N, v_{\mathcal{E}}^{n+1})_{\Gamma^N} = a_D(\underline{u}_D^\infty, \underline{v}_D^{n+1}),$$

whence, by bilinearity of $a_D$, we infer

$$(u_M^{n+1} - u_M^n, v_M^{n+1})_{\Omega} = -\Delta t a_D (\underline{u}_D^{n+1} - \underline{u}_D^\infty, \underline{v}_D^{n+1}) = -\Delta t \mathbb{D}^{n+1}.$$

Combined to (1.66), this proves the entropy/dissipation relation (1.63). Now, since the advection field $V^\phi$ satisfies (1.21), we can invoke (1.22) from Proposition 2 to infer that

$$\mathbb{D}^{n+1} \geq \kappa |\underline{u}_D^{n+1} - \underline{u}_D^\infty|_{1,D}^2,$$

where $\kappa > 0$ only depends on $\Lambda, \beta, \Omega, d, \Gamma^D$, and $\theta_D$. Combining this estimate with a discrete Poincaré inequality from Proposition 7 (applied to $\underline{u}_D^{n+1} - \underline{u}_D^\infty \in \underline{V}_{D,0}$), and with the definition of the discrete (relative) entropy, yields

$$\mathbb{D}^{n+1} \geq \frac{\kappa}{C_p^2} \|u_M^{n+1} - u_M^\infty\|_{L^2(\Omega)}^2 = \frac{2\kappa}{C_p^2} \mathbb{E}^{n+1},$$

where $C_p$ is either equal to $C_{PW}$ if $|\Gamma^D| = 0$ or to $C_{P,\Gamma^D}$ otherwise. This last inequality, combined with the entropy/dissipation relation (1.63), implies the entropy decay (1.64). The inequality (1.65) is then a straightforward consequence of the definition of $\mathbb{E}^n$. $\square$

The result of Theorem 2 does not use the fact that $V^\phi$ is related to the gradient of a potential, it thus extends to general advection fields.

## 1.4.2 Asymptotic behaviour of the exponential fitting scheme

We recall that $\rho = \frac{u}{\omega}$, with $\omega = e^{-\phi}$, is the solution to Problem (1.34), and that $\rho^\infty$ is the corresponding steady-state, solution to Problem (1.28). We consider the following $\omega$-weighted definitions of the relative entropy and dissipation:

$$\mathbb{E}_{\omega}(t) = \frac{1}{2} \|\rho(t) - \rho^{\infty}\|_{L_{\omega}^{2}(\Omega)}^{2}, \quad \mathbb{D}_{\omega}(t) = \int_{\Omega} \omega \Lambda \nabla(\rho(t) - \rho^{\infty}) \cdot \nabla(\rho(t) - \rho^{\infty}).$$

It can be easily verified that the following entropy/dissipation relation holds at the continuous level:

$$\frac{d}{dt} \mathbb{E}_{\omega}(t) + \mathbb{D}_{\omega}(t) = 0.$$

At the discrete level, let us recall that $(\underline{\rho}_D^n \in \underline{V}_{D,0}^\omega)_{n \geq 1}$ is the solution to Problem (1.35). We then set $\underline{u}_D^n = \underline{\omega}_D \times \underline{\rho}_D^n$ with $\underline{\omega}_D$ defined by (1.31), in such a way that $\underline{u}_D^n \in \underline{V}_{D,0}$. Similarly, $\underline{\rho}_D^\infty \in \underline{V}_{D,0}^\omega$ is the corresponding steady-state, solution to Problem (1.30), and we set $\underline{u}_D^\infty = \underline{\omega}_D \times \underline{\rho}_D^\infty \in \underline{V}_{D,0}$. We consider the following equivalents of the $\omega$-weighted (relative) entropy and
dissipation: for all $n \in \mathbb{N}^*$,

$$ \mathbb{E}_{\omega}^{n} = \frac{1}{2} \|\rho_{M}^{n} - \rho_{M}^{\infty}\|_{L_{\omega}^{2}(\Omega)}^{2}, \quad \mathbb{D}_{\omega}^{n} = a_{D}^{\omega}(\underline{\rho}_{D}^{n} - \underline{\rho}_{D}^{\infty}, \underline{\rho}_{D}^{n} - \underline{\rho}_{D}^{\infty}), $$

where the discrete bilinear form $a_D^\omega$ is defined (locally) by (1.29). The definition of the relative entropy is seamlessly extended to the case $n = 0$. Our main result on the exponential fitting HFV scheme is the following, whose proof is very similar to the one of Theorem 2.

**Theorem 3 (Asymptotic stability).** The following discrete entropy/dissipation relation holds true:

$$ \forall n \in \mathbb{N}, \quad \frac{\mathbb{E}_{\omega}^{n+1} - \mathbb{E}_{\omega}^{n}}{\Delta t} + \mathbb{D}_{\omega}^{n+1} \le 0. \qquad (1.67) $$

Furthermore, the discrete entropy decays exponentially fast in time: there is $v_\omega = \frac{2\omega_b\lambda_b\alpha_b}{\omega_\#C_P^2} > 0$, where $\lambda_b\alpha_b$ is the coercivity constant of (1.14) (only depending on $\Lambda$, $\Omega$, $d$, and $\theta_D$), (we recall that) the bounds $\omega_b, \omega_\#$ only depend on $\phi$ and $\Omega$, and $C_P$ is either equal to $2C_{PW}$ if $|\Gamma^D| = 0$ or to $C_{P,\Gamma^D}$ otherwise (where $C_{PW}, C_{P,\Gamma^D}$ are the Poincaré constants of Proposition 7), such that

$$ \forall n \in \mathbb{N}, \quad \mathbb{E}_{\omega}^{n+1} \le (1 + v_{\omega} \Delta t)^{-1} \mathbb{E}_{\omega}^{n}. \qquad (1.68) $$

Consequently, the discrete solution converges exponentially fast in time towards its associated discrete steady-state: for all $n \in \mathbb{N}^*$,

$$ \|u_M^n - u_M^\infty\|_{L^2(\Omega)} \le \sqrt{\frac{\omega_\#}{\omega_b}} (1 + v_\omega \Delta t)^{-\frac{n}{2}} \|u_M^0 - u_M^\infty\|_{L^2(\Omega)}. \qquad (1.69) $$

*Proof.* Let $n \in \mathbb{N}$. Reasoning as in the proof of Theorem 2, we infer

$$ \mathbb{E}_{\omega}^{n+1} - \mathbb{E}_{\omega}^{n} \le -\Delta t a_{D}^{\omega}(\underline{\rho}_{D}^{n+1} - \underline{\rho}_{D}^{\infty}, \underline{\rho}_{D}^{n+1} - \underline{\rho}_{D}^{\infty}) = -\Delta t \mathbb{D}_{\omega}^{n+1}, $$

which proves (1.67). By the coercivity estimate (1.33), we have that $\mathbb{D}_{\omega}^{n+1} \ge \omega_b\lambda_b\alpha_b|\underline{\rho}_{D}^{n+1} - \underline{\rho}_{D}^{\infty}|_{1,D}^2$ where $\alpha_b > 0$ from Section 1.2.3 only depends on $\Omega$, $d$, and $\theta_D$, and $\omega_b > 0$ only depends on $\phi$ and $\Omega$. Reasoning as in the proof of Proposition 3, and using a discrete Poincaré inequality from Proposition 7 (combined with [51, Lemma 5.2] in the case $|\Gamma^D| = 0$), we also infer that

$$ \|\rho_M^{n+1} - \rho_M^\infty\|_{L^2(\Omega)} \le \sqrt{\omega_\# C_P} |\rho_D^{n+1} - \rho_D^\infty|_{1,D}, $$

where $\omega_\# > 0$ only depends on $\phi$ and $\Omega$, and $C_P$ is either equal to $2C_{PW}$ if $|\Gamma^D| = 0$ or to $C_{P,\Gamma^D}$ otherwise. Thus, we finally get that

$$ \mathbb{D}_{\omega}^{n+1} \ge \frac{2\omega_b\lambda_b\alpha_b}{\omega_\# C_P^2} \mathbb{E}_{\omega}^{n+1}. $$

Combined to (1.67), this yields (1.68). Deriving (1.69) is then straightforward. $\square$

<figure><img src="image_16.png" /></figure>

### 1.4.3 Asymptotic behaviour of the nonlinear scheme

Recall that $u > 0$ is the solution to Problem (1.1) endowed with pure Neumann boundary conditions ($|\Gamma^D| = 0$), and data $f = 0, g^N = 0$, and $u^{in} \ge 0$ with $M = \int_{\Omega} u^{in} > 0$ (see Appendix 1.B
for the case of mixed Dirichlet-Neumann boundary conditions), and that $u^\infty > 0$, solution to Problem (1.2) with same data, is the thermal equilibrium $u_{th}^\infty$ given by (1.25) with $\bar{\rho} = \rho^M$. The analysis of the nonlinear scheme relies on the entropy/dissipation structure (1.36)-(1.37) introduced in Section 1.3.3. Notice that the relative dissipation (or relative Fisher information) of (1.37) can be equivalently rewritten

$$ \mathbb{D}(t) = \int_{\Omega} u(t) \Lambda \nabla \log \left( \frac{u(t)}{u^{\infty}} \right) \cdot \nabla \log \left( \frac{u(t)}{u^{\infty}} \right) = 4 \int_{\Omega} u^{\infty} \Lambda \nabla \sqrt{\frac{u(t)}{u^{\infty}}} \cdot \nabla \sqrt{\frac{u(t)}{u^{\infty}}}. $$

At the discrete level, recalling that $(\underline{u}_D^n \in \underline{V}_D)_{n \ge 1}$ is a (positive) solution to Problem (1.45), and that $\underline{u}_D^n \in \underline{V}_D$ is the corresponding steady-state, solution to Problem (1.46), that is equal to the HFV interpolate of $u_{th}^\infty$, we consider the discrete entropy $\mathbb{E}^n$ and dissipation $\mathbb{D}^n$ defined by (1.49), and we define a discrete counterpart of the relative dissipation written in root-form: for all $n \ge 1$,

$$ \hat{\mathbb{D}}^n = 4 \sum_{K \in \mathcal{M}} u_{K,b}^{\infty} \int_K \Lambda \nabla_K \underline{\xi}_K^n \cdot \nabla_K \underline{\xi}_K^n = 4 \sum_{K \in \mathcal{M}} u_{K,b}^{\infty} \delta_K \underline{\xi}_K^n \cdot \mathbb{A}_K \delta_K \underline{\xi}_K^n, \quad (1.70) $$

where, for all $K \in \mathcal{M}$, $u_{K,b}^\infty = \min(u_K^\infty, \min_{\sigma \in \mathcal{E}_K} u_\sigma^\infty)$ and the matrix $\mathbb{A}_K$ is defined by (1.10), and $\underline{\xi}_D^n$ is the element of $\underline{V}_D$ such that

$$ \xi_K^n = \sqrt{\frac{u_K^n}{u_K^\infty}} \quad \forall K \in \mathcal{M}, \quad \xi_\sigma^n = \sqrt{\frac{u_\sigma^n}{u_\sigma^\infty}} \quad \forall \sigma \in \mathcal{E}. \quad (1.71) $$

At the discrete level, and as opposed to the continuous level, the quantities $\mathbb{D}^n$ and $\hat{\mathbb{D}}^n$ are not equal, therefore we need to compare them. The definition of $u_{K,b}^\infty$ results from the following observation: according to the structures of $\mathbb{D}^n$ and $\hat{\mathbb{D}}^n$, locally, we expect to have to compare $u_{K,b}^\infty$ with $r_K(\underline{u}_K^n)$, which depends on $u_K^n$ and on the $(u_\sigma^n)_{\sigma \in \mathcal{E}_K}$.

**Proposition 6** (Fisher information). There is $C_F > 0$, only depending on $\Lambda$, $\Omega$, $d$, and $\theta_D$ such that

$$ \forall n \ge 1, \quad \hat{\mathbb{D}}^n \le C_F \mathbb{D}^n. \quad (1.72) $$

*Proof.* Let $n \in \mathbb{N}^*$, and $K \in \mathcal{M}$. By (1.96) from Lemma 4, we first have that

$$ 4u_{K,b}^{\infty} \delta_K \underline{\xi}_K^n \cdot \mathbb{A}_K \delta_K \underline{\xi}_K^n \le 4u_{K,b}^{\infty} \delta_K \underline{\xi}_K^n \cdot \mathbb{B}_K \delta_K \underline{\xi}_K^n, \quad (1.73) $$

where the matrix $\mathbb{B}_K$ is the diagonal matrix defined by (1.95). Since for all $(x,y) \in (\mathbb{R}_+^*)^2$, $\int_x^y \frac{dz}{\sqrt{z}} = 2(\sqrt{y} - \sqrt{x})$, the Cauchy-Schwarz inequality yields

$$ 4(\sqrt{y} - \sqrt{x})^2 \le (y - x)(\log(y) - \log(x)). $$

By the property (1.42d) of the function $m$, we then get that

$$ \forall (x, y) \in (\mathbb{R}_+^*)^2, \quad 4(\sqrt{y} - \sqrt{x})^2 \le m(x, y)(\log(y) - \log(x))^2. \quad (1.74) $$
Since $\mathbb{B}_K$ is diagonal, the combination of (1.73), (1.71), and (1.74), yields

$$ 
\begin{aligned}
4u_{K,b}^{\infty} \delta_K \xi_K^n \cdot \mathbb{A}_K \delta_K \xi_K^n &\le \sum_{\sigma \in \mathcal{E}_K} 4u_{K,b}^{\infty} B_K^{\sigma\sigma} \left( \sqrt{\frac{u_K^n}{u_K^{\infty}}} - \sqrt{\frac{u_{\sigma}^n}{u_{\sigma}^{\infty}}} \right)^2 \\
&\le \sum_{\sigma \in \mathcal{E}_K} u_{K,b}^{\infty} B_K^{\sigma\sigma} m \left( \frac{u_K^n}{u_K^{\infty}}, \frac{u_{\sigma}^n}{u_{\sigma}^{\infty}} \right) \left( \log \left( \frac{u_K^n}{u_K^{\infty}} \right) - \log \left( \frac{u_{\sigma}^n}{u_{\sigma}^{\infty}} \right) \right)^2.
\end{aligned}
$$

By definition of $u_{K,b}^{\infty}$, and monotonicity (1.42a) and homogeneity (1.42c) of $m$, we infer that, for all $\sigma \in \mathcal{E}_K$,

$$ u_{K,b}^{\infty} m \left( \frac{u_K^n}{u_K^{\infty}}, \frac{u_{\sigma}^n}{u_{\sigma}^{\infty}} \right) \le u_{K,b}^{\infty} m \left( \frac{u_K^n}{u_{K,b}^{\infty}}, \frac{u_{\sigma}^n}{u_{K,b}^{\infty}} \right) = m(u_K^n, u_{\sigma}^n). $$

By definition (1.43) of $f|_{\mathcal{E}_K}$, and the bound (1.4) on $|\mathcal{E}_K|$, we then have

$$ \max_{\sigma \in \mathcal{E}_K} u_{K,b}^{\infty} m \left( \frac{u_K^n}{u_K^{\infty}}, \frac{u_{\sigma}^n}{u_{\sigma}^{\infty}} \right) \le \max_{\sigma \in \mathcal{E}_K} m \left( u_K^n, u_{\sigma}^n \right) \le |\mathcal{E}_K| r_K(\underline{u}_K^n) \le d\theta_D^2 r_K(\underline{u}_K^n). $$

We deduce that

$$ 
\begin{aligned}
4u_{K,b}^{\infty} \delta_K \xi_K^n \cdot \mathbb{A}_K \delta_K \xi_K^n &\le d\theta_D^2 r_K(\underline{u}_K^n) \sum_{\sigma \in \mathcal{E}_K} B_K^{\sigma\sigma} \left( \log \left( \frac{u_K^n}{u_K^{\infty}} \right) - \log \left( \frac{u_{\sigma}^n}{u_{\sigma}^{\infty}} \right) \right)^2 \\
&= d\theta_D^2 r_K(\underline{u}_K^n) \delta_K \underline{w}_K^n \cdot \mathbb{B}_K \delta_K \underline{w}_K^n,
\end{aligned}
$$

where $\underline{w}_K^n \in \underline{V}_K$ is such that $\underline{u}_K^n = \underline{u}_K^{\infty} \times \exp(\underline{w}_K^n)$. Using again (1.96) from Lemma 4, we finally infer that

$$ 4u_{K,b}^{\infty} \delta_K \xi_K^n \cdot \mathbb{A}_K \delta_K \xi_K^n \le d\theta_D^2 C_B r_K(\underline{u}_K^n) \delta_K \underline{w}_K^n \cdot \mathbb{A}_K \delta_K \underline{w}_K^n, $$

with $C_B > 0$ only depending on $\Lambda$, $\Omega$, $d$, and $\theta_D$. Summing over $K \in \mathcal{M}$, and recalling the definitions (1.70) of $\hat{\mathbb{D}}^n$, and (1.49) of $\mathbb{D}^n$, eventually yields (1.72) with $C_F = d\theta_D^2 C_B$. $\square$

The long-time behaviour of the nonlinear HFV scheme is studied in the following result.

**Theorem 4 (Asymptotic stability).** Recall the discrete entropy/dissipation relation of Proposition 5. The discrete entropy decays exponentially fast in time: there is $v_{nl} = \frac{4u_b^{\infty} \lambda_b \alpha_b}{C_F C_{LS,\infty}^2} > 0$, depending on $\Lambda$, $\phi$, $M$, $\Omega$, $d$, and $\theta_D$ such that

$$ \forall n \in \mathbb{N}, \quad \mathbb{E}^{n+1} \le (1 + v_{nl} \Delta t)^{-1} \mathbb{E}^n. \qquad (1.75) $$

Consequently, the discrete solution converges exponentially fast in time towards its associated discrete steady-state: for all $n \in \mathbb{N}^*$,

$$ \|u_M^n - u_M^\infty\|_{L^1(\Omega)} \le \sqrt{2M\mathbb{E}^0} (1 + v_{nl}\Delta t)^{-\frac{n}{2}}. \qquad (1.76) $$

**Proof.** Let $n \in \mathbb{N}$. By definition (1.70) of $\hat{\mathbb{D}}^n$, and from the coercivity estimate (1.14), we first infer that

$$ \hat{\mathbb{D}}^{n+1} = 4 \sum_{K \in \mathcal{M}} u_{K,b}^{\infty} a_K^{\Lambda}(\underline{\xi}_K^{n+1}, \underline{\xi}_K^{n+1}) \ge 4u_b^{\infty} a_D^{\Lambda}(\underline{\xi}_D^{n+1}, \underline{\xi}_D^{n+1}) \ge 4u_b^{\infty} \lambda_b \alpha_b |\underline{\xi}_D^{n+1}|_{1,D}^2 $$
which, combined with (1.72), implies that

$$ \mathbb{D}^{n+1} \geq \frac{1}{C_F} \hat{\mathbb{D}}^{n+1} \geq \frac{4u_b^{\infty} \lambda_b \alpha_b}{C_F} |\xi_D^{n+1}|_{1,D}^2 $$

In order to compare this quantity with the entropy, we use the discrete log-Sobolev inequality (1.86) from Proposition 9, applied to the couple $(u_D^{n+1}, u_D^{\infty})$ (which satisfies the mass condition owing to Proposition 4). We get

$$ \mathbb{E}^{n+1} = \int_{\Omega} u_{\mathcal{M}}^{\infty} \Phi_1 \left( \frac{u_{\mathcal{M}}^{n+1}}{u_{\mathcal{M}}^{\infty}} \right) \leq C_{LS,\infty}^2 |\xi_D^{n+1}|_{1,D}^2 $$

which, combined with the previous estimate, yields

$$ \mathbb{D}^{n+1} \geq \frac{4u_b^{\infty} \lambda_b \alpha_b}{C_F C_{LS,\infty}^2} \mathbb{E}^{n+1}. $$

Combined with (1.50) from Proposition 5, this shows (1.75). The $L^1$-norm estimate (1.76) is then a direct consequence of (1.75) and of the Csiszár-Kullback lemma (cf., e.g., [51, Lemma 5.6]) applied to the probability measure $\mu(x)dx = u_{\mathcal{M}}^{\infty}(x) \frac{dx}{M}$ and to the function $g = \frac{u_{\mathcal{M}}^{n}}{u_{\mathcal{M}}^{\infty}}$ such that $\int_{\Omega} g d\mu = 1$, which yields $\|u_{\mathcal{M}}^{n} - u_{\mathcal{M}}^{\infty}\|_{L^1(\Omega)} \leq \sqrt{2M\mathbb{E}^{n}}$ for all $n \geq 1$. $\square$

**Remark 9** (Norms and long-time behaviour). Notice that Theorem 4 states an exponential decay in $L^1$-norm, whereas Theorems 2 and 3 assert a convergence in $L^2$, and thus a convergence in $L^p$ for any $p \in [1, 2]$. This is reminiscent of the fact that the natural topologies for the linear and nonlinear problems differ.

# 1.5 Numerical results

## 1.5.1 Implementation

In this section, we discuss some practical aspects concerning the implementation of the schemes described in this paper. In all the test-cases presented below, the two-dimensional domain is taken to be $\Omega = (0, 1)^2$. The meshes used for the numerical tests, presented on Figure 1.2, are the classical Cartesian, triangular, and Kershaw meshes from the FVCA V benchmark (see [154]), as well as a tilted hexagonal-dominant mesh (cf. [101]). These meshes have convex

<figure><img src="image_17.png" /><figcaption>Figure 1.2 – Implementation. Coarsest meshes of each family used in the numerical tests.</figcaption></figure>

cells, hence we always choose $x_K$ to be the barycentre of $K \in \mathcal{M}$. In our implementation,
we compute the meshsize as $\tilde{h}_D = \max_{K \in \mathcal{M}} \frac{|K|}{|\partial K|}$. Observe that $\tilde{h}_D/h_D$ is framed by constants only depending on the mesh regularity. Notice also that the Kershaw mesh family is not uniformly regular in the sense defined in Section 1.2.1. In practice, we use a fixed value $\sqrt{2} < \eta = 1.5 < 2$ of the stabilisation parameter (see Remark 1). In the sequel, we denote by HMM the classical HFV linear scheme for advection-diffusion, and we restrict our attention to the Scharfetter-Gummel discretisation of the flux (1.18), namely to the function $A(s) = \frac{s}{e^{s-1}} - 1$, extended by continuity to 0 at $s = 0$.

## Linear systems and static condensation

The two linear (HMM and exponential fitting) schemes are implemented in the same way. To fix ideas, we consider the evolution problem with pure Neumann boundary conditions, of unknown solution $u_D \in V_D$. We denote by $U_M \in \mathbb{R}^{|M|}$ and $U_E \in \mathbb{R}^{|\mathcal{E}|}$ the unknown vectors $(u_K)_{K \in \mathcal{M}}$ and $(u_\sigma)_{\sigma \in \mathcal{E}}$. The linear schemes result in the following block system:

$$ \left( \begin{array}{cc} M_M & M_{M,\mathcal{E}} \\ M_{\mathcal{E},M} & M_{\mathcal{E}} \end{array} \right) \left( \begin{array}{c} U_M \\ U_{\mathcal{E}} \end{array} \right) = \left( \begin{array}{c} S_M \\ S_{\mathcal{E}} \end{array} \right), \qquad (1.77) $$

where $M_M \in \mathbb{R}^{|M|\times|M|}$, $M_{M,\mathcal{E}} \in \mathbb{R}^{|M|\times|\mathcal{E}|}$, $M_{\mathcal{E},M} \in \mathbb{R}^{|\mathcal{E}|\times|M|}$, $M_{\mathcal{E}} \in \mathbb{R}^{|\mathcal{E}|\times|\mathcal{E}|}$, and $S_M \in \mathbb{R}^{|M|}$ and $S_{\mathcal{E}} \in \mathbb{R}^{|\mathcal{E}|}$ stem from the loading term and the boundary conditions. By construction, the matrix $M_M$ is diagonal with non-zero diagonal entries, and can therefore be inverted at a very low computational cost. Thus, one can eliminate the cell unknowns, noticing that

$$ U_M = M_M^{-1} (S_M - M_{M,\mathcal{E}} U_{\mathcal{E}}). \qquad (1.78) $$

Using this relation, one infers that $U_E$ is the solution to the following linear system:

$$ (M_{\mathcal{E}} - M_{\mathcal{E},M} M_M^{-1} M_{M,\mathcal{E}}) U_{\mathcal{E}} = S_{\mathcal{E}} - M_{\mathcal{E},M} M_M^{-1} S_M, \qquad (1.79) $$

where $M_D = M_{\mathcal{E}} - M_{\mathcal{E},M} M_M^{-1} M_{M,\mathcal{E}}$, the so-called Schur complement of the matrix $M_M$, is an invertible matrix of size $|\mathcal{E}| \times |\mathcal{E}|$. In practice, we solve the linear system (1.79) using an LU factorisation algorithm, and we use the solution $U_E$ to reconstruct $U_M$ from (1.78). This method, called *static condensation*, allows one to replace a system of size $|M| + |\mathcal{E}|$ by a system of size $|\mathcal{E}|$ without additional fill-in. In the case of mixed Dirichlet-Neumann boundary conditions, the Dirichlet face unknowns are eliminated from the global linear system.

## Exponential fitting scheme: choice of unknown and harmonic averaging

The exponential fitting scheme can be expressed in either the $u$ or the $\rho = u e^{\phi}$ variable. In the $\rho$ variable, the resulting linear system is symmetric. One can then use, e.g., Cholesky factorisation or a conjugate gradient method. However, the formulation in $\rho$ is ill-conditioned. In our numerical experiments, the ratio between the condition numbers of the linear systems in $\rho$ and in $u$ often exceeds $10^5$. Because of this, we chose and we recommend to solve the linear system in the unknown $u$. Notice that solving the system in $u$ is equivalent to right pre-condition the system in $\rho$ with the inverse of the diagonal matrix with entries the coordinates of the interpolate of $\omega = e^{-\phi}$.

In order to implement the exponential fitting scheme, one needs to evaluate averages of the diffusion tensor $\frac{1}{|P_{K,\sigma}|} \int_{P_{K,\sigma}} \omega \Lambda$. Observe that $\omega(x)/\omega(x_K)$ is of order $e^{h_K \|\nabla\phi\|_{L^\infty(p_{K,\sigma})}}$ in $P_{K,\sigma}$. Therefore, with large advection fields, the diffusion problem (1.28) becomes strongly heterogeneous. It
is pointed out in [43] that an (empirical) solution to improve robustness to this heterogeneity is to use harmonic averages to approximate integrals of the diffusion tensor. In the numerical tests of the following subsections, we compare the “classical” exponential fitting scheme (for which the integral is approximated by a standard - second order - quadrature) with the “harmonic” one, in which case we choose to use the following approximation:

$$ \frac{1}{|P_{K,\sigma}|} \int_{P_{K,\sigma}} \omega \Lambda \approx 3 \left( \sum_{F \in \mathcal{E}_{P_{K,\sigma}}} \frac{1}{\omega(\bar{x}_F)} \right)^{-1} \Lambda(x_K), $$

where $\mathcal{E}_{P_{K,\sigma}}$ denotes the set of edges of the triangle $P_{K,\sigma}$ (recall that $d = 2$ in our experiments), and $\bar{x}_F$ is the barycentre of $F \in \mathcal{E}_{P_{K,\sigma}}$.

## Nonlinear scheme and Newton's method

The implementation of the nonlinear scheme relies on the following formulation: given $u_D^{n-1} \in \underline{V}_D$ positive, we want to solve the nonlinear system $\underline{G}_D^{n,\delta t}(u_D^n) = \underline{0}_D$, where $\underline{G}_D^{n,\delta t}$ is defined as in (1.55) but with a time step $\delta t$ instead of $\Delta t$. The resolution of this system relies on Newton's method.

First, one initialises the method with $\delta t = \Delta t$, and initial guess $\max(u_D^{n-1}, \epsilon \underline{1}_D) \in \underline{V}_D$ (where the maximum is taken coordinate by coordinate), in order to avoid potential problems due to the singularity of the log near 0. The successive linear systems to compute the residue have the same structure as (1.77). We thus perform static condensation at each Newton iteration. As a stopping criterion, we compare the $l^\infty$ relative norm of the residue with a threshold tol. If the method does not converge after $i_{max}$ iterations, we divide the time step by 2, and we restart the resolution. When the method converges, one can proceed with the approximation of $u_D^{n+1}$, with an initial time step of $\min(\Delta t, 2\delta t)$. In practice, we use $\epsilon = 10^{-11}$, $i_{max} = 50$, and tol $= 10^{-11}$.

The implementation of the nonlinear scheme relies on the computation of $\log(\omega_K)$ and $\log(\omega_\sigma)$. Since we have chosen $x_K$ to be the barycentre of $K$, we choose to approximate $\frac{1}{|K|} \int_K e^{-\phi}$ by $e^{-\phi(x_K)}$. Therefore, $\log(\omega_K)$ is computed as $\log(e^{-\phi(x_K)}) = -\phi(x_K)$. The same holds true for $\log(\omega_\sigma)$.

In the simulations shown below, we use arithmetic means for the functions $m$ and $f_{|\mathcal{E}_K|}$ of the reconstruction $r_K$ defined by (1.41)-(1.43). For all $K \in \mathcal{M}$, and all $\underline{u}_K \in \underline{V}_K$, we thus consider

$$ r_K(\underline{u}_K) = \frac{1}{2} \left( u_K + \frac{1}{|\mathcal{E}_K|} \sum_{\sigma \in \mathcal{E}_K} u_\sigma \right). $$

This choice is close to the one advocated in [49, Eq. (58)]. For a discussion on other choices of reconstructions, we refer to [51, Section 6.2].

## 1.5.2 Long-time behaviour of discrete solutions

In this section, we present some numerical illustration of the long-time behaviour of discrete solutions. We focus on a test-case from [56, 52, 51]. We consider homogeneous pure Neumann boundary conditions ($\Gamma^D = \emptyset$ and $g^N = 0$), and zero loading term ($f = 0$). The advective potential and diffusion tensor are set to $\phi(x, y) = -x$ and $\Lambda = \begin{pmatrix} l_x & 0 \\ 0 & 1 \end{pmatrix}$ for $l_x > 0$. The exact solution is given
by

$$u(t, x, y) = C_1 e^{-\alpha t + \frac{\beta}{2}} (2\pi \cos(\pi x) + \sin(\pi x)) + 2C_1 \pi e^{x - \frac{1}{2}},$$

where $C_1 > 0$ and $\alpha = l_x(\frac{1}{4} + \pi^2)$. Note that $u^{in}$ vanishes on $\{x = 1\}$, but for any $t > 0$, $u(t, \cdot) > 0$. The associated steady-state is

$$u^{\infty}(x, y) = 2C_1 \pi e^{x - \frac{1}{2}}.$$

Our experiments are performed using the following values:

$$l_x = 10^{-2} \quad \text{and} \quad C_1 = 10^{-1}.$$

We compute the solution on the time interval $[0, T_f]$, and we denote by $(u_D^n)_{1 \le n \le N_f}$ the corresponding approximate solution. Note that the number of time steps $N_f$ may differ between the linear and nonlinear schemes, because of the adaptive time step refinement procedure used for the nonlinear scheme.

We set $T_f = 350$ in order to see the complete evolution. Since the long-time behaviour of the schemes does not depend on the size of the discretisation, we can explore the evolution using a large time step $\Delta t = 10^{-1}$. We perform the numerical experiments on two Kershaw meshes (see Figure 1.2c) of sizes 0.02 and 0.006. In Figure 1.3, we depict, as a function of time, the $L^1$

<figure><img src="image_18.png" /><figcaption>Figure 1.3 – Long-time behaviour of discrete solutions. Comparison of the long-time behaviour on Kershaw meshes for $T_f = 350$ and $\Delta t = 0.1$.</figcaption></figure>

distance between $u_D^n$ and $u^\infty$ (the exact steady-state) computed as

$$\sum_{K \in \mathcal{M}} |K||u_K^n - u^\infty(x_K)|.$$

We observe the exponential decay towards the steady-state, until some precision is reached. Note that for the HMM scheme, some saturation occurs at precision of magnitude $10^{-6}$ and $10^{-7}$: the scheme does not preserve the thermal equilibrium (see Remark 4). This saturation corresponds to the accuracy of the stationary scheme (see Section 1.5.4), so the threshold is lower on the refined mesh. Note that one could also consider the error measure $\sum_{K \in \mathcal{M}} |K||u_K^n - u^\infty|$
between the discrete solution and the discrete steady-state: this quantity decays exponentially, with a lower saturation of magnitude $10^{-12}$, corresponding to machine precision. The other (nonlinear and exponential fitting) schemes have the same decay rate, and the saturation occurs at machine precision. For the four schemes, the rates of convergence are similar to the real one $\alpha$. In particular, the use of harmonic averages in the exponential fitting scheme does not have any impact on the long-time behaviour.

### 1.5.3 Positivity of discrete solutions

We are now interested in the positivity of the discrete solutions. We use the following test-case with anisotropic diffusion and homogeneous pure Neumann boundary conditions. We set $\Gamma^D = \emptyset$, $f = 0$, $g^N = 0$,

$$ \phi(x, y) = -(x - 0.4)^2 + (y - 0.6)^2 \quad \text{and} \quad \Lambda = \begin{pmatrix} 0.8 & 0 \\ 0 & 1 \end{pmatrix}. $$

For the initial datum, we take

$$ u^{in} = 10^{-3} \mathbb{1}_{B} + \mathbb{1}_{\Omega \setminus B}, $$

where $B$ is the Euclidean ball $\{(x, y) \in \mathbb{R}^2 | (x - 0.5)^2 + (y - 0.5)^2 \le 0.2^2\}$. These data ensure that the solution $u$ is positive on $\mathbb{R}_+ \times \Omega$. The experiment is performed on a tilted hexagonal-dominant mesh (see Figure 1.2d) of size $4.3 \cdot 10^{-3}$, made up of 4192 cells and 12512 edges. Since we deal with a diffusive phenomenon, the smallest values of $u$ are expected for small time, hence we perform the simulation with a relatively small final time $T_f = 5.10^{-4}$, alongside with a time step of $\Delta t = 10^{-5}$.

The results are collected in Table 1.1. The cost is defined as the number of linear systems solved in order to compute the solution $(u_D^n)_{1 \le n \le N_f}$, and the minimum values $\min\_cells$ and $\min\_edges$ are defined by

$$ \min\{u_K^n | 1 \le n \le N_f, K \in \mathcal{M}\} \quad \text{and} \quad \min\{u_\sigma^n | 1 \le n \le N_f, \sigma \in \mathcal{E}\}. $$

The indicated number of negative unknowns is for the whole simulation. Here, $N_f = 50$ also for the nonlinear scheme (no sub-division of the time step was needed in Newton's method). As

<table><caption>Table 1.1 – Positivity of discrete solutions. Numerical results for $T_f = 5.10^{-4}$ and $\Delta t = 10^{-5}$ on a tilted hexagonal-dominant mesh. At each time step, there are 4192 cell unknowns and 12512 edge unknowns.</caption><thead><tr><th></th><th>cost</th><th>min</th><th>min_edges</th><th># negative unknowns</th></tr></thead><tbody><tr><td>Nonlinear</td><td>175</td><td>9.93e-04</td><td>7.36e-04</td><td>0</td></tr><tr><td>HMM</td><td>50</td><td>-5e-03</td><td>-7.74e-02</td><td>593</td></tr><tr><td>ExpF</td><td>50</td><td>-4.98e-03</td><td>-7.72e-02</td><td>590</td></tr><tr><td>ExpF (harmonic)</td><td>50</td><td>-4.98e-03</td><td>-7.74e-02</td><td>588</td></tr></tbody></table>

expected, the nonlinear scheme has positive discrete solutions, whereas the linear ones exhibit a violation of positivity (the value of $\eta$ can have some influence on positivity; see [122]). Note that the use of harmonic averages for the exponential fitting scheme has no impact on the undershoots.

We observe that the nonlinear scheme requires approximately 3.5 times more linear system inversions than the linear schemes. However, this value depends strongly on the final time
of simulation $T_f$. Indeed, the number of linear systems solved at step $n$ decreases when $n$ increases. The first time step costs 9 resolutions, but this number rapidly decreases as the solution approaches the steady-state (the second and the third time steps respectively cost 5 and 4 resolutions).

### 1.5.4 Accuracy of stationary solutions

In this section, we aim at comparing the accuracy of the different schemes for the stationary problem. To do so, we define the discrete $L^2$-norm and $H^1$-seminorm errors as (i) the $L^2$-norm of the difference $u_M - \Pi_M(u)$, and (ii) the $|\cdot|_{1,D}$-seminorm of the difference $\underline{u}_D - \bar{\Pi}_D(u)$, where $\underline{u}_D$ is the discrete solution, and $\bar{\Pi}_D(u)$ is the HFV interpolate of the continuous solution $u$ (computed as $\bar{\Pi}_D(u) \approx ((u(x_K))_{K \in M}, (u(\bar{x}_\sigma))_{\sigma \in \mathcal{E})))$. In what follows, we reason in relative errors.

The nonlinear scheme is extended to a more general setting, in order to take into account a loading term $f \ge 0$ and mixed Dirichlet-Neumann boundary conditions ($|\Gamma^D| > 0$) with $g^D > 0$ and $g^N \ge 0$. The scheme writes:

$$ \text{Find } \underline{u}_D \in \underline{V}_D \text{ positive such that } \underline{G}_D(\underline{u}_D) = \underline{0}_D, \qquad (1.80) $$

where $\underline{G}_D : (\underline{V}_D)^* \to \underline{V}_D$ is the vector field defined by

$$ \underline{G}_K(\underline{u}_D) = \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{K,\sigma}^{\text{nl}}(\underline{u}_K) - \int_K f \qquad \forall K \in \mathcal{M}, \qquad (1.81\text{a}) $$

$$ \underline{G}_{\sigma}(\underline{u}_{D}) = - (\underline{\mathcal{F}}_{K,\sigma}^{\text{nl}}(\underline{u}_{K}) + \underline{\mathcal{F}}_{L,\sigma}^{\text{nl}}(\underline{u}_{L})) \qquad \forall \sigma = K \mid L \in \mathcal{E}_{\text{int}}, \qquad (1.81\text{b}) $$

$$ \underline{G}_{\sigma}(\underline{u}_D) = - \int_{\sigma} g^N - \mathcal{F}_{K,\sigma}^{\text{nl}}(\underline{u}_K) \qquad \forall \sigma \in \mathcal{E}_{\text{ext}}^N \text{ with } \mathcal{M}_{\sigma} = \{K\}, \qquad (1.81c) $$

$$ \underline{G}_{\sigma}(\underline{u}_D) = \frac{1}{|\sigma|} \int_{\sigma} g^D - u_{\sigma} \qquad \forall \sigma \in \mathcal{E}_{\text{ext}}^D, \qquad (1.81\text{d}) $$

$$ \underline{\mathcal{F}}_{K,\sigma}^{\text{nl}}(\underline{u}_K) = r^K(\underline{u}_K) \sum_{\sigma' \in \mathcal{E}_K} A_K^{\sigma\sigma'} (\log(u_K) + \phi(x_K) - \log(u_{\sigma'}) - \phi(\bar{x}_{\sigma'})) \qquad \forall K \in \mathcal{M}, \forall \sigma \in \mathcal{E}_K. \qquad (1.81\text{e}) $$

The implementation of this scheme still relies on a Newton method similar to the one used for the evolution scheme. It is here initialised with $\underline{\Pi}_D(\omega) = \underline{\omega}_D$.

The first test-case we consider is the same as in [19]. It is an isotropic problem, with $\Lambda = I_2$, $\Gamma^N = \emptyset$, $\phi(x, y) = -(2x + 3y)$, and exact solution

$$ u(x, y) = (x - e^{2(x-1)}) (y^2 - e^{3(y-1)}), $$

the other data $f$ and $g^D$ being set accordingly. Note that for this test-case, the diffusion and advection terms are of the same order of magnitude. The numerical experiments are performed on the triangular mesh family (see Figure 1.2b). The convergence results are depicted in Figure 1.4. As expected, the two linear schemes are of order two in $L^2$-norm, and one in $H^1$-seminorm. The same holds for the nonlinear scheme, whose accuracy is rather the same as the classical HMM scheme, one order of magnitude better than the exponential fitting schemes. On this test-case, the use of harmonic averages for the exponential fitting scheme does not have a significant impact.

The second test-case is an advection-dominated problem, with anisotropic diffusion and
<figure><img src="image_19.png" /><figcaption>Figure 1.4 – Accuracy of stationary solutions. Relative errors in discrete L²-norm and H¹-seminorm for the first test-case on triangular meshes.</figcaption></figure>

mixed Dirichlet-Neumann boundary conditions. We set $\Gamma^D = (\{0\} \times (0,1)) \cup (\{1\} \times (0,1))$, $\Gamma^N = ((0,1) \times \{0\}) \cup ((0,1) \times \{1\})$, $g^D = 1$, $g^N = 0$, and $f = 0$. The diffusion tensor and the potential are defined by the following expressions:

$$ \Lambda = \begin{pmatrix} 1 & 0 \\ 0 & l_y \end{pmatrix} \quad \text{and} \quad \phi(x, y) = \log\left(\frac{1}{v} + x\right), $$

with $l_y, v > 0$. Note that the advection field $V^\phi = -\left(\frac{v}{1+vx}\right)$ has a magnitude of order $v$ when $x$ is small. Thus, near the boundary $\{0\} \times [0,1]$, the problem is advection-dominated if $v$ is large enough. Moreover, $\text{div}(V^\phi) = \frac{v^2}{(1+vx)^2} > 0$ and $V^\phi \cdot n = 0$ on $\Gamma^N$, so the problem is coercive. The exact solution is given by

$$ u(x, y) = \frac{v}{1+vx} \left( \frac{2vx}{2+v} \left( \frac{1}{v} + \frac{x}{2} \right) + \frac{1}{v} \right). $$

We perform our numerical experiments on the Cartesian mesh family (see Figure 1.2a), with

$$ l_y = 100 \quad \text{and} \quad v = 200. $$

The results are depicted in Figure 1.5. They show that the HMM scheme suffers, most probably because of the fact that the advective term predominates over the diffusive term, at least in some part of the domain. The order of convergence of the HMM scheme is less than one in $H^1$-seminorm, and than two in $L^2$-norm. The other schemes converge with order one in $H^1$-seminorm, and two in $L^2$-norm. Moreover, on this test-case, their accuracy is better than that of the HMM scheme. Notice that $\omega = e^{-\phi} = \frac{v}{1+vx}$ has small variations (i.e., not exponential) in the cells, even if $v$ is large. Therefore, the diffusion tensor $\omega\Lambda$ of the problem in the $\rho$ unknown for the exponential fitting schemes is not that heterogeneous (locally). It could explain the good performances of the exponential fitting schemes in this case. Moreover, on this test-case, using harmonic averages in the exponential fitting scheme gives a substantial gain in accuracy for both
<figure><img src="image_20.png" /><figcaption>Figure 1.5 – Accuracy of stationary solutions. Relative errors in discrete L²-norm and H¹-seminorm for the second test-case on Cartesian meshes.</figcaption></figure>

the L² and H¹ relative errors, of magnitude 10¹.

## 1.6 Conclusion

In this paper, by means of discrete entropy methods, we have analysed the long-time behaviour of three hybrid finite volume schemes for linear advection-diffusion equations. We have proved that the solutions to all schemes converge exponentially fast in time towards the associated discrete steady-states. Two schemes among the three are new, that are the (linear) exponential fitting scheme (adapting known ideas to the HFV context) and the nonlinear scheme, for which we have proved the existence of solutions. All schemes can handle anisotropy and general meshes. The two linear schemes can deal with general data and mixed Dirichlet-Neumann boundary conditions, however they do not preserve the positivity of solutions. On the other hand, the nonlinear scheme preserves positivity and can be used in practice with general boundary conditions. However, at the moment, its asymptotic analysis is limited to systems that converge in time towards the thermal equilibrium, restricting the admissible data. We have finally validated our theoretical findings on different numerical tests, assessing long-time behaviour, positivity, and spatial accuracy of the schemes.

### 1.A Functional inequalities

#### 1.A.1 Discrete Poincaré inequalities

We recall the following hybrid discrete Poincaré inequalities (cf. [111, Lemmas B.25 and B.32, p = 2]).

**Proposition 7 (Discrete Poincaré inequalities).** Let $\mathcal{D}$ be a given discretisation of $\Omega$, with regularity parameter $\theta_D$. There exists $C_{PW} > 0$, only depending on $\Omega$, $d$, and $\theta_D$ such that

$$ \forall \underline{v}_D \in \underline{V}_{D,0}^N, \quad \|v_M\|_{L^2(\Omega)} \le C_{PW}|\underline{v}_D|_{1,D}. \qquad (1.82) $$
Assume that $|\Gamma^D| > 0$. Then, there exists $C_{P,\Gamma^D} > 0$, only depending on $\Omega$, $d$, $\Gamma^D$, and $\theta_D$ such that

$$ \forall \underline{v}_D \in \underline{V}_{D,0}^D, \quad \|v_M\|_{L^2(\Omega)} \le C_{P,\Gamma^D} |v_D|_{1,D}. \qquad (1.83) $$

### 1.A.2 Logarithmic Sobolev inequalities

In this section, we derive logarithmic Sobolev inequalities on a bounded domain, in the continuous setting. The intermediate results of Proposition 8 below will be useful in the discrete setting. In the following, $\mu$ is a probability measure on the bounded domain $\Omega$, and the space $L_{\mu}^q(\Omega)$ denotes the space endowed with the norm $\|f\|_{L_{\mu}^q(\Omega)}^q = \int_{\Omega} |f|^q d\mu$. We start with a preliminary lemma, which is an adaptation of part of the proof of [96, Theorem 6.1.22] (see also [153]). We recall that $\Phi_1(s) = s \log(s) - s + 1$.

**Lemma 3.** For all $t \in \mathbb{R}$ and $\psi \in L_{\mu}^2(\Omega)$ such that $\|\psi\|_{L_{\mu}^2(\Omega)} = 1$, one has

$$ \int_{\Omega} \Phi_1((1+t\psi)^2) d\mu \le t^2 \int_{\Omega} \psi^2 \log(\psi^2) d\mu + (1+t^2) \log(1+t^2) + (1+|\langle\psi\rangle_{\mu}|)t^2, $$

where $\langle\psi\rangle_{\mu} = \int_{\Omega} \psi d\mu$.

*Proof.* Let us define, for $\delta > 0$,

$$ f_{\delta}(t) = \int_{\Omega} \Phi_{1}((1+t\psi)^2 + \delta) d\mu - t^2 \int_{\Omega} \psi^2 \log(\psi^2) d\mu - (1+t^2) \log(1+t^2). $$

Differentiating $f_{\delta}$ yields

$$ f'_{\delta}(t) = 2 \int_{\Omega} (1+t\psi)\psi \log((1+t\psi)^2 + \delta) d\mu - 2t \int_{\Omega} \psi^2 \log(\psi^2) d\mu - 2t \log(1+t^2) - 2t. $$

In particular, $f'_{\delta}(0) = 2\log(1+\delta)\langle\psi\rangle_{\mu}$. Differentiating once more, and using that $\|\psi\|_{L_{\mu}^2(\Omega)}^2 = 1$, we obtain

$$ f_{\delta}''(t) = 2 \int_{\Omega} \psi^2 \log \left( \frac{(1+t\psi)^2 + \delta}{(1+t^2)\psi^2} \right) d\mu + 4 \int_{\Omega} \psi^2 \frac{(1+t\psi)^2}{\delta + (1+t\psi)^2} d\mu - \frac{4t^2}{1+t^2} - 2. $$

Therefore, using that $\log(x) \le x - 1$ in the first term, that $\frac{x}{\delta+x} \le 1$ in the second, together with the fact that $\mu$ is a probability measure and that $\|\psi\|_{L_{\mu}^2(\Omega)}^2 = 1$, one gets

$$ f_{\delta}''(t) \le \frac{2\delta}{1+t^2} + \frac{4t}{1+t^2}\langle\psi\rangle_{\mu} + 4 - \frac{4t^2}{1+t^2} - 2 \le 2\delta + 2|\langle\psi\rangle_{\mu}| + 2. $$

One concludes by integrating this inequality twice between 0 and $t$, using that $f_{\delta}(0) = \Phi_1(1+\delta)$, and letting $\delta \to 0$. $\square$

**Proposition 8.** For any $\xi \in L_{\mu}^q(\Omega)$ with $q > 2$, one has

$$ \int_{\Omega} \xi^2 \log \left( \frac{\xi^2}{\|\xi\|_{L_{\mu}^2(\Omega)}^2} \right) d\mu \le \frac{q}{q-2} \|\xi - \langle\xi\rangle_{\mu}\|_{L_{\mu}^2(\Omega)}^2 + \frac{q-4}{q-2} \|\xi - \langle\xi\rangle_{\mu}\|_{L_{\mu}^2(\Omega)}^2, \quad (1.84) $$
where $\langle \xi \rangle_{\mu} = \int_{\Omega} \xi d\mu$. Besides, one also has

$$ \int_{\Omega} \Phi_1(\xi^2) d\mu \le \frac{q}{q-2} \|\xi - 1\|_{L_{\mu}^2(\Omega)}^2 + \frac{2q-6}{q-2} \|\xi - 1\|_{L_{\mu}^2(\Omega)}^2 + \Phi_1 \left( 1 + \|\xi - 1\|_{L_{\mu}^2(\Omega)}^2 \right). \quad (1.85) $$

**Proof.** i) Assume that $\langle \xi \rangle_{\mu} \ne 0$, and take $t$ and $\psi$ such that $\xi = \langle \xi \rangle_{\mu}(1+t\psi)$ and $\|\psi\|_{L_{\mu}^2(\Omega)} = 1$.

In particular, $\langle \psi \rangle_{\mu} = 0$, and $1+t^2 = \frac{\|\xi\|_{L_{\mu}^2(\Omega)}^2}{\langle \xi \rangle_{\mu}^2}$. Using Lemma 3, a somewhat tedious but straightforward computation yields

$$ \int_{\Omega} \xi^2 \log \left( \frac{\xi^2}{\|\xi\|_{L_{\mu}^2(\Omega)}^2} \right) d\mu \le \int_{\Omega} (\xi - \langle \xi \rangle_{\mu})^2 \log \left( \frac{(\xi - \langle \xi \rangle_{\mu})^2}{\|\xi - \langle \xi \rangle_{\mu}\|_{L_{\mu}^2(\Omega)}^2} \right) d\mu + 2\|\xi - \langle \xi \rangle_{\mu}\|_{L_{\mu}^2(\Omega)}^2. $$

Observe that the last inequality also holds if $\langle \xi \rangle_{\mu} = 0$. Let $\phi = \xi - \langle \xi \rangle_{\mu}$. Then,

$$ \int_{\Omega} \xi^2 \log \left( \frac{\xi^2}{\|\xi\|_{L_{\mu}^2(\Omega)}^2} \right) d\mu \le \frac{2}{q-2} \|\phi\|_{L_{\mu}^2(\Omega)}^2 \int_{\Omega} \frac{\phi^2}{\|\phi\|_{L_{\mu}^2(\Omega)}^2} \log \left( \frac{\phi^{q-2}}{\|\phi\|_{L_{\mu}^2(\Omega)}^2} \right) d\mu + 2\|\phi\|_{L_{\mu}^2(\Omega)}^2. $$

Therefore, by Jensen's inequality for the probability measure $\frac{\phi^2}{\|\phi\|_{L_{\mu}^2(\Omega)}^2} d\mu$ applied to the concave function log, one obtains

$$ \begin{aligned} \int_{\Omega} \xi^2 \log \left( \frac{\xi^2}{\|\xi\|_{L_{\mu}^2(\Omega)}^2} \right) d\mu &\le \frac{2}{q-2} \|\phi\|_{L_{\mu}^2(\Omega)}^2 \log \left( \frac{\|\phi\|_{L_{\mu}^2(\Omega)}^q}{\|\phi\|_{L_{\mu}^2(\Omega)}^q} \right) + 2\|\phi\|_{L_{\mu}^2(\Omega)}^2 \\ &= \frac{q}{q-2} \|\phi\|_{L_{\mu}^2(\Omega)}^2 \log \left( \frac{\|\phi\|_{L_{\mu}^2(\Omega)}^q}{\|\phi\|_{L_{\mu}^2(\Omega)}^q} \right) + 2\|\phi\|_{L_{\mu}^2(\Omega)}^2, \end{aligned} $$

and one concludes using that $\log(x) \le x - 1$.

ii) Take $t$ and $\psi$ such that $\xi = 1+t\psi$ and $\|\psi\|_{L_{\mu}^2(\Omega)} = 1$. Remark that $t = \|\xi - 1\|_{L_{\mu}^2(\Omega)}$. Using that $|\langle \psi \rangle_{\mu}| \le \|\psi\|_{L_{\mu}^2(\Omega)} = 1$, Lemma 3 yields

$$ \int_{\Omega} \Phi_1(\xi^2) d\mu - \Phi_1 \left( 1 + \|\xi - 1\|_{L_{\mu}^2(\Omega)}^2 \right) \le \int_{\Omega} (\xi - 1)^2 \log \left( \frac{(\xi - 1)^2}{\|\xi - 1\|_{L_{\mu}^2(\Omega)}^2} \right) d\mu + 3\|\xi - 1\|_{L_{\mu}^2(\Omega)}^2. $$

Letting $\phi = \xi - 1$, the proof goes on as for i). □

<figure><img src="image_21.png" /></figure>

From there, logarithmic Sobolev inequalities are immediate consequences of Poincaré-Sobolev inequalities, of [51, Lemma 5.2], and of the fact that $\Phi_1(1+s) \le s \log(1+s)$.

**Corollary 1** (Logarithmic Sobolev inequalities). Assume that $\mu$ has a density (still denoted by $\mu$) with respect to the Lebesgue measure such that $0 < \mu_b \le \mu(x) \le \mu_\#$ for a.e. $x \in \Omega$. Then, for any
$\xi \in H^1(\Omega)$, one has

$$ \int_{\Omega} \xi^2 \log \left( \frac{\xi^2}{\|\xi\|_{L_{\mu}^2(\Omega)}^2} \right) d\mu \le C(\Omega, d, \mu_b, \mu_\#) \|\nabla \xi\|_{L_{\mu}^2(\Omega; \mathbb{R}^d)}^2. $$

Besides, if $|\Gamma^D| > 0$ and $\xi - 1 \in H_0^{1,D}(\Omega) = \{v \in H^1(\Omega) | v = 0 \text{ on } \Gamma^D\}$, then

$$ \int_{\Omega} \Phi_1(\xi^2) d\mu \le C'(\Omega, d, \mu_b, \mu_\#) \left(1 + \log\left(1 + \|\xi - 1\|_{L_{\mu}^2(\Omega)}^2\right)\right) \|\nabla \xi\|_{L_{\mu}^2(\Omega; \mathbb{R}^d)}^2. $$

### 1.A.3 Discrete logarithmic Sobolev inequalities

Similarly to what was done in [51], one can derive discrete logarithmic Sobolev inequalities adapted to the hybrid setting.

**Proposition 9** (Discrete logarithmic Sobolev inequality, Neumann case). Let $\mathcal{D}$ be a given discretisation of $\Omega$, with regularity parameter $\theta_D$. Let $\underline{v}_D, \overline{v}_D \in \underline{V}_D$ be two positive vectors of unknowns such that

$$ \int_{\Omega} v_M = \int_{\Omega} v_M^\infty = M, $$

and set $v_{M,\#}^\infty = \sup_{K \in \mathcal{M}} v_K^\infty$. Define $\underline{\xi}_D$ as the element of $\underline{V}_D$ such that

$$ \underline{\xi}_K = \sqrt{\frac{v_K}{v_K^\infty}} \quad \forall K \in \mathcal{M}, \quad \overline{\xi}_\sigma = \sqrt{\frac{v_\sigma}{v_\sigma^\infty}} \quad \forall \sigma \in \mathcal{E}. $$

Then, there exists $C_{LS,\infty} > 0$, only depending on $M$, $v_{M,\#}^\infty$, $\Omega$, $d$, and $\theta_D$ such that

$$ \int_{\Omega} v_M^\infty \Phi_1(\xi_M^2) \le C_{LS,\infty} |\underline{\xi}_D|_{1,D}^2. \quad (1.86) $$

*Proof.* By (1.84) and [51, Lemma 5.2] applied to the probability measure $\mu(x)dx = v_M^\infty(x) \frac{dx}{M}$ and to the function $\underline{\xi}_M = \sqrt{\frac{v_M}{v_M^\infty}}$, we first infer that, for $q > 2$,

$$ \int_{\Omega} v_M \log(\xi_M^2) \le C(M, v_{M,\#}^\infty, q) \left( \|\underline{\xi}_M - \overline{\xi}_M\|_{L^q(\Omega)}^2 + \|\underline{\xi}_M - \overline{\xi}_M\|_{L^2(\Omega)}^2 \right), $$

where we let $\overline{\xi}_M = \frac{1}{|\Omega|} \int_{\Omega} \underline{\xi}_M dx$. The conclusion then falls in two steps. On the one hand, since $\int_{\Omega} v_M = \int_{\Omega} v_M^\infty$, we remark that

$$ \int_{\Omega} v_M^\infty \Phi_1(\xi_M^2) = \int_{\Omega} v_M \log(\xi_M^2). $$

On the other hand, we invoke (1.82) and the discrete Poincaré–Sobolev inequality of [111, Lemma B.25, p = 2] for $2 < q < \frac{2d}{d-2}$:

$$ \forall \underline{w}_D \in \underline{V}_{D,0}^N, \quad \|w_M\|_{L^q(\Omega)} \le C_{PS}|\underline{w}_D|_{1,D}, $$
where $C_{PS} > 0$ only depends on $\Omega$, $d$, and $\theta_D$, that we apply to $w_D = \underline{\xi}_D - \overline{\xi_M} \underline{1}_D \in \underline{V}_{D,0}^N$. This
proves (1.86).

<figure><img src="image_22.png" /></figure>

Starting from (1.85), a similar proof yields the following result. The relevant discrete
Poincaré–Sobolev inequality in this case is given in [111, Lemma B.32, p = 2].

**Proposition 10** (Discrete logarithmic Sobolev inequality, Dirichlet case). Assume that $|\Gamma^D| > 0$. Let $\mathcal{D}$ be a given discretisation of $\Omega$, with regularity parameter $\theta_D$. Let $v_D, v_D^\infty \in \underline{V}_D$ be two positive vectors of unknowns such that

$$
\underline{v}_D - \underline{v}_D^\infty \in \underline{V}_{D,0}^D,
$$

and set $v_{M,\#}^{\infty} = \sup_{K \in \mathcal{M}} v_K^{\infty}$ and $M^{\infty} := \int_{\Omega} v_M^{\infty}$. Define $\underline{\xi}_D$ as the element of $\underline{V}_D$ such that

$$
\xi_K = \sqrt{\frac{v_K}{v_K^\infty}} \quad \forall K \in \mathcal{M}, \quad \xi_\sigma = \sqrt{\frac{v_\sigma}{v_\sigma^\infty}} \quad \forall \sigma \in \mathcal{E}.
$$

Then, letting $\mu = \frac{v_M^\infty}{M^\infty}$, there exists $C_{LS,\Gamma^D,\infty} > 0$, only depending on $M^\infty$, $v_{M,\#}^\infty$, $\Omega$, $d$, $\Gamma^D$, and $\theta_D$ such that

$$
\int_{\Omega} v_{\mathcal{M}}^{\infty} \Phi_1 (\xi_{\mathcal{M}}^2) \le C_{\mathrm{LS}, \Gamma^D, \infty}^2 \left( 1 + \log \left( 1 + \| \xi_{\mathcal{M}} - 1 \|_{L^2(\Omega)}^2 \right) \right) |\underline{\xi}_D|^2 |\underline{1}_{\mathrm{D}}|^2 \quad (1.87)
$$

1.B Nonlinear scheme for mixed Dirichlet-Neumann bound-
ary conditions

In this appendix, we introduce and analyse a version of the nonlinear scheme for the evolution
problem (1.1) when $|\Gamma^D| > 0$. In order to perform the asymptotic analysis, we need to assume
that the data are compatible with the thermal equilibrium:

$$
f = 0, \quad g^N = 0, \quad \text{and there exists } \rho^D > 0 \text{ such that } g^D = \rho^D e^{-\phi} = \rho^D \omega.
$$

For such data, given $u^{in} \ge 0$, the solution $u$ to (1.1) is positive for $t > 0$, and converges towards $u^{\infty} = \rho^D e^{-\phi}$ when $t \to \infty$.

1.B.1 Scheme and well-posedness

Accordingly to this setting, we define $u_D^\infty = \rho^D \underline{\omega}_D$. One has $u_b^\infty \underline{1}_D \le u_D^\infty \le u_\#^\infty \underline{1}_D$, where
$0 < u_b^\infty \le u_\#^\infty$ only depend on $\rho^D$, $\phi$, and $\Omega$. Remind that, as in (1.47), given a positive $u_D \in \underline{V}_D$,
one defines $w_D(u_D) \in \underline{V}_D$ as

$$
w_K = \log\left(\frac{u_K}{u_K^\infty}\right) \quad \forall K \in \mathcal{M}, \quad w_\sigma = \log\left(\frac{u_\sigma}{u_\sigma^\infty}\right) \quad \forall \sigma \in \mathcal{E}.
$$
For mixed boundary conditions, the discrete problem reads: Find $(\underline{u}_D^n \in \underline{V}_D)_{n \ge 1}$ positive such that

$$ \left\{ \begin{array}{ll} \frac{1}{\Delta t}(u_M^n - u_M^{n-1}, v_M)_\Omega + T_D(\underline{u}_D^n, \underline{w}_D(\underline{u}_D^n), \underline{v}_D) = 0 & \forall \underline{v}_D \in \underline{V}_{D,0}^D, \quad (1.88a) \\ \underline{w}_D(\underline{u}_D^n) \in \underline{V}_{D,0}^D, & (1.88b) \end{array} \right. $$

$$ \underline{w}_D(\underline{u}_D^n) \in \underline{V}_{D,0}^D, \quad (1.88b) $$

$$ u_K^0 = \frac{1}{|K|} \int_K u^{in} \quad \forall K \in \mathcal{M}. \quad (1.88c) $$

Notice that, since for all $\sigma \in \mathcal{E}$, $w_\sigma^n = \log(\frac{u_\sigma^n}{u_\sigma^\infty})$, the equation (1.88b) only means that, for all $\sigma \in \mathcal{E}_{ext}^D$, $u_\sigma^n = u_\sigma^\infty$, which enforces strongly the Dirichlet boundary condition on $\Gamma^D$. One can show the following existence result.

**Theorem 5** (Existence of positive solutions and entropy/entropy dissipation relation). Let $u^{in} \in L^2(\Omega)$ be a non-negative function. There exists at least one positive solution $(\underline{u}_D^n \in \underline{V}_D)_{n \ge 1}$ to the nonlinear scheme (1.88). It satisfies the following entropy/dissipation relation:

$$ \forall n \in \mathbb{N}, \quad \frac{\mathbb{E}^{n+1} - \mathbb{E}^n}{\Delta t} + \mathbb{D}^{n+1} \le 0, \quad (1.89) $$

where $\mathbb{E}^n$ and $\mathbb{D}^n$ are, respectively, the discrete relative entropy and dissipation defined in (1.49). Moreover, there exists $\varepsilon > 0$, depending on $\Lambda$, $\phi$, $u^{in}$, $\rho^D$, $\Omega$, $d$, $\Delta t$, and $\mathcal{D}$ such that, for all $n \ge 1$, $u_K^n \ge \varepsilon$ for all $K \in \mathcal{M}$ and $u_\sigma^n \ge \varepsilon$ for all $\sigma \in \mathcal{E}$.

The proof of this theorem relies on the same arguments as the one of Theorem 1 for (homo- geneous) pure Neumann boundary conditions. The major difference lies in the counterpart of Lemma 2, which is no longer based on the positivity of the mass, but on the prescribed (zero) value on the Dirichlet faces.

### 1.B.2 Long-time behaviour

In the next theorem, we state the long-time behaviour of the discrete solutions to the nonlinear scheme (1.88).

**Theorem 6** (Asymptotic stability). If $(\underline{u}_D^n \in \underline{V}_D)_{n \ge 1}$ is a (positive) solution to (1.88), then the discrete entropy decays exponentially fast in time: there is $v_{nl,\Gamma^D} > 0$, depending on $\Lambda$, $\phi$, $\Gamma^D$, $\rho^D$, $u^{in}$, $\Omega$, $d$, and $\theta_D$ such that

$$ \forall n \in \mathbb{N}, \quad \mathbb{E}^{n+1} \le (1 + v_{nl,\Gamma^D} \Delta t)^{-1} \mathbb{E}^n. \quad (1.90) $$

Consequently, the discrete solution converges exponentially fast in time towards its associated discrete steady-state.

*Proof.* Let $n \in \mathbb{N}^*$. As in Section 1.4.3, one has $\mathbb{D}^n \ge \frac{1}{C_F} \hat{\mathbb{D}}^n \ge \frac{4u_b^\infty \lambda_b \alpha_b}{C_F} |\xi_D^n|_{1,D}^2$, where $C_F > 0$ depends on the data. Using the discrete log-Sobolev inequality (1.87) from Proposition 10, we
get

$$
\begin{align}
\mathbb{E}^n \le C_{LS,\Gamma^D,\infty} (1 + \log(1 + \|\xi_M^n - 1\|_{L^2_\mu(\Omega)}^2)) |\xi_D^n|_{1,D}^2 & \nonumber \\
& \le \frac{C_{LS,\Gamma^D,\infty} C_F}{4u_b^\infty \lambda_b \alpha_b} (1 + \log(1 + \|\xi_M^n - 1\|_{L^2_\mu(\Omega)}^2)) \mathbb{D}^n. \tag{1.91}
\end{align}
$$

Then, there is $C > 0$ such that (recall that $\underline{\xi}_D^n$ is positive)

$$
\|\xi_M^n - 1\|_{L^2_\mu(\Omega)}^2 \le \|(\xi_M^n)^2\|_{L^1_\mu(\Omega)} + 1 \le \|\Phi_1((\xi_M^n)^2)\|_{L^1_\mu(\Omega)} + C = (M^\infty)^{-1}\mathbb{E}^n + C,
$$

where the last inequality is an application of the Fenchel-Young inequality $x \le \Phi_1(x) + \Phi_1^*(1)$,
where $\Phi_1^*$ is the convex conjugate of $\Phi_1$ and $x = (\xi_M^n)^2$. But, since the entropy/dissipation relation
(1.89) holds, the discrete entropy decays and $\mathbb{E}^n \le \mathbb{E}^0$. Therefore, one has

$$
\|\xi_M^n - 1\|_{L^2_\mu(\Omega)}^2 \le (M^\infty)^{-1} \mathbb{E}^0 + C.
$$

Combining this estimate with (1.91), we deduce that there exists $v_{nl,\Gamma^D} > 0$, depending on $\Lambda$, $\phi$, $\Gamma^D$, $\rho^D$, $u^{in}$, $\Omega$, $d$, and $\theta_D$ such that $\mathbb{E}^n \le v_{nl,\Gamma^D} \mathbb{D}^n$. Then, using the entropy/dissipation relation (1.89), we get (1.90). $\square$

1.C Proofs of technical results

1.C.1 Discrete boundedness by mass and dissipation

We prove Lemma 2 from Section 1.3.3. To ease the reading, we first recall the result.

**Lemma 2.** Let $w_D \in V_D$, and assume that there exist $C_\# > 0$, and $M_\# \ge M_b > 0$ such that

$$
M_b \le \sum_{K \in \mathcal{M}} |K| u_K^\infty e^{w_K} \le M_\# \quad \text{and} \quad \mathbb{D}(w_D) \le C_\# \tag{1.53}
$$

Then, there exists $C > 0$, depending on $\Lambda$, $u_b^\infty$, $u_\#^\infty$, $M_b$, $M_\#$, $C_\#$, $\Omega$, $d$, and $\mathcal{D}$ such that

$$
|w_K| \le C \quad \forall K \in \mathcal{M} \quad \text{and} \quad |w_\sigma| \le C \quad \forall \sigma \in \mathcal{E}.
$$

*Proof.* For $K \in \mathcal{M}$, using (1.12) and (1.13), we first infer that

$$
\delta_K \underline{w}_K \cdot \mathbb{A}_K \delta_K \underline{w}_K = a_K^\Lambda (\underline{w}_K, \underline{w}_K) \ge \lambda_b \alpha_b |\underline{w}_K|_{1,K}^2 = \lambda_b \alpha_b \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (\underline{w}_K - \underline{w}_\sigma)^2.
$$

By definition (1.3) of the regularity parameter $\theta_D$, we have that $|\frac{\sigma|}{d_{K,\sigma}}| \ge \frac{h_K^{d-2}}{\theta_D}$ for all $\sigma \in \mathcal{E}_K$, so that

$$
\delta_K \underline{w}_K \cdot \mathbb{A}_K \delta_K \underline{w}_K \geq \frac{\lambda_b \alpha_b}{\theta_D} h_K^{d-2} |\delta_K \underline{w}_K|^2. \quad (1.92)
$$
By the expression (1.52) of $\mathbb{D}(\underline{w}_D)$, and the local lower bound (1.92), we thus get

$$
\begin{align*}
\mathbb{D}(\underline{w}_D) &= \sum_{K \in \mathcal{M}} r_K(\underline{u}_K^\infty \times \exp(\underline{w}_K)) \delta_K \underline{w}_K \cdot \mathbb{A}_K \delta_K \underline{w}_K \\
&\geq \frac{\lambda_b \alpha_b}{\theta_D} \sum_{K \in \mathcal{M}} h_K^{d-2} r_K(\underline{u}_K^\infty \times \exp(\underline{w}_K)) |\delta_K \underline{w}_K|^2 \\
&= \frac{\lambda_b \alpha_b}{\theta_D} \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{E}_K} h_K^{d-2} r_K(\underline{u}_K^\infty \times \exp(\underline{w}_K))(w_K - w_\sigma)^2.
\end{align*}
$$

Let $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}_K$ be fixed. Using, successively, the definition (1.41) of $r_K$ combined with the definition (1.43) of $f_{\mathcal{E}_K}$, the combination of (1.48) with assumptions (1.42a) and (1.42c), and the assumptions (1.42b) and (1.42d) combined with the bound (1.4) on $|\mathcal{E}_K|$, we infer, for $w_\sigma \neq w_K$,

$$
\begin{align*}
r_K(\underline{u}_K^\infty \times \exp(\underline{w}_K))(w_K - w_\sigma)^2 &\ge \frac{1}{|\mathcal{E}_K|} m(u_K^\infty e^{w_K}, u_\sigma^\infty e^{w_\sigma}) (w_K - w_\sigma)^2 \\
&\ge \frac{u_b^\infty}{|\mathcal{E}_K|} m(e^{w_K}, e^{w_\sigma}) (w_K - w_\sigma)^2 \\
&\ge \frac{u_b^\infty}{d\theta_D^2} (e^{w_K} - e^{w_\sigma}) (w_K - w_\sigma) \ge 0,
\end{align*}
$$

and we verify that this inequality still holds when $w_\sigma = w_K$. Since $\mathbb{D}(\underline{w}_D) \le C_\#$ by (1.53), for all $K \in \mathcal{M}$, and all $\sigma \in \mathcal{E}_K$, we have

$$
0 \le (e^{w_K} - e^{w_\sigma})(w_K - w_\sigma) \le \zeta h_K^{2-d}, \quad (1.93)
$$

with $\zeta = \frac{dC_\#\theta_D^3}{\lambda_b \alpha_b u_b^\infty} > 0$ (recall that $\alpha_b$ depends on $\Omega$, $d$, and $\theta_D$). Besides, since $\sum_{K \in \mathcal{M}} |K| u_K^\infty e^{w_K} \le M_\#$
again by (1.53), we have $|K| u_K^\infty e^{w_K} \le M_\#$ for all $K \in \mathcal{M}$. Similarly, since $\sum_{K \in \mathcal{M}} |K| u_K^\infty e^{w_K} \ge M_\#$,
there exists $K_0 \in \mathcal{M}$ such that $|\Omega| u_{K_0}^\infty e^{w_{K_0}} \ge M_\# > 0$. Combining these bounds, we infer that there
exists $K_0 \in \mathcal{M}$ such that

$$
\log\left(\frac{M_b}{|\Omega|u_{\#}^{\infty}}\right) \le w_{K_0} \le \log\left(\frac{M_{\#}}{|K_0|u_{b}^{\infty}}\right). \quad (1.94)
$$

Now, let us show that we can similarly frame all the other components of $\underline{w}_D$.

For $a, x \in \mathbb{R}$, let us define $E(a, x) = (e^x - e^a)(x-a) \ge 0$. Observe that $E(a, y+a)e^{-a} = (e^y-1)y =:$
$\xi(y)$ and that $\xi$ is continuous, strictly decreasing for $y < 0$, strictly increasing for $y > 0$, $\xi(0) = 0$,
and $\xi(y) \to +\infty$ when $y \to \pm\infty$. Let $b, a_\# > 0$, and take $|a| \le a_\#$. By the properties of $\xi$, if $E(a, x) \le b$,
then $|x| \le \kappa_b(a_\#) := a_\# + \max\{|y| \text{ s.t. } \xi(y) = b e^{a_\#}\}$. We can thus infer that if $(x_k)_{k=0,...,m}$ is a finite
sequence of real numbers such that $E(x_k, x_{k+1}) \le b$ and $|x_0| \le a_\#$, then $|x_m| \le \kappa_b^{(m)}(a_\#)$ where $\kappa_b^{(m)}$ is
$m$ compositions of $\kappa_b$. In particular, the bound only depends on $a_\#$, $m$ and $b$.

Now we can conclude the proof. Because of the connectivity of the mesh, for any cell
K (respectively, face σ) there is a finite sequence of components of w_D, denoted (x_k)_{k=0,...,m},
starting at x_0 = w_{K_0} and finishing at x_m = w_K (respectively, x_m = w_σ) such that, by (1.93),
E(x_k, x_{k+1}) ≤ b := ζh_D^{2-d}. The inequality (1.94) yields the initial bound on |x_0|, and one concludes
by the above argument.
□
1.C.2 A local comparison result

We prove a local comparison result between the matrices A_K and some (local) diagonal matrices. The proof relies on arguments that are similar to those advocated in [56] to analyse the VAG scheme.

**Lemma 4.** For $K \in \mathcal{M}$, let $A_K \in \mathbb{R}^{|\mathcal{E}_K| \times |\mathcal{E}_K|}$ be the matrix defined by (1.10). The matrices $A_K$ are symmetric positive-definite, and there exists $C_A > 0$, only depending on $\Lambda, \Omega, d$, and $\theta_D$ such that

$$ \forall K \in \mathcal{M}, \text{ Cond}_2(A_K) \leq C_A, $$

where $\text{Cond}_2(A_K) = \|A_K^{-1}\|_2 \|A_K\|_2$ is the condition number of the matrix $A_K$. Moreover, letting for $K \in \mathcal{M}$, $B_K \in \mathbb{R}^{|\mathcal{E}_K| \times |\mathcal{E}_K|}$ be the diagonal matrix with entries

$$ B_K^{\sigma\sigma} = \sum_{\sigma' \in \mathcal{E}_K} |A_K^{\sigma\sigma'}| \quad \text{for all } \sigma \in \mathcal{E}_K, \qquad (1.95) $$

there exists $C_B > 0$, only depending on $\Lambda, \Omega, d$, and $\theta_D$ such that

$$ \forall K \in \mathcal{M}, \forall w \in \mathbb{R}^{|\mathcal{E}_K|}, \quad w \cdot A_K w \le w \cdot B_K w \le C_B w \cdot A_K w. \qquad (1.96) $$

*Proof.* Let $K \in \mathcal{M}$ and $k = |\mathcal{E}_K|$. As a direct consequence of its definition (1.10), the matrix $A_K \in \mathbb{R}^{k \times k}$ is symmetric and positive semi-definite. Now, let $w = (w_\sigma)_{\sigma \in \mathcal{E}_K} \in \mathbb{R}^k$, and define $\underline{v}_K \in \underline{V}_K$ such that

$$ v_K = 0 \quad \text{and} \quad v_\sigma = -w_\sigma \quad \text{for all } \sigma \in \mathcal{E}_K. $$

Then, $\delta_K \underline{v}_K = (v_K - v_\sigma)_{\sigma \in \mathcal{E}_K} = w$. By (1.92), we immediately get that

$$ w \cdot A_K w \ge \frac{\lambda_b \alpha_b}{\theta_D} h_K^{d-2} |w|^2, $$

which implies, since $w \in \mathbb{R}^k$ is arbitrary, that $A_K$ is invertible, and gives us a lower bound on its smallest eigenvalue. By the same arguments advocated to prove (1.92), noticing that $\frac{|\sigma|}{d_{K,\sigma}} \le \theta_D h_K^{d-2}$ for all $\sigma \in \mathcal{E}_K$, we infer that

$$ w \cdot A_K w \le \lambda_\# \alpha_\# \theta_D h_K^{d-2} |w|^2. $$

We eventually get, using the estimates on the eigenvalues of $A_K$, that

$$ \text{Cond}_2(A_K) \le \frac{\lambda_\# \alpha_\#}{\lambda_b \alpha_b} \theta_D^2 = C_A, \qquad (1.97) $$

with $C_A > 0$ only depending on $\Lambda, \Omega, d$, and $\theta_D$. Now, by (1.95), since $A_K$ is symmetric, we have

$$ w \cdot B_K w = \sum_{\sigma \in \mathcal{E}_K} \sum_{\sigma' \in \mathcal{E}_K} |A_K^{\sigma\sigma'}| w_\sigma^2 = \sum_{\sigma \in \mathcal{E}_K} \sum_{\sigma' \in \mathcal{E}_K} |A_K^{\sigma\sigma'}| w_{\sigma'}^2, $$

and we can use the half-sum to get

$$ w \cdot B_K w = \sum_{\sigma \in \mathcal{E}_K} \sum_{\sigma' \in \mathcal{E}_K} |A_K^{\sigma\sigma'}| \frac{w_\sigma^2 + w_{\sigma'}^2}{2}. $$
Using Young's inequality, we infer

$$ 
\begin{aligned}
w \cdot A_K w &= \sum_{\sigma \in \mathcal{E}_K} \sum_{\sigma' \in \mathcal{E}_K} A_K^{\sigma\sigma'} w_\sigma w_{\sigma'} \le \sum_{\sigma \in \mathcal{E}_K} \sum_{\sigma' \in \mathcal{E}_K} |A_K^{\sigma\sigma'} ||w_\sigma|| w_{\sigma'}| \\
&\le \sum_{\sigma \in \mathcal{E}_K} \sum_{\sigma' \in \mathcal{E}_K} |A_K^{\sigma\sigma'}| \frac{w_\sigma^2 + w_{\sigma'}^2}{2} = w \cdot B_K w.
\end{aligned}
$$

For the second inequality, by symmetry of $A_K$, we have

$$ w \cdot B_K w = \sum_{\sigma \in \mathcal{E}_K} B_K^{\sigma\sigma} w_\sigma^2 \le \max_{\sigma \in \mathcal{E}_K} (B_K^{\sigma\sigma}) \sum_{\sigma' \in \mathcal{E}_K} w_\sigma^2 = \max_{\sigma \in \mathcal{E}_K} \left( \sum_{\sigma' \in \mathcal{E}_K} |A_K^{\sigma'\sigma}| \right) |w|^2 = \|A_K\|_1 |w|^2. $$

The space $\mathbb{R}^{k \times k}$ being of finite dimension, the norms $\| \cdot \|_1$ and $\| \cdot \|_2$ are equivalent, and there exists $\gamma_k > 0$ such that $\| \cdot \|_1 \le \gamma_k \| \cdot \|_2$. Moreover, since $A_K$ is symmetric positive-definite, the following inequality holds:

$$ w \cdot A_K w \ge \frac{\|A_K\|_2}{\mathrm{Cond}_2(A_K)} |w|^2. $$

From the previous estimates and (1.97), we deduce that

$$ w \cdot B_K w \le \gamma_k \mathrm{Cond}_2(A_K) w \cdot A_K w \le \gamma_k C_A w \cdot A_K w. $$

But, according to (1.4), we have $\max_{K \in \mathcal{M}} \gamma_k \le \max_{(d+1) \le l \le d \theta_D^2} \gamma_l$, therefore

$$ w \cdot B_K w \le C_B w \cdot A_K w, $$

where $C_B = C_A \max_{(d+1) \le l \le d \theta_D^2} \gamma_l$ is a positive constant only depending on $\Lambda$, $\Omega$, $d$, and $\theta_D$. This completes the proof of the comparison result (1.96). $\square$
Chapter 2

A structure preserving hybrid finite
volume scheme for semiconductor
models with magnetic field on
general meshes

## Outline of the current chapter

<table><tr><td>2.1</td><td>Introduction</td><td>78</td></tr><tr><td>2.2</td><td>Discrete setting and schemes</td><td>83</td></tr><tr><td>2.2.1</td><td>Mesh, discrete unknowns and boundary data</td><td>83</td></tr><tr><td>2.2.2</td><td>Foundations of the hybrid finite volume method</td><td>86</td></tr><tr><td>2.2.3</td><td>Description of the schemes and main theorem</td><td>88</td></tr><tr><td>2.3</td><td>Analysis of the stationary scheme</td><td>91</td></tr><tr><td>2.3.1</td><td>Well-posedness</td><td>91</td></tr><tr><td>2.3.2</td><td>Correspondence between discrete densities and discrete potentials</td><td>92</td></tr><tr><td>2.4</td><td>Analysis of the transient scheme</td><td>94</td></tr><tr><td>2.4.1</td><td>Discrete entropy structure</td><td>94</td></tr><tr><td>2.4.2</td><td>Existence of solutions</td><td>97</td></tr><tr><td>2.4.3</td><td>Long-time behaviour</td><td>99</td></tr><tr><td>2.5</td><td>Numerical results</td><td>101</td></tr><tr><td>2.5.1</td><td>Implementation</td><td>102</td></tr><tr><td>2.5.2</td><td>Proof of concept</td><td>105</td></tr><tr><td>2.5.3</td><td>Long-time behaviour of the discrete solutions</td><td>110</td></tr><tr><td>2.6</td><td>Conclusion</td><td>112</td></tr><tr><td>2.A</td><td>Discrete boundedness by entropy and dissipation</td><td>113</td></tr></table>

This chapter content is the article [213] published in ESAIM: Mathematical Modelling and Numerical Analysis.
We are interested in the discretisation of a drift-diffusion system in the frame-
work of hybrid finite volume (HFV) methods on general polygonal/polyhedral
meshes. The system under study is composed of two anisotropic and nonlinear
convection-diffusion equations with nonsymmetric tensors, coupled with a Pois-
son equation and describes in particular semiconductor devices immersed in a
magnetic field. We introduce a new scheme based on an entropy-dissipation re-
lation and prove that the scheme admits solutions with values in admissible sets
- especially, the computed densities remain positive. Moreover, we show that the
discrete solutions to the scheme converge exponentially fast in time towards the
associated discrete thermal equilibrium. Several numerical tests confirm our theo-
retical results. Up to our knowledge, this scheme is the first one able to discretise
anisotropic drift-diffusion systems while preserving the bounds on the densities.

2.1 Introduction

We are interested in the numerical approximation of a generalised anisotropic Van Roos-
broeck's drift-diffusion system, inspired by the model introduced in [143] by Gajewski and
Gärtner to describe semiconductors devices immersed in exterior magnetic fields. Such a
model differs from the classical drift-diffusion system from its anisotropic nature, induced by
the magnetic field. From a numerical point of view, this anisotropy leads to several difficul-
ties, including the preservation of bounds on the solution (see the numerical results of [143]).
The main originality of our work is to address this difficulty using a nonlinear Hybrid Finite
Volume method [123, 70], specifically designed to handle anisotropy while preserving the
bounds. Hence, up to our knowledge, the scheme introduced here is the first bound-preserving
scheme for anisotropic drift-diffusion systems. From a general perspective, the Van Roosbroeck's
drift-diffusion system, initially introduced in [258], is one of the fundamental models for the
description and the simulation of semiconductor devices. It is a macroscopic model taking into
account the densities of the charge carriers (the electrons, negatively charged, and the holes, of
positive charge) alongside with an electrostatic potential induced by the spatial inhomogeneity of
the charges. Mathematically speaking, the system consists of two parabolic convection-diffusion
equations coupled with a Poisson equation on the electrostatic potential. Among the other
different generalisations of the initial model proposed in [258], one can mention the incorpo-
ration of nonlinear terms of recombination-generation (see for example [208]), as well as the
use of non-Boltzmann statistics [144, 125] to describe some specific physical situations (high
carrier densities) or devices (for example, organic semiconductors [257]), which leads to the
modification of the diffusive term into a nonlinear one. Another generalisation of the system, as
discussed above, is presented in [143]: the authors consider that the semiconductor is immersed
into an exterior magnetic field. Such a situation leads to the introduction of anisotropy tensors
(which are nonsymmetric) in the convection-diffusion equations, related to the magnetic field.

In the present work, we consider a general model that encompasses the different features
mentioned above. More precisely, let $\Omega$ be an open, bounded, connected polytopal subsets
of $\mathbb{R}^d$, with $d \in \{1, 2, 3\}$ whose boundary $\partial\Omega$ is divided into two disjoint subset $\partial\Omega = \bar{\Gamma}^D \cup \bar{\Gamma}^N$,
with $|\bar{\Gamma}^D| > 0$. We are interested in the following problem, where the unknowns $N, P$ and $\phi$ are
functions from $\mathbb{R}_+ \times \Omega$ to $\mathbb{R}$:

$$ egin{cases} 
\partial_t N - \operatorname{div}(N \Lambda_N \nabla(h(N) - \phi)) = -R(N, P) & \text{in } \mathbb{R}_+ \times \Omega \\ 
\partial_t P - \operatorname{div}(P \Lambda_P \nabla(h(P) + \phi)) = -R(N, P) & \text{in } \mathbb{R}_+ \times \Omega \\ 
-\operatorname{div}(\Lambda_\phi \nabla\phi) = C + P - N & \text{in } \mathbb{R}_+ \times \Omega \\ 
N = N^D, P = P^D \text{ and } \phi = \phi^D & \text{on } \mathbb{R}_+ \times \Gamma^D \\ 
N \Lambda_N \nabla(h(N) - \phi) \cdot n = P \Lambda_P \nabla(h(P) + \phi) \cdot n = \Lambda_\phi \nabla\phi \cdot n = 0 & \text{on } \mathbb{R}_+ \times \Gamma^N \\ 
N(0, \cdot) = N^{in} \text{ and } P(0, \cdot) = P^{in} & \text{in } \Omega, 
\end{cases} 
\quad (2.1) $$

where $n$ denotes the unit normal vector to $\partial\Omega$ pointing outward $\Omega$.

Let us give some insight and explanations about this system, and make precise assumptions on the data. First, the unknowns $N$ and $P$ refer respectively to the densities of electrons and holes, and $\phi$ refers to the electrostatic potential. The densities $N$ and $P$ take values in the set of admissible densities $I_h = ]0, a[, \text{ where } a \in ]0, +\infty[)$ is the upper bound on the density. $I_h$ is in fact the definition domain of the function $h: I_h \to \mathbb{R}$, which depends on the statistics used to describe the relation between the densities and the chemical potential. We assume that $h$ is a $C^1$ function, such that $h' > 0$ on $I_h$, with limits $\lim_{0}h = -\infty$ and $\lim_{a}h = +\infty$. We denote by $g = h^{-1}: \mathbb{R} \to I_h$ the inverse function of $h$. Note that $g$ is positive, and that $g' > 0$ on $\mathbb{R}$. Moreover, we assume that there exists $g_0 > 0$ such that for any real $s$, $g'(s) \le g_0 g(s)$. We refer the reader to [144, 125, 127] for more details about these statistics functions, but let us emphasise here three classical cases that fall under the scope of our assumptions:

(i) the Boltzmann statistics, for which $g = \exp$, $a = +\infty$ and $h = \log$;

(ii) the Blakemore statistics, for which $g(s) = (\gamma + e^{-s})^{-1}$, $a = \frac{1}{\gamma}$ and $h(s) = \log(s/(1 - \gamma s))$, with $\gamma > 0$ (a relevant physical value is $\gamma = 0.27$);

(iii) the Fermi-Dirac statistics of order $1/2$, for which $g(s) = \frac{2}{\sqrt{\pi}} \int_{0}^{+\infty} \frac{\sqrt{s}}{1+e^{-2s}} ds$ and $a = +\infty$.

The statistics (ii) can be interpreted as an approximation of (iii) in low density situations (see the seminal paper of Blakemore [31]), and similarly (ii) can be approximated by (i) if the carrier densities are small enough. For the purpose of future analysis, we define $H: x \to \int_{1}^{x} h(s) ds$ and $G: x \mapsto \int_{-\infty}^{x} g(s) ds$. Notice that $H$ and $G$ are strictly convex, and that $G$ is positive. The function $C \in L^\infty(\Omega)$ is the doping profile of the semiconductor, which characterises the device under study. In practice, $C$ is a discontinuous function.

The term $R(N, P)$ corresponds to the recombination-generation rate, which can be interpreted as a reaction term between electrons and holes. We assume that this term has the following form:

$$ R(N, P) = r(N, P) (e^{h(N)+h(P)} - 1), \quad (2.2) $$

where $r$ is a continuous non-negative function, which can be nonconstant in space. Pertinent choices of $r$ are for example (see [207, 208, 128]) the Auger recombination $r(N, P) = c_N N + c_P P$, the Shokley-Read-Hall (SRH) term $r(N, P) = (\tau_N N + \tau_P P + \tau_C)^{-1}$, or no recombination term $r = 0$. The tensor $\Lambda_\phi \in L^\infty(\Omega, \mathbb{R}^{d \times d})$ is the (rescaled) permittivity of the medium. We assume that it is a symmetric and uniformly elliptic tensor, in the sense that there exists $\lambda_\#^\phi \ge \lambda_b^\phi > 0$ such that $\forall \xi \in \mathbb{R}^d$,

$$ \lambda_b^\phi |\xi|^2 \le \xi \cdot \Lambda_\phi \xi \le \lambda_\#^\phi |\xi|^2, $$

In practical situations, the permittivity is often assumed isotropic, leading to a tensor of the
form $\Lambda_{\phi} = \lambda^2 \epsilon I_d$, where $\lambda$ is the rescaled Debye length of the system, which accounts for the nondimensionalisation, and $\epsilon : \Omega \to R_+^*$ is a uniformly positive function corresponding to the dielectric permittivity of the material. Note that relevant values of the Debye length can be very small, inducing some stiff behaviours. Moreover, since the devices are often made of different materials, the function $\epsilon$ can be non-regular, and exhibits discontinuities at junctions between different materials.

The tensors $\Lambda_N$ and $\Lambda_P$ are diffusion tensors in $L^{\infty}(\Omega, \mathbb{R}^{d \times d})$, related to the exterior magnetic field. We refer to [143] for a detailed description of the semiconductor models with magnetic field. A typical example is the case $d = 2$, where the magnetic field $\vec{B}$ is orthogonal to the device: letting $b = |\vec{B}|$, the tensors write

$$ \Lambda_N = \frac{\mu_N}{1+b^2} \begin{pmatrix} 1 & b \\ -b & 1 \end{pmatrix} \quad \text{and} \quad \Lambda_P = \frac{\mu_P}{1+b^2} \begin{pmatrix} 1 & -b \\ b & 1 \end{pmatrix}, $$

where $\mu_N : \Omega \to \mathbb{R}_+^*$ and $\mu_P : \Omega \to \mathbb{R}_+^*$ are the rescaled mobilities of the electrons and holes. For our analysis, we will consider general tensors, and only assume that the tensors are uniformly elliptic and bounded, in the sense that there exist $\lambda_{\#} \ge \lambda_b > 0$ such that for any $\xi \in \mathbb{R}^d$,

$$ \lambda_b |\xi|^2 \le \xi \cdot \Lambda_N \xi, \quad \lambda_b |\xi|^2 \le \xi \cdot \Lambda_P \xi \quad \text{and} \quad |\Lambda_N \xi| \le \lambda_{\#} |\xi|, |\Lambda_P \xi| \le \lambda_{\#} |\xi|. $$

Note that these tensors are nonsymmetric in general (if $b \ne 0$).

Concerning the boundary data, we consider mixed Dirichlet-Neumann conditions. On $\Gamma^D$, which corresponds to the ohmic contacts, we impose that the values of the densities and the electrostatic potential are equal to $N^D$, $P^D$ and $\phi^D$ which are the traces of $H^1$ functions, denoted the same way. On $\Gamma^N$, which corresponds to insulated contacts, we impose a null flux. Concerning the initial conditions, we impose the initial densities $N^{in}$ and $P^{in}$, which are in $L^{\infty}(\Omega)$. We assume that the initial and boundary data are uniformly bounded in $I_h$: there exist $0 < m < M < a$ such that, almost everywhere on $\Omega$,

$$ m \le N^{in}, P^{in}, N^D, P^D \le M. \qquad (2.3) $$

In the sequel, we will use the generic term of “data” or “physical data” to refer to the objects previously introduced, namely $h, r, \Lambda_{\phi}, \Lambda_N, \Lambda_P, C, N^D, P^D, \phi^D, N^{in}$ and $P^{in}$.

The existence and uniqueness of the global solutions to the drift-diffusion system (2.1) have been originally investigated in the case of Boltzmann statistics ($h = \log$), see for example [214, 142, 208]. Concerning models with general statistics described above, we refer to [144, 149]. One of the main properties of the solutions is the fact that the carrier densities $N$ and $P$ take values in $I_h$: in the sequel we will say that the densities are $I_h$-valued. Another natural concern lies in the long-time behaviour of the system, that is to say the asymptotics of the solutions when $t \to +\infty$: under some assumptions on the data, the solutions are shown to converge exponentially fast towards some steady-state called thermal equilibrium. The thermal equilibrium is a specific steady-state ($N^e, P^e, \phi^e$) for which the electrons and holes currents cancel. Under the hypothesis on the function $h$, according to [209], there is no vacuum in this state, which hence satisfies $\nabla(h(N^e) - \phi^e) = 0$ and $\nabla(h(P^e) + \phi^e) = 0$ on $\Omega$. Therefore, there exist some constants $\alpha_N$ and $\alpha_P$ such that $h(N^e) - \phi^e = \alpha_N$ and $h(P^e) + \phi^e = \alpha_P$ in $\Omega$ (recall that $\Omega$ is connected). In particular, this implies the following compatibility condition on the boundary data:

$$ h(N^D) - \phi^D = \alpha_N \quad \text{and} \quad h(P^D) + \phi^D = \alpha_P \quad \text{on} \quad \Gamma^D. \qquad (2.4) $$

Moreover, at the thermal equilibrium, there is no more recombination nor generation, that is to
say $R(N^e, P^e) = 0$. This implies, according to the definition of $R$, that $h(N^e) + h(P^e) = 0$ if $r \neq 0$. Combined with the previous constraint, this imposes another compatibility condition on the constants if $r \neq 0$:

$$ \alpha_N + \alpha_P = 0. \qquad (2.5) $$

Under these conditions, we get that $N^e = g(\alpha_N + \phi^e)$ and $P^e = g(\alpha_P - \phi^e)$. Thus, the electrostatic potential at the equilibrium $\phi^e$ is a solution to the following nonlinear Poisson equation (called Poisson-Boltzmann in case of Boltzmann statistics, i.e. if $g = \exp$):

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\Lambda_{\phi} \nabla \phi^e) = C + g(\alpha_P - \phi^e) - g(\alpha_N + \phi^e) & \text{in } \Omega, \\ \phi^e = \phi^D & \text{on } \Gamma^D, \\ \Lambda_{\phi} \nabla \phi^e \cdot n = 0 & \text{on } \Gamma^N. \end{array} \right. \qquad (2.6) $$

In the sequel, we will assume that the boundary data are compliant with the thermal equilibrium and satisfy the compatibility conditions (2.4) and (2.5) (if $r \neq 0$). Notice that (2.4) and (2.3) imply that $\phi^D$ is bounded on $\Gamma^D$. Hence, we assume that its lifting on $\Omega$, also denoted by $\phi^D$, is in $H^1(\Omega) \cap L^\infty(\Omega)$. The analysis of the long-time behaviour of the drift-diffusion system was carried out in [143] in the case of the Boltzmann statistics with a magnetic field, and in [144] for other statistics (without magnetic field). It relies on the use of the entropy method, which consists in using a Lyapunov functional -the entropy- decaying with time. Such a method was initially used in kinetic theory and then extended to the study of various dissipative systems, as explained in [8, 174]. Let us sketch here the main lines of this method: letting

$$ \mathbb{E}(t) = \int_{\Omega} H(N) - H(N^e) - h(N^e)(N - N^e) + \int_{\Omega} H(P) - H(P^e) - h(P^e)(P - P^e) \\ + \frac{1}{2} \int_{\Omega} \nabla(\phi - \phi^e) \cdot \Lambda_{\phi} \nabla(\phi - \phi^e), $$

$$ \mathbb{D}(t) = \int_{\Omega} N \nabla(h(N) - \phi) \cdot \Lambda_N \nabla(h(N) - \phi) + \int_{\Omega} P \nabla(h(P) + \phi) \cdot \Lambda_P \nabla(h(P) + \phi) \\ + \int_{\Omega} R(N, P)(h(N) + h(P)), $$

one has the following entropy-entropy dissipation relation:

$$ \frac{d\mathbb{E}}{dt}(t) = -\mathbb{D}(t). \qquad (2.7) $$

From the assumptions on $R$ and $h$, one deduces that the entropy $\mathbb{E}$ and the dissipation $\mathbb{D}$ are non-negative quantities. Therefore, the entropy decreases, and this implies the convergence of $(N(t), P(t), \phi(t))$ towards $(N^e, P^e, \phi^e)$ when $t \to \infty$.

Note that the system (2.1) can also be interpreted as a Poisson-Nernst-Planck (PNP) system, which describes the evolution of charged particles, typically ions in a solution. Especially, the results presented in this article can be easily generalised to systems with more than two charge carriers.

From a numerical point of view, approximating the solutions to system (2.1) proves to be a challenging task due to various factors, including the nonlinear coupling between the equations, the discontinuity of the doping profile $C$, the stiffness induced by small values of the Debye
length $\lambda$, the anisotropy in the convection-diffusion equations or the non-symmetric nature of
the tensors. Moreover, it is essential to design numerical schemes that preserve the qualitative
physical properties of the continuous model. Here, the fact that the densities take values in
$I_h$ as well as the long-time behaviour are characteristic features of the continuous system that
should be preserved at the discrete level. Another challenge lies in the variety of meshes the
scheme can handle. Indeed, semiconductor devices are subject to boundary layers phenomena,
which often require local mesh refinement to be performed: matching simplicial meshes are
strongly constrained, and not very suitable for this purpose (see [207, Chapter 5.4]). Among
the various numerical methods that have been proposed for solving drift-diffusion systems, the
seminal work of Scharfetter and Gummel [239] in one dimension with Boltzmann statistics
presents the keystone idea of numerous schemes: to use the exponential relation between
chemical potentials and densities to ensure a preservation of the steady-state as well as a good
discrete long-time behaviour. This idea was generalised for different statistics (of the form
$h(s) = s^\alpha$, which do not match the assumptions presented above and allow vacuum) in [26] by
Bessemoulin-Chatard. This scheme was proved to have a good long-time behaviour in [27, 28],
following ideas from Chainais-Hilairet and Filbet introduced in [68] to analyse the long-time
behaviour of an upwind TPFA scheme. The analysis strongly relies on the adaptation of the
entropy method at the discrete level. The schemes discussed above are essentially part of the
two-points flux approximation (TPFA) [121], which suffers from limitations: the mesh has to
satisfy some orthogonality conditions, and the diffusion tensor has to be isotropic. In order to
overcome these strong constraints, various finite volume methods have been introduced in the
framework of the Poisson equation (see [109]), using auxiliary unknowns. Regrettably, these
schemes do not enjoy monotonicity properties, and especially they can have negative solutions.
In [55, 56], Cancès and Guichard proposed a solution to compensate this limitation, introducing
nonlinear schemes based on the entropic structure of the problem. The specific question of
approximating drift-diffusion systems on general meshes was already investigated in the past.
One can cite the discrete duality finite volume (DDFV, see [157, 105] for a presentation and
analysis of these schemes applied to Poisson equations) scheme of Chainais-Hilairet [64], which
can handle general polygonal meshes with $d=2$ in the framework of Boltzmann statistics. More
recently, Su and Tang designed and analysed a scheme for Poisson–Nernst–Planck systems
on general meshes in [250]. This scheme is based on two discretisation methods: a virtual
element method [20] for the Poisson equation alongside with a positive nonlinear finite volume
method [32, 48] for the convection-diffusion ones. The scheme is shown to admit solutions with
positive densities and to have an entropic structure. However, it is restricted to the Boltzmann
statistics with pure homogeneous Neumann boundary conditions and the Debye length is
assumed to be 1. Apart from these finite volume schemes, many finite elements schemes were
also designed for this type of drift-diffusion systems. Among them, one can cite the mixed and
hybrid exponential fitting schemes [42, 43] of Brezzi, Marini and Pietra, which are proposed as
finite elements generalisations of the Scharfetter–Gummel scheme [239], for linear diffusion
on simplicial meshes. These schemes were then generalised to other statistics, especially in
dimension 2, by Jüngel and Pietra in [175]. For other references about these schemes and their
extensions, we refer the reader to [41]. Up to our knowledge, the only existing scheme which
can handle the presence of a magnetic field in the model (anisotropic tensors $\Lambda_N$ and $\Lambda_p$) is the
one introduced and analysed by Gajewski and Gärtner in [143]. It can be seen as a modified
Scharfetter–Gummel scheme, and is hence restricted to Boltzmann statistics ($h = \log$). Moreover,
the scheme is restricted to triangular meshes, and bound to some strong constraints between
the mesh geometry and the magnetic field intensity. Last, the scheme does not preserve the
positivity of the densities in the presence of a strong magnetic field.

In this article, we are concerned with the design and the analysis of a numerical scheme
for (2.1) that preserves the continuous features of the system (long-time behaviour and $I_h$-valuation of the densities). One of our main concerns is the ability of the scheme to handle the large variety of possible data: statistics function $h$, recombination-generation term $R$, physical data (doping $C$, Debye length $\lambda$, initial and boundary data, magnetic field) and discretisation data (mesh, time step). Note that among the different motivations, for the support of general meshes is the desire to make local mesh adaptation much simpler in order to capture the boundary layers. Such an approach based on a local adaptation has already been investigated in 1D [58], but should become much more complicated or even impossible to use in higher dimension on constrained meshes. Hence, the scheme presented here could be a possible way to use local mesh refinement in 2D or 3D. To devise such a scheme, we use the hybrid finite volume (HFV) method, introduced in [123] by Eymard et al. in the framework of stationary diffusion problems. This method entails several interesting features: it can handle anisotropic and heterogeneous diffusion tensors alongside with very general polytopal meshes, and it benefits from a unified 1D/2D/3D formulation. The scheme introduced here and its analysis are based on the entropy structure (2.7) of the system, following the ideas of [29, 27] (in the framework of drift-diffusion systems with TPFA schemes) and [56, 52, 70] (in the framework of advection-diffusion equations, with positivity-preserving schemes supporting anisotropy and general meshes). In particular, the nonlinearity induced by the function $h$ is discretised along the principles of the nonlinear HFV scheme introduced and analysed in [70] for advection-diffusion equations. Using this approach, we get a unified scheme for the system (2.1), robust with respect to the various data of the problem. The main result of this paper, stated in Theorem 7, is the existence of solutions (with $I_h$-valued discrete densities) to the nonlinear scheme. In Theorem 9, we establish the exponential decay of the discrete solutions towards the associated thermal equilibrium. These theoretical results are validated by various numerical simulations.

The article is organised as follows. In Section 2.2, we introduce the HFV framework (mesh, discrete unknowns, discrete operators), present the schemes for the equations (2.1) and a generalisation of (2.6) (namely, the semi-linear Poisson equation (2.30)), and state our main result. In Section 2.3, we show that the stationary scheme is well posed, and discuss an important consequence, namely the correspondence between discrete densities and discrete quasi-Fermi potentials. In Section 2.4, we analyse the scheme for the transient problem, showing that it has an entropy structure and admits solutions with $I_h$-valued densities whose long-time behaviour mimics that of the continuous solutions. In Section 2.5, we discuss the implementation of the schemes, and give some numerical evidences of our theoretical results. Finally, in Appendix 2.A, we state and prove a technical result that is instrumental to show the existence of solutions.

## 2.2 Discrete setting and schemes

The aim of this section is to recall the HFV framework for diffusive problems, and introduce the schemes for the discretisation of the drift-diffusion system (2.1) and the nonlinear Poisson equation (2.30). Here, for the purpose of analysis, we present the schemes from the finite element viewpoint. We refer the reader to Section 2.5.1 for a presentation of the schemes within the framework of finite volume methods (especially for the purpose of implementation). For a detailed presentation of the HFV method in the framework of a steady variable diffusion problem with symmetric tensor, we refer the reader to [123].

### 2.2.1 Mesh, discrete unknowns and boundary data

The definitions and notation we adopt for the discretisation are essentially the same as in [123]. A discretisation of the (open, bounded, connected) polytopal set $\Omega \subset \mathbb{R}^d$, $d \in \{1, 2, 3\}$, is
defined as a triplet $\mathcal{D} = (\mathcal{M}, \mathcal{E}, \mathcal{P})$, where:

- $\mathcal{M}$ (the mesh) is a partition of $\Omega$, i.e., a finite family of nonempty disjoint (open, connected) polytopal subsets $K$ of $\Omega$ (the mesh cells) such that (i) for all $K \in \mathcal{M}$, $|K| > 0$, and (ii) $\bar{\Omega} = \bigcup_{K \in \mathcal{M}} \bar{K}$.

- $\mathcal{E}$ (the set of faces) is a partition of the mesh skeleton $\bigcup_{K \in \mathcal{M}} \partial K$, i.e., a finite family of nonempty disjoint (open, connected) subsets $\sigma$ of $\bar{\Omega}$ (the mesh faces, or mesh edges if $d=2$) such that for all $\sigma \in \mathcal{E}$, $|\sigma| > 0$ and there exists $\mathcal{H}_{\sigma}$ affine hyperplane of $\mathbb{R}^d$ such that $\sigma \subset \mathcal{H}_{\sigma}$, and $\bigcup_{K \in \mathcal{M}} \partial K = \bigcup_{\sigma \in \mathcal{E}} \bar{\sigma}$. We assume that, for all $K \in \mathcal{M}$, there exists $\mathcal{E}_K \subset \mathcal{E}$ (the set of faces of $K$) such that $\partial K = \bigcup_{\sigma \in \mathcal{E}_K} \bar{\sigma}$. For $\sigma \in \mathcal{E}$, we let $\mathcal{M}_{\sigma} = \{K \in \mathcal{M} \mid \sigma \in \mathcal{E}_K\}$ be the set of cells whose $\sigma$ is a face. Then, for all $\sigma \in \mathcal{E}$, either $\mathcal{M}_{\sigma} = \{K\}$ for a cell $K \in \mathcal{M}$, in which case $\sigma$ is a boundary face ($\sigma \subset \partial \Omega$) and we note $\sigma \in \mathcal{E}_{\text{ext}}$, or $\mathcal{M}_{\sigma} = \{K, L\}$ for two cells $K, L \in \mathcal{M}$, in which case $\sigma$ is an interface and we note $\sigma = K|L \in \mathcal{E}_{\text{int}}$.

- $\mathcal{P}$ (the set of cell centres) is a finite family $\{x_K\}_{K \in \mathcal{M}}$ of points of $\Omega$ such that, for all $K \in \mathcal{M}$, (i) $x_K \in K$, and (ii) $K$ is star-shaped with respect to $x_K$. Moreover, we assume that the Euclidean (orthogonal) distance $d_{K,\sigma}$ between $x_K$ and the affine hyperplane $\mathcal{H}_{\sigma}$ containing $\sigma$ is positive (equivalently, the cell $K$ is strictly star-shaped with respect to $x_K$).

For a given discretisation $\mathcal{D}$, we denote by $h_D > 0$ the size of the discretisation (the meshsize), defined by $h_D = \sup_{K \in \mathcal{M}} h_K$ where, for all $K \in \mathcal{M}$, $h_K = \sup_{x,y \in \bar{K}} |x-y|$ is the diameter of the cell $K$. For all $\sigma \in \mathcal{E}$, we let $\bar{x}_{\sigma} \in \sigma$ be the barycentre of $\sigma$. Finally, for all $K \in \mathcal{M}$, and all $\sigma \in \mathcal{E}_K$, we let $n_{K,\sigma} \in \mathbb{R}^d$ be the unit normal vector to $\sigma$ pointing outward $K$, and $P_{K,\sigma}$ be the (open) pyramid of base $\sigma$ and apex $x_K$ (notice that, when $d=2$, $P_{K,\sigma}$ is always a triangle). Since $|\sigma|$ and $d_{K,\sigma}$ are positive, we have $|P_{K,\sigma}| = \frac{|\sigma|d_{K,\sigma}}{d} > 0$. We depict on Figure 2.1 an example of discretisation. Notice that the mesh cells are not assumed to be convex, neither $x_K$ is assumed to be the barycentre of $K \in \mathcal{M}$. Moreover, hanging nodes are seamlessly handled with our assumptions, so that meshes with non-conforming cells are allowed (see the orange cross in Figure 2.1; the cell $K$ therein is treated as an hexagon).

<figure><img src="image_23.png" /><figcaption>Figure 2.1 – Two-dimensional discretisation and corresponding notations.</figcaption></figure>
We consider the following measure of regularity for the discretisation (which is the same as in [70, Eq. (2.1)]):

$$ \theta_D = \max \left( \max_{K \in \mathcal{M}, \sigma \in \mathcal{E}_K} \frac{h_K}{d_{K,\sigma}}, \max_{\sigma \in \mathcal{E}, K \in \mathcal{M}_\sigma} \frac{h_K^{d-1}}{|\sigma|} \right). \quad (2.8) $$

We now introduce the set of (hybrid, cell- and face-based) discrete unknowns:

$$ \underline{V}_D = \{ \underline{v}_D = ( (v_K)_{K \in \mathcal{M}}, (v_\sigma)_{\sigma \in \mathcal{E}} ) \mid v_K \in \mathbb{R}, \forall K \in \mathcal{M}, v_\sigma \in \mathbb{R}, \forall \sigma \in \mathcal{E} \}. $$

Given a mesh cell $K \in \mathcal{M}$, we let $\underline{V}_K = \mathbb{R} \times \mathbb{R}^{|\mathcal{E}_K|}$ be the restriction of $\underline{V}_D$ to $K$, and $\underline{v}_K = (v_K, (v_\sigma)_{\sigma \in \mathcal{E}_K}) \in \underline{V}_K$ be the restriction of a generic element $\underline{v}_D \in \underline{V}_D$ to $K$. Also, for $\underline{v}_D \in \underline{V}_D$, we let $v_M : \Omega \to \mathbb{R}$ and $v_\varepsilon : \bigcup_{K \in \mathcal{M}} \partial K \to \mathbb{R}$ be the piecewise constant functions such that

$$ v_{M|K} = v_K \text{ for all } K \in \mathcal{M}, \quad \text{and} \quad v_{\varepsilon|\sigma} = v_\sigma \text{ for all } \sigma \in \mathcal{E}. $$

For further use, we let $\underline{1}_D$ (respectively $\underline{0}_D$) denote the element of $\underline{V}_D$ with all coordinates equal to 1 (respectively 0). Also, given a function $f : \mathbb{R} \to \mathbb{R}$, and with a slight abuse in notation, we denote by $f(\underline{v}_D)$ the element of $\underline{V}_D$ whose coordinates are the $(f(v_K))_{K \in \mathcal{M}}$ and the $(f(v_\sigma))_{\sigma \in \mathcal{E}}$. Given $\underline{u}_D$ and $\underline{v}_D$ two elements of $\underline{V}_D$, we say that $\underline{u}_D \le \underline{v}_D$ (respectively $<, \ge,>$) if and only if for any $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}$, $\underline{u}_K \le \underline{v}_K$ and $\underline{u}_\sigma \le \underline{v}_\sigma$ (respectively $<, \ge,>$). In particular, a vector of unknowns $\underline{v}_D \in \underline{V}_D$ is $I_h$-valued if and only if $\underline{0}_D < \underline{v}_D < a\underline{1}_D$. Note that $I_h$-valued vectors have therefore positive coordinates.

We assume that the discretisation $\mathcal{D}$ is compliant with the partition $\partial\Omega = \overline{\Gamma^D} \cup \overline{\Gamma^N}$ of the boundary of the domain, in the sense that the set $\mathcal{E}_{\text{ext}}$ can be split into two disjoint subsets $\mathcal{E}_{\text{ext}}^D = \{\sigma \in \mathcal{E}_{\text{ext}} | \sigma \subset \Gamma^D\}$ and $\mathcal{E}_{\text{ext}}^N = \{\sigma \in \mathcal{E}_{\text{ext}} | \sigma \subset \Gamma^N\}$ such that $\mathcal{E}_{\text{ext}} = \mathcal{E}_{\text{ext}}^D \cup \mathcal{E}_{\text{ext}}^N$. Notice that, since $|\Gamma^D| > 0$, one has $|\mathcal{E}_{\text{ext}}^D| \ge 1$. We define the following subspace of $\underline{V}_D$, enforcing strongly a homogeneous Dirichlet boundary condition on $\Gamma^D$:

$$ \underline{V}_{D,0} = \{ \underline{v}_D \in \underline{V}_D | \forall \sigma \in \mathcal{E}_{\text{ext}}^D, v_\sigma = 0 \}. $$

In view of the upcoming analysis, we define a discrete counterpart of the $H^1$ semi-norm. Locally to any cell $K \in \mathcal{M}$, we let, for any $\underline{v}_K \in \underline{V}_K$, $|\underline{v}_K|^2_{1,K} = \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (v_K - v_\sigma)^2$. At the global level, for any $\underline{v}_D \in \underline{V}_D$, we let

$$ |\underline{v}_D|_{1,D} = \sqrt{\sum_{K \in \mathcal{M}} |\underline{v}_K|^2_{1,K}}. $$

Notice that $|\cdot|_{1,D}$ does not define a norm on $\underline{V}_D$, but if $|\underline{v}_D|_{1,D} = 0$, then there is $c \in \mathbb{R}$ such that $\underline{v}_D = c\underline{1}_D$ ($\underline{v}_D$ is constant). Thus, $|\cdot|_{1,D}$ defines a norm on the space $\underline{V}_{D,0}$.

We finally introduce discrete hybrid (Dirichlet) boundary data. The definitions below are motivated by the need to preserve the compatibility condition (2.4) at the discrete level. First, we define the discrete liftings of the densities $\underline{N}_D^D$ and $\underline{P}_D^D$ as the elements of $\underline{V}_D$ such that for any $K \in \mathcal{M}$ and any $\sigma \in \mathcal{E}$,

$$ N_{\sigma}^{D} = g\left(\frac{1}{|\sigma|} \int_{\sigma} h(N^{D})\right) \quad \text{and} \quad N_{K}^{D} = g\left(\frac{1}{|K|} \int_{K} h(N^{D})\right), \qquad (2.9) $$

$$ P_{\sigma}^{D} = g\left(\frac{1}{|\sigma|} \int_{\sigma} h(P^{D})\right) \quad \text{and} \quad P_{K}^{D} = g\left(\frac{1}{|K|} \int_{K} h(P^{D})\right). \qquad (2.10) $$
We also define the hybrid interpolate $\underline{\phi}_D^D \in \underline{V}_D$ of the lifting $\phi^D$, defined by $\phi_K^D = \frac{1}{|K|} \int_K \phi^D$ for any $K \in \mathcal{M}$ and $\phi_\sigma^D = \frac{1}{|\sigma|} \int_\sigma \phi^D$ for all $\sigma \in \mathcal{E}$. One has (see [111, Proposition B.7]) the following stability results:

$$|\underline{\phi}_D^D|_{1,D} \le C_{l,\Gamma^D} \|\phi^D\|_{H^{1/2}(\Gamma^D)} \quad \text{and} \quad -\|\phi^D\|_{L^\infty(\Omega)} \underline{1}_D \le \underline{\phi}_D^D \le \|\phi^D\|_{L^\infty(\Omega)} \underline{1}_D, \qquad (2.11)$$

where $C_{l,\Gamma^D}$ is a positive constant only depending on $\theta_D$, $\Omega$ and $\Gamma^D$. By (2.4) and (2.9)-(2.10), the following compatibility condition, discrete counterpart of (2.4), holds:

$$\forall \sigma \in \mathcal{E}_{\text{ext}}, h(N_{\sigma}^{D}) - \phi_{\sigma}^{D} = \alpha_{N} \text{ and } h(P_{\sigma}^{D}) + \phi_{\sigma}^{D} = \alpha_{p}. \qquad (2.12)$$

**Remark 10** (Discrete liftings). The definition of the discrete boundary densities (2.9) and (2.10) is driven by the desire to obtain the discrete compatibility condition (2.12). Such definition boils down to discretise quantities carrying a physical meaning, homogeneous to the quasi-Fermi potentials. Moreover, this definition ensures that $\underline{N}_D^D$ and $\underline{P}_D^D$ are $I_h$-valued. Notice that in practice, if $N^D$ is regular enough, one can choose to approximate $\frac{1}{|\sigma|} \int_\sigma h(N^D)$ by $h(N^D(\bar{x}_\sigma))$ (respectively, $\frac{1}{|K|} \int_K h(N^D)$ by $h(N^D(x_K))$ if $x_K$ is chosen to be the barycenter of $K$), and therefore find a classical expression for the discrete liftings, namely $N_\sigma^D = N^D(\bar{x}_\sigma)$ and $N_K^D = N^D(x_K)$. The same holds true for $\underline{P}_D^D$.

### 2.2.2 Foundations of the hybrid finite volume method

The HFV method hinges on the definition of a discrete gradient operator $\nabla_D$, that maps any element $\underline{v}_D \in \underline{V}_D$ to a piecewise constant $\mathbb{R}^d$-valued function on the pyramidal submesh of $\mathcal{M}$ formed by all the $P_{K,\sigma}$'s, for $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}_K$. More precisely, for all $K \in \mathcal{M}$, and all $\sigma \in \mathcal{E}_K$,

$$\nabla_D \underline{v}_D|_K = \nabla_K \underline{v}_K \quad \text{with} \quad \nabla_K \underline{v}_K|_{P_{K,\sigma}} = \nabla_{K,\sigma} \underline{v}_K = G_K \underline{v}_K + S_{K,\sigma} \underline{v}_K \in \mathbb{R}^d,$$

where $G_K \underline{v}_K$ is the consistent part of the gradient given by

$$G_K \underline{v}_K = \frac{1}{|K|} \sum_{\sigma' \in \mathcal{E}_K} |\sigma'|(v_{\sigma'} - v_K) n_{K,\sigma'} = \frac{1}{|K|} \sum_{\sigma' \in \mathcal{E}_K} |\sigma'| v_{\sigma'} n_{K,\sigma'}$$

and $S_{K,\sigma} \underline{v}_K$ is a stabilisation, given, for some parameter $\eta > 0$, by

$$S_{K,\sigma} \underline{v}_K = \frac{\eta}{d_{K,\sigma}} (v_\sigma - v_K - G_K \underline{v}_K \cdot (\bar{x}_\sigma - x_K)) n_{K,\sigma}. \qquad (2.13)$$

Let us consider a generic tensor $\Lambda \in L^\infty(\Omega, \mathbb{R}^{d \times d})$, such that for any $\xi \in \mathbb{R}^d$, $\lambda_b|\xi|^2 \le \xi \cdot \Lambda\xi$ and $|\Lambda\xi| \le \lambda_\sharp|\xi|$. Locally to any cell $K \in \mathcal{M}$, we introduce the discrete bilinear form $a_K^\Lambda: \underline{V}_K \times \underline{V}_K \to \mathbb{R}$ such that, for all $\underline{u}_K, \underline{v}_K \in \underline{V}_K$,

$$a_K^\Lambda(\underline{u}_K, \underline{v}_K) = \sum_{\sigma \in \mathcal{E}_K} |P_{K,\sigma}| \nabla_{K,\sigma} \underline{v}_K \cdot \Lambda_{K,\sigma} \nabla_{K,\sigma} \underline{u}_K = \int_K \nabla_K \underline{v}_K \cdot \Lambda \nabla_K \underline{u}_K \qquad (2.14)$$
where we set $\Lambda_{K,\sigma} = \frac{1}{|P_{K,\sigma}|} \int_{P_{K,\sigma}} \Lambda$. At the global level, we let $a_D^\Lambda: \underline{V}_D \times \underline{V}_D \to \mathbb{R}$ be the discrete bilinear form such that, for all $\underline{u}_D, \underline{v}_D \in \underline{V}_D$,

$$a_D^\Lambda(\underline{u}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} a_K^\Lambda(\underline{u}_K, \underline{v}_K) = \int_{\Omega} \nabla_D \underline{v}_D \cdot \Lambda \nabla_D \underline{u}_D.$$

As for the continuous case, the well-posedness of HFV methods for diffusion problems relies on a coercivity argument. Let $K \in \mathcal{M}$, and reason locally. By definition (2.14) of the local discrete bilinear form $a_K^\Lambda$, and from the bounds on the diffusion coefficient, we have $\lambda_b \|\nabla_K \underline{v}_K\|_{L^2(K;\mathbb{R}^d)}^2 \le a_K^\Lambda(\underline{v}_K, \underline{v}_K) \le \lambda_\# \|\nabla_K \underline{v}_K\|_{L^2(K;\mathbb{R}^d)}^2$ for all $\underline{v}_K \in \underline{V}_K$. Furthermore, the following comparison result holds (cf. [111, Lemma 13.11]): there exist $\alpha_b, \alpha_\# > 0$ with $0 < \alpha_b \le \alpha_\# < \infty$, only depending on $\Omega, d$, and $\theta_D$ such that $\alpha_b |\underline{v}_K|_{1,K}^2 \le \|\nabla_K \underline{v}_K\|_{L^2(K;\mathbb{R}^d)}^2 \le \alpha_\# |\underline{v}_K|_{1,K}^2$ for all $\underline{v}_K \in \underline{V}_K$. Combining both estimates, we infer a local coercivity and boundedness result:

$$\forall \underline{v}_K \in \underline{V}_K, \quad \lambda_b \alpha_b |\underline{v}_K|_{1,K}^2 \le a_K^\Lambda(\underline{v}_K, \underline{v}_K) \le \lambda_\# \alpha_\# |\underline{v}_K|_{1,K}^2. \qquad (2.15)$$

Summing over $K \in \mathcal{M}$, we get the following global estimates:

$$\forall \underline{v}_D \in \underline{V}_D, \quad \lambda_b \alpha_b |\underline{v}_D|_{1,D}^2 \le a_D^\Lambda(\underline{v}_D, \underline{v}_D) \le \lambda_\# \alpha_\# |\underline{v}_D|_{1,D}^2. \qquad (2.16)$$

Last, we recall the hybrid discrete Poincaré inequality (cf. [111, Lemma B.32, p = 2]): there exists $C_{p,\Gamma^D} > 0$, only depending on $\Omega, d, \Gamma^D$, and $\theta_D$ such that

$$\forall \underline{v}_D \in \underline{V}_{D,0}, \quad \|v_M\|_{L^2(\Omega)} \le C_{p,\Gamma^D} |\underline{v}_D|_{1,D}. \qquad (2.17)$$

In the sequel, we will omit the tensor in superscript and use the following notations instead:

$$a_D^\phi = a_D^{\Lambda_\phi}, \quad a_D^P = a_D^{\Lambda_P} \text{ and } a_D^N = a_D^{\Lambda_N}. \qquad (2.18)$$

One can note that, since $\Lambda_\phi$ is symmetric, one has

$$\forall (\underline{u}_D, \underline{v}_D) \in \underline{V}_{D,0}^2, \quad |a_D^\phi(\underline{u}_D, \underline{v}_D)| \le \alpha_\# \lambda_\#^\phi |\underline{u}_D|_{1,D} |\underline{v}_D|_{1,D}. \qquad (2.19)$$

For further use, we introduce the linear form $L_D^\phi: \underline{V}_D \to \mathbb{R}$ such that

$$\forall \underline{v}_D \in \underline{V}_D, \quad L_D^\phi(\underline{v}_D) = \int_{\Omega} C v_M - a_D^\phi(\underline{\phi}_D^D, \underline{v}_D). \qquad (2.20)$$

Note that, according to (2.19), the Poincaré inequality (2.17) and the estimates (2.11), one has

$$\forall \underline{v}_D \in \underline{V}_{D,0}, \quad |L_D^\phi(\underline{v}_D)| \le \|C\|_{L^2(\Omega)} \|v_M\|_{L^2(\Omega)} + \alpha_\# \lambda_\# |\underline{\phi}_D^D|_{1,D} |\underline{v}_D|_{1,D} \le L_\#^\phi |\underline{v}_D|_{1,D}, \qquad (2.21)$$

where there exists $C_L > 0$, only depending on $\Omega, \Gamma^D, d, \Lambda_\phi$ and $\theta_D$ such that

$$L_\#^\phi \le C_L (\|C\|_{L^2(\Omega)} + \|\phi^D\|_{H^{1/2}(\Gamma^D)}).$$
### 2.2.3 Description of the schemes and main theorem

In this section, we introduce an HFV scheme for Problem (2.1). The scheme is designed so as to preserve the entropic structure presented in the introduction. In the sequel, we consider a fixed spatial discretisation $\mathcal{D}$ of $\Omega$ which satisfies the assumptions of Section 2.2.1, and a fixed time step $\Delta t > 0$ for the time discretisation. For any $n \in \mathbb{N}$, we let $t^n = n\Delta t$.

Our HFV scheme relies on the hybrid nonlinear discretisation introduced in [70] for linear advection-diffusion equations, which was in turn inspired from the VAG and DDFV schemes of [56, 52]. We refer to [70] for a detailed description of the discretisation (with $h = \log$). The common idea of these schemes is to discretise the flux in a nonlinear way: formally, given a function $u$ with values in $I_h$ and a function $\phi$, the problem $\text{div}(-u\Lambda\nabla(h(u) + \phi)) = 0$ can be expressed in term of the flux

$$J = -u\Lambda\nabla(h(u) + \phi) = -u\Lambda\nabla w,$$

where $w = h(u) + \phi$ is a quasi-Fermi potential. Therefore, locally to any cell $K \in \mathcal{M}$, we define an approximation of

$$(u, w, v) \mapsto - \int_K J \cdot \nabla v = \int_K u \Lambda \nabla w \cdot \nabla v$$

under the form

$$T_K^\Lambda(\underline{u}_K, \underline{w}_K, \underline{v}_K) = \int_K r_K(\underline{u}_K) \Lambda \nabla_K \underline{w}_K \cdot \nabla_K \underline{v}_K$$

for all $I_h$-valued $\underline{u}_K \in \underline{V}_K$ and all $(\underline{w}_K, \underline{v}_K) \in \underline{V}_K^2$, where $r_K: \underline{V}_K \to I_h$, called local reconstruction operator, maps $I_h$-valued elements of $\underline{V}_K$ to $I_h$. Since $r_K(\underline{u}_K)$ is a (positive) constant on $K$, we have

$$T_K^\Lambda(\underline{u}_K, \underline{w}_K, \underline{v}_K) = r_K(\underline{u}_K) a_K^\Lambda(\underline{w}_K, \underline{v}_K), \quad (2.22)$$

where $a_K^\Lambda$ is defined by (2.14). As already pointed out in previous works on similar nonlinear schemes [49, 52, 51, 70], the definition of the local reconstruction operator is crucial to guarantee the existence of solutions as well as the long-time behaviour of the discrete solutions. Especially, we use a reconstruction which embeds information from both the local cell and face unknowns. It is a natural generalisation of the one introduced in [70]. For $\underline{u}_K$ an $I_h$-valued element of $\underline{V}_K$, we let

$$r_K(\underline{u}_K) = \frac{1}{|\mathcal{E}_K|} \sum_{\sigma \in \mathcal{E}_K} m(u_K, u_\sigma) \quad (2.23)$$

where $m: I_h^2 \to I_h$ is a continuous (mean-)function satisfying

$$\forall (x, y) \in I_h^2, m_h(x, y) \le m(x, y) \le \max(x, y), \quad (2.24)$$

and the function $m_h$ is defined by

$$m_h(x, y) = \frac{[xh(x) - H(x)] - [yh(y) - H(y)]}{h(x) - h(y)} \quad \text{if } x \ne y \text{ and } m_h(x, x) = x. \quad (2.25)$$

This function $m_h$ was used in [131, Definition 3.3] to analyse an entropic scheme for nonlinear advection-diffusion equations. One can note that in the case of the Boltzmann statistics ($h = \log$), $m_h(x, y)$ is the classical logarithmic average of $x$ and $y$, which coincides with the assumptions of [70]. Moreover, for any $(x, y) \in \mathbb{R}^2$, one has $m_h(g(x), g(y)) = \frac{xg(x) - H(g(x)) - [yg(y) - H(g(y))]}{x-y}$. Since the
derivative of $x \mapsto xg(x) - H(g(x))$ is $g$, the following simple expression holds:

$$ \forall (x, y) \in \mathbb{R}^2, m_h(g(x), g(y)) = \frac{G(x) - G(y)}{x - y}. \qquad (2.26) $$

Note that, for all $(x, y) \in I_h^2$, one has $m_h(x, y) \le \frac{x+y}{2} \le \max(x, y)$, and each expression of the previous sequence is a mean function $m$ satisfying (2.24). Heuristically, $r_K(\underline{u}_K)$ computes an average of the unknowns attached to the cell $K$, especially it contains information about all the local face unknowns. The property (2.24) will be instrumental to prove Lemma 7 and the existence of solutions to the scheme. Finally, we let $T_D^\Lambda$ be such that, for all $I_h$-valued $\underline{u}_D \in \underline{V}_D$, and all $(\underline{w}_D, \underline{v}_D) \in \underline{V}_D^2$,

$$ T_D^\Lambda(\underline{u}_D, \underline{w}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} T_K^\Lambda(\underline{u}_K, \underline{w}_K, \underline{v}_K), \qquad (2.27) $$

where the local contributions $T_K^\Lambda$ are defined by (2.22). In the sequel, we will use the following notation

$$ T_D^N = T_D^{\Lambda N} \text{ and } T_D^P = T_D^{\Lambda P} = (\text{resp., } T_K^N = T_K^{\Lambda N} \text{ and } T_K^P = T_K^{\Lambda P}). $$

Using a backward Euler discretisation in time, and the nonlinear HFV discretisation previously described in space, we consider the following scheme for the drift-diffusion system (2.1): find $((\underline{N}_D^n, \underline{P}_D^n, \underline{\phi}_D^n))_{n \in \mathbb{N}} \in \underline{V}_D^3$ such that

$$ \left\{ \begin{array}{l} \underline{w}_D^{N,n+1} = h(\underline{N}_D^{n+1}) - \underline{\phi}_D^{n+1} - \alpha_N \underline{1}_D \in \underline{V}_{D,0}, \qquad (2.28a) \\ \underline{w}_D^{P,n+1} = h(\underline{P}_D^{n+1}) + \underline{\phi}_D^{n+1} - \alpha_P \underline{1}_D \in \underline{V}_{D,0}, \qquad (2.28b) \\ \underline{\psi}_D^{n+1} = \underline{\phi}_D^{n+1} - \underline{\phi}_D^D \in \underline{V}_{D,0}, \qquad (2.28c) \\ \begin{aligned} \int_{\Omega} \frac{\underline{N}_M^{n+1} - \underline{N}_M^n}{\Delta t} v_M + T_D^N(\underline{N}_D^{n+1}, \underline{w}_D^{N,n+1}, \underline{v}_D) &= - \int_{\Omega} R(\underline{N}_M^{n+1}, \underline{P}_M^{n+1}) v_M, & \forall \underline{v}_D \in \underline{V}_{D,0}, \quad (2.28d) \\ \int_{\Omega} \frac{\underline{P}_M^{n+1} - \underline{P}_M^n}{\Delta t} v_M + T_D^P(\underline{P}_D^{n+1}, \underline{w}_D^{P,n+1}, \underline{v}_D) &= - \int_{\Omega} R(\underline{N}_M^{n+1}, \underline{P}_M^{n+1}) v_M, & \forall \underline{v}_D \in \underline{V}_{D,0}, \quad (2.28e) \\ \underline{a}_D^\phi(\underline{\psi}_D^{n+1}, \underline{v}_D) - \underline{L}_D^\phi(\underline{v}_D) &= \int_{\Omega} (\underline{P}_M^{n+1} - \underline{N}_M^{n+1}) v_M, & \forall \underline{v}_D \in \underline{V}_{D,0}, \quad (2.28f) \\ N_K^0 &= \frac{1}{|K|} \int_K N^{in} \text{ and } P_K^0 = \frac{1}{|K|} \int_K P^{in}, & \forall K \in \mathcal{M}. \quad (2.28g) \end{aligned} \right. $$

Notice that the equations (2.28a) and (2.28b) implicitly assume that, for any $n \ge 1$, the potentials $h(\underline{N}_D^n)$ and $h(\underline{P}_D^n)$ are defined, meaning that the discrete densities $\underline{N}_D^n$ and $\underline{P}_D^n$ are $I_h$-valued. Therefore, a solution to (2.28) has $I_h$-valued discrete densities by definition, which can be expressed as functions of the discrete quasi-Fermi potentials: for any $n \ge 1$,

$$ \underline{N}_D^n = g(\underline{w}_D^{N,n} + \underline{\phi}_D^n + \alpha_N \underline{1}_D) \text{ and } \underline{P}_D^n = g(\underline{w}_D^{P,n} - \underline{\phi}_D^n + \alpha_P \underline{1}_D). \qquad (2.29) $$

**Remark 11** (Discrete potentials and boundary values). The definitions of the discrete quasi-Fermi potentials (2.28a) and (2.28b) are used to enforce the boundary conditions on the discrete unknowns. Indeed, since $\underline{w}_D^{N,n} \in \underline{V}_{D,0}$, from (2.29) we get that for any $\sigma \in \mathcal{E}_{\text{ext}}^D$, $N_\sigma^n = g(0 + \underline{\phi}_\sigma^D + \alpha_N) = N_\sigma^D$
according to the compatibility condition (2.12). A similar result holds for $\underline{P}_D^n$. In fact, by definition of the trilinear forms, one can rewrite in a more natural way the equations (2.28d) and (2.28e) noticing that

$$ T_D^N(\underline{N}_D^{n+1}, \underline{w}_D^{N,n+1}, \underline{v}_D^N) = T_D^N(\underline{N}_D^{n+1}, h(\underline{N}_D^{n+1}) - \underline{\phi}_D^{n+1}, \underline{v}_D^N), \text{ and} \\ T_D^P(\underline{P}_D^{n+1}, \underline{w}_D^{P,n+1}, \underline{v}_D^P) = T_D^P(\underline{P}_D^{n+1}, h(\underline{P}_D^{n+1}) + \underline{\phi}_D^{n+1}, \underline{v}_D^P). $$

When considering general boundary conditions (i.e., if $\alpha_N$ and $\alpha_P$ are no longer constant on $\Gamma^D$, and therefore with liftings no longer constant in $\Omega$), the scheme will be the same, except for the equations (2.28d) and (2.28e) where arguments of the trilinear forms will be

$$ T_D^N(\underline{N}_D^{n+1}, \underline{w}_D^{N,n+1} + \underline{\alpha}_D^N, \underline{v}_D^N) \text{ and } T_D^P(\underline{P}_D^{n+1}, \underline{w}_D^{P,n+1} + \underline{\alpha}_D^P, \underline{v}_D^P). $$

The analysis of such a scheme is more sophisticated than the one presented here, and shall be the subject of a future work.

Plugging the expression of the discrete carrier densities (2.29) into Equation (2.28f), we get an elliptic semi-linear equation on $\underline{\psi}_D^{n+1}$. This stationary equation gives a correspondence between densities and potentials at a given time, and will be of great importance in the analysis. At the continuous level, it corresponds to the following semi-linear Poisson equation:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\Lambda_{\phi} \nabla \phi) = C + g(z^P - \phi) - g(z^N + \phi) & \text{in } \Omega, \\ \phi = \phi^D & \text{on } \Gamma^D, \\ \Lambda_{\phi} \nabla \phi \cdot n = 0 & \text{on } \Gamma^N, \end{array} \right. \qquad (2.30) $$

where $z^P$ and $z^N$ are given functions in $L^{\infty}(\Omega)$.

**Remark 12** (Electrostatic potential at thermal equilibrium). If $z^P$ and $z^N$ are constant with respective values $\alpha_P$ and $\alpha_N$ in $\Omega$, then the solution to (2.30) is the electrostatic potential at thermal equilibrium $\phi^e$ solution to (2.6).

We introduce an HFV scheme to approximate (2.30): find $\underline{\phi}_D = \underline{\psi}_D + \underline{\phi}_D^D \in \underline{V}_D$ where $\underline{\psi}_D \in \underline{V}_{D,0}$ is such that

$$ \forall \underline{v}_D \in \underline{V}_{D,0}, \quad a_D^\phi(\underline{\psi}_D, \underline{v}_D) = L_D^\phi(\underline{v}_D) + \int_{\Omega} [g(z_M^P - \phi_M^D - \psi_M) - g(z_M^N + \phi_M^D + \psi_M)] \underline{v}_M, \quad (2.31) $$

where $z_M^P$ and $z_M^N$ are the piecewise constant functions such that, for any cell $K \in \mathcal{M}$, $z_M^P = \frac{1}{|K|} \int_K z^P$ and $z_M^N = \frac{1}{|K|} \int_K z^N$ on $K$. We will prove in Theorem 8 that this scheme is well-posed. Especially, in the light of Remark 12, one can define the discrete thermal equilibrium as a triplet $(\underline{N}_D^e, \underline{P}_D^e, \underline{\phi}_D^e) \in \underline{V}_D^3$ of discrete unknowns satisfying:

$$ \left\{ \begin{array}{ll} \underline{\phi}_D^e = \underline{\phi}_D^D + \underline{\psi}_D^e, & \text{where } \underline{\psi}_D^e \in \underline{V}_{D,0} \text{ is the solution to (2.31) with } (z^N, z^P) = (\alpha_N, \alpha_P), \\ \underline{N}_D^e = g(\alpha_N \underline{1}_D + \underline{\phi}_D^e), & \\ \underline{P}_D^e = g(\alpha_P \underline{1}_D - \underline{\phi}_D^e). & \end{array} \right. \qquad (2.32) $$

Note that, by definition, the discrete densities at thermal equilibrium $\underline{N}_D^e$ and $\underline{P}_D^e$ are $I_h$-valued. We remind the reader that alternative formulations (in the spirit of finite volume methods) of the schemes (2.28) and (2.31) are presented in Section 2.5.1.
We are now in position to state the main result of this paper, whose proof is the subject of Section 2.4.

**Theorem 7** (Existence of discrete $I_h$-valued solutions to (2.28)). Let $N^{in}$ and $P^{in}$ two positive functions in $L^{\infty}(\Omega)$ satisfying the condition (2.3). There exists at least one solution $(\underline{N}_D^n, \underline{P}_D^n, \underline{\phi}_D^n)_{n \in \mathbb{N}} \in \underline{V}_D^3$ to the coupled scheme (2.28). It satisfies the following entropy-dissipation relation:

$$ \forall n \in \mathbb{N}, \quad \frac{\underline{E}^{n+1} - \underline{E}^n}{\Delta t} + \underline{D}^{n+1} \le 0, \qquad (2.33) $$

where $\underline{E}^n$ and $\underline{D}^n$ are, respectively, the discrete relative entropy and dissipation at time $t^n = n\Delta t$, defined by

$$ \underline{E}^n = \int_{\Omega} H(N_M^n) - H(N_M^e) - h(N_M^e)(N_M^n - N_M^e) \\ + \int_{\Omega} H(P_M^n) - H(P_M^e) - h(P_M^e)(P_M^n - P_M^e) + \frac{1}{2} a_D^{\phi}(\underline{\phi}_D^n - \underline{\phi}_D^e, \underline{\phi}_D^n - \underline{\phi}_D^e), \quad (2.34) $$

$$ \text{and } \underline{D}^n = T_D^N(\underline{N}_D^n, \underline{w}_D^{N,n}, \underline{w}_D^{N,n}) + T_D^P(\underline{P}_D^n, \underline{w}_D^{P,n}, \underline{w}_D^{P,n}) + \int_{\Omega} R(N_M^n, P_M^n)(h(N_M^n) + h(P_M^n)). \quad (2.35) $$

Moreover, there exist $0 < M_b \le M_\# < a$, depending on the physical data, $\Delta t$, and $\mathcal{D}$ such that

$$ \forall n \ge 1, \quad M_b \underline{1}_D \le \underline{N}_D^n \le M_\# \underline{1}_D \text{ and } M_b \underline{1}_D \le \underline{P}_D^n \le M_\# \underline{1}_D. \quad (2.36) $$

**Remark 13** ($L^\infty$ bounds). Notice that the bounds $M_b$ and $M_\#$ on the discrete densities depend on the time step $\Delta t$ and on the spatial discretisation $\mathcal{D}$ (and therefore, on the meshsize). Thus, our result is weaker than the one proved in [28, Theorem 1] in the framework of TPFA schemes. In fact, we believe that the discrete Moser iteration process used in [28] to get uniform bounds cannot be adapted in the hybrid framework, because the method is intrinsically non-monotone.

## 2.3 Analysis of the stationary scheme

In this section, we analyse the scheme (2.31) for the approximation of the semi-linear elliptic equation (2.30). First, we show that the scheme is well-posed. Then, we use this result to show that there is a correspondence between the discrete densities and the discrete quasi-Fermi potentials, which will be very useful for the following.

### 2.3.1 Well-posedness

In the next theorem, we show that the scheme (2.31) admits a unique solution. The proof relies on an energy minimisation approach.

**Theorem 8** (Well-posedness for (2.31)). Let $z^N$ and $z^P$ be two functions in $L^\infty(\Omega)$. There exists a unique $\underline{\psi}_D \in \underline{V}_{D,0}$ solution to (2.31), and there exists a positive constant $C_{Poisson}$, only depending on $\Lambda_\phi, \Omega, \Gamma^D, d$ and $\theta_D$ such that the following stability estimate holds:

$$ |\underline{\psi}_D|^2_{1,D} \le C_{Poisson} \left[ G(\|z^P\|_{L^\infty(\Omega)} + \|z^N\|_{L^\infty(\Omega)} + \|\phi^D\|_{L^\infty(\Omega)}) + \|C\|_{L^2(\Omega)}^2 + \|\phi^D\|_{H^{1/2}(\Gamma^D)}^2 \right]. \quad (2.37) $$
Moreover, the application $(z^N, z^P) \mapsto \underline{\psi}_D$ is continuous from $L^\infty(\Omega) \times L^\infty(\Omega)$ to $\underline{V}_{D,0}$.

*Proof.* We define the discrete energy functional $J_D: \underline{V}_{D,0} \to \mathbb{R}$ by

$$J_D(\underline{u}_D) = \frac{1}{2} a_D^\phi(\underline{u}_D, \underline{u}_D) - L_D^\phi(\underline{u}_D) + \int_{\Omega} G(z_M^P - \phi_M^D - u_M) + G(z_M^N + \phi_M^D + u_M). \quad (2.38)$$

Since $G$ is $C^2$ on $\mathbb{R}$, $J_D$ is $C^2$ on $\underline{V}_{D,0}$, and one can compute its differential at a point $\underline{u}_D \in \underline{V}_{D,0}$:

$$\forall \underline{v}_D \in \underline{V}_{D,0}, \quad dJ_D(\underline{u}_D) \cdot \underline{v}_D = a_D^\phi(\underline{u}_D, \underline{v}_D) - L_D^\phi(\underline{v}_D) - \int_{\Omega} [g(z_M^P - \phi_M^D - u_M) - g(z_M^N + \phi_M^D + u_M)] v_M.$$

Hence, a solution $\underline{\psi}_D \in \underline{V}_{D,0}$ to (2.31) is an element of $\underline{V}_{D,0}$ such that $dJ_D(\underline{\psi}_D) = 0$. Now, one can notice that $J_D$ is strongly convex on $\underline{V}_{D,0}$. Indeed, since $a_D^\phi$ is coercive and $L_D^\phi$ is a continuous linear form, $\underline{u}_D \mapsto \frac{1}{2} a_D^\phi(\underline{u}_D, \underline{u}_D) - L_D^\phi(\underline{u}_D)$ is strongly convex. Moreover, $G$ is a convex function (since $G' = g$ is strictly increasing), so $\underline{u}_D \mapsto \int_{\Omega} G(z_M^P - \phi_M^D - u_M) + G(z_M^N + \phi_M^D + u_M)$ is also convex. Therefore, $J_D$ has a unique minimum $\underline{\psi}_D \in \underline{V}_{D,0}$, which is a solution to (2.31). On the other hand, by convexity of $J_D$, any critical point of $J_D$ in $\underline{V}_{D,0}$ is a minimum, therefore the solution to (2.31) is unique. The continuity of $(z_M^N, z_M^P) \mapsto \underline{\psi}_D$ follows from the regularity of $J_D$ and the implicit function theorem. To obtain (2.37), one can use the monotonicity of $G$ and (2.11) to get

$$\begin{aligned} J_D(\underline{\psi}_D) \le J_D(\underline{0}_D) &= \int_{\Omega} G(z_M^P - \phi_M^D) + G(z_M^N + \phi_M^D) \\ &\le \int_{\Omega} G(\|z^P\|_{L^\infty(\Omega)} + \|\phi^D\|_{L^\infty(\Omega)}) + G(\|z^N\|_{L^\infty(\Omega)} + \|\phi^D\|_{L^\infty(\Omega)}) \\ &\le 2|\Omega|G(\|z^N\|_{L^\infty(\Omega)} + \|z^P\|_{L^\infty(\Omega)} + \|\phi^D\|_{L^\infty(\Omega)}). \end{aligned}$$

Furthermore, using the fact that $G$ is positive alongside with (2.16) and (2.21), and using Young's inequality on $L^\phi_{\#} |\underline{\psi}|_{1,D}$, one has the following lower bound:

$$J_D(\underline{\psi}_D) \ge \frac{1}{2} a_D^\phi(\underline{\psi}_D, \underline{\psi}_D) - L_D^\phi(\underline{\psi}_D) \ge \frac{\alpha_b \lambda_b}{2} |\underline{\psi}_D|_{1,D}^2 - L_\#^\phi |\underline{\psi}_D|_{1,D} \ge \frac{\alpha_b \lambda_b}{4} |\underline{\psi}_D|_{1,D}^2 - \frac{1}{\alpha_b \lambda_b} L_\#^\phi.$$

We conclude by combining the previous estimates and using the bound on $L_\#^\phi$. $\square$

<figure><img src="image_24.png" /></figure>

Notice that this result ensures the existence and uniqueness of the discrete thermal equilibrium $(\underline{N}_D^e, \underline{P}_D^e, \underline{\phi}_D^e)$ described in equation (2.32).

**Remark 14 (Uniform bounds).** The stability result (2.37) gives a uniform bound (with respect to the meshsize $h_D$) on $\underline{\psi}_D$ in energy norm $|\cdot|_{1,D}$. There is however no clear $L^\infty$ bound (uniform in $D$) on the potential. It is rather different from some previous results in the TPFA framework (see [68]), where $L^\infty$ bounds were obtained alongside with the $H^1$ estimate.

### 2.3.2 Correspondence between discrete densities and discrete potentials

In this section, we discuss the link between discrete densities and discrete potentials.
For a given couple $(\underline{N}_D, \underline{P}_D) \in \underline{V}_D^2$ of $I_h$-valued vectors, one associates a unique discrete electrostatic potential $\underline{\phi}_D = \underline{\psi}_D + \underline{\phi}_D^D \in \underline{V}_D$ such that $\underline{\psi}_D \in \underline{V}_{D,0}$ is the solution to

$$ \forall \underline{v}_D \in \underline{V}_{D,0}, \quad a_D^\phi(\underline{\psi}_D, \underline{v}_D) = L_D^\phi(\underline{v}_D) + \int_{\Omega} (P_M - N_M) v_M. \qquad (2.39) $$

Since $a_D^\phi$ is coercive on $\underline{V}_{D,0}$ and $L_D^\phi$ is a linear form, (2.39) is well-posed and $\underline{\psi}_D$ is uniquely defined, so is $\underline{\phi}_D$. Then, one associates to $(\underline{N}_D, \underline{P}_D)$ a unique couple of discrete quasi-Fermi potentials $(\underline{w}_D^N, \underline{w}_D^P) \in \underline{V}_D^2$ defined by

$$ \underline{w}_D^N = h(\underline{N}_D) - \underline{\phi}_D - \alpha_N \underline{1}_D \quad \text{and} \quad \underline{w}_D^P = h(\underline{P}_D) + \underline{\phi}_D - \alpha_P \underline{1}_D. \qquad (2.40) $$

We call these vectors the discrete electrostatic $(\underline{\phi}_D)$ and quasi-Fermi potentials $((\underline{w}_D^N, \underline{w}_D^P))$ associated to the discrete densities $(\underline{N}_D, \underline{P}_D)$.

Note that, as in Remark 11, the discrete quasi-Fermi potentials associated to $(\underline{N}_D, \underline{P}_D)$ are in $\underline{V}_{D,0}$ if and only if the densities satisfy the following discrete boundary condition:

$$ \forall \sigma \in \mathcal{E}_{\text{ext}}^D, \quad N_\sigma = N_\sigma^D \quad \text{and} \quad P_\sigma = P_\sigma^D. \qquad (2.41) $$

Conversely, for a given couple $(\underline{w}_D^N, \underline{w}_D^P) \in \underline{V}_D^2$, one associates a unique electrostatic potential $\underline{\phi}_D = \underline{\psi}_D + \underline{\phi}_D^D \in \underline{V}_D$ such that $\underline{\psi}_D \in \underline{V}_{D,0}$ is the solution to

$$ \forall \underline{v}_D \in \underline{V}_{D,0}, \quad a_D^\phi(\underline{\psi}_D, \underline{v}_D) = L_D^\phi(\underline{v}_D) + \int_{\Omega} \left( g(w_M^P + \alpha_P - \phi_M^D - \psi_M) - g(w_M^N + \alpha_N + \phi_M^D + \psi_M) \right) v_M. \qquad (2.42) $$

Indeed, since the functions $w_M^P + \alpha_P$ and $w_M^N + \alpha_N$ are in $L^\infty(\Omega)$, according to Theorem 8 applied to $(z^N, z^P) = (w_M^P + \alpha_P, w_M^N + \alpha_N)$, there exists a unique solution to problem (2.42). Then, one associates to $(\underline{w}_D^N, \underline{w}_D^P)$ a couple of $I_h$-valued discrete carrier densities $(\underline{N}_D, \underline{P}_D) \in \underline{V}_D$ by

$$ \underline{N}_D = g(\underline{w}_D^N + \underline{\phi}_D + \alpha_N \underline{1}_D) \quad \text{and} \quad \underline{P}_D = g(\underline{w}_D^P - \underline{\phi}_D + \alpha_P \underline{1}_D). $$

Notice that with this process, according to (2.42) and the definition of $\underline{N}_D$ and $\underline{P}_D$, one has that

$$ \forall \underline{v}_D \in \underline{V}_{D,0}, \quad a_D^\phi(\underline{\psi}_D, \underline{v}_D) = L_D^\phi(\underline{v}_D) + \int_{\Omega} (P_M - N_M) v_M. $$

Therefore, it means that $(\underline{w}_D^N, \underline{w}_D^P)$ and $\underline{\phi}_D$ are the discrete potentials associated to the discrete densities $(\underline{N}_D, \underline{P}_D)$ in the sense of (2.39)-(2.40).

Moreover, $(\underline{w}_D^N, \underline{w}_D^P) \in \underline{V}_{D,0}^2$ if and only if $(\underline{N}_D, \underline{P}_D)$ satisfy the discrete boundary condition (2.41).

Thus, there is a bijective correspondence between $(I_h)$-valued discrete densities and discrete quasi-Fermi potentials (in $\underline{V}_D$).

**Remark 15** (Thermal equilibrium and quasi-Fermi potentials). From the definition (2.32) of the discrete thermal equilibrium, one can notice that $(\underline{N}_D^e, \underline{P}_D^e, \underline{\phi}_D^e)$ are the discrete densities and electrostatic potential associated to the quasi-Fermi potentials at thermal equilibrium, which are the null vectors $(\underline{w}_D^{N,e}, \underline{w}_D^{P,e}) = (\underline{0}_D, \underline{0}_D)$.
## 2.4 Analysis of the transient scheme

In this section, we perform the analysis of the scheme (2.28) for the drift-diffusion system. The first two subsections are dedicated to the proof of Theorem 7. We also state and prove a qualitative property about the discrete solutions, namely the long-time behaviour in Theorem 9, in line with the work of [27].

### 2.4.1 Discrete entropy structure

The scheme (2.28) is designed so as to mimic the continuous entropic structure presented in introduction. In this part, we are interested in the entropic property of the scheme (2.28).

First, we define a discrete counterpart of the entropy and entropy dissipation. Given $(\underline{N}_D, \underline{P}_D)$ a couple of $I_h$-valued vectors, one can define the relative discrete entropy as

$$ \mathbb{E}(\underline{N}_D, \underline{P}_D) = \int_{\Omega} H(N_M) - H(N_M^e) - h(N_M^e)(N_M - N_M^e) \\ + \int_{\Omega} H(P_M) - H(P_M^e) - h(P_M^e)(P_M - P_M^e) + \frac{1}{2} a_D^{\phi}(\underline{\phi}_D - \underline{\phi}_D^e, \underline{\phi}_D - \underline{\phi}_D^e), \quad (2.43) $$

where $(\underline{N}_D^e, \underline{P}_D^e, \underline{\phi}_D^e)$ are defined in (2.32), and $\underline{\phi}_D$ is the discrete electrostatic potential associated to $(\underline{N}_D, \underline{P}_D)$ in the sense of Section 2.3.2. Notice that, by convexity of $H$, $H(N_M) - H(N_M^e) - h(N_M^e)(N_M - N_M^e)$ and the analogous term in $P_M$ are non-negative. Therefore, the discrete entropy is non-negative. Similarly, one defines the discrete dissipation by

$$ \mathbb{D}(\underline{N}_D, \underline{P}_D) = T_D^N(\underline{N}_D, \underline{w}_D^N, \underline{w}_D^P) + T_D^P(\underline{P}_D, \underline{w}_D^P, \underline{w}_D^P) + \int_{\Omega} R(N_M, P_M)(h(N_M) + h(P_M)). \quad (2.44) $$

By definition of $T_D^N$ and $T_D^P$, the first two terms are non-negative. One can also notice that the integrand of the third term writes

$$ R(N_M, P_M)(h(N_M) + h(P_M)) = r(x, N_M, P_M)(\exp(h(N_M) + h(P_M)) - 1)(h(N_M) + h(P_M)), $$

so since $r$ is non-negative and $x(e^x - 1) \ge 0$ for any real $x$, the discrete dissipation is also non-negative. Note that, given a couple $(\underline{w}_D^N, \underline{w}_D^P) \in \underline{V}_D^2$ of discrete quasi-Fermi potentials, one can also define the corresponding relative entropy and dissipation by $\mathbb{E}(\underline{w}_D^N, \underline{w}_D^P) = \mathbb{E}(\underline{N}_D, \underline{P}_D)$ and $\mathbb{D}(\underline{w}_D^N, \underline{w}_D^P) = \mathbb{D}(\underline{N}_D, \underline{P}_D)$, where $(\underline{N}_D, \underline{P}_D)$ are the discrete densities associated to $(\underline{w}_D^N, \underline{w}_D^P)$.

At the continuous level, the analysis of the drift-diffusion system relies on the entropic structure. At the discrete level, we use analogous results, based on the following lemma, which will be used in the proofs of Theorems 7 and 9.

**Lemma 5** (Discrete differentiation of the entropy). Let $\underline{N}_D^n, \underline{P}_D^n, \underline{N}_D$ and $\underline{P}_D$ be four $I_h$-valued elements of $\underline{V}_D$, and $(\underline{w}_D^N, \underline{w}_D^P) \in \underline{V}_D^2$ be the discrete quasi-Fermi potentials associated to the discrete densities $(\underline{N}_D, \underline{P}_D)$. Then, the following inequality holds:

$$ \mathbb{E}(\underline{N}_D, \underline{P}_D) - \mathbb{E}(\underline{N}_D^n, \underline{P}_D^n) \le \int_{\Omega} (N_M - N_M^n) \underline{w}_M^N + \int_{\Omega} (P_M - P_M^n) \underline{w}_M^P. \quad (2.45) $$
*Proof.* By definition of the discrete entropy, one has

$$ \mathbb{E}(\underline{N}_D, \underline{P}_D) - \mathbb{E}(\underline{N}_D^n, \underline{P}_D^n) = \int_{\Omega} H(N_M) - H(N_M^n) - h(N_M^e) (N_M - N_M^n) \\ + \int_{\Omega} H(P_M) - H(P_M^n) - h(P_M^e) (P_M - P_M^n) \\ + \frac{1}{2} (a_D^{\phi}(\underline{\phi}_D - \underline{\phi}_D^e, \underline{\phi}_D - \underline{\phi}_D^e) - a_D^{\phi}(\underline{\phi}_D^n - \underline{\phi}_D^e, \underline{\phi}_D^n - \underline{\phi}_D^e)). \quad (2.46) $$

Let us first consider the first term (in N) of (2.46): by convexity of H and definition of $\underline{N}_D^e$, we get

$$ \int_{\Omega} H(N_M) - H(N_M^n) - h(N_M^e) (N_M - N_M^n) \leq \int_{\Omega} (h(N_M) - h(N_M^e)) (N_M - N_M^n) \\ = \int_{\Omega} (h(N_M) - \phi_M^e - \alpha_N) (N_M - N_M^n). $$

For the second term of (2.46), one has an analogous result, namely

$$ \int_{\Omega} H(P_M) - H(P_M^n) - h(P_M^e) (P_M - P_M^n) \leq \int_{\Omega} (h(P_M) + \phi_M^e - \alpha_P) (P_M - P_M^n). $$

Now, notice that since $a_D^\phi$ is a symmetric and positive semi-definite bilinear form on $\underline{V}_D$, the quadratic form $\underline{u}_D \mapsto a_D^\phi(\underline{u}_D, \underline{v}_D)$ is convex, and therefore for any $(\underline{u}_D, \underline{v}_D) \in \underline{V}_D^2$, one has

$$ a_D^\phi(\underline{u}_D, \underline{u}_D) - a_D^\phi(\underline{v}_D, \underline{v}_D) \leq 2a_D^\phi(\underline{u}_D - \underline{v}_D, \underline{u}_D). $$

Using this relation with $\underline{u}_D = \underline{\phi}_D - \underline{\phi}_D^e$ and $\underline{v}_D = \underline{\phi}_D^n - \underline{\phi}_D^e \in \underline{V}_D$, we obtain the following estimate on the third term of (2.46):

$$ \frac{1}{2} (a_D^\phi(\underline{\phi}_D - \underline{\phi}_D^e, \underline{\phi}_D - \underline{\phi}_D^e) - a_D^\phi(\underline{\phi}_D^n - \underline{\phi}_D^e, \underline{\phi}_D^n - \underline{\phi}_D^e)) \leq a_D^\phi(\underline{\phi}_D - \underline{\phi}_D^n, \underline{\phi}_D - \underline{\phi}_D^e). $$

We recall that the electrostatic potentials are defined by (2.39), and that $\underline{\phi}_D - \underline{\phi}_D^n = \underline{\psi}_D - \underline{\psi}_D^n$.

Hence, letting $\underline{v}_D = \underline{\phi}_D - \underline{\phi}_D^e \in \underline{V}_{D,0}$ and using the linearity of $a_D^\phi$ with respect to the first variable, one gets

$$ \begin{aligned} a_D^\phi(\underline{\phi}_D - \underline{\phi}_D^n, \underline{v}_D) &= a_D^\phi(\underline{\psi}_D, \underline{v}_D) - a_D^\phi(\underline{\psi}_D^n, \underline{v}_D) \\ &= L_D^\phi(\underline{v}_D) + \int_{\Omega} (P_M - N_M)\underline{v}_M - L_D^\phi(\underline{v}_D) - \int_{\Omega} (P_M^n - N_M^n)\underline{v}_M \\ &= \int_{\Omega} (P_M - P_M^n)\underline{v}_M - \int_{\Omega} (N_M - N_M^n)\underline{v}_M. \end{aligned} $$
Combining the previous estimates, we obtain

$$ \mathbb{E}(\underline{N}_D, \underline{P}_D) - \mathbb{E}(\underline{N}_D^n, \underline{P}_D^n) \le \int_{\Omega} (h(N_M) - \phi_M^e - \alpha_N - v_M)(N_M - N_M^n) \\ + \int_{\Omega} (h(P_M) + \phi_M^e - \alpha_P + v_M)(P_M - P_M^n). $$

But, by definition of $\underline{v}_D$, we have $h(N_M) - \phi_M^e - \alpha_N - v_M = h(N_M) - \phi_M - \alpha_N = w_M^N$ and $h(P_M) + \phi_M^e - \alpha_P + v_M = w_M^P$, therefore (2.45) holds. $\square$

One can now state the following a priori result about the dissipation of entropy, which ensures that (2.33) holds.

**Proposition 11 (Entropy dissipation).** Assume that $(\underline{\underline{N}}_D^n, \underline{\underline{P}}_D^n, \underline{\underline{\phi}}_D^n)_{n \in \mathbb{N}} \in \underline{\underline{V}}_D^3$ is a solution to (2.28). Then, the following entropy-entropy dissipation relation holds:

$$ \forall n \ge 0, \quad \frac{\mathbb{E}(\underline{\underline{N}}_D^{n+1}, \underline{\underline{P}}_D^{n+1}) - \mathbb{E}(\underline{\underline{N}}_D^n, \underline{\underline{P}}_D^n)}{\Delta t} \le -\mathbb{D}(\underline{\underline{N}}_D^{n+1}, \underline{\underline{P}}_D^{n+1}). $$

*Proof.* Since the discrete quasi-Fermi potentials $w_D^{N,n+1}$ and $\underline{\underline{w}}_D^{P,n+1}$ are in $\underline{\underline{V}}_{D,0}$ by (2.28a) and (2.28b), one can use them as test functions in (2.28d) and (2.28e): summing the two identities, and noticing that $w_M^{N,n+1} + w_M^{P,n+1} = h(N_M^{n+1}) + h(P_M^{n+1})$ because of (2.5), one gets that

$$ \int_{\Omega} \frac{N_M^{n+1} - N_M^n}{\Delta t} w_M^{N,n+1} + \int_{\Omega} \frac{P_M^{n+1} - P_M^n}{\Delta t} w_M^{P,n+1} = -\mathbb{D}(\underline{\underline{N}}_D^{n+1}, \underline{\underline{P}}_D^{n+1}). $$

One concludes using Lemma 5. $\square$

From this entropic structure follows an important preservation property of the scheme (2.28).

**Proposition 12 (Thermodynamic consistency).** Let $(\underline{\underline{N}}_D^\infty, \underline{\underline{P}}_D^\infty, \underline{\underline{\phi}}_D^\infty) \in \underline{\underline{V}}_D^3$ be a discrete stationary state of the scheme (2.28), in the sense that there exist $(\underline{\underline{N}}_D^n, \underline{\underline{P}}_D^n, \underline{\underline{\phi}}_D^n)_{n \in \mathbb{N}} \in \underline{\underline{V}}_D^3$ a solution to (2.28) and $n_\infty \in \mathbb{N}$ such that,

$$ \forall n \ge n_\infty, \quad (\underline{\underline{N}}_D^n, \underline{\underline{P}}_D^n, \underline{\underline{\phi}}_D^n) = (\underline{\underline{N}}_D^\infty, \underline{\underline{P}}_D^\infty, \underline{\underline{\phi}}_D^\infty). $$

Then the discrete stationary state coincides with the discrete thermal equilibrium:

$$ (\underline{\underline{N}}_D^\infty, \underline{\underline{P}}_D^\infty, \underline{\underline{\phi}}_D^\infty) = (\underline{\underline{N}}_D^e, \underline{\underline{P}}_D^e, \underline{\underline{\phi}}_D^e). \quad (2.47) $$

*Proof.* Let $n \ge n_\infty$, we have

$$ \mathbb{E}(\underline{\underline{N}}_D^{n+1}, \underline{\underline{P}}_D^{n+1}) - \mathbb{E}(\underline{\underline{N}}_D^n, \underline{\underline{P}}_D^n) = \mathbb{E}(\underline{\underline{N}}_D^\infty, \underline{\underline{P}}_D^\infty) - \mathbb{E}(\underline{\underline{N}}_D^\infty, \underline{\underline{P}}_D^\infty) = 0. $$

According to Proposition 11, one has by positivity of the dissipation

$$ \mathbb{D}(\underline{\underline{N}}_D^\infty, \underline{\underline{P}}_D^\infty) = \mathbb{D}(\underline{\underline{N}}_D^{n+1}, \underline{\underline{P}}_D^{n+1}) = 0. $$

But, since all the terms of the dissipation are non-negative (see (2.44)), it means that

$$ T_D^N(\underline{\underline{N}}_D^\infty, \underline{\underline{w}}_D^{N,\infty}, \underline{\underline{w}}_D^{N,\infty}) = T_D^P(\underline{\underline{P}}_D^\infty, \underline{\underline{w}}_D^{P,\infty}, \underline{\underline{w}}_D^{P,\infty}) = 0, $$
where $(w_D^{N,\infty}, w_D^{P,\infty}) \in \underline{V}_{D,0}^2$ are the quasi-Fermi potentials associated to the discrete stationary densities $(\underline{N}_D^\infty, \underline{P}_D^\infty)$. Now, let $K \in \mathcal{M}$. Since all the coordinates of $\underline{N}_K^\infty$ are positive, $r_K(\underline{N}_K^\infty)$ is also positive, and by definition (2.22) of $T_K^N$, we get that

$$a_K^N(\underline{w}_K^{N,\infty}, \underline{w}_K^{N,\infty}) = 0.$$

Using the local coercivity estimates (2.15), we infer that $|\underline{w}_K^{N,\infty}|_{1,K}^2 = 0$. This holds for any cell $K \in \mathcal{M}$, so we have $|\underline{w}_D^{N,\infty}|_{1,D} = 0$, but $\underline{w}_D^{N,\infty} \in \underline{V}_{D,0}$, therefore $\underline{w}_D^{N,\infty} = \underline{0}_D$. A similar result holds for $\underline{w}_D^{P,\infty}$. Thus, one has $\underline{w}_D^{N,\infty} = \underline{w}_D^{P,\infty} = \underline{0}_D$ and, by Remark 12, it means that $(\underline{w}_D^{N,\infty}, \underline{w}_D^{P,\infty})$ are equal to the quasi-Fermi potentials at thermal equilibrium. By correspondence between densities and potential discussed in Section 2.3.2, we get (2.47). $\square$

**Remark 16 (Preservation of the thermal equilibrium).** The statement of Proposition 12 means that the scheme (2.28) preserves the thermal equilibrium: the only admissible discrete stationary state is the discrete thermal equilibrium. This remarkable feature is sometimes called thermodynamic consistency in the context of TPFA schemes (see [127]). In such a framework, this result is expressed in the following way: if the numerical fluxes of the convection-diffusion equations vanish, then the discrete quasi-Fermi potentials are constant (equal to zero). For our HFV scheme, this statement still holds, using the numerical fluxes (2.5.1) defined in Section 2.5.1. In fact, both statements in terms of fluxes and stationary states are equivalent. This notion of thermodynamic consistency can indeed be reformulated using the entropic property of the schemes: if the discrete dissipation associated to a discrete solution vanishes, then the solution is the discrete thermal equilibrium.

### 2.4.2 Existence of solutions

The goal of this section is to prove the existence result and the estimates (2.36) of Theorem 7. The proof proposed here is an adaptation of the one introduced in [70] for similar schemes in the context of linear advection-diffusion equations. According to the discussion in Section 2.3.2, it is equivalent to seek discrete solutions in terms of densities or in terms of quasi-Fermi potentials. Therefore, we will seek discrete quasi-Fermi potentials in the whole space $\underline{V}_{D,0} \times \underline{V}_{D,0}$ instead of discrete densities in the set of $I_h$-valued elements of $\underline{V}_D \times \underline{V}_D$ (which is not a vector space).

In the sequel, we consider the vector space $\underline{V}_{D,0}^2$, and denote by $\underline{w}_D = (\underline{w}_D^N, \underline{w}_D^P) \in \underline{V}_{D,0}^2$ a typical element. We endow $\underline{V}_{D,0}^2$ with the following natural inner product and norm:

$$\langle \underline{w}_D, \underline{v}_D \rangle = \sum_{K \in \mathcal{M}} (w_K^N v_K^N + w_K^P v_K^P) + \sum_{\sigma \in \mathcal{E}_{int} \cup \mathcal{E}_{ext}^N} (w_\sigma^N v_\sigma^N + w_\sigma^P v_\sigma^P) \text{ and } \|\underline{w}_D\| = \sqrt{\langle \underline{w}_D, \underline{w}_D \rangle}.$$

One can also endow $\underline{V}_{D,0}^2$ with the $I^\infty$ norm $\|\cdot\|_\infty$ defined by

$$\|\underline{w}_D\|_\infty = \max_{K \in \mathcal{M}, \sigma \in \mathcal{E}_{int} \cup \mathcal{E}_{ext}^N} (|w_K^N|, |w_K^P|, |w_\sigma^N|, |w_\sigma^P|).$$

Since $\underline{V}_{D,0}^2$ is a finite-dimensional vector space of dimension $2(|\mathcal{M}| + |\mathcal{E}_{int}| + |\mathcal{E}_{ext}^N|)$ there exists a constant $c_{dim}$ only depending on $\dim(\underline{V}_{D,0})$ such that

$$\forall \underline{w}_D \in \underline{V}_{D,0}^2, c_{dim} \|\underline{w}_D\| \le \|\underline{w}_D\|_\infty. \qquad (2.48)$$

We can now establish the existence result by induction on $n$. Let $\underline{N}_D^n$ and $\underline{P}_D^n$ be two $I_h$-valued vectors, we want to show that there exists a solution to (2.28a)-(2.28f). Given $\underline{w}_D =$
$(w_D^N, w_D^P) \in \underline{V}_{D,0}^2$, (and $(\underline{N}_D, \underline{P}_D, \underline{\psi}_D)$ the associated discrete densities and electrostatic potential),
the application

$$
\underline{w}_D \mapsto \int_{\Omega} \frac{N_M - N_M^n}{\Delta t} v_M^N + T_D^N(\underline{N}_D, \underline{w}_D^N, \underline{v}_D^N) + \int_{\Omega} R(N_M, P_M) v_M^N \\
+ \int_{\Omega} \frac{P_M - P_M^n}{\Delta t} v_M^P + T_D^P(\underline{P}_D, \underline{w}_D^P, \underline{v}_D^P) + \int_{\Omega} R(N_M, P_M) v_M^P
$$

is a continuous linear form on $\underline{V}_{D,0}^2$. Hence, there exists a unique $\mathcal{G}(\underline{w}_D) \in \underline{V}_{D,0}^2$ such that

$$
\forall \underline{w}_D \in \underline{V}_{D,0}^2, \langle \mathcal{G}(\underline{w}_D), \underline{w}_D \rangle = \int_{\Omega} \frac{N_M - N_M^n}{\Delta t} v_M^N + T_D^N(\underline{N}_D, \underline{w}_D^N, \underline{v}_D^N) + \int_{\Omega} R(N_M, P_M) v_M^N \\
+ \int_{\Omega} \frac{P_M - P_M^n}{\Delta t} v_M^P + T_D^P(\underline{P}_D, \underline{w}_D^P, \underline{v}_D^P) + \int_{\Omega} R(N_M, P_M) v_M^P, \quad (2.49)
$$

and it is straightforward to see that $\mathcal{G}(\underline{w}_D) = (\underline{0}_D, \underline{0}_D)$ if and only if $(\underline{N}_D, \underline{P}_D, \underline{\phi}_D)$ is a solution to (2.28a)-(2.28f). Moreover, by the continuity result of Theorem 8, the vector field $\mathcal{G}: \underline{V}_{D,0}^2 \to \underline{V}_{D,0}^2$ is continuous. As in [70], our proof relies on two key results. The first one, that can be found, e.g., in [119, Section 9.1], is a corollary of Brouwer's fixed-point theorem.

**Lemma 6.** Let $k$ be a positive integer and let $P: \mathbb{R}^k \to \mathbb{R}^k$ be a continuous vector field. Assume that there is $r > 0$ such that

$$
P(x) \cdot x \ge 0 \quad \text{if } |x| = r.
$$

Then, there exists a point $x_0 \in \mathbb{R}^k$ such that $P(x_0) = 0$ and $|x_0| \le r$.

The second lemma, whose proof is postponed until Appendix 2.A, establishes that bounds on
the discrete entropy and dissipation terms imply bounds on the discrete quasi-Fermi potentials.

**Lemma 7.** Let $(\underline{w}_D^N, \underline{w}_D^P) \in \underline{V}_{D,0}^2$, and assume that there exists $B_\# \ge 0$ such that

$$
\mathbb{E}(\underline{w}_D^N, \underline{w}_D^P) \le B_\# \quad \text{and} \quad \mathbb{D}(\underline{w}_D^N, \underline{w}_D^P) \le B_\# \tag{2.50}
$$

Then, there exists $C_{\#} > 0$, depending on the data, $B_{\#}$ and $\mathcal{D}$ such that

$$
-C_{\#} \underline{1}_D \le \underline{w}_D^N \le C_{\#} \underline{1}_D, \quad -C_{\#} \underline{1}_D \le \underline{w}_D^P \le C_{\#} \underline{1}_D \text{ and } -C_{\#} \underline{1}_D \le \underline{\phi}_D \le C_{\#} \underline{1}_D,
$$

where $\underline{\phi}_D$ is the electrostatic potential associated to $(\underline{w}_D^N, \underline{w}_D^P)$.

We are now in position to prove the existence of solutions to the scheme and the estimates on
these solutions.

*Proof of existence - Theorem 7.* In order to use the result of Lemma 6, notice that by definition of *G*, and because $w_M^N + w_M^P = h(N_M) + h(P_M)$, we have

$$
\forall \underline{w}_D \in \underline{V}_{D,0}^2, \langle \mathcal{G}(\underline{w}_D), \underline{w}_D \rangle = \int_{\Omega} \frac{N_M - N_M^n}{\Delta t} w_M^N + \int_{\Omega} \frac{P_M - P_M^n}{\Delta t} w_M^P + \mathbb{D}(\underline{w}_D^N, \underline{w}_D^P).
$$

Since $(\underline{w}_D^N, \underline{w}_D^P)$ are the quasi-Fermi potentials associated to $(\underline{N}_D, \underline{P}_D)$, the result of Lemma 5
ensures that

$$ \forall \underline{w}_D \in \underline{V}_{D,0}^2, \langle \mathcal{G}(\underline{w}_D), \underline{w}_D \rangle \ge \frac{\mathbb{E}(\underline{w}_D^N, \underline{w}_D^P) - \mathbb{E}(\underline{N}_D^n, \underline{P}_D^n)}{\Delta t} + \mathbb{D}(\underline{w}_D^N, \underline{w}_D^P). \quad (2.51) $$

Let $B^n = \mathbb{E}(\underline{N}_D^n, \underline{P}_D^n)$. According to Lemma 7, there exists $C^n > 0$ depending on the data, $\mathcal{D}$ and $\Delta t$ such that

$$ \text{if } \mathbb{E}(\underline{w}_D^N, \underline{w}_D^P) + \Delta t \mathbb{D}(\underline{w}_D^N, \underline{w}_D^P) \le B^n, \text{ then } \|\underline{w}_D\|_{\infty} \le C^n. \quad (2.52) $$

Now, letting $r^n = 2C^n/c_{\text{dim}}$, one can notice that for any $\underline{w}_D \in \underline{V}_{D,0}^2$ such that $\|\underline{w}_D\| = r^n$, we have $C^n < 2C^n = c_{\text{dim}}\|\underline{w}_D\| \le \|\underline{w}_D\|_{\infty}$. Therefore, by (2.51) and contraposition of (2.52), if $\|\underline{w}_D\| = r^n$, one has

$$ \Delta t \langle \mathcal{G}(\underline{w}_D), \underline{w}_D \rangle \ge \mathbb{E}(\underline{w}_D^N, \underline{w}_D^P) + \Delta t \mathbb{D}(\underline{w}_D^N, \underline{w}_D^P) - B^n > B^n - B^n \ge 0. $$

Thus, we can use the result of Lemma 6 applied to the vector field $\mathcal{G}$, and there exists at least one solution $(\underline{N}_D^{n+1}, \underline{P}_D^{n+1}, \underline{\phi}_D^{n+1})$ to (2.28a)-(2.28f).

To prove the bounds (2.36), it suffices to note that the dissipation relation (2.33) implies that

$$ \mathbb{E}^{n+1} \le \mathbb{E}^n \le \mathbb{E}^0 \text{ and } \mathbb{D}^{n+1} \le \frac{\mathbb{E}^0}{\Delta t}. $$

Hence, by Lemma 7, there exists $C_{\#}$ depending on the data, $\mathcal{D}$, $\Delta t$ and $\mathbb{E}^0$ (but not on $n$) such that

$$ -C_{\#} \underline{1}_D \le \underline{w}_D^{N,n+1} \le C_{\#} \underline{1}_D, \quad -C_{\#} \underline{1}_D \le \underline{w}_D^{P,n+1} \le C_{\#} \underline{1}_D \text{ and } -C_{\#} \underline{1}_D \le \underline{\phi}_D^{n+1} \le C_{\#} \underline{1}_D. $$

Therefore, by (2.29), $g(-2C_{\#} + \alpha_N)\underline{1}_D \le \underline{N}_D^{n+1} \le g(2C_{\#} + \alpha_N)\underline{1}_D$, and an analogous estimate holds for $\underline{P}_D^{n+1}$. $\square$

### 2.4.3 Long-time behaviour

In this section, we analyse the long-time behaviour of the scheme (2.28). The analysis proposed here relies on an adaptation of the arguments presented in [27] in the framework of TPFA schemes.

**Theorem 9 (Discrete long-time behaviour).** Let $(\underline{N}_D^n, \underline{P}_D^n, \underline{\phi}_D^n)_{n \in \mathbb{N}} \in \underline{V}_D^3$ be a solution to the coupled scheme (2.28). There exists $v_{D,\Delta t} > 0$ depending on the data, $\mathcal{D}$ and $\Delta t$ such that

$$ \forall n \in \mathbb{N}, \mathbb{E}^{n+1} \le (1 + v_{D,\Delta t} \Delta t)^{-1} \mathbb{E}^n. \quad (2.53) $$

Moreover, the discrete solution converges geometrically fast towards the discrete thermal equilibrium: there exists a positive constant $c_{D,\Delta t}$ depending on the data, $\mathcal{D}$ and $\Delta t$ such that

$$ \forall n \in \mathbb{N}, \|N_M^n - N_M^e\|_{L^2(\Omega)}^2 + \|P_M^n - P_M^e\|_{L^2(\Omega)}^2 + \|\phi_M^n - \phi_M^e\|_{L^2(\Omega)}^2 \le c_{D,\Delta t} \mathbb{E}^0(1 + v_{D,\Delta t} \Delta t)^{-n}. \quad (2.54) $$

*Proof.* In this proof, we will denote by $c$ a generic constant depending on the data, $\mathcal{D}$ and $\Delta t$ (but not on $n$). First, recall that the solutions satisfy the bounds (2.36). One can notice (using a Taylor expansion of $H$) that there exists two positive constants $\tilde{c}_1$ and $\tilde{c}_2$ such that for any $(x, y) \in [M_b, M_{\#}]$,

$$ \tilde{c}_1(x-y)^2 \le H(x) - H(y) - h(y)(x-y) \le \tilde{c}_2(x-y)^2. \quad (2.55) $$

Therefore, letting $\hat{\mathbb{E}}^n = \|N_M^n - N_M^e\|_{L^2(\Omega)}^2 + \|P_M^n - P_M^e\|_{L^2(\Omega)}^2 + \frac{1}{2}a_D^\phi(\underline{\phi}_D - \underline{\phi}_D^e, \underline{\phi}_D - \underline{\phi}_D^e)$ and using
(2.36) and (2.55), we deduce that there exist $c_1$ and $c_2$ positive such that

$$c_1 \hat{\mathbb{E}}^n \le \mathbb{E}^n \le c_2 \hat{\mathbb{E}}^n. \qquad (2.56)$$

Since $\mathbb{E}(\underline{N}_D^e, \underline{P}_D^e) = 0$, by Lemma 5 (used with $(\underline{N}_D, \underline{P}_D) = (\underline{N}_D^n, \underline{P}_D^n)$ and $(\underline{N}_D^n, \underline{P}_D^n) = (\underline{N}_D^e, \underline{P}_D^e)$), one gets

$$c_1 \hat{\mathbb{E}}^n \le \mathbb{E}^n \le \int_{\Omega} (N_M^n - N_M^e) w_M^{N,n} + \int_{\Omega} (P_M^n - P_M^e) w_M^{P,n}.$$

Then, we use Young inequality (with scaling $c_1$) to obtain

$$c_1 \hat{\mathbb{E}}^n \le \frac{c_1}{2} \left( \|N_M^n - N_M^e\|_{L^2(\Omega)}^2 + \|P_M^n - P_M^e\|_{L^2(\Omega)}^2 \right) + \frac{1}{2c_1} \left( \|w_M^{N,n}\|_{L^2(\Omega)}^2 + \|w_M^{P,n}\|_{L^2(\Omega)}^2 \right).$$

Notice that $\|N_M^n - N_M^e\|_{L^2(\Omega)}^2 + \|P_M^n - P_M^e\|_{L^2(\Omega)}^2 \le \hat{\mathbb{E}}^n$, therefore we have

$$\frac{c_1}{2} \hat{\mathbb{E}}^n \le \frac{1}{2c_1} \left( \|w_M^{N,n}\|_{L^2(\Omega)}^2 + \|w_M^{P,n}\|_{L^2(\Omega)}^2 \right).$$

Combining this estimate with (2.56), we deduce that

$$\frac{c_1}{c_2} \mathbb{E}^n \le c_1 \hat{\mathbb{E}}^n \le \frac{1}{c_1} \left( \|w_M^{N,n}\|_{L^2(\Omega)}^2 + \|w_M^{P,n}\|_{L^2(\Omega)}^2 \right). \qquad (2.57)$$

On the other hand, one has $\mathbb{D}^n \ge T_D^N(\underline{N}_D^n, \underline{w}_D^{N,n}, \underline{w}_D^{P,n}) + T_D^P(\underline{P}_D^n, \underline{w}_D^{P,n}, \underline{w}_D^{P,n})$, and, for any $K \in \mathcal{M}$, using (2.22) alongside with the positivity of $r_K(\underline{N}_K)$ and the local coercivity (2.15) of $a_K^n$, one gets that

$$T_K^N(\underline{N}_K^n, \underline{w}_K^{N,n}, \underline{w}_K^{P,n}) \ge r_K(\underline{N}_K^n) a_K^N(\underline{w}_K^{N,n}, \underline{w}_K^{P,n}) \ge \alpha_b \lambda_b r_K(\underline{N}_K^n) |\underline{w}_K^{N,n}|_{1,K}^2.$$

Moreover, by definition of $r_K$, properties on $m$ and $m_h$ and the bounds (2.36) on $\underline{N}_D^n$, we deduce that

$$r_K(\underline{N}_K^n) = \sum_{\sigma \in \mathcal{E}_K} \frac{m(N_K^n, N_\sigma^n)}{|\mathcal{E}_K|} \ge \sum_{\sigma \in \mathcal{E}_K} \frac{m_h(N_K^n, N_\sigma^n)}{|\mathcal{E}_K|} \ge \sum_{\sigma \in \mathcal{E}_K} \frac{m_h(M_b, M_b)}{|\mathcal{E}_K|} \ge M_b.$$

The same holds for the term $T_D^P(\underline{P}_D^n, \underline{w}_D^{P,n}, \underline{w}_D^{P,n})$, therefore by the discrete Poincaré inequality (2.17) one has

$$\mathbb{D}^n \ge \alpha_b \lambda_b M_b \left( |\underline{w}_D^{N,n}|_{1,D}^2 + |\underline{w}_D^{P,n}|_{1,D}^2 \right) \ge \frac{\alpha_b \lambda_b M_b}{C_{P,\Gamma D}^2} \left( \|w_M^{N,n}\|_{L^2(\Omega)}^2 + \|w_M^{P,n}\|_{L^2(\Omega)}^2 \right). \qquad (2.58)$$

Combining (2.57) and (2.58), we get the following control on the entropy:

$$\nu_{D,\Delta t} \mathbb{E}^n \le \mathbb{D}^n \text{ with } \nu_{D,\Delta t} = \frac{\alpha_b \lambda_b M_b c_1^2}{c_2 C_{P,\Gamma D}^2}.$$

From this inequality and the entropy dissipation relation (2.33), we deduce (2.53). To get (2.54), notice that $c_1 \hat{\mathbb{E}}^n \le \mathbb{E}^n$, so using the discrete Poincaré inequality alongside with the coercivity of
$a_D^\phi$, we get

$$ \|N_M^n - N_M^e\|_{L^2(\Omega)}^2 + \|P_M^n - P_M^e\|_{L^2(\Omega)}^2 + \frac{\alpha_b \lambda_b^\phi}{2C_{P,\Gamma^D}^2} \|\phi_M^n - \phi_M^e\|_{L^2(\Omega)}^2 \le \hat{E}^n \le \frac{1}{c_1} \mathbb{E}^n. $$

We conclude by using (2.53).

<figure><img src="image_25.png" /></figure>

**Remark 17** (Non-uniformity with respect to the mesh). The result of Theorem 9 states a geometric convergence of the solutions to the scheme towards the discrete thermal equilibrium. However, the constants involved in (2.53) and (2.54) may depend strongly on the mesh (and in particular on $h_D$), which yields a weaker result than the one of [27, 28] for TPFA schemes. Such a dependency is a reminiscence of the non-monotonicity of the HFV schemes (see Remark 13). In practice, the long-time behaviour of the scheme seems not to depend on the meshsize: see Section 2.5.3.

## 2.5 Numerical results

In this section, we give some numerical evidences of the good behaviour of the scheme (2.28). We use test-cases inspired by the 2D PN-junction studied in [27], whose geometry is described in Figure 2.2. The domain $\Omega$ is the unit square $[0, 1]^2$. For the boundary conditions, we split $\Gamma^D = \Gamma_0^D \cup \Gamma_1^D$ with $\Gamma_0^D = [0, 1] \times \{0\}$ and $\Gamma_1^D = [0, 0.25] \times \{1\}$. For $i \in \{0, 1\}$, we let

$$ N^D = N_i^D, P^D = P_i^D \text{ and } \phi^D = \frac{h(N_i^D) - h(P_i^D)}{2} \text{ on } \Gamma_i^D. $$

To be consistent with the compatibility condition (2.4) we assume that there exists a constant $\alpha_0$ such that $h(N^D) + h(P^D) = \alpha_0$, therefore for given $N^D$ and $\alpha_0$ we set $P^D = g(\alpha_0 - h(N^D))$ on $\Gamma^D$. Thus, we get $\alpha_N = \alpha_P = \frac{\alpha_0}{2}$. If $r \ne 0$, we finally impose that $\alpha_0 = 0$ to satisfy (2.5). We use the

<figure><img src="image_26.png" /><figcaption>Figure 2.2 – The geometry of the PN diode.</figcaption></figure>
following initial condition:

$$N_0(x, y) = N_1^D + (N_0^D - N_1^D)(1 - \sqrt{y}) \text{ and } P_0(x, y) = P_1^D + (P_0^D - P_1^D)(1 - \sqrt{y}).$$

Finally, we use a piecewise constant doping profile C, equal to -1 in the P-region and 1 in the N-region (see Figure 2.2). Concerning the tensors, we assume that the permittivity is isotropic, of the form $\Lambda_{\phi} = \lambda^2 I_2$, where $\lambda > 0$ is the rescaled Debye length. We also assume that the magnetic field is constant, of magnitude $b \ge 0$, and that the rescaled mobilities are equal to 1 ($\mu_N = \mu_P = 1$), therefore the tensors for the convection-diffusion equations are

$$\Lambda_N = \frac{1}{1+b^2} \begin{pmatrix} 1 & b \\ -b & 1 \end{pmatrix} \text{ and } \Lambda_P = \frac{1}{1+b^2} \begin{pmatrix} 1 & -b \\ b & 1 \end{pmatrix}. \qquad (2.59)$$

### 2.5.1 Implementation

In this section, we discuss some practical details concerning the implementation of the schemes introduced in this paper. The mesh families used for the numerical tests are classical Cartesian and triangular families - see for example [154] - and a tilted hexagonal-dominant mesh family, depicted on Figure 2.3. These meshes have convex cells, hence we always choose $x_K$ to be

<figure><img src="image_27.png" /><figcaption>Figure 2.3 – Implementation. Coarsest meshes of each family used in the numerical tests.</figcaption></figure>

the barycentre of K. Moreover, we use a fixed stabilisation parameter $\eta = 1.5$ (see (2.13)). In the simulations showed below, we use the arithmetic mean as an $m$ function for the reconstruction defined in (2.23), therefore for $u_D \in V_D$ and $K \in \mathcal{M}$,

$$r_K(u_K) = \frac{1}{2} \left( u_K + \frac{1}{|\mathcal{E}_K|} \sum_{\sigma \in \mathcal{E}_K} u_{\sigma} \right).$$

For a discussion on other choices of reconstructions we refer to [51, Section 6.2].

#### Finite volume formulation

The hybrid schemes described in this paper define finite volume methods, in the sense that they can be equivalently rewritten under a conservative forms, with local mass balances, flux equilibration at interfaces, and boundary conditions. For more details about this formulation in the framework of a linear Poisson equation, we refer the reader to [123]. Let us detail succinctly this formulation for our schemes. Let $\Lambda$ be a generic uniformly elliptic tensor. For all $K \in \mathcal{M}$,
and all $\sigma \in \mathcal{E}_K$, in the framework of HFV methods, the normal diffusive flux $-\int_{\sigma} \Lambda \nabla u \cdot n_{K,\sigma}$ is approximated by the following numerical flux:

$$F_{K,\sigma}^{\Lambda}(\underline{u}_K) = \sum_{\sigma' \in \mathcal{E}_K} A_K^{\sigma\sigma'}(\underline{u}_K - \underline{u}_{\sigma'}), \qquad (2.60)$$

where the $A_K^{\sigma\sigma'}$ are defined by

$$A_K^{\sigma\sigma'} = \sum_{\sigma'' \in \mathcal{E}_K} |P_{K,\sigma''}| y_K^{\sigma''\sigma} \cdot \Lambda_{K,\sigma''} y_K^{\sigma''\sigma'}, \qquad (2.61)$$

and the $y_K^{\sigma\sigma'} \in \mathbb{R}^d$ only depend on the geometry of the discretisation $\mathcal{D}$ (see, for example, [123, Eq. (2.22)] for an exact definition with $\eta = \sqrt{d}$). For all $K \in \mathcal{M}$, one can express the local discrete bilinear form $a_K^{\Lambda}$ in terms of the local fluxes $(F_{K,\sigma}^{\Lambda})_{\sigma \in \mathcal{E}_K}$: for all $(\underline{u}_K, \underline{v}_K) \in \underline{V}_K^2$,

$$a_K^{\Lambda}(\underline{u}_K, \underline{v}_K) = \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}^{\Lambda}(\underline{u}_K)(\underline{v}_K - \underline{v}_{\sigma}).$$

With these considerations, the scheme (2.31) for the Poisson equation (2.30) writes under the following form:

$$\left\{ \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}^{\Lambda_{\phi}}(\underline{\phi}_K) = |K|(C_K + g(z_K^P - \phi_K) - g(z_K^N + \phi_K)) \quad \forall K \in \mathcal{M}, \right. \qquad (2.62a)$$

$$F_{K,\sigma}^{\Lambda_{\phi}}(\underline{\phi}_K) = -F_{L,\sigma}^{\Lambda_{\phi}}(\underline{\phi}_L) \quad \forall \sigma = K|L \in \mathcal{E}_{\text{int}}, \qquad (2.62b)$$

$$\phi_{\sigma} = \phi_{\sigma}^{D} \quad \forall \sigma \in \mathcal{E}_{\text{ext}}^{D}, \qquad (2.62c)$$

$$F_{K,\sigma}^{\Lambda_{\phi}}(\underline{\phi}_K) = 0 \quad \forall \sigma \in \mathcal{E}_{\text{ext}}^{N} \text{ with } \mathcal{M}_{\sigma} = \{K\}, \qquad (2.62d)$$

where the fluxes are defined by (2.60), and $C_K = \frac{1}{|K|} \int_K C$. The first equation corresponds to local balance, the second imposes the local conservativity of the fluxes at interfaces and the last one enforces the boundary conditions.

Following an analogous approach, we define the nonlinear flux for the advection-diffusion: for any $I_h$-valued $\underline{u}_K \in \underline{V}_K$ and $\underline{\phi}_K \in \underline{V}_K$, we let

$$\mathcal{F}_{K,\sigma}^{\Lambda}(\underline{u}_K, \underline{\phi}_K) = r_K(\underline{u}_K) F_{K,\sigma}^{\Lambda} (h(\underline{u}_K) + \underline{\phi}_K). \qquad (2.63)$$

Therefore, letting $\underline{w}_K = h(\underline{u}_K) + \underline{\phi}_K - \alpha \underline{1}_K$ (with $\alpha \in \mathbb{R}$), one can write the local trilinear form $T_K^{\Lambda}$ in terms of the nonlinear flux: for any $\underline{v}_K \in \underline{V}_K$,

$$T_K^{\Lambda}(\underline{u}_K, \underline{w}_K, \underline{v}_K) = \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{K,\sigma}^{\Lambda}(\underline{u}_K, \underline{\phi}_K)(\underline{v}_K - \underline{v}_{\sigma}).$$
The scheme (2.28) for the drift-diffusion system then writes under the following form:

$$ \forall K \in \mathcal{M}, \qquad |K| \frac{N_K^{n+1} - N_K^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{K,\sigma}^{\Lambda_N}(\underline{N}_K^{n+1}, -\underline{\phi}_K^{n+1}) = -|K|R(N_K^{n+1}, P_K^{n+1}), \quad (2.64a) $$

$$ \forall K \in \mathcal{M}, \qquad |K| \frac{P_K^{n+1} - P_K^n}{\Delta t} + \sum_{\sigma \in \mathcal{E}_K} \mathcal{F}_{K,\sigma}^{\Lambda_P}(\underline{P}_K^{n+1}, \underline{\phi}_K^{n+1}) = -|K|R(N_K^{n+1}, P_K^{n+1}), \quad (2.64b) $$

$$ \forall K \in \mathcal{M}, \qquad \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}^{\Lambda_{\phi}}(\underline{\phi}_K^{n+1}) = |K|(C_K + P_K^{n+1} - N_K^{n+1}), \quad (2.64c) $$

$$ \forall \sigma = K|L \in \mathcal{E}_{int}, \qquad \mathcal{F}_{K,\sigma}^{\Lambda_N}(\underline{N}_K^{n+1}, -\underline{\phi}_K^{n+1}) + \mathcal{F}_{L,\sigma}^{\Lambda_N}(\underline{N}_L^{n+1}, -\underline{\phi}_L^{n+1}) = 0, \quad (2.64d) $$

$$ \forall \sigma = K|L \in \mathcal{E}_{int}, \qquad \mathcal{F}_{K,\sigma}^{\Lambda_P}(\underline{P}_K^{n+1}, \underline{\phi}_K^{n+1}) + \mathcal{F}_{L,\sigma}^{\Lambda_P}(\underline{P}_L^{n+1}, \underline{\phi}_L^{n+1}) = 0, \quad (2.64e) $$

$$ \forall \sigma = K|L \in \mathcal{E}_{int}, \qquad F_{K,\sigma}^{\Lambda_{\phi}}(\underline{\phi}_K^{n+1}) + F_{L,\sigma}^{\Lambda_{\phi}}(\underline{\phi}_L^{n+1}) = 0, \quad (2.64f) $$

$$ \forall \sigma \in \mathcal{E}_{ext}^D, \qquad N_{\sigma}^{n+1} = N_{\sigma}^D, P_{\sigma}^{n+1} = P_{\sigma}^D \text{ and } \underline{\phi}_{\sigma}^{n+1} = \underline{\phi}_{\sigma}^D, \quad (2.64g) $$

$$ \forall \sigma \in \mathcal{E}_{ext}^N, \qquad \mathcal{F}_{K,\sigma}^{\Lambda_N}(\underline{N}_K^{n+1}, -\underline{\phi}_K^{n+1}) = \mathcal{F}_{K,\sigma}^{\Lambda_P}(\underline{P}_K^{n+1}, \underline{\phi}_K^{n+1}) = F_{K,\sigma}^{\Lambda_{\phi}}(\underline{\phi}_K^{n+1}) = 0, \quad (2.64h) $$

where the nonlinear fluxes are defined by (2.63), the initial data $(N_K^0, P_K^0)_{K \in \mathcal{M}}$ are defined as in (2.28g) and the cell $K$ in (2.64h) is such that $\mathcal{M}_{\sigma} = \{K\}$.

Note that these formulations yield nonlinear systems. Thus, we can introduce natural functions $G^{\text{sta}}: \underline{V}_D \to \underline{V}_D$ and $G^{n,\Delta t}: \underline{V}_D^3 \to \underline{V}_D^3$ such that (2.62) rewrites $G^{\text{sta}}(\underline{\phi}_D) = \underline{0}_D$ and (2.64) writes $G^{n,\Delta t}(\underline{N}_D^{n+1}, \underline{P}_D^{n+1}, \underline{\phi}_D^{n+1}) = (\underline{0}_D, \underline{0}_D, \underline{0}_D)$. The two functions are regular on their domains.

### Newton's method and static condensation

The implementation of the nonlinear schemes relies on their finite volume formulation. To fix ideas, we consider the case of the transient scheme (2.64). Given $(\underline{N}_D^n, \underline{P}_D^n) \in \underline{V}_D^2$, we want to solve the nonlinear system $G^{n,\delta t}(\underline{N}_D^{n+1}, \underline{P}_D^{n+1}, \underline{\phi}_D^{n+1}) = (\underline{0}_D, \underline{0}_D, \underline{0}_D)$ (with a time step $\delta t$) instead of $\Delta t$. The resolution of this system relies on a Newton method with time step adaptation.

First, one initialises the method with initial guess $(\underline{N}_D^n, \underline{P}_D^n) \in \underline{V}_D^2$, where the coordinates of $\underline{N}_D^n$ (respectively $\underline{P}_D^n$) are the projections of the coordinates of $\underline{N}_D^n$ (resp. $\underline{P}_D^n$) on $[\epsilon, a - \epsilon]$ (if $a = +\infty$, we project on $[\epsilon, +\infty[$) in order to avoid potential problems due to the singularity of $h$ near 0 and $a$.

The computation of the residues follows the process described below: let us denote by $R_M \in \mathbb{R}^{3|M|}$ and $R_E \in \mathbb{R}^{3|\mathcal{E}|}$ the residue vectors $(r_K^N, r_K^P, r_K^\phi)_{K \in \mathcal{M}}$ and $(r_\sigma^N, r_\sigma^P, r_\sigma^\phi)_{\sigma \in \mathcal{E}}$. They are solution to the following linear block system:

$$ \left( \begin{array}{cc} \mathbb{M}_M & \mathbb{M}_{M,\mathcal{E}} \\ \mathbb{M}_{\mathcal{E},M} & \mathbb{M}_{\mathcal{E}} \end{array} \right) \left( \begin{array}{c} R_M \\ R_{\mathcal{E}} \end{array} \right) = \left( \begin{array}{c} S_M \\ S_{\mathcal{E}} \end{array} \right), \quad (2.65) $$

where $\mathbb{M}_M \in \mathbb{R}^{3|M| \times 3|M|}$, $\mathbb{M}_{M,\mathcal{E}} \in \mathbb{R}^{3|M| \times 3|\mathcal{E}|}$, $\mathbb{M}_{\mathcal{E},M} \in \mathbb{R}^{3|\mathcal{E}| \times 3|M|}$, $\mathbb{M}_{\mathcal{E}} \in \mathbb{R}^{3|\mathcal{E}| \times 3|\mathcal{E}|}$, and $S_M$ and $S_{\mathcal{E}}$ are vectors of size $3|M|$ and $3|\mathcal{E}|$ issued from the previous iteration. By construction, the matrix $\mathbb{M}_M$ is block diagonal with $3 \times 3$ diagonal blocks, which are expected to be invertible. Therefore, this matrix can be inverted at a very low computational cost, inverting only small matrices. Thus, we
can eliminate the cell unknowns, noticing that

$$R_M = M_M^{-1}(S_M - M_{M,\varepsilon}R_\varepsilon). \qquad (2.66)$$

Using this relation, one shows that $R_\varepsilon$ is the solution to the following linear system:

$$(M_\varepsilon - M_{\varepsilon,M} M_M^{-1} M_{M,\varepsilon}) R_\varepsilon = S_\varepsilon - M_{\varepsilon,M} M_M^{-1} S_M, \qquad (2.67)$$

where $M_D = M_\varepsilon - M_{\varepsilon,M}M_M^{-1}M_{M,\varepsilon}$, the Schur complement of the block $M_M$, is an invertible matrix of size $3|\varepsilon| \times 3|\varepsilon|$. In practice, we solve the linear system (2.67) using an LU factorization algorithm, and we use the solution $R_\varepsilon$ in order to compute $R_M$ from (2.66). This technique, called static condensation, allows one to replace a system of size $3(|\varepsilon| + |M|)$ by a system of size $3|\varepsilon|$ without any additional fill-in. As a stopping criterion for the Newton iterations, we compare the $l^\infty$ relative norm of the residue with a threshold tol. If the method does not converge after $i_{max}$ iterations, we divide the time step by 2 and restart the resolution. When the method converges, one can proceed to the approximation of $(\underline{N}_D^{n+2}, \underline{P}_D^{n+2}, \underline{\phi}_D^{n+2})$, with an initial time step of $\min(\Delta t, 1.4 \times \delta t)$.

The initial time step (used to compute $(\underline{N}_D^1, \underline{P}_D^1, \underline{\phi}_D^1)$) is $\Delta t$. In practice, we use $\varepsilon = 10^{-9}$, $i_{max} = 50$ and $\text{tol} = 10^{-10}$.

For the discrete thermal equilibrium, we solve the nonlinear system $G^{\text{sta}}(\underline{\phi}_D) = \underline{0}_D$ (with $(z^P, z^N) = (\alpha_N, \alpha_P)$) using a similar approach, based on a Newton's method alongside with a continuation method. Moreover, note that the counterpart in this case of the matrix $M_M$ is diagonal with non-zero entries, so its inversion is straightforward.

### 2.5.2 Proof of concept

In this section, we are interested in qualitative and quantitative properties of the discrete solutions, and we present the profiles of some computed discrete solutions.

<figure><img src="image_28.png" /><figcaption>Figure 2.4 – Test-case 1. Evolution of the discrete density of holes</figcaption></figure>

For the test-case 1, we use a set of data introduced in [68], with Boltzmann statistics ($h = \log$), no recombination ($r = 0$), no magnetic field ($b = 0$), $\lambda = 1$, boundary values $N_0^D = 0.9$, $N_1^D = 0.1$ and $\alpha_0 = \log(N_0^D \times N_1^D)$. We perform a simulation on a tilted hexagonal-dominant mesh constituted of 280 cells, with $\Delta t = 0.1$. On Figure 2.4, we show the profiles of $P_M^n$ for different
values of $n$. The evolution is in accordance with results obtained in the TPFA context (see [68], and notice that we use a different initial condition). Note that the discrete density remains positive, as expected. In fact, even on coarser meshes, the positivity of the densities is preserved (both for the cell and edge values).

Now, we want to assert the robustness of the scheme with respect to the $h$ function and the anisotropy. We consider the test-case 2, with Blakemore statistics ($h(s) = \log(\frac{s}{1-\gamma s})$, $\gamma = 0.27$),

<figure><img src="image_29.png" /><figcaption>Figure 2.5 – Test-case 2 ($b = 1$). Evolution of the discrete density of holes (note that the scale varies from a figure to the other)</figcaption></figure>

a SRH recombination term ($r(N, P) = \frac{10}{1+N+P}$) and a strong magnetic field $b = 1$. We also use a (realistic) small Debye length $\lambda = 0.05$. In order to check the fact that the discrete densities are $I_h$-valued (here, $a = 1/\gamma \approx 3.7$), we consider boundary values close to the maximum admissible densities: $N_0^D = 3.5$ and $N_1^D = 1.5$. Since $r \neq 0$, we let $\alpha_0 = 0$. We perform a simulation on a refined triangular mesh (57 344 cells), with a time step $\Delta t = 0.1$. We show the profile of $P_M^n$ for different values of $n$ in Figure 2.5. Notice that this test-case is subject to boundary layers (essentially because $\lambda$ is small and $b$ is big, see the discussion below and Table 2.1), therefore the scheme performs numerous adaptations of the time step at the beginning of the simulation (the first admissible time step is $\delta t \approx 10^{-8}$). Moreover, one can see a rotation movement for the density of the holes, which is in accordance with the expected physical effect of the magnetic field.

To give more quantitative informations, we show in Figure 2.6 the evolution of the minimal
<figure><img src="image_30.png" /><figcaption>Figure 2.6 – Test-case 2 (b = 1). Evolution of the discrete extremal values, time step and cost</figcaption></figure>

and maximal values, along with the time step and the number of Newton’s iterations needed to compute the solutions at a given time. First, note that we display each of the steps (from $t^n$ to $t^{n+1}$) used to compute the solutions on the time interval $[0, 1]$, there are 63 of them for this simulation. The extremal values account for the cell and edge unknowns: we let

$$ \min \underline{N}_D^n = \min \left( \min_{K \in \mathcal{M}} N_K^n, \min_{\sigma \in \mathcal{E}} N_\sigma^n \right) \text{ and } \max \underline{N}_D^n = \max \left( \max_{K \in \mathcal{M}} N_K^n, \max_{\sigma \in \mathcal{E}} N_\sigma^n \right), $$

as well as analogous definitions for the holes density. In order to observe the extremal values on relevant scales, we look at the minimum and the distance between $a$ (the maximal density allowed by the model) and the maximal densities computed, and print them on a log scale. As expected from Theorem 7, the values of the discrete densities stay in $I_h = ]0, a[$. In particular, our scheme is not subject to the same lack of positivity as the scheme of [143] when the magnetic field is intense. Moreover, it is remarkable to note that the scheme seems very robust: the densities are close to the limit values, at a distance reaching 7.56e-8 in the most difficult situation. Note that
this value mean that in practice, we do not perform the projection step described in Section 2.5.1 to initialise the Newton's methods, since the projection threshold $\epsilon = 10^{-9}$ is smaller than 7.56e-8. Using larger values for $\epsilon$ could perhaps improve the convergence of Newton's methods, but the impact of such a modification has not been investigated here. As previously announced, one can see that the first effective time step needed to compute $(\underline{N}_D^1, \underline{P}_D^1, \underline{\phi}_D^1)$ is relatively small: $0.1 \times 2^{-22} \approx 2.38e-8$. It means that the time step adaptation procedure needed to perform 22 time step reductions before managing to compute a solution. On the other hand, we can see that after this initial reduction, there are only five others time step reductions (one around 1.5e-5, one at 2e-5, two around 1.5e-4, and a last one at 9e-4). At the end of the simulation, around $t = 0.3$, the time step reach its maximal value $\Delta t = 0.1$ and there is no more time step adaptation. On the third graph, we show the number of Newton iterations needed to compute the discrete solution from one time to another. Note that we do not take into account the iterations used in non-convergent Newton's methods (i.e. methods that leads to a time step reduction). For all the time step reductions of this simulation, the Newton's methods do not converge because at least one of the computed discrete densities was not $I_h$-valued. On average, the convergent Newton's methods converge in 10.5 iterations. Finally, we can see on Figure 2.6 a clear correlation between the extremal values reached by the discrete densities and the number of Newton iterations needed to compute the solution. This can be explained by at least one fact: since $h$ is singular in 0 and $a$, the values of its derivative near these limit values blow up. Therefore, the Jacobian matrix used in the Newton's method tends to be ill-conditioned, which induces numerical instabilities. We can give a last remark on this test-case: the most extreme values observed in Figure 2.6 are located near to the boundary and appear at the beginning of the simulation (see for example the minimal values on Figure 2.5b, on the left and right sides of the square). Hence, the difficulties are essentially due to the presence of the (strong) magnetic field, which induces rotation of the charges and creates boundary layers.

<table><caption>Table 2.1 – Test-case 2. Comparison of the extremal values and costs for different magnetic fields</caption><thead><tr><th>Magnetic field intensity</th><th>b = 0</th><th>b = 0.5</th><th>b = 1</th></tr></thead><tbody><tr><td>min {min <u>N</u><sub>D</sub><sup>n</sup> | 0 ≤ <u>t</u><sup>n</sup> ≤ 1}</td><td>3.20e-1</td><td>5.36e-3</td><td>1.41e-4</td></tr><tr><td>min {min <u>P</u><sub>D</sub><sup>n</sup> | 0 ≤ <u>t</u><sup>n</sup> ≤ 1}</td><td>1.56e-2</td><td>2.53e-3</td><td>7.56e-8</td></tr><tr><td>min {a - max <u>N</u><sub>D</sub><sup>n</sup> | 0 ≤ <u>t</u><sup>n</sup> ≤ 1}</td><td>2.04e-1</td><td>4.17e-3</td><td>2.18e-6</td></tr><tr><td>min {a - max <u>P</u><sub>D</sub><sup>n</sup> | 0 ≤ <u>t</u><sup>n</sup> ≤ 1}</td><td>2.31</td><td>9.39e-1</td><td>4.23e-2</td></tr><tr><td>Number of steps</td><td>18</td><td>40</td><td>63</td></tr><tr><td>Minimal time step</td><td>3.13e-3</td><td>1.53e-6</td><td>2.38e-8</td></tr><tr><td>Number of initial time step reductions</td><td>5</td><td>16</td><td>22</td></tr><tr><td>Total number of time step reductions</td><td>5</td><td>16</td><td>27</td></tr><tr><td>Maximal number of Newton iterations</td><td>7</td><td>7</td><td>35</td></tr><tr><td>Average number of Newton iterations</td><td>5.05</td><td>5.45</td><td>10.54</td></tr><tr><td>Total cost</td><td>96</td><td>234</td><td>692</td></tr></tbody></table>

To confirm this statement, we perform simulations with the same parameters (still on the time interval $[0, 1]$) except that we consider situations with smaller magnetic field intensities ($b = 0$ and $b = 0.5$). The results are presented in Table 2.1. As expected, it seems that the extremal values are strongly related to the intensity of the magnetic field. It follows that the number of Newton iterations, and hence the global computational cost of the simulation increase with the intensity of the magnetic field. We can also notice that the time step reductions occur only at the
first time step for the moderate magnetic fields ($b = 0$ and $b = 0.5$): the computational difficulties lie in the boundary layers appearing at small times. As previously, note that the “number of iterations” mentioned on the table do not take into account the iterations of non-convergent Newton’s methods. In the last line, we give the “Total cost” of the simulation, that is to say the number of linear systems solved during the simulation, including the iterations of non-convergent Newton’s methods. Taking into account every iterations performed, the simulation with $b = 1$ is basically 7 times more time consuming than the one without magnetic field.

In order to assert the robustness of the method with respect to the mesh, we consider the test case 3, which is characterised by the same physical parameters as the previous one: $h(s) = \log\left(\frac{s}{1-\gamma s}\right)$, $\gamma = 0.27$), $r(N, P) = \frac{10}{1+N+P}$, $b = 1$, $\lambda = 0.05$, $N_0^D = 3.5$, $N_1^D = 1.5$ and $\alpha_0 = 0$. Contrarily to the previous test, we perform a simulation on a tilted hexagonal mesh, constituted of 4192 cells, with a time step $\Delta t = 0.1$. One has to notice that the spatial mesh used here is a “general polygonal” one, in the sense that it is not a admissible mesh for the TPFA scheme. Moreover, the geometry of this mesh is not well-suited with respect to the geometry of the device, since the junction of the PN-diode (which corresponds to the discontinuity of the doping profile C) crosses some cells. On Figure 2.7, we show the profile of the density of holes $P_M$ computed

<figure><img src="image_31.png" /><figcaption>Figure 2.7 – Test-case 3. Evolution of the discrete density of holes, on a tilted hexagonal mesh.</figcaption></figure>

at different times. These profiles are to be compared with these of the Figure 2.5: even if the tilted hexagonal does not fit the geometry of the problem perfectly and is much coarser than the triangular mesh used for the Figure 2.5, the profiles look quite the same.

We give quantitative information for this test case on Figure 2.8, which is the counterpart of Figure 2.6: we show the evolution of the bounds of the discrete densities, as well as the time step and the number of Newton’s iterations used for a given time step. First, note that the first admissible time step ($\delta t \approx 10^{-6}$) is bigger than the one on the refined triangular mesh. Moreover, the extremal values are much farther from the bounds of $I_h$ than on the triangular mesh. This behaviour can be explained by the fineness of the triangular mesh, which enable the scheme to capture in a very accurate way the boundary layers: the continuous solution take very small/big values near to the boundary of the domain, but if one average these values on relatively big cells (such as the cells of the tilted hexagonal mesh), then the mean values are not so extreme. We can also remark that the number of iterations needed to compute one time step peaks around $t \approx 3.10^{-4}$, as in Figure 2.6: it corresponds to the time when the densities are closest to the bounds of $I_h$. Such fact enforces the hypothesis formulated before: the values of the densities have an important impact on the Jacobian entries and hence on the convergence of the Newton’s method.

Overall, the behaviour of the scheme does not depend on the geometry of the mesh, and the
<figure><img src="image_32.png" /><figcaption>Figure 2.8 – **Test-case 3** ($b = 1$, tilted hexagonal mesh). Evolution of the discrete extremal values, time step and cost.</figcaption></figure>

fact that the PN-junction crosses some cells does not have a noticeable impact.

### 2.5.3 Long-time behaviour of the discrete solutions

We are now interested in the long-time behaviour of the solutions computed with the scheme. For the test-case 4, we consider a test-case from [27] with Boltzmann statistics ($h = \log$), no recombination ($r = 0$), no magnetic field ($b = 0$), $\lambda = 1$ and boundary values $N_0^D = e$, $N_1^D = 1$ and $\alpha_0 = 1$. We are interested in the evolution of the discrete relative entropy and the $L^2$ distance to the equilibrium (namely, $\sqrt{\|N_M^n - N_M^e\|_{L^2(\Omega)}^2 + \|P_M^n - P_M^e\|_{L^2(\Omega)}^2 + \|\phi_M^n - \phi_M^e\|_{L^2(\Omega)}^2}$) with respect to time.

In Figure 2.9, we show the evolution of these quantities for simulations performed on the family of triangular meshes (the coarsest mesh has 56 cells, and a size $h_0$), with $\Delta t = 0.01$. First, one can note that the evolutions are exponentially fast, as expected from Theorem 9. One can also see a saturation phenomenon when the machine precision is reached. Moreover, as announced in Remark 17, the decay rate is not strongly impacted by the meshsize. Last but not least, the quantitative values of the decay rate are in accordance with those obtained in [27, Figure 2] with a Scharfetter–Gummel TPFA scheme.

In Figure 2.10, we show the evolution of the discrete relative entropy and the $L^2$ distance to the equilibrium for simulations computed with $\Delta t = 0.01$ on meshes with different geometry: a Cartesian one (64 cells, $h_D = 3.12 \ 10^{-2}$), a triangular one (56 cells, $h_D = 3.07 \ 10^{-2}$) and a tilted hexagonal one (76 cells, $h_D = 3.42 \ 10^{-2}$). Note that the meshes under consideration have a similar meshsize. As for the previous results, the expected exponential decay of entropy is