{'fraction_non_alphanumeric': 0.05200856679786988, 'fraction_numerical': 0.014413058578374624, 'mean_word_length': 4.319938414164742, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hypocoercivity-compatible finite element methods for the long-time computation of Kolmogorov's equation

This work is concerned with the development of a family of Galerkin finite element methods for the classical Kolmogorov's equation. Kolmogorov's equation serves as a sufficiently rich, for our purposes, model problem for kinetic-type equations and is characterised by diffusion in one of the two (or three) spatial directions only. Nonetheless, its solution admits typically decay properties to some long time equilibrium, depending on closure by suitable boundary/decay-at-infinity conditions. A key attribute of the proposed family of methods is that they also admit similar decay properties at the (semi)discrete level for very general families of triangulations. The method construction uses ideas by the general theory of hypocoercivity developed by Villani [23], along with judicious choice of numerical flux functions. These developments turn out to be sufficient to imply that the proposed finite element methods admit a priori error bounds with constants independent of the final time, despite Kolmogorov equation's degenerate diffusion nature. Thus, the new methods provably allow for robust error analysis for final times tending to infinity. The extension to three spatial dimensions is also briefly discussed.


Introduction
Degenerate parabolic problems in kinetic modelling are often characterised by the explicit presence of diffusion/dissipation in some of the spatial directions only, yet may still admit decay properties to some long time equilibrium. The development of Galerkin finite element methods for such problems satisfying provably also such decay properties is, remarkably, an unexplored area. This is despite the potential attraction that general, possibly highly nonuniform or adaptive, grid approximations can bring to kinetic simulations, especially those concerning expansive long-time simulation phenomena. Aiming to contribute in this direction, we shall develop and analyse a new class of Galerkin finite element methods for a simple kinetic equation, equipped with suitable closures by known boundary/periodicity conditions from the literature. To this end, we shall be predominantly concerned with the classical Kolmogorov's partial differential equation (PDE) for suitably smooth forcing f ∈ L 2 (R + ; H 1 (Ω)), and Cauchy initial data u(0, ·) = u 0 ∈ L 2 (Ω). The notational conventions used above are somewhat non-standard in kinetic theory, whereby L may be written as Lf := f t − f vv + vf x with v denoting the particle velocity variable, x the displacement/position and f the respective probability density function. In Section 6, we shall also briefly discuss the extension of the developments in this work to the three-dimensional version of (1), namely Kolmogorov showed already in 1934 that (1) admits a smooth fundamental solution outside the pole and, hence, it is hypoelliptic, i.e., it admits a smooth solution in (0, ∞) × R 2 for any smooth initial and forcing data u 0 , f , respectively [16]. This may be a surprising assertion at first sight, since there is no explicit dissipation built into the PDE in the y-variable. The advection field (0, x) T , however, is "appropriately non-constant" so that it suffices to propagate the built-in dissipation in the x-variable onto the entire spatial domain. In 1967, Hörmander, in the celebrated work [14], gave a sufficient condition for hypoellipticity for a wide class of 2nd order operators with non-negative characteristic form, which are nowadays often referred to as Hörmander sum-of-squares operators; interestingly, a motivating example for the developments in [14] has, indeed, been Kolmogorov's equation (1). Hence, although (1) may appear to be a rather special equation at first sight, its significance as a model problem is paramount, for it encompasses a number of pertinent structural properties of large classes of kinetic models; we refer, e.g., to [23] for very instructive expositions. The proof of trend to equilibrium for the directly related inhomogeneous Fokker-Planck equation has been given by Hérau & Nier in [13] upon realising that certain entropies admitting mixed derivatives (u xy or, respectively, f vx ) give rise to full gradients in certain weighted Sobolev spaces, from which Poincaré-Friedrichs inequalities (also known as spectral gap properties in this context) are sufficient to prove decay to an equilibrium distribution. Related ideas have also been used by Eckmann & Hairer in [7] to study the spectral properties of certain hypoelliptic operators involving degenerate diffusions and by Mouhot & Neumann [18] in the study of kinetic models with integral-type collision operators, among others. The idea of using entropies involving mixed derivatives was elevated to a general framework in proving decay to equilibrium for kinetic equations by Villani [23] via the introduction of the concept of hypocoercivity. Roughly speaking, hypocoercivity is the property of certain degenerate elliptic or parabolic differential operators to yield dissipation of the solution also in the directions where no diffusion is explicitly present. Astonishingly, sufficient conditions for hypocoercivity are given by the rank spanned by Hörmander vector fields and their commutators with respect to the vector field related to 1st order term [23]. Thus, sufficient conditions for hypocoercivity are closely related, but not identical, to the classical Hörmander's rank condition for hypoellipticity.
To fix ideas, denoting by (·, ·) ≡ (·, ·) L 2 (Ω) the inner product of L 2 (Ω), and respective norm · 2 := (·, ·), and assuming decay of u as |(x, y)| → ∞, the classical energy-type analysis for (1) equation yields for almost all t ∈ (0, t f ], for some final time t f ∈ R + , upon observing that the skew-symmetric term (xu y , u) vanishes. A combination of the Cauchy-Schwarz' and Young's inequalities on the last term on the right-hand side (3), therefore, yields for any δ > 0. Therefore, Grönwall's Lemma implies Since the right-hand side of (5) is always greater than or equal to u 0 2 for every f ∈ L 2 (0, ∞; L 2 (Ω)), we conclude that the basic energy estimate may be severely pessimistic as t f → ∞ for any fixed δ > 0: for any finite final time t f > 0, the above implies that the righthand side of (5) will be positive away from zero for any non-trivial initial condition. Notice that when f = 0, we can take δ → 0, yielding the stability estimate u(t f ) ≤ u 0 . This is at odds with the fast decay to 0 observed for this problem by the solution of Kolmogorov's equation for the respective homogeneous problem, at least when Ω = R 2 . Evidently the key culprit in the above analysis is the absence of control in a time-integrated norm on the left-hand side of (5). The situation does not appear to improve if different versions of Grönwall-type inequalities are applied instead.
On the other hand, for the non-degenerate parabolic equatioñ with the same forcing and initial conditions, we can easily deduce the energy identity Assuming now the validity of a Poincaré-Friedrichs inequality of the form u 2 ≤ C P F ∇u 2 , and working as before, we can arrive at the estimate from which, Grönwall's Lemma implies whose right-hand side clearly decays exponentially as t f → ∞ when f = 0, for instance. This is in sharp contrast with (5). The above is obviously not a satisfactory state of affairs not only at the PDE level, but also under the prism of the numerical analysis for these problems. Indeed, almost all provably convergent, general numerical methods which can be potentially applied to this class of PDEs use such standard energy-type arguments for their error analysis, resulting to possibly severely pessimistic a priori and/or a posteriori error bounds for large final times t f . Indeed, this is the state of affairs discussed in various works on Galerkin-type numerical methods for (degenerate) second order evolution problems of diffusion type; we non-exhaustively refer to [1,15,8,3] as representatives of different approaches.
To the best of our knowledge, there exist only few, yet quite inspiring, exceptions which employ hypocoercivity ideas in the context of numerical methods. In particular, the work of Porretta & Zuazua [20] proves the hypocoercivity of classical finite difference operators discretising Kolmogorov's equation for f = 0 on uniform spatial grids over [0, ∞) ×R 2 . The recent manuscript by Dujardin, Herau, & Lafitte [6] takes the approach of [20] further by establishing discrete versions of weighted Poincaré inequalities for difference operators, thereby showing decay to equilibrium for finite difference approximations of the inhomogeneous two-dimensional Fokker-Planck equation, which is, of course, closely related to Kolmogorov's equation. Both these approaches apply directly the abstract semigroup result of Villani [23,Theorem 18] to the finite difference operators used by proving the necessary commutator properties required for its validity. With regard to practical approaches, Foster, Lohéac & Tran [10] discuss the development 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 [2] present design principles for the construction of equilibrium-preserving finite volume methods for nonlinear degenerate parabolic problems.
For practical computations, one is typically forced to confine the spatial computational domain into an open and bounded one, say Ω ⊂ R d d = 2, 3. If the objective is to compute equilibrium states for the whole of R d , the truncation of the computational domain and the imposition of boundary conditions should be performed carefully, for it may alter the rate of decay to equilibrium compared to the non-confined problem [10]. If one is interested instead in computations on bounded domains as such, then the modelling often dictates the imposition of additional boundary conditions/constraints, such as periodicity or far field decay. As we shall see below, the imposition of boundary conditions in a fashion allowing concepts of hypocoercivity to be ported to bounded spatial domains is a challenge both in the PDE and the numerical analysis contexts. 2 This work is concerned with the development of a Galerkin-type finite element method which is compatible with hypocoercivity concepts and, therefore, allows for the proof of error bounds whose error constants do not grow as the final time t f → ∞, such as the ones in (5) do. Thus, as we shall show below, there will be no dependence of the error constant on the final time t f . This, in turn, supports the potential advantages of the numerical methods presented below for long time computations. The construction of the finite element method is based on a variational interpretation of the modified entropies involving mixed derivatives of Hérau & Nier [13] and Villani [23], which are able to reveal the implicit diffusion properties in the y-variable also. The derived modified variational problem is based on the PDE and higher order laws stemming from the PDE and involving up to 4th order partial derivatives, when viewed formally in strong form. The finite element method is constructed via a Galerkin projection onto a suitable C 0 -finite element space subordinate to a non-uniform, in general, triangulation. The proposed method is, nevertheless, non-conforming due to the 4th order symbol of the variational problem. To recover, therefore, the consistency of the finite element method without compromising on the hypocoercivity properties, the variational form has to be enriched by carefully constructed consistent numerical fluxes over the skeleton of the triangulation. The proposed finite element method is proven to decay to equilibrium when f = 0, when we equip the PDE with additional natural boundary conditions; this is, to the best of our knowledge, the first manifestation of hypocoercivity for a Galerkin scheme in the literature. Moreover, we prove a priori error bounds in the natural norm determined by the hypocoercivity result.
We view the developments presented below as proof of concept that Galerkin finite element methods over general non-uniform grids can be constructed in a compatible fashion to retain hypocoercivity properties at discrete level. This, in turn, leads to a priori error bounds which are robust with respect to their dependence on the final time t f . Kolmogorov's equation has the role of a sufficiently rich model problem to highlight a new finite element methodology for hypocoercive PDEs on bounded computational domains with view to retaining the longtime behaviour of the exact initial/boundary value problems. Apart from the mathematical challenge, the development of numerical methods for kinetic equations on general, possibly nonquasi-uniform meshes has important practical ramifications. Indeed, modern applications of kinetic equations involve the movement of large particles for which Lagrangian schemes appear to be preferential by practitioners which greatly benefits from the ability to discretise over unstructured, possibly moving, meshes.
The remainder of this work is structured as follows. In Section 2, we introduce a variational formulation motivated by design concepts of hypocoercivity [23], which is shown to be coercive with respect to an H 1 -equivalent norm in Section 3. A discrete bilinear form using continuous, yet non-conforming, elements with appropriately constructed numerical fluxes ensuring the coercivity in the discrete setting also is proposed in Section 4 and is shown to admit completely analogous decay to equilibrium properties as the respective continuous problem in the nonexternal forcing scenario. The error analysis of the proposed finite element method is given in Section 5 along with a discussion on its properties. To highlight the potential of the proposed framework for the class of kinetic equations with degenerate diffusion, we briefly discuss the extension of the developments in the first five sections to the case of the three dimensional (in space)version of Kolmogorov's equation; this is the content of Section 6.

A special weak formulation
The problem of closing a degenerate parabolic PDE with suitable boundary conditions is well understood via the, now classical, theory of linear second order equations with non-negative characteristic form [9,19]. In particular, with n(·) := (n 1 (·), n 2 (·)) T denoting the unit outward normal vector at each point of ∂Ω, we define the elliptic boundary along with the inflow and outflow parts of the non-elliptic boundary: Introducing, now, the notation Lu ≡ −u xx + xu y , we consider the initial/boundary-value problem: with ∂ −,0 Ω := ∂ − Ω ∪ ∂ 0 Ω, for t f > 0 and for f ∈ H 1 (Ω); note the non-standard regularity assumption for f . The well-posedness of the above problem is assured upon assuming that ∂ −,0 Ω has positive one-dimensional Hausdorff measure [19]. It is also possible to impose Neumanntype boundary conditions on parts of ∂ 0 Ω; this is not done here in the interest of simplicity of the presentation only. For the three-dimensional counterpart of the above model problem we refer to Section 6.
Let A := α β β γ a symmetric and non-negative definite matrix with α, β, γ non-negative parameters, and set the fact that A is non-negative definite implies that 0 ≤ δ ≤ 1. Assuming sufficient regularity for the exact solution u for the moment, so that the following calculations are well defined, (1) implies ∇u t + ∇Lu = ∇f, which, tested against A∇v for any v ∈ H 1 (Ω), results to (∇u t , A∇v) + (∇Lu, A∇v) = (∇f, A∇v); here and below, we denote (·, ·) ≡ (·, ·) L 2 (Ω) and · ≡ · L 2 (Ω) for brevity. Resorting to (1) once more, after an integration by parts, we deduce the variational form . The latter will be the basis of constructing a finite element method with the sought-after properties. To highlight the idea behind this construction in the context of hypocoercivity, we shall now discuss the properties of the spatial part of the operator on the left-hand side of (12), viz., 3 The hypocoercivity of a(·, ·) We begin by somewhat modifying the approach presented in Villani [23] arriving, nevertheless, to an effectively similar result. For accessibility, we shall first revert to (1), rather than (10), i.e., we have Ω = R 2 , assuming that |u(x, y)| → 0 as |(x, y)| → ∞; the case Ω ⊂ R 2 will be discussed below. When the domain is clear by the context, we shall drop the subscript in the L 2 (Ω)-norm, viz., · L 2 (Ω) ≡ · , for brevity. Setting v = u into (13) gives which, upon integration by parts, already yields, respectively, noting that (xu xy , u x ) = 0 = (xu y , u yy ) from divergence theorem and the decay properties of u as |(x, y)| → ∞. Young's and Cauchy-Schwarz' inequalities now imply Selecting α > 0, β = 4α 2 /9, γ = α 3 /3, and ǫ = 3/4, for instance, we deduce Another instructive choice is α > 0, β = α 2 , γ = α 3 and ǫ = 1/2, giving Thus, remarkably, using the non-standard test function v − ∇ · (A∇v), results into L providing coercivity with respect to u y also! This is a manifestation of the concept of hypocoercivity. As briefly discussed above, the mechanics of this idea can be traced back to the celebrated work of Hörmander [14] for sufficient conditions for hypoellipticity for second order linear operators: the first order term is "appropriately non-constant" with respect to the independent variables t, x, y, allowing for control also in the commutator vector field [∂ x , x∂ y ] = ∂ y . Informally speaking, the first order vector field in conjunction with the direction of ∂ x is able to "propagate" the x-variable dissipation of L to the whole of R 2 , resulting to the property of hypocoercivity [23].
Notice that the choice of α, β, γ yielding (14) corresponds to v − ∇ · (A∇v) being an elliptic differential operator, while the respective choice for (15) corresponds to the parabolic operator v − ∇ · (A∇v). A comment on the different scalings present on each term on the right-hand side of (14) is also in order. It is possible to introduce dependence of the tensor A on the time variable t for the study of the singularities near t = 0; we refer to [12,20] for details. Since our interest lies in the long time provable computation robustness for long times, we shall not pursue this direction further here.
On the other hand, kinetic PDE models posed on bounded domains Ω ⊂ R d , are often complemented with additional boundary /compatibility conditions to convey the physical properties in the vicinity of boundaries, e.g., reflection, diffusion, periodicity, etc. In the case of the classical Kolmogorov's equation considered here, the x-variable is typically associated with particle velocity, while the y-variable models position/displacement. Performing the above calculations on a bounded domain Ω instead, integrations by parts yield A sufficient condition, therefore, to arrive again at (14) or (15) is to assert that or, at least, that part of the left-hand side of (17) can be controlled by the non-negative terms on the right-hand side of (16) and the remaining to satisfy a respective positivity condition. For example, assuming that Ω : Therefore, assuming H 2 (Ω)-smoothness of the solution, imposing additional periodic boundary conditions on ∂ 0 Ω, we deduce since xn 2 = 0 on ∂ 0 Ω. Supplementing also with a no-flux boundary conditions on the inflow boundary ∂ − Ω, viz., we conclude since ∇u = (u x , u y ) T = 0 and n 1 = 0 on ∂ ± Ω. The periodicity in the x-variable and (19), therefore, implies We refer, e.g., to [21,6] for further discussion on additional boundary conditions. Different combinations of assumptions are also possible to yield (17), such as vanishing second derivative conditions across the elliptic boundary, etc. For the remaining of this work, we shall assume the validity of the additional boundary conditions: 1 2 xn 2 ∇u − n 1 ∇u x ∂Ω\∂ + Ω = 0 and n 1 ∇u x ∂ + Ω = 0, which trivially implies (17). We do so, instead of prescribing particular additional boundary conditions; we shall clearly state if the validity of (20) is assumed or not for the proof of various results below. We stress that, under minor modifications only, the theory presented below can still be valid for other suitable assumptions also by careful inspection of the proofs. This is not pursued here in the interest of a unified and clearer presentation.

A finite element method
Motivated by the above discussion, we return to the variational representation (12) of problem (10) on a bounded domain Ω ⊂ R 2 . For simplicity of the presentation we, henceforth, assume that Ω is polygonal also; extension to curved domains is possible via isoparametric versions of the finite element spaces discussed below. As is standard in discretisation by a Galerkin type finite element method, we introduce a family of triangulation of Ω, say T , consisting of mutually disjoint open triangular elements T ∈ T , whose closures coverΩ exactly. Let also h : ∪ T ∈T T → R + be the local meshsize function defined elementwise by h| T := h T := diam(T ). For simplicity, are further assume that T to be shape-regular, in the sense that the radius ρ T of the largest inscribed circle of each T ∈ T is bounded with respect to each element's diameter h| T , uniformly as h L∞(Ω) → 0 under mesh refinement. Also, we assume that T is locally quasi-uniform in the sense that the diameters of adjacent elements are uniformly bounded from above and below. Finally, let Γ := ∪ T ∈T ∂T denote the mesh skeleton, and Γ int := Γ\∂Ω. We define the finite element space subordinate to T by Finally, by ∇ T we shall denote the element-wise broken gradient. Starting from (12), we consider the spatially discrete finite element method: for every t ∈ (0, t f ], find U ≡ U(t) ∈ V h , such that for all V ∈ V h , whereby [v] 1 | e := (v| e ) + n + 1 + (v| e ) − n − 1 , for each internal face e and n 1 denoting the first component of the unit normal vector to e and [v] 1 | e := (v| e ) + n + 1 , for e ⊂ ∂Ω, for some τ : Γ → R nonnegative function to be determined precisely below and user-defined parameters κ, λ ≥ 0; we also set U(0) := Πu 0 ∈ V h , with Π : L 2 (Ω) → V h denoting the orthogonal L 2 projection onto the finite element space. For brevity, we shall use the notation: We also observe that the first term of s h vanishes for (x, y) T ∈ ∂Ω\(∂ − Ω ∩ ∂ + Ω) since xn 2 = 0 there. As we shall see below the stabilisation s h (·, ·) is constructed so that in the method (21) is consistent with respect to (12) upon assuming the additional boundary conditions (20).
for C τ > 0 large enough but independent of h and of p. Then, there exists a positive constant c 0 , independent of the approximate solution U ∈ V h and of t f , such that for any 0 < ǫ, ζ < 1, we have when f ≡ 0, we can trivially select ζ = 0.
Proof. Setting V = U in (21), standard arguments give upon observing that (xU y , U) = √ xn 2 U 2 ∂ + Ω /2 as U = 0 on ∂ − Ω. We calculate Observing now the identities Then, summing over all T ∈ T , we deduce Using (24) (25) Focusing on the fifth and sixth terms on the left-hand side of (25), recalling (11), the Cauchy-Schwarz' and Young's inequalities along with elementary manipulations imply respectively, upon defining the diagonal matrix B = diag(2(1 − ǫ), 2β − α 2 /(2ǫ)), for any 0 < ǫ < 1. We, then, have where λ min B−A denotes the smallest eigenvalue of B−A. Since U = 0 on ∂ − Ω which was assumed to be of positive one-dimensional Hausdorff measure, the validity of a Poincaré-Friedrichs/spectral gap inequality U ≤ C P F ∇U with C P F ≡ C P F (Ω) positive constant, independent of U, is assured, thereby allowing us to conclude the lower bound Also, the symmetry of A along with the Cauchy-Schwarz inequality and a standard inverse estimate yield, respectively, using the standard inverse inequality v 2 ∂T ≤ C inv p 2 /h T v 2 T on (29), and selecting ϕ 2 = τ ≥ 2C inv p 2 /{ {h} } and resorting to the local quasiuniformity of the mesh. Inserting the last two bounds into (25) and, using the Cauchy-Schwarz inequality, setting c 0 := C −1 P F λ min B−A , we deduce the result for any 0 < ζ < 1; when f ≡ 0 we can trivially select ζ = 0. Remark 4.2. Setting β = α 2 γ = α 3 , and ǫ = 1/2, we have det(B − A) = α 2 (1 − 2α). A numerical investigation reveals that max 0≤α≤1/2 λ min B−A ≈ 0.054429 attained for a value α ≈ 0.35060, noting that B − A has positive determinant for 0 < α < 1/2. It is possible to optimise further the constant c 0 by working as follows: starting from √ B∇U 2 , we have and selecting c 0 = λ min B−νA C −1 P F for suitable ν > 0. We will revisit this idea in Section 6. The last proof highlights that the validity of the Poincaré-Friedrichs inequality is of central importance here, and in general in the study of trend to equilibrium for kinetic equations [21,13,23]. In the present context of bounded domain spatial Ω, we have taken the viewpoint of using standard (non-weighted) integral norms, whose measure is present in the constant c 0 . The extension of the above framework to general Fokker-Planck equations on function spaces weighted by equilibrium distributions is an important question and will be discussed elsewhere. The availability of such a Poincaré-Friedrichs inequality (also known as spectral gap property in kinetic theory) facilitates the crucial feature of the proposed method: the absence of an exponential term of the form exp(ct f ) multiplying the data terms in the above stability estimate. As a consequence, (21) also immediately implies the decay of numerical solutions for the respective homogeneous problem.
for all t f > 0.
Proof. From (22) we deduce from which an application of Grönwall's Lemma already implies the result.

Lemma 5.1 ((Hypo)coercivity). Let
For the remaining of this work, unless explicitly stated, we assume that ǫ is away from 1 and that β ∼ α 2 , so that B = c B diag(1, α 2 ) for a given constant 0 < c B ≤ 1. Thus, we can take C 0 = (2c B ) −1 . Note that this also covers the case of A being singular, i.e., the operator ∇ · A∇(·) being parabolic. Next, we establish the consistency of the bilinear form B.

Lemma 5.3 ((hypo)elliptic projection). Assume that
for all V ∈ V h . Moreover, assuming further that v ∈ H k T (T ) for T ∈ T , k T ≥ 3, we have the approximation estimate: with s T := min{p + 1, k T } − 2 and x T := max x∈T {1, |x| 2 }, T ∈ T .
Proof. Uniqueness (and, therefore, existence due to linearity) follows from the coercivity of B shown in Lemma 5.1. Setting, now η := v − Pv and ξ := Pv − πv, with P denoting an optimal projection operator from H 1 −,0 (Ω) onto V h , whose particular properties apart from approximation are irrelevant at this point, we have and (∇ T (xη y ), A∇ξ) = T ∈T (η y + xη xy , αξ x + βξ y ) T + (xη yy , βξ x + γξ y ) T = α(η y , ξ x ) + β(η y , ξ y ) − T ∈T (η x , x(αξ x + βξ y ) y ) T + (η y , x(βξ x + γξ y ) y ) T + T ∈T ∂T xn 2 αη x ξ x + βη x ξ y + βη y ξ x + γη y ξ y ds. (36) Observe that the last term on the right-hand side of (36) is equal to and we note that the first term of (37) cancels with the first term of s h (η, ξ). Using now (35), (36) and (37) into (34), along with standard Cauchy-Schwarz' and Young's inequality arguments, and working as in (29) gives for any µ, κ, λ > 0, for some C > 1, independent of the, relevant to the argument, parameters. Alternatively, when κ, λ = 0, we work as follows: using the inverse estimate ξ x T and, similarly for ξ y , we further estimate the second term on the right-hand side of (37) to arrive instead at with |·| 2 denoting the matrix-2-norm. Applying the inverse estimate ∇v 2 T ≤ Cp 4 /h 2 T v T , for v ∈ P p (T ), on the terms containing µξ in the last two alternative estimates, selecting µ = h/p 2 , and performing standard manipulations, we arrive at the combined estimate which holds for any κ, λ ≥ 0. Choosing P to be an hp-optimal projection operator onto the finite element space V h , we can assume approximation estimates of the form for k > m, with u ∈ H k (T ), T ∈ T (see, e.g., [22,3] for examples of such operators). These, together with the trace estimate w 2 ∂T ≤ C w T w H 1 (T ) and the triangle inequality |ρ| ≤ |η| + |ξ| , already imply (33).
We now show the (potentially super-)approximation of for the left-hand side of (31). Such results are typical to finite element method for parabolic problems. We shall show that this is the case also for (21) discretising the, degenerate, Kolmogorov's equation.
for any t ∈ (0, t f ], for some positive constant C, depending only on the shape-regularity and the local quasi-uniformity of the mesh.
Proof. Integrating in time between 0 and t ∈ (0, t f ] multiplying by 2 and using standard arguments, gives recalling the notation c 0 = C −1 P F λ min B−A . Using now (28), we deduce the result. The above developments yield an a priori error bound in the norm appearing on the left-hand side of (38).
Theorem 5.5. Assume that u 0 , u, u t ∈ H 1 −,0 (Ω) ∩ H k T (T ) for k T ≥ 3, T ∈ T , the latter two for almost all t ∈ (0, t f ]. Assume also that (20) are satisfied. Then, under the assumptions of Lemmata 4.1 and 5.3, the error e := u − U of the finite element method (21) satisfies the bound: where with s T := min{p + 1, k T } − 2, and the constant C(A, κ, λ) > 0 is bounded away from +∞ for 0 < α ≤ 1/2 and all κ, λ ≥ 0 and independent of the mesh parameters and of u. Moreover, Proof. Working as for (28), we have Similarly, we also have The last two estimates can be further bounded from above by (33) and hp-approximation bounds for the orthogonal L 2 -projection [17,4], respectively. Also, from (33), we have Combining the above bounds with the triangle inequality recalling that u − u h = ρ + ϑ, the result already follows.
The crucial property of the above error estimate, i.e., the independence of the respective constant C(A, κ, λ) from the final time t f , is a direct consequence of the hypocoercivity-compatible discretisation (21).
Assumption u 0 ∈ H 1 −,0 (Ω) is a natural compatibility condition in this context. Since a typical setting in kinetic simulation concerns initial profiles with compact support within a computational domain Ω, u 0 ∈ H 1 −,0 (Ω) is of significant practical relevance also. Nonetheless, we shall now discuss the dependence of the error committed by the numerical method as t f → ∞ and we shall see that the effect of the initial condition error u 0 − Πu 0 diminishes exponentially as t f → ∞. Thus, possibly incompatible initial conditions u 0 ∈ H 1 (Ω)\H 1 −,0 (Ω) will have an exponentially diminishing effect in the accuracy of the method with respect to t f . where for s T := min{p + 1, k T } − 2, with C(A, κ, λ) independent of t f , of u, and of the mesh parameters (cf. Theorem 5.5).
Proof. Starting from (31), we use Lemmata 5.1 and 5.3 and, subsequently, apply (28) to arrive at upon multiplication by 2. Grönwall's Lemma, thus, implies As before, we also have The last bound combined with (43) and the completely analogous estimate for already imply the result.
Remark 5.7. The first term on the right-hand side of (43) vanishes upon altering the finite element method (21) at t = 0 from U(0) = Πu 0 to U(0) = πu 0 . In practical terms, this results to a computational overhead of a stiffness matrix solve as opposed to a mass matrix one.
The a priori error estimate derived in Theorem 5.5 is optimal with respect to mesh-size h, upon observing that the bilinear form B is the weak form of a formally 4th order differential operator. At the same time, the highest order terms involving A of B(w, v) play the role of "stabilisation" in this non-standard Galerkin context, resulting to increase of the spatial operator order from second to fourth. The estimate (40) is slightly suboptimal with respect to the polynomial degree p, as is typical in error analyses of hp-version interior penalty procedures, such as (21), involving inverse estimates [11]. Correspondingly, the a priori error bound in Theorem 5.6 is formally suboptimal with respect to the mesh size h by one order; the errors are measured in weaker norms that the ones appearing on the right-hand side of (42). We view this as a reasonable price to pay for the long time robustness of method when t f → ∞. Indeed, in many practical scenarios, such as long time computations of decay to equilibrium distributions, the favourable exponential dependence with respect to t f may be preferable to a slight suboptimality in the rate of convergence with respect to the mesh size h. Moreover, given the typically high smoothness of the exact solution u, the use of high order finite element spaces may diminish further the practical significance of this slight h-suboptimality.
One may be tempted, therefore, to scale the stabilisation terms involving A by appropriate powers of the mesh-size h to reduce the stiffness matrix scaling to one of a second order operator; this is a standard practice in classical Galerkin contexts, e.g., for streamline upwinded Petrov-Galerkin methods for convection-dominated problems. In the present context, however, such scaling would introduce new challenges, most important of which is that the spectral gap constant c 0 would be proportional to a negative power of h. Although such dependence may not turn out to be catastrophic in certain scenarios, e.g., when u t decays very fast (cf., (42) in Theorem 5.6), this is a challenging question that requires further research.
Another possibility to retrieve optimality in the rate of convergence for the method (21) in weaker norms is the incorporation of Aubin-Nitsche type duality techniques allowing for h-optimal error control of the (hypo)elliptic projection error u − πu with respect to L 2 -and/or H 1 -equivalent norms. Although this is classical in parabolic problems, it appears to be a challenge in the present context of hypocoercive operators. This is due to the non-closedness of the hypocoercivity property with respect to the adjoint operation in the present variational context (12). This constitutes an interesting direction for future research.
Finally, we remark on the practical performance of the proposed method. Preliminary computations show agreement with the theoretical developments presented above. Since this is, however, an interesting topic in its own right, it will be discussed in detail elsewhere [5].