The role of stabilization in the virtual element method: a survey

The virtual element method was introduced 10 years ago and it has generated a large number of theoretical results and applications ever since. Here, we overview the main mathematical results concerning the stabilization term of the method as an introduction for newcomers in the field. In particular, we summarize the proofs of some results for two dimensional ``nodal'' conforming and nonconforming virtual element spaces to pinpoint the essential tools used in the stability analysis. We discuss their extensions to several other virtual elements. Finally, we show several ways to prove interpolation estimates, including a recent one that is based on employing the stability bounds.


Introduction
The virtual element method was introduced ten years ago [9] as a generalization of the finite element method to polytopal meshes. Local virtual element spaces consist of solutions to local problems with polynomial data. Therefore, virtual element functions are not known in closed form; in the spirit of the mimetic finite differences [15,48], only the evaluation of their degrees of freedom is required in the design of the method.
Consequently, the bilinear forms appearing in the variational formulation of given partial differential equations are not computable and have to be approximated, based on two main ingredients: projections from local virtual element spaces onto polynomial spaces; bilinear forms that stabilize the scheme.
This entails that the error analysis for virtual elements has the form of a Strang-type result, where several variational crimes have to be taken into account. In particular, one has to cope with certain stability bounds and interpolation estimates in virtual element spaces.
The first work on the virtual element method contains the following statement concerning stability estimates [9,Section 4.6]: "In general, the choice of the bilinear form S K (·, ·) [the local virtual element stabilization] would depend on the problem and on the degrees of freedom. From (4.20) it is clear that S K (·, ·) must scale like a K (·, ·) [the "grad-grad bilinear form] on the kernel of Π ∇ p [an H 1 polynomial projector]. Choosing then the canonical basis ϕ 1 , . . . , ϕ NK as [χ i is the i-th local degree of freedom] the local stiffness matrix is given by In our case it is easy to check that, on a "reasonable" polygon, a K (ϕ i , ϕ j ) ≈ 1. Note that this holds true for all i = 1, 2, . . . , N K [N K is the dimension of the local discrete space] since we defined the local degrees of freedom suitably. [. . . ] However, several types of misbehaviour can occur for awkwardly-shaped polygons, in particular if two or more vertices tend to coalesce, although, in our numerical experiments, the method appears to be quite robust in this respect." So, some facts were made clear since the very inception of the method: • the stabilization is not required to have approximation properties: it only has to scale as the corresponding "continuous" bilinear form; • a reasonable choice for the stabilization is the "dofi-dofi" one given by S K (ϕ i , ϕ j ) = δ i,j for all i, j = 1, . . . , N K ; • the "dofi-dofi" stabilization needs to be carefully tuned/changed in presence of "awkwardlyshaped polygons".
Similar considerations are contained in other works tracing back to the early years of the virtual elements literature, say, the period 2013-2016: there, one can find heuristic motivations but no proofs of the stability bounds, i.e., bounds of the form for suitable discrete functions v h and positive constants α * < α * .
In 2017, the first paper [16] on the theoretical aspects of the stabilization in virtual elements was published: stability properties for two dimensional nodal conforming virtual elements were investigated on rather general geometries. Ever since, several related contributions have been proposed; amongst them we mention the other three pioneering works [28,29,39]. Most of the literature on this topic is concerned with nodal conforming virtual elements; fewer works can be found on other types of virtual elements, such as face, edge, Stokes, immersed-like virtual elements.
In light of this, the virtual element literature on the theoretical aspects of the stabilization might be classified into three periods: • the early years (2013-2016), when the properties of the stabilization were introduced and motivated heuristically; • the pioneering years (2017-2018), when the first papers [16,28,29,39] on the theoretical aspects of the stabilization were published; • the consolidation years (2019-2023), when several other works on the topic were written, the pioneering analysis was generalized, and other types of virtual elements were considered.
We deemed useful to collect and review all contributions that we are aware of about the theoretical aspects of the stabilization. We hope that this work could represent a starting guide in this area. More details and the topics we are going to cover can be found at the end of this section.
Disclaimer. As we are interested here only in reviewing the literature on theoretical aspects of the virtual element stabilization, we do not review contributions on more "practical" issues as in [14,45,49,52,54,55].
Notation and assumptions. Given D ⊂ R d , d = 1, 2, 3, a Lipschitz domain with measure |D|, we introduce the Sobolev space H s (D), s ≥ 0, and endow it with the usual norm · s,D , seminorm |·| s,D , and bilinear form (·, ·) s,D . If s = 0, we let the Sobolev space H 0 (D) be the Lebesgue space L 2 (D). Negative order Sobolev spaces are defined by duality.
We set the spaces P p (K) and P p (e) of polynomials of maximum degree p over a polygon K and an interval e. We use the notation P −1 (K) = {0}. For the space P p (K), given x K = (x K , y K ) and h K denote the centroid and diameter of K, it is convenient to introduce the basis {m α } of scaled and centred monomials (1) Analogously, we can define the basis {m e α }, α = 1, . . . , p, of shifted and scaled monomials for P p (e). Given two quantities a and b, we write a b meaning that there exists c depending on the shape of K, but not on its size, such that a ≤ c b.
Outline. We prove the stability bounds for two dimensional nodal conforming virtual elements in Section 2: here, we focus on two paradigmatic stabilizations; derive corresponding stability bounds; underline what the main tools in the proof are; review the relevant literature. We discuss the generalization to nonconforming "nodal" elements and arbitrarily regular "elliptic-type" elements in Section 3. "Nonelliptic-type" elements, such as face, edge, Stokes-like, and immersed virtual elements, are overviewed in Section 4. Comments on "p-version" spaces, and the role of the stabilization in relability and efficiency bounds for residual error estimators are given in Section 5. In Section 6, we show that the stability bound easily implies interpolation estimates. We present some conclusions and perspectives for the virtual element method in Section 7.
2 Basic results: stability in nodal conforming virtual elements In this section, we review the role of the stabilization in two dimensional nodal conforming virtual elements. After designing local virtual elements with a set of unisolvent degrees of freedom, we describe computable polynomial projectors and discrete bilinear forms in Section 2.1. We introduce two paradigmatic stabilizations and provide an easy proof of stability bounds in Section 2.2. In Section 2.3, we review the essential literature concerning the theory behind the stabilization in nodal conforming virtual elements.

Nodal conforming virtual elements
Following [9], given a polygonal element K and a positive integer number p, we define the nodal conforming virtual element We endow the space V h (K) with the following set of unisolvent degrees of freedom: • the point values of v h at the vertices of K; • on each edge e of K, the point values of v h at the p − 1 internal Gauß-Lobatto nodes; • given the scaled monomial basis {m α } of P p−2 (K) as in (1), the scaled moments The second set of degrees of freedom can be replaced by suitably scaled edge moments or uniformly spaced nodes. We collect the above degrees of freedom in the set The degrees of freedom of V h (K) allow for the computation of several polynomial projections [9]. In what follows, we need the operator Π ∇ p : and the operator Π 0 p−2 : L 2 (K) → P p−2 (K) defined as Given the degrees of freedom of a function v h in V h (K), we can compute Π ∇ p v h and Π 0 p−2 v h ; see [9]. Other options to fix the constant part of Π ∇ p , see the second condition in (3), can be found in the literature. For instance, it is alternatively possible to fix the average in the bulk or the arithmetic average of the values at the vertices of the polygon.
The standard discretization of the bilinear form a K (·, ·) := (∇·, ∇·) 0,K in nodal conforming virtual elements is given by The bilinear form is typically required to be computable via the degrees of freedom and satisfies the following stability bounds:

Two stabilizations for nodal conforming virtual elements
We introduce two stabilizations that are common in the literature of nodal conforming virtual elements. The first one is known as the "dofi-dofi" stabilization [9] and is given by The second one, which we shall refer to as "projected" stabilization, is [13] S Both stabilizations are computable using the degrees of freedom of V h (K).
Proof. We prove the assertion for the "projected" stabilization (8) only; the details for the "dofidofi" stabilization (7) are similar but slightly more involved. We point out that the assertion for the stabilization in (7) follows combining the bounds for the stabilization in (8) and the techniques, e.g., in [16,29,39].
The lower bound. We recall an inverse estimate for functions with polynomial; see, e.g., [32,Lemma 10]: Integrating by parts, recalling that ∆v h belongs to P p−2 (K), applying the Cauchy-Schwarz inequality, invoking the definition of negative Sobolev norms, and using the inverse estimate (9) and a polynomial inverse estimate on ∂K, the lower bound in (6) follows: The upper bound. Using the stability of Π 0 p−2 in the L 2 norm, the trace inequality, and the Poincaré inequality (recall we are assuming v h belongs to ker(Π ∇ p )), the upper bound in (6) follows: From the proof of Lemma 2.1, we can see that • the lower bound in (6) follows from (polynomial and virtual) inverse estimates, but no zero average condition (based on the fact that v h belongs also to ker(Π ∇ p )) is used; • the upper bound in (6) is proven based only on "direct estimates", such as the Poincaré inequality and the trace inequality, and therefore is valid for H 1 functions with zero average.
For these reasons, Lemma 2.1 immediately generalizes as follows.
Corollary 2.2. The bilinear forms S K (·, ·) in (7) and (8) satisfy the following stability bounds: where D is any subset with nonzero measure of either K or ∂K.
The stability bounds in Theorem 2.1 and Corollary 2.2 can be easily extended to the case of the so-called enhanced virtual element spaces introduced in [1]. More precisely, define This space can be endowed with the same degrees of freedom of the standard space virtual element space in (2). Such degrees of freedom allow to compute the projectors Π ∇ p and Π 0 p ; see (3) and (4). It is immediate to check that the stabilization satisfies Lemma 2.1 and Corollary 2.2

Early essential literature and later contributions
The first work containing mathematical proofs on the stabilization is [16]. This and other three works constitute the early backbone of the stability analysis of the virtual element method. In particular, we pinpoint and rapidly describe these four works, which were published in 2017-2018: • [16] contains the first analysis of the stabilization in two dimensional nodal conforming virtual elements; an arbitrary number of edges is allowed; stability bounds are there derived for different stabilizations (including the "dofi-dofi" one); • [28,29] contain the stability analysis for two and three dimensional nodal conforming virtual elements; arbitrary numbers of edges and faces are allowed; stability bounds are derived for different stabilizations (including the "dofi-dofi" one); • [39] contains the stability analysis for two dimensional nodal conforming virtual elements; more standard geometries are employed.
It is our opinion that the presentation in [39] is the simplest one and is therefore recommended to virtual elements novices; the presentation in [16] and [28,29] is more technical and applies to more general geometries, whence these contributions are recommended to "more expert readers".
After 2018, other works were devoted to additional mathematical aspects of the stability in nodal conforming virtual elements. Amongst others, we mention nodal conforming virtual elements on curved domains [24]; robustness with respect to anisotropic elements [33,34]; virtual elements stabilized by means of higher-order polynomial energy projection terms [25,26,37,47]; lack of robustness with respect to the degree of accuracy of the method [13]; nodal serendipity virtual elements [21]; presence of stabilization terms in residual error estimators [19,32]. We shall elaborate more on the two last topics in Section 5 below.
The literature stability analysis for other types of elements is more limited. The reason for this is that nodal conforming virtual element functions have a "second order elliptic structure". In other virtual elements, local spaces consist of functions solving local problems (with polynomial data) that are in general more technical to handle.

Stability in other "elliptic-type" virtual elements
In this section, we investigate the role of the stabilization for other "elliptic-type" virtual elements that still have an "elliptic structure". In particular, we focus on "second order" nonconforming virtual elements in Section 3.1; investigate their stability properties in Section 3.2; review conforming and nonconforming arbitrarily regular virtual element spaces in Section 3.3.

"Nodal" nonconforming virtual elements
Following [8], given a polygonal element K, E K its set of edges, and a positive integer number p, we define the "nodal" nonconforming virtual element We endow the space V h (K) with the following set of unisolvent degrees of freedom: • on each edge e of K, given the scaled monomial basis {m e α } of P p−1 (e) discussed in Section 1, the scaled moments 1 |e| e m e α v h ; • given the scaled monomial basis {m α } of P p−2 (K) as in (1), the scaled moments We collect the above degrees of freedom in the set {dof j } NK j=1 , N K being the dimension of V h (K). The degrees of freedom of V h (K) allow for the computation of several polynomial projections [8].
Given the degrees of freedom of a function v h in V h (K), we can compute Π ∇ p v h , Π 0 p−2 v h , and Π 0,e p−1 v h for all edges e; see [8]. The standard discretization of the bilinear form a K (·, ·) := (∇·, ∇·) 0,K in "nodal" nonconforming virtual elements is as in (5).
Also in the "nodal" nonconforming virtual element setting, the stabilization S K (·, ·) : V h (K) × V h (K) is a bilinear form that is computable via the degrees of freedom and is such that there exist 0 < α * ≤ α * independent of h K such that the stability bounds in (6) are valid.

Stability bounds in "nodal" nonconforming virtual elements
We introduce two stabilizations that are common in the literature on "nodal" nonconforming virtual elements. The first one is the "dofi-dofi" stabilization as in (7). The second one, which we shall refer to as "projected" stabilization, is Both stabilization are computable using the degrees of freedom of V h (K).
Proof. We prove only the assertion for the "projected" stabilization (12); to this aim, we follow the guidelines in [50, Theorem 3.2]. The details for the "dofi-dofi" stabilization (7) are similar but slightly more involved.
The lower bound. Integrating by parts, recalling the properties of functions in "nodal" nonconforming virtual element spaces in (11), applying the Cauchy-Schwarz inequality twice, using a polynomial inverse inequality on ∂K, invoking the Neumann trace inequality [53,Theorem A.33], and recalling the inverse estimate (9), we can write The upper bound. Using the stability of Π 0 p−2 and Π 0,e p−1 in the L 2 (K) and L 2 (e) (for all edges e of K) norm, respectively, the trace inequality, and the Poincaré inequality, the upper bound in (6) follows: As in the case of nodal conforming virtual elements, Lemma 3.1 can be generalized in the sense of Corollary 2.2. Other theoretical results on the stabilization for "nodal" nonconforming virtual elements can be found in [27].

Arbitrarily regular virtual element spaces
Arbitrarily regular virtual element spaces were one of the first generalization of the nodal conforming virtual element method detailed in Section 2.1 and trace back to almost the inception of the method; see [18,31]. Over the years, several conforming and nonconforming variants have been proposed; see, e.g., [4,5,56].
Differently from standard C 0 elements, arbitrarily regular virtual elements are designed so as global C k spaces, k > 0, can be constructed. This is accomplished by a suitable enrichment of interface and bulk degrees of freedom. In particular, local virtual element spaces are the set of solutions to suitable polyharmonic problems with polynomial data. So, the "elliptic structure" of the spaces can be still employed while deriving the stability bounds. Thus, the analysis is quite similar to that for standard "nodal" conforming and nonconforming virtual elements. We refer to [40] and [38] for the stability bounds in nonconforming and conforming arbitrarily regular virtual elements. It is quite interesting the fact that both references are rather recent (2020 and 2022, respectively). This is probably motivated by the fact that (i) stability bounds are derived in any dimension; (ii) even though the structure of the arbitrarily regular virtual element spaces is still "elliptic", the interaction of high order derivatives and the presence of multiple boundary conditions render the analysis of the stability bounds quite involved.

Face virtual elements
Face virtual elements were designed in [10,11,30] and generalize the Raviart-Thomas elements to polytopal meshes. Both in two and three dimensions, local face spaces consist of functions solving local "div-rot" problems and with polynomial normal components on faces. The degrees of freedom consist of bubble-like moments in the element and moments of normal components over faces.
The first stability analysis (on 2D curved elements) can be found in [42]. Lowest order two and three dimensional face spaces were analyzed in [20], general order two and three dimensional (standard and serendipity) face spaces in [22], and three dimensional face spaces on curved polyhedra in [43]. Related results are discussed in [57].
Compared to the elliptic case, the stability analysis is here complicated by the "div-rot" structure of the spaces. So, the analysis needs the design of certain Helmholtz-like decompositions, as well as employing more sophisticated results from the theory of vector potential and div-rot systems; see, e.g., [2,41].

Edge virtual elements
Edge virtual elements were first designed in [10] and generalize the Nédélec element to polytopal elements. In three dimensions, local edge spaces consist of functions solving local "div-curlcurl" problems inside the element; local "div-rot" problems on faces; polynomial tangential traces over edges. The degrees of freedom consist of bubble-like moments in the element and on faces, and moments of tangential components over edges.
The first stability analysis (lowest order, in two and three dimensions) can be found in [20]; the general order (standard and serendipity) case was later tackled in [22].
Compared to the the case of face elements, the analysis of edge elements in three dimensions is further complicated by the fact that virtual element functions are solutions to different types of problems inside the element ("div-curlcurl" problems) and on faces ("div-rot" problems), which require a different treatment to establish suitable stability estimates.

Stokes-like virtual elements
Lowest order Stokes-like virtual elements appeared in [3] and were later generalized to the general order case in [17]. Local spaces consist of functions solving local Stokes-like with polynomial data; in particular, the divergence of Stokes-like virtual element functions is polynomial. This allows for the design of methods that are automatically divergence free. In two dimensions, the degrees of freedom consist of bubble-like moments in the element and point values on the boundary; in three dimensions, of element and face moments, and point values on the edges.
Stokes-like virtual elements are very successful and were the object of a plethora of articles. However, the first and only paper dealing with the stability analysis of Stokes-like virtual element spaces is [23]. The stability bounds are essentially derived based on the stable inf-sup structure of the local Stokes problems, and polynomial and (novel) virtual element inverse estimates.

Immersed-like virtual elements
More recently, immersed-like virtual elements have been introduced in two and three dimensions for elliptic and Maxwell problems; see [35,36]. Such spaces are particularly effective when handling problems with nonsmooth coefficients. The irregular behaviour of the exact solution is captured by local immersed-like virtual element functions in a Trefftz fashion, as they are solutions to local problems involving the same nonsmooth coefficients as in the continuous problem.

Dependence on the degree of accuracy and error estimators
Two situations where the presence of the stabilization might create some issues are the p-and hp-versions of the method, and when proving the reliability and efficiency of error estimators. We elaborate on these two settings in Sections 5.1 and 5.2, respectively.

Stability bounds depending on the degree of accuracy
The p-and hp-versions of a Galerkin method aim at achieving convergence by increasing the dimension of the local approximation spaces (while keeping the mesh fixed), and by mesh refinement and p-version at once, respectively. Stability bounds with explicit dependence on the degree of accuracy p were investigated in several works. Amongst them, we recall the following contributions: p-explicit stability bounds were first derived in [13] for the "projected" stabilization in (8); this result was extended [6] to the "dofi-dofi" stabilization in (7); the p-version nonconforming virtual elements was addressed in [50].
In all cases, at least one of the constants α * and α * in (6) depend on p. This is not surprising as inverse estimates are employed to derive the lower bound in (6); see the lower bounds in 2.1 and 3.1. However, the suboptimality with respect to the degree of accuracy is much milder in practice; see [13, Section 4.1].

The stability in residual error estimators
In finite elements, an error estimator η given by the combination of local error estimators η K is a quantity that can be computed by means of the discrete solution u h to the method and has to be designed so as to be comparable to the error of the method.
For the Poisson problem, standard reliability and efficiency estimates read where ω K is the patch of elements around K and HOT D (f ) is a high-order oscillation term on a domain D involving the right-hand side f . Residual error estimators in virtual elements [19,32] satisfy analogous bounds that differ from (13) in two aspects: (i) the discrete solution u h is not available in closed form, so it must be replaced by a polynomial projection; (ii) the stabilization appears in the error estimator and error estimates.
So, for a computable virtual element error estimator η given by the combination of local virtual element error estimators η K , the virtual element counterpart of (13) has the following form: The presence of the extra stabilization term in the bounds above might be inauspicious. On the one hand, such term is not robust with respect to the degree of accuracy as discussed in Section 5.1.
On the other hand, adaptive mesh refinements may lead to polygonal meshes with several hanging nodes and possibly very small facets; in that case, one should be extremely careful in designing stabilization that are robust with respect to badly-shaped elements. A partial breakthrough facing these issues is contained in [12], where the stabilization term is removed in the reliability and efficiency estimates. However, this has been proved under restrictive assumptions: meshes must consist of elements with triangular shape; the 3D version is not covered; lowest order elements are employed. So, general reliability and efficiency estimates that are robust with respect to the stabilization are still unknown.

The stability bounds may imply interpolation estimates
Interpolation estimates by means of functions in the virtual element space are an important ingredient in the error analysis of the method. In this section, we first report standard ways of proving interpolation estimates in conforming and nonconforming virtual elements; see Sections 6.1 and 6.2, respectively. Next, in Section 6.3, we show a recent (faster and easier) strategy to derive interpolation estimates based on the stability bounds covering the conforming and nonconforming cases at once.

Standard proof for interpolation estimates: conforming elements
The first proof of interpolation estimates for nodal conforming virtual elements traces back to 2015; see [51]. Several variants have been proposed, most of them extending the ideas therein.
To explain the main ingredients to derive such estimates in the conforming case, we provide details on the two dimensional nodal conforming virtual elements; see [51,Proposition 4.2]. Given an element K split into a shape-regular triangulation T (K), v in H 1 (K), q p in P p (K), and v C the Clément quasi-interpolant of v over T (K), we define v I over K as the weak solution to A minimum energy argument entails whence interpolation estimates follow from polynomial quasi-interpolation and approximation estimates. So, deriving interpolation estimates for conforming virtual elements strongly hinges upon using the structure of the local virtual element spaces, i.e., the stability of the continuous formulation. When considering other types of virtual elements, this may result in interpolation estimates that are much harder to derive. Prototypical examples are interpolation estimates in Stokes-like [17]; Hellinger-Reissner-like [7,44]; face and edge [20,22] virtual elements.

Standard proof for interpolation estimates: nonconforming elements
The first proof of interpolation estimates for "nodal" nonconforming virtual elements traces back to 2018; see [50]. Several variants have been proposed, most of them extending the ideas therein.
To explain the main ingredients in deriving such estimates in the nonconforming case, we provide details on the two dimensional "nodal" nonconforming virtual elements; see [50, Proposition 3.1] or the later work [46,Corollary 4.1].
Given K in T h and v in H 1 (K), we define v I in the "nodal" nonconforming space V h (K) defined in (11) as the degrees of freedom interpolant of v. Notably, v I is identified by For any q K p in P p (K), recalling that ∆v I belongs to P p−2 (K) and n K · ∇v I |e belongs to P p−1 (e) for all the edges e of K, we readily deduce Extending this approach to other types of virtual elements may be not trivial.

A more recent proof for interpolation estimates
There is a recent alternative approach [23] to derive interpolation estimates, which applies to conforming and nonconforming virtual elements at once and is based on having stability bounds at hand.
Again, to simplify the presentation we stick to the two dimensional nodal conforming virtual element case, even though the following steps can be easily generalized to several other elements and dimensions. Given an element K, let S K (·, ·) be either defined in (7) or (8); v belong to H s (K), s > 1, for (7) and s ≥ 1 for (7); q p be any function in P p (K) with the same average of v over K. We denote the degrees of freedom interpolant of any sufficiently smooth function by adding a subscript · I to such function. Observe that q p = (q p ) I . Then, recalling Corollary 2.2 and using that S K (·, ·) is computable via the degrees of freedom, whence S K (v I , v I ) = S K (v, v), we readily get |v − v I | 1,K ≤ |v − q p | 1,K + |v I − q p | 1,K ≤ |v − q p | 1,K + α −1 * S K (v I − q p , v I − q p ) Thus, we proved that, upon having at disposal stability bounds of the form (6), the best interpolation error is bounded, up to constants, by the best polynomial error. The ideas above are very general and can be extended to several other types of elements. Notably, it is apparent that the same proof carries over identically to "nodal" nonconforming virtual elements.

Conclusions and outlook
The first paper on the virtual element method was published ten years ago. The analysis of the stabilization of the method is more recent and is still a current area of research. Several works on the stability in nodal conforming virtual elements have been published, fewer on other types of virtual elements. In this contribution, we made an overview on the literature about the role and theoretical analysis of the stabilization in the virtual element method. We presented paradigmatic proofs for "nodal" conforming and nonconforming virtual elements; we only reviewed the main ideas for other types of elements. We further underlined the fact that stability bounds imply interpolation estimates.