Quanling Deng
email: quanling.deng@math.wisc.edu

Alexandre Ern
email: alexandre.ern@enpc.fr

SoftFEM: revisiting the spectral finite element approximation of second-order elliptic operators

Keywords: finite element method (FEM), Laplacian, spectral approximation, eigenvalues, stiffness, gradient-jump penalty

We propose, analyze mathematically, and study numerically a novel approach for the finite element approximation of the spectrum of secondorder elliptic operators. The main idea is to reduce the stiffness of the problem by subtracting a least-squares penalty on the gradient jumps across the mesh interfaces from the standard stiffness bilinear form. This penalty bilinear form is similar to the known technique used to stabilize finite element approximations in various contexts. The penalty term is designed to dampen the high frequencies in the spectrum and so it is weighted here by a negative coefficient. The resulting approximation technique is called softFEM since it reduces the stiffness of the problem. The two key advantages of softFEM over the standard Galerkin FEM are to improve the approximation of the eigenvalues in the upper part of the discrete spectrum and to reduce the condition number of the stiffness matrix. We derive a sharp upper bound on the softness parameter weighting the stabilization bilinear form so as to maintain coercivity for the softFEM bilinear form. Then we prove that softFEM delivers the same optimal convergence rates as the standard Galerkin FEM approximation for the eigenvalues and the eigenvectors. We next compare the discrete eigenvalues obtained when using Galerkin FEM and softFEM. Finally, a detailed analysis of linear softFEM for the 1D Laplace eigenvalue problem delivers a sensible choice for the softness parameter. With this choice, the stiffness reduction ratio scales linearly with the polynomial degree. Various numerical experiments illustrate the benefits of using softFEM over Galerkin FEM.

Introduction

The optimal approximation of eigenvalues and eigenfunctions from second-order elliptic spectral problems by means of Galerkin finite element methods (FEM) is well-established. We refer the reader to the seminal contributions in Vainikko [START_REF] Vainikko | Asymptotic error bounds for projection methods in the eigenvalue problem[END_REF][START_REF] Vainikko | Rapidity of convergence of approximation methods in eigenvalue problems[END_REF], Bramble and Osborn [START_REF] Bramble | Rate of convergence estimates for nonselfadjoint eigenvalue approximations[END_REF], Strang and Fix [START_REF] Strang | An analysis of the finite element method[END_REF], Osborn [START_REF] Osborn | Spectral approximation for compact operators[END_REF], Descloux et al. [START_REF] Descloux | On spectral approximation. I. The problem of convergence[END_REF][START_REF] Descloux | On spectral approximation. II. Error estimates for the Galerkin method[END_REF], Babuška and Osborn [START_REF] Babuška | Eigenvalue problems[END_REF], and to the more recent reviews in [START_REF] Boffi | Finite element approximation of eigenvalue problems[END_REF][START_REF] Ern | Finite elements. II. Galerkin Approximation[END_REF]. The approximation of elliptic spectral problems has also been studied by means of mixed finite element methods [START_REF] Canuto | Eigenvalue approximations by mixed methods[END_REF][START_REF] Mercier | Eigenvalue approximation via non-conforming and hybrid finite element methods[END_REF][START_REF] Mercier | Eigenvalue approximation by mixed and hybrid methods[END_REF], discontinuous Galerkin methods [START_REF] Antonietti | Discontinuous Galerkin approximation of the Laplace eigenproblem[END_REF][START_REF] Giani | hp-adaptive composite discontinuous Galerkin methods for elliptic eigenvalue problems on complicated domains[END_REF], hybridizable discontinuous Galerkin methods [START_REF] Cockburn | Hybridization and postprocessing techniques for mixed eigenfunctions[END_REF][START_REF] Gopalakrishnan | Spectral approximations by the HDG method[END_REF], hybrid high-order methods [START_REF] Calo | Spectral approximation of elliptic operators by the hybrid high-order method[END_REF][START_REF] Carstensen | Guaranteed lower bounds on eigenvalues of elliptic operators with a hybrid high-order method[END_REF], and virtual element methods [START_REF] Gardini | Virtual element method for second-order elliptic eigenvalue problems[END_REF]. All of these methods deliver optimally convergent approximations. Since the eigenfunctions become more and more oscillatory in the upper part of the spectrum, their approximation is accurate only in the lower part of the spectrum. In contrast, isogeometric analysis [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] delivers a more accurate approximation in the upper part of the spectrum (see also [START_REF] Deng | Dispersion-minimized mass for isogeometric analysis[END_REF][START_REF] Calo | Dispersion optimized quadratures for isogeometric analysis[END_REF][START_REF] Deng | A boundary penalization technique to remove outliers from isogeometric analysis on tensor-product meshes[END_REF] for some recent improvements on the subject).

The goal of this work is to improve on the Galerkin FEM spectral approximation so as to increase the accuracy in the upper part of the spectrum. This goal is achieved by reducing the stiffness of the discrete spectral problem. With this in mind, we refer to the newly coined method as softFEM. The idea is to subtract a least-squares penalty on the gradient jumps across the mesh interfaces from the standard stiffness bilinear form. Thus, the softFEM bilinear form is defined as â(⋅, ⋅) ∶= a(⋅, ⋅) -ηs(⋅, ⋅),

where a(⋅, ⋅) is the standard Galerkin FEM stiffness bilinear form, η is the socalled softness parameter, and s(⋅, ⋅) is the bilinear form penalizing the gradient jumps across the mesh interfaces. The idea behind softFEM shares some common ground with isogeometric analysis where the basis functions have at least C 1 -smoothness. In softFEM, the same basis functions are used as in Galerkin FEM so that the smoothness is only C 0 . However, by considering the bilinear form â(⋅, ⋅) instead of a(⋅, ⋅), one reduces the amount of energy stored in the gradient jumps of eigenfunctions associated with the large eigenvalues in the spectrum. This change is not needed for eigenfunctions associated with the lower part of the spectrum since those eigenfunctions are smooth and can be accurately approximated on a given mesh. We notice that the bilinear form s(⋅, ⋅) has been considered for the purpose of stabilization (i.e., leading to a positive contribution and not to a negative one as in the present work) in various contexts related, in particular, to advection-dominated advection-diffusion equations and to the Stokes equations [START_REF] Burman | A unified analysis for conforming and nonconforming stabilized finite element methods using interior penalty[END_REF][START_REF] Burman | Edge stabilization for the generalized Stokes problem: a continuous interior penalty method[END_REF][START_REF] Burman | Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations[END_REF]. In the context of the Helmholtz equation, the bilinear form s(⋅, ⋅) is weighted by a coefficient with positive imaginary part to ensure coercivity [START_REF] Wu | Pre-asymptotic error analysis of CIP-FEM and FEM for the Helmholtz equation with high wave number. Part I: linear version[END_REF]. In addition, the possibility of using a weighting coefficient with negative real part has been considered in [START_REF] Burman | Linear continuous interior penalty finite element method for Helmholtz equation with high wave number: onedimensional analysis[END_REF][START_REF] Du | Preasymptotic error analysis of higher order FEM and CIP-FEM for Helmholtz equation with high wave number[END_REF] to improve the phase error. Incidentally, we mention that the term softFEM has been used recently in [START_REF] Peña | SoftFEM: the soft finite element method[END_REF] in a completely different context related to heuristic optimization and soft computing for solid mechanics.

To give the reader a first view on the benefits of softFEM over Galerkin FEM, we present in Figure 1 the relative eigenvalue and eigenfunction errors for the 1D Laplace eigenvalue problem (with Dirichlet boundary conditions) using Galerkin FEM and softFEM. We use a uniform mesh composed of N h = 100 elements and a polynomial degree p ∈ {1, 2, 3}. The total number of discrete eigenpairs is N h p ∶= pN h -1. The benefit of using softFEM is evident when looking at the upper part of the spectrum. Another salient advantage of softFEM with respect to Galerkin FEM is that softFEM tempers the condition number Laplace eigenvalue problem when using Galerkin FEM and softFEM with N h = 100 uniform elements and polynomial degrees p ∈ {1, 2, 3}. Upper row: p = 1; middle row: p = 2; bottom row: p = 3. The eigenfunction errors for linear Galerkin FEM and linear softFEM are the same as both discretization methods give the same eigenvectors (but not the same eigenvalues).

of the stiffness matrix. This can have practically important consequences in the context of explicit time-marching schemes for time-dependent PDEs by reducing the CFL constraint on the time step. In many situations we observe that the stiffness reduction ratio scales linearly with p and is of the order of 1 + p 2 . The main mathematical results of this work can be summarized as follows.

In Theorem 1 we show that in order to maintain the coercivity of the softFEM bilinear form â(⋅, ⋅), the softness parameter can be chosen so that η ∈ [0, η max ), where the limit value depends on the polynomial degree p and the type of mesh (tensor-product or simplicial). Specifically η max = 1 2p(p+1) on tensor-product meshes and η max = 1 2p(p+d-1) on simplicial meshes (here d ≥ 2 denotes the space dimension). This result is established by means of discrete trace inequalities with sharp constants. In Theorem 2 we establish that softFEM maintains the same optimal convergence rates as Galerkin FEM. In Theorem 4 we prove for the 1D Laplace eigenvalue problem approximated by linear softFEM (i.e., p = 1), that the choice η = 1 2(p+1)(p+2) = 1 12 leads to superconvergence of the eigenvalue errors (quartic convergence rate instead of quadratic). We retain this choice for the value of the softness parameter in the rest of this work and notice that it is compatible with the maximum value η max obtained in Theorem 1. Finally, in Theorem 3 we establish lower and upper bounds on the discrete softFEM eigenvalues by those approximated by Galerkin FEM. In particular, the lower bound shows that the optimal value for the stiffness reduction ratio should be 1 + p 2 on tensor-product meshes and 1 + p 4-d on simplicial meshes with d ∈ {2, 3}. Both values are close to those observed in our numerical experiments.

The rest of this paper is organized as follows. Section 2 presents the exact spectral problem, its Galerkin FEM discretization, the softFEM approximation, as well as the following salient results concerning softFEM: coercivity (Theorem 1), error estimates (Theorem 2), and lower and upper bounds on the discrete eigenvalues (Theorem 3). Theorem 2 and Theorem 3 are proved in Section 2, but the proof of Theorem 1 is postponed to Section 5. Section 3 is concerned with the softFEM approximation of the 1D Laplace eigenvalue problem on uniform meshes. It contains the superconvergence result for softFEM (Theorem 4) motivating the choice η = 1 2(p+1)(p+2) for the softness parameter, and numerical experiments for various polynomial degrees illustrating the benefits of using soft-FEM with respect to Galerkin FEM both for the accuracy of the upper part of the spectrum and for the stiffness reduction. Section 4 collects more challenging numerical examples (Laplace eigenvalue problem in multiple dimensions, elliptic eigenvalue problem and non-uniform meshes for the 1D Laplace eigenvalue problem, and the use of simplicial meshes on the unit square and the L-shaped domain still for the Laplace eigenvalue problem) which corroborate the positive conclusions drawn on softFEM in Section 3. In Section 5, we first study discrete trace inequalities with sharp constants and then use these inequalities to prove Theorem 1. Concluding remarks are presented in Section 6.

Main idea and results

In this section, we state the elliptic eigenvalue problem and describe its approximation by means of Galerkin FEM and softFEM. We then state the main results concerning softFEM.

Problem statement

Let Ω be a bounded, open subset of R d , d ≥ 1, with Lipschitz boundary ∂Ω. For simplicity, we assume in what follows that Ω is a polyhedron. We use standard notation for the Lebesgue and Sobolev spaces. For any measurable subset S ⊆ Ω, we denote the L 2 -inner product and norm as (⋅, ⋅) S and ⋅ S , respectively, and the same notation is used for vector-valued fields. For any integer m ≥ 1, we denote the H m -norm and H m -seminorm as ⋅ H m (S) and ⋅ H m (S) , respectively.

We consider the following second-order elliptic eigenvalue problem with homogeneous Dirichlet boundary conditions: Find an eigenpair

(λ, u) ∈ R + × H 1 0 (Ω) such that u Ω = 1 and -∇ ⋅ (κ∇u) = λu in Ω, u = 0 on ∂Ω, (2.1)
with the diffusion coefficient κ ∈ L ∞ (Ω) uniformly bounded from below away from zero, and we set κ min ∶= ess inf x∈Ω κ(x) > 0. For κ = 1, the problem (2.1) reduces to the Laplace (Dirichlet) eigenvalue problem. The variational formulation of (2.1) is

a(u, w) = λb(u, w), ∀w ∈ H 1 0 (Ω), (2.2)
with the bilinear forms

a(v, w) ∶= (κ∇v, ∇w) Ω , b(v, w) ∶= (v, w) Ω . (2.3)
The eigenvalue problem (2.1) has a countable set of eigenvalues λ j ∈ R + (see, for example, [32, Sec. 9.8])

0 < λ 1 < λ 2 ≤ λ 3 ≤ . . .
and an associated set of L 2 -orthonormal eigenfunctions u j , that is, (u j , u k ) Ω = δ jk , where δ jk is the Kronecker delta. With (2.2) in mind, the normalized eigenfunctions are also orthogonal in the energy inner product since we have

a(u j , u k ) = λ j b(u j , u k ) = λ j δ jk .
In what follows we always sort the eigenvalues in ascending order counted with their order of algebraic multiplicity.

Galerkin FEM

Let (T h ) h>0 be a shape-regular sequence of meshes of Ω. A generic mesh element is denoted τ , its diameter h τ , and its outward unit normal n τ . We set h ∶= max τ ∈T h h τ . To stay general, we consider both tensor-product meshes where the mesh elements are cuboids (and so is the domain Ω), and simplicial meshes where the mesh elements are simplices (triangles if d = 2, tetrahedra if d = 3). Let p ≥ 1 be the polynomial degree. Let P p (τ ) (resp., Q p (τ )) be the space composed of the restriction to τ of polynomials of total degree at most p (resp., of degree at most p in each variable). For tensor-product meshes, the Galerkin finite element approximation space is defined as

V h p ∶= {v h ∈ C 0 (Ω) ∶ v h ∂Ω = 0, ∀τ ∈ T h , v h τ ∈ Q p (τ )}, (2.4)
whereas for simplicial meshes, it is defined as

V h p ∶= {v h ∈ C 0 (Ω) ∶ v h ∂Ω = 0, ∀τ ∈ T h , v h τ ∈ P p (τ )}. (2.5)
It is well-known that in both cases V h p ⊂ H 1 0 (Ω). The Galerkin FEM approximation of (2.1) seeks

(λ h , u h ) ∈ R + × V h p such that u h Ω = 1 and a(u h , w h ) = λ h b(u h , w h ), ∀w h ∈ V h p . (2.6)
The algebraic realization of (2.6) follows by choosing basis functions {φ h j } j∈{1,...,N h p } of V h p with N h p ∶= dim(V h p ) (typically, one considers nodal basis functions). This leads to the following generalized matrix eigenvalue problem (GMEVP):

KU = λ h MU, (2.7)
where

K kl ∶= a(φ h k , φ h l ) and M kl ∶= b(φ h k , φ h l )
, for all k, l ∈ {1, . . . , N h p }, are the entries of the stiffness and mass matrices, respectively, and U ∈ R N h p is the eigenvector collecting the components of u h in the chosen basis.

SoftFEM

For all τ ∈ T h , we define h 0 τ to be the length of the smallest edge of τ if τ is a cuboid, whereas we set h 0 τ ∶= d τ ∂τ if τ is a simplex. Let F i h be the collection of the mesh interfaces. For all F ∈ F i h , we have F = ∂τ 1 ∩ ∂τ 2 for two distinct mesh elements τ 1 , τ 2 ∈ T h . We then set

h F ∶= min(h 0 τ1 , h 0 τ2 ), κ F ∶= min(κ τ1 , κ τ2 ), (2.8)
with κ τ ∶= ess inf x∈τ κ(x) (i.e., κ F is the smallest value of κ on the two elements that share the interface F ). Moreover, for any function v h ∈ V h p , we define the jump of its normal derivative across F as

⟦∇v h ⋅ n⟧ F ∶= ∇v h τ1 ⋅ n τ1 + ∇v h τ2 ⋅ n τ2 .
(2.9)

We drop the subscript F when the context is unambiguous. The softFEM approximation of (2.1) seeks ( λh , ûh

) ∈ R + × V h p such that ûh Ω = 1 and â(û h , w h ) = λh b(û h , w h ), ∀w h ∈ V h p , (2.10)
where for all

v h , w h ∈ V h p , â(⋅, ⋅) ∶= a(⋅, ⋅)-ηs(⋅, ⋅) with s(v h , w h ) ∶= F ∈F i h κ F h F (⟦∇v h ⋅n⟧, ⟦∇w h ⋅n⟧) F , (2.11
) and η ≥ 0 is a parameter to be specified below. The terminology softFEM is motivated by the fact that the term -ηs(⋅, ⋅) reduces the stiffness of the system. We refer to η as the softness parameter. We will see below that one can take η ∈ [0, η max ) for some η max depending on the polynomial degree p and the type of mesh elements so that the bilinear form â(⋅, ⋅) remains coercive. When η = 0, softFEM reduces to FEM.

Similarly to Galerkin FEM, the algebraic realization of the softFEM approximation (2.10) leads to the GMEVP K Û = λh M Û, (2.12) where K ∶= K -ηS with S kl ∶= s(φ h k , φ h l ), K and M are respectively the stiffness and mass matrices as in (2.7), and Û is the eigenvector collecting the components of ûh in the chosen basis {φ h j } j∈{1,...,N h p } of V h p .

Remark 1 (Variants). For p ≥ 2, the stiffness can be further reduced by imposing least-squares penalties on higher-order derivative jumps. However, these additional terms increase the computational costs while our numerical experiments (not shown for brevity) indicate only a further marginal improvement in terms of spectral errors. We also mention the recent work [START_REF] Deng | A boundary penalization technique to remove outliers from isogeometric analysis on tensor-product meshes[END_REF] which penalizes both the higher-order derivatives as well as the mass bilinear form near the boundary to eliminate the so-called outliers in isogeometric spectral approximations.

Main results on softFEM

In this section, we present our main results on softFEM. We first derive an upper bound on the softness parameter to ensure coercivity of the bilinear form â(⋅, ⋅). To improve readability, the proof is postponed to Section 5.

Theorem 1 (Coercivity). Let â(⋅, ⋅) be defined in (2.11). Set η max ∶=

1 2p(p+1)
for tensor-product meshes with d ≥ 1 and η max ∶= 1 2p(p+d-1) for simplicial meshes with d ≥ 2. Assume that the softness parameter η ∈ [0, η max ). The following holds:

β 1 w h 2 H 1 (Ω) ≤ â(w h , w h ), ∀w h ∈ V h p , (2.13)
with β 1 ∶= κ min (1 -η ηmax ) > 0. Let us now consider the convergence of eigenvalues and eigenfunctions for softFEM. We define the solution operator T ∶ L 2 (Ω) → H 1 0 (Ω) ⊂ L 2 (Ω) such that for all φ ∈ L 2 (Ω), a(T (φ), w) = b(φ, w), ∀w ∈ H 1 0 (Ω).

(2.14)

We notice that T is selfadjoint and compact, and the elliptic regularity theory implies that there is s ∈ ( 1 2 , 1] such that T maps boundedly from L 2 (Ω) into H 1+s (Ω). Moreover (λ, u) is an eigenpair of (2.2) if and only if (µ, u) is an eigenpair of T with µ = λ -1 .

Theorem 2 (Eigenvalue and eigenfunction errors

). Let (λ j , u j ) ∈ R + × H 1 0 (Ω) solve (2.2) and let ( λh j , ûh j ) ∈ R + × V h p solve (2.10) with the normalizations u j Ω = 1 and ûh j Ω = 1. Let s ∈ ( 1 2 , 1]
be the index of elliptic regularity. Assume that there is t ∈ [s, p] and a constant C t such that one has the following smoothness property:

φ H 1+t (Ω) + T (φ) H 1+t (Ω) ≤ C t φ Ω for all φ ∈ G j ∶= ker(µ j I -T ) with µ j ∶= λ -1 j .
Then, the following holds:

λh j -λ j ≤ Ch 2t , u j -ûh j H 1 (Ω) ≤ Ch t , (2.15)
where C is a positive constant independent of the mesh-size h. The convergence rates are optimal whenever t = p.

Proof. We cannot apply directly the classical theory for error analysis derived in [START_REF] Babuška | Eigenvalue problems[END_REF]Thm. 7.2 & 7.4] since the softFEM bilinear form â(⋅, ⋅) differs from a(⋅, ⋅). Instead, we can apply the extension of this theory presented in [START_REF] Ern | Finite elements. II. Galerkin Approximation[END_REF]Chap. 48] to finite element approximations with so-called variational crimes. We can work on the extended space Y h ∶= V h p + H 1+s (Ω) and establish the boundedness of â on Y h × Y h using the H 1 -seminorm augmented by s(⋅, ⋅) 1 2 . Optimal approximation properties in this norm are readily derived for smooth functions. Moreover, consistency holds true since we have s(u j , y) = 0 for all y ∈ Y h because s > 1 2 . This implies the above error estimates.

Remark 2 (Pythagorean identity). A classical identity relating the eigenvalue and eigenfunction errors (see, e.g., [START_REF] Strang | An analysis of the finite element method[END_REF]Chap. 6]) is

u j -ûh j 2 E = λ j u j -ûh j 2 Ω + λh j -λ j ,
where

⋅ 2 E ∶= â(⋅, ⋅) ≥ β 1 ⋅ 2 H 1 (Ω)
owing to Lemma 1. Our third main result quantifies the stiffness reduction by softFEM for one particular choice of the softness parameter η that is further motivated in Section 3 (see, in particular, Theorem 4), namely η = 1 2(p+1)(p+2) . Notice that η < η max = 1 2p(p+1) for tensor-product meshes and that η < η max = 1 2p(p+d-1) with d ∈ {2, 3} on simplicial meshes.

Theorem 3 (Eigenvalue lower and upper bounds). Assume that

η = 1 2(p+1)(p+2) . Assume that d ∈ {2, 3} if simplicial meshes are used. Let j ∈ N, let (λ h j , u h j ) ∈ R×V h p solve (2.6
), and let ( λh j , ûh j ) ∈ R×V h p solve (2.10) with the normalizations u h j Ω = 1 and ûh j Ω = 1. The following holds:

γ p λ h j ≤ λh j < λ h j , (2.16)
with γ p ∶= 2 p+2 on tensor-product meshes and γ p ∶= 4-d p+4-d on simplicial meshes. Proof. For all v h ∈ V h p ∖ {0}, let us define the Rayleigh quotients

R(v h ) ∶= a(v h , v h ) b(v h , v h ) , R(v h ) ∶= â(v h , v h ) b(v h , v h ) .
As shown in Section 5 (see (5.8)), we have

(1 -2p(p + 1)η)a(v h , v h ) ≤ â(v h , v h ) < a(v h , v h )
on tensor-product meshes and

(1 -2p(p + d -1)η)a(v h , v h ) ≤ â(v h , v h ) < a(v h , v h ) on simplicial meshes. With the choice η = 1 2(p+1)(p+2) , a direct calculation shows that γ p a(v h , v h ) ≤ â(v h , v h ) < a(v h , v h
), with γ p defined in the assertion, which readily implies that

γ p R(v h ) ≤ R(v h ) < R(v h
).

(2.17) Let V j denote the set of the subspaces of V h p of dimension j ≥ 1. Classical results on the Rayleigh quotient imply that

λ h j = min Ej ∈Vj max v h ∈Ej R(v h ), λh j = min Ej ∈Vj max v h ∈Ej R(v h ).
The bounds in (2.16) then readily follow from (2.17).

Since the stiffness matrices K and K are symmetric, their condition numbers are given by

σ ∶= λ h max λ h min , σ ∶= λh max λh min , (2.18)
where λ h max , λh max are the largest eigenvalues and λ h min , λh min are the smallest eigenvalue of the GMEVPs (2.7) and (2.12) that are associated with (2.6) and (2.10), respectively. We define the stiffness reduction ratio of softFEM with respect to Galerkin FEM as

ρ ∶= σ σ = λ h max λh max ⋅ λh min λ h min . (2.19)
In general, for Galerkin FEM and softFEM with sufficient elements (i.e., as h → 0), one has λ h min ≈ λh min . Thus, the stiffness reduction ratio depends only on the largest eigenvalues for both methods. Since softFEM leads to a smaller largest eigenvalue, softFEM lowers the condition number of the stiffness matrix, i.e., ρ ≥ 1. We define the asymptotic stiffness reduction ratio of softFEM with respect to Galerkin FEM as

ρ ∞ ∶= lim h→0 λ h max λh max .
(2.20)

Theorem 3 shows that for η = 1 2(p+1)(p+2) , the best possible asymptotic stiffness reduction ratio is 1+ p 2 on tensor-product meshes and 1+ p 4-d on simplicial meshes with d ∈ {2, 3}. Notice that for both types of meshes, this value grows linearly with p. Our numerical experiments reported in Section 3.2 for the 1D Laplace eigenvalue problem show that the asymptotic stiffness reduction ratio is indeed

ρ ∞ = 1 + p 2 .
Moreover, the values of the asymptotic stiffness reduction ratio observed in the more general situations studied in Section 4 are also close to the predictions of Theorem 3. Finally, we define the stiffness reduction percentage of softFEM with respect to Galerkin FEM as

= 100 σ - σ σ % = 100(1 -ρ -1 ) %, (2.21)
and the asymptotic stiffness reduction percentage as ∞ ∶= 100(1 -ρ -1 ∞ ) %, respectively.

Remark 3 (SoftFEM eigenvalues). It is well-known that for Galerkin FEM, one has λ j ≤ λ h j for all j ≥ 1, but this is not necessarily the case for softFEM. Our numerical experiments indicate that softFEM approximates the exact eigenvalues from above in the low-frequency region and from below in the high-frequency region.

Laplace eigenvalue problem in 1D

In this section, we focus on the spectral problem (2.1) with Ω ∶= (0, 1) and κ ∶= 1, that is, on the 1D Laplace eigenvalue problem. In this case, the problem (2.1) has exact eigenvalues and L 2 -normalized eigenfunctions

λ j = j 2 π 2 and u j (x) = √ 2 sin(jπx), j = 1, 2, . . . , (3.1)
respectively. We partition the interval Ω = (0, 1) into N h uniform elements so that the mesh size h = 1 N h . We first focus on the case of linear finite elements (p = 1) and derive some analytical results showing that in this case the optimal choice for the softness parameter is η = 1 12 , that is, η = 1 2(p+1)(p+2) for p = 1. Then we present numerical experiments for this choice of the softness parameter and various polynomial degrees.

Analytical results for linear softFEM

The advantage of using linear elements is that it is possible to compute analytically the eigenvalues and eigenvectors for Galerkin FEM and softFEM. Firstly, it is well-known that the bilinear forms a(⋅, ⋅) and b(⋅, ⋅) with p = 1 lead to the following stiffness and mass matrices:

K = 1 h ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2 -1 -1 2 -1 ⋱ ⋱ ⋱ -1 2 -1 -1 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , M = h ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2 3 1 6 1 6 2 3 1 6 ⋱ ⋱ ⋱ 1 6 2 3 1 6 1 6 2 3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (3.2)
which are of order (N h -1) × (N h -1). The bilinear form s(⋅, ⋅) leads to the matrix

S = 1 h ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5 -4 1 -4 6 -4 1 1 -4 6 -4 1 ⋱ ⋱ ⋱ ⋱ ⋱ 1 -4 6 -4 1 1 -4 6 -4 1 -4 5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (3.3)
which is also of order (N h -1) × (N h -1). Recall that we then have K ∶= K -ηS, and that according to Theorem 1, we must take the softness parameter η ∈ [0, η max ) with η max = 1 2p(p+1) = 1 6 since p = 1 here. Lemma 1 (Analytical eigenvalues and eigenvectors). (i) Galerkin FEM approximation: The GMEVP KU = λ h MU has eigenpairs (λ h j , U j ) for all j ∈ {1, . . . , N h -1} with

λ h j = 6 h 2 1 -cos(t j ) 2 + cos(t j ) , U j = c j sin(kt j ) k∈{1,...,N h -1} , (3.4)
with t j ∶= jπh and some normalization constant c j > 0. (ii) SoftFEM approximation: The GMEVP K Û = λh M Û has eigenpairs ( λh j , Ûj ) for all j ∈ {1, . . . , N h -1} with An interesting consequence of (3.4)-(3.5) is that for linear softFEM, the stiffness reduction ratio and the asymptotic stiffness reduction ratio are

λh j = 6 h 2 1 + 3η -(1 + 4η) cos(t j ) + η cos(2t j ) 2 + cos(t j ) , Ûj = U j . ( 3
ρ = λ h max λh max ⋅ λh min λ h min = 5 + cos(πh) 5 -cos(πh) , ρ ∞ = lim h→0 5 + cos(πh) 5 -cos(πh) = 3 2 . (3.6)
Thus, asymptotically, linear softFEM reduces the stiffness of Galerkin FEM by about 33.3%. For all η ∈ [0, η max ) with η max = 1 2p(p+1) = 1 6 , Theorem 2 shows that one should expect a quadratic convergence rate for the discrete eigenvalues. We now show that for the specific choice η = 1 2(p+1)(p+2) = 1 12 , one obtains a quartic convergence rate, uniformly for all the discrete eigenvalues.

Theorem 4 (Eigenvalue superconvergence). Let λ j be the j-th exact eigenvalue of (2.1) and let λh j be the j-th approximate eigenvalue using linear softFEM. Assume that η = 1 2(p+1)(p+2) = 1 12 . The following holds:

λh j -λ j λ j < 1 360 (jπh) 4 , ∀j ∈ {1, . . . , N h -1}. (3.7)
Proof. The exact eigenvalues λ j are given in (3.1), and the approximate eigenvalues λh j are given in (3.5). To motivate the result of Theorem 4, we observe that applying a Taylor expansion to λh j , we obtain (recall that t j ∶= jπh)

λh j -λ j λ j = 1 -12η 12 t 2 j + 1 360 t 4 j - 17 -84η 60480 t 6 j + O(t 8 j ),
showing that the choice η = 1 12 leads to a cancellation of the dominant term in the expansion and that the sixth-order term has a negative coefficient. More rigorously, using (3.1), (3.5), and algebraic manipulations, we infer that

λh j -λ j λ j = 9 -(2 + cos(t j )) 2 t 2 j (2 + cos(t j )) -1 .
Since t j samples the interval (0, π), we can consider a continuous variable t ∈ (0, π) and prove more generally that

9 -(2 + cos(t)) 2 t 2 (2 + cos(t)) -1 < 1 360 t 4 ,
or, equivalently, that

-t 6 (2 + cos(t)) < 3240 -360(2 + cos(t)) 2 -360t 2 (2 + cos(t)) < t 6 (2 + cos(t)),
for all t ∈ (0, π). For the first inequality, we notice that the function

f (t) ∶= t 6 (2 + cos(t)) + 3240 -360(2 + cos(t)) 2 -360t 2 (2 + cos(t))
is increasing on (0, t 0 ) and decreasing on (t 0 , π) with t 0 ≈ 2.79911, that f (0) = 0 and f (π) ≈ 0.800921. The minimum value of f in (0, π) is thus f (0) = 0. For the second inequality, we notice that the function

g(t) ∶= t 6 (2 + cos(t)) -3240 + 360(2 + cos(t)) 2 + 360t 2 (2 + cos(t))
is increasing on (0, π) and that g(0) = 0. This completes the proof.

Remark 4 (Literature). The same value for the penalty parameter to achieve superconvergence is obtained in [START_REF] Burman | Linear continuous interior penalty finite element method for Helmholtz equation with high wave number: onedimensional analysis[END_REF] for the Helmholtz problem under Robin boundary conditions, still for p = 1. Values of the penalty parameter for p ∈ {2, 3} are derived in [START_REF] Du | Preasymptotic error analysis of higher order FEM and CIP-FEM for Helmholtz equation with high wave number[END_REF].

Numerical results for arbitrary-order softFEM in 1D

In this section, we explore numerically softFEM for various polynomial degrees p ≥ 1 using in all cases the softness parameter η = Recall that Figure 1 shows the relative eigenvalue and eigenfunction errors for Galerkin FEM and softFEM with N h = 100 uniform elements and polynomial orders p ∈ {1, 2, 3}. Notice that there are N h p ∶= pN h -1 eigenpairs both for Galerkin FEM and for softFEM. We refer the reader to Section 3.3 for a brief discussion on the structure of the discrete spectrum for Galerkin FEM, including the notions of acoustic/optical branches and stopping bands. The improvement offered by softFEM over Galerkin FEM for the eigenvalues is clearly visible in Figure 1 over the whole spectrum. For the eigenfunctions, there is no difference for p = 1 (see Lemma 1), whereas the improvement of softFEM over Galerkin FEM for p ∈ {2, 3} is salient around the stopping bands (that is, around j = N h for p = 2 and around j ∈ {N h , 2N h } for p = 3). Incidentally, we notice that for the H 1 -seminorm, the errors in the low-frequency region are slightly larger with softFEM than with Galerkin FEM, although the convergence order for softFEM remains optimal. This is expected since in the low-frequency region, best-approximation errors in the finite element space decay optimally, and the softFEM approximation leads to an additional optimally-converging contribution due to the interface jump penalty on the normal gradient. Table 1 reports the errors for the first and sixth eigenpairs using softFEM and polynomial degrees p ∈ {1, . . . , 4}. We observe that in all the cases, the convergence rates match well the predictions of Theorem 2 (and of Theorem 4 for p = 1). To motivate the choice of the softness parameter η = 1 2(p+1)(p+2) = 1 24 for p = 2, we show in Figure 2 the softFEM discrete spectra using various values for the softness parameter η. In this experiment, we increase the mesh resolution to N h = 1000 elements. In the left panel of Figure 2, for the sake of illustration, we actually increase the stiffness, i.e., we set â ∶= a + ηs. As expected, increasing η merely worsens the results. Instead, in the right panel of Figure 2, we return to softFEM and consider â ∶= aηs. We observe that the choice η = 1 24 appears to deliver the best overall result concerning the accuracy of the discrete eigenvalues over the whole spectrum. In the high-frequency region, the accuracy of the discrete eigenvalues is sensitive to the value of the softness parameter. For reference, we also display the results for η = η max = 1 2p(p+1) = 1 12 which show that the limit value on the softness parameter derived in Theorem 1 is indeed sharp.

1 2(p+1)(p+2) . p N h λh 1 -λ1 λ1 u 1 -ûh 1 H 1 u 1 -ûh 1 L 2 λh 6 -λ6 λ6 u 6 -ûh 6 H 1 u 6 -ûh 6 L 2 8
In Figure 3, we present the ratio ηs(û h j , ûh j ) a(û h j , ûh j ) for softFEM eigenfunctions. The mesh is composed of 240, 120, 80, 60 uniform elements for p ∈ {1, . . . , 4}, respectively, so that the number of eigenpairs is always the same. As predicted by Theorem 1, this ratio is always lower than one. We see that the amount of stiffness removed by softFEM is more substantial in the highfrequency region.

Table 2 shows the minimal and maximal eigenvalues, the condition numbers, the stiffness reduction ratios, and the percentages for Galerkin FEM and .) We observe that the stiffness reduction ratio increases with the polynomial degree, starting at ρ = 1.5 for p = 1 up to ρ = 3.2476 for p = 5. Thus, the benefit of using softFEM in tempering the condition number of the stiffness matrix becomes more pronounced as p is increased. We also notice that the computed value for the stiffness reduction ratio ρ is quite close to the optimal value 1 + p 2 resulting from Theorem 3 (see the lower bound in (2.16)).

Discrete spectrum for Galerkin FEM

The goal of this section is to briefly outline some basic facts about the spectrum of Galerkin FEM for the 1D Laplace eigenvalue problem. We explore the polynomial degrees p ∈ {1, 2, 3}. For p = 1, all the degrees of freedom (dofs) in V h p are attached to the N h p mesh vertices. Letting Λ ∶= λh 2 , solving the GMEVP leads us to look for nonzero vectors in the kernel of the following matrix of order

(N h -1) × (N h -1): A vv ∶= Λ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4 1 1 4 1 ⋱ ⋱ ⋱ 1 4 1 1 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ -6 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2 -1 -1 2 -1 ⋱ ⋱ ⋱ -1 2 -1 -1 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 3.8)
For p = 2, there are N h 2 = 2N h -1 dofs. It is interesting to order first the N h -1 dofs associated with the mesh vertices and then the N h dofs associated with the mesh elements. The basis functions associated with these N h dofs are bubble functions supported in a single mesh element. Solving the GMEVP problem leads us to look for nonzero vectors in the kernel of the following matrix whose block decomposition reflects the above partition into vertex and bubble dofs:

A vv 0 A bv A bb . (3.9)
It turns out that there is one vector in the kernel of A bb whose bubble dofs oscillate from one cell to the next one, and the corresponding eigenvalue is λ b = 10h -2 . The other vectors are obtained by considering the kernel of the block A vv which admits the following structure:

A vv ∶=Λ 2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 6 -1 -1 6 -1 ⋱ ⋱ ⋱ -1 6 -1 -1 6 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ -16Λ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 13 1 1 13 1 ⋱ ⋱ ⋱ 1 13 1 1 13 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + 240 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2 -1 -1 2 -1 ⋱ ⋱ ⋱ -1 2 -1 -1 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (3.10)
Finally, for p = 3, there are N h 3 = 3N h -1 dofs. We order first the N h -1 dofs associated with the mesh vertices and then the 2N h dofs associated with the mesh elements. The basis functions associated with these 2N h dofs are bubble functions supported in a single mesh element (2 per element). Solving the GMEVP problem leads us to look for nonzero vectors in the kernel of a matrix with the same block-structure as in (3.9), but this time the block A bb is two times larger. The kernel of A bb is two-dimensional and the corresponding eigenfunctions are thus composed only of bubble functions. Moreover, we have

A vv ∶=Λ 3 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 8 1 1 8 1 ⋱ ⋱ ⋱ 1 8 1 1 8 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ -30Λ 2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 36 -1 -1 36 -1 ⋱ ⋱ ⋱ -1 36 -1 -1 36 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + 360Λ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 64 3 3 64 3 ⋱ ⋱ ⋱ 3 64 3 3 64 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ -25200 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2 -1 -1 2 -1 ⋱ ⋱ ⋱ -1 2 -1 -1 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .
(3.11) For all p ∈ {1, 2, 3}, one can readily verify that the matrix A vv has a non-trivial kernel if and only if Λ is a root of the following polynomials (the subscript refers to the polynomial degree):

f 1 (Λ) = (2 + ζ j )Λ -6(1 -ζ j ), f 2 (Λ) = 2(3 -ζ j )Λ 2 -16(13 + 2ζ j )Λ + 480(1 -ζ j ), f 3 (Λ) = (4 + ζ j )Λ 3 -30(18 -ζ j )Λ 2 + 360(32 + 3ζ j )Λ -25200(1 -ζ j ), (3.12)
where ζ j ∶= cos(πt j ), t j ∶= jh and j ∈ {1, . . . , N h -1}. By replacing ζ j by the continuous variable ζ ∶= cos(πt) with t ∈ (0, 1), one obtains one branch of eigenvalues for p = 1, two branches of eigenvalues for p = 2, and three branches of eigenvalues for p = 3. Each branch contains N h -1 eigenvalues. For p ∈ {2, 3}, the spectrum is completed by the one or two eigenvalues associated with the eigenfunction(s) composed of bubble functions only.

In the literature, one refers to these latter eigenvalues as stopping band(s), whereas the branch associated with the lowest eigenvalues is called acoustical branch and the other branches are called optical branches. For instance, [START_REF] Brillouin | Wave propagation in periodic structures[END_REF] reported that quadratic finite elements for the 1D Laplace eigenvalue problem delivered an acoustical branch (low-frequency region) and an optical branch (high-frequency region) separated by one stopping band. We refer the reader to the left plots in Figure 1 for an illustration of these notions. We also observe that the notions of acoustical and optical branches as well as stopping bands depend on the sorting of the eigenvalues and that some overlap between the branches can happen in multiple dimensions; see Figure 4 for an illustration in 2D.

SoftFEM on more challenging numerical examples

In this section, we present more challenging numerical tests to illustrate the performances of softFEM. We consider Laplace eigenvalue problems on tensorproduct meshes in Section 4.1, elliptic eigenvalue problems and non-uniform meshes in 1D in Section 4. and 4.3, we use a higher-order method with a large number of elements to produce reference eigenpairs so as to quantify the approximation errors.

Laplace eigenvalue problems on tensor-product meshes

We consider the spectral problem (2.1) posed on Ω = (0, 1) d , d ∈ {2, 3}, with κ = 1. For d = 2, the exact eigenvalues and eigenfunctions are respectively for all i, j = 1, 2, . . .,

λ ij = (i 2 + j 2 )π 2 , u ij (x, y) = c ij sin(iπx) sin(jπy),
for some normalization constant c ij > 0, whereas for d = 3, the exact eigenvalues and eigenfunctions are respectively for all k, l, m = 1, 2, . . .,

λ klm = (k 2 + l 2 + m 2 )π 2 , u klm (x, y, z) = c klm sin(kπx) sin(lπy) sin(mπz),
for some normalization constant c klm > 0. For the Galerkin FEM and softFEM approximation, we use uniform tensor-product meshes. Theorem 1 shows that admissible values for the softness parameter are η ∈ [0, η max ) with η max = 1 2p(p+1) . Motivated by the 1D numerical experiments reported Section 3, we take again η = 1 2(p+1)(p+2) . Figures 5 and6 show the relative eigenvalue errors when using quadratic and cubic Galerkin FEM and softFEM in 2D. For quadratic elements, we use a uniform mesh with 40 × 40 elements, whereas for cubic elements, we use a uniform mesh with 20 × 20 elements. Figure 7 shows the relative eigenvalue errors for the 3D problem with 20 × 20 × 20 elements and p ∈ {2, 3, 4}. We observe in these plots that softFEM significantly improves the accuracy in the high-frequency region. Moreover, the plots using the log-log scale indicate that the spectral accuracy is maintained for quadratic elements and even improved for cubic elements in the low-frequency region. The convergence rates for the errors are optimal, and we omit them for brevity.

Figure 8 shows the ratio ηs(û h j , ûh j ) a(û h j , ûh j ) for the softFEM eigenfunctions in both the 2D and 3D settings. In 2D, there are 48 × 48, 24 × 24, 16 × 16, and 12 × 12 uniform elements for p ∈ {1, . . . , 4}, respectively, whereas in 3D, there are 24 × 24 × 24, 12 × 12 × 12, 8 × 8 × 8, and 6 × 6 × 6 uniform elements for p ∈ {1, . . . , 4}, respectively. These results essentially show how much stiffness is removed from the eigenfunctions by means of softFEM. The fact that the ratio ηs(û h j , ûh j ) a(û h j , ûh j ) is more pronounced in the high-frequency region corroborates the reduction of the spectral errors in this region. Finally, we mention that the stiffness reduction ratios and percentages are quite close to those reported in 1D, that is, ρ ≈ 1 + p 2 and ≈ 100 p p+2 % for p ∈ {1, . . . , 4} in both 2D and 3D.

Elliptic eigenvalue problems and non-uniform meshes in 1D

We now consider the 1D elliptic eigenvalue problem (2.1) with κ(x) ∶= e x sin(2πx) , so that κ max ≈ 1.34 and κ min ≈ 0.46. The exact eigenpairs are approximated using Galerkin FEM with C 6 septic B-spline basis functions and a mesh composed of N h = 1000 elements. Figure 9 compares the relative eigenvalue errors for Galerkin FEM and softFEM on a uniform mesh composed of N h = 200 elements and polynomial degrees p ∈ {2, 3, 4, 5}. We observe that softFEM reduces the spectrum errors, especially in the high-frequency region. The convergence rates for the errors are optimal, and we omit them here for brevity. Table 3 shows the smallest and largest eigenvalues, the condition numbers, the stiffness reduction ratios, and the percentages for Galerkin FEM and soft-FEM. In all cases, we observe that softFEM leads to smaller largest eigenvalues and hence to smaller condition numbers. The stiffness reduction ratio is about ρ = σ σ ≈ 1 + p 2 while the percentage is about ≈ 100 p p+2 %; this is consistent with the 1D results reported in Section 3.2.

Figure 10 compares the relative eigenvalue errors for the 1D Laplace eigenvalue problem when using Galerkin FEM and softFEM with p ∈ {2, 3, 4, 5} on a non-uniform mesh composed of N h = 10 elements. The mesh nodes have been randomly set to {0, 0.1, 0.18, 0.29, 0.41, 0.5, 0.59, 0.66, 0.81, 0.92, 1}. Refer- ence eigenvalues to evaluate the errors are computed as above. We observe that the improvement offered by softFEM over Galerkin FEM is similar to the one observed on uniform meshes. Table 4 reports the smallest and largest eigenvalues, the condition numbers, the stiffness reduction ratios, and the percentages. We observe that the stiffness reduction ratios are slightly larger than when using uniform meshes (compare with Table 3).

Simplicial meshes and L-shaped domain

In this section, we consider the 2D Laplace eigenvalue problem posed on the unit square domain or on the L-shaped domain, and we use simplicial meshes (triangulations) as depicted in Figure 11. Theorem 1 shows that admissible values for the softness parameter on simplicial meshes are η ∈ [0, η max ) with

η max = 1 2p(p+d-1) = 1 2p(p+1) if d = 2
. Motivated by the 1D numerical experiments reported in the previous section, we take again η = 1 2(p+1)(p+2) . Figures 12 and13 compare the relative eigenvalue errors for the 2D Laplace eigenvalue problem on the unit square domain and the L-shaped domain, respectively, when using Galerkin FEM and softFEM with p ∈ {1, 2, 3} and an unstructured mesh (triangulation). As observed in the previous numerical experiments, softFEM leads to smaller spectral errors than Galerkin FEM especially in the high-frequency region. Tables 5 and6 report the smallest and largest eigenvalues, the condition numbers, the stiffness reduction ratios and the percentages for the 2D Laplace eigenvalue problem on the unit square domain and the L-shaped domain, respectively. We use Galerkin FEM and softFEM with p ∈ {1, 2, 3} and an unstructured mesh. Once again we observe that softFEM is capable to reduce significantly the stiffness of the resulting matrix on unstructured meshes as well.

Proof of Theorem 1

In this section, we prove Theorem 1 which establishes the coercivity of the bilinear form â(⋅, ⋅) under the condition that the softness parameter η ∈ [0, η max ) for some real number η max depending on the polynomial degree p and the type of mesh. To this purpose, we first establish some useful discrete trace inequalities.

Discrete trace inequalities

For a natural number m ∈ N, we define the sets I m ∶= {0, . . . , m}, ∂I m ∶= {0, m}, and I 0 m ∶= I m ∖ ∂I m . Let p ∈ N be the polynomial degree. We are going to consider the Gauss-Lobatto rule with (p + 2) points (see, for example, [START_REF] Quarteroni | Numerical mathematics[END_REF]), which is exact for polynomials of degree at most (2p + 1).

, j = 2 (p + 1)(p + 2)(L p+1 (ξ j )) 2 , ∀j ∈ I 0 p+1 , (5.1
) where L p+1 is the Legendre polynomial of degree (p + 1). For a univariate function v that is k-times differentiable, we denote its k-th derivative as v (k) .

v (k) ∂τ ≤ C 1 (k, p)h -1 2 τ v (k) τ , ∀v ∈ P p (τ ), (5.2)
where C 1 (k, p) ∶= (p -k + 1)(pk + 2). Moreover, the constant C 1 (k, p) is sharp. In particular, for p ≥ 1 and k = 1, we have Proof. It is clear that it suffices to prove (5.2) for k = 0. Let v ∈ P p (τ ). We observe that v 2 ∂τ = v(a)

v ′ ∂τ ≤ p(p + 1)h -1 2 τ v ′ τ , ∀v ∈ P p (τ ). ( 5
+ j∈I 0 p+1 j v 2 (g j ) = b -a (p + 1)(p + 2) v 2 ∂τ + b -a 2 j∈I 0 p+1 j v 2 (g j ) ≥ C 1 (0, p) -2 h τ v 2 ∂τ ,
where we used that g 0 = a and g p+1 = b, the definition of C 1 (0, p) and h τ , and the fact that the weights j are non-negative for all j ∈ I 0 p+1 . This proves (5.1) for k = 0. Finally, that the inequality is sharp follows from the fact that it is possible to find a nonzero polynomial in P p (τ ) that vanishes at all the points g j for all j ∈ I 0 p+1 . Let us now turn to the multi-dimensional case. We consider first the tensorproduct case. For simplicity, we focus on bounding the normal derivative on the boundary of a cuboid cell. For a different result bounding any partial derivative on the boundary, we refer the reader to Remark 5.

Lemma 3 (Discrete trace inequality, cuboid)

. Let τ ∶= [a 1 , b 1 ] × . . . × [a d , b d ] ⊂ R d ,
with b j > a j for all j ∈ {1, . . . , d}, be a cuboid with boundary ∂τ and outward normal n τ . Recall that h 0 τ ∶= min i∈{1,...,d} (b i -a i ) is the length of the smallest edge of τ . Let p ≥ 1. The following holds:

∇v ⋅ n τ ∂τ ≤ p(p + 1)(h 0 τ ) -1 2 ∇v τ , ∀v ∈ Q p (τ ). (5.4)
Moreover, the constant is sharp. Notice that (5.4) coincides with (5.3) for d = 1.

Proof. We present the proof in the 2D case (d = 2); the general case is treated similarly. Let v ∈ Q p (τ ). One can write v(x, y) = ∑ jx,jy∈Ip α jxjy ψ jx (x)ψ jy (y), where {ψ j } j∈Ip are basis functions of the univariate polynomial space of degree at most p. Moreover, we have ∂τ = F x ∪ F y . F x contains two faces (located at x = a 1 , b 1 ) and so does F y (located at y = a 2 , b 2 ). We consider the linear mappings

x ∶ [-1, 1] → [a 1 , b 1 ] and y ∶ [-1, 1] → [a 2 , b 2 ].
Let us first consider the two faces in F x . Since we are integrating the partial derivative of v with respect to x, we consider a Gauss-Lobatto quadrature in τ obtained as the tensorproduct of a Gauss-Lobatto quadrature with (p + 1) points in the x variable and a Gauss-Lobatto quadrature with (p + 2) points in the y variable. We use a superscript for the weights and nodes to indicate the number of points in the quadrature, and we set g x j ∶= x(ξ p+1 j ) for all j ∈ I p and g y j ∶= y(ξ p+2 j ) for all j ∈ I p+1 . Using the same arguments as in the proof of Lemma 2, we obtain

∇v ⋅ n τ 2 Fx = b2 a2 (∂ x v x=a1 ) 2 dy + b2 a2 (∂ x v x=b1 ) 2 dy = b 2 -
= p(p + 1) b 1 -a 1 ∂ x v 2 τ .
Similarly, we have

∇v ⋅ n τ 2 Fy ≤ p(p + 1) b 2 -a 2 ∂ y v 2 τ .
Summing the above two inequalities and recalling the definition of h 0 τ gives

∇v ⋅ n τ 2 ∂τ ≤ p(p + 1) b 1 -a 1 ∂ x v 2 τ + p(p + 1) b 2 -a 2 ∂ y v 2 τ ≤ p(p + 1) h 0 τ ∇v 2 τ . (5.5)
Taking square roots completes the proof for d = 2. Finally, the constant is sharp since the upper bound in (5.4) can be attained by univariate functions.

Finally, we consider the case of a simplex.

∇v ⋅ n τ ∂τ ≤ p(p + d -1)(h 0 τ ) -1 2 ∇v τ , ∀v ∈ P p (τ ). (5.6)
Proof. Let v ∈ P p (τ ), let F be a face of τ and set n F ∶= n τ F . Then w ∶= ∇v⋅n F ∈ P p-1 (τ ). Applying the discrete trace inequality from [START_REF] Warburton | On the constants in hp-finite element trace inverse inequalities[END_REF] yields

w 2 L 2 (F ) ≤ p(p + d -1) F d τ w 2 L 2 (τ ) .
Since ∇v ⋅ n F ≤ ∇v 2 (the Euclidean norm of ∇v), we infer that

∇v ⋅ n F 2 L 2 (F ) ≤ p(p + d -1) F d τ ∇v 2 L 2 (τ ) .
Summing over the faces of τ , taking square roots, and recalling the definition of h 0 τ conclude the proof. Remark 5 (Lemma 2). Using the Gauss-Lobatto nodes and their tensor-products to prove discrete trace inequalities is a known technique. The result of Lemma 2 however slightly differs from previous results from the literature and provides a sharper constant. For instance, for p ≥ 1, d = 1 and k = 0, [37, Thm. 2] leads to the constant 2(p + 1) 2 and [27, Lemma 3.1] to the constant (p + 1)(2p + 1), which are both less sharp than C 1 (0, p) in (5.2). Notice also that (5.6) with

d = 1 leads to v ′ ∂τ ≤ 2p 2 h -1 2 τ
v ′ τ which is again less sharp that (5.3) for p ≥ 2. Finally, we have the following multidimensional extension of Lemma 2 in a cuboid; the proof is omitted for brevity and follows arguments similar to those above

. Let τ ∶= [a 1 , b 1 ] × . . . × [a d , b d ] ⊂ R d ,
with b j > a j for all j ∈ {1, . . . , d}, be a cuboid with boundary ∂τ . Let p ≥ 1. For any multi-index (k) = (k 1 , . . . , k d ) with k j ∈ I p for all j ∈ {1, . . . , d}, denoting the k-th partial derivative of v as v (k) , the following holds:

v (k) ∂τ ≤ C d (k, p, τ ) v (k) τ , ∀v ∈ Q p (τ ), (5.7)
with C d (k, p, τ ) ∶= ∑ j∈{1,...,d}

(p-kj +1)(p-kj +2) bj -aj . Moreover, the constant C d (k, p, τ )
is sharp.

Coercivity proof

We can now give the proof of Theorem 1.

Proof of Theorem 1. (i) Tensor-product meshes. For all F ∈ F i h , let T F be the set collecting the two mesh elements sharing F . For all w h ∈ V h p , we have s(w h , w h ) =

F ∈F i h κ F h F ⟦∇w h ⋅ n⟧ 2 F ≤ 2 F ∈F i h τ ∈T F κ F h F ∇w h τ ⋅ n τ 2 F .
Since h F = min τ ∈T F h 0 τ and κ F = min τ ∈T F κ τ (see (2.8)) and exchanging the order of the two summations, we infer that

s(w h , w h ) ≤ 2 τ ∈T h κ τ h 0 τ ∇w h ⋅ n τ 2 ∂τ .
Applying Lemma 3 yields s(w h , w h ) ≤ 2p(p + 1)

τ ∈T h κ τ ∇w h 2 τ ≤ 2p(p + 1)a(w h , w h ).

Recalling that â(⋅, ⋅) = a(⋅, ⋅) -ηs(⋅, ⋅) with η > 0, we conclude that â(w h , w h ) ≥ (1 -2p(p + 1)η)a(w h , w h ).

(5.8)

(ii) Simplicial meshes. The proof is similar but we now invoke Lemma 4 instead of Lemma 3.

Concluding remarks

In this work, we have shown by mathematical analysis and numerical experiments the benefits of tempering the stiffness of the Galerkin FEM approximation of second-order elliptic spectral problems. The idea is to subtract a least-squares penalty on the gradient jumps across the mesh interfaces from the stiffness bilinear form. This novel approximation technique has been named softFEM since it reduces the stiffness of the problem. SoftFEM is formulated in terms of one softness parameter for which we provided an admissible range of values to maintain coercivity on both tensor-product and simplicial meshes. We also gave a practical choice of the softness parameter that leads to superconvergence for linear softFEM in 1D and to attractive numerical performances in more general situations. The main feature of softFEM is that it preserves the optimal accuracy of the eigenvalues in the low-frequency region, while at the same time improving significantly the accuracy in the high-frequency region. The main explanation for this improvement is, as illustrated numerically in our experiments, that in the high-frequency region the standard Galerkin FEM approximation tends to store a substantial amount of energy for the eigenfunctions in the form of gradient jumps across the mesh interfaces. Another very important advantage of softFEM that we illustrated in several settings is its ability to offer a sizable reduction of the conditioning of the stiffness matrix. The optimal value of the asymptotic stiffness reduction ratio increases linearly with the polynomial degree and fairly close values to those predicted theoretically are recovered in our various numerical experiments. As for future work, a first possible direction is the generalization to other differential operators, such as the biharmonic operator. Two possible approaches for the FEM spectral approximation are the mixed (see, e.g., [START_REF] Andreev | Postprocessing and higher order convergence of the mixed finite element approximations of biharmonic eigenvalue problems[END_REF][START_REF] Deng | Optimal spectral approximation of 2n-order differential operators by mixed isogeometric analysis[END_REF]) and the primal (see, e.g., [START_REF] Brenner | C 0 interior penalty Galerkin method for biharmonic eigenvalue problems[END_REF]) formulations. Figure 14 shows the FEM and softFEM spectral approximations of the 1D biharmonic eigenvalue problem: Find an eigenpair (λ, u) such that ∆ 2 u = u (4) = λu in Ω = (0, 1) with the simply supported plate boundary conditions u = u ′′ = 0 on ∂Ω. For the mixed formulation, we decompose u (4) = λu as v ′′ = λu and u ′′ -v = 0. We then apply FEM and soft-FEM, as developed in Section 2, to the decomposed problem with C 0 quadratic elements (notice that softFEM is employed for both equations). The left plot in Figure 14 shows that softFEM maintains the same advantageous features of softFEM as for the second-order operator. In particular, softFEM reduces significantly the high-frequency spectral errors. For the primal formulation, we consider the bilinear form â(v, w) = (v ′′ , w ′′ ) Ω -ηs(v, w) with the softness bilinear form s(v, w) ∶= ∑ F ∈F i h h F (⟦v ′′ ⟧, ⟦w ′′ ⟧) F . The right plot of Figure 14 shows the comparison of FEM and softFEM with various softness parameters when using C 1 cubic splines. The optimal choice for the softness parameter is η = 1 24 , that is, η = 1 2p(p+1) with p = 3. Notice that this optimal value is different from the one found for the second-order elliptic operator ( 1 2(p+1)(p+2) ). Further analysis of the optimality parameter along with the error analysis is postponed to future work.

Another future work direction is the generalization to other discretization methods, such as isogeometric analysis (leading to softIGA) and discontinuous Galerkin methods. Preliminary numerical tests indicate that softIGA has the same features as softFEM: it reduces the stiffness and condition numbers and it improves high-frequency spectral accuracy. More details will be reported in future work. Finally, the stiffness reduction by softFEM lends itself naturally to tempering the CFL condition in explicit time-marching schemes applied to time-dependent PDEs.

Figure 1 :

 1 Figure1: Relative eigenvalue (left) and eigenfunction (right) errors for the 1D Laplace eigenvalue problem when using Galerkin FEM and softFEM with N h = 100 uniform elements and polynomial degrees p ∈ {1, 2, 3}. Upper row: p = 1; middle row: p = 2; bottom row: p = 3. The eigenfunction errors for linear Galerkin FEM and linear softFEM are the same as both discretization methods give the same eigenvectors (but not the same eigenvalues).

. 5 )

 5 Proof. The result (3.4) is well-known; see, for example, [9, Sec. 2] or [33, Sec. 4], whereas the result (3.5) follows for instance from an application of [34, Thm. 2.1].

Figure 2 :

 2 Figure 2: Quadratic softFEM spectra in 1D with N h = 1000 elements using various softness parameters η. Left: â = a + ηs; Right: â = aηs.

2 ,Figure 4 :

 24 Figure 4: Quadratic FEM approximate spectrum with N h = 40×40 elements for the 2D Laplace eigenvalue problem. Left: eigenvalues sorted in each dimension. Right: eigenvalues sorted in 2D.

Figure 5 :Figure 6 :

 56 Figure 5: Relative eigenvalue errors for the 2D Laplace eigenvalue problem when using quadratic Galerkin FEM and softFEM with 40 2 elements.

Figure 7 :

 7 Figure 7: Relative eigenvalue errors for the 3D Laplace eigenvalue problem when using FEM and softFEM with p = 2, 3, 4.

Figure 8 :

 8 Figure 8: Ratio ηs(û h j , ûh j ) a(û h j , ûh j ) for the softFEM eigenfunctions for the Laplace eigenvalue problem in 2D (left) and 3D (right).

Figure 9 :

 9 Figure 9: Relative eigenvalue errors for the elliptic eigenvalue problem (2.1) in 1D with κ(x) ∶= e x sin(2πx) when using Galerkin FEM and softFEM with p ∈ {2, 3, 4, 5}. The mesh has N h = 200 uniform elements.

-Figure 10 :

 10 Figure10: Relative eigenvalue errors for the 1D Laplace eigenvalue problem when using Galerkin FEM and softFEM with p ∈ {2, 3, 4, 5} on a non-uniform mesh composed of N h = 10 elements. The mesh nodes have been randomly set to {0, 0.1, 0.18, 0.29, 0.41, 0.5, 0.59, 0.66, 0.81, 0.92, 1}.

Lemma 2 (

 2 Discrete trace inequality, 1D). Let τ ∶= [a, b] with b > a and ∂τ = 9653 1.2631e3 8.0985e2 1.2675e2 8.1267e1 1.5597 35.88% 2 9.8698 7.2767e3 3.2585e3 7.3727e2 3.3014e2 2.2332 55.22% 3 9.8696 2.1782e4 7.6596e3 2.2070e3 7.7608e2 2.8438 64.84% 4 9.8696 5.0056e4 1.5948e4 5.0717e3 1.6159e3 3.1387 68.14% 5 9.8696 9.9119e4 2.9618e4 1.0043e4 3.0009e3 3.3466 70.12% Table4: Minimal and maximal eigenvalues, condition numbers, stiffness reduction ratios, and percentages for the 1D Laplace eigenvalue problem when using Galerkin FEM and softFEM on a non-uniform mesh composed of N h = 10 elements. The mesh nodes have been randomly set to {0, 0.1, 0.18, 0.29, 0.41, 0.5, 0.59, 0.66, 0.81, 0.92, 1}.

Figure 11 :

 11 Figure 11: Unstructured meshes for the unit square domain (left) and the Lshaped domain (right).

Figure 12 :Figure 13 :

 1213 Figure12: Relative eigenvalue errors for the 2D Laplace eigenvalue problem on the unit square domain when using Galerkin FEM and softFEM with p ∈ {1, 2, 3} and an unstructured mesh.

2 2v 2

 22 at most 2p, it is integrated exactly by the Gauss-Lobatto quadrature with (p+2) points. Considering the linear mapping from ξ ∈ [-1, 1] to [a, b] with x(ξ) ∶= b-a 2 ξ + a+b 2 and setting g j ∶= x(ξ j ) for all j ∈ I p+1 , we have (a) + 2v 2 (b) (p + 1)(p + 2)

Figure 14 :

 14 Figure14: FEM and softFEM spectral approximations of the 1D biharmonic eigenvalue problem using mixed (left) and primal (right) formulations on a uniform mesh composed of 100 elements.

Table 1 :

 1 Errors and convergence rates for the first and sixth eigenpairs using softFEM and polynomial degrees p ∈ {1, . . . , 4}.

			6.54e-5	3.58e-1	5.85e-3	2.10e-2	1.40e1	3.56e-1
		16	4.12e-6	1.78e-1	1.44e-3	4.80e-3	6.63	6.06e-2
	1	32	2.58e-7	8.91e-2	3.60e-4	3.27e-4	3.23	1.35e-2
		64	1.61e-8	4.45e-2	8.98e-5	2.08e-5	1.61	3.27e-3
		rate	4.00	1.00	2.01	3.38	1.04	2.25
		4	4.38e-4	7.57e-2	2.54e-3	3.08e-2	1.37e1	2.82e-1
		8	3.15e-5	1.84e-2	3.40e-4	1.11e-2	3.95	4.47e-2
	2	16	2.04e-6	4.53e-3	4.33e-5	1.80e-3	1.04	7.78e-3
		32	1.29e-7	1.13e-3	5.43e-6	1.50e-4	2.52e-1	1.11e-3
		64	8.06e-9	2.82e-4	6.80e-7	1.02e-5	6.15e-2	1.45e-4
		rate	3.94	2.02	2.97	2.93	1.96	2.72
		4	1.16e-7	5.82e-3	8.08e-5	4.32e-2	5.24	1.04e-1
		8	4.47e-10	7.19e-4	4.80e-6	7.64e-4	9.12e-1	9.29e-3
	3	16	2.08e-12	8.96e-5	2.96e-7	3.02e-6	1.20e-1	4.41e-4
		32	4.04e-13	1.12e-5	1.85e-8	1.15e-8	1.46e-2	2.48e-5
		rate	6.21	3.01	4.03	7.35	2.84	4.05
		4	4.55e-9	2.71e-4	4.54e-6	2.29e-4	2.12	2.39e-2
	4	8	2.09e-11	1.55e-5	1.47e-7	6.70e-6	1.38e-1	7.88e-4
		16	1.25e-13	9.38e-7	4.65e-9	9.01e-8	8.72e-3	3.24e-5
		rate	7.58	4.09	4.97	5.65	3.96	4.77

Table 2 :

 2

			0.8
				linear softFEM
				quadratic softFEM
				cubic softFEM
			0.6	quartic softFEM
			0.4
			0.2
			0 0	100	j	200	250
	Figure 3: Ratio ηs(û h j , ûh j ) a(û h j , ûh j ) for softFEM eigenfunctions. The mesh is
	composed of 240, 120, 80, 60 uniform elements for p ∈ {1, . . . , 4}, respectively.
	p	λ h min	λ h max	λh max	σ	σ	ρ
	1 9.8698 4.7991e5 3.1995e5 4.8624e4 3.2417e4 1.5000 33.33%
	2 9.8696 2.3998e6 1.2000e6 2.4315e5 1.2158e5 1.9999 50.00%
	3 9.8696 6.8046e6 2.7255e6 6.8945e5 2.7615e5 2.4967 59.95%
	4 9.8696 1.5209e7 5.1587e6 1.5410e6 5.2269e5 2.9482 66.08%
	5 9.8696 2.9555e7 9.1006e6 2.9946e6 9.2208e5 3.2476 69.21%

Minimal and maximal eigenvalues, condition numbers, stiffness reduction ratios, and percentages when using Galerkin FEM and softFEM for a mesh composed of N h = 200 uniform elements and polynomial degrees p ∈ {1, . . . , 5}. softFEM for a mesh composed of N h = 200 uniform elements and polynomial degrees p ∈ {1, . . . , 5}. (Recall that λh min ≈ λ h min so that we only show λ h min in the table

Table 3 :

 3 Minimal and maximal eigenvalues, condition numbers, stiffness reduction ratios, and percentages for the 1D elliptic eigenvalue problem with κ(x) ∶= e x sin(2πx) when using Galerkin FEM and softFEM on a uniform mesh composed of N h = 200 elements.

	The weights are

Table 5 :

 5 2 + v(b)2 . Moreover, since v 2 is a polynomial of degree Minimal and maximal eigenvalues, condition numbers, stiffness reduction ratios, and percentages for the 2D Laplace eigenvalue problem on the unit square domain when using Galerkin FEM and softFEM with p ∈ {1, 2, 3} and an unstructured mesh.

	p	λ h min	λ h max	λh max	σ	σ	ρ
	1 2.0020e1 2.9992e3 9.8013e2 1.4981e2 4.8957e1 3.0600 67.32%
	2 1.9740e1 1.5224e4 4.1819e3 7.7122e2 2.1185e2 3.6404 72.53%
	3 1.9739e1 4.0719e4 1.2356e4 2.0628e3 6.2598e2 3.2954 69.65%
	p	λ h min	λ h max	λh max	σ	σ	ρ
	1 4.0162e1 4.7287e3 1.9228e3 1.1774e2 4.7875e1 2.4593 59.34%
	2 3.8707e1 2.5394e4 9.2240e3 6.5605e2 2.3830e2 2.7530 63.68%
	3 3.8619e1 7.1172e4 2.7611e4 1.8429e3 7.1496e2 2.5777 61.21%

Table 6 :

 6 Minimal and maximal eigenvalues, condition numbers, stiffness reduction ratios, and percentages for the 2D Laplace eigenvalue problem on the L-shaped domain when using Galerkin FEM and softFEM with p ∈ {1, 2, 3} and an unstructured mesh.

  Lemma 4 (Discrete trace inequality, simplex). Let τ be a simplex in R d , d ≥ 2, with boundary ∂τ and outward normal n τ . Recall that h 0 τ ∶= d τ ∂τ . Let p ≥ 1. The following holds: