{'fraction_non_alphanumeric': 0.04658560513720011, 'fraction_numerical': 0.013824062788001308, 'mean_word_length': 4.3546641197071, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The virtual element method for eigenvalue problems with potential terms on polytopal meshes

We extend the conforming virtual element method to the numerical resolution of eigenvalue problems with potential terms on a polytopal mesh. An important application is that of the Schrodinger equation with a pseudopotential term. This model is a fundamental element in the numerical resolution of more complex problems from the Density Functional Theory. The VEM is based on the construction of the discrete bilinear forms of the variational formulation through certain polynomial projection operators that are directly computable from the degrees of freedom. The method shows a great flexibility with respect to the meshes and provide a correct spectral approximation with optimal convergence rates. This point is discussed from both the theoretical and the numerical viewpoint. The performance of the method is numerically investigated by solving the Quantum Harmonic Oscillator problem with the harmonic potential and a singular eigenvalue problem with zero potential for the first eigenvalues.


Introduction
The numerical treatment of the Schrödinger equation with local pseudopotentials is one of the most expensive step in solving electronic-structures in large-scale Density Functional calculations, e.g., see [10,60,43]. These calculations make it possible to determine properties of materials from quantummechanical first principles (ab initio), hence without the need of adaptable parameters. A widely used approach for solving the Schrödinger equation in large-scale quantum-mechanical physical systems is provided by the the plane wave (PW) pseudopotential method [55] and its many variants. PW methods are spectral methods based on an expansion on Fourier basis functions (the plane waves). Such methods are generally accurate, but their computer implementation may be inefficient as it normally relies on the fast Fourier transform (FFT), whose nonlocal communication pattern compromises the method's scalability on parallel architectures. Moreover, the strictly uniform resolution of a plane waves expansion makes resolution adaptation infeasible, thus requiring to consider a big number of PWs to capture the highly oscillatory behavior in the atomic region. Such a high resolution is unnecessary in transition regions between atoms where the solution to the Schrödinger equation is normally much smoother. This fact may eventually lead to computationally inefficient discretizations [53].
An alternative to the PW approach has been offered in recent years by real-space methods such as finite differences, finite volumes, and finite elements. Such methods are based on the direct approximation of the solution of the Schrödinger equation on a computational grid. In particular, the finite element method (FEM) is a variational method based on the expansion of the solution in shape basis functions, usually piecewise polynomials that are strictly locally defined in each mesh element. As noted in [51,52], the strictly local nature of the shape functions has several important consequences. First, the FEM produces sparse matrices that can efficiently be treated by standard iterative methods (preconditioned Krylov schemes) [29]; the computational mesh can be refined near the atom locations where the eigenfunctions are expected to vary the most in order to increase the efficiency of the representation; highly scalable implementations are possible on parallel machines. Moreover, the accuracy of the method can easily be improved by increasing the polynomial degree of the shape functions and systematically enhanced by adding other nonpolynomial functions that incorporates in the local approximation some physical insight from the eigenfunction behavior [56,54].
A very recent and important development in the field of the FEM consists in the virtual element method (VEM), which was introduced in [11] as a variational reformulation of the Mimetic Finite Difference (MFD) method [14,45,13]. The VEM is a very successful approach for the construction of numerical approximation of any order of accuracy and regularity on general polygonal/polyhedral meshes. Despite its relative recentness (the first paper was published in 2013), the VEM has been developed successfully for a large range of mathematical and engineering problems [35,22,4,15,32,20,5,25,59,36,57,12,7,58], extended to curved edges [24], and three dimensional problems [3,21,37], using also mixed spaces [28] and nonconforming spaces [8,46,6,31,33]. High-order and higher-order continuity schemes have been presented in [35] and [27,23,5], respectively.
The VEM is indeed a finite element method, so all good properties of the FEM when applied to the Schrodinger equation still hold. In addition, we can exploit the great flexibility of the method, which comes from the fact that the shape functions used in the variational formulation are not known in a closed form, but are defined as the solution of suitable differential problems. This fact is also the motivation of the name "virtual".
The construction of the method and its practical implementations relies on the special choice of the degrees of freedom rather than the explicit knowledge of the local shape functions. The degrees of freedom allow the calculation of certain projection operators from local virtual element spaces into polynomial subspaces. Using such operators, we can properly construct the discrete bilinear forms that approximate the continuous ones of the variational formulation in the virtual element framework.
The present work is the first instance of a long term project that aims to extend the virtual element approach to the real-space numerical approximation of the equations of the Density Functional Theory. We start here by considering the Schrödinger equation with local pseudopotentials and Dirichlet/Neumann boundary conditions. Despite its simplicity, this model allows us to compute the spectrum of the classical Quantum Harmonic Oscillator. We emphasize that the zero potential case, which corresponds to the standard eigenvalue problem for the Laplace operator, also provides the classical "atom in a box" model. Previous works investigating the VEM for eigenvalue problems regard the approximation of the Steklov eigenvalue problem [47,48], the Laplace eigenvalue problem [41,40] with conforming and nonconforming virtual elements, respectively, the acoustic vibration problem [16], the vibration problem of Kirchhoff plates [49], the transmission eigenvalue problem [50] whereas [30] deals with the Mimetic Finite Difference approximation of the eigenvalue problem in mixed form.
The outline of the paper is as follows. In Section 2 we recall the eigenvalue problem under investigation, introducing the classical variational formulation and the necessary notations. Section 3 details the proposed discretization procedure. The approximation spaces and all the bilinear forms that define the discrete problem, are introduced and described. Section 4 deals with the theoretical analysis, which leads to the optimal error estimates of Theorems 4.5, 4.6, and 4.7. In Section 5 we present several numerical tests, which highlight the actual performance of our approach. Finally, in Section 6 we offer our final remarks and conclusions.
2 The continuous eigenvalue problem

Technicalities and definitions
We use the standard definition and notation of Sobolev spaces, norms and seminorms as given in [1]. Hence, the Sobolev space H s (ω) consists of functions defined on the open bounded connected subset ω of R d , d = 1, 2, 3, that are square integrable and whose weak derivatives up to order s are also square integrable. As usual, if s = 0, we prefer the notation L 2 (ω). Norm and seminorm in H s (ω) are denoted by || · || s,ω and | · | s,ω , respectively, and (·, ·) ω denote the L 2 -inner product. The subscript ω may be omitted when ω is the whole computational domain Ω.
If ≥ 0 is an integer number, P (ω) is the space of polynomials of degree up to defined on ω, with the convention that P −1 (ω) = {0}. The L 2 -orthogonal projection onto the polynomial space P (ω) is denoted by Π 0,ω : L 2 (ω) → P (ω). Space P (ω) is the span of the finite set of scaled monomials of degree up to , that are given by where • x ω denotes the center of gravity of ω and h ω its characteristic length, as, for instance, the edge length, the face diameter, or the cell diameter for d = 1, 2, 3; We will also use the set of scaled monomials of degree exactly equal to , denoted by M * (ω) and obtained by setting |α| = in the definition above. Finally, we use the letter C in the estimates throughout the paper to denote a strictly positive constant that is independent of the mesh size h, but may depend on the problem constants, like the coercivity and continuity constants, or other discretization constants like the mesh regularity constant, the stability constants, etc. Constant C generally has a different value at each occurrence.

The continuous model
Let Ω ⊆ R d for d = 2, 3 be the computational domain and let Γ be the boundary of Ω. We are interested in the numerical approximation of the eigenvalues ε ∈ R and the eigenfunctions u : Ω → R, u = 0, solving the following problem in strong form: where V (x) is a smooth real-valued scalar potential function. In the context of atomic and molecular quantum theory, ε and u are the energy level and the corresponding wavefunction of the Hamilton operator H := − 1 2 ∆ + V (x), and potential V (x) collects all local and nonlocal terms from the Density Functional Theory [10,60,43]. For a proper mathematical formulation, problem (2) is supplemented with suitable boundary conditions that, depending on the problem, can be of Dirichlet, Neumann, or periodic type (if Ω is a parallelepiped). In the following we consider for sake of simplicity homogeneous Dirichlet boundary conditions. The other cases easily follow the same construction. Furthermore we assume that V (x) is uniformly bounded from below and above, i.e., there exist two strictly positive constant V * and V * such that V * ≤ V (x) ≤ V * for almost every x ∈ Ω.
The variational formulation of (2) reads as: Find ε ∈ R and u ∈ H 1 0 (Ω), ||u|| L 2 (Ω) = 1, such that where the bilinear form a : and the bilinear form b : Remark 2.1. From the standard eigenvalue theory, see [26,9], we know that (i) problem (3) admits a discrete infinite set of eigenvalues forming a positive increasing divergent sequence; (ii) the corresponding eigenfunctions are an orthonormal basis of H 1 0 (Ω) with respect to the L 2 -inner product and the scalar product associated with the bilinear form a(·, ·); (iii) the eigenvalues may have multiplicity bigger than one, but in such a case the corresponding eigenspace must have finite dimension.
We also consider the source problem with homogeneous Dirichlet boundary conditions: where we assume that f ∈ L 2 (Ω). Well-posedness of problem (6), i.e., existence and uniqueness of its solution, is proved by using the Lax-Milgram Lemma [39] since, due to the boundedness assumption on the potential field V , the bilinear form a in (4) is coercive and the bilinear form b in (5) is continue. Moreover, due to regularity result [2,42], there exists r > 0, depending only on Ω, such that u s ∈ H 1+r (Ω). Eventually, the following stability estimate holds 3 The virtual element method We are interested in developing the virtual element approximation of the eigenvalue problem in variational form (3). To this end, we first discuss which meshes can be used for the numerical formulation and introduce a proper set of regularity assumptions. Then, we define the local and global virtual element spaces, the degrees of freedom and the bilinear forms a h and b h that approximate a and b. Finally, we review the estimate of the convergence rate for the related VEM approximation of the source problem.

Mesh definition and regularity assumptions
Let Ω h denote a decomposition of the computational domain Ω into a finite set of polytopal elements P. As usual, the subindex h that labels the mesh Ω h is the maximum of the diameters h P = sup x,y∈P |x − y| of the elements of the mesh. We assume that the elements are nonoverlapping and for each element P we denote its (d − 1)-dimensional nonintersecting boundary by ∂P; its center of gravity by x P ; its d-dimensional measure by |P|. The boundary of P is formed by straight edges when d = 2 and flat faces when d = 3. On 3D polyhedral meshes, we denote the midpoint and length of each mesh edge e by x e and h e , respectively, and the center of gravity, diameter and area of each face f are denoted by x f , h f , and |f |, respectively. In the 2D case, we do not make any special distinction between the terms "edge" and "face", which we consider as synonyms. To unify the notation we may use the symbol σ instead of e or f and, for example, refer to the geometric objects forming the elemental boundary ∂P by the term side instead of edge/face. According to such notation, we denote the center of gravity, diameter, and measure of side σ by x σ , h σ , and |σ|, respectively.
Consider the set T = {Ω h } h formed by the decompositions of Ω for h → 0. The convergence analysis of the conforming VEM we want to consider in this work requires some regularity assumptions that must be satisfied by all the members of mesh family T = {Ω h } h . For completeness we state these assumptions for both d = 2 and d = 3 case, although those for d = 2 can be derived from those for d = 3 by reducing the spatial dimension.
• d = 3. There exists a positive constant independent of h (and, hence, of Ω h ) such that for every polyhedral element P ∈ Ω h it holds that (i) P is star-shaped with respect to a ball with radius ≥ h P ; (ii) every face f ∈ P is star-shaped with respect to a disk with radius ≥ h f ; (iii) for every edge e ∈ ∂f of every face f ∈ ∂P it holds that h e ≥ h f ≥ 2 h P ; There exists a positive constant independent of h (and, hence, of Ω h ) such that for every polygonal element P ∈ Ω h it holds that (i) P is star-shaped with respect to a disk with radius ≥ h P ; (ii) for every edge e ∈ ∂P it holds that h e ≥ h P ; The scaling assumption implies that the number of edges and faces in each elemental boundary is uniformly bounded over the whole mesh family {Ω h }. The star-shapedness property implies that elements and faces are simply connected subsets of R d and R d−1 , respectively.

The conforming virtual element space
We construct the local conforming virtual element space by resorting to the so-called enhancement strategy introduced in [3]. The construction of the conforming virtual element space in the multidimensional case for d ≥ 3 is recursive. We discuss here only the more general case for d = 3, while the case for d = 2 follows from a simple dimensional reduction.
To this end, on every polygonal face f of the boundary ∂P and for any integer number k ≥ 1 we first define the finite element space It is worth noting that the space of polynomials of degree up to k defined on f is a subspace of V h k (f ).

Then, we introduce the set of continuous linear functionals from
(D3) the moments of v h of order up to k − 2 on each two-dimensional face f : Finally, we introduce the elliptic projection operator Π ∇,f together with the additional conditions: As proved in [11,34], the polynomial projection The local conforming virtual element space of order k on the polygonal face f is the subspace of Space V h k (f ) has the two important properties that we outline below: (i) it still contains the space of polynomials of degree at most k; (ii) the values provided by the set of continuous linear functionals (D1)-(D3) uniquely determine every function v h of V h k (f ) and can be taken as the degrees of freedom of v h .
Property (i) above is a direct consequence of the space definition, while property (ii) follows from the unisolvency of the degrees of freedom (D1)-(D3) that was proved in [3].
To define the conforming virtual element space on the polyhedral cell P, we first need to introduce the "extended" virtual element space: (D4) the moments of v h of order up to k − 2 on P: Then we introduce the H 1 -seminorm projection operator Π ∇,P coupled with the conditions: The polynomial projection Π ∇,P k v h can be computed in terms of the values of the linear functionals (D1)-(D4). Finally, Π ∇,P k is polynomial-preserving, i.e., Π ∇,P k q = q for every q ∈ P k (P).
We are now ready to introduce the local conforming virtual element space of order k on the polytopal element P, which is the subspace of V h k (P) defined as follow: We recall that, by construction, the local space V h k (P) enjoys the following fundamental properties (see [3,21]): (i) it still contains the space of polynomials of degree at most k; (ii) the values provided by the set of continuous linear functionals (D1)-(D4) uniquely determine every function v h of V h k (P) and can be taken as the degrees of freedom of v h .
(iii) we can define an interpolation operator in V h k (P) with optimal approximation properties so that for every v ∈ H r+1 (P) with 1 ≤ r ≤ k the interpolant v I satisfies the inequality: for some positive constant C independent of h.
As for the 2D case, the Finally, the global conforming virtual element space V h k of order k ≥ 1 subordinate to the mesh Ω h is obtained by gluing together the elemental spaces V h k (P) to form a subspace of the conforming space H 1 (Ω). The formal definition reads as: A set of degrees of freedom for V h k is given by collecting the values from the linear functionals (D1)-(D4) for all the mesh elements. The unisolvence of such degrees of freedom is an immediate consequence of their unisolvence on each local space V h k (P).

The VEM for the eigenvalue problem
The next step in the construction of our method is to define a discrete version of the bilinear forms a and b given in (4) and (5). First of all we split the bilinear form a into the sum of local terms: and we note that for an arbitrary pair (u, v) ∈ V h k × V h k the quantity a(u, v) is not computable. Then, following a standard procedure in the VEM framework, we consider a computable discrete local bilinear form a h (·, ·) given by the sum of elemental contributions where we define S P (·, ·) being the stabilization term that will be discussed in the following. The bilinear form a P h depends on the orthogonal projections Π 0,P k−1 ∇u h and Π 0,P k−1 ∇v h , which are computable from the degrees of freedom of u h and v h , respectively [3]. In fact, starting from the definition of the orthogonal projection, an integration by parts yields: where n P,f denotes the unit outward normal to f . The first integral on the last right-hand side is computable from the degrees of freedom (D4) as it is the moment of u h against a polynomial of degree k − 2 over P. The face integrals above are also computable since and the L 2 -orthogonal projection Π 0,f k−1 u h , as we have seen, is computable exactly using only the degrees of freedom of u h , c.f. Remark 3.1.
The discrete form a P h (·, ·) must satisfy the two fundamental properties: k-consistency: for all v h ∈ V h k and for all q ∈ P k (P) it holds stability: there exists two positive constants α * , α * , independent of h and of P, such that Stability is ensured by adding the bilinear form S P , which can be any symmetric positive definite bilinear form on the element P for which there exist two positive constants c * and c * such that Note that S P (·, ·) must scale like a P (·, ·), namely S P (·, ·) h d−2 P (see also Section 5). Following [41,40], we consider two different discretizations of the eigenvalue problem (3) that are obtained by considering two possible choices for the discretization of the bilinear form b (cf. (5)). We split b into the local contributions In the first choice we consider an approximated bilinear form b h , which satisfies the k-consistency property but not the stability property (extending to b h the definitions above). Therefore we simply take The second possible choice consists in considering a discrete bilinear formb h (·, ·) which, as done for the discrete form a h (·, ·), enjoys both the k-consistency property and the stability property. In particular we defineb whereS P is any positive definite bilinear form on the element P such that there exist two uniform positive constants β * and β * such that Remark 3.2. In analogy with the condition on the form S P (·, ·), we require that the formS P (·, ·) scales like b P (·, ·), that isS P (·, ·) h d .
The resulting virtual element scheme read as: if we adopt the first choice b h for the approximation of b. The second virtual element formulation reads where we consider the stabilized bilinear form b h .
Finally, in what follows, we will also need the discrete source problem corresponding to both discrete formulations (36) and (37), which reads respectively as: and find u s The well-posedness of the discrete formulations (38) and (39)  We finally observe that both bilinear forms are fully computable for any couple of functions (u h , v h ) ∈ V h k , since the enhancement technique implies that Π 0,P k u h (resp. Π 0,P k v h ) can be computed using only the degrees of freedom of u h (resp. v h ).
The following convergence estimate theorem holds for the approximation of the source problem [3].

Spectral approximation for compact operators
In this section, we briefly recall the results of the spectral approximation theory for compact operators. For more general results, we refer to the original papers [9,26,44].
We introduce a natural compact operator associated with problem (3) and its discrete counterpart, and we recall their connection with the eigenmode convergence.
Let T ∈ L(L 2 (Ω)) be the solution operator associated with problem (3). T is the bounded linear operator T : L 2 (Ω) → L 2 (Ω) which maps the forcing term f to u =: T f : Operator T is self-adjoint and positive definite with respect to the bilinear forms a(·, ·) and b(·, ·) on H 1 (Ω), and compact due to the compact embedding of H 1 (Ω) in L 2 (Ω).
Similarly, let T h ∈ L(L 2 (Ω)) and T h ∈ L(L 2 (Ω)) be the discrete solution operators associated with the stabilized and non-stabilized discrete source problems. The former is the bounded linear operator mapping the forcing term f to u h =: T h f and satisfies: The latter is the bounded linear operator mapping the forcing term f to u h =: T h f and satisfies: Both operators T h and T h are self-adjoint and positive definite with respect to the discrete bilinear form a h (·, ·), b h (·, ·) and a h (·, ·), b h (·, ·). They are also compact since their ranges are finite dimensional.
The eigensolutions of the continuous problem (3) and the discrete problems (36) and (37) , u) is an eigenpair for the operator T , and analogously for problems (36) and (37) and operators T h and T h . Thanks to this correspondence, the convergence analysis can be derived from the spectral approximation theory for compact operators. In the rest of this section we refer only to operators T and T h . Identical considerations hold for operators T and T h and we omit them for brevity.
A sufficient condition for the correct spectral approximation of a compact operator T is the uniform convergence to T of the family of discrete operators { T h } h (see [26,Proposition 7.4], cf. also [9]): or, equivalently, with ρ(h) tending to zero as h goes to zero. Condition (43) usually follows by a-priori estimates with no additional regularity assumption on f . Besides the convergence of the eigenmodes, condition (42), or the equivalent condition (43), implies that no spurious eigenvalues may pollute the spectrum. In fact, each discrete eigenvalue approximates a continuous eigenvalue and each continuous eigenvalue is approximated by a number of discrete eigenvalues (counted with their multiplicity) that corresponds exactly to its multiplicity. We now report the main results about the spectral approximation for compact operators.   1,h , . . . , µ m,h (repeated according to their respective multiplicities) converges to µ. Moreover, let E µ,h be the direct sum of the eigenspaces corresponding to the discrete eigenvalues µ 1,h , · · · , µ m,h converging to µ. Then with where, in general,δ (U, W ) = sup u∈U,||u|| 0 =1 inf w∈W ||u − w|| 0 denotes the gap between U , W ⊆ L 2 (Ω).
Concerning the eigenvalue approximation error, we recall the following result.
where µ 1,h , . . . , µ m,h are the m discrete eigenvalues converging to µ repeated according to their multiplicities.

Convergence analysis for the stabilized formulation
In this section we study the convergence of the discrete eigenmodes provided by the VEM approximation to the continuous ones. We will consider the stabilized discrete formulation (37). The analysis can be easily applied to the non-stabilized one (36).
Proof. The proof follows the same lines as the one of Theorem 6.4 in [41]. We recall it here for the convenience of the reader. Let u s and u s h be the solutions to the continuous and the discrete source problems (6) and (38), respectively. From the L 2 -estimate of Theorem 3.3 with f ∈ L 2 (Ω) and the stability estimate in (7) we have that ||u s − u s h || 0 ≤ Ch min(t+1,1) ||f || 0 with t = min(k, r), k ≥ 1 being the order of the method and r being the regularity index of the solution u s ∈ H 1+r (Ω) to the continuous source problem. Then it follows that We remark that if f ∈ E µ then, thanks to the L 2 a priori error estimate in Remark 3.4, it holds Putting together Theorem 4.1, Theorem 4.3, and the above observation, we can state the following result.
Theorem 4.4. Let µ be an eigenvalue of T , with multiplicity m, and denote the corresponding eigenspace by E µ . Then, exactly m discrete eigenvalues µ 1,h , . . . , µ m,h (repeated according to their respective multiplicities) converges to µ. Moreover, let E µ,h be the direct sum of the eigenspaces corresponding to the discrete eigenvalues µ 1,h , · · · , µ m,h converging to µ. Then A direct consequence of the previous result (cf. [9,26]) is the following one.
where t = min{k, r}, being k the order of the method and r the regularity index of u.
We now state the usual double order convergence of the eigenvalues. Theorem 4.6. Let ε be an eigenvalue of problem (3) with multiplicity m, and denote by ε 1,h , · · · , ε m,h the m discrete eigenvalues of problem (37) converging towards ε. Then the following optimal double order convergence holds: with t = min{k, r}, being k the order of the method and r the regularity index of the eigenfunction corresponding to ε.
Proof. The proof follows the guidelines of Theorem 6.4 in [41] and Theorems 4.2 and 4.3 in [49]. For an alternative proof see also Theorem 6.6. in [40], taking into account that in our case the term N h (·, ·), relative to the conformity error, vanishes.
Eventually, we state the optimal error estimate for the eigenfunctions in the energy norm.
Theorem 4.7. With the same notation as in Theorem 4.5, we have where t = min(k, r), k being the order of the method and r the regularity index of u ∈ H 1+r (Ω).
Proof. The proof of this result is similar to the one for the finite element method. We briefly report it here for the sake of completeness. It holds that The first term at the right-hand side of the previous equation is of order h 2t , while the second one is of order h t . Finally, for the last term, using (31), the continuity of the operator T h , and Theorem 4.5, we obtain

Numerical experiments
In this section, we investigate the behavior of our virtual element method for the numerical treatment of the eigenvalue problem (3). In particular, we present the performance of the conforming VEM applied to the eigenvalue problem on a two-dimensional square domains. We use the "diagonal" stabilization [18] for the bilinear form a P h (·, ·) (cf. (27)) and b P h (·, ·) (cf. (35), which reads as follows: where v h , w h denote the vectors containing the values of the local degrees of freedom associated to v h , w h ∈ V h k (P) and the parameters σ P and τ P are two positive h-independent constants. In the numerical tests we choose σ P as the mean value of the eigenvalues of the matrix stemming from the consistency term (Π 0,P k−1 ∇·, Π 0,P k−1 ∇·) P for the grad-grad form (see (27)). In the same way we pick τ P as the mean value of the eigenvalues of the matrix resulting from the term 1 h 2 P (Π 0,P k ·, Π 0,P k ·) P for the mass matrix (see (35)).
The convergence of the numerical approximation is shown through the relative error quantity where ε denotes an eigenvalue of the continuous problem and ε h its virtual element approximation.

Test 1.
In this test case, we numerically solve the 2D Quantum Harmonic Oscillator problem that corresponds to the Schrodinger equation with the harmonic potential V (x, y) = (1/2)(x 2 + y 2 ). The eigenvalues are a suitable combinations of the eigenvalues of the one dimensional problem and are given by the natural numbers n = 1, 2, 3, . . ., each one with multiplicity n. The eigenfunctions of such a problem are obtained through the two-dimensional tensor product of one-dimensional Hermite functions, which are given by the Hermite polynomials multiplied by the Gaussian function w(x, y) = exp −(x 2 +y 2 ) . As these eigenfunctions are rapidly decreasing to zero for x, y tending to infinity due to the Gaussian term, we can assume homogeneous Dirichlet boundary conditions if the computational domain is sufficiently large. For such reason, we solve the eigenvalue problem on the square domain Ω =] − 10, 10[×] − 10, 10[. On this domain, we consider four different mesh sequences, hereafter denoted by: • Mesh 1, mainly hexagonal mesh with continuously distorted cells; • Mesh 2, nonconvex octagonal mesh; • Mesh 3, randomized quadrilateral mesh; • Mesh 4, central Voronoi tessellation.
The first mesh of each sequence is shown in Figure 1. These mesh sequences have been widely used in the mimetic finite difference and virtual element literature, and a detailed description of their construction can easily be found elsewhere, for example, see [13].
The convergence curves for the four mesh sequences above are reported in Figures 2, 3, 4, and 5.
The expected rate of convergence is shown in each panel by the triangle closed to the error curve and indicated by an explicit label. For these calculations, we used the VEM approximation based on the conforming virtual element space V h k , k = 1, 2, 3, 4, and the VEM formulation (36) using the nonstabilized bilinear form b h (·, ·). As already observed in [41] for the conforming VEM approximation of the Laplace eigenvalue problem, the same computations using formulation (37) and the stabilized bilinear     b h (·, ·) produce almost identical results, which, for this reason, are not shown here. These plots confirm that the conforming VEM formulations proposed in this work provide a numerical approximation with optimal convergence rate on a set of representative mesh sequences, including deformed and nonconvex cells, of the Schrodinger equation problem, i.e., the standard eigenvalue problem with a regular potential term in the Hamilton operator at left hand-side.

Test 2 (piecewise constant diffusivity tensor)
The present test problem is taken from the benchmark singular solution set in [38]. We here consider the square domain Ω = (−1, 1) 2 split into two subdomains Ω δ and Ω 1 (see the left plot in Figure 6), and we study the eigenvalue problem on the square with discontinuous diffusivity tensor and zero potential V coupled with Neumann homogeneous boundary conditions i.e. we consider the following problem in strong form: −∇ · (K(x)∇u(x)) = εu(x) in Ω and ∂u ∂n = 0 on Γ, Therefore the continuous bilinear form associated to the eigenvalue problem is a P K (u, v) := P K∇u · ∇v dx whose virtual approximation (see [19,17]) is given by to be used in place of a P h (u h , v h ) (cf. (27)) in Problem (37), where K = K ∞,P . We consider K |Ω1 = I and K |Ω δ = δ −1 I with four different values of δ, namely δ = 0.50, 0.10, 0.01, 1e − 8.
We apply the Virtual Element method (37) using a sequence of Voronoi meshes with mesh diameter h = 1/2, 1/4, 1/8, 1/16 (see the right plot in Figure 6 for an example of the adopted meshes). We show the plot of the convergence for the first eight computed eigenvalues in Figures 7 ,8,9,10. We compute the relative errors by comparing our results with the values given in [38].
In accordance with Theorem 4.6, we obtain different rates of convergence that are determined by the polynomial order of the method and by the regularity of the corresponding exact eigenfunctions [38]. Taking this into account, we show that the method is overall optimal, and thus stable with respect to discontinuities in the diffusivity tensor.

Conclusions
We have discussed the application of the conforming virtual element method to the numerical resolution of eigenvalue problems with potential terms on polytopal meshes. The most notable case is that of    the Schrödinger equation with a suitable pseudopotential, which is fundamental in the numerical treatment of more complex problems in the Density Functional Theory. The VEM approximation of such problem was discussed from both the theoretical and the numerical viewpoint, proving that the method provides a correct spectral approximation with optimal rates of convergence. The performance of the method was shown by computing the first eigenvalues of the Quantum Harmonic Oscillator provided by the harmonic potential and a singular eigenvalue problem with zero potential.