<figure><img src="image_33.png" /><figcaption>Figure 2.9 – Test-case 4. Influence of the meshsize: entropy and $L^2$ distance to equilibrium</figcaption></figure>

observed. The decay rates computed are almost the same for the three mesh geometries, which reinforces the previous observations: in practice, the decay rate is not impacted by the mesh used for the simulation (either by its geometry or its size). It is also remarkable to note that on the

<figure><img src="image_34.png" /><figcaption>Figure 2.10 – Test-case 5. Influence of the mesh geometry on the long-time behaviour: entropy and $L^2$ distance to equilibrium.</figcaption></figure>

tilted hexagonal mesh, which is not adapted to the geometry of the semiconductor device, the long-time behaviour is essentially similar to the one on meshes with adapted geometry. Again, this indicates the robustness of the scheme with respect to the mesh used.

Last, we investigate the influence of the magnetic field over the long-time behaviour. We consider the test-case 6, with Blakemore statistics ($h(s) = \log\left(\frac{s}{1-\gamma s}\right)$), no recombination ($r = 0$), $\lambda = 1$ and boundary values $N_0^D = e$, $N_1^D = 1$ and $\alpha_0 = 1$. We perform our simulations on a Cartesian mesh, with a time step $\Delta t = 0.01$, with different values of $b$. The results are presented
<figure><img src="image_35.png" /><figcaption>Figure 2.11 – **Test-case 6.** Long-time behaviour for the Blakemore statistics and influence of the magnetic field</figcaption></figure>

in Figure 2.11a. One can see the relative entropy decreases exponentially fast in time. Moreover, the presence of a magnetic field tends to attenuate the dissipative effects, and to slow down the evolution. This was expected from a physical point of view, since the magnetic field induces a rotation of the charge carriers which delays the natural relaxation towards equilibrium. On Figure 2.11b, we show the evolution of the entropy decay rate with respect to the intensity of the magnetic field *b*. We plot the values of the decay rate for different values of *b*, as well as a reference rate $v(b) = \frac{v_0}{\sqrt{1+b^2}}$, where $v_0$ is the numerical value of the decay rate computed with no magnetic field ($b=0$). The computed values seem to fit very well the reference rate. Note that the decay rate seems therefore to be proportional to $\frac{1}{\sqrt{1+b^2}}$, which happens to be the modulus of the eigenvalues of the diffusion tensors $\Lambda_N$ and $\Lambda_p$.

## 2.6 Conclusion

In this article, we have designed and analysed a scheme for general anisotropic drift-diffusion systems on general polytopal meshes. The scheme is based on the Hybrid Finite Volume method. We have proved that the scheme has a discrete entropic structure, and used it to show the existence of solutions (with $I_h$-valued densities). As a by-product of this structure, we have also proved that the solutions to the scheme converge exponentially fast in time towards the associated discrete thermal equilibrium. The results are established for general statistics functions $h$, general diffusion tensors (potentially anisotropic and nonsymmetric) and general recombination terms, under the hypothesis that the boundary data are compatible with the thermal equilibrium (conditions (2.4) and (2.5)). Finally, we have validated our scheme on different numerical tests, highlighting the $I_h$-valuation of the discrete densities, the ability to withstand intense magnetic fields and the long-time behaviour. Numerical experiments (not presented in this paper) suggest that the scheme introduced here also works in situations where the compatibility condition with the thermal equilibrium does not hold. Hence, a future
direction would be to analyse the scheme in such a situation.

## 2.A Discrete boundedness by entropy and dissipation

In this appendix, we prove Lemma 7. To do so, we need a technical result stated below.

**Lemma 8.** Let $E: (x,y) \mapsto (G(x)-G(y))(x-y)$, and $(x,y) \in \mathbb{R}^2$.

If there exist two positive constants $M$ and $C_1$ such that $|y| \le M$ and $0 \le E(x,y) \le C_1$, then, there exists a constant $C_{M,C_1}$ only depending on $M$ and $C_1$ such that $|x| \le C_{M,C_1}$.

*Proof.* First, we define $E_y: \delta \mapsto (G(y+\delta)-G(y))\delta$, therefore letting $\delta = x-y$, one has $E(x,y) = E_y(\delta) \le C_1$. The inequality $|x| \le M + |\delta|$ holds, so it suffices to get a bound on $\delta$ to conclude:

(i) if $\delta \ge 2M+1 > 0$, then since $G$ is increasing, one has

$$G(y+\delta) \ge G(y+2M+1) \ge G(-M+2M+1) = G(M+1) \text{ and } G(y) \le G(M),$$

so $C_1 \ge E_y(\delta) \ge (G(M+1)-G(M))\delta$, and we get $2M+1 \le \delta \le C_1 (G(M+1)-G(M))^{-1}$;

(ii) if $\delta \le -2M-1$, a similar computation shows that $-2M-1 \ge \delta \ge C_1 (G(-M-1)-G(-M))^{-1}$.

Hence, one has $|\delta| \le \max(2M+1, \frac{C_1}{G(M+1)-G(M)}, \frac{C_1}{G(-M)-G(-M-1)})$, which concludes the proof. □

We can now prove Lemma 7. The proof is similar to [70, Lemma 2], except for the use of the unknowns on $\mathcal{E}_{ext}^D$.

**Lemma 7.** Let $(\underline{w}_D^N, \underline{w}_D^P) \in \mathcal{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_\# \qquad (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)$.

*Proof.* The proof is divided into two steps. First, we prove bounds on the discrete electrostatic potential $\underline{\phi}_D$ thanks to the bound on the entropy. Then, we use these bounds to estimate $(\underline{w}_D^N, \underline{w}_D^P)$.

Since $\mathbb{E}(\underline{w}_D^N, \underline{w}_D^P) \le B_\#$, by (2.16) one has $\frac{\alpha_b \lambda_b^\phi}{2} |\underline{\phi}_D - \underline{\phi}_D^e|_{1,D}^2 \le \frac{1}{2} a_D^\phi (\underline{\phi}_D - \underline{\phi}_D^e, \underline{\phi}_D - \underline{\phi}_D^e) \le B_\#$. Therefore, letting $c_1 = |\underline{\phi}_D^e|_{1,D} + \sqrt{\frac{2B_\#}{\alpha_b \lambda_b^\phi}}$, we get that $|\underline{\phi}_D|_{1,D} \le c_1$. By definition of $|\cdot|_{1,D}$, we deduce that, for any $K \in \mathcal{M}$ and any $\sigma \in \mathcal{E}_K$, $\frac{|\sigma|}{d_{K,\sigma}} (\phi_K - \phi_\sigma)^2 \le c_1^2$. Letting $c_2 = \max_{K \in \mathcal{M}, \sigma \in \mathcal{E}_K} c_1 \sqrt{\frac{d_{K,\sigma}}{|\sigma|}}$, one has $\forall K \in \mathcal{M}, \forall \sigma \in \mathcal{E}_K, |\phi_K - \phi_\sigma| \le c_2$.
On the other hand, there exists $\sigma_0 \in \mathcal{E}_{\text{ext}}^D$ such that $\phi_{\sigma_0} = \phi_{\sigma_0}^D$ and so, by (2.11), we have

$$|\phi_{\sigma_0}| \le \|\phi^D\|_{L^\infty(\Omega)}.$$

Now, one can use the connectivity of the mesh: for any cell $K$ (resp. face $\sigma$) there is a finite sequence of components of $\phi_D$, denoted $(x_k)_{0 \le k \le l}$, starting at $x_0 = \phi_{\sigma_0}$ and finishing at $x_l = \phi_K$ (resp. $x_l = \phi_\sigma$) such that, for any $k$ in $\{0, \dots, l-1\}$, one has $|x_{k+1} - x_k| \le c_2$. Therefore, one concludes that

$$|x_l| \le lc_2 + |x_0| \le 2|\mathcal{M}|c_2 + \|\phi^D\|_{L^\infty(\Omega)}.$$

Thus there exists $\phi_\#$ positive depending on the data and $\mathcal{D}$ such that

$$-\phi_\# \mathbb{1}_D \le \underline{\phi}_D \le \phi_\# \mathbb{1}_D.$$

Now, note that since $\mathbb{D}(\underline{w}_D^N, \underline{w}_D^P) \le B_\#$, one has $T_D^N(\underline{N}_D, \underline{w}_D^N, \underline{w}_D^N) \le B_\#$ and $T_D^P(\underline{P}_D, \underline{w}_D^P, \underline{w}_D^P) \le B_\#$, where we recall that the discrete densities are defined by $\underline{N}_D = g(\underline{w}_D^N + \underline{\phi}_D + \alpha_N \mathbb{1}_D)$ and $\underline{P}_D = g(\underline{w}_D^P - \underline{\phi}_D + \alpha_P \mathbb{1}_D)$. For $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

$$B_\# \ge r_K(\underline{N}_K) a_K^N(\underline{w}_D^N, \underline{w}_D^N) \ge \alpha_b \lambda_b \sum_{\sigma \in \mathcal{E}_K} r_K(\underline{N}_K) \frac{|\sigma|}{d_{K,\sigma}} (w_K^N - w_\sigma^N)^2. \quad (2.68)$$

Given $\sigma \in \mathcal{E}_K$, by definition of $r_K$ and $\underline{N}_K$, using (2.24) and the positivity of $m$, we obtain

$$r_K(\underline{N}_K) = \frac{1}{|\mathcal{E}_K|} \sum_{\sigma' \in \mathcal{E}_K} m(N_K, N_{\sigma'}) \ge \frac{1}{|\mathcal{E}_K|} m_h \left( g(w_K^N + \phi_K + \alpha_N), g(w_\sigma^N + \phi_\sigma + \alpha_N) \right).$$

Recall that formula (2.26) holds, so

$$r_K(\underline{N}_K) \ge \frac{1}{|\mathcal{E}_K|} \frac{G(w_K^N + \phi_K + \alpha_N) - G(w_\sigma^N + \phi_\sigma + \alpha_N)}{(w_K^N + \phi_K + \alpha_N) - (w_\sigma^N + \phi_\sigma + \alpha_N)}.$$

Moreover, since $G$ is convex, $(x, y) \mapsto \frac{G(x)-G(y)}{x-y}$ is non-decreasing with respect to both its variables, therefore using the bound on $\underline{\phi}_D$ proved previously we get

$$r_K(\underline{N}_K) \ge \frac{1}{|\mathcal{E}_K|} \frac{G(w_K^N - \phi_\# + \alpha_N) - G(w_\sigma^N - \phi_\# + \alpha_N)}{(w_K^N - \phi_\# + \alpha_N) - (w_\sigma^N - \phi_\# + \alpha_N)} \\ = \frac{1}{|\mathcal{E}_K|} \frac{G(w_K^N - \phi_\# + \alpha_N) - G(w_\sigma^N - \phi_\# + \alpha_N)}{w_K^N - w_\sigma^N}.$$

Using (2.68), one gets that for any $K \in \mathcal{M}$ and for any $\sigma \in \mathcal{E}_K$,

$$0 \le (G(w_K^N - \phi_\# + \alpha_N) - G(w_\sigma^N - \phi_\# + \alpha_N))(w_K^N - w_\sigma^N) \le \zeta$$

where $\zeta = \frac{B_\#}{\alpha_b \lambda_b} \max_{K \in \mathcal{M}, \sigma \in \mathcal{E}_K} \frac{d_{K,\sigma}|\mathcal{E}_K|}{|\sigma|}$. Now, let $\underline{u}_D = \underline{w}_D^N + (-\phi_\# + \alpha_N)\mathbb{1}_D \in \underline{V}_D$. It is clear that it
suffices to get bounds on $\underline{u}_D$ to bound $\underline{w}_D^N$. The previous relation writes

$$ \forall K \in \mathcal{M}, \forall \sigma \in \mathcal{E}_K, 0 \le E(u_K, u_\sigma) \le \zeta, \qquad (2.69) $$

where $E: \mathbb{R}^2 \to \mathbb{R}_+$ is the function defined in Lemma 8. Now, notice that since $\underline{w}_D^N \in \underline{V}_{D,0}$ there exists $\sigma_0 \in \mathcal{E}_{ext}^D$ such that $u_{\sigma_0} = -\phi_\# + \alpha_N$. One can then use the connectivity of the mesh to conclude: for any cell $K$ (resp. face $\sigma$) there is a finite sequence of components of $\underline{u}_D$, denoted $(x_k)_{0 \le k \le l}$, starting at $x_0 = u_{\sigma_0}$ and finishing at $x_l = u_K$ (resp. $x_l = u_\sigma$) such that for any $k$ in $\{0, \dots, l-1\}$, one has $E(x_{k+1}, x_k) \le \zeta$. Therefore, by $l$ successive applications of Lemma 8, we get the existence of some $c^N > 0$, depending on $\zeta$, $l$, and $x_0 = -\phi_\# + \alpha_N$ such that $|x_l| \le c^N$. Thus, there exists $C_\#^N > 0$ depending on the data, $\mathcal{D}$ and $B_\#^N$ such that

$$ \forall K \in \mathcal{M}, \forall \sigma \in \mathcal{E}, \quad |w_K^N| \le C_\#^N \quad \text{and} \quad |w_\sigma^N| \le C_\#^N. $$

Using the same strategy, one gets bounds for $\underline{w}_D^P$, which concludes the proof. $\square$

<figure><img src="image_36.png" /></figure>
Chapter 3

# A comparison of structure-preserving schemes for drift-diffusion systems on general meshes: DDFV vs HFV

## Outline of the current chapter

<table><tr><td>3.1 Motivation</td><td>118</td></tr><tr><td>3.2 Description of the schemes</td><td>119</td></tr><tr><td>3.2.1 The DDFV scheme</td><td>120</td></tr><tr><td>3.2.2 The HFV scheme</td><td>121</td></tr><tr><td>3.2.3 Some structural differences between schemes</td><td>123</td></tr><tr><td>3.3 Numerical experiments</td><td>123</td></tr><tr><td>3.3.1 Description of the test case</td><td>123</td></tr><tr><td>3.3.2 Positivity</td><td>124</td></tr><tr><td>3.3.3 Long-time behaviour</td><td>125</td></tr><tr><td>3.4 Conclusion</td><td>126</td></tr></table>

This chapter is a work in collaboration with Stella Krell, accepted as a proceeding of the FVCA X conference [187]. For the purpose of unified presentation and comparison, the notations used in this chapter are slightly different from those of the two previous chapters.

We made a comparison between a Discrete Duality Finite Volume (DDFV) scheme and a Hybrid Finite Volume (HFV) scheme for a drift-diffusion model with mixed boundary conditions on general meshes. Both schemes are based on a nonlinear discretisation of the convection-diffusion fluxes, which ensures the positivity of the discrete densities. We investigate the behaviours of the schemes on various numerical test cases.
## 3.1 Motivation

We are interested in the numerical discretization of a drift-diffusion model on general meshes. Let $\Omega$ be a polygonal connected open bounded subset of $\mathbb{R}^2$, whose boundary $\Gamma = \partial\Omega$ is divided into two parts $\Gamma = \Gamma^D \cup \Gamma^N$ with $m(\Gamma^D) > 0$. The problem writes:

$$ egin{cases} 
\partial_t N - \operatorname{div}(\nabla N - N \nabla \phi) = 0 & \text{in } \mathbb{R}_+ \times \Omega, \\
\partial_t P - \operatorname{div}(\nabla P + P \nabla \phi) = 0 & \text{in } \mathbb{R}_+ \times \Omega, \\
-\lambda^2 \operatorname{div}(\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, \\
(\nabla N - N \nabla \phi) \cdot n = (\nabla P + P \nabla \phi) \cdot n = \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 (3.1) $$

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

(i) the parameter $\lambda > 0$ is the rescaled Debye length of the system, which accounts for the nondimensionalisation (relevant values of this parameter can be very small, inducing some stiff behaviours),

(ii) the initial conditions $N^{in}$ and $P^{in}$ belong to $L^\infty(\Omega)$ and are positive,

(iii) the doping profile $C$ is in $L^\infty(\Omega)$, and characterises the semiconductor device used.

In the following, we also assume that the boundary conditions are the trace of some $H^1$ function on $\Omega$, such that the following relation holds:

$$ \log(N^D) - \phi^D = \alpha_N \text{ and } \log(P^D) + \phi^D = \alpha_P \text{ on } \Gamma^D, \quad (3.2) $$

where $\alpha_N$ and $\alpha_P$ are two real constants. It follows that $N^D$ and $P^D$ are positive.

The solution to (3.1) enjoys some natural physical properties: the densities $N$ and $P$ are positive for all time, and the solution converges exponentially fast towards some thermal equilibrium $(N^e, P^e, \phi^e)$ as the time goes to infinity. This thermal equilibrium, which is a stationary solution to (3.1), is defined such that

$$ N^e = e^{\alpha_N + \phi^e} \quad \text{and} \quad P^e = e^{\alpha_P - \phi^e}, \quad (3.3) $$

where the electrostatic potential at the equilibrium $\phi^e$ is the solution to the Poisson-Boltzmann equation

$$ egin{cases} 
-\lambda^2 \operatorname{div}(\nabla \phi^e) = C + \exp(\alpha_P - \phi^e) - \exp(\alpha_N + \phi^e) & \text{in } \Omega, \\
\phi^e = \phi^D & \text{on } \Gamma^D, \\
\nabla \phi^e \cdot n = 0 & \text{on } \Gamma^N. 
\end{cases} 
\quad (3.4) $$

The relation (3.2) is a compatibility condition in order to ensure that the only stationary state of (3.1) is the thermal equilibrium defined by (3.4).

When designing numerical schemes for (3.1), it is crucial to ensure that the scheme preserves these properties at the discrete level. This structure preserving features are ensured by classical Two-Point Flux Approximation (TPFA) schemes on admissible orthogonal meshes (see [27]). Unfortunately, these schemes cannot be used on general meshes, and their use is essentially restricted to isotropic problems. On the other hand, it appears that the use of adaptive mesh refinement could be a relevant solution in order to handle the stiff behaviour of these problems, as presented in [58]. Therefore, a first step towards this kind of strategy lies in the development
of structure preserving scheme on general meshes. However, the main issue of schemes handling general meshes for diffusive problems is the lack of discrete positivity, both at the theoretical and practical level (see [109]). One of the possible solution in order tackle this problem is to use a nonlinear approach based on the entropy structure of the model. Such idea was originally introduced in [56] with a Vertex Approximate Gradient (VAG) scheme for advection-diffusion problems. Following the ideas of this seminal work, a nonlinear positivity preserving Discrete Duality Finite Volume (DDFV) scheme for linear advection-diffusion equations has been introduced in [52]. The good long-time behaviour of this scheme has then been established in [51]. In the spirit of these works, a nonlinear structure preserving Hybride Finite Volume (HFV) scheme was introduced and analysed in [213]: this scheme is proved to preserve the positivity of the densities as well as the long-time behaviour of the solutions. These three schemes are based on different spatial discretisation methods handling general meshes (VAG, DDFV and HFV), but they share the same spirit in terms of conception and analysis.

The aim of this paper is therefore to compare two of these methods for a realistic coupled problem, from a numerical point of view. To do so, we to introduce a nonlinear structure preserving DDFV scheme for (3.1) based on the scheme of [52] and we compare it numerically with the HFV scheme of [213].

## 3.2 Description of the schemes

The schemes used here are based on the same nonlinear strategy, introduced in [56], consisting in the reformulation of the convection-diffusion fluxes:

$$ \nabla N - N\nabla\phi = N\nabla(\log(N) - \phi) \text{ and } \nabla P + P\nabla\phi = P\nabla(\log(P) + \phi). $$

The quantities $\log(N) - \phi$ and $\log(P) + \phi$ are often called quasi-Fermi potentials in the literature associated to semiconductor models, and have an important physical meaning. At the discrete level, both schemes rely on discrete gradients operators to approximate the continuous gradients. The major issue lies in the discretisation of the prefactors $P$ and $N$, which will be handled by local reconstruction operators. This discretisation strategy is a way of ensuring, both at the theoretical and practical levels, the positivity of the discrete densities. We refer to [213, Theorem 1] (HFV scheme for drift-diffusion system) and [52, Theorem 2.1] (DDFV scheme for a single advection-diffusion equation) for proofs of this statement. We also refer the reader to these proofs for more insight about the reconstruction operators.

Both schemes are based on a backward Euler discretisation in time. To fix ideas, we will use a constant time step $\Delta t > 0$. For more precise descriptions and statements about the schemes and the meshes, we refer to [52] (DDFV) and [213] (HFV).

**Remark 18 (Generalisation to anisotropic models).** In this paper, we consider isotropic convection-diffusion equations for the charges carriers for the sake of brevity. One could add anisotropic diffusion tensors and consider the framework described in [213]. Thus, the presentation made here can easily be extended to models with external magnetic fields (see [143]).

Both schemes rely on a spatial discretisation (or mesh) of the domain $\Omega$. The (primal interior) mesh $\mathcal{M}$ is a partition of $\Omega$ in polygonal control volumes (or cells). We let $\partial\mathcal{M}$ be the set of boundary edges, seen either as degenerate control volumes (DDFV framework) or as edges (HFV framework). The primal mesh $\overline{\mathcal{M}}$ is defined as the reunion of $\mathcal{M}$ and $\partial\mathcal{M}$. Given a cell $K \in \overline{\mathcal{M}}$, we fix a point $x_K \in K$, called the center of $K$. For all neighboring primal cells $K$ and $L$, we assume that $\partial K \cap \partial L$ is a segment, corresponding to an internal edge of the mesh $\mathcal{M}$, denoted by $\sigma = K|L$ and we let $\mathcal{E}_{int}$ be the set of such edges. We denotes by $\mathcal{E} = \mathcal{E}_{int} \cup \partial\mathcal{M}$ the set of all (internal and
exterior) edges of the mesh, and define $\mathcal{E}_K$ the set of edges of the cell $K \in \mathbb{N}$. For any $K \in \mathbb{N}$ and $\sigma \in \mathcal{E}_K$, we define $\mathbf{n}_{\sigma K}$ as the unit normal to $\sigma$ outward $K$.

Given any measurable $X \subset \mathbb{R}^2$, we denote by $m_X$ the measure of the object X.

### 3.2.1 The DDFV scheme

In order to define the DDFV scheme, we need to introduce two other meshes: the dual mesh denoted $\overline{\mathbb{N}}^*$ and the diamond mesh denoted $\mathcal{D}$ (see [52] for more details). The dual mesh $\overline{\mathbb{N}}^*$ is also composed of interior dual mesh $\mathbb{N}^*$ (corresponding of cells around vertex in $\Omega$) and of boundary dual mesh $\partial\mathbb{N}^*$ (corresponding of cells around vertex on $\partial\Omega$). For any vertex $x_{K^*}$ of the primal mesh satisfying $x_{K^*} \in \Omega$, we define a polygonal control volume $K^*$ by connecting all the centers of the primal cells sharing $x_{K^*}$ as vertex. For any vertex $x_{K^*} \in \partial\Omega$, we define a polygonal control volume $K^*$ by connecting the centers $x_K$ of the interior primal cells and the midpoints of the boundary edges sharing $x_{K^*}$ as vertex and $x_{K^*}$. We define the set $\mathcal{E}_{int}^*$ of internal edges of the dual mesh similarly as $\mathcal{E}_{int}$. We denote by $\mathbf{n}_{\sigma^* K^*}$ the unit normal to $\sigma^*$ outward $K^*$. For each couple $(\sigma, \sigma^*) \in \mathcal{E} \times \mathcal{E}_{int}^*$ such that $\sigma = [x_{K^*}, x_{L^*}]$ and $\sigma^* = K^*|L^*$, we define the quadrilateral diamond $\mathcal{D}_{\sigma, \sigma^*}$ whose diagonals are $\sigma$ and $\sigma^*$ (if $\sigma \subset \partial\Omega$, it degenerates into a triangle). The set of the diamonds defines the diamond mesh $\mathcal{D}$, which is a partition of $\Omega$. Finally, the DDFV

<figure><img src="image_37.png" /><figcaption>Figure 3.1 – Definition of the diamonds $\mathcal{D}_{\sigma, \sigma^*}$ and related notations.</figcaption></figure>

mesh is made of $\mathcal{T} = (\overline{\mathbb{N}}, \overline{\mathbb{N}}^*)$ and $\mathcal{D}$.

We now introduce the space of scalar fields which are associated to each primal and dual cell $\mathbb{R}^T$, and space of vector fields constant on the diamonds $(\mathbb{R}^2)^\mathcal{D}$:

$$u_T \in \mathbb{R}^T \iff u_T = (u_K)_{K \in \overline{\mathbb{N}}}, (u_{K^*})_{K^* \in \overline{\mathbb{N}}^*} \text{ and } \xi_D \in (\mathbb{R}^2)^\mathcal{D} \iff \xi_D = (\xi_D)_{D \in \mathcal{D}}.$$

To enforce Dirichlet boundary conditions, we introduce the set of Dirichlet boundary primal and dual cells: $\partial\mathbb{N}_D = \{K \in \partial\mathbb{N} : K \subset \Gamma_D\}$ and $\partial\mathbb{N}_D^* = \{K^* \in \partial\mathbb{N}^* : x_{K^*} \in \overline{\Gamma}_D\}$, and, for a given $v \in C(\Gamma^D)$, we define

$$E_v^{\Gamma_D} = \{u_T \in \mathbb{R}^T \mid \forall K \in \partial\mathbb{N}_D, u_K = v(x_K) \text{ and } \forall K^* \in \partial\mathbb{N}_D^*, u_{K^*} = v(x_{K^*})\}.$$

These sets can be generalised to less regular boundary data $v$ by substituting punctual evaluation
by averages. We also define discrete bilinear forms on $\mathbb{R}^T$ and $(\mathbb{R}^2)^D$ by

$$ \|v_T, u_T\|_T = \frac{1}{2} \left( \sum_{K \in \mathfrak{M}} m_K u_K v_K + \sum_{K^* \in \mathfrak{M}^*} m_{K^*} u_{K^*} v_{K^*} \right), \quad \forall (u_T, v_T) \in (\mathbb{R}^T)^2, $$

$$ (\xi_D, \varphi_D)_D = \sum_{D \in \mathcal{D}} m_D \xi_D \cdot \varphi_D, \quad \forall (\xi_D, \varphi_D) \in ((\mathbb{R}^2)^D)^2. $$

The DDFV method is based on the definition of a discrete gradient operator $\nabla^D : \mathbb{R}^T \to (\mathbb{R}^2)^D$, defined by $\nabla^D u_T = (\nabla^D u_T)_{D \in \mathcal{D}}$, where

$$ \nabla^D u_T = \frac{1}{2m_D} (m_\sigma(u_L - u_K)\mathbf{n}_{\sigma K} + m_{\sigma^*}(u_{L^*} - u_{K^*})\mathbf{n}_{\sigma^* K^*}) \quad \forall D \in \mathcal{D}. \quad (3.5) $$

Finally, we introduce a reconstruction operator on diamonds $r^D$. It is a mapping from $\mathbb{R}^T$ to $\mathbb{R}^D$ defined for all $u_T \in \mathbb{R}^T$ by $r^D u_T = (r^D u_T)_{D \in \mathcal{D}}$, where for $D \in \mathcal{D}$, whose vertices are $x_K$, $x_L$, $x_{K^*}$, $x_{L^*}$, $r^D u_T = \frac{1}{4}(u_K + u_L + u_{K^*} + u_{L^*})$. One can now introduce a DDFV discretisation of $(u, w, v) \mapsto \int_{\Omega} u \nabla w \cdot \nabla v$, defined by

$$ T_D : (u_T, w_T, v_T) \mapsto \sum_{D \in \mathcal{D}} m_D r^D u_T \nabla^D w_T \cdot \nabla^D v_T. $$

Now, we first discretise the data by taking the mean values of $N^{in}$, $P^{in}$ and C on the primal and dual cells, which define $N_T^0$, $P_T^0$ and $C_T$. Then, for all $n \ge 0$, we look for $(N_T^{n+1}, P_T^{n+1}, \phi_T^{n+1}) \in E_{N^D}^{\Gamma_D} \times E_{P^D}^{\Gamma_D} \times E_{\phi^D}^{\Gamma_D}$ solution to:

$$ \left\| \frac{N_T^{n+1} - N_T^n}{\Delta t}, v_T \right\|_T + T_D(N_T^{n+1}, \log(N_T^{n+1}) - \phi_T^{n+1}, v_T) = 0 \quad \forall v_T \in E_0^{\Gamma_D}, \quad (3.6a) $$

$$ \left\| \frac{P_T^{n+1} - P_T^n}{\Delta t}, v_T \right\|_T + T_D(P_T^{n+1}, \log(P_T^{n+1}) + \phi_T^{n+1}, v_T) = 0 \quad \forall v_T \in E_0^{\Gamma_D}, \quad (3.6b) $$

$$ \lambda^2 (\nabla^D \phi_T^{n+1}, \nabla^D v_T)_D = \left\| C_T + P_T^{n+1} - N_T^{n+1}, v_T \right\|_T \quad \forall v_T \in E_0^{\Gamma_D}. \quad (3.6c) $$

In (3.6a) and (3.6b), we use the notation $\log(u_T) = ((\log(u_K))_{K \in \mathfrak{M}}, (\log(u_{K^*}))_{K^* \in \mathfrak{M}^*})$.

### 3.2.2 The HFV scheme

In order to define the HFV scheme, we need to introduce a pyramidal submesh. To do so, one has to assume that each cell $K \in \mathfrak{M}$ is star-shaped with respect to its center $x_K$ (we recall that $x_K$ is not necessarily the barycentre of $K$). We then define $P_{K,\sigma}$ as the pyramid (triangle) of base $\sigma$ and apex $x_K$. Given any $\sigma \in \mathcal{E}$, we denote by $\bar{x}_\sigma$ the barycentre of $\sigma$, and by $d_{K,\sigma}$ the euclidean distance between $\sigma$ and $x_K$. Finally, we define the hybrid discretisation (or mesh) as $\mathcal{D} = (\mathfrak{M}, \mathcal{E})$. We now introduce the space of discrete (scalar) hybrid unknowns $\underline{u}_\mathcal{D}$:

$$ \underline{u}_{\mathcal{D}} \in \underline{V}_{\mathcal{D}} \iff \underline{u}_{\mathcal{D}} = ((u_K)_{K \in \mathfrak{M}}, (u_\sigma)_{\sigma \in \mathcal{E}}), $$

where the $u_K \in \mathbb{R}$ are the cell unknowns and the $u_\sigma \in \mathbb{R}$ are the edges unknowns (approximation
of the trace of the solutions on the edges). To enforce Dirichlet boundary conditions, for a given $v \in C(\Gamma^D)$, we define

$$ \underline{V}_{\mathcal{D},v}^{\Gamma_D} = \{ \underline{u}_{\mathcal{D}} \in \underline{V}_{\mathcal{D}} \mid \forall \sigma \in \partial \mathfrak{M}_{\mathcal{D}}, u_{\sigma} = v(\bar{x}_{\sigma}) \}. $$

As for the DDFV framework, we define a bilinear form on $\underline{V}_{\mathcal{D}}$, discrete counterpart of the inner product on $L^2(\Omega)$ as

$$ \llbracket \underline{u}_{\mathcal{D}}, \underline{v}_{\mathcal{D}} \rrbracket_{\mathfrak{M}} = \sum_{K \in \mathfrak{M}} m_K u_K v_K, \quad \forall (\underline{u}_{\mathcal{D}}, \underline{v}_{\mathcal{D}}) \in \underline{V}_{\mathcal{D}}^2. $$

Note that this discrete inner product only takes into account the values of the cell unknowns.

The HFV method is based on the definition of a discrete gradient operator $\nabla_{\mathcal{D}}: \underline{V}_{\mathcal{D}} \to (\mathbb{R}^2)^{\Omega}$ which maps discrete hybrid unknowns onto piecewise constant functions on the pyramidal submesh. More precisely, given $\underline{v}_{\mathcal{D}} \in \underline{V}_{\mathcal{D}}$, $K \in \mathfrak{M}$ and $\sigma \in \mathcal{E}_K$ we let

$$ \nabla_{\mathcal{D}} \underline{v}_{\mathcal{D}}|_{P_{K,\sigma}} = G_K \underline{v}_{\mathcal{D}} + S_{K,\sigma} \underline{v}_{\mathcal{D}}, $$

where the consistent part of the gradient $G_K \underline{v}_{\mathcal{D}}$ and the stabilisation part $S_{K,\sigma} \underline{v}_{\mathcal{D}}$ are defined as

$$ G_K \underline{v}_{\mathcal{D}} = \frac{1}{m_K} \sum_{\sigma' \in \mathcal{E}_K} m_{\sigma'} v_{\sigma'} n_{K,\sigma'} \quad \text{and} \quad S_{K,\sigma} \underline{v}_{\mathcal{D}} = \frac{\eta}{d_{K,\sigma}} (v_{\sigma} - v_K - G_K \underline{v}_K \cdot (\bar{x}_{\sigma} - x_K)) n_{K,\sigma}, $$

where $\eta$ is a given positive stabilisation parameter. One can now define the discrete counterpart of $(u,v) \mapsto \int_{\Omega} \nabla u \cdot \nabla v$ as

$$ a_{\mathcal{D}} : (\underline{u}_{\mathcal{D}}, \underline{v}_{\mathcal{D}}) \mapsto \int_{\Omega} \nabla_{\mathcal{D}} \underline{u}_{\mathcal{D}} \cdot \nabla_{\mathcal{D}} \underline{v}_{\mathcal{D}}. $$

Finally, we introduce as previously local reconstruction operators on cells $r^K: \underline{V}_{\mathcal{D}} \to \mathbb{R}$, such that for any $u_{\mathcal{D}} \in \underline{V}_{\mathcal{D}}$, $r^K(\underline{u}_{\mathcal{D}}) = \frac{1}{|\mathcal{E}_K|} \sum_{\sigma \in \mathcal{E}_K} \frac{u_K + u_{\sigma}}{2}$, where $|\mathcal{E}_K|$ is the cardinal of the finite set $\mathcal{E}_K$.

One can now introduce a HFV discretisation of $(u,w,v) \mapsto \int_{\Omega} u \nabla w \cdot \nabla v$, defined by

$$ T_{\mathcal{D}} : (\underline{u}_{\mathcal{D}}, \underline{w}_{\mathcal{D}}, \underline{v}_{\mathcal{D}}) \mapsto \sum_{K \in \mathfrak{M}} r^K(\underline{u}_{\mathcal{D}}) \int_K \nabla_{\mathcal{D}} \underline{w}_{\mathcal{D}} \cdot \nabla_{\mathcal{D}} \underline{v}_{\mathcal{D}}. $$

We now discretise the data by taking the mean values of $N^{in}$, $P^{in}$ and $C$ on the cells and edges, which define $\underline{P}_{\mathcal{D}}^0$, $\underline{N}_{\mathcal{D}}^0$ and $\underline{C}_{\mathcal{D}}$. Then, for all $n \ge 0$, we look for $(\underline{N}_{\mathcal{D}}^{n+1}, \underline{P}_{\mathcal{D}}^{n+1}, \underline{\phi}_{\mathcal{D}}^{n+1}) \in \underline{V}_{\mathcal{D},N^D}^{\Gamma_D} \times \underline{V}_{\mathcal{D},P^D}^{\Gamma_D} \times \underline{V}_{\mathcal{D},\phi^D}^{\Gamma_D}$ solution to:

$$ \left[ \left\| \frac{\underline{N}_{\mathcal{D}}^{n+1} - \underline{N}_{\mathcal{D}}^n}{\Delta t}, \underline{v}_{\mathcal{D}} \right\|_{\mathfrak{M}} + T_{\mathcal{D}}(\underline{N}_{\mathcal{D}}^{n+1}, \log(\underline{N}_{\mathcal{D}}^{n+1}) - \underline{\phi}_{\mathcal{D}}^{n+1}, \underline{v}_{\mathcal{D}}) = 0 \quad \forall \underline{v}_{\mathcal{D}} \in \underline{V}_{\mathcal{D},0}^{\Gamma_D}, \quad (3.7a) $$

$$ \left[ \left\| \frac{\underline{P}_{\mathcal{D}}^{n+1} - \underline{P}_{\mathcal{D}}^n}{\Delta t}, \underline{v}_{\mathcal{D}} \right\|_{\mathfrak{M}} + T_{\mathcal{D}}(\underline{P}_{\mathcal{D}}^{n+1}, \log(\underline{P}_{\mathcal{D}}^{n+1}) + \underline{\phi}_{\mathcal{D}}^{n+1}, \underline{v}_{\mathcal{D}}) = 0 \quad \forall \underline{v}_{\mathcal{D}} \in \underline{V}_{\mathcal{D},0}^{\Gamma_D}, \quad (3.7b) $$

$$ \lambda^2 a_{\mathcal{D}}(\underline{\phi}_{\mathcal{D}}^{n+1}, \underline{v}_{\mathcal{D}}) = \left[ \underline{C}_{\mathcal{D}} + \underline{P}_{\mathcal{D}}^{n+1} - \underline{N}_{\mathcal{D}}^{n+1}, \underline{v}_{\mathcal{D}} \right]_{\mathfrak{M}} \quad \forall \underline{v}_{\mathcal{D}} \in \underline{V}_{\mathcal{D},0}^{\Gamma_D}. \quad (3.7c) $$

As previously, we use the notation $\log(\underline{u}_{\mathcal{D}}) = (\log(u_K))_{K \in \mathfrak{M}}, (\log(u_{\sigma}))_{\sigma \in \mathcal{E}}$.
### 3.2.3 Some structural differences between schemes

As highlighted by the unified presentation above, both schemes are very similar and rely on the same features. Note that both local reconstruction operators $r^D$ and $r^K$ take into account all the local unknowns of the geometric entity considered (diamond or cells), this property is the key point of the analysis of this kind of schemes, see [52, 213].

However, the schemes exhibit differences, some of which are listed below:

- the discrete HFV gradient $\nabla_D$ includes a stabilisation term for the sake of coercivity and the stabilisation parameter $\eta$ has to be chosen a priori, whereas the DDFV one is simpler and do not need any choice of parameter;

- the DDFV unknowns are all "volumic", in the sense that there are associated to geometric entities with non-zero two-dimensional measures, whereas the faces unknowns of the HFV method have no mass and have no influence on the discrete time derivative terms $\llbracket \underline{N}_D^{n+1} - \underline{N}_D^n, \underline{v}_D \rrbracket_M$ and $\llbracket \underline{P}_D^{n+1} - \underline{P}_D^n, \underline{v}_D \rrbracket_M$;

- the cells unknowns of the HFV scheme can be eliminated before solving linear systems, using a static condensation procedure (see [213, Section 5.1.2.]), whereas one has to solve a system including all primal and dual unknowns for DDFV;

- the HFV scheme can be used in 3D without any modification (the edges become faces), whereas using a DDFV method in 3D requires more sophisticated changes (see [90]).

## 3.3 Numerical experiments

The two numerical schemes described here are nonlinear, hence their algebraic realisations boil down to the resolution of nonlinear systems of equations. To solve these systems, we use Newton method, with an adaptative time stepping strategies: if the Newton method does not converge, we try to compute the solution for a smaller time step $0.5 \times \Delta t$. If the method converges, we use a bigger time step $1.4 \times \Delta t$. We fix the initial time step as well as the maximal time step allowed, and denote these quantities by $\Delta t_{\text{init}}$ and $\Delta t_{\text{max}}$.

We use the same stopping criterion for both scheme, namely the size of the $l^2$ norm of the goal function (with the same threshold). Note that even if we use the same algebraic criterion, this does not mean in practice that the criterion has the same level of exigence for both schemes. In particular, since the HFV scheme has two different type of unknowns (cell and face), understanding how to create equivalent stopping criteria for the two schemes is unclear at this point. The results presented below indicate that, indeed, the criterion used leads to slightly different behaviours.

For the HFV scheme, at each system resolution, a static condensation is used to eliminate the cell unknowns (see [213, Section 5.1.2.]), and we use $\eta = 1.5$.

Note that we use $N$, $P$ and $\phi$ as discrete unknowns in the schemes. One could also solve the system using the quasi-Fermi potentials and the electrostatic potential as discrete unknowns. Such a choice corresponds at the level of the implementation to a preconditioning. It could maybe leads to better performances in practice, but it seems difficult to find relevant criterion in order to know when switch between the two sets of unknowns.

#### 3.3.1 Description of the test case

The test cases used below follow the framework used in [213] to describe a 2D PN-junction, whose geometry is described in Figure 3.2. The domain $\Omega$ is the unit square $[0, 1]^2$. For the
<figure><img src="image_38.png" /><figcaption>Figure 3.2 – PN diode geometry.</figcaption></figure>

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 (3.2) we assume that there exists a constant $\alpha_0$ such that $\log(N^D \times P^D) = \alpha_0$. Therefore for given $N^D$ and $\alpha_0$ we set

$$P^D = \frac{e^{\alpha_0}}{N^D} \text{ on } \Gamma^D.$$

Thus, one has $\alpha_N = \alpha_P = \frac{\alpha_0}{2}$. The doping profile C is piecewise constant, equal to -1 in the P-region and 1 in the N-region (see Figure 3.2). Last, we use the following smooth initial conditions:

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

### 3.3.2 Positivity

In this section, we compare the discrete positivity preservation of the schemes. The test used here corresponds to the following values:

$$\lambda = 0.05, \quad N_0^D = 0.1, \quad N_1^D = 1 \quad \text{and} \quad \alpha_0 = -4.$$

We perform a test on a distorted quadrangle mesh (mesh_quad_6 of the FVCA 8 Benchmark), with $\Delta t_{ini} = 1.410^{-3}$ and $\Delta t_{max} = 0.1$. We show in Figure 3.3 the evolution of the minimal values of P and N, along with the time step and the number of Newton's iterations needed to compute the solutions at a given time for each time step. The minimal values are taken on every unknowns (primal and dual cells for the DDFV scheme, cells and faces for the HFV one). One can see that both schemes compute, as expected by the theoretical results, positive densities. The minimal values computed are of the same order for both schemes. Moreover, both computations proceed without the need of a time step reduction. Regarding the cost, it appears that the HFV scheme
<figure><img src="image_39.png" /><figcaption>Figure 3.3 – Positivity. Evolution of the discrete minimal values, time step and cost</figcaption></figure>

needs more Newton iterations than the DDFV one (90 vs 63). For both schemes, the number of iteration decay as the time increases, since the solutions converge exponentially fast towards the equilibrium.

### 3.3.3 Long-time behaviour

Here, we investigate the long-time behaviour of the schemes. At the continuous level, one usually quantify the distance between the solution $(N, P, \phi)$ and the equilibrium $(N^e, P^e, \phi^e)$ by looking at the relative entropy, defined as

$$ \mathbb{E}(t) = \int_{\Omega} N^e H\left(\frac{N}{N^e}\right) + \int_{\Omega} P^e H\left(\frac{P}{P^e}\right) + \frac{\lambda^2}{2} \|\nabla(\phi - \phi^e)\|_{L^2(\Omega)}^2, $$

with $H: s \mapsto s \log(s) - s + 1$. One can check that $(N, P, \phi)$ coincides with the equilibrium if and only if the relative entropy cancels. In the following, we are interested in the evolution of the
discrete counterparts of this quantities, defined as

$$ \mathbb{E}_{\mathcal{D}}^{n} = \left[ \left\| \underline{N}_{\mathcal{D}}^{e} H \left( \frac{\underline{N}_{\mathcal{D}}^{n}}{\underline{N}_{\mathcal{D}}^{e}} \right), \underline{1}_{\mathcal{D}} \right\|_{\mathfrak{M}} + \left\| \underline{P}_{\mathcal{D}}^{e} H \left( \frac{\underline{P}_{\mathcal{D}}^{n}}{\underline{P}_{\mathcal{D}}^{e}} \right), \underline{1}_{\mathcal{D}} \right\|_{\mathfrak{M}} + \frac{\lambda^2}{2} a_{\mathcal{D}} \left( \underline{\phi}_{\mathcal{D}}^{n} - \underline{\phi}_{\mathcal{D}}^{e}, \underline{\phi}_{\mathcal{D}}^{n} - \underline{\phi}_{\mathcal{D}}^{e} \right) \right] $$

for the HFV scheme (where $\underline{1}_{\mathcal{D}}$ is the discrete elements whose coordinates are 1, and the product, quotient and functions are applied coordinate-wise) and similar definition for the DDFV scheme. Note that the HFV entropy does not take into account the edge unknowns of the discrete densities. To compute the discrete equilibrium, we use a nonlinear scheme for (3.4) and get $\underline{\phi}_{\mathcal{D}}^{e}$, then we defined the associated densities following the continuous relations $N^{e} = e^{\alpha_{N}+\phi^{e}}$ and $P^{e} = e^{\alpha_{P}-\phi^{e}}$ (defined as the discrete level coordinate by coordinate).

We consider a test case with physical data $N_{0}^{D} = e$, $N_{1}^{D} = 1$ and $\alpha_{0} = 0$. We also use two different values of the Debye length $\lambda$, respectively 1 and 0.01. We perform simulations on a triangular mesh, with a $\Delta t_{ini} = \Delta t_{max} = 0.1$. On Figure 3.4, we show the evolutions of the

<figure><img src="image_40.png" /><figcaption>Figure 3.4 – Long-time behaviour. Evolution of the discrete relative entropies.</figcaption></figure>

discrete relative entropies along time, for the two values of $\lambda$ and both schemes. As expected, the convergence towards the equilibrium is exponentially fast, as in the continuous framework. Moreover, it is remarkable to notice that the decay rates are almost the same for both schemes. Moreover, with the small Debye length (Figure 3.4b), both schemes are able to capture the behaviour with a very fast evolution far from the equilibrium, then slower once close to it.

## 3.4 Conclusion

We have compared numerically two schemes for drift-diffusion systems. Both schemes handle seamlessly general polygonal meshes, and are designed in order to preserve the positivity of the density as well as the entropy structure of the continuous problem. It appears that the two schemes present similar behaviours on practical test cases, both from a qualitative and quantitative viewpoint.
Chapter 4

# High-order polytopal schemes for advection-diffusion equations: linear and nonlinear approaches

## Outline of the current chapter

<table><tr><td>4.1 Motivations and context</td><td>128</td></tr><tr><td>4.2 Discrete setting and schemes</td><td>131</td></tr><tr><td>4.2.1 Mesh</td><td>131</td></tr><tr><td>4.2.2 Polynomials, discrete unknowns and discrete operators</td><td>131</td></tr><tr><td>4.2.3 Exponential fitting scheme (linear)</td><td>134</td></tr><tr><td>4.2.4 Nonlinear scheme</td><td>135</td></tr><tr><td>4.3 Main features of the schemes</td><td>137</td></tr><tr><td>4.3.1 Exponential fitting scheme</td><td>138</td></tr><tr><td>4.3.2 Nonlinear scheme</td><td>141</td></tr><tr><td>4.4 Numerical results</td><td>149</td></tr><tr><td>4.4.1 Implementation</td><td>150</td></tr><tr><td>4.4.2 Positivity</td><td>153</td></tr><tr><td>4.4.3 Convergence, accuracy and efficiency</td><td>155</td></tr><tr><td>4.4.4 Discrete long-time behaviour</td><td>158</td></tr><tr><td>4.5 Conclusion</td><td>160</td></tr></table>

This chapter is an extended version of [212], which is an accepted proceeding of the FVCA X conference.

We are interested in the high-order approximation of anisotropic advection-diffusion equations on general meshes. We introduce two schemes, based on the Hybrid High-Order methodology. The first one is a linear scheme, while the second one is nonlinear. Both schemes admit solutions and possess a discrete entropy structure,
which ensures that the long-time behaviour of the discrete solutions mimics the continuous one. Moreover, the nonlinear scheme preserves the positivity of the solution at the discrete level. On the contrary, we present numerical evidence indicating that the linear scheme does not preserve positivity, independently of the order considered. Finally, we show on numerical experiments that the nonlinear scheme has optimal order of accuracy, and is less demanding, for a given accuracy, in terms of computational resources than a low-order Hybrid Finite Volume nonlinear method.

## 4.1 Motivations and context

We are interested in the discretisation of a linear advection-diffusion equation on general meshes with a high-order scheme. Our goal is to compare nonlinear structure-preserving high-order methods with similar low-order methods and linear high-order methods. More precisely, let $\Omega$ be an open, bounded, connected polytopal subset of $\mathbb{R}^d$, $d \in \{2, 3\}$, with Lipschitz boundary. We consider the following linear advection-diffusion problem with homogeneous Neumann boundary conditions: find $u: \mathbb{R}_+ \times \Omega \to \mathbb{R}$ solution to

$$ egin{cases} 
\partial_t u - \operatorname{div}(\Lambda(\nabla u + u \nabla \phi)) = 0 & \text{in } \mathbb{R}_+ \times \Omega, \\
\Lambda(\nabla u + u \nabla \phi) \cdot n = 0 & \text{on } \mathbb{R}_+ \times \partial \Omega, \\
u(0, \cdot) = u^{in} & \text{in } \Omega,
\end{cases} 
\quad (4.1) $$

where $n$ is the unit normal vector to $\partial\Omega$ pointing outward from $\Omega$. We assume that the data satisfy:

(i) $\Lambda \in W^{1,\infty}(\Omega; \mathbb{R}^{d \times d})$ is a uniformly elliptic diffusion tensor: there exists $\lambda_b > 0$ such that, for a.e. $x$ in $\Omega$, $\Lambda(x)\xi \cdot \xi \ge \lambda_b |\xi|^2$ for all $\xi \in \mathbb{R}^d$;

(ii) $\phi \in C^1(\bar{\Omega})$ is a regular potential;

(iii) $u^{in} \in L^1(\Omega)$ is a non-negative initial datum, such that $\int_{\Omega} u^{in} \log(u^{in}) < \infty$.

The solutions to (4.1) enjoy some specific and well-known properties. First the mass of the solution is preserved along time, i.e. for almost every $t > 0$,

$$ \int_{\Omega} u(t) = \int_{\Omega} u^{in} = M \quad (4.2) $$

where $M > 0$ is the initial mass. Second, the solution is positive:

$$ \text{for } t > 0, \quad u(t, \cdot) > 0 \quad \text{almost everywhere on } \Omega. \quad (4.3) $$

Last, the solution has a specific long-time behaviour: it converges exponentially fast when $t \to \infty$ towards the thermal equilibrium $u^{\infty}$, solution to the stationary problem associated to (4.1), defined as

$$ u^{\infty} = \frac{M}{\int_{\Omega} e^{-\phi}} e^{-\phi}. \quad (4.4) $$

In order to get a reliable numerical approximation of such problems, one has to preserve these structural properties at the discrete level. It is well-known that Two-Point Flux Approximation (TPFA) finite volume methods are structure-preserving (see [69] for the long-time behaviour), but these methods can only be used on meshes satisfying some orthogonality conditions (with
respect to the inner product induced by $\Lambda$), which essentially restricts the use of these methods to problems with isotropic tensors. On the other hand, a number of finite volume methods using auxiliary unknowns has been introduced for anisotropic problems on general meshes within the past twenty years. On can cite Discrete Duality Finite Volume (DDFV) methods, with additional unknowns on a dual mesh [157, 105], Vertex Approximate Gradient (VAG) methods, with auxiliary unknowns at the vertices of the mesh [124] or Hybrid Finite Volume (HFV) and Mimetic Finite Difference (MFD) methods, with auxiliary unknowns attached to the faces of the mesh [40, 123]. Such methods have proven to be relevant solutions to the anisotropy issue, but none of these linear methods preserves the positivity of the solutions (see [109]). A possible alternative was proposed in [56], with the introduction and analysis of a nonlinear positivity-preserving VAG scheme. The design and analysis of this scheme, as well as of its DDFV and HFV counterparts of [52, 70], relies on the entropy structure of (4.1): there exists some physically motivated quantity, called entropy, which decays along time. Reproducing this structure at the discrete level is a key point to get stability and convergence results.
One can also cite the methods introduced and analysed in [245, 113, 32], which rely on the introduction of nonlinearities within the schemes at the discrete level. These schemes exhibit a good robustness with respect to anisotropy and meshes for linear stationary diffusion. However they rely on nonlinearities which do not correspond to a particular equation at the continuous level. Therefore, it seems difficult to show that they enjoy discrete entropy structures (and associated long-time behaviour).

All the schemes discussed above are at most of order two in space in $L^2$-norm. A natural extension of these works is to develop similar methods with high-order accuracy. Such a task is not easy, because high-order polynomials are well-known to produce oscillations, which makes it difficult to ensure that the discrete solutions stay positive on the whole domain. There exist various positivity-preserving schemes in the literature, based on Discontinuous Galerkin (DG) methods [13, 9, 10]. To our knowledge, the first high-order positive DG scheme for diffusive problems is the one introduced in [198] for linear advection-diffusion. The idea of this scheme is to choose some stabilisation parameters (depending on the degree of the polynomial unknowns and on the dimension of the domain) so that the positivity of the unknowns is ensured in average: the mean value of each cell unknown is positive, even though the (polynomial) unknowns can take negative values inside the cells. The main drawback of such a scheme lies in the difficulty to analyse it: since the solution can become negative on some area of the domain, it is not possible to adapt the analysis techniques used in the continuous framework, such as the use of the Boltzmann entropy.

Another positivity-preserving DG scheme was introduced in [36] for a Fisher-KPP equation $\partial_t u - \Delta u = u(1-u)$. This nonlinear scheme is developed so as to preserve the entropy structure of the continuous problem, and relies on a transformation of the equation based on the variable $\lambda = \log(u)$. The scheme provides positive densities, defined as $u = e^\lambda$. Compared to the high-order DG schemes discussed above, the main improvement lies in the fact that the unknowns are positive everywhere. Such a feature allows the authors to fully analyse the scheme, including existence, long-time behaviour, and convergence towards a semi-discretised solution. The analysis is based on the properties of a well-chosen stabilisation function whose expression implies $L^\infty$-norms of the polynomial unknowns along the faces of the mesh. One can also cite the space-time DG discretisation [37] for cross-diffusion systems.

The aim of this chapter is to introduce a high-order scheme for anisotropic problems on general meshes, preserving the three structural properties (4.2)-(4.3)-(4.4) discussed above. Since the HFV method can be seen as the low-order version of the Hybrid High-Order (HHO) scheme introduced in [100] for stationary diffusion, it is rather natural to try to adapt the nonlinear HFV scheme of [70] to the HHO framework. Our goal here is twofold. First, we want to design
a structure-preserving high-order scheme, and assert from a numerical point of view that the scheme is indeed of high order accuracy. On the other hand, we also want to compare this scheme with other approaches, based on two criteria:

(i) in terms of efficiency, with respect to the nonlinear low-order schemes of Chapter 1;

(ii) in terms of structure preservation, with respect to linear high-order schemes.

For the second criterion we mainly focus on the positivity of the discrete solutions. Indeed, with high-order schemes, one can expect to compute a numerical solution which is very close to the continuous one. Therefore, one could expect that the use of linear high-order schemes (with sufficiently high order) is already a solution in itself to preserve the positivity, at least in practice.

In order to perform this comparison, we introduce two HHO schemes. The first one is linear, and relies on the exponential fitting strategy [42] already used in Chapter 1. The key idea of this scheme is the change of unknown $\rho = ue^{\phi}$, which allows one to reformulate (4.1) into an unconditionally coercive problem in $\rho$. As a by-product of this reformulation, we also get a scheme which preserves the thermal equilibrium. The second scheme is a nonlinear one, which can be seen as a high-order generalisation of the nonlinear HFV scheme introduced in Chapter 1. This scheme is designed so as to preserve the positivity of the discrete solutions. For the sake of simplicity, both schemes rely on a mixed-order HHO space: given an integer $k \ge 0$, the methods hinge on face unknowns of degree $k$ and enriched cell unknowns of degree $k+1$. The main interest of such a discrepancy in the degree between face and cell unknowns is the simplification of the stabilisation term. In the meantime, such a choice preserves optimal accuracy (order $k+2$ in $L^2$-norm) and frugality (the linear system only takes into account the face unknowns, of degree $k$). This kind of method is sometimes referred to as HDG+ in the literature, and is associated to the so-called Lehrenfeld-Schöberl stabilisation [190, 191].

Our first theoretical results stated in Propositions 13 and 14 concern the well-posedness of the exponential fitting scheme and its long-time behaviour. Regarding the nonlinear scheme, we prove the existence of (positive) solutions in Theorem 10. Last, these discrete solutions are proved to converge (in time) towards an associated discrete thermal equilibrium in Proposition 16.

Note that one could also compare the schemes introduced in this chapter with the linear scheme for advection-diffusion of [98], based on separate discretisations of the advection and diffusion. This scheme is a high-order generalisation of the HMM scheme [19] which was considered in Chapter 1. As observed in [70], the low-order version of this scheme has two major drawbacks, namely the fact that it does not preserve the thermal equilibrium, and the need for some coercivity assumptions (which constrain the variety of potentials that can be considered for the analysis). We expect the high-order scheme to inherit these features. Moreover, in order to perform a fair comparison between the schemes, one should also consider a mixed-order version of the scheme of [98], but such a variation has, to our knowledge, not been studied yet. Thus, we will not consider this scheme here, but a numerical comparison including this scheme could be the subject of a future work.

The rest of this chapter is organised as follows. In Section 4.2, we first present the discrete framework, and describe the two schemes under consideration. Then, in Section 4.3 we discuss the main properties of both schemes, and provide some elements of analysis regarding the well-posedness and discrete long-time behaviours. Last, in Section 4.4 we discuss the implementation of the nonlinear scheme and we assert the behaviour of the schemes on some test-cases. These results show that the nonlinear high-order scheme offers a better efficiency in terms of computational cost than low-order schemes, while preserving the positivity of the solution, contrary to linear schemes.
## 4.2 Discrete setting and schemes

In this section, we introduce the discrete framework. Since our focus is mainly on the numerical aspects of the schemes, the description is concise. We refer the reader to [97, Section 1.1] for a more detailed description.

Here, as in the rest of this manuscript, we focus on space discretisation. Hence, both schemes are based on a backward Euler discretisation in time. For the sake of readability, we will consider a constant positive time step $\Delta t$. The time discretisation is then defined as $(t^n)_{n \in \mathbb{N}}$, where $t^n = n\Delta t$. Note that it is rather straightforward to generalise the discussion below to variable time steps.

### 4.2.1 Mesh

We define a discretisation of $\Omega$ as a couple $\mathcal{D} = (\mathcal{M}, \mathcal{E})$, where:

- the mesh $\mathcal{M}$ is a partition of $\Omega$, i.e. $\mathcal{M}$ is a finite collection of disjoint, open, Lipschitz polytopes $K \subset \Omega$ with positive measure $|K| > 0$ (the cells) such that $\bar{\Omega} = \bigcup_{K \in \mathcal{M}} \bar{K}$;

- the set of faces $\mathcal{E}$ is a partition of the mesh skeleton $\bigcup_{K \in \mathcal{M}} \partial K$, i.e. $\mathcal{E}$ is a finite collection of disjoint, connected, relatively open subsets $\sigma$ (the faces) of $\bar{\Omega}$ with positive measure $|\sigma| > 0$ such that $\bigcup_{K \in \mathcal{M}} \partial K = \bigcup_{\sigma \in \mathcal{E}} \bar{\sigma}$. It is assumed that, for all $\sigma \in \mathcal{E}$, $\sigma$ is a Lipschitz polytopal subset of an hyperplane. We denote by $\mathcal{E}_K$ the set of faces of the cell $K$, and we let $n_{K,\sigma} \in \mathbb{R}^d$ be the (constant) unit normal vector to $\sigma \in \mathcal{E}_K$ pointing outward from $K$.

The diameter of a subset $X \subset \bar{\Omega}$ is denoted by

$$h_X = \sup \{ |x - y| \mid (x, y) \in X^2 \}.$$

We define the size of $\mathcal{D}$ as $h_D = \sup\{h_K \mid K \in \mathcal{M}\}$.

When studying asymptotic behaviours with respect to the meshsize, one has to adopt a measure of regularity for (refined) families of discretisations. We classically follow [97, Definition 1.9], in which regularity for a refined family of discretisations is quantified by a uniform (with respect to the meshsize) parameter $0 < \theta < 1$, called regularity parameter. This parameter takes into account the ratios between the different characteristic scales of the mesh (typically, the size of the cells and the size of the faces).

### 4.2.2 Polynomials, discrete unknowns and discrete operators

In the following, $k$ is a fixed non-negative integer. First, we introduce polynomial spaces on a subset $X \subset \bar{\Omega}: \mathbb{P}^k(X)$ and $\mathbb{P}^k(X)^d$ denote respectively the spaces of polynomial functions $X \to \mathbb{R}$ and polynomial vector fields $X \to \mathbb{R}^d$ of total degree at most $k$. If $X$ is a geometric object of Hausdorff dimension $l \in \{1, 2, 3\}$, the dimension of the space $\mathbb{P}^k(X)$ is $\binom{l+k}{k}$. We also define the $L^2(X)$-orthogonal projector $\Pi_X^k: L^1(X) \to \mathbb{P}^k(X)$ such that, given any $v \in L^1(X)$, $\Pi_X^k(v)$ is the only element of $\mathbb{P}^k(X)$ which satisfies

$$\forall w \in \mathbb{P}^k(X), \quad \int_X \Pi_X^k(v)w = \int_X vw.$$

Given any $K \in \mathcal{M}$, $\sigma \in \mathcal{E}_K$ and $v \in W^{1,1}(K)$, we also introduce the shortcut notation

$$\Pi_{\sigma}^{k}(v) = \Pi_{\sigma}^{k}(v|_{\sigma}).$$
We now introduce the set of discrete unknowns corresponding to the mixed-order HHO method [83, 97], with face unknowns of degree $k$ and (enriched) cell unknowns of degree $k+1$:

$$ \underline{V}_D^k = \left\{ \underline{v}_D = \left( (v_K)_{K \in \mathcal{M}}, (v_\sigma)_{\sigma \in \mathcal{E}} \right) \left| \begin{array}{ll} \forall K \in \mathcal{M}, & v_K \in \mathbb{P}^{k+1}(K) \\ \forall \sigma \in \mathcal{E}, & v_\sigma \in \mathbb{P}^k(\sigma) \end{array} \right. \right\}. $$

Given a cell $K \in \mathcal{M}$, we let

$$ \underline{V}_K^k = \mathbb{P}^{k+1}(K) \times \prod_{\sigma \in \mathcal{E}_K} \mathbb{P}^k(\sigma) $$

be the restriction of $\underline{V}_D^k$ to $K$, and for any generic discrete unknown $\underline{v}_D \in \underline{V}_D^k$ we denote by $\underline{v}_K = (v_K, (v_\sigma)_{\sigma \in \mathcal{E}_K}) \in \underline{V}_K^k$ its local restriction to the cell $K$. Given any $\underline{v}_D \in \underline{V}_D^k$, we associate two piecewise polynomial functions $v_M: \Omega \to \mathbb{R}$ and $v_\mathcal{E}: \bigcup_{K \in \mathcal{M}} \partial K \to \mathbb{R}$ such that

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

We also introduce $\underline{1}_D \in \underline{V}_D^k$ the discrete element such that $\underline{1}_K = 1$ for any cell $K \in \mathcal{M}$ and $\underline{1}_\sigma = 1$ for any face $\sigma \in \mathcal{E}$. Last, given a cell $K \in \mathcal{M}$, we define the local interpolator $\underline{I}_K^k: W^{1,1}(K) \to \underline{V}_K^k$ such that, for any $v \in W^{1,1}(K)$,

$$ \underline{I}_K^k(v) = (\Pi_K^{k+1}(v), (\Pi_\sigma^k(v))_{\sigma \in \mathcal{E}_K}). $$

Similarly, the global interpolator $\underline{I}_D^k: W^{1,1}(\Omega) \to \underline{V}_D^k$ is defined as

$$ \underline{I}_D^k(v) = ((\Pi_K^{k+1}(v))_{K \in \mathcal{M}}, (\Pi_\sigma^k(v))_{\sigma \in \mathcal{E}}) \quad \forall v \in W^{1,1}(\Omega). $$

In order to analyse the schemes, one needs to introduce an $H^1$-like discrete semi-norm. To do so, given a cell $K \in \mathcal{M}$, we first introduce a local semi-norm $|\cdot|_{1,K}$ on $\underline{V}_K^k$, defined as

$$ |\underline{v}_K|^2_{1,K} = \|\nabla \underline{v}_K\|_{L^2(K)}^2 + \sum_{\sigma \in \mathcal{E}_K} \frac{1}{h_\sigma} \|\underline{v}_\sigma - \underline{v}_K\|_{L^2(\sigma)}^2. $$

At the global level, we define the discrete $H^1$ semi-norm $|\cdot|_{1,D}$ on $\underline{V}_D^k$ as

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

Note that $|\cdot|_{1,D}$ is not a norm on $\underline{V}_D^k$, but any $\underline{v}_D \in \underline{V}_D^k$ satisfying $|\underline{v}_D|_{1,D} = 0$ is a constant discrete element: there exists $c \in \mathbb{R}$ such that $\underline{v}_D = c \underline{1}_D$. In particular, this implies that $|\cdot|_{1,D}$ is a norm on the space $\underline{V}_{D,0}^k$ of null-mass discrete unknowns, defined as

$$ \underline{V}_{D,0}^k = \left\{ \underline{v}_D \in \underline{V}_D^k \mid \int_{\Omega} v_M = 0 \right\}. $$

The HHO schemes are based on a local approach, which ensures the local conservativity of the method. Given a cell $K \in \mathcal{M}$, we introduce local operators defined on $\underline{V}_K^k$, which are split into a consistent part and a stabilisation part. We first define a local discrete gradient operator
$G_K^K : \underline{V}_K^K \to \mathbb{P}^k(K)^d$ such that, for any $\underline{v}_K \in \underline{V}_K^K$, $G_K^K(\underline{v}_K)$ satisfies

$$ \int_K G_K^K(\underline{v}_K) \cdot \tau = \int_K \nabla v_K \cdot \tau + \sum_{\sigma \in \mathcal{E}_K} \int_{\sigma} (v_{\sigma} - v_K) \tau \cdot n_{K,\sigma} \quad \forall \tau \in \mathbb{P}^k(K)^d. \quad (4.6) $$

This operator is a discrete counterpart of the continuous gradient and is consistent. It satisfies the following commutation property:

$$ \forall v \in W^{1,1}(K), G_K^K \circ \underline{I}_K^K(v) = \Pi_K^K(\nabla v), $$

where $\Pi_K^K$ acts component-wise. However, this operator suffers from a lack of coercivity. To ensure stability, we need to add a stabilisation, which controls the jumps between the cell unknown (or, to be more precise, the trace of this unknown on $\partial K$) and the face unknowns. To do so, given a face $\sigma \in \mathcal{E}_K$, we define the jump operator $J_{K,\sigma} : \underline{V}_K^K \to \mathbb{P}^k(\sigma)$ by

$$ J_{K,\sigma}(\underline{v}_K) = \Pi_{\sigma}^{k}(v_K) - v_{\sigma}. \quad (4.7) $$

One can now introduce the classical (see [83, Section 3.2.1] and [97, Sections 4.2 and 5.1]) discrete counterpart of the local bilinear form $(w, v) \mapsto \int_K \Lambda \nabla w \cdot \nabla v$ as the bilinear form $a_K : \underline{V}_K^K \times \underline{V}_K^K \to \mathbb{R}$ such that

$$ a_K : (\underline{w}_K, \underline{v}_K) \mapsto \int_K \Lambda G_K^K(\underline{w}_K) \cdot G_K^K(\underline{v}_K) + \eta_l \sum_{\sigma \in \mathcal{E}_K} \frac{\Lambda_{K,\sigma}}{h_{\sigma}} \int_{\sigma} J_{K,\sigma}(\underline{w}_K) J_{K,\sigma}(\underline{v}_K), \quad (4.8) $$

where $\Lambda_{K,\sigma} = \|\Lambda|_K n_{K,\sigma} \cdot n_{K,\sigma}\|_{L^\infty(\sigma)}$ and $\eta_l > 0$ is a stabilisation parameter ($l$ stands for linear). The stabilisation used here, based on the jump operators $J_{K,\sigma}$, is the main specificity of the mixed-order HHO method with respect to the equal-order one [100]. Note that such a strategy used in mixed-form (formulation on the flux) is sometimes called HDG+ methods, and the corresponding stabilisation is called Lehrenfeld-Schöberl stabilisation, as it was first introduced in [190, 191]. Here, as in the rest of this manuscript, the HHO scheme is written on primal form.

One can now define a global bilinear form $a_D : \underline{V}_D^K \times \underline{V}_D^K \to \mathbb{R}$, discretisation of the continuous one $(w, v) \mapsto \int_{\Omega} \Lambda \nabla w \cdot \nabla v$, by summing the local contributions:

$$ a_D(\underline{w}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} a_K(\underline{w}_K, \underline{v}_K). \quad (4.9) $$

In terms of analysis, the definition of $a_D$ implies (see for example [97, Section 5.1.6]) that the following stability estimate holds: there exists a positive constant $\alpha_b$, depending only on $\Omega$, $k$ and $\theta$ (the regularity of the mesh) such that

$$ \forall \underline{v}_D \in \underline{V}_D^K, a_D(\underline{v}_D, \underline{v}_D) \geq \lambda_b \alpha_b |\underline{v}_D|^2_{1,D}. \quad (4.10) $$

In particular, since $|\cdot|_{1,D}$ is a norm on $\underline{V}_{D,0}^k$, this estimate implies that $a_D$ is a coercive bilinear form on $\underline{V}_{D,0}^k$.

In order to perform the analysis of the linear scheme, as in the continuous setting, we need a control of the $L^2$-norm of $v_M$ by its energy-norm. This is the purpose of the discrete Poincaré-Wirtinger inequality: there exists a positive constant $C_P$, depending only on $\Omega$, $k$ and $\theta$ such that

$$ \forall \underline{v}_D \in \underline{V}_{D,0}^k, \|v_M\|_{L^2(\Omega)} \leq C_P |\underline{v}_D|_{1,D}. \quad (4.11) $$
We refer to [97, Theorem 6.5, $p = q = 2$] for a demonstration of this result (in the equal-order case).

### 4.2.3 Exponential fitting scheme (linear)

The first scheme introduced herein is a linear one, which follows the exponential fitting strategy. Such an idea has, to our knowledge, initially been introduced at the numerical level in [43] in the context of mixed finite elements. In Chapter 1, we introduce and analyse an HFV scheme based on this strategy. In a nutshell, the exponential fitting approach hinges on the following reformulation of the continuous flux: letting $\omega = e^{-\phi}$, we introduce the Slotboom variable $\rho = \frac{u}{\omega}$. Then, at least formally, one has the following relation:

$$ \nabla u + u \nabla \phi = \omega \nabla \rho, \qquad (4.12) $$

which allows one to transform an advection-diffusion equation into a purely diffusive problem. The main goal of this procedure is to ensure the coercivity of the new problem (since, by regularity of $\phi$, $\omega$ is uniformly bounded away from zero). The exponential fitting scheme is based on this formulation. We discretise the Slotboom variable as a polynomial, that is to say we approximate $\rho$ as a discrete unknown in $\underline{V}_D^k$. Then, we reconstruct the discrete density by mimicking the continuous relation $u = \omega \times \rho$. Thus, a solution $(\rho_D^n)_{n \ge 1}$ to the scheme (4.14) described below corresponds to an approximation of the Slotboom variable associated to the density $u$.

More precisely, for a given discrete Slotboom variable $\rho_D \in \underline{V}_D^k$, we associate a discrete density

$$ \underline{u}_D = ((u_K^\omega)_{K \in \mathcal{M}}, (u_\sigma^\omega)_{\sigma \in \mathcal{E}}) $$

defined as a collection of (non-polynomial) regular functions:

$$ \forall K \in \mathcal{M}, u_K^\omega = \omega_{|K} \rho_K \quad \text{and} \quad \forall \sigma \in \mathcal{E}, u_\sigma^\omega = \omega_{|\sigma} \rho_\sigma. $$

To this discrete density, we associate two piecewise continuous functions: $u_M^\omega : \Omega \to \mathbb{R}$, the discrete density on the cells and $u_\varepsilon^\omega : \bigcup_{K \in \mathcal{M}} \partial K \to \mathbb{R}$, the discrete density on the faces. These reconstructed densities are defined as

$$ u_{M|K}^{\omega} = u_{K}^{\omega} \quad \text{for all } K \in \mathcal{M} \quad \text{and} \quad u_{\varepsilon|\sigma}^{\omega} = u_{\sigma}^{\omega} \quad \text{for all } \sigma \in \mathcal{E}. \qquad (4.13) $$

Note that $\underline{u}_D$ is not a collection of polynomials in general (as highlighted by the use of a double underline).

**Remark 19 (Non-polynomial reconstruction).** Here, we choose to reconstruct the density as a non-polynomial function, since we locally multiply the polynomial unknown $\rho_K$ (or $\rho_\sigma$ for the faces) by the continuous (and, a priori non-polynomial) $\omega$. One could also think of reconstructing a polynomial unknown, by multiplying (component by component) $\underline{\rho}_D$ by $\underline{I}_D^k(\omega)$. Such a polynomial reconstruction is the choice made for the HFV exponential fitting scheme of Chapter 1, in the sense that we use an interpolation of $\omega$ to define the discrete density. But for the low-order scheme of Chapter 1, both polynomials are constants, and therefore the reconstructed density is also constant. For the arbitrary-order scheme presented here, we believe that a polynomial reconstruction of the density is not so relevant. Indeed, the main issue of using $\underline{I}_D^k(\omega)$ lies in the lack of positivity: there is no way of ensuring that the polynomials constituting $\underline{I}_D^k(\omega)$ are positive on their definition domain. Therefore, the use of a
polynomial reconstruction should lead to larger positivity issues. Such issues do not exist in Chapter 1, since the interpolation is positivity preserving.

In order to define a mixed-order HHO exponential fitting scheme, one needs to introduce a discrete counterpart of $(\rho, v) \mapsto \int_{\Omega} \omega \Lambda \nabla \rho \cdot \nabla v$. To do so, given $K \in \mathcal{M}$, leveraging the definition (4.8) of $a_K$, we define at the local level $a_K^\omega: \underline{V}_K^k \times \underline{V}_K^k \to \mathbb{R}$ as the bilinear form

$$a_K^\omega: (\underline{\rho}_K, \underline{v}_K) \mapsto \int_K \omega \Lambda G_K^k(\underline{\rho}_K) \cdot G_K^k(\underline{v}_K) + \eta_l \sum_{\sigma \in \mathcal{E}_K} \frac{\Lambda_{K\sigma}^\omega}{h_\sigma} \int_\sigma J_{K,\sigma}(\underline{\rho}_K) J_{K,\sigma}(\underline{v}_K),$$

where $\Lambda_{K\sigma}^\omega = \|\omega\Lambda_{|K}n_{K,\sigma} \cdot n_{K,\sigma}\|_{L^\infty(\sigma)}$. At the global level, as previously, we construct the bilinear form $a_D^\omega: \underline{V}_D^k \times \underline{V}_D^k \to \mathbb{R}$ by summing the local contributions:

$$a_D^\omega: (\underline{\rho}_D, \underline{v}_D) \mapsto \sum_{K \in \mathcal{M}} a_K^\omega(\underline{\rho}_K, \underline{v}_K).$$

We can now introduce the exponential fitting scheme for Problem (4.1): find $(\underline{\rho}_D^n)_{n \ge 1} \in (\underline{V}_D^k)^{\mathbb{N}^*}$ such that, for all $n \in \mathbb{N}$,

$$\left\{ \int_{\Omega} \frac{u_{M}^{\omega, n+1} - u_{M}^{\omega, n}}{\Delta t} v_{M} = -a_{D}^{\omega}(\rho_{D}^{n+1}, \underline{v}_{D}) \quad \forall \underline{v}_{D} \in \underline{V}_{D}^{k}, \quad (4.14a) \right. $$

$$u_{M|K}^{\omega, n+1} = \omega_{|K} \rho_K^{n+1} \quad \forall K \in \mathcal{M}, \quad (4.14b)$$

$$u_{M|K}^{\omega, 0} = u_{|K}^{in} \quad \forall K \in \mathcal{M}. \quad (4.14c)$$

Given a solution $(\rho_D^n)_{n \ge 1}$ to (4.14), as described above, we associate a sequence of reconstructed densities $(\underline{u}_D^n)_{n \ge 1}$.

One can check that, for any $c \in \mathbb{R}$, the discrete element $\underline{\rho}_D = c1_D \in \underline{V}_D^k$ is a stationary solution to the exponential fitting scheme (4.14) (in Slotboom variable). Therefore, we define the discrete Slotboom variable at equilibrium $\underline{\rho}_D^\infty \in \underline{V}_D^k$ as

$$\underline{\rho}_D^\infty = c_\omega^M 1_D, \text{ where } c_\omega^M = \frac{M}{\int_{\Omega} \omega}.$$

It is naturally associated to the the discrete thermal equilibrium $\underline{u}_D^\infty = ((u_K^{\omega, \infty})_{K \in \mathcal{M}}, (u_\sigma^{\omega, \infty})_{\sigma \in \mathcal{E}})$. It is the only stationary solution to (4.14) associated to a density of mass M.

### 4.2.4 Nonlinear scheme

Following the ideas from [56, 70], our nonlinear scheme relies on a nonlinear reformulation of Problem (4.1). To do so, we introduce the logarithm potential $\ell = \log(u)$ and the quasi-Fermi potential $w = \ell + \phi$. At least formally, if $u$ is positive, one has the following relation:

$$\nabla u + u \nabla \phi = u \nabla (\log(u) + \phi) = e^\ell \nabla w. \quad (4.15)$$

We choose to discretise the potentials as polynomials, i.e. to approximate $\ell$ and $w$ as discrete unknowns in $\underline{V}_D^k$. Then, mimicking the relation $u = e^\ell$, we reconstruct the density as the
exponential of polynomial functions, thus ensuring its positivity. Therefore, a solution $(\underline{\ell}_D^n)_{n \ge 1}$ to the nonlinear scheme (4.20) described below corresponds to an approximation of the logarithms of the solution $u$ (density).

More specifically, for a given discretisation $\underline{\ell}_D \in \underline{V}_D^k$ of the potential $\ell$, one associates a discrete density

$$u_D = ((u_K)_{K \in \mathcal{M}}, (u_\sigma)_{\sigma \in \mathcal{E}}),$$

defined as a collection of positive smooth functions

$$\forall K \in \mathcal{M}, u_K = \exp(\ell_K) \quad \text{and} \quad \forall \sigma \in \mathcal{E}, u_\sigma = \exp(\ell_\sigma).$$

As previously, we associate to this discrete density a couple of piecewise smooth functions $u_M : \Omega \to \mathbb{R}$ (corresponding to the cell unknowns) and $u_\varepsilon : \bigcup_{K \in \mathcal{M}} \partial K \to \mathbb{R}$ (corresponding to the face unknowns), defined as

$$u_M = \exp(\ell_M) \text{ and } u_\varepsilon = \exp(\ell_\varepsilon). \qquad (4.16)$$

Note that a discrete density $u_D$ is not a collection of polynomials (which is highlighted by the use of the wave under $u$), but it enjoys positivity, both on the cells and faces, since it is defined as the exponential of real functions.

As previously, this scheme is based on local cell contributions, split into a consistent term and a stabilisation term. Given a cell $K \in \mathcal{M}$, we define a discretisation of $(\ell, w, v) \mapsto \int_K e^\ell \Lambda \nabla w \cdot \nabla v$ as a sum of nonlinear consistent (4.17a) and stabilising (4.17b) contributions:

$$C_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) = \int_K e^{\ell_K} \Lambda G_K^k(\underline{w}_K) \cdot G_K^k(\underline{v}_K), \qquad (4.17a)$$

$$S_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) = \eta_{nl} \sum_{\sigma \in \mathcal{E}_K} \frac{\Lambda_{K\sigma}}{h_\sigma} \int_\sigma \frac{e^{\Pi_\sigma^k(\ell_K)} + e^{\ell_\sigma}}{2} J_{K,\sigma}(\underline{w}_K) J_{K,\sigma}(\underline{v}_K), \qquad (4.17b)$$

where $\eta_{nl} > 0$ is a stabilisation parameter (nl stands for nonlinear). We can now define a local application $\mathcal{T}_K : \underline{V}_K^k \times \underline{V}_K^k \times \underline{V}_K^k \to \mathbb{R}$ by

$$\mathcal{T}_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) = C_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) + S_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K) + \varepsilon h_K^{k+2} a_K(\underline{w}_K, \underline{v}_K), \qquad (4.18)$$

where $\varepsilon$ is a non-negative parameter and $a_K$ is the bilinear form defined in (4.8)-(4.9). At the global level, we define $\mathcal{T}_D : \underline{V}_D^k \times \underline{V}_D^k \times \underline{V}_D^k \to \mathbb{R}$ by summing the local contributions:

$$\mathcal{T}_D(\underline{\ell}_D, \underline{w}_D, \underline{v}_D) = \sum_{K \in \mathcal{M}} \mathcal{T}_K(\underline{\ell}_K, \underline{w}_K, \underline{v}_K). \qquad (4.19)$$

**Remark 20** (Parameter $\varepsilon$). The operator $\mathcal{T}_D$ is to be understood as a discretisation of $(\ell, w, v) \mapsto \int_\Omega (e^\ell + \varepsilon) \Lambda \nabla w \cdot \nabla v$, with $\varepsilon \sim \varepsilon h_D^{k+2}$ a small parameter. The $\varepsilon$ perturbation is used in order to show the existence result of Proposition 10 and can be seen as a kind of stabilisation. The scaling factor $h_K^{k+2}$ in (4.18) is used to get the expected order of convergence. Although the choice made here gives the right order of convergence, the use of larger parameters (such as $h_K^{k+1}$) has not yet been studied. In practice, numerical results for $\varepsilon = 1$ and $\varepsilon = 0$ are almost the same. The influence of this term will be further investigated in future works.

We now let $\underline{\phi}_D = \underline{I}_D^k(\underline{\phi}) \in \underline{V}_D^k$ be the interpolate of $\underline{\phi}$. Now, using a backward Euler dis-
cretisation in time with time step $\Delta t > 0$, we introduce the following scheme for (4.1): find $(\underline{\ell}_D^n)_{n \ge 1} \in (\underline{V}_D^k)^{N^*}$ such that, for all $n \ge 0$,

$$ \left\{ \int_{\Omega} \frac{u_M^{n+1} - u_M^n}{\Delta t} v_M = -T_D(\underline{\ell}_D^{n+1}, \underline{\ell}_D^{n+1} + \underline{\phi}_D, \underline{v}_D) \quad \forall \underline{v}_D \in \underline{V}_D^k, \quad (4.20a) $$

$$ u_K^{n+1} = e^{\underline{\ell}_K^{n+1}} \quad \forall K \in \mathcal{M}, \quad (4.20b) $$

$$ u_K^0 = u^{in}_{|K} \quad \forall K \in \mathcal{M}. \quad (4.20c) $$

Given a solution $(\underline{\ell}_D^n)_{n \ge 1}$ to the scheme (4.20) (which are approximations of the logarithm of the density at times $t^n = n\Delta t$), as discussed above, we associate a sequence of positive discrete densities $(u_D^n)_{n \ge 1}$.

**Remark 21** (Initial condition). Note that the scheme (4.20) does not rely on discrete initial data, as described in (4.20c). This strategy allows one not to define $\underline{\ell}_D^0$ (as the interpolate of $\underline{\ell}^{in} = \log(u^{in})$), which is not possible in regions where $u^{in}$ vanishes.

The question of defining an initial discrete datum is a major difficulty when it comes to the numerical implementation and the initialisation of the Newton method (see Section 4.4.1).

We define the discrete thermal equilibrium as $u_D^\infty = ((u_K^\infty)_{K \in \mathcal{M}}, (u_\sigma^\infty)_{\sigma \in \mathcal{E}})$ where

$$ \forall K \in \mathcal{M}, \quad u_K^\infty = c_{nl}^M e^{-\phi_K} \quad \text{and} \quad \forall \sigma \in \mathcal{E}, \quad u_\sigma^\infty = c_{nl}^M e^{-\phi_\sigma}, $$

with $c_{nl}^M = M / \int_{\Omega} e^{-\phi_M}$. One can show that $u_D^\infty$ (and the associated logarithm potential at equilibrium $\underline{\ell}_D^\infty = \log(c_{nl}^M) \underline{1}_D - \underline{\phi}_D \in \underline{V}_D^k$) is the only stationary solution to (4.20) with mass $M$.

## 4.3 Main features of the schemes

In this section, we present some theoretical results about the two schemes introduced above. We especially focus on the existence (and stability) of solutions, as well as on questions related to the long-time behaviour. The results presented below are generalisations of the analysis for low-order schemes presented in Chapter 1. In particular, the analysis strongly relies on the entropy structure of both schemes.

**Remark 22** (Low-order versions of the schemes ($k=0$)). Note that the low-order versions of the two schemes introduced above do not coincide with the exponential fitting and nonlinear HFV schemes of Chapter 1. Indeed, even if the equal-order HHO method (face and cell unknowns of the same degree $k$) for $k=0$ is essentially equivalent to the HFV method for diffusion problems [100], such a result does not hold for the mixed-order HHO method. In fact, the low-order versions of the schemes (4.14) and (4.20) have enriched cell unknowns, which are affine functions (and not constant). Therefore, the results presented here are new, even for $k=0$.

Before presenting individual results for each scheme, it should be noted that both schemes exhibits a similar important property: the preservation of the thermal equilibrium. Indeed, both scheme admit an unique stationary solution with mass $M$, namely $\underline{\rho}_D^\infty = c_\omega^M \underline{1}_D$ and $\underline{\ell}_D^\infty = \log(c_{nl}^M) \underline{1}_D - \underline{\phi}_D$, which are associated to the following discrete cell densities:

$$ u_M^{\omega, \infty} = \frac{M}{\int_{\Omega} e^{-\phi}} e^{-\phi} \quad \text{and} \quad u_M^\infty = \frac{M}{\int_{\Omega} e^{-\phi_M}} e^{-\phi_M}. $$
One can notice that $u_M^{\omega,\infty}$ coincides with the continuous thermal equilibrium $u^\infty$. This is a consequence of the choice of a non-polynomial reconstruction of the density in the exponential fitting scheme (see Remark 19). For the nonlinear scheme, we do not preserve the exact thermal equilibrium, since the discrete densities for this scheme are defined as exponentials of polynomial functions (so, contrarily to the scheme of Chapter 1, the discrete equilibrium is not $I_D^k(u^\infty)$). However, the discrete thermal equilibrium $u_M^\infty$ is a reasonable discretisation of $u^\infty$, in the sense that the logarithm potential of the discrete thermal equilibrium and the interpolate of the continuous equilibrium logarithm potential are equal up to an additive constant:

$$ \forall K \in \mathcal{M}, \log(u_K^\infty) - \Pi_K^{k+1}(\log(u^\infty)) = \log\left(\int_{\Omega} e^{-\phi}\right) - \log\left(\int_{\Omega} e^{-\phi_M}\right). $$

In fact, at the continuous level, the thermodynamic consistency is sometimes expressed as follows: if the advection-diffusion flux $J = \nabla u + u\nabla\phi$ vanishes on $\Omega$, then the quasi-Fermi potential $\log(u) + \phi$ is constant in $\Omega$. Both schemes satisfy discrete counterparts of this principle: for the exponential fitting scheme, $\log(u_M^\omega) + \phi$ and $\log(u_E^\omega) + \phi|_{\partial\mathcal{M}}$ are constant; while for the nonlinear scheme, $\underline{\ell}_D + \underline{\phi}_D$ is a constant discrete element (i.e. proportional to $\underline{1}_D$).

### 4.3.1 Exponential fitting scheme

We present here the main properties of the exponential fitting scheme (4.14), and give detailed proofs of the results.

As a preliminary remark, note that since $\phi$ is $C^1$ on $\bar{\Omega}$, there exist positive constants $\omega_b$ and $\omega_\#^n$ such that

$$ \forall x \in \Omega, \omega_b \le \omega(x) \le \omega_\#^n. $$

As a consequence, the diffusion tensor $\omega\Lambda$ is uniformly elliptic. This implies a stability estimate on the bilinear form $a_D^\omega$ on $\underline{V}_D^k$: there exists a positive constant $\alpha_b^\omega = \omega_b\lambda_b a_b$ (which therefore depends only on $\Omega$, $k$, $\phi$, $\Lambda$ and the regularity of the mesh $\Theta$) such that

$$ \forall \underline{v}_D \in \underline{V}_D^k, a_D^\omega(\underline{v}_D, \underline{v}_D) \ge \alpha_b^\omega |\underline{v}_D|_{1,D}^2. \qquad (4.21) $$

We first state an existence result, which is mainly a consequence of the previous stability estimate.

**Proposition 13** (Well-posedness of the exponential fitting scheme). *The exponential fitting scheme (4.14) admits an unique solution $(\underline{\rho}_D^n)_{n \ge 1}$. Moreover, the corresponding discrete densities $(\underline{u}_D^n)_{n \ge 1}$ have a mass equal to M:*

$$ \forall n \ge 1, \int_{\Omega} u_M^{\omega,n} = \int_{\Omega} u^{in} = M. \qquad (4.22) $$

*Proof.* Let $n \ge 0$, and assume that $u_M^{\omega,n}$ is defined. We want to show that equation (4.14b) admits a unique solution. To do so, we first define, for any $\underline{\rho}_D \in \underline{V}_D^k$,

$$ \|\underline{\rho}_D\|_{\omega, \Delta t, D} = \sqrt{\Delta t |\underline{\rho}_D|_{1, D}^2 + \int_{\Omega} \omega \underline{\rho}_M^2}. $$

Since $|\cdot|_{1,D}$ is a semi-norm on $\underline{V}_D^k$, it follows that $\|\cdot\|_{\omega,\Delta t,D}$ is an Euclidean norm on $\underline{V}_D^k$. Thus,
by (4.21), the bilinear form $A : (\underline{\rho}_D, \underline{v}_D) \mapsto \int_{\Omega} \omega \rho_M v_M + \Delta t a_D^\omega(\underline{\rho}_D, \underline{v}_D)$ satisfies the following coercivity estimate:

$$ \forall \underline{\rho}_D \in \underline{V}_D^k, A(\underline{\rho}_D, \underline{\rho}_D) \geq \min(1, \alpha_b^\omega) \|\underline{\rho}_D\|_{\omega, \Delta t, D}^2. $$

Therefore, by Lax-Milgram theorem, equation (4.14b) admits an unique solution $\underline{\rho}_D^{n+1}$ in $\underline{V}_D^k$. To get the mass preservation identity (4.22), we test (4.14a) by $\underline{1}_D$ to get

$$ \int_{\Omega} \frac{u_M^{\omega, n+1} - u_M^{\omega, n}}{\Delta t} = 0. $$

We conclude by noticing that $\int_{\Omega} u_M^{\omega, 0} = \int_{\Omega} u^{in}$ according to (4.14c). $\square$

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

We now state our main result about the exponential fitting scheme, which ensures that the solution to (4.14) has similar long-time behaviour as the continuous solution. As usual with the entropy method, the main idea is to get a control of the entropy by its dissipation. Here, such an estimate is a consequence of the discrete Poincaré inequality (4.11).

**Proposition 14** (Long-time behaviour of the exponential fitting scheme). Assume that $u^{in} \in L^2(\Omega)$, and let $(\underline{\rho}_D^n)_{n \ge 1}$ be the solution to the exponential fitting scheme (4.14). Then, the following discrete entropy relation holds:

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

where the discrete quadratic entropy is defined as

$$ \mathbb{E}_{\omega}^{n} = \frac{1}{2} \int_{\Omega} \omega |\rho_{M}^{n} - \rho_{M}^{\infty}|^{2} \text{ with } \rho_{M}^{0} = \frac{u^{in}}{\omega} \in L^{2}(\Omega) $$

and the discrete dissipation is defined by

$$ \mathbb{D}_{\omega}^{n} = a_{D}^{\omega}(\underline{\rho}_{D}^{n} - \underline{\rho}_{D}^{\infty}, \underline{\rho}_{D}^{n} - \underline{\rho}_{D}^{\infty}) \forall n \ge 1. $$

As a consequence, the discrete density converges exponentially fast in time towards the corresponding discrete density at equilibrium: there exist positive constants $C_\omega$ and $v_\omega$ only depending on the data and on the mesh regularity parameter $\theta$ such that

$$ \forall n \in \mathbb{N}, \|u_M^{\omega, n} - u_M^{\omega, \infty}\|_{L^2(\Omega)} \le C_\omega (1 + v_\omega \Delta t)^{-\frac{n}{2}} \|u^{in} - u_M^{\omega, \infty}\|_{L^2(\Omega)}. \quad (4.24) $$

*Proof.* Let $n \in \mathbb{N}$. By convexity of $x \mapsto x^2$ on $\mathbb{R}$, one has

$$ \mathbb{E}_{\omega}^{n+1} - \mathbb{E}_{\omega}^{n} = \frac{1}{2} \int_{\Omega} \omega (|\rho_{M}^{n+1} - \rho_{M}^{\infty}|^2 - |\rho_{M}^{n} - \rho_{M}^{\infty}|^2) \le \int_{\Omega} \omega (\rho_{M}^{n+1} - \rho_{M}^{n})(\rho_{M}^{n+1} - \rho_{M}^{\infty}). $$

Therefore, testing (4.14a) against $\underline{\rho}_D^{n+1} - \underline{\rho}_D^\infty \in \underline{V}_D^k$, we get

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

Note that this estimate holds true also for $n = 0$ (using the definition of $\rho_M^0$). On the other hand,
by definition of $\rho_D^\infty$ (proportional to $\underline{v}_D$) alongside with the bilinearity of $a_D^\omega$, one has

$$a_D^\omega(\rho_D^\infty, \rho_D^{n+1} - \rho_D^\infty) = 0 \text{ and } a_D^\omega(\rho_D^{n+1}, \rho_D^{n+1} - \rho_D^\infty) = \mathbb{D}_\omega^{n+1},$$

which yields the entropy relation (4.23). To get the exponential decay, one needs to compare $\mathbb{D}_\omega^{n+1}$ with $\mathbb{E}_\omega^{n+1}$. To do so, we let $\underline{v}_D = \rho_D^{n+1} - \rho_D^\infty \in \underline{V}_D^k$, and we define the probability measure $\mu = \frac{\omega}{M_\omega}dx$ on $\Omega$, where $M_\omega = \int_\Omega \omega$. We define $\mu v_M$ as the mass of $v_M$ for the measure $\mu$, i.e.

$$\mu v_M = \int_{\Omega} v_M d\mu.$$

The mass preservation identity (4.22) implies that

$$\mu v_M = \int_{\Omega} v_M d\mu = \frac{1}{M_\omega} \int_{\Omega} \omega v_M = \frac{1}{M_\omega} \int_{\Omega} (u_M^{\omega, n+1} - u_M^{\omega, \infty}) = \frac{M-M}{M_\omega} = 0.$$

Therefore, letting $\bar{v} = \frac{1}{|\Omega|} \int_{\Omega} v_M$, we apply [51, Lemma 5.2, q = 2 and $\mu$ as defined above] to get

$$\int_{\Omega} \frac{\omega}{M_{\omega}} v_{M}^{2} = \int_{\Omega} (v_{M} - \mu v_{M})^{2} d\mu \leq 4 \int_{\Omega} (v_{M} - \bar{v})^{2} d\mu.$$

Using the upper bound on $\omega$ and the definition of $\mu$, we infer that

$$\int_{\Omega} \omega v_M^2 \leq 4\omega_{\#} \int_{\Omega} (v_M - \bar{v})^2 = 4\omega_{\#} \|v_M - \bar{v}\|_{L^2(\Omega)}^2.$$

By definition of $\bar{v}$, one has $\underline{v}_D - \bar{v} \underline{v}_D \in \underline{V}_{D,0}^k$ so we can apply the discrete Poincaré inequality (4.11) to infer that

$$\|v_M - \bar{v}\|_{L^2(\Omega)}^2 \leq C_P^2 |\underline{v}_D - \bar{v} \underline{v}_D|^2_{1,D} = C_P^2 |\underline{v}_D|^2_{1,D}.$$

Combining the two previous estimates, we get

$$\mathbb{E}_{\omega}^{n+1} \leq 2\omega_{\#} C_P^2 |\underline{v}_D|^2_{1,D}.$$

Now, one can use the stability estimate (4.21) on $a_D^\omega$ to deduce that

$$\mathbb{D}_{\omega}^{n+1} = a_{D}^{\omega}(\underline{v}_{D}, \underline{v}_{D}) \geq \alpha_{b}^{\omega} |\underline{v}_{D}|_{1, D}^{2}.$$

Therefore, by the last two estimates, one has the following relation between the entropy and its dissipation:

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

Plugging this estimate into the entropy relation (4.23), we deduce that

$$(1 + v_{\omega} \Delta t) \mathbb{E}_{\omega}^{n+1} \leq \mathbb{E}_{\omega}^{n}, \text{ with } v_{\omega} = \frac{\alpha_{b}^{\omega}}{2\omega_{\#} C_P^2}.$$

This implies the exponential decay of the entropy:

$$\forall n \geq 0, \mathbb{E}_{\omega}^{n} \leq (1 + v_{\omega} \Delta t)^{-n} \mathbb{E}_{\omega}^{0}.$$
To conclude, we just use the bounds on $\omega$ to infer that

$$ \omega_b \|u_M^{\omega,n} - u_M^{\omega,\infty}\|_{L^2(\Omega)}^2 \le 2\mathbb{E}_{\omega}^n \le \omega_{\#} \|u_M^{\omega,n} - u_M^{\omega,\infty}\|_{L^2(\Omega)}^2 $$

which yields (4.24) with $C_{\omega} = \sqrt{\frac{\omega_{\#}}{\omega_b}}$. $\square$

**Remark 23** (Regularity of the initial datum and topology of the convergence). Notice that in Proposition 14 we made the extra assumption that $u^{in} \in L^2(\Omega)$. Indeed, the long-time analysis of the exponential fitting scheme relies on the decay of the quadratic entropy (on the Slotboom variable)

$$ \mathbb{E}_{\omega}(t) = \frac{1}{2} \int_{\Omega} \omega |\rho(t) - \rho^{\infty}|^2. $$

Therefore, even at the continuous level, one needs to have a finite initial quadratic entropy to perform the analysis, and assuming that the initial datum is in $L^2(\Omega)$ is a safe choice. At the end, the exponential fitting scheme gives a convergence (in time) result in the $L^2$ topology, as a reminiscence of its linearity. On the other hand, the nonlinear scheme will yield convergence in a weaker norm (typically $L^1$), but can be used to deal with less regular initial data, which are in $L\log(L)$ only.

### 4.3.2 Nonlinear scheme

We present here some results regarding the analysis of the nonlinear scheme (4.20). Since we deal with a nonlinear scheme, unlike the exponential fitting scheme, the question of the existence of solutions is the main difficulty here. As often for this type of scheme, we start by establishing some a priori estimates, and then give an existence result.

For the purpose of analysis, given $\underline{\ell}_D \in \underline{V}_D^k$ a discrete logarithm potential, we associate a discrete quasi-Fermi potential defined as

$$ \underline{w}_D = \underline{\ell}_D + \underline{\phi}_D - \log(c_{nl}^M) \underline{1}_D. $$

By definition of $c_{nl}^M$, one has $w_M = \log(\frac{u_M}{u_M^{\infty}})$. Note that, on the other hand, for any $(\underline{\ell}_D, \underline{v}_D) \in \underline{V}_D^k \times \underline{V}_D^k$, we have

$$ \mathcal{T}_D(\underline{\ell}_D, \underline{\ell}_D + \underline{\phi}_D, \underline{v}_D) = \mathcal{T}_D(\underline{\ell}_D, \underline{w}_D, \underline{v}_D), $$

since $\underline{\ell}_D + \underline{\phi}_D$ and $\underline{w}_D$ only differ from an additive constant element. As in Chapter 2, the discrete quasi-Fermi potentials are the key variables in order to perform the existence analysis of the scheme.

As a last remark before analysing the scheme, note that by stability of the interpolator (see [97, Proposition 2.2]), there exists $\phi_{\#} > 0$ depending only on $\phi$, $\Omega$, $k$ and $\theta$ such that

$$ |\underline{\phi}_D|_{1,D} \le \phi_{\#}. \qquad (4.25) $$

We now state our fundamental a priori results. As for the exponential fitting scheme, the discrete entropy structure of the nonlinear scheme is mainly a consequence of the convexity of the entropy.

**Proposition 15** (Fundamental a priori relations). Let $(\underline{\ell}_D^n)_{n \ge 1}$ be a solution to the nonlinear scheme (4.20), and $(u_D^n)_{n \ge 1}$ be the associated reconstructed discrete density. Then, the following a priori results hold:
(i) the mass is preserved along time:

$$ \forall n \in \mathbb{N}^*, \quad \int_{\Omega} u_M^n = \int_{\Omega} u^{in} = M, \qquad (4.26) $$

(ii) a discrete entropy/dissipation relation holds true:

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

where the discrete entropy and dissipation are non-negative quantities defined by

$$ \mathbb{E}^n = \int_{\Omega} u_M^{\infty} \Phi_1 \left( \frac{u_M^n}{u_M^{\infty}} \right) \ge 0 \quad \text{and} \quad \mathbb{D}^n = \mathcal{T}_D(\underline{\ell}_D^n, \underline{w}_D^n, \overline{w}_D^n) \ge 0 \text{ for } n \ge 1 $$

with $\Phi_1: s \mapsto s \log(s) - s + 1$ (and $\Phi_1(0) = 1$).

*Proof.* Let $n \ge 0$. Using $\underline{1}_D$ as a test function in (4.20a), we get that the mass is preserved:

$$ \int_{\Omega} u_M^{n+1} = \int_{\Omega} u_M^n. $$

Therefore, by (4.20c), we get the identity (4.26).

To get the entropy relation, we use the convexity of $\Phi_1$ and notice that

$$ \mathbb{E}^{n+1} - \mathbb{E}^n \le \int_{\Omega} u_M^{\infty} \Phi_1' \left( \frac{u_M^{n+1}}{u_M^{\infty}} \right) \frac{u_M^{n+1} - u_M^n}{u_M^{\infty}}. $$

Then, since $w_M^{n+1} = \log(\frac{u_M^{n+1}}{u_M^n})$ and $\Phi_1' = \log$, one has

$$ \mathbb{E}^{n+1} - \mathbb{E}^n \le \int_{\Omega} w_M^{n+1} (u_M^{n+1} - u_M^n). \qquad (4.28) $$

On the other hand, testing (4.20a) with $\underline{w}_D^{n+1} \in \underline{V}_D^k$, we get by definition of the discrete dissipation

$$ \int_{\Omega} w_M^{n+1} (u_M^{n+1} - u_M^n) = -\Delta t \mathcal{T}_D(\underline{\ell}_D^{n+1}, \underline{\ell}_D^{n+1} + \underline{\phi}_D, \overline{w}_D^{n+1}) = -\Delta t \mathcal{T}_D(\underline{\ell}_D^{n+1}, \overline{w}_D^{n+1}, \overline{w}_D^{n+1}), $$

which finally yields (4.27).

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

Note that the previous results hold for any $\varepsilon \ge 0$. Following the ideas of [56, 70], the entropy/dissipation relation should allow one to show the exponential convergence of the discrete solutions towards the discrete thermal equilibrium and to get convergence results (when $(h_D, \Delta t) \to (0, 0)$). These aspects will be the topics of future works.

In the remainder of this section, we focus on the existence of solutions and on the long-time behaviour, for $\varepsilon > 0$. The two results stated rely on a discrete a priori estimate, which is a high-order counterpart of Lemma 2.

In order to perform the analysis, we first introduce an inner product $\langle \cdot, \cdot \rangle$ on $\underline{V}_D^k$ such that

$$ \forall (\underline{w}_D, \underline{v}_D) \in \underline{V}_D^k \times \underline{V}_D^k, \langle \underline{w}_D, \underline{v}_D \rangle = \int_{\Omega} w_M v_M + \sum_{\sigma \in \mathcal{E}} \int_{\sigma} w_{\sigma} v_{\sigma}. $$
We denote by $\|\cdot\|$ the corresponding Euclidean norm:

$$ \forall \underline{v}_D \in \underline{V}_D^k, \|\underline{v}_D\| = \sqrt{\sum_{K \in \mathcal{M}} \|v_K\|_{L^2(K)}^2 + \sum_{\sigma \in \mathcal{E}} \|v_\sigma\|_{L^2(\sigma)}^2}. $$

We now state the fundamental discrete a priori result. Note that the proof proposed here adapts easily to non-Boltzmann statistics, since we do not use the morphism properties of the exponential function (but only the fact that it is non-decreasing).

**Lemma 9** (Discrete boundedness by mass and energy semi-norm). Let $\underline{\ell}_D \in \underline{V}_D^k$ and assume that there exist $C_\# > 0$ and $M_\# \ge M_b > 0$ such that

$$ M_b \le \int_{\Omega} e^{\underline{\ell}_M} \le M_\# \quad \text{and} \quad |\underline{\ell}_D|^2_{1,D} \le C_\#. \qquad (4.29) $$

Then, there exists a positive constant $C$ depending only on $M_b, M_\#, C_\#, k, \Omega$ and $\mathcal{D}$ such that

$$ \|\underline{\ell}_D\| \le C. $$

*Proof.* In order to prove this lemma, we need two preliminary results.

First, we recall that in a finite-dimensional vector space, all the norms are equivalent. Therefore, there exist some positive constants $c_\infty$ and $c_2$ (depending only on $\mathcal{D}$ and $k$) such that

$$ \forall K \in \mathcal{M}, \forall v \in \mathbb{P}^{k+1}(K), \left\{ \begin{array}{ll} \|v\|_{L^\infty(K)} & \le c_\infty \|v\|_{L^2(K)}, \\ \|v\|_{L^2(K)} & \le c_2 \|v\|_{L^\infty(K)}, \\ \|v\|_{L^2(\sigma)} & \le c_2 \|v\|_{L^\infty(\sigma)} \quad \forall \sigma \in \mathcal{E}_K, \\ \|v\|_{L^\infty(\sigma)} & \le c_\infty \|v\|_{L^2(\sigma)} \quad \forall \sigma \in \mathcal{E}_K, \\ \|\nabla v\|_{L^\infty(K)} & \le c_\infty \|\nabla v\|_{L^2(K)}, \end{array} \right. $$

and

$$ \forall \sigma \in \mathcal{E}, \forall v \in \mathbb{P}^k(\sigma), \left\{ \begin{array}{ll} \|v\|_{L^\infty(\sigma)} & \le c_\infty \|v\|_{L^2(\sigma)}, \\ \|v\|_{L^2(\sigma)} & \le c_2 \|v\|_{L^\infty(\sigma)}. \end{array} \right. $$

Note that these inequalities essentially correspond to reverse Lebesgue and discrete trace inequalities, for which the precise dependency of the constant with respect to the considered cell is known [97, Section 1.2].

Second, let $K \in \mathcal{M}$. By the mean value theorem (see [118, Remark 2.12]), one shows that there exists some positive constant $c_K$ depending on $K$ (essentially its size and shape) such that

$$ |\inf_K \ell_K - \sup_K \ell_K| \le c_K \|\nabla \ell_K\|_{L^\infty(K)} \le c_K c_\infty \|\nabla \ell_K\|_{L^2(K)} \le c_K c_\infty \sqrt{C_\#}. $$

Therefore, defining the maximal oscillation over the cells as

$$ o_D = c_\infty \sqrt{C_\#} \max_{K \in \mathcal{M}} c_K, $$

we get that

$$ \forall K \in \mathcal{M}, |\inf_K \ell_K - \sup_K \ell_K| \le o_D. \qquad (4.30) $$

We now proceed in three steps. First, we get an $L^\infty$-control over one of the cell unknowns. Then, we propagate this bound to the neighbouring unknowns. Last, we use the connectivity of
the mesh to get an estimate on $||\ell_D||$.

Since $M_b \le \int_{\Omega} e^{\ell_M} = \sum_{K \in \mathcal{M}} \int_K e^{\ell_K} \le M_\#$, there exists a cell $K_0 \in \mathcal{M}$ such that

$$ \frac{|K_0|M_b}{|\Omega|} \le \int_{K_0} e^{\ell_{K_0}} \le M_\# $$

Let $\ell_b = \inf_{K_0} \ell_{K_0}$ and $\ell_\# = \sup_{K_0} \ell_{K_0}$ be the extremal values of $\ell_{K_0}$ on $K_0$. According to (4.30), one has

$$ \ell_\# - o_D \le \ell_b. $$

Thus, integrating this inequality we deduce that

$$ \int_{K_0} e^{\ell_\# - o_D} \le \int_{K_0} e^{\ell_{K_0}} \le M_\# \quad \text{so} \quad \ell_\# \le \log\left(\frac{M_\#}{|K_0|}\right) + o_D. $$

For the other bound, one has

$$ \frac{|K_0|M_b}{|\Omega|} \le \int_{K_0} e^{\ell_{K_0}} \le \int_{K_0} e^{\ell_b + o_D} \quad \text{so} \quad \ell_b \ge \log\left(\frac{M_b}{|\Omega|}\right) - o_D. $$

Letting $C_{K_0} = \max(\log(\frac{M_\#}{|K_0|}) + o_D, o_D - \log(\frac{M_b}{|\Omega|}))$, we get a uniform control of the cell unknown $\ell_{K_0}$:

$$ \|\ell_{K_0}\|_{L^\infty(K_0)} \le C_{K_0}. $$

Now, we can estimate the neighbouring face unknowns. Let $\sigma \in \mathcal{E}_{K_0}$ be a fixed face. The bound on $|\underline{\ell}_D|^2_{1,D}$ implies that $\|\ell_{K_0} - \ell_\sigma\|_{L^2(\sigma)}^2 \le h_\sigma C_\#$. Since $\ell_{K_0}$ is polynomial, one has

$$ \|\ell_\sigma\|_{L^2(\sigma)} \le \|\ell_{K_0}\|_{L^2(\sigma)} + \sqrt{h_\sigma C_\#} \le c_2 \|\ell_{K_0}\|_{L^\infty(\sigma)} + \sqrt{h_\sigma C_\#} \le c_2 C_{K_0} + \sqrt{h_\sigma C_\#}. $$

Therefore, letting $C_\sigma = c_\infty(c_2 C_{K_0} + \sqrt{h_D C_\#})$ we get a uniform control of the face unknown:

$$ \|\ell_\sigma\|_{L^\infty(\sigma)} \le C_\sigma. $$

On the other hand, let $L \in \mathcal{M}$ be a cell such that $\sigma \in \mathcal{E}_L$. As previously, one has $\|\ell_L - \ell_\sigma\|_{L^2(\sigma)}^2 \le h_\sigma C_\# \le h_D C_\#$, therefore

$$ \|\ell_L\|_{L^\infty(\sigma)} \le c_\infty \|\ell_L\|_{L^2(\sigma)} \le c_\infty \left( \|\ell_\sigma\|_{L^2(\sigma)} + \sqrt{h_D C_\#} \right) \le c_\infty \left( c_2 \|\ell_\sigma\|_{L^\infty(\sigma)} + \sqrt{h_D C_\#} \right). $$

Thus, assuming that $\|\ell_\sigma\|_{L^\infty(\sigma)}$ is controlled, one gets a uniform estimate on the trace of $\ell_L$ on $\sigma$. But, using the control of the oscillations (4.30), we deduce that

$$ \|\ell_L\|_{L^\infty(L)} \le \|\ell_L\|_{L^\infty(\sigma)} + o_D. $$

All in all, we can get a uniform estimate on any cell or face unknown. Because of the connectivity of the mesh, for any cell $K$ (respectively, face $\sigma$) there is a finite sequence of components of $\ell_D$, denoted $(\ell_k)_{k=0,...,m}$, starting at $\ell_0 = \ell_{K_0}$ and finishing at $\ell_m = \ell_K$ (respectively, $\ell_m = \ell_\sigma$) so that $\ell_K$ (resp. $\ell_\sigma$) can be connected to $\ell_{K_0}$. Given $i$, we denote by $C_i$ the $L^\infty$-norm of $\ell_i$ (either on a face or on a cell, depending on the nature of $\ell_i$). According to the previous
estimates, if $i$ is even (and corresponds to a cell unknown), one has

$$C_{i+1} \leq c_{\infty} \left( c_2 C_i + \sqrt{h_D C_{\#}} \right),$$

and else, if $i$ is odd (and corresponds to a face unknown), we have

$$C_{i+1} \leq c_{\infty} \left( c_2 C_i + \sqrt{h_D C_{\#}} \right) + o_D.$$

At the end, one can deduce that for any $0 \le i \le m-1$,

$$C_{i+1} \le o_D + c_{\infty} \sqrt{h_D C_{\#}} + c_{\infty} c_2 C_i.$$

Using a comparison with an arithmetico-geometric sequence, we infer that

$$C_m \le (c_\infty c_2)^m |C_0 - \gamma| + |\gamma| \quad \text{with} \quad \gamma = \frac{o_D + c_\infty \sqrt{h_D C_{\#}}}{1 - c_\infty c_2}. \qquad (4.31)$$

Therefore, $C_m$ depends only on $C_{K_0}$, $o_D$, $h_D$, $c_2$, $c_\infty$, $C_\#$ and $m$. Iterating this process for each unknown, since there is a finite number of cells and faces (at most $|M| + |\mathcal{E}|$), we deduce that there exists a positive constant $C$ depending on $C_{K_0}$, $o_D$, $h_D$, $c_2$, $c_\infty$, $C_\#$ and $|M| + |\mathcal{E}|$ such that

$$\forall K \in \mathcal{M}, \quad \|\ell_K\|_{L^\infty(K)} \le C \quad \text{and} \quad \forall \sigma \in \mathcal{E}, \quad \|\ell_\sigma\|_{L^\infty(\sigma)} \le C.$$

To conclude and get a bound on $\|\ell_D\|$, we just use the fact that the local (cell and face) $L^2$-norms are dominated by the local $L^\infty$-norms. $\square$

**Remark 24** (Non-uniform bound). Note that the estimate obtained in Lemma 9 involves quantities which strongly depend on the mesh, such as the meshsize $h_D$ or the number of cells and faces. Therefore, this estimation is not uniform with respect to the mesh.

We now state an existence result, which holds only for positive $\varepsilon$. The proof strongly relies on the methodology developed in Chapter 1 in order to show the existence result of Theorem 1 for the nonlinear HFV scheme.

**Theorem 10** (Existence of solutions to the nonlinear scheme (4.20)). Assume that the stabilisation parameter $\varepsilon$ in (4.18) is positive. Then, there exists at least one solution $(\underline{u}_D^n)_{n \ge 1}$ to the scheme (4.20). The associated densities $(\underline{u}_D^n)_{n \ge 1}$ are positive.

*Proof.* The proof proceeds by induction. Let $n \in \mathbb{N}$, and assume that $u_M^n$ is well defined, following the definition (4.20b) (if $n \ge 1$) or (4.20c) (if $n = 0$). We now prove the existence of a solution $\underline{\ell}_D^{n+1} \in \underline{V}_D^k$ to (4.20a). In order to simplify the computations, instead of looking for the logarithmic potential, we will seek for the discrete quasi-Fermi potential $\underline{w}_D^{n+1} = \underline{\ell}_D^{n+1} + \underline{\phi}_D - \log(c_{nl}^M) \underline{1}_D$.

To do so, we first notice that, given any $\underline{w}_D \in \underline{V}_D^k$ (in the following, we denote by $\underline{\ell}_D = \underline{w}_D - \underline{\phi}_D + \log(c_{nl}^M) \underline{1}_D$ the discrete logarithm potential and by $\underline{u}_D$ the associated density), the application

$$\underline{v}_D \mapsto \int_{\Omega} \frac{u_M - u_M^n}{\Delta t} v_M + T_D(\underline{\ell}_D, \underline{w}_D, \underline{v}_D)$$

is a continuous linear form on $\underline{V}_D^k$. Therefore, by Riesz representation theorem, there exists a
unique $\underline{G}_D(\underline{w}_D) \in \underline{V}_D^k$ such that

$$ \forall \underline{v}_D \in \underline{V}_D^k, \langle \underline{G}_D(\underline{w}_D), \underline{v}_D \rangle = \int_{\Omega} \frac{u_M - u_M^n}{\Delta t} \underline{v}_M + \mathcal{T}_D(\underline{\ell}_D, \underline{w}_D, \underline{v}_D). $$

On the other hand, since $\underline{V}_D^k$ is a finite-dimensional vector space, one can show that $\underline{w}_D \mapsto \underline{G}_D(\underline{w}_D)$ is a continuous mapping on $\underline{V}_D^k$. Notice that, for any discrete quasi-Fermi potential $\underline{w}_D \in \underline{V}_D^k$ satisfying $\underline{G}_D(\underline{w}_D) = \underline{0}_D$, the corresponding discrete logarithm potential $\underline{\ell}_D = \underline{w}_D - \underline{\phi}_D + \log(c_{hl}^M)\underline{1}_D$ is a solution to (4.20a). Thus, our goal from now on is to prove that $\underline{G}_D$ vanishes on $\underline{V}_D^k$.

We introduce a regularisation of $\underline{G}_D$: given any $\mu > 0$, we define an application $\underline{G}_D^\mu$ as

$$ \underline{G}_D^\mu : \underline{w}_D \mapsto \underline{G}_D(\underline{w}_D) + \mu \underline{w}_D. $$

By definition of $\underline{G}_D$, one has

$$ \begin{aligned} \langle \underline{G}_D^\mu(\underline{w}_D), \underline{w}_D \rangle &= \langle \underline{G}_D(\underline{w}_D), \underline{w}_D \rangle + \mu \|\underline{w}_D\|^2 \\ &= \int_{\Omega} \frac{u_M - u_M^n}{\Delta t} \underline{w}_M + \mathcal{T}_D(\underline{\ell}_D, \underline{w}_D, \underline{w}_D) + \mu \|\underline{w}_D\|^2. \end{aligned} $$

As it was already shown in the proof of Proposition 15 (equation (4.28)), by convexity of $\Phi_1$ one has

$$ \int_{\Omega} \frac{u_M - u_M^n}{\Delta t} \underline{w}_M \geq \frac{\underline{E}(\underline{w}_D) - \underline{E}^n}{\Delta t}, $$

where the relative entropies are defined by

$$ \underline{E}(\underline{w}_D) = \int_{\Omega} u_M^{\infty} \Phi_1 \left( \frac{u_M}{u_M^{\infty}} \right) \quad \text{and} \quad \underline{E}^n = \int_{\Omega} u_M^{\infty} \Phi_1 \left( \frac{u_M^n}{u_M^{\infty}} \right). $$

As already mentioned, since $\Phi_1$ is a non-negative function, these two quantities are non-negative. Moreover, since $u_M^n$ is piecewise $C^\infty$ if $n \ge 1$ and $\int_{\Omega} u_M^0 \log(u_M^0) = \int_{\Omega} u_M^{in} \log(u_M^{in}) < \infty$, one has $\underline{E}^n < \infty$. Note that it may occur that $\underline{E}^n = 0$ (which is equivalent to $u_M^n = u_M^\infty$), in which case $\underline{\ell}_D = \underline{\ell}_D^\infty$ is the unique solution to (4.20a) by mass preservation and entropy relation (4.27). In the following, we therefore assume that $\underline{E}^n > 0$. The previous identities and the non-negativity of the dissipation and entropy imply that

$$ \begin{aligned} \langle \underline{G}_D^\mu(\underline{w}_D), \underline{w}_D \rangle &\ge \frac{\underline{E}(\underline{w}_D) - \underline{E}^n}{\Delta t} + \mathcal{T}_D(\underline{\ell}_D, \underline{w}_D, \underline{w}_D) + \mu \|\underline{w}_D\|^2 \\ &\ge \mu \|\underline{w}_D\|^2 - \frac{\underline{E}^n}{\Delta t}. \end{aligned} \quad (4.32) $$

Thus, for any $\underline{w}_D \in \underline{V}_D^k$ such that $\|\underline{w}_D\| \ge \sqrt{\frac{\underline{E}^n}{\mu\Delta t}}$, one has $\langle \underline{G}_D^\mu(\underline{w}_D), \underline{w}_D \rangle \ge 0$. Therefore, according to Lemma 1 (see also [119, Section 9.1]), there exists a least one $\underline{w}_D^\mu \in \underline{V}_D^k$ such that

$$ \underline{G}_D^\mu(\underline{w}_D^\mu) = \underline{0}_D \quad \text{and} \quad \|\underline{w}_D^\mu\| \le \sqrt{\frac{\underline{E}^n}{\mu\Delta t}}. \quad (4.33) $$
Now, using (4.32) with $\underline{w}_D^\mu$, we get

$$ \frac{\mathbb{E}(\underline{w}_D^\mu)}{\Delta t} + \mathcal{T}_D(\underline{\ell}_D^\mu, \underline{w}_D^\mu, \underline{w}_D^\mu) + \mu \|\underline{w}_D^\mu\|^2 \le \frac{\mathbb{E}^n}{\Delta t}, $$

so $\mathcal{T}_D(\underline{\ell}_D^\mu, \underline{w}_D^\mu, \underline{w}_D^\mu) \le \frac{\mathbb{E}^n}{\Delta t}$. Therefore, recalling the definition (4.18)-(4.19) of $\mathcal{T}_D$, alongside with the fact that $\varepsilon > 0$ and the global stability result (4.10) for $a_D$, we then infer that

$$ \varepsilon \alpha_b \lambda_b h_b^{k+2} |\underline{w}_D^\mu|^2_{1,D} \le \frac{\mathbb{E}^n}{\Delta t} \quad \text{with } h_b = \min_{K \in \mathcal{M}} h_K. $$

On the one hand, by definition of the discrete logarithm potential $\underline{\ell}_D^\mu$ associated to $\underline{w}_D^\mu$ and the estimate (4.25) on $|\underline{\phi}_D|_{1,D}$, it holds

$$ |\underline{\ell}_D^\mu|_{1,D} = |\underline{w}_D^\mu - \underline{\phi}_D + \log(c_{nl}^M) \underline{1}_D|_{1,D} \le \sqrt{\frac{\mathbb{E}^n}{\varepsilon \alpha_b \lambda_b h_b^{k+2} \Delta t}} + \phi_\# \quad (4.34) $$

On the other hand, by definition of $\underline{\mathcal{G}}_D^\mu$, one first infers

$$ \begin{aligned} 0 = \langle \underline{\mathcal{G}}_D^\mu(\underline{w}_D^\mu), \underline{1}_D \rangle &= \int_{\Omega} \frac{u_M^\mu - u_M^n}{\Delta t} + \mathcal{T}_D(\underline{\ell}_D^\mu, \underline{w}_D^\mu, \underline{1}_D) + \mu \langle \underline{w}_D^\mu, \underline{1}_D \rangle \\ &= \int_{\Omega} \frac{u_M^\mu - u_M^n}{\Delta t} + \mu \langle \underline{w}_D^\mu, \underline{1}_D \rangle. \end{aligned} $$

Second, recalling the estimate (4.33) on $\underline{w}_D^\mu$ and using the Cauchy-Schwarz inequality, one gets

$$ \left| \int_{\Omega} (u_M^\mu - u_M^n) \right| \le \mu \Delta t \| \underline{w}_D^\mu \| \| \underline{1}_D \| \le \sqrt{\mu} \sqrt{\Delta t \mathbb{E}^n \| \underline{1}_D \|}. $$

Thus, letting $M^n = \int_{\Omega} u_M^n > 0$ and $\mu_n = \frac{(M^n)^2}{4\Delta t \mathbb{E}^n \|\underline{1}_D\|} > 0$, for all $0 < \mu < \mu_n$ one has

$$ \frac{M^n}{2} \le \int_{\Omega} e^{\underline{\ell}_M^\mu} \le \frac{3M^n}{2}. \quad (4.35) $$

Therefore, by estimates (4.34) and (4.35), one can apply Lemma 9 to $\underline{\ell}_D^\mu$ with $M_b = \frac{M^n}{2}$, $M_\# = \frac{3M^n}{2}$ and $C_\# = \left( \sqrt{\frac{\mathbb{E}^n}{\varepsilon \alpha_b \lambda_b h_b^{k+2} \Delta t}} + \phi_\# \right)^2$ (note that these three constants do not depend on $\mu$): there exists a constant $C > 0$ which does not depend on $\mu$ such that

$$ \forall \mu \in ]0, \mu_n[, \quad \|\underline{\ell}_D^\mu\| \le C. $$

Then, by compactness, there exists $\underline{\ell}_D^{n+1} \in \underline{V}_D^k$ such that, up to extraction, $\underline{\ell}_D^\mu \to \underline{\ell}_D^{n+1}$ when $\mu \to 0$. On the other hand, $\underline{\mathcal{G}}_D^\mu$ tends to $\underline{\mathcal{G}}_D$ as $\mu$ tends to 0. Therefore, letting $\underline{w}_D^{n+1} = \underline{\ell}_D^{n+1} + \underline{\phi}_D - \log(c_{nl}^M) \underline{1}_D$, we have $\underline{0}_D = \underline{\mathcal{G}}_D^\mu(\underline{w}_D^\mu) \to \underline{\mathcal{G}}_D(\underline{w}_D^{n+1})$ as $\mu \to 0$, which implies that

$$ \underline{\mathcal{G}}_D(\underline{w}_D^{n+1}) = \underline{0}_D. $$
It follows that $\underline{\ell}_{\underline{D}}^{n+1}$ is a solution to (4.20a). $\square$

**Remark 25** (Uniqueness of the solution). As for the nonlinear HFV scheme of Chapter 1, the existence result only asserts that solutions exist. However, the question of uniqueness is still open currently. A possible approach to show such a result could be to consider the relative entropy of a solution with respect to another solution, and show that this quantity vanishes.

Last, we state a result which indicates that the nonlinear scheme has a good long-time behaviour, in the sense that the solution converges in time towards the discrete thermal equilibrium. Note that this result is relatively weaker compared to the one of Proposition 14, since it does not give any indication about the speed of convergence.

**Proposition 16** (Long-time behaviour of the nonlinear scheme). Assume that the stabilisation parameter $\epsilon$ in (4.18) is positive, and let $(\underline{\ell}_{\underline{D}}^n)_{n \ge 1}$ be a solution to the scheme (4.20). Then, the solution converges towards the discrete logarithm potential at equilibrium as time tends to $+\infty$:

$$ \underline{\ell}_{\underline{D}}^n \underset{n \to \infty}{\longrightarrow} \underline{\ell}_{\underline{D}}^{\infty} \text{ in } \underline{V}_{\underline{D}}^k. \qquad (4.36) $$

Consequently, the corresponding discrete density $(u_M^n)_{n \ge 1}$ converges to $u_M^\infty$ in $L^\infty(\Omega)$.

*Proof.* First, note that by the entropy relation (4.27), one has

$$ \sum_{n \ge 1} \mathbb{D}^n \le \sum_{n \ge 1} \frac{\mathbb{E}^{n-1} - \mathbb{E}^n}{\Delta t} \le \frac{\mathbb{E}^0}{\Delta t}. $$

Thus, according to the definition of the discrete dissipation $\mathbb{D}^n$, alongside with the definition (4.19) of $T_D$, the fact that $\epsilon > 0$ and the global stability result (4.10) on $a_D$, we infer that

$$ \sum_{n \ge 1} |\underline{w}_{\underline{D}}^n|_{1,D}^2 \le \frac{\mathbb{E}^0}{\epsilon h_b^{k+2} \alpha_b \lambda_b \Delta t}, $$

where $h_b = \min_{K \in \mathcal{M}} h_K > 0$. This implies in particular that

$$ \forall n \ge 1, \quad |\underline{w}_{\underline{D}}^n|_{1,D}^2 \le \frac{\mathbb{E}^0}{\epsilon h_b^{k+2} \alpha_b \lambda_b \Delta t} \quad \text{and} \quad |\underline{w}_{\underline{D}}^n|_{1,D} \underset{n \to \infty}{\longrightarrow} 0. \qquad (4.37) $$

Let $n \ge 1$. By definition of $\underline{\ell}_{\underline{D}}^n$ and the bound on $\underline{\phi}_{\underline{D}}$, one has $|\underline{\ell}_{\underline{D}}^n|_{1,D} \le \underline{\phi}_{\#} + \sqrt{\frac{\mathbb{E}^0}{\epsilon h_b^{k+2} \alpha_b \lambda_b \Delta t}}$. On the other hand, by the mass preservation identity (4.26), we have $\int_{\Omega} e^{\underline{\ell}_{\underline{M}}^n} = M > 0$. Therefore, one can apply Lemma 9 to $\underline{\ell}_{\underline{D}}^n$: there exists a positive constant $C$ (which depends on $\epsilon$, $D$, $k$, $\Delta t$ and $M$, but not on $n$) such that

$$ \forall n \ge 1, \quad \|\underline{\ell}_{\underline{D}}^n\| \le C. $$

Therefore, by compactness, there exists $\underline{\ell}_{\underline{D}} \in \underline{V}_{\underline{D}}^k$ such that, up to extraction,

$$ \lim_{n \to \infty} \underline{\ell}_{\underline{D}}^n = \underline{\ell}_{\underline{D}} \text{ in } \underline{V}_{\underline{D}}^k. $$

On the other hand, by continuity of $|\cdot|_{1,D}$ on $\underline{V}_{\underline{D}}^k$, by definition of $\underline{w}_{\underline{D}}^n$ and (4.37), we have

$$ |\underline{\ell}_{\underline{D}} + \underline{\phi}_{\underline{D}}|_{1,D} = 0. $$
This means that there exists $a \in \mathbb{R}$ such that $\underline{\ell}_D + \underline{\phi}_D = a \underline{1}_D$. Moreover, thanks to mass preservation, it holds

$$M = \int_{\Omega} e^{\underline{\ell}_{M}^{n}} \underset{n \to \infty}{\longrightarrow} \int_{\Omega} e^{\underline{\ell}_{M}} \quad \text{so that} \quad \int_{\Omega} e^{\underline{\ell}_{M}} = M.$$

Hence, one has $a = \log(c_{nl}^M)$, which implies that $\underline{\ell}_D = \log(c_{nl}^M) \underline{1}_D - \underline{\phi}_D = \underline{\ell}_D^\infty$. By uniqueness of the limit, we get the convergence of the whole sequence $(\underline{\ell}_D^n)_{n \ge 1}$ towards $\underline{\ell}_D^\infty$. This implies in particular that $\underline{\ell}_M^n \underset{n \to \infty}{\longrightarrow} \underline{\ell}_M^\infty$ in $L^\infty(\Omega)$ by norm equivalence in finite-dimensional vector spaces. Then, by the mean value theorem we deduce that

$$\|u_M^n - u_M^\infty\|_{L^\infty(\Omega)} = \|e^{\underline{\ell}_M^n} - e^{\underline{\ell}_M^\infty}\|_{L^\infty(\Omega)} \le e^{\max(\|e_M^n\|_{L^\infty(\Omega)}, \|e_M^\infty\|_{L^\infty(\Omega)})} \|e_M^n - e_M^\infty\|_{L^\infty(\Omega)},$$

which implies the convergence of the densities. $\square$

**Remark 26** (Estimation of the rate of convergence). The previous result states a convergence of the potential towards the corresponding equilibrium potential. It follows that $u_M^n \underset{n \to \infty}{\longrightarrow} u_M^\infty$ in $L^\infty(\Omega)$. This result is relatively weaker than other long-time behaviour results presented in this manuscript. Indeed, there is no indication on the speed of convergence, even on a fixed mesh. However, note that the numerical results of Section 4.4.4 indicate that the convergence is as expected exponentially fast, and that the rate of convergence is similar to the continuous one. In order to state exponential convergence towards the equilibrium, one has to establish that the discrete dissipation controls the discrete entropy. Such an estimate (with a dependence on the meshsize) can be obtained for Dirichlet boundary conditions, using the same trick (comparison of the Boltzmann entropy with the $H^1$-norm of the quasi-Fermi potential, see (2.58)) as in Chapter 2. In order to get a uniform rate of convergence (with respect to the meshsize), one has to use further a discrete Logarithmic-Sobolev inequality. This strategy, already advocated for the low-order schemes of Chapter 1, is the subject of an ongoing work.

## 4.4 Numerical results

In this section, we assert the properties of the schemes on some test-cases. We also aim at investigating the relevance of using (nonlinear) high-order schemes in terms of computational efficiency (accuracy for a given computational cost). The tests considered below are all set in the domain $\Omega = ]0, 1]^2$, and are the same tests as in Chapter 1, to which we refer for more detailed descriptions.

Given a (face) degree $k$, the exponential fitting scheme (4.14) will be referred to as expf_k, the nonlinear scheme (4.20) as nlhho_k, whereas the HFV scheme of Chapter 1 will be denoted nlhfv. Note that nlhho_0 hinges on affine cell unknowns, whereas the cell unknowns for nlhfv are constant: these two schemes hence do not coincide (see Remark 22), and nlhho_0 is expected to be more costly.

In the sequel, we use the following stabilisation parameters: $\eta_{nl} = \eta_l = 1$. Regarding the parameter $\varepsilon$ we will use the value $\varepsilon = 1$. However, in some situations we will compare the two values $\varepsilon = 1$ and $\varepsilon = 0$. The nonlinear scheme with $\varepsilon = 0$ will then be denoted nlhho_k_0 for better readability.

The computations presented in this chapter were run on a laptop equipped with an Intel Core i7-9850H processor clocked at 2.60GHz and 32Gb of RAM.
### 4.4.1 Implementation

Here, we discuss some implementation aspects of the schemes. For a detailed introduction to the implementation of mixed-order HHO methods for standard linear diffusion problems (with a different consistent part), we refer the reader to [83, Chapter 8]. A different presentation (for equal-order methods) can be found in [97, Appendix B].

#### General setting

The HHO implementation here makes use of monomial basis functions for both the cell and face unknowns. Note that such a choice is known to introduce some numerical instabilities for large values of $k$. In the sequel, we will thus restrict the schemes to $k \le 3$ in order to avoid as most as possible these instabilities. Keep in mind that the use of orthonormal basis functions could improve the results presented below (and especially the convergence of the Newton method used for the nonlinear scheme), and should be one of the next future improvements. The quadrature formulas used in our implementation are based on the Dunavant rules [115], which implies to subdivide the cells into triangles. For the tests below, we use mesh families with cells star-shaped with respect to their barycenters, therefore the subdivisions correspond to the pyramidal submesh introduced in the HFV context (see Chapter 1). To deal with integrals of non-polynomial functions (for both schemes), we use quadrature formulas of order $2k + 5$. We performed a few tests (not shown here) with higher-order formulas and did not observe significant changes. Last but not least, the local computations are performed sequentially. One could expect a significant gain in terms of computational time parallelising these latter.

#### Exponential fitting scheme

For the linear exponential fitting scheme, the implementation follows the classical HHO strategy for linear diffusion problems. We directly solve for the discrete Slotboom variables $(\rho_D^n)_{n \ge 1}$. As standard for skeletal methods, we do not solve the full linear system but perform static condensation, which allows one to locally eliminate the cell unknowns (see [97, Appendix B.3.2]). Since the scheme relies on the same LHS matrix at each time step, we perform once and for all an LU decomposition of the matrix at the beginning of the computation. At each time step, the resolution is then inexpensive (the RHS has to be updated, but only through a matrix-vector product).

The implementation detailed here is relatively naive. In particular, we do not address the main questions which were highlighted in Chapter 1 about the (harmonic) averaging of $\omega$ (which is linked to the question of the quadrature formulas for the high-order scheme) and the preconditioning of the system (which was equivalent to choose to solve the system in the density variable $u$ for the HFV scheme). These questions should be addressed in future works. However, in view of the results obtained in Chapter 1 for the HFV exponential fitting scheme, we believe these potential improvements will have no effect on the lack of positivity of this scheme.

#### Nonlinear scheme

The numerical scheme (4.20) requires to solve a nonlinear system of equations at each time step. For a given $n \in \mathbb{N}$, one wants to find $\underline{\ell}_D^{n+1} \in \underline{V}_D^k$ solution to the nonlinear scheme (4.20): this scheme can be written as the equation

$$ \underline{G}_D^{n, \Delta t}(\underline{\ell}_D^{n+1}) = \underline{0}_D, $$
where $\underline{G}^{n,\Delta t}_D : \underline{V}_D^k \to \underline{V}_D^k$ is a regular (nonlinear) vector field. To find a zero of $\underline{G}^{n,\Delta t}_D$, we use a Newton method.

In practice, the use of a naive method without any adaptation proved not to be enough to compute a solution in general. In order to get a robust implementation, which can handle various data and meshes, we use a few techniques. For further use, we denote by $\|\underline{\ell}_D\|_{l^\infty}$ the $l^\infty$-norm of the coefficients of $\underline{\ell}_D$ in the polynomial bases (cells and faces) used to implement the scheme. This $\|\cdot\|_{l^\infty}$ is a norm on $\underline{V}_D^k$ which is easily (and at very low cost) computable in practice. To fix the ideas and notation, the Newton method is defined as follows: given an initialisation $\underline{\ell}_{D,0} \in \underline{V}_D^k$ and a time step $0 < \delta t \le \Delta t$, one defines a sequence $\underline{\ell}_{D,i}$ of elements of $\underline{V}_D^k$ such that

$$J_{\underline{\ell}_{D,i}}^{n,\delta t}(\underline{\ell}_{D,i+1} - \underline{\ell}_{D,i}) = -\underline{G}_{D}^{n,\delta t}(\underline{\ell}_{D,i}), \quad (4.38)$$

where $\underline{G}^{n,\delta t}_D$ is the vector field associated to the nonlinear scheme (4.20) with time step $\delta t$ instead of $\Delta t$, and $J_{\underline{\ell}_{D,i}}^{n,\delta t}$ is the differential (Jacobian in practice in the implementation) of $\underline{G}^{n,\delta t}_D$ at the point $\underline{\ell}_{D,i}$. Note that in practice, we do not solve this linear system, but perform static condensation (see [97, Appendix B.3.2]) in order to eliminate (locally) the cell unknowns. The resulting linear system is called "condensed system" in what follows. We discuss below the main tricks deployed to reach robustness in the implementation of the Newton method.

**Stopping criterion.** In order to have information about the convergence of the Newton method, we consider two quantities: the relative norm of the residual $r_{i+1}$ and the norm of the objective function $g_i$ defined as

$$r_{i+1} = \frac{\|\underline{\ell}_{D,i+1} - \underline{\ell}_{D,i}\|_{l^\infty}}{\|\underline{\ell}_{D,i}\|_{l^\infty}} \quad \text{and} \quad g_i = \|\underline{G}_D^{n,\delta t}(\underline{\ell}_{D,i})\|_{l^\infty}.$$

We consider that the Newton method has converged when either

$$(r_{i+1} \le 0.1 \times \text{tol}) \text{ or } (r_{i+1} \le \text{tol} \text{ and } g_i \le \text{tol})$$

with $\text{tol} = 5.10^{-10}$, in which case we set $\underline{\ell}_{D}^{n+1} = \underline{\ell}_{D,i+1}$. On the other hand, if this criterion is not met at $i = 50$, the method is considered as non-convergent (and we then proceed to a time step reduction, see below). In practice for the test of this chapter, we never reached $i = 50$, either because the method converged or because of a loop break, as described below.

**Loop break for large $\underline{\ell}_D$.** The computation of $\underline{G}^{n,\delta t}_D(\underline{\ell}_{D,i})$ and $J_{\underline{\ell}_{D,i}}^{n,\delta t}$ implies punctual evaluations of $e^{\underline{\ell}_{K,i}}$ and $e^{\underline{\ell}_{\sigma,i}}$ (for $K$ a cell and $\sigma$ a face) in the quadrature formulas. Such computations can lead to severe numerical issues if the values at the quadrature nodes are too large. Therefore, we declare that $\underline{\ell}_{D,i}$ is too large for the computations if there exists a cell quadrature node $x_{K,q} \in \bar{K}$ or a face quadrature node $x_{\sigma,q} \in \bar{\sigma}$ such that

$$|\underline{\ell}_{K,i}(x_{K,q})| \ge 100 \quad \text{or} \quad |\underline{\ell}_{\sigma,i}(x_{\sigma,q})| \ge 100.$$

In such a case, the method is immediately considered as non-convergent, and we proceed with a time step reduction (see below). Note that the choice of the value 100 allows one to compute densities $u_D$ on a range from $10^{-43}$ to $10^{43}$, and hence should not be a significant restriction in practice. In the numerical simulations presented in Sections 4.4.2-4.4.4, the use of this loop-breaking procedure is absolutely necessary in order to avoid some "explosion" of the method which leads to evaluations of too large quantities and crash of the code. Moreover, we also add a loop break if the linear solver does not
perform a successful resolution of the condensed linear system associated to (4.38), which corresponds to situations where either $J_{\ell_{D,i}}^{n,\delta t}$ or its condensed counterpart are not invertible. Such situations occur in practice, essentially on very coarse meshes.

**Adaptive time stepping.** The previous strategies can lead to resolution failures for a given time step $\delta t$. If the Newton method does not converge, we try to compute the solution for a smaller time step $0.5 \times \delta t$. On the other hand, if the method converges we use for the subsequent time step the value $2 \times \delta t$. The maximal time step allowed is the initial time step, denoted by $\Delta t$. In practice, the scheme may perform numerous time step reductions at the beginning (early times) of the computation.

**Initialisation by truncation and filtration.** As for any Newton method, the question of the initialisation is fundamental in order to get a robust implementation. It appears that, for $n \ge 1$, the natural initialisation $\underline{\ell}_{D,0} = \underline{\ell}_D^n$ is satisfactory when used with the adaptative time stepping strategy. However, for $n = 0$, such a choice is not possible, since $\underline{\ell}_D^0$ does not exist in general if $u^{in}$ vanishes locally or is too small (see Remark 21). A first way of tackling this problem (in the spirit of the techniques already used for the low-order schemes of the previous Chapters) is to define a truncated initial logarithm potential $\tilde{\ell}^{in}$ as

$$ \tilde{\ell}^{in} = \log(\max(u^{in}, 10^{-8})), $$

and to initialise the Newton method with a polynomial interpolation $\tilde{\ell}_D^0 \in \underline{V}_D^k$ of $\tilde{\ell}^{in}$. In fact, such a strategy induces another issue: $\tilde{\ell}_D^0$ exhibits strong oscillations in the regions where the truncation is performed (this also holds true if $u^{in}$ is not continuous on $\Omega$ near its discontinuities, as in Section 4.4.2). These oscillations usually make the method diverge, even with extremely small time step. Therefore, we eventually initialise the method with a “filtered” discrete unknown $\underline{\ell}_D^0 \in \underline{V}_D^k$ with no oscillation, which corresponds to a projection of $\tilde{\ell}_D^0$ to piecewise constants (on cells and faces):

$$ \forall K \in \mathcal{M}, \ell_K^0 = \Pi_K^0(\tilde{\ell}_K^0) = \frac{1}{|K|} \int_K \tilde{\ell}_K^0 \quad \text{and} \quad \forall \sigma \in \mathcal{E}, \ell_\sigma^0 = \Pi_\sigma^0(\tilde{\ell}_\sigma^0) = \frac{1}{|\sigma|} \int_\sigma \tilde{\ell}_\sigma^0. $$

In practice, using $\underline{\ell}_{D,0} = \underline{\ell}_D^0$ as the first initialisation (when $n = 0$) yields convergent Newton methods for all tests presented below. In particular, the use of filtered initial discrete data seems crucial for high-order schemes ($k \ge 1$). For the lowest order version of the scheme ($k = 0$), the use of $\tilde{\ell}_D^0$ as the first initialisation often yields convergent Newton methods (up to time step reduction). Note that in order to define the initialisation on the faces, one has to assume more regularity on $u^{in}$ than the initial hypothesis.

Of course, the chosen values for the stopping criterion and the thresholds are arbitrary and could be modified. However, the set of values advocated here makes the scheme robust enough so as to compute solutions for all the test-cases in this chapter.

**Remark 27 (Potentials vs. densities).** One of the main differences between the present scheme and the low-order ones from Chapter (1) lies in the fact that we use the potential $\ell$ as our unknown, whereas the density $u$ was used in the low-order schemes. As a consequence, our stopping criterion only provides information on $\ell$, while we are interested in the reconstructed density. Moreover, we only take into account the coefficients of the polynomials (through the norm $\|\cdot\|_{\infty}$), but such a norm does not give much information about the effective behaviour of the unknown. A more relevant stopping criterion
could be to consider the residual in terms of densities

$$ \frac{\|u_{M,i+1} - u_{M,i}\|_{L^2(\Omega)}}{\|u_{M,i}\|_{L^2(\Omega)}} $$

(and analogous definition for the face unknowns) in order to ensure a satisfying accuracy on $u$. However, the main issue of such a criterion is its evaluation cost. Thus, in this work we choose to use a purely algebraic stopping criterion on $\ell$, whose cost is marginal. The testing of other stopping criteria will be the subject of future investigations.

Last, recall that for the HFV schemes, we used a loop-breaking strategy: when the computed density had non-positive values, we performed a time step reduction. Here, such a situation cannot occur since $\ell$ takes values in $\mathbb{R}$, but this apparent absence of negativity on the density is in fact pernicious. Indeed, situations where $\ell_D$ takes very negative values are the counterpart of a negative $u$ for the low-order schemes, and lead to divergent Newton methods. The main difficulty here lies in the design of a relevant criterion in order to avoid this divergence.

## 4.4.2 Positivity

This first section is dedicated to assessing discrete positivity preservation. For the test considered here, we set the advection field and the diffusion tensor to

$$ \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$. We perform the simulation on the time interval $[0, 5.10^{-4}]$ with $\Delta t = 10^{-5}$ on a refined tilted hexagonal-dominant mesh (4192 cells and 12512 edges). The computed discrete densities are denoted by $(u_D^n)_{1 \le n \le N_f}$ and $(u_D^{\omega,n})_{1 \le n \le 50}$. Note that the situation $N_f > 50$ may occur if the nonlinear scheme has to perform time step adaptation. In this section, we also compare our results with the linear HMM scheme of [19] which was studied in Chapter 1.

In Table 4.1, we collect the minimal values reached by the discrete solutions. The values of “mincells” are defined by

$$ \min \left\{ \frac{1}{|K|} \int_K u_M^n | K \in \mathcal{M}, 1 \le n \le N_f \right\} \quad \text{and} \quad \min \left\{ \frac{1}{|K|} \int_K u_M^{\omega,n} | K \in \mathcal{M}, 1 \le n \le 50 \right\}, $$

for, respectively, the nonlinear scheme and the exponential fitting scheme. The values of “mincellQN” are the minimal values taken by the densities at the cell quadrature nodes. Analogous definitions hold for the faces. The values of “#resol” correspond to the number of linear systems solved during the computation. Note that the size of these systems depends on the value of $k$, so it is not a relevant information to compare the cost of the schemes for different values of $k$. Last, the “computing time” is the total time needed to compute the discrete solution (it includes the precomputational steps, such as the computation of the matrices representing $G_K^k$). Recall
<table><caption>Table 4.1 – Positivity of discrete solutions.</caption><thead><tr><th>scheme</th><th>computing time (in s)</th><th>#resol</th><th>mincells</th><th>minfaces</th><th>mincellQN</th><th>minfaceQN</th></tr></thead><tbody><tr><td>nlhfv</td><td>1.77e+01</td><td>175</td><td>9.93e-04</td><td>7.36e-04</td><td>9.93e-04</td><td>7.36e-04</td></tr><tr><td>nlhho_0</td><td>7.17e+01</td><td>224</td><td>1.00e-03</td><td>1.01e-03</td><td>2.41e-06</td><td>1.01e-03</td></tr><tr><td>nlhho_1</td><td>4.13e+02</td><td>248</td><td>6.65e-04</td><td>2.05e-05</td><td>1.78e-04</td><td>3.57e-08</td></tr><tr><td>nlhho_2</td><td>1.45e+03</td><td>251</td><td>9.50e-04</td><td>5.99e-04</td><td>2.67e-07</td><td>1.06e-05</td></tr><tr><td>nlhho_3</td><td>3.87e+03</td><td>254</td><td>9.85e-04</td><td>8.58e-04</td><td>1.10e-05</td><td>1.79e-05</td></tr><tr><td>HMM</td><td>2.20e-01</td><td>50</td><td>-5.00e-03</td><td>-7.74e-02</td><td>-5.00e-03</td><td>-7.74e-02</td></tr><tr><td>expf_0</td><td>5.66e-01</td><td>50</td><td>1.02e-03</td><td>1.89e-03</td><td>-3.78e-01</td><td>1.89e-03</td></tr><tr><td>expf_1</td><td>2.23e+00</td><td>50</td><td>-1.29e-02</td><td>-2.40e-01</td><td>-4.91e-01</td><td>-3.71e-01</td></tr><tr><td>expf_2</td><td>6.34e+00</td><td>50</td><td>-6.14e-03</td><td>-1.02e-01</td><td>-5.08e-01</td><td>-5.35e-01</td></tr><tr><td>expf_3</td><td>1.53e+01</td><td>50</td><td>-3.24e-04</td><td>-1.02e-02</td><td>-5.52e-01</td><td>-4.05e-01</td></tr></tbody></table>

that the HMM and the exponential fitting schemes are linear (with corresponding matrices not depending on time), whence their extremely low cost compared to the nonlinear schemes. However, when an LU decomposition is unaffordable and an iterative solver has to be used instead, nlhho_k is “only” five times more costly than expf_k.

The results of Table 4.1 first indicate that, as expected, all nonlinear schemes preserve the positivity of the discrete solutions. On the other hand, none of the linear schemes preserves positivity on the whole domain $\Omega$. In fact, except for expf_0, all linear schemes studied here do not even preserve the average positivity on each cell, in the sense that there exists $K_0 \in \mathcal{M}$ and $n_0$ an integer such that

$$ \int_{K_0} u_{M}^{\omega, n_0} < 0 \quad (\text{or } u_{K_0}^{n_0} < 0 \text{ for the HMM scheme}). $$

Moreover, it is interesting to note that the positivity violation peak (which can be approximated as $|mincellQN|$ and $|minfaceQN|$) increases as $k$ increases, while in average (values of mincells and minfaces) the lack of positivity becomes smaller as the order increases.

At this stage, it is worth pointing out the fact that quantifying the negativity of the solution is much more difficult for high-order schemes, since it is not possible to “count” the number of negative values (which corresponds to the degrees of freedom for low-order schemes). While the “mincellQN” value gives information about the minimum value reached on the whole domain, it does not give any indication about the measure of the set $\{x \in \bar{\Omega} | u_{M}^{\omega,n}(x) < 0\}$ where the discrete cell unknown takes negative values. The same remark applies to “mincells”. As an attempt to provide an idea of the size of this set, we display in Table 4.2 the number of cells with negative average over the whole simulation, defined as the cardinal of the set

$$ \{(K, n) \in \mathcal{M} \times [[1; 50]] | \int_K u_{M}^{\omega, n} < 0\}. $$

These data reveal that the higher the order, the smaller the size of the negative-values set (except for expf_0).

The previous observations seem to indicate a competition between two phenomena for the linear methods. As $k$ increases, the accuracy is improved, and therefore the discrete solution becomes closer to the continuous one. Hence, in average, high-order schemes compute solutions with smaller area of negativity. However, high values of $k$ induce larger oscillations for the
<table><caption>Table 4.2 – Number of negative cell averages</caption><thead><tr><th>scheme</th><th>HMM</th><th>expf_0</th><th>expf_1</th><th>expf_2</th><th>expf_3</th></tr></thead><tbody><tr><td># cells with negative average</td><td>90</td><td>0</td><td>824</td><td>136</td><td>1</td></tr></tbody></table>

polynomial solution: the computed solution takes negative values on smaller sets but the undershoots become bigger as $k$ increases. At the end, it seems that there is no hope to get a positive discrete solution on the whole domain $\Omega$ with a linear method on general meshes.

**Remark 28** (An accuracy criterion taking into account positivity). The previous observations suggest that, for applications where preserving the positivity of the solution is an essential feature, the accuracy of the scheme should not simply be defined as an $LP$ distance between the discrete solution $(u_M)$ and the exact one $(u)$. We believe that a relevant criterion in order to take into account both “classical accuracy” (distance between $u_M$ and $u$) and positivity is to look at the relative Boltzmann entropy (or other kinds of “relative $\Phi$-entropies” as defined in [33]) with respect to the exact solution, that is

$$Err(u_M) = \int_{\Omega} u \Phi_1 \left( \frac{u_M}{u} \right), \qquad (4.39)$$

where $\Phi_1(s) = s \log(s) - s + 1$ for $s > 0$, $\Phi_1(0) = 1$ and $\Phi(s)$ with large values for $s < 0$. The interest of such a definition is twofold. First, the negativity of $u_M$ is penalized. Second, if $u_M$ is positive and $\int_{\Omega} u_M = \int_{\Omega} u$ (which is the case in practice for problems with homogeneous Neumann boundary conditions), by Csiszár-Kullback inequality (see for example [51, Lemma 5.6]), one has

$$\|u_M - u\|_{L^1(\Omega)} \le \sqrt{\|u\|_{L^1(\Omega)} Err(u_M)}.$$

### 4.4.3 Convergence, accuracy and efficiency

Here, we are interested in the convergence as $(h_D, \Delta t) \to (0, 0)$ of the nonlinear schemes.

We consider a test-case with known exact solution. We set the advective potential and diffusion tensor to

$$\phi(x, y) = -x \quad \text{and} \quad \Lambda = \begin{pmatrix} l_x & 0 \\ 0 & 1 \end{pmatrix}$$

for $l_x > 0$. The exact solution is then given by

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

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$. Here, our experiments are performed using $l_x = 1$ and $C_1 = 10^{-1}$.

To get information about the accuracy of the schemes, we compute the discrete solution on the time interval $[0, 0.1]$, and we denote by $(u_D^n)_{1 \le n \le N_f}$ the corresponding discrete density. Then, we compute the relative $L_t^2(L_x^2)$ error on the solution and on the gradient of the solution, defined as

$$\frac{\sqrt{\sum_{n=1}^{N_f} \delta t^n \|u_M^n - u(t^n, \cdot)\|_{L^2(\Omega)}^2}}{\|u\|_{L_t^2(L_x^2)}} \quad \text{and} \quad \frac{\sqrt{\sum_{n=1}^{N_f} \delta t^n \|\mathbb{G}_M(u_D^n) - \nabla u(t^n, \cdot)\|_{L^2(\Omega)}^2}}{\|\nabla u\|_{L_t^2(L_x^2)}}$$

where $\delta t^n = t^n - t^{n-1}$ (note that $\sum_{1 \le n \le N_f} \delta t^n = 0.1$) and $\mathbb{G}_M(u_D^n)$ is the discrete gradient of the
densities. This discrete operator is defined by mimicking the continuous relation $\nabla u = e^\ell \nabla \ell$ at the discrete level: $G_M(u_D^n)$ is a piecewise continuous vector field satisfying, for any cell $K \in \mathcal{M}$,

$$G_M(u_D^n)|_K = e^{\ell K} G_K^K(\ell_K) \text{ in } K.$$

As already mentioned, $L^2$-norms (in space) are computed using quadrature formulas of order $2k + 5$, and we use the approximation

$$\|u\|_{L_i^2(L_x^2)} \approx \sqrt{\sum_{n=1}^{N_f} \delta t^n \|u(t^n, \cdot)\|_{L^2(\Omega)}^2}$$

for the exact solution norms (with an analogous definition for $\|\nabla u\|_{L_i^2(L_x^2)}$). Notice that, with the chosen definitions, we do not take into account the time $t = 0$. To plot the error graphs, we perform our simulations on a triangular mesh family ($D_i$)$_{1 \le i \le 5}$, such that $h_{D_i}/h_{D_{i+1}} = 2$. Since the time discretisation is of order one, on the i-th mesh of the family, we use an initial (and maximal) time step of

$$\Delta t_i = \frac{\Delta t_0}{2^{i(k+2)}}$$

where $\Delta t_0 = 0.05$. The goal of this choice is to fit the time step according to the expected order of convergence of the scheme (in space).

<figure><img src="image_43.png" /><figcaption>Figure 4.1 – Accuracy of transient solutions. Relative errors on triangular meshes.</figcaption></figure>

In Figure 4.1, we plot the errors as functions of the meshsize $h_D$ for different values of $k$. We see that, as one could expect, the scheme $\text{nlhho}_k$ converges at order $k + 1$ in energy-norm, and $k + 2$ in $L^2$ norm of the density.

We now study efficiency, that is to say accuracy for a given computational cost. In Figure 4.2, we plot the errors as functions of the computing walltime. It is remarkable to see that, even with a low-order discretisation in time, significant efficiency gains can be obtained using a large value of $k$. The gain is expected to be even bigger after parallelising the local computations. Of course, the use of higher-order time-stepping methods should also lead to significant gains, and this
<figure><img src="image_44.png" /><figcaption>Figure 4.2 – Accuracy vs. computational cost. Relative errors on triangular meshes.</figcaption></figure>

will be investigated in future works.

**Remark 29** (High-order schemes in time and space). The extension of the nonlinear scheme (4.20) to arbitrary orders in time and space is a rather natural goal in order to achieve optimal efficiency. However, even with a discretisation of order 2 (like for example BDF2), there is currently no successful approach retaining the discrete entropy structure. Since this structure is the cornerstone of the analysis (including the existence of solutions), it is important to preserve it. Some numerical investigations on nonlinear entropic TPFA schemes for diffusive problems with order two time discretisation (BDF2) have been performed in [89, Chapter 3], and indicate that such a time discretisation could be a good candidate, even with respect to long-time behaviour (see [89, Chapter 3.4.4]).

As a last accuracy test, we compare the nonlinear HHO scheme with the HFV scheme of Chapter 1. For the sake of completeness, we also include the scheme `nlhho_2_0` (i.e. with $\epsilon = 0$ in (4.18)) in our comparison. For the `nlhfv` scheme, we use the same definition of the errors as in Chapter 1. The results are presented on Figure 4.3. First, we can see that `nlhho_2` and `nlhho_2_0` have the same behaviour. Tests with other values of $k$, not shown here, indicate that the influence of $\epsilon$ (0 or 1) on the accuracy of the scheme is not noticeable. Regarding the efficiency, the interest of using high-order schemes seems clear (in particular when one is interested in the accuracy on the gradient). For the accuracy on the density, the trend is quite clear and shows that for (relative) errors smaller than $10^{-5}$, the high-order schemes are more efficient. Note that the comparison presented here between the HFV and HHO schemes is not completely fair. Indeed, the HFV scheme has been implemented within an independent code, where a lot of optimisations (including precomputations and static condensation with explicit inverse matrices) are performed. Moreover, the question of the stopping criterion for the Newton method (see Remark 27) implies that it is difficult to compare "similar" versions.

Overall, the tests presented here confirm that the nonlinear high-order schemes allow for a consequent gain of efficiency compared with low-order schemes.
<figure><img src="image_45.png" /><figcaption>Figure 4.3 – Accuracy: nlhfv vs. nlhho. Relative errors on triangular meshes.</figcaption></figure>

### 4.4.4 Discrete long-time behaviour

We are now interested in the long-time behaviour of discrete solutions. We use the same test-case as before, but with an anisotropic tensor: we set $l_x = 10^{-2}$. The corresponding steady-state is

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

We compute the solution on the time interval $[0, 350]$, with $\Delta t = 10^{-1}$, on two Kershaw meshes of sizes 0.02 and 0.006. In Figure 4.4, we show the evolution along time of the $L^1$ distance between $u_D^n$ and $u^\infty$, computed as

$$\int_{\Omega} |u_M^n - u^\infty|,$$

for the nlhho_k schemes, for $k \in \{0, 1, 2\}$. We observe the exponential convergence towards the thermal equilibrium for each value of k, until some precision is reached. For $k \ge 1$, the rates of convergence are similar to the exact one ($\alpha$), and do not depend on the size of the mesh. For $k = 0$, the rate of convergence towards equilibrium differs a bit from $\alpha$ on the coarsest mesh,
<figure><img src="image_46.png" /><figcaption>Figure 4.4 – Long-time behaviour of discrete solutions. Comparison of the long-time behaviour of the solutions to the nonlinear HHO schemes on Kershaw meshes for $T_f = 350$ and $\Delta t = 0.1$.</figcaption></figure>

but these two rates seem to coincide on sufficiently refined meshes. Note that for the test-case considered, $\phi \in \mathbb{P}^1(\Omega)$, therefore $\phi_M = \phi$. It follows that $u_M^\infty = u^\infty$, hence the very low saturation threshold.

For the exponential fitting scheme, the distance between $u_D^n$ and $u^\infty$ is computed as

$$ \int_{\Omega} |u_{M}^{\omega,n} - u^{\infty}| $$

To assess the behaviour of the exponential fitting scheme, we last compare the results for three

<figure><img src="image_47.png" /><figcaption>Figure 4.5 – Long-time behaviour of discrete solutions. Comparison of the long-time behaviour of the solutions to the nonlinear HHO schemes on Kershaw meshes for $T_f = 350$ and $\Delta t = 0.1$.</figcaption></figure>
schemes with face unknowns of degree $k = 1$: nlhho_1, nlhho_1_0 and expf_1.
The results are presented in Figure 4.5. As previously, one can note that nlhho_1 and nlhho_1_0
exhibit an extremely similar behaviour. Regarding the exponential fitting scheme, as stated
in Proposition 14, the discrete density converges exponentially fast towards the equilibrium.
Moreover, as for the nonlinear schemes, the decay rate is similar to the continuous one.

4.5 Conclusion

In this chapter, we have introduced two arbitrary-order schemes for linear anisotropic advection-diffusion equations on general meshes. The first one is a linear scheme which follows the exponential fitting strategy, whereas the second is a nonlinear scheme which follows the low-order philosophy [56, 52, 70]. We have proved that both schemes admit solutions, have an entropy structure and preserve the long-time behaviour of the continuous problem. Moreover, the solutions to the nonlinear scheme are positive by construction. We have validated these theoretical results on numerical test-cases. We have also highlighted the lack of positivity of the linear method, which justifies the use of (more costly) nonlinear methods. In the meantime, the use of nonlinear schemes with polynomial unknowns of sufficient degree proves to improve the efficiency (accuracy at given computational cost) with respect to low-order nonlinear methods. These results confirm the benefits of using high-order nonlinear schemes in order to get reliable approximation of diffusive problems. Future directions of this work include a deeper analysis of the nonlinear scheme, as well as the development of similar schemes for more complex and nonlinear problems.
# Conclusion and perspectives

## Outline of the current chapter

<table><tr><td><b>Further analysis of the schemes</b></td><td><b>161</b></td></tr><tr><td><b>More efficient and robust implementations</b></td><td><b>162</b></td></tr><tr><td><b>Applications to other complex problems</b></td><td><b>163</b></td></tr><tr><td><b>Numerical comparisons with existing methods</b></td><td><b>164</b></td></tr><tr><td><b>High-order discretisations in time</b></td><td><b>165</b></td></tr></table>

The work presented in this thesis shows that high-order nonlinear structure-preserving methods for anisotropic diffusive problems on general meshes can be used successfully. The study performed in Chapter 1 shows that non-linear methods are needed in order to preserve positivity of solutions to anisotropic problems on general meshes, and introduces a nonlinear structure-preserving HFV method. This method is then used on complex and nonlinear anisotropic problems in Chapter 2 and 3, and allows one to approximate solutions to the drift-diffusion system (3) while preserving the structural bounds (implying in particular positivity) on the densities and the long-time behaviour. These encouraging results motivate the design of high-order schemes with similar structure-preserving features. Using the quasi-Fermi potentials as key variables, we successfully adapt the nonlinear HFV discretisation to an HHO one in Chapter 4. The scheme introduced exhibits an optimal order of accuracy, while preserving the positivity and the long-time behaviour of the solutions. The implementation proposed enjoys a robustness similar to the HFV one, and the use of sufficiently high-order schemes gives a consequent gain of computational efficiency with respect to the HFV scheme.

Thus, this work paves the way to various extensions and applications to other models. We discuss some of them below, regarding the analysis, the implementation and the applications.

## Further analysis of the schemes

Regarding the analysis of the schemes, the most interesting challenges lie in the way of getting uniform (with respect to the meshsize) estimates needed in order to get long-time behaviour and convergence results.

The first natural question about the analysis of the HHO scheme is to understand how to get uniform (with respect to the meshsize) estimates. Indeed, such estimates are needed in order to prove both convergence of the scheme and uniform long-time behaviour. Performing an analysis of the long-time behaviour of the nonlinear HHO scheme similar to the one performed in Chapter 1 is a first important step and an ongoing work. The main issues is to get an high-order
counterpart of Proposition 6. The main difference when it comes to the high-order scheme is the difference between potential (polynomials) and densities: one has to deal with both in order to get the appropriate estimates. Another important task then is to adapt discrete functional inequalities of Section 1.A to the case of high-order unknowns.

The convergence analysis of the non-linear schemes for advection-diffusion equations is an important perspective of our work, for both HFV and HHO schemes. The analysis of the HFV scheme was not proved in this thesis, however the long-time behaviour analysis gives some uniform estimates on $H^1$ discrete semi-norms in Proposition 6. Such estimates, combined with the compactness toolbox [6] are the key ideas to show the convergence of similar DDFV and VAG schemes [56, 52]. The adaptation of these arguments to the HFV framework should be possible. For the high-order scheme, one needs to get uniform estimates (see above) and handle non-polynomial density. It is not known whether the toolbox [6] can be used directly.

The convergence analysis of schemes for drift-diffusion systems are more difficult. Indeed, the coupling between the unknowns makes it difficult to get uniform estimates on the relevant quantities. Up to our knowledge, the convergence analysis of existing schemes for drift-diffusion systems relies on $L^\infty$ bounds on the densities and/or on the electrostatic potential. On the other hand, getting such estimates (uniform with respect to the meshsize) seems rather difficult for the methods considered here. A potential approach could be to adapt arguments used for the recent analysis of a continuous memristor model [170], which is essentially a drift-diffusion systems with three densities instead of two. The main interest of this analysis is the fact that it does not rely on $L^\infty$ estimates on the solution (which are not available because there are three species instead of two) but rather on a pure entropy/energy method. We believe that such a method could be used at the discrete level with our schemes. Note that the analysis of long-time behaviour (with an uniform rate of decay with respect to the meshsize) of HFV schemes for the drift-diffusion system should also be performed using similar ideas.

Another important task is related to the analysis of the nonlinear schemes in situations out of equilibrium. Such analysis is performed for a nonlinear $P^1$ finite element scheme [57] which enjoys the same structure as the nonlinear HFV scheme. The adaptation of this analysis to the nonlinear HFV scheme seems possible, but rather technical since one has to deal with non-conformal method. A related work could be the analysis of positivity-preserving schemes for stationary problems. The accuracy study of Chapter 1 suggests that the stationary version of the nonlinear HFV scheme is suited in order to approximate solutions to stationary problem, but up to now the analysis this scheme (or similar VAG or DDFV schemes) was not performed. A relevant analysis approach should be to use the relative entropy with respect to a lifting of the boundary condition in order to get entropy estimations. However, the main difficulty in order to analyse such a stationary scheme is related to the way of getting estimates on the entropy dissipation without the help of the time derivative term. Such estimates could be consequences of refined Log-Sobolev inequalities.

## More efficient and robust implementations

In terms of implementation, we are working on a couple of ways to improve the numerical behaviour of the schemes.

A first interesting question about the numerical robustness is the implementation of the Newton method. Indeed, some recent work [38, 24] suggests that the use of parametrisation (or change of unknowns) in the Newton method could strongly improve its robustness. For the HFV schemes, such strategy would mean be to solve the problem using either densities or quasi-Fermi potentials as primal unknowns. Note that the use of quasi-Fermi potentials as
unknowns is the solution implemented in the TPFA code **ddfemi** [104]. On the other hand, for the HHO scheme, since we can only use the coefficients of the polynomials representing the quasi-Fermi potentials, it is not straightforward to switch between densities and potentials. Another possibility could be to consider more subtle resolution method, using for example line search techniques. Some numerical test performed during the phase development of **hho_nl** with relaxed Newton steps indicated that a better robustness can be achieved with such strategy, but at the cost of many iterations. Of course, such numerical work could also be associated to a theoretical study, in order to ensure the behaviour of the Newton method (invertible Jacobian matrices, global convergence of the methods).

Another important feature which should be investigated is the question of mesh adaptation. In [143], it is suggested that an equidistribution of the entropy dissipation over the domain (which means, the same dissipation on every cells) leads to a better repartition of the errors, and therefore to a better efficiency. We believe that refining the mesh in areas where the dissipation is too big can be a way of improving the convergence of the Newton methods. Such questions are the topic of an ongoing work on TPFA schemes with collaborators from WIAS and Tyndall National Institute for semiconductors models studied in Appendices: we aim at understanding how to use the entropy and its dissipation in order to get coarser mesh without losing too much accuracy. Our ultimate goal is to simulate bigger devices than the ones studied in Appendix D. Of course, the use of HFV and HHO schemes should be a better solution when it comes to mesh adaptation, given the flexibility afforded by general meshes.

Last but not least, an important question which was not addressed during this PhD is the parallelisation of the computations for the hybrid schemes. Indeed, these methods are based on local contribution, and a lot of computations can be performed in parallel (static condensation, computation of every local quantity). Such implementation should allow one to get a way better efficiency for the scheme implemented, and could be part of the development of the ParaSkel++ [18] platform. Natural related development in this framework includes methods handling three dimensional domains and the use orthogonal polynomial basis instead of monomial ones.

## Applications to other complex problems

In term of non-academic problems, my work essentially focused on semiconductor models, and the main goal was to approximate (3). In view of the results of Chapter 4, an extension of the work performed in Chapters 2 and 4 to develop and analyse an high-order scheme for semiconductor models is a short-term objective in order to get more efficient structure-preserving schemes for this problem.

Of course, the methodology developed in this thesis is intended to be used in a variety of different settings and models. Here, we discuss three situations of interest. As above, the schemes developed could be either HFV ones (as a preliminary work) or HHO ones.

A first natural application is related to porous media problems, for which there are challenges related to anisotropy. In particular, there exists some HMM, DDFV and HHO based schemes [65, 71, 72, 5] for Peaceman problem. These schemes are shown to be convergent, but there is no theoretical guarantees that the computed concentrations are non-negative. Therefore, the use of the non-linear methods of this PhD could be a solution in order to get a more reliable approximation of the solutions.

Another field of application could be the wide family of cross-diffusion systems. These are systems of multiples unknowns, in which the diffusion coefficients for each unknown depends on the other unknowns. As for the drift-diffusion systems (which can be seen as simple prototypes
of general cross-diffusion systems), the developpement and analysis of TPFA schemes for cross-
diffusion is an active field of research [54, 147, 177]. On the other hand, there exists some
anisotropic cross-diffusion models obtained as rigorous limits of stochastic models [230, 117].
For such models, the use of the previous TPFA schemes is not possible: the use of HFV and
HHO schemes based on the work of this thesis should be a natural solution to get reliable
approximation methods. Of course, the use of high-order schemes should also be an interesting
solution to approximate such systems in an efficient way. Since these systems are generally
analysed thanks to boundedness by entropy methods [47, 173], the methodology used to analyse
our scheme should be adaptable without too much difficulty to cross-diffusion.

Another interesting question is the generalisation of the results presented in this thesis to
convection fields, which are not the gradient of potentials. Indeed, in the study performed here,
the only scheme which is able to handle such a general field is the HMM scheme [19] studied
in Chapter 1. However, this scheme does not preserve positivity, neither thermal equilibrium.
Moreover, its analysis relies on coercivity assumptions. From a more general point of view, the
main difficulty of such a situation in order to develop structure-preserving methods is that
the quasi-Fermi potentials, keystone of our structure preserving methods and the associated
analysis, do not exist anymore. A first direction in order to solve this issue could be to use a Hodge
decomposition of the field $V = \nabla\phi + V_0$, where $V_0$ is divergence-free, and treat the irrotational
part in a non-linear way by using a partial quasi-Fermi potential $\log(u) + \phi$. The advective (and
coercive) part induced by $V_0$ could then be discretised using the scheme of [19]. Computations
at the continuous level indicate that this strategy -which is in a way similar to the one used to
analyse non-coercive advection-diffusion problems [106]- allows one to get some estimates on
the entropy dissipation. On the other hand, at the discrete level the way of getting this estimate
is not straightforward and one has to understand how to perform the Hodge decomposition in
an appropriate way. Note that this situation is not purely academic, and can be used in practical
situations. In particular, there exists a series of work by Jochmann [163, 164, 165, 167, 166, 168]
about semiconductor models in which the electrostatic filed does not solve a Poisson equation.
Such a model correspond to a situation where the internal magnetic effects are not neglected
anymore: the electrostatic field is not the gradient of a potential anymore, but is given by Maxwell
equations. It could be interesting to develop positivity preserving high-order methods for such a
system, using the methodology introduced in this thesis for the convection-diffusion equation
alongside with some recent HHO methods for Maxwell equations [76]. Another interesting
direction for such work could be to investigate the asymptotic preserving features of schemes
with respect to the asymptotic analysis done by Jochman for the continuous model (long-time
behaviour and singular limit when the magnetic susceptibility tends to zero).

## Numerical comparisons with existing methods

As far as possible, we gave numerical comparison between different methods in order to
identify and exhibit the strengths of each method introduced.

The results of Chapter 4 indicate that the high-order structure-preserving method developed
in this thesis have a better efficiency that the nonlinear HFV method of Chapter 1. However,
keeping in mind that the methods developed here address the need to handle anisotropy and
substitute for TPFA schemes, a question of great importance is the comparison of the hybrid
methods with the TPFA ones. Indeed, in situations where both classes of methods are available
(isotropic problems on orthogonal meshes), the efficiency of the high-order schemes could be
a strong argument in favour of HHO schemes. In order to understand the performances of
each scheme, one need to perform tests on various nonlinear problems. In order to avoid the
limitations induced by the (first order accurate) time discretisation and give a fair comparison,
the ideal situation should be to use stationary problems.

Another interesting work could be an extension of the results presented in Chapter 3. Indeed,
we only focused on an isotropic drift-diffusion system, but one could ask if both HFV and DDFV
schemes have the same behaviour for the anisotropic problem (3) with general statistics. Ongoing
numerical investigations indicate that both schemes handle anisotropy for Boltzmann statistics.
Moreover, the influence of the intensity of the magnetic field on the rate of exponential decay
towards equilibrium (results in the same spirit as Figure 2.11b) is rather similar for both schemes.
In view of the results of Chapter 4, it would be interesting to also compare the low-order schemes
of Chapter 3 with a nonlinear structure-preserving HHO scheme for anisotropic drift-diffusion
systems. Such comparative work could be a first step in order to investigate the behaviour of
high-order structure-preserving schemes on complex problems and get reference results. On the
other hand, such a work could also be the opportunity to investigate the question of efficiency
(and the expected gain for the high-order methods) of the schemes on nonlinear problems.

**High-order discretisations in time**

A last important perspective of this work is related to the time-stepping method used. Indeed,
all schemes presented in this manuscript are based on a backward Euler discretisation in time.
Such a choice allows one to get discrete entropy structures relatively easily as a consequence
of the convexity of the entropies. On the other hand, since every schemes are of order at least
two in space, which leads to a rather large unbalanced space-time accuracy. This is even worse
for high-order schemes: the use of high-order methods for time discretisation should provide a
scheme higher efficiency.

The main difficulty when it comes to high-order methods in time is the way of getting the
entropy relation, which is the base of all the nonlinear schemes developed here. A first step
should be to test different numerical strategies in order to check both stability and efficiency of
the schemes. Such a comparative work was performed for TPFA schemes in [89], and shows that
BDF2 schemes provide good stability properties for entropic schemes. The analysis of such a
scheme could rely on the methodology used in [146] to analyse an entropic scheme using BDF2.
Another interesting direction could be the use of Runge-Kutta schemes following the work [176]
on semi-discrete problems. A possible interesting solution could be the use of Strong Stability
Preserving Methods, which are often used in the context of hyperbolic equations [151, 152].
Note that for the linear schemes used in this thesis, the use of a Crank–Nicolson discretisation
in time leads to an exact entropy-dissipation relation at the discrete level (with an equality as
in the continuous framework). However, getting a control of the entropy by the dissipation
(which then depends on two times) in this situation seems rather difficult. Moreover, such a time
discretisation is known to generate spurious oscillations which should degrade the positivity of
solutions.
# Semiconductor models with varying band edge energies: an overview

## Outline of the current chapter

<table><tr><td>A.1 Motivations and context</td><td>167</td></tr><tr><td>A.2 A comparison between numerical simulations and experimental results (Appendix B)</td><td>169</td></tr><tr><td>A.3 The question of thermodynamic consistency (Appendix C)</td><td>169</td></tr><tr><td>A.4 Application to the design of efficient LEDs (Appendix D)</td><td>170</td></tr><tr><td>A.5 Some mathematical issues raised by these works</td><td>171</td></tr></table>

## A.1 Motivations and context

In the Appendices, we focus on models with irregular convection fields. Such models are particularly useful in describing realistic devices, where the structure induces local fluctuations in internal energies. The prototype equation for modeling such a framework can be reduced to the linear problem:

$$ \partial_t u - \text{div}(\nabla u + u \nabla(\phi + E)) = 0, \qquad (A.1) $$

where $\phi$ is a regular potential (solving a Poisson equation in the context of drift-diffusion systems) and $E$ is a function known as the Band Edge Energy (BEE), which depends on the material (alloy) used in the semiconductor. In classical models (as in the drift-diffusions systems studied in Chapter 2 and 3), this energy is assumed to be constant. However, there are practical situations in which the BEEs are not constant and can even be irregular. Two cases of interest are:

- regions of the alloy with strong local fluctuations, referred to as quantum wells, where the BEEs exhibit extreme irregularity (essentially, a function taking different values at each atom);

- junctions of two layers made of different materials, known as heterojunctions, where the BEEs are piecewise constant.
Random alloy fluctuations significantly affect the electronic, optical, and transport properties of (In,Ga)N-based optoelectronic devices. To bridge the gap between macroscale drift-diffusion simulations and atomistic band-edge fluctuations, recently a multiscale framework was developed to integrate the macroscopic and microscopic worlds [221].

To understand how we can appropriately approximate the solutions of such problem, it is necessary to establish a framework for analysing the continuous equation (A.1). To the best of our knowledge, the most comparable situations in the literature are related to “distributional drifts” in stochastic problems [136, 137, 238, 161, 162]. In these works, the problems considered are driven by advection fields that lie in low-regularity Sobolev spaces, and the solutions are defined as limits of solutions to regularised problems. However, these advections are too regular for our study, as we aim to consider drifts that are singular measures, such as gradients of piecewise constant functions.

On the other hand, while the advection considered in (A.1) can be highly irregular, it has the specific structure of a gradient. By formally defining $\rho = \frac{u}{\omega}$, with $\omega = e^{-\phi-E} \in L^{\infty}(\Omega)$, we can express an equivalent formulation of (A.1) in terms of $\rho$, namely

$$ \omega \partial_t \rho - \text{div}(\omega \nabla \rho) = 0. \qquad (A.2) $$

Such an equation is classically well-posed, as the diffusion coefficient $\omega$ is bounded and uniformly positive on $\Omega$. Hence, one can define the concept of solution to (A.1) by employing the exponential fitting formulation (A.2). Interestingly, this notion of solution can also be obtained as a limit of regularised solutions.

This approach can be generalised to problems with non-linear statistics (written here for the stationary problem to fix ideas)

$$ -\text{div}(u \nabla(h(u) + \phi + E)) = f, $$

where $h$ is a function corresponding to the considered statistics (see, for example, the introduction of Chapter 2). By defining $\rho = \frac{e^{h(u)}}{\omega}$, one formally obtains an exponential fitting formulation:

$$ -\text{div}\left(\frac{g(\log(\rho) - \phi - E)}{\rho}\nabla\rho\right) = f, \qquad (A.3) $$

where $g$ is a positive function, the inverse of $h$. In this formulation, similar to the previous discussion, the gradient of the irregular term $E$ no longer appears.

At the numerical level, we propose a strategy based on the exponential fitting formulations and Two Point Flux Approximation (TPFA) finite volume methods. As for the HFV exponential fitting developed in Chapter 1, one of the main challenges lies in averaging the diffusion coefficient $\omega$ (or $g(\log(\rho) - \phi - E)$ for general statistics), which is irregular. It is worth noting that this issue has been addressed in the context of heterogeneous diffusion problems for TPFA schemes [121, Section 2.3.1]. Our numerical results (see Figure A.1) demonstrate the suitability of the exponential fitting for this type of problem, provided that an appropriate averaging technique is employed to handle the irregular BEE. In practice, it is necessary to use a harmonic mean for the averaging process.

The remainder of this appendix is organized as follows. First, we present the results of Appendices B-C-D, introducing the physical context of the problems and the issues addressed. Then, in Section A.5, we summarise various perspectives linked to this work, from both physical and mathematical points of view.
<figure><img src="image_48.png" /><figcaption>Figure A.1 – Accuracy of schemes. Boltzmann statistics and heterojunctions.</figcaption></figure>

## A.2 A comparison between numerical simulations and experimental results (Appendix B)

In Appendix B, we compare the numerical results obtained using our strategy with those obtained from commercial packages and experimental data.

Specifically, we consider an (In,Ga)N/GaN semiconductor system with quantum wells and introduce a model that incorporates random fluctuations in the alloy. This is achieved by computing appropriate band edge energies using the Localisation Landscape Theory (LLT) approach [11, 132, 224]. The LLT approach allows for the prediction of charge localisation regions and their corresponding energies. Our method differs from classical strategies implemented in commercial software, which employ a Schrödinger equation to incorporate quantum corrections. However, these Schrödinger-Poisson solvers have a significant computational cost, and currently, these commercial packages only support 1D simulations.

To assess the validity of our approach, we compare our simulations with results obtained using the commercial package nextnano [30], which employs a full Schrödinger-Poisson equation, as well as experimental data. In particular, we investigate the impact of the well order on device behaviour. Our results demonstrate good agreement between our simulations and those presented in [30]. Moreover, considering the random fluctuations is necessary to accurately reproduce the observed trends in the experiment.

## A.3 The question of thermodynamic consistency (Appendix C)

In Appendix C, we investigate certain issues related to the approximation of models with non-Boltzmann statistics. These discussions are motivated by observations resulting from the use of incorrect discretisation methods for solving problems with varying BEEs. However, they also apply to classical models with constant BEE.

These issues are well-known within the mathematical community and have led to numerous studies, ranging from the development of schemes handling general statistics to the comparison of such methods [265, 120, 26, 140, 50, 125, 127, 50]. However, in more applied scientific
communities, these questions are relatively under-documented.

The work presented in Appendix C aims to provide physics-based evidence that using inconsistent discretisations of the flux can significantly impact the results of simulations. The main novelty of this work lies in investigating the influence of discretisation on important physical behaviours, namely the recombination process and current-voltage (I-V) curves. The former is of major importance in LED applications, while the latter serves as a basis for engineers to model the device behaviour in an electrical circuit. The study conducted provides strong evidence that employing an inconsistent Scharfetter–Gummel flux in situations where non-Boltzmann statistics are used leads to non-physical results.

As a by-product of this study, we also introduce a new definition of thermodynamically consistent numerical flux (C.7). This notion is based on the relationship between the monotonicity of the quasi-Fermi potential and the sign of the flux. The introduced definition corresponds to a discrete version of this continuous feature. The main advantage of this notion is that, unlike the classical one<fnref n="1" />, it holds meaning even outside of thermal equilibrium. It appears that this notion is linked to the methodology used in [50] to analyse TPFA schemes for Blakemore statistics: the numerical fluxes are expressed as the product of a positive interface term by the discrete derivative of the quasi-Fermi potentials.

## A.4 Application to the design of efficient LEDs (Appendix D)

The last appendix focuses on a practical question: how to design efficient LEDs. We are specifically interested in LEDs operating in the deep-UV-C wavelength range (below 280 nm), which can be used for applications such as sterilisation or water purification filters. Aluminium gallium nitride ((Al,Ga)N) based LEDs are considered ideal candidates for such devices, as their emission wavelengths can be tuned across the entire UV spectrum without the need for toxic mercury. However, the use of (Al,Ga)N alloy is relatively new, and there is limited literature available compared to more classical (In,Ga)N-based optoelectronic systems. The main goal of our work is to understand the impact of random fluctuations in the alloy on LED efficiency.

To understand the behaviour and efficiency of LEDs, it is important to examine the recombination processes. In the situation under study, two main types are of interest:

- Radiative recombination, which involves electron transitions from the conduction band to the valence band accompanied by photon emission. This process is the primary source of light emission.

- Non-radiative Auger-Meitner recombination, which involves internal transfer of electrons (from the valence band to the conduction band) or of holes. This internal process can be viewed as an energy leak associated with the "emission" of phonons. Experimental results [148, 233] indicate that at high densities of charge carriers, this recombination process can physically degrade the device.

As a preliminary approach, we can consider that the ideal device will exhibit high radiative recombination and low non-radiative recombination. Our study focuses on these two types of recombination, comparing the device behaviour with and without random fluctuations.

Ultimately, our results suggest that random fluctuations increase both radiative and non-radiative recombination. To improve efficiency, careful device design is necessary to control the increase in non-radiative recombination. This can be achieved, for example, by using wider wells or increasing the number of wells in relevant areas.

<footnote n="1">1. If the numerical fluxes vanish, then the quasi-Fermi potentials are constant.</footnote>
## A.5 Some mathematical issues raised by these works

Regarding the work performed on the toy model, the numerical investigation suggests that it is possible to develop accurate schemes even for irregular BEEs. A useful first step before analysing schemes for full drift-diffusion systems could be to prove the convergence of the scheme for the toy problem equation. Convergence analysis for Boltzmann statistics is a classical approach using the compactness method since it only involves an equation with a heterogeneous diffusion coefficient. However, the case of non-Boltzmann statistics is more challenging. In fact, even for the continuous modified equation (A.3), we do not yet know if solutions exist. One possible first step could be to study the continuous equation and show, as in the linear case, that an equivalence holds between regularised solutions and exponentially fitted ones. Another possibility could be to directly prove the convergence of the numerical schemes using compactness. Such a result would imply, in particular, that a solution to the continuous problem exists. The techniques used in [50] could be relied upon for the analysis.

When it comes to simulating realistic devices, one of our medium-term objectives is to compare the numerical results obtained using our framework with experimental data that will soon be produced at the National Tyndall Institute. This comparative work will allow us to confirm the relevance of the approach used (LLT combined with the drift-diffusion system) in modelling complex devices. Once this confirmation has been obtained, it will be possible to apply our simulations to study optimal designs for LEDs.

The main hindrance to perform these comparisons is the cost of the simulations. Since we use an extremely refined mesh in the wells, our computations are not affordable for realistic-sized devices with a large number of wells. Furthermore, the adaptation between refined and coarse areas is challenging due to the orthogonality constraint on the (TPFA) admissible meshes. To address this issue, we are currently working on ways to locally coarsen the mesh while maintaining the same device behaviour. One of our goals is to find good heuristic criteria for mesh adaptation. A related question could be linked to homogenisation. Instead of using the exponential fitting approach to define a continuous solution, we could view the BEEs as oscillating functions and employ more classical homogenisation methods to obtain a well-posed limit equation. Such homogenised equation could be easier to discretise and would allow us to simulate larger devices with more quantum wells.

All the previous observations suggest that, according the promising results of Chapter 4, an interesting perspective could be to use HHO structure-preserving methods in order to simulate the devices under consideration. Regarding the irregularity of the BEEs, we believe that the non-linear strategies of the hybrid schemes developed in this thesis are adapted to handle it. Indeed, returning at the initial flux

$$J = u \nabla(\log(u) + \phi + E),$$

one can introduce a (partial) quasi-Fermi potential $w = \log(u) + E$, and expresses the flux as

$$J = e^{w-E} \nabla(w + \phi).$$

Here, the irregular BEE $E$ only appears in the diffusion coefficient. Hence, we can approximate the corresponding PDE using the same strategy as in Chapter 4, by discretising the potentials $w$ and $\phi$ as polynomials. From a numerical perspective, the main difficulty caused by the irregularity of the BBE is the heterogeneity of the diffusion coefficient. However, this feature is not an issue in theory when it comes to the HHO method used (since it can handle heterogeneous and anisotropic diffusion). Of course, from a practical point of view, as usual with exponential fitting, the main challenge lies in approximating the integrals $\int_K e^{-E}$ with irregular and oscillating
integrands.
# Impact of random alloy fluctuations on the carrier distribution in multi-color (In,Ga)N/GaN quantum well systems

## Outline of the current chapter

<table><tr><td>B.1 Introduction</td><td>174</td></tr><tr><td>B.2 Model MQW structures and literature experimental findings</td><td>175</td></tr><tr><td>B.3 Theoretical framework</td><td>178</td></tr><tr><td>B.3.1 Tight-binding energy landscape</td><td>178</td></tr><tr><td>B.3.2 Device simulation</td><td>179</td></tr><tr><td>B.4 Results</td><td>182</td></tr><tr><td>B.4.1 Continuum-based simulations of the carrier transport in (In,Ga)N-based LEDs</td><td>182</td></tr><tr><td>B.4.2 Impact of a random alloy fluctuations on the carrier transport in (In,Ga)N/GaN MQWs</td><td>185</td></tr><tr><td>B.5 Conclusions</td><td>186</td></tr></table>

This chapter corresponds to a submitted work [219], in collaboration with Michael O'Donovan, Patricio Farrell, Timo Streckenbach, Thomas Koprucki and Stefan Schulz.

In this work, we study the impact that random alloy fluctuations have on the distribution of electrons and holes across the active region of a (In,Ga)N/GaN multi-quantum well based light emitting diode (LED). To do so, an atomistic tight-binding model is employed to account for alloy fluctuations on a microscopic level and the resulting tight-binding energy landscape forms input to a quantum corrected drift-diffusion model. Here, quantum corrections are introduced via localization landscape theory
and we show that when neglecting alloy disorder our established theoretical frame-
work yields results very similar to commercial software packages that employ a
self-consistent Schrödinger-Poisson-drift-diffusion solver; this provides validation of
the developed quantum corrected transport model. Similar to experimental studies
in the literature, we have focused on a multi-quantum well system where two of the
three wells have the same In content while the third well differs in In content. By
changing the order of wells in this ‘multi-color’ quantum well structure and looking at
the relative radiative recombination rates of the different emitted wavelengths, we (i)
gain insight into the distribution of carriers in such a system and (ii) can compare our
findings to trends observed in experiment. Our results indicate that the distribution
of carriers depends significantly on the treatment of the quantum well microstructure.
For instance, when including random alloy fluctuations and quantum corrections in
the simulations, the calculated trends in the relative radiative recombination rates as a
function of the well ordering are consistent with previous experimental studies. How-
ever, the results from the widely employed virtual crystal approximation contradict
the experimental data. Our calculations clearly demonstrate that when accounting for
random alloy fluctuations in the simulations, no further ad-hoc modifications to the
transport model are required, in contrast to previous studies neglecting alloy disor-
der. Overall, our work highlights the importance of a careful and detailed theoretical
description of the carrier transport in an (In,Ga)N/GaN multi-quantum well system
to ultimately guide the design of the active region of III-N-based LED structures.

B.1 Introduction

At the heart of modern light emitting diodes (LEDs) operating in the blue to violet spectral region are (In,Ga)N/GaN multi-quantum well (MQW) systems [158]. While the efficiency of these LEDs is and can be very high, further efficiency gains will still directly reduce the cost of operating such LEDs. Moreover, extending efficient operation of (In,Ga)N-based LEDs into the green to red spectral range is a topic of current research interest [201, 240, 12, 116, 159]. To achieve all this, understanding the fundamental electronic, optical and transport properties of (In,Ga)N-based MQW systems is of central importance to guide design of (In,Ga)N-based LEDs with new and improved capabilities. While experimental and theoretical studies on the electronic [59, 75, 252] and optical [169, 95] properties of such systems have already revealed that these properties are significantly impacted by alloy fluctuation induced carrier localization effects, the impact of alloy disorder on the carrier transport has only been targeted recently [220, 221, 222, 203, 44, 194]. Here, studies are ranging from fully atomistic quantum mechanical approaches [220] up to modified continuum-based models [221, 222, 203, 44, 194].

Recently, we have developed a three-dimensional (3-D) multiscale simulation framework that connects atomistic tight-binding theory with a modified, quantum corrected drift-diffusion (DD) solver. The framework has been employed to investigate uni-polar carrier transport and it was found that alloy fluctuations result in an increase in carrier transport of electrons (in an *n*-i-*n* system) [221], but decreases transport in the case of holes (in a *p*-i-*p* systems) [222]. In the present work we extend this scheme to investigate the active region of (In,Ga)N-based MQW LED structures (thus *p*-i-*n* systems) and study how carriers distribute across the active region. In general, understanding the carrier distribution can help to guide maximizing the efficiency in an LED, since ideally the carriers shall be distributed evenly across the entire MQW region so that all QWs will contribute to emission [266]. However, previous experimental studies on
carrier distribution in (In,Ga)N/GaN MQW systems have indicated that mainly the well closest to the *p*-doped contact side contributes to the light emission process [145, 92, 199, 266]. These samples were specifically designed to gain insight into the carrier distribution inside the active region of an LED.

Overall, this has been attributed to a sequential filling of the QWs, resulting in a high hole density only in the *p*-side QW. To establish accurate carrier transport models the trends found in the experimental studies of Refs. [145, 92, 199, 266] need to be captured. Previous theoretical studies have reproduced the experimentally observed behaviour, however this required (i) treating bound carriers in a quantum mechanical picture, (ii) softening of the QW barrier interface to account for tunneling effects, (iii) distinguishing between continuum and bound carriers in the carrier transport model (multi-population model), and (iv) allowing for scattering between the different populations [236]. But, the impact of alloy disorder is basically neglected in this advanced but also complex carrier transport model.

In this paper we show that when employing our quantum corrected 3-D simulation framework that accounts for random alloy fluctuations, the experimentally observed trends are captured, without introducing for instance a multi-population scheme. This highlights that our developed solver presents an ideal starting point for future device design studies.

To highlight clearly the impact that random alloy fluctuations have on the carrier distribution in the active region of an (In,Ga)N-based LED, we use as a reference point a virtual crystal approximation (VCA) which effectively can be described by a 1-D model. The benefit of this is twofold. Firstly, this enables us to compare directly the outcomes of our quantum corrected model with results from 1-D commercial software simulations; commercial software packages often employ a standard Schrödinger-Poisson-DD solver, which is numerically very costly and therefore unfeasible in large 3-D transport simulations. This motivates the need for an alternative implementation of quantum corrections. Secondly, and building on this benchmark, alloy fluctuations can be included in the calculations, revealing clearly their impact on the results. Our studies show, and when using the same input parameter set, only the model accounting for random alloy fluctuations produces trends that are consistent with the experimental data. The widely employed VCA yields results that are in contradiction with the experimental data, thus indicating that radiative recombination stems mainly from the well *furthest away* from the *p*-side. Overall, this highlights (i) that alloy fluctuations are essential to achieve an accurate description of the carrier transport and (ii) have to be taken into account when theoretically guiding the design of energy efficient III-N light emitters.

The paper is organized as follows: in Section B.2 we outline the model structure used for calculations and briefly summarise some of the literature experimental data from Ref. [145]. The theoretical framework which we use is summarized in Section B.3. Our results are discussed in Section B.4. Finally Section B.5 presents our conclusions.

## B.2 Model MQW structures and literature experimental findings

To investigate the carrier distribution in (In,Ga)N/GaN MQW systems we proceed similar to experimental studies in the literature [145, 199] and target MQW systems where one of the wells in the MQW stack has a slightly higher In content compared to the remaining wells. In our case, we study MQW systems with three (In,Ga)N/GaN wells. Here two are $In_{0.1}Ga_{0.9}N$ ("shallow") wells and one is an $In_{0.125}Ga_{0.875}N$ ("deep") QW. These QWs are 3 nm wide and separated by 5 nm GaN barriers. The band edge profile of such a system along the transport (*c*)- direction, using a VCA, is shown in Fig. B.1 at a current density of 50 A/cm².
<figure><img src="image_49.png" /><figcaption>Figure B.1 – Conduction and valence band edges (black) along with the quasi-Fermi energies for electrons and holes (grey) in an (In,Ga)N/GaN multi-quantum well system described in virtual crystal approximation. The band edge profile and the quasi Fermi levels are shown at a current density of 50 A/cm². The leftmost (In,Ga)N quantum well contains 12.5% indium while the other two (In,Ga)N wells (centre and right) contain 10% indium.</figcaption></figure>
<figure><img src="image_50.png" /><figcaption>Figure B.2 – Schematic illustration of multi-quantum well system. The *n*-doped region is shown in cyan, the *p*-doped is in red and undoped regions are in grey. The quantum wells are numbered starting from the *n*-side.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure B.3 – Profile of (a) valence band edge energy, (b) conduction band edge energy, and (c) radiative recombination rate in the growth plane (c-plane) of an In$_{0.1}$Ga$_{0.9}$N quantum well; the current density is 50 A/cm$^2$ in all depicted figures. The slice displayed is the through the center well. The data are shown in all cases on a linear scale.</figcaption></figure>

In the following, we investigate the carrier transport properties in two settings: (i) on an atomistic level accounting for random alloy fluctuations and (ii) in the frame of a VCA thus neglecting alloy fluctuations. In the latter VCA case (ii), at a given z-position (along the c-direction), there is no variation in material properties within the growth plane (c-plane). This assumption is also made in the widely used 1-D transport simulations on (In,Ga)N MQWs.

To study the carrier distribution in MQW systems using the simulation settings (i) and (ii), we follow again the experimental approach e.g. presented in Ref. [145] and the deep QW is moved from the *n*-side (position 1 (*n*-side) in Fig. B.2) to the *p*-side (position 3 (*p*-side) in Fig. B.2). In the case of the random alloy structures, the same microscopic configuration is kept for each well and only the ordering is changed.

For each of these systems the ratio of radiative recombination from the shallow wells to the deep well is calculated using:

$$ \rho = \frac{R_{\Omega_S}^{RAD}}{R_{\Omega_D}^{RAD}} \qquad (B.1) $$
where

$$R_{\Omega_i}^{\text{RAD}} = \int_{\Omega_i} R^{\text{RAD}}(\mathbf{r})dV, \qquad (B.2)$$

is the total radiative recombination from the region $\Omega_i$. Here, $\Omega_D$ is the region containing the deep QW, $\Omega_S$ is the region containing the shallow wells (as there are two shallow QWs this is the union of the two shallow QW regions). The radiative recombination rate at position $\mathbf{r}$, $R_{\text{RAD}}(\mathbf{r})$, is discussed in further detail in section B.3.2. Since we are studying a system with three QWs, an even distribution of carriers across the MQWs would result in a ratio of $\rho = 2$. Previous experimental work on a similar system by Galler et. al. [145] found that $\rho$ was small (i.e. emission is dominated by the deep QW) only when the deep well was closest to the p-doped side of the MQW system (thus position 3 (p-side) in Fig.B.2). The authors conclude that holes are responsible for this behavior, and argue that they are mainly found in the p-side QW and not in wells further away from the p-side. As a consequence, the overall emission from the (In,Ga)N/GaN MQW system is dominated by the emission from this well closest to the p-doped region. In line with Ref. [145], we calculate $\rho$ at a current density of 50 A/cm$^2$, which allows us to compare the here predicted trends with the trends found in the experimental studies. The theoretical framework employed to gain insight into $\rho$ is discussed in the following section.

## B.3 Theoretical framework

In this section, we introduce the underlying (microscopic) theory of our multiscale simulations. We start in Section B.3.1 with the electronic structure model, an atomistic tight-binding (TB) model, and discuss the drift-diffusion approach in Section B.3.2. Since all these ingredients have been discussed in detail in Refs. [221, 222], we here give only a brief summary.

### B.3.1 Tight-binding energy landscape

In order to model the electronic structure of the above described (In,Ga)N MQW systems on an atomistic level, we employ a nearest-neighbour, sp$^3$ TB model [242]. Here, local strain and polarization effects are included using a valence force field model and local polarization theory, respectively [60].

To connect the TB model to a DD solver we proceed as follows. Firstly, the TB model is used to extract a potential energy landscape describing the MQW region of the device using an atomistic framework. To do so, at each atomic site in the three dimensional (3-D) supercell, a local TB Hamiltonian is constructed from the full TB Hamiltonian [75]. Subsequently, only the local TB Hamiltonian is diagonalized, yielding the conduction and valence band edge energy at each lattice site. These band edge energies include now already effects arising from alloy fluctuations and connected fluctuations in strain and built-in polarization field. The obtained 3-D confining energy landscape, after employing a Gaussian softening, forms the basis for our DD calculations. In previous studies we have investigated and discussed in detail the influence of the Gaussian softening on transport calculations for electrons and holes [221, 222]. Here, we choose a Gaussian broadening on the order of the GaN lattice constant, $\sigma_{c,v} = a^{\text{GaN}} = 0.3189$ nm, in all calculations. This value is large enough to average over a number of neighboring sites, while also small enough to retain fluctuations in the energy landscape.

To obtain an accurate description of carrier transport in (In,Ga)N-based LEDs, the DD equations, which will be discussed below, are often coupled with solving the Schrödinger equation to account for quantum corrections. Such a Schrödinger-Poisson solver is widely available in commercial software packages. However, in these packages it is largely restricted