# A hypocoercivity-exploiting stabilised finite element method for Kolmogorov equation

Zhaonan Dong, Emmanuil H Georgoulis, Philip J. Herbert

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

## To cite this version:

Zhaonan Dong, Emmanuil H Georgoulis, Philip J. Herbert. A hypocoercivity-exploiting stabilised finite element method for Kolmogorov equation. SIAM Journal on Numerical Analysis, 2025, 63 (3), pp.1105-1127. 10.1137/24M163373X. hal-04413333v2

HAL Id: hal-04413333

https://inria.hal.science/hal-04413333v2

Submitted on 27 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
A HYPOCOERCIVITY-EXPLOITING STABILISED

FINITE ELEMENT METHOD FOR KOLMOGOROV EQUATION*

ZHAONAN DONG†, EMMANUIL H. GEORGOULIS‡, AND PHILIP J. HERBERT§

**Abstract.** We propose a new stabilised finite element method for the classical Kolmogorov equation. The latter serves as a basic model problem for large classes of kinetic-type equations and, crucially, is characterised by degenerate diffusion. The stabilisation is constructed so that the resulting method admits a numerical hypocoercivity property, analogous to the corresponding property of the PDE problem. More specifically, the stabilisation is constructed so that spectral gap is possible in the resulting “stronger-than-energy” stabilisation norm, despite the degenerate nature of the diffusion in Kolmogorov, thereby the method has a provably robust behaviour as the “time” variable goes to infinity. We consider both a spatially discrete version of the stabilised finite element method and a fully discrete version, with the time discretisation realised by discontinuous Galerkin timestepping. Both stability and a priori error bounds are proven in all cases. Numerical experiments verify the theoretical findings.

**Key words.** Kolmogorov equation, hypocoercivity, Galerkin methods, finite element methods, stabilised finite elements

MSC codes. 65N30

**1. Introduction.** In the study of kinetic modelling, degenerate evolution problems are abundant. These may arise from the integration of stochastic processes, which upon integration typically give rise to Kolmogorov/Fokker-Planck type partial differential equations (PDEs). In many examples, diffusion (stemming from Brownian motion/randomness) is present in some of the spatial directions in the respective Kolmogorov/Fokker-Planck PDE. Quite interestingly, despite this, many such PDEs with degenerate diffusion give rise to initial/boundary value problems that are, nonetheless, characterised by convergence to long time equilibria. Villani in [18] coined the term hypocoercivity to signify the subtle property of certain degenerate evolution (integro)differenetal operators to yield dissipation in directions where no diffusion is explicitly present. Hypocoercivity has been formulated as a general framework in proving decay to equilibrium for kinetic equations in the seminal work [18] by Villani, who conceptualised and expanded key ideas proposed in the seminal works of Hérau & Nier in [10], Eckmann & Hairer in [6], Mouhot & Neumann [13] and others.

To highlight the role of degeneracy of diffusion with respect to long-time behaviour, we consider the classical Kolmogorov equation:

$$ (1.1) \qquad \mathcal{L}u \equiv u_t - u_{xx} + xu_y = f, \quad \text{in } (0, t_f] \times \Omega; $$

for some final time $t_f > 0$, with $\Omega \subset \mathbb{R}^2$, and a suitably smooth forcing $f: (0, \infty) \times \Omega \to \mathbb{R}$, along with appropriate initial and boundary data. (The above notation

<footnote>*Submitted to the editors DATE.</footnote>

<footnote>Funding: This work was funded by the EPSRC (grant number EP/W005840/1).</footnote>

<footnote>†1) Inria, 48 rue Barrault, 75647 Paris, France, 2) CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée, France (zhaonan.dong@inria.fr).</footnote>

<footnote>‡1) Maxwell Institute for Mathematical Sciences and Department of Mathematics, School of Mathematical and Computer Sciences, Heriot-Watt University, Edinburgh EH14 4AS, United Kingdom (E.Georgoulis@hw.ac.uk), 2) Department of Mathematics, School of Applied Mathematical and Physical Sciences, National Technical University of Athens, Zografou 15780, Greece & 3) IACMFORTH, Greece.</footnote>

<footnote>§Department of Mathematics, University of Sussex, Brighton BN1 9RF, United Kingdom. (P.Herbert@sussex.ac.uk)</footnote>
used above is typical in PDE literature; in kinetic applications, the operator $\mathcal{L}$ may appear as $\mathcal{L}f := f_t - f_{vv} + v f_x$ with $v$ denoting the particle velocity variable, $x$ the displacement/position and $f$ the respective probability density function.) Let $(\cdot, \cdot)_{L_2(\Omega)}$ the standard $L_2(\Omega)$-inner product, for some 'spatial' domain $\Omega \subset \mathbb{R}^2$, along with the respective norm $\|\cdot\|_{L_2(\Omega)}$. Assuming suitable boundary conditions, such that the skew-symmetric term $(xu_y, u)$ vanishes, the standard energy argument for (1.1) yields

$$ (1.2) \qquad \frac{1}{2} \frac{d}{dt} \|u(t)\|_{L_2(\Omega)}^2 + \|u_x\|_{L_2(\Omega)}^2 = (f, u)_{L_2(\Omega)} \le \frac{1}{2\epsilon} \|f\|_{L_2(\Omega)}^2 + \frac{\epsilon}{2} \|u\|_{L_2(\Omega)}^2, $$

for almost all $t \in (0, t_f]$ and all $\epsilon > 0$, where $t_f \in \mathbb{R}_+$ is some final time. Observing that a Poincaré-Friedrichs/spectral gap estimate of the form $\|u\|_{L_2(\Omega)} \le C\|\nabla u\|_{L_2(\Omega)}$ is not applicable, and since an estimate of the type $\|u\|_{L_2(\Omega)} \le C\|u_x\|_{L_2(\Omega)}$ may not hold, standard arguments involving Grönwall's Lemma give

$$ (1.3) \qquad \|u(t_f)\|_{L_2(\Omega)}^2 \le e^{\epsilon t_f} (\epsilon^{-1} \|f\|_{L_2(0,t_f;L_2(\Omega))}^2 + \|u_0\|_{L_2(\Omega)}^2), $$

for any $\epsilon > 0$, where we have used standard notation for Bochner spaces. The presence of the term $e^{\epsilon t_f}$ mean that the stability estimate (1.3) is not robust with respect to the final time.

The hypocoercivity of (1.1), however, can be manifested through non-standard energy arguments. More specifically, by considering different inner products of form $(\cdot, \cdot)_{\text{hc}} = (\cdot, \cdot)_{L_2(\Omega)} + (A\nabla\cdot, \nabla\cdot)_{L_2(\Omega)}$, where $A \in \mathbb{R}_{\text{symm}}^{2 \times 2}$ is a suitable positive semi-definite matrix, and the associated norm $\|\cdot\|_{\text{hc}} := \sqrt{(\cdot, \cdot)}_{\text{hc}}$, one arrives at an energy estimate of the form

$$ \frac{1}{2} \frac{d}{dt} \|u(t)\|_{\text{hc}}^2 + c (\|\nabla u\|_{L_2(\Omega)}^2 + \|u_x\|_{\text{hc}}^2) \le (f, u)_{\text{hc}}, $$

for some constant $c > 0$. Crucially, if $u$ is such that it satisfies a Poincaré-Friedrichs inequality, (e.g., satisfying suitable boundary conditions,) we can deduce

$$ \frac{1}{2} \frac{d}{dt} \|u(t)\|_{\text{hc}}^2 + c_{\text{hc}} \|u\|_{\text{hc}}^2 \le (f, u)_{\text{hc}}, $$

for some constant $c_{\text{hc}} > 0$ which depends on $A$, so that $\frac{d}{dt} \|u(t)\|_{\text{hc}}^2 + c_{\text{hc}} \|u\|_{\text{hc}}^2 \le c_{\text{hc}}^{-1} \|f\|_{\text{hc}}^2$, using standard arguments. Thus, Grönwall's Lemma gives

$$ (1.4) \qquad \|u(t_f)\|_{\text{hc}}^2 \le e^{-c_{\text{hc}}t_f} \|u_0\|_{\text{hc}}^2 + c_{\text{hc}}^{-1} \int_0^{t_f} e^{-c_{\text{hc}}(t_f-s)} \|f(s)\|_{\text{hc}}^2 ds. $$

This is reminiscent of the standard parabolic decay estimates. In particular, for $f = 0$, this provides decay as $t \to \infty$.

Retaining hypocoercive structures upon discretisation of respective PDE problems is desirable. This is because numerical methods admitting such structures would give rise to robust performance for long time computations. So far, only few works have been concerned with such numerical hypocoercivity properties. We mention [15] for a study of hypocoercivity for a central difference method discretising the Kolmogorov equation (1.1), using the abstract framework from [18, Theorem 18], and [5] where the approach from [15] was taken further to discuss the availability of spectral gap for central difference discretisations of the respective Fokker-Planck equation. Also, Foster, Lohéac & Tran [8] discuss the development of a Lagrangian-type splitting
method based on a carefully constructed similarity transformation and linear finite
elements over quasiuniform meshes, although no proof of decay to equilibrium is given
for the numerical method itself. Also, Bessemoulin-Chatard & Filbet [1] present design
principles for the construction of equilibrium-preserving finite volume methods for
nonlinear degenerate parabolic problems. Recently, Bessemoulin-Chatard, Herda &
Rey [2] presented an asymptotic-preserving in the diffusive limit finite volume scheme
for a class of one-dimensional kinetic equations, using the hypocoercivity framework
of Dolbeault, Mouhot & Schmeiser [4].

Developing Galerkin finite element methods which preserve/admit hypocoercivity structures is a largely unexplored area. This is despite the potential advantages of Galerkin methods, allowing for high order discretisations, over unstructured, possibly locally adapted meshes. We are only aware of [9], whereby a hypocoercivity-preserving family of spatially discrete non-conforming finite element methods for (1.1) is constructed and analysed. The key idea in the design of the methods in [9] is the construction of certain consistent stabilisations based on carefully selected numerical fluxes on the skeleton of the mesh, along with imposition of additional boundary conditions in the PDE problem. The modified inner product (·, ·)$_{hc}$, containing derivatives of order two or higher, results into the family of methods introduced in [9] to scale like a fourth order discrete operator in space. This may result to higher computational cost if used in conjunction with a fully discrete implicit time-stepping fashion and iterative linear solvers for each time-step.

This work is concerned with constructing a new, hypocoercivity-exploiting, stabilised finite element method, which address three key issues that arose in [9]: 1) the new method scales like a second order discrete operator, thereby reinstating the expected computational complexity of discretising a second order differential operator, 2) both spatially discrete and fully discrete methods are presented, with the time discretisation realised by an arbitrary order discontinuous Galerkin timestepping approach, and 3) no additional boundary conditions other than standard Cauchy/Dirichlet-Neumann ones are required for the analysis of the new method. On the other hand, however, the three aforementioned developments are possible at the cost of a discretization-parameter dependent, negative Grönwall exponent when compared to the one in [9]. In particular, the available spectral gap degenerates when increasing the dimension of the space-time approximation spaces (i.e., reducing meshsize/timestep and/or increasing polynomial degree). Nevertheless, for fixed discretisation parameters, hypocoercivity is proven and used, in turn, for the proof of a priori error bounds which are robust with respect to the final time $t_f$.

The remainder of this work is structured as follows. In Section 2, we discuss the special weak formulation of (1.1) which utilizes the concept of hypocoercivity in order to generate a contractive semigroup. Following this, we discuss two semi-discretisations. The first semi-discretisation is discussed in Section 3 comes in the form of a streamline-upwinded Petrov-Galerkin method which does not itself immediately achieve an optimal rate of convergence. This is followed in Section 4 by a discretisation which makes use of both the properties discussed in Section 2 and the streamline-upwinded Petrov-Galerkin method. Both of these semi-discretisations make use of continuous finite elements in space. Once the semi-discrete case has been treated, we move to Section 5 which considers a full discretisation of (1.1) by use of discontinuous Galerkin methods in time.

**Notation.** To simplify notation, we abbreviate the $L_2(\omega)$-inner product and $L_2(\omega)$-norm for a Lebesgue-measurable subset $\omega \subset \mathbb{R}^d$ as $(\cdot, \cdot)_\omega$ and $\|\cdot\|_\omega$, respectively.
tively. Moreover, when $\omega = \Omega$, we will further compress the notation to $(\cdot, \cdot) \equiv (\cdot, \cdot)_\Omega$ and $\|\cdot\| \equiv \|\cdot\|_\Omega$. The standard notation $H^k(\omega)$ for Hilbertian Sobolev spaces, $k \in \mathbb{R}$ will be used. In addition, given an interval $J \subset \mathbb{R}$ and a Banach space $V$, we use the standard notation Bochner spaces that are denoted by $L_p(J; V)$.

For the rest of this work, we restrict to $\Omega$ being a polygonal domain and introduce a family of triangulations of $\Omega$, say $\mathcal{T}$, consisting of mutually disjoint open triangular elements $T \in \mathcal{T}$, whose closures cover $\bar{\Omega}$ exactly. Also, let $h: \cup_{T \in \mathcal{T}} T \to \mathbb{R}_+$ be the local meshsize function defined elementwise by $h|_T := h_T := \text{diam}(T)$. For simplicity, we further assume that $\mathcal{T}$ is shape-regular, in the sense that the radius $\rho_T$ of the largest inscribed circle of each $T \in \mathcal{T}$ is bounded from below with respect to each element's diameter $h|_T$, uniformly as $\|h\|_{L_\infty(\Omega)} \to 0$ under mesh refinement. Also, we assume that $\mathcal{T}$ is locally quasi-uniform in the sense that the diameters of adjacent elements are uniformly bounded from above and below. Finally, let $\Gamma := \cup_{T \in \mathcal{T}} \partial T$ denote the mesh skeleton, and $\Gamma_{\text{int}} := \Gamma \setminus \partial \Omega$.

The family of finite element spaces subordinate to $\mathcal{T}$ we consider are defined by

$$V_h \equiv V_h^p := \{V \in H_0^1(\Omega) : V|_T \in \mathbb{P}_p(T), T \in \mathcal{T}\},$$

with $\mathbb{P}_p(\omega)$, $\omega \subset \mathbb{R}^d$, denoting the space of polynomials of total degree at most $p$ over $\omega$, $p = 1, 2, 3, \dots$. Notice that this choice of $V_h$ is only $H^1$-conforming despite the higher order derivatives which will later appear in the bilinear form of interest. Further, we define the broken Sobolev spaces $H^r(\Omega, \mathcal{T}) := \{v|_T \in H^r(T), T \in \mathcal{T}\}$, and the broken gradient $\nabla_{\mathcal{T}}$ to be the element-wise gradient operator with $\nabla_{\mathcal{T}} w|_T = \nabla w|_T$, $T \in \mathcal{T}$, for $w \in H^1(\Omega, \mathcal{T})$.

**2. A special weak formulation.** The problem of closing a degenerate parabolic PDE with suitable boundary conditions is well understood via the classical theory of linear second order equations with non-negative characteristic form [7, 14]. In particular, with $n(\cdot) := (n_1(\cdot), n_2(\cdot))^T$ denoting the unit outward normal vector at almost every point of $\partial\Omega$, which is assumed to be piecewise smooth and Lipschitz, we first define the elliptic portion of the boundary

$$\partial_0 \Omega := \{(x, y) \in \partial \Omega : n_1(x, y) \neq 0\}.$$

On the non-elliptic portion of the boundary $\partial\Omega \setminus \partial\Omega_0$, we define the inflow and outflow boundaries:

$$\partial_{-\Omega} := \{(x, y) \in \partial\Omega \setminus \partial_0\Omega : xn_2(x, y) < 0\}, \quad \partial_{+\Omega} := \{\xi \in \partial\Omega \setminus \partial_0\Omega : xn_2(x, y) \ge 0\}.$$

For notational brevity, any characteristic portions of the boundary $\{(x, y) \in \partial\Omega \setminus \partial_0\Omega : xn_2(x, y) = 0\}$ have been included into the outflow part, since their treatment is identical in what follows. Notice that $\partial\Omega = \partial_{+\Omega} \cup \partial_{-\Omega} \cup \partial_0\Omega$. Introducing the notation $Lu \equiv -u_{xx} + xu_y$, we consider the initial/boundary-value problem:

$$\begin{array}{ll}
u_t + Lu = u_t - u_{xx} + xu_y = f, & \text{in } (0, t_f] \times \Omega, \\
& \\
u = u_0, & \text{on } \{0\} \times \Omega, \\
& \\
u = 0, & \text{on } (0, t_f] \times \partial_{-\Omega}, \\
& \\
n_1 u_x = 0, & \text{on } (0, t_f] \times \partial_0\Omega,
\end{array} \tag{2.1}$$

for $t_f > 0$ and for $f \in L^2((0, t_f); H^1(\Omega))$, noting the stronger, non-standard regularity assumption for $f$. The well-posedness of the above problem is assured upon assuming
that $\partial_{-\Omega}$ has positive one-dimensional Hausdorff measure [14]. We refer to [14] for further details on the admissible boundary conditions for ultraparabolic PDEs and their physical meaning. We note that the developments below (with minor modifications) still apply upon replacing the non-flux condition by a Dirichlet one. We opt for the no-flux condition in the present work for its relevance in kinetic-type systems. For brevity, denote the time interval by $I := (0, t_f]$. We note that, in contrast to the previous work [9], in what follows, we do not prescribe any additional boundary conditions.

Let $A$ be a symmetric and non-negative definite matrix of the form

$$ (2.2) \qquad A := \begin{pmatrix} \alpha & \beta \\ \beta & \gamma \end{pmatrix} $$

where $\alpha, \beta, \gamma$ are non-negative parameters. Denote by $((\cdot, \cdot))_A$ the inner product

$$ ((w, v))_A := (w, v) + (\nabla w, A\nabla v), $$

for $w, v \in H^1(\Omega)$ with induced inner product $\|\cdot\|_A := ((\cdot, \cdot))_A^{\frac{1}{2}}$. This norm will play the role of $\|\cdot\|_{hc}$ from the Introduction.

Assuming sufficient regularity for the exact solution $u$, so that the following calculations are well defined, (1.1) implies

$$ (2.3) \qquad ((u_t, v))_A + ((Lu, v))_A = ((f, v))_A, $$

for all $v \in H^1_-(\Omega) := \{v \in H^1(\Omega) : v|_{\partial_{-\Omega}} = 0\}$, and set

$$ (2.4) \qquad a(u, v) := (u_x, v_x) + (xu_y, v) + (\nabla L u, A \nabla v). $$

**2.1. The hypocoercive structure of $a(\cdot, \cdot)$**. Setting $v = u$ into (2.4) gives

$$ a(u, u) = \|u_x\|^2 + (xu_y, u) + (-u_{xxx} + (xu_y)_x, \alpha u_x + \beta u_y) + (-u_{xxy} + xu_{yy}, \beta u_x + \gamma u_y), $$

which, upon integration by parts, gives

$$ (2.5) \qquad \begin{aligned} a(u, u) = \|u_x\|^2 &+ \frac{1}{2}\|\sqrt{xn_2}u\|_{\partial_+\Omega}^2 + \|\sqrt{A}\nabla u_x\|^2 + \beta\|u_y\|^2 + \alpha(u_y, u_x) \\ &+ (\frac{1}{2}xn_2\nabla u - n_1\nabla u_x, A\nabla u)_{\partial\Omega}. \end{aligned} $$

In [9] additional boundary conditions are considered, so that the last two terms on the right-hand side of (2.5) are non-negative and, thus, they can be omitted. After applying Young's and Cauchy-Schwarz' inequalities and electing $\alpha > 0$, $\beta = 4\alpha^2/9$, and $\gamma = \alpha^3/3$, for instance, we deduce

$$ (2.6) \qquad a(u, u) \ge \frac{1}{4}\|u_x\|^2 + \frac{\alpha}{3}\|u_{xx}\|^2 + \frac{\alpha^2}{9}\|u_y\|^2 + \frac{\alpha^3}{27}\|u_{xy}\|^2. $$

Thus, remarkably, using the non-standard test function $v - \nabla \cdot (A\nabla v)$ for well-chosen $A$, results in $L$ providing coercivity with respect to $u_y$ also: this is a manifestation of the concept of hypocoercivity, as $\partial_y = [x\partial_y, \partial_x]$; we refer to [18, Theorem 24] for details.
Remark 2.1. To emphasise the potential generality of the approach, one may consider a $(d+1)$-dimensional version of (1.1): for $\Omega \subset \mathbb{R}^{d+1}$ open, find $u: (0, t_f] \times \Omega \to \mathbb{R}$, such that

$$ (2.7) \qquad u_t + L_{d+1}u = u_t - u_{x_1x_1} + \sum_{i=1}^{d} x_i u_{x_{i+1}} = f, \quad \text{in } (0, t_f] \times \Omega, $$

for $t_f > 0$ and $f \in H^1(\Omega)$, along with suitable initial/boundary conditions. Interestingly, (2.7) is smoothing, despite possessing explicit diffusion in one spatial dimension only. Indeed, we have, respectively, $[\partial_{x_i}, \mathbf{b} \cdot \nabla] = \partial_{x_{i+1}}$ thereby, Hörmander's rank condition is satisfied [11], implying that (2.7) is, in fact, hypoelliptic! Moreover, since full rank is achieved via commutators involving the skew-symmetric 1st order part $\mathbf{b} \cdot \nabla u$, the PDE in (2.7) satisfies also the commutator hypotheses of [18, Theorem 24]. The developments discussed in this work can be transferred to this problem also with minor modifications only.

### 3. A streamline-upwinded Petrov-Galerkin formulation.
For equations with non-negative characteristic form, such as (2.1), a streamline-upwinded Petrov-Galerkin formulation (SUPG) was proposed [12] to achieve optimal rates of convergence in the presence of degenerate diffusion for continuous finite element spaces like $V_h^p$. In the present setting, the method in [12] reads: find $U \equiv U(t) \in V_h$, $t \in (0, t_f]$, such that

$$ (3.1) \qquad (U_t, V) + a_{h,\text{SUPG}}(U, V) = (f, V + \tau(V_t + xV_y)), $$

for all $V \in V_h$, with

$$ (3.2) \qquad a_{h,\text{SUPG}}(U,V) := (U_x, V_x) + (xU_y, V) + \sum_{T \in \mathcal{T}} (U_t + LU, \tau(V_t + xV_y))_T, $$

for an SUPG parameter $\tau: \Omega \to \mathbb{R}$ defined element-wise precisely below. We also set $U(0) := \pi u_0 \in V_h$, with $\pi: L_2(\Omega) \to V_h$ denoting the orthogonal $L_2$-projection onto the finite element space. The intuition behind testing the PDE against the non-standard term $\tau(V_t + xV_y)$ is to “symmetrize” the bilinear form with respect to the first order terms.

#### 3.1. Coercivity in SUPG-norm.
Setting $U = V = w \in V_h^p$ in (3.2), we have the coercivity relation:

$$ (3.3) \qquad \begin{aligned} a_{h,\text{SUPG}}(w,w) &= \|w_x\|^2 + \frac{1}{2} \|\sqrt{xn_2}w\|_{\partial_+\Omega}^2 + \|\sqrt{\tau}(w_t + xw_y)\|^2 \\ &\quad - \sum_{T \in \mathcal{T}} (w_{xx}, \tau(w_t + xw_y))_T \\ &\ge \|w_x\|^2 + \frac{1}{2} \|\sqrt{xn_2}w\|_{\partial_+\Omega}^2 + \|\sqrt{\tau}(w_t + xw_y)\|^2 \\ &\quad - \|C_{\text{INV}} p^2 h^{-1} \sqrt{\tau} w_x\| \|\sqrt{\tau}(w_t + xw_y)\| \ge \frac{1}{2} \|w\|_{\text{SUPG}}^2, \end{aligned} $$

for

$$ \|w\|_{\text{SUPG}} = \left( \|w_x\|^2 + \|\sqrt{xn_2}w\|_{\partial_+\Omega}^2 + \|\sqrt{\tau}(w_t + xw_y)\|^2 \right)^{\frac{1}{2}}, $$

whereby, for each $T \in \mathcal{T}$, we set

$$ (3.4) \qquad \tau|_T := C_{\text{INV}}^{-2} p^{-4} h^2, $$
with $C_{INV}$ the constant of the inverse estimate $\|\nabla v\|_T \le C_{INV} p^2 h_T^{-1} \|v\|_T$, valid for $v \in \mathbb{P}_p(T)$, depending only on the shape-regularity of the mesh.

For $f = 0$, starting from (3.1), we set $V = U$ and, employing (3.3) along with standard arguments, we arrive at

$$ \frac{d}{dt} \|U\|^2 + \|U\|_{\text{SUPG}}^2 \le 0, $$

which implies the expected stability:

$$ (3.5) \qquad \|U(t_f)\| \le \|U(0)\|, $$

for all $\tau \ge 0$, i.e., we have stability of the SUPG scheme with respect to the PDE initial data.

### 3.2. Error estimate. Set now

$$ \rho := u - \pi u, \quad \vartheta := \pi u - U, $$

so that $e = \rho + \vartheta$, whereby, following [3], $\pi: L_2(\Omega) \to V_h^p$ is the orthogonal $L_2$-projection operator onto the finite element space, giving $(\rho_t, \vartheta) = 0$.

Employing (3.1) with $V = \vartheta$ gives the error equation

$$ \begin{align} \tag{3.6} & \frac{1}{2} \frac{d}{dt} \|\vartheta\|^2 + a_{h,\text{SUPG}}(\vartheta, \vartheta) \\ & = ((\pi u)_t, \vartheta) + a_{h,\text{SUPG}}(\pi u, \vartheta) - (f, \vartheta_t + \tau(\vartheta_t + x\vartheta_y)), \\ & = -(\rho_t, \vartheta) - a_{h,\text{SUPG}}(\rho, \vartheta) = -a_{h,\text{SUPG}}(\rho, \vartheta). \end{align} $$

Using $(\rho, \vartheta_t) = 0$, we have $(\rho, x\vartheta_y) = (\rho, x\vartheta_y + \vartheta_t)$. Then, we have,

$$ \begin{align} a_{h,\text{SUPG}}(\rho, \vartheta) &= (\rho_x, \vartheta_x) - (\rho, x\vartheta_y) + \int_{\partial_+\Omega} xn_2\rho\vartheta \,ds \nonumber \\ & \quad + (\rho_t + x\rho_y, \tau(\vartheta_t + x\vartheta_y)) - \sum_{T \in \mathcal{T}} (\rho_{xx}, \tau(\vartheta_t + x\vartheta_y))_T \nonumber \\ \tag{3.7} & \le \left( \|\rho_x\| + \|\tau^{-\frac{1}{2}}\rho\| + \|\sqrt{xn_2}\rho\|_{\partial_+\Omega} + \|\sqrt{\tau}(\rho_t + x\rho_y)\| \right. \nonumber \\ & \qquad \left. + \left( \sum_{T \in \mathcal{T}} \|\sqrt{\tau}\rho_{xx}\|_T^2 \right)^{\frac{1}{2}} \right) \nonumber \\ & \qquad \times \left( \|\vartheta_x\| + \|\sqrt{xn_2}\vartheta\|_{\partial_+\Omega} + 2\|\sqrt{\tau}(\vartheta_t + x\vartheta_y)\| \right) \nonumber \\ & \le 4 \left( \|\rho\|_{\text{SUPG}} + \|\tau^{-\frac{1}{2}}\rho\| + \left( \sum_{T \in \mathcal{T}} \|\sqrt{\tau}\rho_{xx}\|_T^2 \right)^{\frac{1}{2}} \right) \|\vartheta\|_{\text{SUPG}}, \end{align} $$

where $\tau$ is defined in (3.4). Using (3.7) along with (3.3) into (3.6), we arrive at

$$ (3.8) \quad \frac{d}{dt} \|\vartheta\|^2 + \frac{1}{2} \|\vartheta\|_{\text{SUPG}}^2 \le 32 \left( \|\rho\|_{\text{SUPG}} + \|\tau^{-\frac{1}{2}}\rho\| + \left( \sum_{T \in \mathcal{T}} \|\sqrt{\tau} \rho_{xx}\|_T^2 \right)^{\frac{1}{2}} \right)^2. $$

We note that $\|\cdot\|_{\text{SUPG}}$ may only be a seminorm, in general, in $V_h^p$, for instance when $(0, y) \in \Omega$ for some $y \in \mathbb{R}$. In particular, it is not necessarily true that $\|\vartheta\| \le$
$C_{PF}$, for some constant $C_{PF}$, independent of $\vartheta$. Therefore, from (3.8), using standard arguments, we can only conclude the a priori error estimate

$$ (3.9) \qquad \begin{aligned} & \|e(t_f)\|^2 + \int_0^{t_f} \|e\|_{\text{SUPG}}^2 \, dt \le \|e(0)\|^2 + 2\|\rho(t_f)\|^2 \\ & \qquad + 128 \int_0^{t_f} \left( \|\rho\|_{\text{SUPG}} + \|\tau^{-\frac{1}{2}}\rho\| + \left( \sum_{T \in \mathcal{T}} \|\sqrt{\tau} \rho_{xx}\|_T^2 \right)^{\frac{1}{2}} \right)^2 \, dt. \end{aligned} $$

Assume, now, optimal best approximation properties for the projection $\pi$, viz., (3.10)

$$ \|h^{\kappa-1}\rho\|^2 + \|h^\kappa \nabla \rho\|^2 + \sum_{T \in \mathcal{T}} (\|h^{\kappa+\frac{1}{2}} \nabla \rho\|_{\partial T}^2 + \|h^{\kappa+1} \nabla^2 \rho\|_T^2) \le C_{app} \sum_{T \in \mathcal{T}} (|h^{\kappa+s} u|_{H^{s+1}(T)}^2), $$

for $u \in H^1(I; H^2(\Omega) \cap H^{s+1}(\Omega, \mathcal{T}))$, $s \ge 1$, and $\kappa \in \mathbb{R}$, with $\nabla^2$ denoting the (weak) Hessian of its argument. Assuming, now, sufficient space-time regularity for the exact solution $u$, the above best approximation estimates yield the error bound (3.11)

$$ \begin{aligned} \|e(t_f)\|^2 + \int_0^{t_f} \|e\|_{\text{SUPG}}^2 \, dt \le \|e(0)\|^2 + C \sum_{T \in \mathcal{T}} \left( \max_{t \in I} |h^s u|_{H^s(T)}^2 \right) \\ + C \sum_{T \in \mathcal{T}} \left( \int_0^{t_f} (|h^s u|_{H^{s+1}(T)}^2 + |h^{s+1} u_t|_{H^{s+1}(T)}^2) \, dt \right), \end{aligned} $$

for some $C > 0$, independent of $t_f$ and of $u$. Thus, if the exact solution were to satisfy $u \in H^1(I; H^2(\Omega) \cap H^{s+1}(\Omega, \mathcal{T}))$ for $1 \le s \le p$, (3.11) implies

$$ (3.12) \qquad \|e(t_f)\|^2 + \int_0^{t_f} \|e\|_{\text{SUPG}}^2 \, dt \le C_{tf} h_{\max}^{2s}, $$

with $h_{\max} = \max_{x \in \Omega} h$ for some $C > 0$ constant depending on $u$.

Remark 3.1 ($\tau \equiv 0$). We remark that, for $\tau \equiv 0$ in (3.1) and employing energy estimates, the resulting ‘plain’ Galerkin semi-discrete method may admit stronger dependence on $t_f$. Indeed, setting $\tau = 0$ and working as above, (3.3) degenerates to:

$$ (3.13) \qquad a_{h,\text{SUPG}}(\vartheta, \vartheta) = \|\vartheta_x\|^2 + \frac{1}{2} \|\sqrt{xn_2}\vartheta\|_{\partial_+ \Omega}^2, $$

while (3.7) gives

$$ \begin{aligned} a_{h,\text{SUPG}}(\rho, \vartheta) &= (\rho_x, \vartheta_x) - (\rho, x\vartheta_y) + \int_{\partial_+ \Omega} xn_2 \rho \vartheta \, ds \\ &\le (\|\rho_x\| + \|C_{INV} p^2 h^{-1} x \rho\| + \|\sqrt{xn_2} \rho\|_{\partial_+ \Omega}) (\|\vartheta_x\| + \|\sqrt{xn_2} \vartheta\|_{\partial_+ \Omega} + \|\vartheta\|), \end{aligned} $$

so that, (3.6) and standard arguments give

$$ \begin{aligned} & \frac{1}{2} \frac{d}{dt} \|\vartheta\|^2 + \frac{1}{4} \|\vartheta_x\|^2 + \frac{1}{4} \|\sqrt{xn_2}\vartheta\|_{\partial_+ \Omega}^2 \\ & \le 2 (\|\rho_x\| + \|C_{INV} p^2 h^{-1} x \rho\| + \|\sqrt{xn_2} \rho\|_{\partial_+ \Omega}) \|\vartheta\|, \end{aligned} $$

or

$$ \frac{d}{dt} \|\vartheta\|^2 \le \frac{16}{\epsilon} (\|\rho_x\| + \|C_{INV} p^2 h^{-1} x \rho\| + \|\sqrt{xn_2} \rho\|_{\partial_+ \Omega})^2 + \epsilon \|\vartheta\|^2, $$
for any $\epsilon > 0$. The triangle inequality, together with (3.10), then, implies

$$ (3.14) \qquad \|e(t_f)\|^2 \le C(\epsilon^{-1})e^{\epsilon t_f} h_{\max}^s, $$

for some constant $C > 0$ depending proportionally on $\epsilon^{-1}$ and on $u$, assuming the same regularity setting as in (3.12).

**4. A hypocoercivity-exploiting SUPG method.** Starting from (2.3), we consider the spatially discrete finite element method: find $U \equiv U(t) \in V_h$, $t \in (0, t_f]$, such that

$$ (4.1) \qquad ((U_t, V))_A + a_h(U, V) = ((f, V))_A + (f, \tau(V_t + xV_y)), $$

for all $V \in V_h$, for

$$ (4.2) \qquad a_h(U, V) := a_{h, \text{SUPG}}(U, V) + \sum_{T \in \mathcal{T}} (\nabla L U, A \nabla V)_T. $$

We set $U(0) := \hat{\pi}u_0 \in V_h$, with $\hat{\pi}: H_\perp^1(\Omega) \to V_h$ defined as the A-orthogonal projection operator onto the finite element space $V_h$.

**Remark 4.1.** The weak form (4.1) is a Petrov-Galerkin formulation arising formally by testing the PDE against $v + \tau(v_t + xv_y) - \nabla \cdot (A\nabla v)$, modulo a boundary term arising from integration by parts.

**4.1. Numerical hypocoercivity.** To study the numerical hypocoercivity, we first define

$$ \|w\| := \left( \frac{1}{2} \|w_x\|^2 + \|\sqrt{\gamma\delta}w_y\|^2 + \frac{1}{2} \|\sqrt{\tau}(w_t + xw_y)\|^2 + \|\sqrt{A\nabla\tau}w_x\|^2 + \|\sqrt{xn_2}w\|_{\partial_+\Omega}^2 + \sum_{T \in \mathcal{T}} \|\sqrt{xn_2A\nabla w}\|_{\partial_+T}^2 \right)^{\frac{1}{2}}, $$

with $A$ as in (2.2) and $\nabla\tau$ is the element-wise gradient operator, defined in the last section of the Introduction. We have the following coercivity result.

**LEMMA 4.2** ((Hypo)coercivity). For all $w \in V_h$, we have

$$ (4.3) \qquad a_h(w, w) \ge \frac{1}{4} \|w\|^2, $$

when $\alpha = (8\delta)^{-1}$, $\beta = (24\delta^2)^{-1}$, $\gamma = (64\delta^3)^{-1}$, whereby for each $T \in \mathcal{T}$,

$$ (4.4) \qquad \delta|_T := \delta_T = \max \left\{ \frac{C_{inv}^2 p^2 \|xn_2\|_{L_\infty(\partial_-T)}}{h_T}, \frac{2\|n_1\|_{L_\infty(\partial T)}^2 C_{inv}^4 p^4}{3h_T^2} \right\} \quad \text{and} \quad \tau|_T := \frac{h_T^2}{4C_{INV}^2 p^4}, $$

where $C_{inv} > 0$ the trace-inverse inequality constant depending only on the shape-regularity of the mesh, viz., $\|v\|_{\partial T} \le C_{inv} p h_T^{-\frac{1}{2}} \|v\|_T$, valid for $v \in \mathbb{P}_p(T)$, $T \in \mathcal{T}$.

*Proof.* Now, working as above, we have

$$ a_h(w, w) = \|w_x\|^2 + \frac{1}{2} \|\sqrt{xn_2}w\|_{\partial_+\Omega}^2 + \|\sqrt{\tau}(w_t + xw_y)\|^2 + \sum_{T \in \mathcal{T}} \left( (\nabla L w, A \nabla w)_T - (w_{xx}, \tau(w_t + xw_y))_T \right), $$
whilst the last term on the right-hand side gives

$$
\begin{align*}
\mathcal{I} &:= \|\sqrt{A}\nabla w_x\|_T^2 + \beta\|w_y\|_T^2 + \alpha(w_y, w_x)_T + \left(\frac{1}{2}xn_2\nabla w - n_1\nabla w_x, A\nabla w\right)_{\partial T} \\
&\quad - (w_{xx}, \tau(w_t + xw_y))_T \\
&\ge \|\sqrt{A}\nabla w_x\|_T^2 + \beta\|w_y\|_T^2 - \frac{\alpha^2}{2\epsilon_1}\|w_y\|_T^2 - \frac{\epsilon_1}{2}\|w_x\|_T^2 - \frac{1}{2}\|\sqrt{|xn_2|A}\nabla w\|_{\partial_-T}^2 \\
&\quad + \frac{1}{2}\|\sqrt{xn_2A}\nabla w\|_{\partial_+T}^2 - \|n_1\|_{L^\infty(\partial T)}\|\sqrt{A}\nabla w_x\|_{\partial T}\|\sqrt{A}\nabla w\|_{\partial T} \\
&\quad - \frac{1}{2}\|C_{INV}p^2h^{-1}\sqrt{\tau}w_x\|_T^2 - \frac{1}{2}\|\sqrt{\tau}(w_t + xw_y)\|_T^2 \\
&\ge (1 - \epsilon_2)\|\sqrt{A}\nabla w_x\|_T^2 + \left(\beta - \frac{\alpha^2}{2\epsilon_1}\right)\|w_y\|_T^2 - \left(\frac{\epsilon_1}{2} + \frac{1}{8}\right)\|w_x\|_T^2 \\
&\quad - \left(\frac{C_{inv}^2p^2\|xn_2\|_{L^\infty(\partial_-T)}}{2h_T} + \frac{\|n_1\|_{L^\infty(\partial T)}C_{inv}^4p^4}{4h_T^2\epsilon_2}\right)\|\sqrt{A}\nabla w\|_T^2 \\
&\quad - \frac{1}{2}\|\sqrt{\tau}(w_t + xw_y)\|_T^2 + \frac{1}{2}\|\sqrt{xn_2A}\nabla w\|_{\partial_+T}^2,
\end{align*}
$$

upon selecting $\tau := C_{INV}^{-2}p^{-4}h^2/4$. We set $\epsilon_2 = \frac{3}{4}$, and recall the definition of $\delta_T$, for which we note that $\delta_T \neq 0$ since $n_1^2 + n_2^2 = 1$ pointwise. Continuing, we have

$$
\begin{align*}
\mathcal{I} &\ge \frac{1}{4} \|\sqrt{A} \nabla w_x\|_T^2 + \left(\beta - \frac{\alpha^2}{2\epsilon_1}\right) \|w_y\|_T^2 - \left(\frac{\epsilon_1}{2} + \frac{1}{8}\right) \|w_x\|_T^2 - \delta_T \alpha \|w_x\|_T^2 \\
&\quad - \delta_T \gamma \|w_y\|_T^2 - 2\delta_T \beta (w_x, w_y)_T + \frac{1}{2} \|\sqrt{xn_2 A} \nabla w\|_{\partial_+ T}^2 - \frac{1}{2} \|\sqrt{\tau} (w_t + xw_y)\|_T^2 \\
&\ge \frac{1}{4} \|\sqrt{A} \nabla w_x\|_T^2 + \left(\beta - \frac{\alpha^2}{2\epsilon_1} - \delta_T \gamma - \delta_T \frac{\beta^2}{\alpha \epsilon_3}\right) \|w_y\|_T^2 \\
&\quad - \left(\frac{\epsilon_1}{2} + \frac{1}{8} + (1 + \epsilon_3) \delta_T \alpha\right) \|w_x\|_T^2 + \frac{1}{2} \|\sqrt{xn_2 A} \nabla w\|_{\partial_+ T}^2 - \frac{1}{2} \|\sqrt{\tau} (w_t + xw_y)\|_T^2.
\end{align*}
$$

Setting $\epsilon_1 = 1$, $\epsilon_3 = 1$, $\alpha = (8\delta_T)^{-1}$, $\beta = (24\delta_T^2)^{-1}$, and $\gamma = \frac{9}{8}\frac{\beta^2}{\alpha} = (64\delta_T^3)^{-1}$, we deduce

$$
\mathcal{I} \geq \frac{1}{4} \|\sqrt{A} \nabla w_x\|_T^2 + \frac{5\delta_T \gamma}{18} \|w_y\|_T^2 - \frac{7}{8} \|w_x\|_T^2 \\
+ \frac{1}{2} \|\sqrt{xn_2 A} \nabla w\|_{\partial_+ T}^2 - \frac{1}{2} \|\sqrt{\tau} (w_t + xw_y)\|_T^2.
$$

By combining the above estimates, we deduce

$$
\begin{align*}
a_h(w, w) \geq{}& \frac{1}{8} \|w_x\|^2 + \frac{5}{18} \|\sqrt{\gamma \delta} w_y\|^2 + \frac{1}{4} \|\sqrt{A} \nabla_{\tau} w_x\|^2 + \frac{1}{2} \|\sqrt{xn_2 w}\|_{\partial_+ \Omega}^2 \\
& + \frac{1}{2} \sum_{T \in \mathcal{T}} \|\sqrt{xn_2 A} \nabla w\|_{\partial_+ T}^2 + \frac{1}{2} \|\sqrt{\tau}(w_t + xw_y)\|^2,
\end{align*}
$$

from which the result follows.

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

Remark 4.3. The choice for $\alpha, \beta, \gamma$ is not unique in the above proof. The particular numbers have been chosen as a typical example of possible values.

The following result provides a mesh-dependent spectral gap.
LEMMA 4.4. Assume that $h_{\max} = \max_{x \in \Omega} h < 1$, and set $h_{\min} := \min_{x \in \Omega} h$. Then, there exists a constant $1 > c_{hc} > 0$, independent of $h$ and $p$, such that, for all $w \in H^1(I; H^1(\Omega) \cap H^2(\Omega, \mathcal{T}))$, we have

$$ \|w\|^2 \geq \kappa \|w\|_A^2, \quad \text{with} \quad \kappa \equiv \kappa(h, p) := c_{hc} h_{\min}^4 p^{-8}. $$

*Proof.* Let $B := \text{diag}(1, (32\delta^2)^{-1})$. Then we have, respectively,

$$ 2\|w\|^2 \geq \|\sqrt{B}\nabla w\|^2 = \|\sqrt{B - \nu A}\nabla w\|^2 + \nu\|\sqrt{A}\nabla w\|^2. $$

Upon defining $\bar{\delta} := \min_{T \in \mathcal{T}} C_{\delta}^{T} p^{4} h_{\min}^{-2}$, we immediately have $\delta \geq \bar{\delta}$. Selecting now $\nu = \bar{\delta}/2$, the eigenvalues $\lambda_i$, $i = 1, 2$ of $B - \nu A$ satisfy $\frac{1}{2} \leq \lambda_1 \leq 1$ and $\frac{h_{\min}^4}{384p^8} \leq \lambda_2 < \frac{1}{2}$, which implies

$$ \|\sqrt{B - \nu A}\nabla w\|^2 \geq \lambda_2 \|\nabla w\|^2 \geq C_{\text{PF}}^{-1} \lambda_2 \|w\|^2, $$

where $C_{\text{PF}} > 0$ is the Poincaré-Friedrichs/spectral gap inequality $\|w\|^2 \leq C_{\text{PF}} \|\nabla w\|^2$ for any $w \in H^1(\Omega) \cap H^2(\Omega, \mathcal{T})$, independent of $w$ and discretisation parameters. Thus the result follows for $c_{hc} = 2^{-1} \min\{(192C_{\text{PF}})^{-1}, \min_{T \in \mathcal{T}} C_{\delta}^{T}\}$. $\square$

On the other hand, for the hypocoercivity-exploiting SUPG method (4.1), for $f = 0$ and $V = U$ we have, respectively,

$$ \frac{1}{2} \frac{d}{dt} \|U\|_A^2 + \frac{\kappa}{4} \|U\|_A^2 \leq \frac{1}{2} \frac{d}{dt} \|U\|_A^2 + \frac{1}{4} \|U\|^2 \leq \frac{1}{2} \frac{d}{dt} \|U\|_A^2 + a_h(U, U) = 0, $$

implying

$$ (4.5) \qquad \|U(t_f)\| \leq e^{-\kappa t_f / 4} \|U(0)\|, $$

which improves upon (3.5). Therefore, the “hypocoercivity-inducing” stabilisation results into favourable dependence of the error with respect to the final time $t_f$ as $t_f \to \infty$. Even though the negative exponential decay may only be effective for $t_f > \kappa^{-1} \sim p^8 h_{\min}^{-4}$, it is nonetheless a manifestation of the concept of hypocoercivity, i.e., convergence to equilibrium despite the degenerate nature of the diffusion in Kolmogorov’s equation. We shall refer to the preservation of the hypocoercivity structure by a numerical method as numerical hypocoercivity.

## 4.2. Error estimate. We split the error as before,

$$ \rho := u - \hat{\pi}u, \quad \vartheta := \hat{\pi}u - U, $$

so that $e = \rho + \vartheta$. Employing (4.1) with $V = \vartheta$ and using (2.3), respectively, gives the error equation

$$ (4.6) \quad \begin{aligned} \frac{1}{2} \frac{d}{dt} \|\vartheta\|_A^2 + a_h(\vartheta, \vartheta) &= (((\hat{\pi}u)_t, \vartheta))_A + a_h(\hat{\pi}u, \vartheta) - ((f, \vartheta))_A - (f, \tau(V_t + xV_y)), \\ &= - ((\rho_t, \vartheta))_A - a_h(\rho, \vartheta) = -a_h(\rho, \vartheta), \end{aligned} $$
since $((\rho_t, \vartheta))_A = 0$, from the definition of the A-orthogonal projection. Similarly, use the relation $((\rho, \vartheta_t))_A = 0$, we have, respectively,

$$
\begin{align*}
a_h(\rho, \vartheta) &= (\rho_x, \vartheta_x) - (x\rho, \vartheta_y) + (xn_2\rho, \vartheta)_{\partial_+\Omega} + (\rho_t + x\rho_y, \tau(\vartheta_t + x\vartheta_y)) \\
&\quad + \sum_{T \in \mathcal{T}} \left( -(\rho_{xx}, \tau(\vartheta_t + x\vartheta_y))_T - (\nabla\rho_{xx}, A\nabla\vartheta)_T + (\nabla(x\rho_y), A\nabla\vartheta)_T \right) \\
&= (\rho_x, \vartheta_x) - (\rho, x\vartheta_y + \vartheta_t) + (xn_2\rho, \vartheta)_{\partial_+\Omega} + (\rho_t + x\rho_y, \tau(\vartheta_t + x\vartheta_y)) \\
&\quad + \sum_{T \in \mathcal{T}} \left( -(\rho_{xx}, \tau(\vartheta_t + x\vartheta_y))_T - (\nabla\rho_{xx}, A\nabla\vartheta)_T \right. \\
&\qquad \left. + (\nabla(x\rho_y), A\nabla\vartheta)_T - (\nabla\rho, A\nabla\vartheta_t)_T \right) \\
&= (\rho_x, \vartheta_x) - (\rho, x\vartheta_y + \vartheta_t) + (xn_2\rho, \vartheta)_{\partial_+\Omega} + (\rho_t + x\rho_y, \tau(\vartheta_t + x\vartheta_y)) \\
&\quad + \sum_{T \in \mathcal{T}} \left( -(\rho_{xx}, \tau(\vartheta_t + x\vartheta_y))_T + (\nabla\rho_x, A\nabla\vartheta_x)_T - (n_1\nabla\rho_x, A\nabla\vartheta)_{\partial T} \right. \\
&\qquad \left. + (\nabla(x\rho_y), A\nabla\vartheta)_T + (\nabla\rho, A\nabla(x\vartheta_y))_T - (\nabla\rho, A\nabla(\vartheta_t + x\vartheta_y))_T \right).
\end{align*}
$$

For the eighth and ninth terms on the right-hand side of the above identity, we have

$$
\begin{align*}
& (\nabla(x\rho_y), A\nabla\vartheta)_T + (\nabla\rho, A\nabla(x\vartheta_y))_T \\
& = (\rho_y(\nabla x), A\nabla\vartheta)_T + (\nabla\rho_y, A\nabla\vartheta y)_T + (\nabla\rho, A\nabla\vartheta y x)_T + (\nabla\rho, A(\nabla x)\vartheta_y)_T \\
& = (\alpha\rho_y, \vartheta_x)_T + 2(\beta\rho_y, \vartheta_y)_T + (\alpha\rho_x, \vartheta_y)_T + (xn_2\nabla\rho, A\nabla\vartheta)_{\partial T}.
\end{align*}
$$

Using the above relation, we have

$$
\begin{align*}
a_h(\rho, \vartheta) &= (\rho_x, \vartheta_x) - (\rho, x\vartheta_y + \vartheta_t) + (xn_2\rho, \vartheta)_{\partial_+\Omega} + (\rho_t + x\rho_y, \tau(\vartheta_t + x\vartheta_y)) \\
&\quad + \sum_{T \in \mathcal{T}} \left( -(\rho_{xx}, \tau(\vartheta_t + x\vartheta_y))_T + (\nabla\rho_x, A\nabla\vartheta_x)_T - (\nabla\rho, A\nabla(\vartheta_t + x\vartheta_y))_T \right. \\
&\qquad \left. + (\alpha\rho_y, \vartheta_x)_T + 2(\beta\rho_y, \vartheta_y)_T + (\alpha\rho_x, \vartheta_y)_T + (xn_2\nabla\rho, A\nabla\vartheta)_{\partial T} - (n_1\nabla\rho_x, A\nabla\vartheta)_{\partial T} \right).
\end{align*}
$$

For the last two terms on the right-hand side of the above identity, we have

$$
\begin{align*}
& (xn_2 \nabla \rho, A \nabla \vartheta)_{\partial T} - (n_1 \nabla \rho_x, A \nabla \vartheta)_{\partial T} = (xn_2 \rho_x, \alpha \vartheta_x + \beta \vartheta_y)_{\partial T} \\
& \qquad + (xn_2 \rho_y, \beta \vartheta_x + \gamma \vartheta_y)_{\partial T} - (n_1 \rho_{xx}, \alpha \vartheta_x + \beta \vartheta_y)_{\partial T} - (n_1 \rho_{xy}, \beta \vartheta_x + \gamma \vartheta_y)_{\partial T} \\
& \le C_{inv} p h_T^{-\frac{1}{2}} (\|\alpha xn_2 \rho_x\|_{\partial T} + \|\beta xn_2 \rho_y\|_{\partial T} + \|\alpha n_1 \rho_{xx}\|_{\partial T} + \|\beta n_1 \rho_{xy}\|_{\partial T}) \|\vartheta_x\|_T \\
& \qquad + C_{inv} p h_T^{-\frac{1}{2}} (\|\beta(\gamma\delta)^{-\frac{1}{2}} xn_2 \rho_x\|_{\partial T} + \|\sqrt{\gamma/\delta} xn_2 \rho_y\|_{\partial T} + \|\beta(\gamma\delta)^{\frac{1}{2}} n_1 \rho_{xx}\|_{\partial T} \\
& \qquad \qquad + \|\sqrt{\gamma/\delta} n_1 \rho_{xy}\|_{\partial T}) \|\sqrt{\gamma\delta}\vartheta_y\|_T
\end{align*}
$$

using the trace-inverse estimate in the last step. In addition, using the inverse inequality, definition of $\tau$ and $A$ in Lemma 4.3, we have

$$
\begin{align*}
-(\nabla\rho, A\nabla(\vartheta_t + x\vartheta_y))_T &\le \|A\|_{L_\infty(T)} \|\nabla\rho\|_T C_{INV} p^2 h_T^{-1} \|(\vartheta_t + x\vartheta_y)\|_T \\
&\le \|(C_{INV} p^2 h^{-1} \|A\|_{L_\infty(T)} \tau^{-\frac{1}{2}}) \nabla\rho\|_T \|\sqrt{\tau}(\vartheta_t + x\vartheta_y)\|_T \\
&\le \|C_{INV}(3C_\delta^T)^{-1} \nabla\rho\|_T \|\sqrt{\tau}(\vartheta_t + x\vartheta_y)\|_T,
\end{align*}
$$
where we have used the bound $||A||_{L_{\infty}(T)} \le (6\delta)^{-1}$, with $\delta$ defined Lemma 4.3. Performing estimation in standard fashion, and using the last bound, we have

$$
\begin{aligned}
a_h(\rho, \vartheta) \le & \|\rho_x\|\|\vartheta_x\| + \|\tau^{-\frac{1}{2}}\rho\|\|\sqrt{\tau}(x\vartheta_y + \vartheta_t)\| + \|\sqrt{xn_2}\rho\|_{\partial_+\Omega}\|\sqrt{xn_2}\vartheta\|_{\partial_+\Omega} \\
& + \|\sqrt{\tau}(\rho_t + x\rho_y)\|\|\sqrt{\tau}(\vartheta_t + x\vartheta_y)\| + \left(\sum_{T \in \mathcal{T}} \|\sqrt{\tau}\rho_{xx}\|_T^2\right)^{\frac{1}{2}} \|\sqrt{\tau}(\vartheta_t + x\vartheta_y)\| \\
& + \|\sqrt{A}\nabla_{\tau}\vartheta_x\|\|\sqrt{A}\nabla_{\tau}\rho_x\| + \left(\sum_{T \in \mathcal{T}} \|C_{INV}(3C_{\delta}^T)^{-1}\nabla\rho\|_T^2\right)^{\frac{1}{2}} \|\sqrt{\tau}(\vartheta_t + x\vartheta_y)\| \\
& + \|\alpha\rho_y\|\|\vartheta_x\| + 2\|\beta(\gamma\delta)^{-\frac{1}{2}}\rho_y\|\|\sqrt{\gamma\delta}\vartheta_y\| + \|\alpha(\gamma\delta)^{-\frac{1}{2}}\rho_x\|\|\sqrt{\gamma\delta}\vartheta_y\| \\
& + \left(\sum_{T \in \mathcal{T}} C_{inv}^2 p^2 h_T^{-1} (\|\alpha xn_2\rho_x\|_{\partial T}^2 + \|\beta xn_2\rho_y\|_{\partial T}^2 + \|\alpha n_1\rho_{xx}\|_{\partial T}^2 \right. \\
& \qquad \left. + \|\beta n_1\rho_{xy}\|_{\partial T}^2)\right)^{\frac{1}{2}} \|\vartheta_x\| + \left(\sum_{T \in \mathcal{T}} C_{inv}^2 p^2 h_T^{-1} (\|\beta(\gamma\delta)^{-\frac{1}{2}}xn_2\rho_x\|_{\partial T}^2 \right. \\
& \qquad \left. + \|\sqrt{\gamma/\delta}xn_2\rho_y\|_{\partial T}^2 + \|\beta(\gamma\delta)^{-\frac{1}{2}}n_1\rho_{xx}\|_{\partial T}^2 + \|\sqrt{\gamma/\delta}n_1\rho_{xy}\|_{\partial T}^2)\right)^{\frac{1}{2}} \|\sqrt{\gamma\delta}\vartheta_y\|.
\end{aligned}
$$

Rearranging terms, we have

$$
\begin{aligned}
a_h(\rho, \vartheta) \le & \|\sqrt{xn_2}\rho\|_{\partial_+\Omega}\|\sqrt{xn_2}\vartheta\|_{\partial_+\Omega} + \|\sqrt{A}\nabla_{\tau}\rho_x\|\|\sqrt{A}\nabla_{\tau}\vartheta_x\| + \left(\sum_{T \in \mathcal{T}} (\|\tau^{-\frac{1}{2}}\rho\|_T^2 \right. \\
& + \|\sqrt{\tau}(\rho_t + x\rho_y)\|_T^2 + \|\sqrt{\tau}\rho_{xx}\|_T^2 + \|C_{INV}(3C_{\delta}^T)^{-1}\nabla\rho\|_T^2) \left. \right)^{\frac{1}{2}} \|\sqrt{\tau}(\vartheta_t + x\vartheta_y)\| \\
& + \left(\sum_{T \in \mathcal{T}} (\|\rho_x\|_T^2 + \|\alpha\rho_y\|_T^2 + C_{inv}^2 p^2 h_T^{-1} (\|\alpha xn_2\rho_x\|_{\partial T}^2 + \|\beta xn_2\rho_y\|_{\partial T}^2 \right. \\
& + \|\alpha n_1\rho_{xx}\|_{\partial T}^2 + \|\beta n_1\rho_{xy}\|_{\partial T}^2) \left. \right)^{\frac{1}{2}} \|\vartheta_x\| \\
& + \left(\sum_{T \in \mathcal{T}} ((\gamma\delta)^{-1}(\|\alpha\rho_x\|_T^2 + 4\|\beta\rho_y\|_T^2) + C_{inv}^2 p^2 h_T^{-1} (\|\beta(\gamma\delta)^{-\frac{1}{2}}xn_2\rho_x\|_{\partial T}^2 \right. \\
& + \|\sqrt{\gamma/\delta}xn_2\rho_y\|_{\partial T}^2 + \|\beta(\gamma\delta)^{-\frac{1}{2}}n_1\rho_{xx}\|_{\partial T}^2 + \|\sqrt{\gamma/\delta}n_1\rho_{xy}\|_{\partial T}^2) \left. \right)^{\frac{1}{2}} \|\sqrt{\gamma\delta}\vartheta_y\|.
\end{aligned}
$$

Then, using the definition of A and $\tau$, and the approximation (3.10), we deduce

$$
\begin{aligned}
\frac{|a_h(\rho, \vartheta)|}{\|\vartheta\|} & \le C_{app} \left( \|\sqrt{hxn_2}\|_{L_\infty(\partial_+\Omega)} + \|\sqrt{A}h^{-1}\|_{L_\infty(\Omega)} + C_{INV} p^2 \right. \\
& \qquad + C_{INV}^{-1} p^{-2} (1 + \|xh\|_{L_\infty(\Omega)}) + C_{INV}(\max_{T \in \mathcal{T}}(C_{\delta}^T)) + 1 + \|\alpha\|_{L_\infty(\Omega)} \\
& \qquad + C_{inv} p (\|\alpha xn_2 h^{-1}\|_{L_\infty(\Gamma)} + \|\beta xn_2 h^{-1}\|_{L_\infty(\Gamma)} + \|\alpha h^{-2} n_1\|_{L_\infty(\Gamma)} \\
& \qquad + \|\beta h^{-2} n_1\|_{L_\infty(\Gamma)}) + \|\alpha(\gamma\delta)^{-\frac{1}{2}}\|_{L_\infty(\Omega)} + \|\beta(\gamma\delta)^{-\frac{1}{2}}\|_{L_\infty(\Omega)} \\
& \qquad \left. C_{inv} p (\|\beta(\gamma\delta)^{-\frac{1}{2}}h^{-1}xn_2\|_{L_\infty(\Gamma)} + \|\sqrt{\gamma/\delta}h^{-1}xn_2\|_{L_\infty(\Gamma)} + \|\beta(\gamma\delta)^{-\frac{1}{2}}h^{-2}n_1\|_{L_\infty(\Gamma)} \right. \\
& \qquad \left. + \|\sqrt{\gamma/\delta}h^{-2}n_1\|_{L_\infty(\Gamma)} \right) \left( \sum_{T \in \mathcal{T}} |h^s u|_{H^{s+1}(T)}^2 \right)^{\frac{1}{2}} + C_{INV}^{-1} p^{-2} \left( \sum_{T \in \mathcal{T}} |h^{s+1} u_t|_{H^{s+1}(T)}^2 \right)^{\frac{1}{2}} \\
& \le C (1 + \|xn_2\|_{L_\infty(\Gamma)} + \|xh\|_{L_\infty(\Omega)} \left( \sum_{T \in \mathcal{T}} |h^s u|_{H^{s+1}(T)}^2 + |h^{s+1} u_t|_{H^{s+1}(T)}^2 \right)^{\frac{1}{2}}
\end{aligned}
$$
for some $C > 0$, independent of $h$, of $t$, of $\|x\|_{L_\infty(\Omega)}$, and of $u$, but dependent, in principle from $p$ and the shape-regularity of the mesh. We note that the lack of dependence on $\|x\|_{L_\infty(\Omega)}$ is due to the presence of $\|xn_2\|_{L_\infty(\partial_-T)}$ in $\delta_T$, in that the ratio $\|x\|_{L_\infty(T)} / \|xn_2\|_{L_\infty(\partial_-T)}$, is uniformly bounded by a constant dependent only on the shape regularity of the mesh (perhaps times $h$).

Therefore, applying (4.3) and (4.7) into (4.6), gives

$$ (4.7) \qquad \frac{1}{2} \frac{d}{dt} \|\vartheta\|_A^2 + \frac{1}{8} \|\vartheta\|^2 \le C \sum_{T \in \mathcal{T}} (|h^s u|_{H^{s+1}(T)}^2 + |h^{s+1} u_t|_{H^{s+1}(T)}^2), $$

upon absorbing the dependence on $hx$ into the generic constant $C$.

Using the calculations above, we are now ready to state an *a priori* error bound.

**THEOREM 4.5.** *Assuming that the exact solution $u \in H^1(I; H^{2.5+\epsilon}(\Omega)) \cap H^1(I; H^{s+1}(\Omega, \mathcal{T}))$, for $1.5 < s \le p, p \ge 2$, then the error of (4.1) satisfies*

$$ (4.8) \qquad \|e(t_f)\|_A \le e^{-\frac{\kappa}{8}t_f} \|e(0)\|_A \\ \qquad \qquad + C \left( \int_0^{t_f} e^{\frac{\kappa}{8}(s-t_f)} \sum_{T \in \mathcal{T}} (|h^s u|_{H^{s+1}(T)}^2 + |h^{s+1} u_t|_{H^{s+1}(T)}^2) \, ds \right)^{\frac{1}{2}}, $$

for $C > 0$, independent of $t_f$, of $u$ and of $U$.

*Proof.* From the definition of the projection $\hat{\pi}$, we have $((\rho_t, \vartheta))_A = ((\vartheta_t, \rho))_A = ((\vartheta, \rho))_A = 0$. Using the latter, along with Lemma 4.4 and (4.7), we have, respectively,

$$ (4.9) \qquad \frac{1}{2} \frac{d}{dt} \|e(t)\|_A^2 + \frac{\kappa}{8} \|e(t)\|_A^2 \\ \qquad \le \frac{1}{2} \frac{d}{dt} \|\vartheta(t)\|_A^2 + \frac{\kappa}{4} \|\vartheta(t)\|_A^2 + ((\rho_t, \rho))_A + \frac{\kappa}{4} \|\rho(t)\|_A^2 \\ \qquad \le C \sum_{T \in \mathcal{T}} (|h^s u|_{H^{s+1}(T)}^2 + |h^{s+1} u_t|_{H^{s+1}(T)}^2 + |h^{s+3} u|_{H^{s+1}(T)}^2), $$

upon absorbing the dependence of $x$ in the generic constant $C$. Grönwall's Lemma already implies the result. □

Comparing (4.8) and (3.11), (or even (3.14),) we observe that the additional "hypocoercivity-inducing" stabilisation results into favourable dependence of the error with respect to the final time $t_f$ as $t_f \to \infty$.

**Remark 4.6** ($p = 1$). The error bound (4.8) only states the result for polynomial order $p \ge 2$. However, we point out that for $p = 1$, the proposed method will converge at an optimal order $O(h)$ under the regularity assumption $u \in H^{2.5+\epsilon}(\Omega)$.

**Remark 4.7.** It is possible to prove *hp*-version *a priori* error estimates upon considering $H^1$-conforming *hp*-version best approximation estimates instead of (3.10). However, available *hp*-version best approximation estimates in the literature for the $L_2$-projection error in the $H^2$-seminorm are significantly suboptimal with respect to the polynomial degree $p$. To avoid dwelling in such issues, we opted for not tracking the convergence rate with respect to $p$.

**5. A fully discrete hypocoercivity-exploiting method.** We shall now modify the spatially discrete method (4.1) above and further discretise the evolutionary derivative by a discontinuous Galerkin (dG) timestepping method.
We begin by subdividing $I = (0, t_f]$, into the family $\mathcal{I} := \{I_n : n = 1, \dots, N\}$, where $I_n := (t_{n-1}, t_n]$ for a strictly increasing sequence $0 = t_0 < t_1 < \dots < t_N = t_f$, and let $k : I \to \mathbb{R}_+$, with $k|_{I_n} := k_n = t_n - t_{n-1}$, the timestep. We define space-time finite element space

$$V_{h,k} \equiv V_{h,k}^{p,q} := \{V \in L_2(I; V_h^p) : V|_{I_n} \in \mathbb{P}_q(I_n), n = 1, \dots, N\}$$

for $q \in \mathbb{N}$. We denote by $[V]_n$ the time-jump across $t_n$, viz.,

$$[V]_n := V(t_n^+) - V(t_n^-) := \lim_{\epsilon \to 0^+} V(t_n + \epsilon) - \lim_{\epsilon \to 0^+} V(t_n - \epsilon)$$

for $n = 1, \dots, N$.

Starting from (4.1), we consider the following dG methods in time variable, hypocoercivity-exploiting SUPG method: find $U \in V_{h,k}$ such that

$$(5.1) \quad B(U, V) = ((U(t_0^-), V(t_0^+)))_A + \int_0^{t_f} \left( ((f, V))_A + (f, \tau(V_t + xV_y)) \right) dt$$

for all $V \in V_{h,k}$, whereby

$$(5.2) \quad B(U, V) := \sum_{n=1}^{N} \int_{I_n} ((U_t, V))_A + a_h(U, V)) \, dt + \sum_{n=2}^{N} (([U]_{n-1}, V(t_{n-1}^+)))_A \\ + ((U(t_0^+), V(t_0^+)))_A,$$

with $U(t_0^-) := \hat{\pi}u_0$, for $\hat{\pi}: H_0^1(\Omega) \to V_h$ A-orthogonal projection operator.

Of course, in practice (5.1) is typically solved on each time-step, viz.,

$$(5.3) \quad \int_{I_n} ((U_t, V))_A + a_h(U, V)) \, dt + ((U(t_{n-1}^+), V(t_{n-1}^+)))_A \\ = ((U(t_{n-1}^-), V(t_{n-1}^+)))_A + \int_{I_n} ((f, V))_A + (f, \tau(V_t + xV_y)) \, dt,$$

for all $V \in V_{h,k}$, $n = 1, \dots, N$. The SUPG parameter $\tau$ is defined as in (4.4) for all $T \in \mathcal{T}$ on $I_n$.

**5.1. Numerical hypocoercivity.** To study the numerical hypocoercivity structure of the space-time method (5.2), we set $V = U$ and $f = 0$ into (5.3) and, following standard arguments, (see, e.g., [17, Chapter 12],) we have

$$\frac{1}{2}\|U(t_n^-)\|_A^2 + \frac{1}{2}\|U(t_{n-1}^+)\|_A^2 - ((U(t_{n-1}^-), U(t_{n-1}^+)))_A + \int_{I_n} a_h(U, U) \, dt = 0,$$

or

$$(5.4) \quad \frac{1}{2}\|U(t_n^-)\|_A^2 + \frac{1}{2}\|[U]_{n-1}\|_A^2 + \int_{I_n} a_h(U, U) \, dt = \frac{1}{2}\|U(t_{n-1}^-)\|_A^2.$$

From Lemmata 4.2 and 4.4, we deduce

$$(5.5) \quad \|U(t_n^-)\|_A^2 + \|[U]_{n-1}\|_A^2 + \frac{\kappa}{2} \int_{I_n} \|U\|_A^2 \, dt \le \|U(t_{n-1}^-)\|_A^2.$$
Next, we manipulate the third term on the left-hand side of the last inequality to ensure a non-trivial spectral gap. Employing the hp-version trace inverse estimate [19], we have

$$ (5.6) \qquad k_n \|U(t_n^-)\|_A^2 \le (q+1)^2 \int_{I_n} \|U\|_A^2 dt, $$

Using the above estimate, we deduce

$$ \frac{\kappa k_n}{2(q+1)^2} \|U(t_n^-)\|_A^2 \le \frac{\kappa}{2} \int_{I_n} \|U\|_A^2 dt, $$

which upon substitution into (5.5) gives

$$ \left(1 + \frac{\kappa k_n}{2(q+1)^2}\right) \|U(t_n^-)\|_A^2 \le \|U(t_{n-1}^-)\|_A^2, $$

we deduce

$$ (5.7) \qquad \|U(t_f)\|_A^2 \le \prod_{n=1}^N \left(1 + \frac{\kappa k_n}{2(q+1)^2}\right)^{-1} \|U(t_0^-)\|_A^2, $$

which is, again, a manifestation of “numerical” hypocoercivity.

Remark 5.1. From (5.7), assuming constant timestep $k_n = k$ for all $n = 1, \dots, N$, we deduce

$$ (5.8) \qquad \|U(t_f)\|_A \le \left(1 + \frac{\kappa k}{2(q+1)^2}\right)^{-\frac{t_f}{2k}} \|U(t_0^-)\|_A \to e^{-\frac{\kappa}{4(q+1)^2} t_f} \|U(t_0^-)\|_A, $$

as $k \to 0$. Comparing with (4.5), the exponent has the same scale of the respective one in the semidiscrete case for $q = 0$.

5.2. Error estimate. We define the space-time norm

$$ (5.9) \qquad \|U\|_{\text{st}}^2 := \frac{1}{2} \left( \sum_{n=1}^{N-1} \|U\|_n^2_A + \|U(t_N^-)\|_A^2 + \|U(t_0^+)\|_A^2 \right) + \sum_{n=1}^N \int_{I_n} \frac{1}{4} \|U\|^2 dt. $$

Set $e := u - U$. The consistency of the method implies the error equation:
(5.10)

$$ \sum_{n=1}^{n} \int_{I_n} ((e_t, V))_A + a_h(e, V)) \, dt + \sum_{n=2}^{n} ((|e|_{n-1}, V(t_{n-1}^{+})))_A + ((e_0^{+}, V(t_0^{+})))_A = 0, $$

for all $V \in V_{h,k}$.

Definition 5.2. We define the space-time projection $\pi_{\text{st}}: H^1(I_n; H^1(\Omega)) \to V_{h,k}^{p,q}$ constructed for each $I_n$, $n = 1, \dots, N$, as the tensor product $\pi_{\text{st}} := \tilde{\pi}^t \circ \tilde{\pi} = \tilde{\pi} \circ \tilde{\pi}^t$, whereby $\tilde{\pi}: H^1(\Omega) \to V_h$ is the A-orthogonal spatial projection, defined by $((\tilde{\pi}v, V))_A = ((v, V))_A$ for all $V \in V_h$, and $\tilde{\pi}^t|_{I_n}: H^1(I_n) \to \mathbb{P}_q(I_n)$, such that $(\tilde{\pi}^t v)(t_n^-) = v(t_n^-)$, and

$$ \int_{I_n} \tilde{\pi}^t vw \, dt = \int_{I_n} vw \, dt, $$

for all $w \in \mathbb{P}_{q-1}(I_n)$.
We comment on the approximation properties and stability of $\tilde{\pi}^t$, confining ourselves to h-version estimates; see [16, 17] for details. For any function $v \in H^r(I_n)$ with $1 \le r \le q + 1$, we have

$$ (5.11) \qquad \|\partial_t^m (v - \tilde{\pi}^t v)\|_{I_n} \le C_{app} \|k^{r-m} \partial_t^r v\|_{I_n}, \quad m = 0,1 $$

and

$$ (5.12) \qquad |(v - \tilde{\pi}^t v)(t_{n-1}^+)| \le C_{app} \|k^{r-\frac{1}{2}} \partial_t^r v\|_{I_n}, $$

and the stability

$$ (5.13) \qquad \|\tilde{\pi}^t v\|_{I_n} \le \|v\|_{I_n} + C_{app} \|k \partial_t v\|_{I_n}, $$

where constant $C_{app}$ is independent of $k$.

Combining (5.11), (5.13) and the approximation results (3.10), we can deduce the following relation:

$$ (5.14) \qquad \int_{I_n} \|\nabla^\lambda (v - \pi_{\text{st}} v)\|_T^2 \, dt \le \tilde{C}_{app} \int_{I_n} \left( |k^r \partial_t^r v|_{H^\lambda(T)}^2 + |h^{s+1-\lambda} v|_{H^{s+1}(T)}^2 \right. \\ \qquad \qquad \qquad \qquad \left. + |h^{s+1-\lambda} k \partial_t v|_{H^{s+1}(T)}^2 \right) \, dt. $$

Similarly, using the trace inequality, we deduce

$$ (5.15) \qquad \int_{I_n} \|h^{\frac{1}{2}} \nabla^\lambda (v - \pi_{\text{st}} v)\|_{\partial T}^2 \, dt \le \tilde{C}_{app} \int_{I_n} \left( |k^r \partial_t^r v|_{H^\lambda(\partial T)}^2 + |h^{s+1-\lambda} v|_{H^{s+1}(T)}^2 \right. \\ \qquad \qquad \qquad \qquad \left. + |h^{s+1-\lambda} k \partial_t v|_{H^{s+1}(T)}^2 \right) \, dt, $$

for any functions $v$ satisfying $v \in H^1(I_n; H^{s+1}(\Omega, \mathcal{T})) \cap H^r(I_n; H^{2.5+\epsilon}(\Omega))$, with arbitrary small $\epsilon > 0$, $1 \le r \le q + 1$, $1.5 < s \le p$, $p \ge 2$ and $\lambda \le 2$.

We now split the error as $e = \eta + \xi$, where $\eta = u - \pi_{\text{st}}u$ and $\xi = \pi_{\text{st}}u - U$. From (5.10), along with completely analogous elementary manipulations to the derivation of (5.4), we arrive at

$$ (5.16) \qquad \frac{1}{2} \left( \|\xi(t_N^-)\|_A^2 + \|\xi(t_0^+)\|_A^2 + \sum_{n=2}^N \|[\xi]_{n-1}\|_A^2 \right) + \sum_{n=1}^N \int_{I_n} a_h(\xi, \xi) \, dt \\ = - \sum_{n=1}^N \int_{I_n} ((\eta_n, \xi))_A + a_h(\eta, \xi)) \, dt - \sum_{n=2}^N ((|\eta|_{n-1}, \xi(t_{n-1}^+)))_A - ((\eta(t_0^+), \xi(t_0^+)))_A, $$

upon setting $V = \xi$.

A key motivation for using the projection $\pi_{\text{st}}$ is the following orthogonality result.

LEMMA 5.3. For all $V \in V_h$ and $n = 1, \dots, N$, it holds

$$ (5.17) \qquad - \int_{I_n} ((\eta_t, V))_A \, dt = ((\eta(t_{n-1}^+), V(t_{n-1}^+)))_A. $$

*Proof.* By construction of $\pi_{\text{st}}$, we have

$$ ((\eta(t_n^-), V(t_n^-)))_A = ((u(t_n^-) - \hat{\pi}u(t_n^-), V(t_n^-)))_A = 0 \quad \text{and} \quad \int_{I_n} ((\eta, W))_A \, dt = 0, $$

whenever $W \in \mathbb{P}_{q-1}(I_n) \times V_h$. Hence, since $V_t \in \mathbb{P}_{q-1}(I_n) \times V_h$, integration by parts on the left-hand side of (5.17) already yields the result. □
Inserting, now, (5.17) into (5.16), we get

$$ \frac{1}{2} ( \|\xi(t_N^-)\|_A^2 + \|\xi(t_0^+)\|_A^2 + \sum_{n=2}^N \|\xi\|_{n-1} \|_A^2 ) + \sum_{n=1}^N \int_{I_n} a_h(\xi, \xi) \, dt = \sum_{n=1}^N \int_{I_n} a_h(\eta, \xi) \, dt, $$

or, using (4.3) and standard manipulations

$$ (5.18) \qquad \|\xi\|_{\text{st}}^2 \le \sum_{n=1}^{N} \int_{I_n} a_h(\eta, \xi) \, dt. $$

Working as in the proof of (4.7), we have

$$ \frac{|a_h(\eta, \xi)|}{\|\xi\|} \le \left( \|\sqrt{xn_2\eta}\|_{\partial_+\Omega}^2 + \|\sqrt{A}\nabla_T\eta_x\|^2 + \|\tau^{-\frac{1}{2}}\eta\|^2 + \|\sqrt{\tau}(\eta_t + x\eta_y)\|^2 \right. \\ + \sum_{T \in \mathcal{T}} \left( \|\sqrt{\tau}\eta_{xx}\|_T^2 + \|C_{\text{INV}}(3C_{\delta}^T)^{-1}\nabla\eta\|_T^2 + \|\eta_x\|_T^2 + \|\alpha\eta_y\|_T^2 \right. \\ + C_{\text{inv}}^2 p^2 h_T^{-1} \left( \|\alpha xn_2\eta_x\|_{\partial T}^2 + \|\beta xn_2\eta_y\|_{\partial T}^2 + \|\alpha n_1\eta_{xx}\|_{\partial T}^2 + \|\beta n_1\eta_{xy}\|_{\partial T}^2 \right) \\ + \|\alpha(\gamma\delta)^{-\frac{1}{2}}\eta_x\|_T^2 + 4\|\beta(\gamma\delta)^{-\frac{1}{2}}\eta_y\|_T^2 \\ + C_{\text{inv}}^2 p^2 h^{-1} \left( \|\beta(\gamma\delta)^{-\frac{1}{2}}xn_2\eta_x\|_{\partial T}^2 + \|\sqrt{\gamma/\delta}xn_2\eta_y\|_{\partial T}^2 \right. \\ \left. + \|\beta(\gamma\delta)^{-\frac{1}{2}}n_1\eta_{xx}\|_{\partial T}^2 + \|\sqrt{\gamma/\delta}n_1\eta_{xy}\|_{\partial T}^2 \right) \right)^{\frac{1}{2}}. $$

Then, using the definition of A, and the approximation (5.14), we deduce

$$ (5.19) \qquad \int_{I_n} |a_h(\eta, \xi)| \, dt \le \tilde{C}_{\text{app}} \left( \int_{I_n} \sum_{T \in \mathcal{T}} (T_T + S_T) \, dt \right)^{\frac{1}{2}} \left( \int_{I_n} \frac{1}{4} \|\xi\|^2 \, dt \right)^{\frac{1}{2}}, $$

where

$$ \begin{aligned} T_T := & \|\sqrt{xn_2}\|_{L_\infty(\partial_T \cap \partial_+ \Omega)}^2 \|k^r \partial_t^r u\|_{\partial_T \cap \partial_+ \Omega}^2 + \|\sqrt{A}\|_{L_\infty(T)}^2 \|k^r \partial_t^r (\nabla u_x)\|_T^2 \\ & + \tau_T^{-1} \|k^r \partial_t^r u\|_T^2 + \tau_T \|k^{r-1} \partial^r u\|_T^2 + \tau_T \|x\|_{L_\infty(T)}^2 \|k^r \partial_t^r u_y\|_T^2 + \tau_T \|k^r \partial_t^r u_{xx}\|_T^2 \\ & + \|C_{\text{INV}}(C_\delta^T)^{-1}\|_{L_\infty(T)}^2 \|k^r \partial_t^r (\nabla u)\|_T^2 + \|k^r \partial_t^r u_x\|_T^2 + \alpha^2 \|k^r \partial_t^r u_y\|_T^2 \\ & + C_{\text{inv}}^2 p^2 h_T^{-1} \left( \|\alpha xn_2\|_{L_\infty(\partial T)}^2 \|k^r \partial_t^r u_x\|_{\partial T}^2 + \|\beta xn_2\|_{L_\infty(\partial T)}^2 \|k^r \partial_t^r u_y\|_{\partial T}^2 \right. \\ & \qquad \left. + \|\alpha n_1\|_{L_\infty(\partial T)}^2 \|k^r \partial_t^r u_{xx}\|_{\partial T}^2 + \|\beta n_1\|_{L_\infty(\partial T)}^2 \|k^r \partial_t^r u_{xy}\|_{\partial T}^2 \right) \\ & + \|\alpha \gamma \delta^{-\frac{1}{2}}\|_{L_\infty(T)}^2 \|k^r \partial_t^r u_x\|_T^2 + \|\beta \gamma \delta^{-\frac{1}{2}}\|_{L_\infty(T)}^2 \|k^r \partial_t^r u_y\|_T^2 \\ & + C_{\text{inv}}^2 p^2 h_T^{-1} \left( \|\beta(\gamma \delta)^{-\frac{1}{2}} xn_2\|_{L_\infty(\partial T)}^2 \|k^r \partial_t^r u_x\|_{\partial T}^2 + \|\sqrt{\gamma/\delta} xn_2\|_{L_\infty(\partial T)}^2 \|k^r \partial_t^r u_y\|_{\partial T}^2 \right. \\ & \qquad \left. + \|\beta(\gamma \delta)^{-\frac{1}{2}} n_1\|_{L_\infty(\partial T)}^2 \|k^r \partial_t^r u_{xx}\|_{\partial T}^2 + \|\sqrt{\gamma/\delta} n_1\|_{L_\infty(\partial T)}^2 \|k^r \partial_t^r u_{xy}\|_{\partial T}^2 \right), \end{aligned} $$
and

$$
\begin{align*}
S_T := (\|\sqrt{xn_2h}\|_{L_\infty(\partial_T \cap \partial_+\Omega)}^2 + \|\sqrt{Ah}^{-1}\|_{L_\infty(T)}^2 + \tau_T^{-1}h_T^2 + \tau_T\|x\|_{L_\infty(T)}^2 \\
+ \tau_T h_T^{-2} + \|C_{INV}(C_\delta^T)^{-1}\|_{L_\infty(T)}^2 + 1 + \alpha^2 + C_{inv}^2 p^2 (\|\alpha xn_2 h^{-1}\|_{L_\infty(\partial T)}^2 \\
+ \|\beta xn_2 h^{-1}\|_{L_\infty(\partial T)}^2 + \|\alpha n_1 h^{-2}\|_{L_\infty(\partial T)}^2 + \|\beta n_1 h^{-2}\|_{L_\infty(\partial T)}^2) \\
+ \|\alpha \gamma \delta^{-\frac{1}{2}}\|_{L_\infty(T)}^2 + \|\beta \gamma \delta^{-\frac{1}{2}}\|_{L_\infty(T)}^2 + C_{inv}^2 p^2 (\|\beta(\gamma \delta)^{-\frac{1}{2}} xn_2 h_T^{-1}\|_{L_\infty(\partial T)}^2 \\
+ \|\sqrt{\gamma/\delta} xn_2 h_T^{-1}\|_{L_\infty(\partial T)}^2 + \|\beta(\gamma \delta)^{-\frac{1}{2}} n_1 h_T^{-2}\|_{L_\infty(\partial T)}^2 \\
+ \|\sqrt{\gamma/\delta} n_1 h_T^{-2}\|_{L_\infty(\partial T)}^2) (|h^s u|_{H^{s+1}(T)}^2 + |h^s k \partial_t u|_{H^{s+1}(T)}^2) + \tau_T |h^{s+1} \partial_t u|_{H^{s+1}(T)}^2.
\end{align*}
$$

Then, using the definition of A, τ, (5.19) and $\sum_{n=1}^{N} \int_{I_n} \frac{1}{4} \|\xi\|^2 dt \le \|\xi\|_{\text{st}}^2$, we deduce

$$
(5.20) \qquad \|\xi\|_{\text{st}}^2 \le C \sum_{n=1}^{N} \varepsilon_n,
$$

with $\mathcal{E}_n$ given by

$$
(5.21) \qquad \mathcal{E}_n := \int_{I_n} \sum_{T \in \mathcal{T}} \left( k_T^{2r} (\|\partial_t^r u\|_{\partial T \cap \partial_+ \Omega}^2 + (1+h_T^{-2} + (h_T/k_n)^2)\|\partial_t^r u\|_T^2 \right. \\
\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \left. + \|\partial_t^r \nabla u\|_T^2 + \|h^{\frac{1}{2}} \partial_t^r \nabla u\|_{\partial T}^2 + \|h \partial_t^r \nabla^2 u\|_T^2 + \|h^{\frac{3}{2}} \partial_t^r \nabla^2 u\|_{\partial T}^2 \right) \\
\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \left. + h_T^{2s} (|u|_{H^{s+1}(T)}^2 + |k \partial_t u|_{H^{s+1}(T)}^2) \right) dt,
$$

where $C > 0$, independent of $h_T$, $k_n$, and of $u$, but dependent on $p$, $q$, the shape-
regularity of the mesh, $\|x\|_{L_\infty(\Omega)}$ and $\|xn_2\|_{L_\infty(\partial_T)}$.

Remark 5.4. The fully discrete error bound (5.20) is optimal in terms of spatial mesh size $h$ in the space-time norm. Moreover, it is optimal in terms of time step $k$ under the condition that $h \approx \mathcal{O}(k^{\frac{1}{2}})$. For the sufficiently regular solution $u$, the optimal error bound with quasi-uniform spatial and temporal meshes are given below with $s = p$ and $r = q + 1$,

$$
\|\xi\|_{\text{st}}^2 \le C(h^{2p} + k^{2q+1}).
$$

Moreover, from the triangle inequality and best approximation, we deduce the follow-
ing (optimal) a priori error estimate

$$
(5.22) \qquad \|u - U\|_{\text{st}}^2 \le C(h^{2p} + k^{2q+1}).
$$

We are now ready to state an *a priori* error bound for the fully discrete scheme.

THEOREM 5.5. Assuming that the exact solution $u|_{I_n} \in H^1(I_n; H^{s+1}(\Omega, \mathcal{T})) \cap H^r(I_n; H^{2.5+\epsilon}(\Omega))$, with $1 \le r \le q+1$, $1.5 < s \le p$, $p \ge 2$, $1 \le n \le N$ and arbitrary small $\epsilon > 0$, then the error of (5.1) satisfies

$$
(5.23) \qquad \|e(t_f)\|_A^2 \le \prod_{n=1}^{N} (1+\mu_n)^{-1} \|u_0 - \hat{\pi} u_0\|_A^2 + C \sum_{n=1}^{N-1} \prod_{m=n+1}^{N} (1+\mu_m)^{-1} \mathcal{E}_n + C \mathcal{E}_N,
$$

for $\mu_n := \kappa k_n / (4(q+1)^2)$, $\kappa = \kappa(h,p) = c_{hc} h_{\min}^4 p^{-8}$ (as above), $\mathcal{E}_n$ defined in (5.21) and a constant $C > 0$, independent of $t_f$, of $u$ and of $U$.
Proof. First, using the dG-method (5.3) on each time interval $I_n$, with $n = 1, \dots, N$, and choosing $V = \xi$, we have the following error equation,

$$ (5.24) \quad \int_{I_n} ((\langle e_t, \xi \rangle)_A + a_h(e, \xi)) \, dt + (\langle |e|_{n-1}, \xi(t_{n-1}^+)) \rangle_A = 0. $$

Then, using the standard manipulation and (5.17) of projection $\pi_{st}$, we have

$$ \frac{1}{2} \|\xi(t_n^-)\|_A^2 + \frac{1}{2} \|\xi|_{n-1}\|_A^2 - \frac{1}{2} \|\xi(t_{n-1}^-)\|_A^2 + \int_{I_n} a_h(\xi, \xi) \, dt = - \int_{I_n} a_h(\eta, \xi) \, dt. $$

Next, applying Lemma 4.3 and Lemma 4.4, we have

$$ \frac{1}{2} \|\xi(t_n^-)\|_A^2 - \frac{1}{2} \|\xi(t_{n-1}^-)\|_A^2 + \frac{1}{4} \int_{I_n} \|\xi\|^2 \, dt \le \int_{I_n} \frac{a_h(\eta, \xi)}{\|\xi\|} \|\xi\| \, dt. $$

Then, standard estimates and Lemma 4.4 give

$$ \frac{1}{2} \|\xi(t_n^-)\|_A^2 - \frac{1}{2} \|\xi(t_{n-1}^-)\|_A^2 + \frac{\kappa}{8} \int_{I_n} \|\xi\|_A^2 \, dt \le C\varepsilon_n. $$

Applying the trace inverse inequality (5.6) and setting $\mu_n := \kappa k_n / (4(q+1)^2)$, we have

$$ (5.25) \quad (1 + \mu_n) \|\xi(t_n^-)\|_A^2 - \|\xi(t_{n-1}^-)\|_A^2 \le C\varepsilon_n. $$

Next, using $((\eta(t_n^-), \xi(t_n^-)))_A = ((\eta(t_{n-1}^-), \xi(t_{n-1}^-)))_A = 0$ along with (5.25), we have

$$ \begin{aligned} & (1 + \mu_n) \|e(t_n^-)\|_A^2 - \|e(t_{n-1}^-)\|_A^2 \\ &= (1 + \mu_n) \|\eta(t_n^-)\|_A^2 - \|\eta(t_{n-1}^-)\|_A^2 + (1 + \mu_n) \|\xi(t_n^-)\|_A^2 - \|\xi(t_{n-1}^-)\|_A^2 \\ &\le (1 + \mu_n) \|\eta(t_n^-)\|_A^2 + C\varepsilon_n. \end{aligned} $$

By the best approximation results (5.15), we have that $\|\eta(t_n^-)\|_A$ can be bounded by $\varepsilon_n$, thereby, concluding

$$ (1 + \mu_n) \|e(t_n^-)\|_A^2 \le \|e(t_{n-1}^-)\|_A^2 + C(1 + \mu_n)\varepsilon_n. $$

The last recurrence relation for $n = 1, \dots, N$, implies

$$ (5.26) \quad \|e(t_N^-)\|_A^2 \le \prod_{n=1}^{N} (1 + \mu_n)^{-1} \|e(t_0^-)\|_A^2 + C \sum_{n=1}^{N-1} \prod_{m=n+1}^{N} (1 + \mu_m)^{-1} \varepsilon_n + C\varepsilon_N. $$

Finally, using the identity $e(t_0^-) = u_0 - \hat{\pi}u_0$, (5.23) is proven. $\square$

**6. Numerical experiments.** We now present a basic numerical experiment which aims to verify the rates predicted by the theory. To begin with, we measure semi-discrete error which was estimated in (4.8). This is followed by measuring the error for an instationary problem which was estimated in (5.23). To complete, we demonstrate the exponential decay of the A-norm of a problem with no forcing which was considered in (5.8).

For the first and second examples of the experiments, we take $\Omega = (0, 1)^2$. This will be discretised by a sequence of successively refined uniform triangular meshes consisting of 32, 128, 512, 2048, and 8192 elements. The convergence rate of the error will be presented as a function of mesh size $h$ and of dof, where dof denotes the total number of space time basis functions.
<figure><img src="image_4.png" /><figcaption>FIG. 6.1. Experiment for the semi-discrete error.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 6.2. Experiment for the fully-discrete error.</figcaption></figure>

**6.1. Semi-discrete error.** For this we do not partition the time interval $(0, t_f)$ and use a single polynomial of order $q = 0$ on this. We measure the error $|||e|||_{st}$ defined in (5.22). By the stationary properties of the problem, this should not be dissimilar to the error $|||e(t_n^-)|||$. We set $t_f = 0$ and choose $f$ such that $u(t, x, y) = \sin^2(\pi x) \sin(\pi y)$. This results in the graph in Figure 6.1. The optimal convergence rate $O(h^p)$ for $p = 1, 2, 3, 4$ is observed.

**6.2. Fully-discrete error.** For this experiment, we consider $k = h^2$ as mentioned in Remark 5.4. The error is again measured by $|||e|||_{st}$. Again we set $t_f = 1$. Here, we choose $f$ such that

$$u(t, x, y) = \exp\left(-\left(x - 0.5\right)^2 - \left(y - 0.5\right)^2\right) \frac{\sin^2(\pi x) \sin(\pi(y - t - 0.5))}{2(t - 1)} + 1.$$

This satisfies a relatively inhomogeneous boundary condition on $y = 0$. The results of the experiment are provided in Figure 6.2. The optimal convergence rate $O(h^p)$ for $p = 1, 2, 3, 4$ is also observed upon selecting polynomial order $q = 0, 1, 2, 2$.

**6.3. Exponential convergence to equilibrium.** For this experiment, we take $\Omega = (-\frac{1}{2}, \frac{1}{2})^2$, $f \equiv 0$ and we set $t_f = 100$. We consider $k = h$. Rather than measuring error, we plot the value of $||U(t^-)||_A$ for the points $t$ which make up the right most end points of each interval, as well as the initial condition. We start with the initial condition $u_0(x, y) = \max(0, \frac{1}{4} - \max(|x|, |y|))$ which corresponds to a hat-type function in the middle of the domain. This results in the graphs in Figure 6.3.
<figure><img src="image_6.png" /><figcaption>FIG. 6.3. Convergence to equilibrium for different values of the discretization parameters.</figcaption></figure>

The exponential convergence observed in this particular example is faster than the one postulated in (5.8), as it does not appear to be adversely dependent on *h* and/or *p*. It remains an open question whether the degeneration of exponential decay for small *h* and/or large *p* can be witnessed in carefully constructed numerical examples and it is a direction of future research.

REFERENCES

[1] M. BESSEMOULIN-CHATARD AND F. FILBET, *A finite volume scheme for nonlinear degenerate parabolic equations*, SIAM J. Sci. Comput., 34 (2012), pp. B559–B583, https://doi.org/10.1137/110853807.

[2] M. BESSEMOULIN-CHATARD, M. HERDA, AND T. REY, *Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations*, Math. Comp., 89 (2020), pp. 1093–1133, https://doi.org/10.1090/mcom/3490.

[3] K. CHRYSAFINOS AND L. S. HOU, *Error estimates for semidiscrete finite element approximations of linear and semilinear parabolic equations under minimal regularity assumptions*, SIAM J. Numer. Anal., 40 (2002), pp. 282–306, https://doi.org/10.1137/S0036142900377991.

[4] J. DOLBEAULT, C. MOUHOT, AND C. SCHMEISER, *Hypocoercivity for linear kinetic equations conserving mass*, Trans. Amer. Math. Soc., 367 (2015), pp. 3807–3828, https://doi.org/10.1090/S0002-9947-2015-06012-7.

[5] G. DUJARDIN, F. HÉRAU, AND P. LAFITTE, *Coercivity, hypocoercivity, exponential time decay and simulations for discrete Fokker-Planck equations*, Numer. Math., 144 (2020), pp. 615–697, https://doi.org/10.1007/s00211-019-01094-y.

[6] J.-P. ECKMANN AND M. HAIRER, *Spectral properties of hypoelliptic operators*, Comm. Math. Phys., 235 (2003), pp. 233–253, https://doi.org/10.1007/s00220-003-0805-9.

[7] G. FICHERA, *Sulle equazioni differenziali lineari ellittico-paraboliche del secondo ordine*, Atti Accad. Naz. Lincei Mem. Cl. Sci. Fis. Mat. Natur. Sez. Ia (8), 5 (1956), pp. 1–30.

[8] E. L. FOSTER, J. LOHÉAC, AND M.-B. TRAN, *A structure preserving scheme for the Kolmogorov-Fokker-Planck equation*, J. Comput. Phys., 330 (2017), pp. 319–339, https://doi.org/10.1016/j.jcp.2016.11.009.

[9] E. H. GEORGOULIS, *Hypocoercivity-compatible finite element methods for the long-time computation of Kolmogorov's equation*, SIAM J. Numer. Anal., 59 (2021), pp. 173–194, https://doi.org/10.1137/19M1296914.

[10] F. HÉRAU AND F. NIER, *Isotropic hypoellipticity and trend to equilibrium for the Fokker-Planck equation with a high-degree potential*, Arch. Ration. Mech. Anal., 171 (2004), pp. 151–218, https://doi.org/10.1007/s00205-003-0276-3.

[11] L. HÖRMANDER, *Hypoelliptic second order differential equations*, Acta Math., 119 (1967), pp. 147–171, https://doi.org/10.1007/BF02392081.

[12] P. HOUSTON AND E. SÜLI, *Stabilised hp-finite element approximation of partial differential equations with nonnegative characteristic form*, vol. 66, 2001, pp. 99–119, https://doi.org/10.1007/s006070170030. Archives for scientific computing. Numerical methods for transport-dominated and related problems (Magdeburg, 1999).
[13] C. MOUHOT AND L. NEUMANN, *Quantitative perturbative study of convergence to equilibrium for collisional kinetic models in the torus*, Nonlinearity, 19 (2006), pp. 969–998, https://doi.org/10.1088/0951-7715/19/4/011.

[14] O. A. OLEĪNIK AND E. V. RADKEVIČ, *Second order equations with nonnegative characteristic form*, Plenum Press, New York-London, 1973. Translated from the Russian by Paul C. Fife.

[15] A. PORRETTA AND E. ZUAZUA, *Numerical hypocoercivity for the Kolmogorov equation*, Math. Comp., 86 (2017), pp. 97–119, https://doi.org/10.1090/mcom/3157.

[16] D. SCHÖTZAU AND C. SCHWAB, *Time discretization of parabolic problems by the hp-version of the discontinuous Galerkin finite element method*, SIAM J. Numer. Anal., 38 (2000), pp. 837–875, https://doi.org/10.1137/S0036142999352394.

[17] V. THOMÉE, *Galerkin finite element methods for parabolic problems*, vol. 25 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, second ed., 2006.

[18] C. VILLANI, *Hypocoercivity*, Mem. Amer. Math. Soc., 202 (2009), pp. iv+141, https://doi.org/10.1090/S0065-9266-09-00567-5.

[19] T. WARBURTON AND J. S. HESTHAVEN, *On the constants in hp-finite element trace inverse inequalities*, Comput. Methods Appl. Mech. Engrg., 192 (2003), pp. 2765–2773, https://doi.org/10.1016/S0045-7825(03)00294-9.