Exponential convergence in H1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1$$\end{document} of hp-FEM for Gevrey regularity with isotropic singularities

For functions u∈H1(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u\in H^1(\Omega )$$\end{document} in a bounded polytope Ω⊂Rd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega \subset {\mathbb {R}}^d$$\end{document}d=1,2,3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1,2,3$$\end{document} with plane sides for d=2,3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=2,3$$\end{document} which are Gevrey regular in Ω¯\S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Omega }\backslash {\mathscr {S}}$$\end{document} with point singularities concentrated at a set S⊂Ω¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {S}}\subset \overline{\Omega }$$\end{document} consisting of a finite number of points in Ω¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Omega }$$\end{document}, we prove exponential rates of convergence of hp-version continuous Galerkin finite element methods on affine families of regular, simplicial meshes in Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document}. The simplicial meshes are geometrically refined towards S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {S}}$$\end{document} but are otherwise unstructured.


Introduction
Many nonlinear PDEs admit solutions which are smooth in a bounded, polytopal domain ⊂ R, but exhibit isolated point singularities at a set S ⊂ . We mention only nonlinear Schrödinger equations with self-focusing, density functional models in electron structure calculations (eg. [3,6,7,17] and the references there), nonlinear parabolic PDEs with critical growth (eg. [22,29] and the references there, or continuum models of crystalline solids with isolated point defects (eg. [25] and the references there).
We prove an exponential convergence result for C 0 -conforming hp-FEM on regular, simplicial mesh families with isotropic, geometric refinement towards the singular point(s) c ∈ S . These meshes are in addition required to be shape-regular. This type Research performed in part while the authors were visiting the Erwin Schrödinger Institute (ESI) in Vienna, Austria, during the ESI thematic term "Numerical Analysis of Complex PDE Models in the Sciences" from June-August 2018. CS was supported in part by the Swiss National Science Foundation. Seminar for Applied Mathematics, ETH Zurich HG G57.1, CH8092 Zurich, Switzerland of mesh arises for example in adaptive bisection-tree refinements. Specifically, for singular solutions u ∈ H 1 ( ) in the bounded domain ⊂ R d , d = 2, 3 which belong, in addition, to a countably normed space with non-homogeneous, radial weights as introduced, for example, in [4,13], and with Gevrey-regular growth of derivatives in \S , we construct a sequence {I hp N } N of continuous, piecewise polynomial (quasi-)interpolation operators on sequences of regular, simplicial partitions that are geometrically refined towards S with exponential convergence in H 1 ( ): for a bounded domain ⊂ R d and for functions u ∈ H 1 ( ) ∩ G δ ( ), a class of δ-Gevrey-regular functions in \S (to be defined in (8) below), there exist constants b, C > 0 which depend on and on u, such that for all N Here, d = 2, 3 denotes the space dimension and N denotes the number of degrees of freedom in the hp-FE approximation, (·) denotes the Gamma function and δ > 0 denotes the Gevrey regularity parameter.
Since the pioneering work [4], Gevrey regularity of solutions has been verified in a number of nonlinear partial differential equations; we mention only incompressible Euler [9] and KdV equations [10].
Note also that δ = 1 corresponds to functions which are analytic in \S which case was considered in [2,12,15,16,[19][20][21]. In the presently considered case δ > 1, i.e. smooth but non-analytic functions u are admitted. This precludes the use of analytic continuation into Bernstein Ellipses and of classical complex variable methods of polynomial approximation to establish exponential convergence rate bounds. We therefore opt here for tensorized univariate hp-projectors with analytic error bounds (being explicit in the polynomial degree and in the regularity parameter). Approximations on tetrahedral elements T are obtained by restricting the resulting local hp-interpolants from a corresponding parallelepiped K T containing T . Rendering the tensorized hpinterpolant well-defined in K T requires K T ⊂ . In Lemma 1, we show that this requirement can be achieved in regular, geometric meshes of triangles resp. tetrahedra with sufficient refinement, in any bounded polytope .
The rate (1) coincides, in the cases d = 1, 2 and for analytic solutions, i.e. when δ = 1, with the exponential convergence rate bounds obtained in [18,19] for corner singularities on structured geometric meshes (consisting of axiparallel quadrilaterals with inserted triangles to remove irregular nodes). In space dimension d = 3, (1) generalizes the hp-approximations in [30,Sec. 5.2.2] in the case of vertex singularities, for meshes of axiparallel hexahedra to unstructured, tetrahedral meshes with geometric refinement towards S . For 0 < δ < 1 we obtain super-exponential convergence rate bounds, which benefit from the higher regularity of u.
The structure of the note is as follows: in Sect. 2, we introduce a model problem, the geometric assumptions on the singularities, and precise the analytic regularity in countably normed, weighted Sobolev spaces with radial weight functions. In Sect. 3, we introduce the hp-version FEM; we specify in particular the assumptions on the simplicial, geometric meshes, on the elemental polynomial degrees, and on the def-inition of the hp FE spaces. Sect. 4 contains a proof of the exponential convergence bound in H 1 ( ) on regular, simplicial geometric mesh families.

Analytic regularity
Analytic regularity is characterized in countably normed weighted Sobolev spaces which have been introduced and used in exponential convergence estimates in a number of references; we only mention [2,13,[18][19][20][21] and the references there. Here, we denote by S ⊂ the set of singular points c; we consider solutions u ∈ H 1 ( ) which are smooth in \S so that the singular support of u coincides with S . We work under the following separation assumption on S .
The singular set S consist of a finite number of isolated points c ∈ . (2) , and allows to partition the set into |S | many disjoint neighborhoods ω c of the singularities c ∈ S . We set denote 0 := \ c∈S ω c .

Weighted Sobolev spaces: Gevrey classes
We characterize analytic regularity of singular solutions by weighted Sobolev spaces. To define these, we introduce distance functions: With c ∈ S we collect all singular exponents β c ∈ R in the "multi-exponent" We assume for d = 3 (β > s and β ± s being understood componentwise for For d = 2, we assume for some ε > 0 that We consider the inhomogeneous, weighted semi-norms |u| N k β ( ) given by (cp. [13, Definition 6.2 and Equation (6.9)], [2] and [20]), We define the inhomogeneous weighted norm u N m β ( ) by u 2 N m Here, |u| H m ( 0 ) signifies the Hilbertian Sobolev semi-norm of integer order m on 0 , and D α denotes the weak partial derivative of order α ∈ N d 0 . The space N m β ( ) is the weighted Sobolev space obtained as the closure of C ∞ 0 ( ) with respect to the norm · N m β ( ) .
The spaces N m β ( ) are closely related to the nonhomogeneous, weighted spaces of type J m γ ( ) which arise in connection with the Mellin transformation of elliptic problems in conical domains. We refer to [12] for a definition and properties of the spaces J m γ ( ). (iv) The spaces N m β ( ) are related to the weighted Sobolev spaces introduced by Babuška and Guo [2,20] in space dimension d = 2 and d = 3, resp. For example, in space dimension d = 3, the weighted seminorm of order k ≥ 2 in (7) coincides with with the vertex-weighted seminorm H k,l β m defined in [20, page 83, top] for l = 2, if we note that for l = 2 and |α| = k ≥ 2 in [20], the vertex-weight function α,l β m (x) = |x| β m +|α|−l = |x| β m +k−2 . Comparing with the weighted seminorm in (7) yields β c = −1 − β m so that the condition β m ∈ (0, 1/2) in [20,Thm. 3.5] translates into the condition β c ∈ (−1, −3/2) in item (i) above.
With N k β ( ) as defined in (7), for δ > 0 we define the δ-Gevrey regular class of solutions with point singularities at S by We mention that in the case δ = 1 and space dimension three, the norm in the definition (8) of the Gevrey class G δ β (S ; ) coincides with the norm for the analytic class B β (S ; ) introduced in [13, Definition 6.9-6.11] in three space dimensions, and, by the observation in Remark 1, (iv), and the Definition of the countably normed spaces B l β m in [20, Page 83,top] with l = 2. In two space dimensions, it equals the weighted analytic classes introduced in [19,20]. All ensuing approximation results in particular apply for this analytic solution class, as has been indicated in [31,32]. Naturally, the present construction parallels earlier constructions in particular cases; for example, the polynomial trace lifting in Sect. 4.2.6 is identical to the analytic case in [32].

Examples of boundary value problems with Gevrey-regular solutions
Large classes of linear and nonlinear elliptic boundary value and eigenvalue problems with analytic input data admit solutions in the analytic class G δ β (S ; ) with δ = 1. We refer to, e.g., [4,13,21] and also [17] for electron structure models, [2,13,21] for elliptic problems in polyhedral domains, and [27] and the references there for nonlinear Schrödinger eigenvalue problems.

Linear elliptic boundary value problems in polygons
In space dimension d = 2, let denote a polygon with straight sides. Consider the model Dirichlet boundary value problem In (9), we assume that A(x) = (a i j (x)) 1≤i, j≤2 and f (x) are analytic in and that the matrix function x → A(x) ∈ R 2×2 sym is uniformly positive definite: there exists α > 0 such that for every ξ ∈ R 2 holds ess inf The unique, weak solution u ∈ V = H 1 0 ( ) of (9) exists by the Lax-Milgram Lemma, and satisfies the weak form of (9): find For a closed subspace V N ⊂ V , approximate solutions u N ∈ V N of (10) are obtained by Galerkin projection: find The approximate solutions u N exist, are unique and quasioptimal: Convergence rates of sequences {u N } N of approximate solutions thus depend on (a) the choice of V N and (b) on the solution regularity. For problem (10), it has been shown in [2] that the solution u ∈ G 1 δ ( ). For {V N } N being a sequence of so-called hp-FE spaces (to be defined in the next section), we recover from (12) and (1) (with δ = 1 and d = 2) the exponential convergence rate exp(−b 3 √ N ) already obtained in [18,31]. Gevrey regularity in conical domains for Gevrey-regular data A and f for solutions u of (10) was first obtained in [4].

Three-dimensional problems
For the analog of (9) in polyhedral domains , the regularity classes G δ β ( ) are not adequate, as even for analytic data A and f , the solutions are locally analytic in , but exhibit apart from corner singularities also so-called edge-singularities. Their precise mathematical description mandates more sophisticated function spaces (see, e.g., [13,20] and the references there and [30] for exponential convergence results for hp-FEM. However, in large classes of applications, solutions are Gevrey regular with point singularities only. We mention only the source problem (10) in domains which exhibit isolated vertices, e.g. conical domains with a smooth (analytic) base, such as circular cones with apex c (see, e.g., [23]).
Another important class of problems arises from mathematical models of quantum chemistry (see, e.g., [6,7] and the references there). For instance, consider the nonlinear Schrödinger EVP: find λ ∈ R and 0 = u ∈ H 1 (R 3 ) such that Here, for analytic potentials V which become singular at a finite set S ⊂ R 3 of isolated points, eigenfunctions u belong to G 1 β ( ) for compact sets ⊂ R 3 containing S in their interior, see [17,26] and [27,Thm. 7]. Quasioptimality in H 1 ( ) of Galerkin-FEM for the EVP (13) can be found, for example, in [6,7].

hp-Finite element spaces
The hierarchies of FE spaces which underlie the hp-FEM are based on two key ingredients: (i) geometric mesh families M κ,σ = {M ( ) } ≥0 and (ii) simultaneous refinement of meshes and polynomial degree distributions. They also exhibit (iii) a layer-structure among the Finite Elements T ∈ M ( ) which we describe next.

Geometric mesh families M Ä,
For two parameters 0 < κ, σ < 1, we consider in the bounded polyhedron geometric mesh families M κ,σ = {M ( ) } ≥1 of geometrically refined, regular simplicial triangulations M ( ) ∈ M κ,σ . Here,the parameter denotes roughly the number of refinements towards the singularities. The meshes M ∈ M κ,σ are regular partitions of the polyhedron into a finite number of open simplices (triangles in space dimen- Here, regular means that for every M ∈ M κ,σ , the intersections of closures of any two distinct T , T ∈ M are either empty, a vertex v, an entire edge e or an entire face f . We assume the family M κ,σ to be uniformly κ-shape regular: for a simplex T ∈ M ( ) , we denote by h T = diam(T ) its diameter and by ρ T = sup{ρ > 0|B ρ ⊂ T }, the radius of the largest ball B ρ that can be inscribed into T . For a regular, simplicial mesh M , the (nondimensional) shape parameter κ( For a regular, simplicial triangulation M of with κ(M ) < ∞, the affine element maps are nondegenerate: the jacobians B T = D F T in (14) are nonsingular, and and assume that K T ⊂ .

Lemma 1
In a bounded polyhedron ⊂ R 3 with plane sides, for any given regular partition M of into simplices, there exists k ∈ N depending only on d ∈ {2, 3} such that k fold uniform newest-vertex bisection refinement of M guarantees that for each T ∈ M there exists K T as described above with K T ⊆ .
Proof We refine each simplex T 0 ∈ M exactly k-times using newest-vertex-bisection (NVB) from [34] resulting in the regular simplicial mesh M k in . Since NVB generates only a bounded number of different shapes in a bounded number of spatial orientations per element T 0 ∈ M (depending only on d [35, Section 5]), we may choose k ∈ N such that for each Obviously, we may construct K T 1 ⊆ K T and it suffices to show that K T ⊆ T . We define F T to map the origin to a vertex of T which is closest to a vertex of T . Thus, application of F −1 T simplifies the situation to T = a T + b ⊂ T for b ∈ R d with | b| ≤ 1/2. Let K denote a parallelepiped for which the set of vertices contains all vertices of T . Denote b a vertex closest to the origin. Then, the distance of every vertex of K to the origin is bounded by Reverting the transformation concludes the proof.

Local polynomial spaces
For T ∈ M the local polynomial approximation space P p (T ) = span{x α : |α| ≤ p} is the span of all multivariate polynomials on T ∈ M whose total degree does not exceed p. The space P p (T ) is invariant under the affine mapping F T , i.e. u ∈ P p (T ) if and only ifû : For each parallelepiped K T associated with a d-simplex T ∈ M with associated affine element mapping F T : K → K T and polynomial degree p ≥ 0, we set For polynomial degree p ≥ 1, and for a family of regular, simplicial triangulations M ( ) ∈ M κ,σ of , we introduce finite element spaces of continuous, piecewise polynomial functions of total degree at most p on each T ∈ M ( ) , i.e.
Typically, hp-FEMs are obtained when the level of geometric mesh refinement is tied to the polynomial degree p.

Mesh layers
A key ingredient in exponential convergence proofs of hp-FEM is geometric mesh refinement towards the set S of singularities. For a parameter 0 < σ < 1, we call a regular, simplicial mesh family We tag members of a σ -geometric family M κ,σ by a subscript σ , i.e. we write M ( ) σ . A simplicial mesh family M κ,σ in a polytopal Lipschitz domain which is σ -geometrically refined towards the singular set S can be generated by the following algorithm: such that each c ∈ S corresponds to a vertex of some T ∈ M (0) and such that for each T ∈ M (0) , T ∩ S is either empty or contains exactly one element of S . For = 0, 1, . . . do: We now show that the output of Algorithm 1 can, for each fixed 0 < σ < 1, be identified with sequences M κ,σ of regular, simplicial meshes which are σ -geometrically refined towards S in .
We start by observing that the mesh completion step 2. in Algorithm 1 from [34] guarantees that #(M ( ) ) · #S , where the hidden constant depends only on M (0) and S . We next observe that the sequence M produced by Algorithm 1 does not depend on σ ∈ (0, 1). Therefore, choosing σ ∈ (0, 1) after executing Algorithm 1 we claim that M can be identified with a σ -geometric family M κ,σ in sense that (18) holds.

Proposition 1
For every fixed 0 < σ < 1, the output M of Algorithm 1 can be identified with a σ -geometric family M κ,σ . In particular, for given 0 < σ < 1, for every > d all elements T ∈ M ( ) can be grouped in mesh-layers: there exists a partition where #T ( ) σ ≤ c T (κ, σ ) and for k log(2)/| log(σ )| and there exists c T > 0 being independent of such that for all holds There exists a constant c(M κ,σ ) ≥ 1 with and such that, for every T ∈ L k and every k ≥ 1, and conv x d , , there is exactly one child T which contains the singularity and also c ∈ E T . This shows that once the refinement edge includes the singularity, this remains the case for all descendants of T which include the singularity. Moreover, it is evident that any vertex of T is part of the refinement edge of one of its descendants T after at least d-bisections. Those two observations imply that for > d, the refinement where T ( ) contains all the remaining elements which include a singularity and k is bounded by k log(2)/| log(σ )|. By shape-regularity, we know that each T ∈ M ( ) which does not include a singularity has a similarly sized neighbor between itself and S . Hence Proof Assertion (i) is already stated in Proposition 1 and repeated just for reference. Property (20) implies that for every T ∈ T ( ) σ , dist(T , S ) ≤ c T (κ, σ )σ . This implies, with the shape regularity of T , that for every T ∈ T ( ) σ holds |T | ≤ c T (κ, σ )σ d . This, in turn, implies assertion (ii). With shape-regularity, we derive (iii) directly from (ii).

Remark 2 (i)
We do not use that fact that the singular supports c ∈ S comprise nodes of some triangulation M ( ) ∈ M κ,σ . This implies, in particular, that the ensuing exponential convergence proofs remain valid for "nearly coalescing" singular supports c, c ∈ S : for c, c ∈ S such that dist(c, c ) < σ p , both c and c are contained in the terminal layers T ( ) σ . There, a low-order quasi interpolant of Clément (resp. Scott-Zhang) type is used, see Sect. 4.2.8 ahead. The constants in the exponential convergence bound are uniform in w.r. to dist(c, c ). (ii) Due to Prop. 2, item (iii), geometric mesh refinement implies that dist(c, c ) is resolved with geometric refinements with ≥ O (| log(dist(c, c ))|) many mesh layers. Here, If, moreover, u| ∂ = 0, then ( p κ,σ u)| ∂ = 0 and (23) holds.

Proof
The proof of the approximation result Theorem 1 is based on constructing the projectors p κ,σ ; our construction will proceed in several steps and we detail it for d = 3, the case d = 2 being a (minor) modification. First, we review from [30, Section 5] a family of univariate hp-projections with error bounds which are explicit in the polynomial degree as well as in the regularity of the functions to be approximated. A corresponding family of polynomial projectors on the unit cube K = (0, 1) 3 with analogous consistency error bounds is then obtained as in [30, Section 5] by tensorization and scaling. We shall use these bounds for a tetrahedron T ∈ O ( ) σ ⊂ M ( ) σ ∈ M κ,σ as follows. By Proposition 1, T ∈ L k for some 1 ≤ k ≤ − 1. The (up to orientation) unique parallelepiped K T = F T ( K ) associated with T ∈ L k has the same scaling properties as T , in particular (22) also holds for K T . For u belonging to the Gevrey class (8) with weight vector satisfying (5), the pullback u T = u| K T • F T satisfies on K the same derivative bounds as u| T • F T on T (with possibly larger constant C u , depending on κ, but independent of and of T ). The tensorized hp interpolation operator from [30] on K is therefore well-defined and allows to construct a polynomial approximation u p T ∈ Q p ( K ) with analytic consistency error bounds on K ; since T ⊂ K , and since Q p ( T ) ⊂ P pd ( T ), the pushforwards of the restrictions u p T | T under the affine mapping F T : T → T will be local polynomial approximations of degree pd with exponential convergence estimates in H 1 (T ). Moreover, since the tensorized interpolant is nodally exact in the vertices of K , and since the set of vertices of T is a subset of the set of vertices of K , the pushforwards of u p T | T under F T are nodally exact in the vertices of T . For elements T ∈ T ( ) σ , we only require a first order approximation property, as the geometric refinement guarantees the necessary convergence rate. We can not use nodal interpolation as functions u ∈ G δ β (S ; ) may not be bounded near a singularity c ∈ S . Thus, we construct a quasi interpolation operator on elements in the terminal layers T ∈ T ( ) σ that interpolates at those vertices of T which are not in S . By the continuity of u ∈ G δ β (S ; ) on \S , the resulting global, piecewise polynomial approximation is nodally exact in all vertices of M ( ) σ except those which coincide with singularities c ∈ S . Particularly, the resulting piecewise polynomial hp-approximation is globally continuous at all vertices of M ( ) σ . However, it still has polynomial jump discontinuities across edges and (in space dimension d = 3) faces of T ∈ M ( ) σ which we remove by polynomial trace liftings, preserving the exponential convergence estimates.

Univariate hp-projectors and hp error bounds
Let I = (−1, 1) be the unit interval. For any k ≥ 1, we write H k (I ) for the usual Sobolev space endowed with norm u H k (I ) . For q ≥ 0, we denote by π q,0 : L 2 (I ) → P q (I ) the L 2 (I )-projection. The following C k−1 -conforming and univariate projector has been constructed in [14,Section 8].
Moreover, there holds: (i) For every k ∈ N, there exists a constant C k > 0 such that (ii) For integers p, k ∈ N with p ≥ 2k − 1, κ = p − k + 1 and for u ∈ H k+s (I ) with any k ≤ s ≤ κ there holds the error bound We refer to [14,Proposition 8.4] and [14,Theorem 8.3], respectively, for proofs, and further references.

Tensor projector on the unit cube
Based on the univariate projectors π p,k , we constructed in [30] polynomial projection operators on I d = (0, 1) d by (a) translation and scaling of the projectors π p,k to (0, 1) and (b) by tensorization, as follows: for integers k ≥ 0 and d > 1, we define where ⊗ denotes the tensor-product of separable Hilbert spaces.
where π (i) p,k denotes the univariate projector in Lemma 2, applied in coordinate 1 ≤ i ≤ d. For d, k ≥ 1 there exists a constant C k,d > 0 such that for all p ≥ 2k − 1 there holds the stability bound and We choose throughout what follows k = 2 as in [30], and obtain from (29), (25) Proposition 3 [30] Assume that the polynomial degree p ≥ 5. Then, for any integers where the constant implied in is independent of s and of p, and where q,r = 2 2(r +3) (q + 1 − r ) Moreover, 3 p,2 v is nodally exact in the vertices of K = (0, 1) 3 :

Transformation formula
For u ∈ H k ( ), and for a simplex T ∈ O ( ) σ , consider the transformationû T = u| T • F T ∈ H k ( T ) for every k ≥ 0. Quantitative bounds on derivatives under affine transformations F T in (14) are provided by the transformation formula (eg. [5, Section II.6.6]).

Element interpolants
For any simplex T ∈ O ( ) σ , the function u ∈ G δ β (S ; ) the polynomial approximation of u| T , u ∈ G δ β (S ; ) is obtained by applying Proposition 3 toû T := u| K T • F T : With u p T as in (34) The bound (20) with c T > 0 sufficiently large, independent of ensures that there exists c(κ, σ ) > 0 such that the associated K T satisfies

Exponential convergence in broken Sobolev norms
Proposition 4 For u ∈ G δ β (S ; ) with (5), there are b, C > 0 (depending on u) such that for every p ≥ 1 and forĨ p from (35) holds with ≥ 1 Here C > 0 depends on u and σ , but is independent of p, and H 1 (O ( ) σ ) denotes the broken H 1 space over O ( ) , with corresponding norm.
Proof Since S consists of finitely many singular points c, by localization and superposition, we may assume wlog. S = {c} and denote by β = β c > −2.
By assumption, K T ⊂ and, by (20), dist(K T , S ) ≥ c T σ k . Then, for u ∈ G δ β (S ; ) and for this T ∈ L k ,û T := u| K T • F T is smooth in K and satisfies, by (33) with G = K T and G = K , We obtain for |u| m,K T using (18) and (22) We define u p T ∈ Q p (T ) ⊂ P pd (T ) as in (34). From (30), for every integer 3 ≤ s ≤ p and with q,r as in (31) and for j = 0, 1, 2, Using the κ-shape regularity of T ∈ L k ⊂ M (k) σ ∈ M κ,σ , we find h T B T F ≤ κh T (eg. [5, (Chap. II, (6.9)]) and, by (22) and (33), that h T κσ k so that for every We obtain for j = 0, 1, 2 the bound Transporting to T = F T ( T ) ∈ L k , we find for β c = −1 − b c and j = 0, 1, 2.
For T ∈ O ( ) , we define the piecewise polynomial interpolantĨ p u| T by (34). Theñ I p u coincides with u in the vertices of all T ∈ O ( ) and is in particular continuous in these vertices; it is, however, in general discontinuous across edges and faces. Using the finite cardinality (21), and summing the bound (38) with j = 0, 1 over layers L 1 , ..., L −1 , we obtain withC := C u κd and We have for s < p with the recursion formula (z + 1) = z (z) that In the case δ ≥ 1, for all p ∈ N choosing s = s( p) > 0 such that there holds with a constant c > 1 independent of p (to be selected below) p = cs δ (this ensures s < p in the upper bounds (40)), we obtain where the constant hidden in is independent of the polynomial degree p.
In the case 0 < δ < 1, we choose s = p/2 in (40) and obtain for some 0 < b < 1. Inserting this bound into (39) completes the proof. ]] F ∈ P pd (F). We build a continuous, piecewise polynomial interpolant by successively lifting these polynomial trace jumps ofĨ p while retaining its consistency, in particular the analytic estimates (30).

Polynomial trace lifting in O
First, we lift jumps on interelement edges E ∈ E T and, second, in dimension d = 3 also for all interelement faces F ∈ F T , for every We recapitulate from [28, Lemma 15, Thm. 1] the required lifting and the stability estimates. Consider the reference simplex T ⊂ R d , d = 2, 3. Given a piecewise polynomial functionĝ p of degree p on each F ∈ F T that is continuous on ∂ T , in [28, Lemma 15, Thm. 1], a polynomial trace liftingv As . With the polynomial inverse inequality (see, e.g., [36]) on each face F ⊂ ∂ T we get (with a possibly different constantĈ > 0 which is independent of p) Squaring this and scaling T to Iterating (42) twice, fromÊ ⊂ ∂F toF ⊂ ∂T toT , we obtain for g p ∈ P p 0 ( E) a polynomial edge lifting L T , E ( g p ) ∈ P p ( T ) on the reference simplex T ⊂ R 3 with Squaring (44) and scaling to Let now d = 3 and let F, F ∈ F T be two distinct faces which share edge E = F ∩ F . Using (42) in dimension d = 2 and scaled to T , we lift g p = [[Ĩ p u]] E ∈ P pd 0 (E) twice, once into F and once into F , resulting in a v p ∈ C 0 (F ∪ F ), v p ∈ P pd (F)∪P pd (F ), and v p | ∂ F∪F = 0 which satisfies (43) with F in place of T . We may therefore extend this continuous, piecewise polynomial function v p from F ∪ F by zero to a functioñ v p ∈ C 0 (∂ T ) which is, on each F ∈ F T , a polynomial of total degree at most pd. There exists a lifting L T ,F (ṽ p ) ∈ P pd (T ) such that for each F ∈ F T we have for which E ∈ E T by the edge-lifting operator By κ shape regularity, #{T ∈ O ( p) σ : E ∈ E T } is bounded independently of p and of the particular edge E by an absolute constant depending only on κ. WithĨ p in (35), we defineȊ Then,Ȋ p u is continuous across edges E ∈ E T for every T ∈ O . (49) The first term was bound in Propostion 4. We bound the second term.
The multiplicative trace inequality implies for a κ-shape regular simplex T ⊂ R d with diameter h T that for every F ∈ F T and for every ϕ ∈ H 1 (T ) reads Iterating this for where the implied constant depends only on κ.
Using (38) and that h T ∼ σ k for T ∈ L k we obtain We estimate further, using the stability of the lifting L T ,F and (51), Recalling (47), we bound for j = 0, 1 We use (38) for the first term, and (53) for the second term to conclude for j = 0, 1 Using again that T ∈ L k satisfies h T ∼ σ k , we insert into (54) and arrive at Inserting this and the bound (53) into (49), we obtain for u − I p u H 1 (O ( ) σ ) exactly once more the bound (39) (with a slightly higher power of p). Absorbing the polynomial factor into the exponential, we conclude the exponential error bounds from Proposition 4, i.e., also for the resulting continuous hp-interpolant I p u defined in (48) in O ( p) σ using again (41).

Enforcement of homogeneous Dirichlet boundary conditions
The preceding polynomial trace liftings allow to obtain interpolation operators {I hp N } N which preserve homogeneous Dirichlet boundary conditions on ∂ . For simplicity, we discuss this only for the case of global homogeneous Dirichlet boundary conditions, i.e., for u| ∂ = 0 (the argument being local, i.e., element-by-element, allows to treat homogeneous Dirichlet boundary conditions also on a proper subset D ⊂ ∂ , as long as D coincides with the closure of a set of boundary faces). In space dimension d = 3, for T ∈ O ( ) σ with F ∈ F T satisfying F ⊂ ∂ , it holds u| F = 0. Hence, on T we may adjust the (nodally exact) hp (quasi-)interpolant (Ĩ p u)| T by lifting its trace (Ĩ p u)| F = −(u −Ȋ p u)| F on the boundary face F ∈ F T ∩ ∂ exactly as in (48), in particular preserving the exponential convergence bound (55). In space dimension d = 2, a corresponding polynomial edge-lifting can like wise be applied. In space dimension d = 1, is a bounded interval on R. The nodal exactness of the hp (quasi-)interpolant (Ĩ p u) implies that is satisfies the zero Dirichlet boundary conditions, so that trace lifting is not necessary in space dimension d = 1.

Approximation in T ( )
Exponential consistency errors for error contributions of the hp-interpolant from the terminal layer will be obtained essentially by a Bramble-Hilbert style scaling ("h-version FEM") argument combined with the exponentially small meshwidth of elements T ∈ T with β m := −1 − β c ∈ (0, 1/2), and 0 < ε < 1/2 − β m = 3/2 + β c ). For ⊂ R 2 , N 2 β ( ) ⊂ H 1+θ ( ) for some θ > 0 (cp. [19]), which implies θ = 2 + β c − ε > 1/2). From Proposition 2 items (i)-(iii), the collections T c := {T ∈ T ( ) σ : T ∈ ω c }, c ∈ S have uniformly bounded (w.r. to ) cardinality and shape regularity. We construct the approximation by use of a Scott-Zhang quasi-interpolating projection operator . This operator is constructed by choosing faces (edges) F z for each vertex of T c via where φ z ∈ S 1 (T c ) denotes the hat function associated with the vertex z and φ z ∈ S 1 (F z ) denotes the unique dual basis function associated with φ z (in the sense that F z φ z φ z dx = δ zz for all nodes z , see e.g. [33]). We define c := ∂( T c ) ∩ (the interface of T c and O ( ) Lemma 3] shows that J c (·)| c is a H 1/2 -stable projection (with constants depending only on shape regularity of T c ) and hence is quasi-optimal in the sense We construct the approximation u c ∈ S 1 (T ( ) c ) by setting u c = J c u on vertices in T c \ c and u c = I p u on the remaining vertices in c . The estimate (56) allows us to bound the difference This leads to By Proposition 2, items (ii) and (iii) it holds that diam( T c ) σ , we obtain with (55) and u S := c∈S u c u − u S H 1 ( T ( ) σ ) ≤ c(κ, σ )σ θ + C exp(−bp 1/δ ). completes the construction of I hp in (1). Choosing p δ concludes the proof for δ ≥ 1.

Concluding remarks
We have proved the exponential convergence rate (23) for continuous hp-FE approximations of κ shape-regular, simplicial meshes with geometric refinement to analytic functions with isolated point singularities at a finite set S in a bounded domain D ⊂ R d , d = 1, 2, 3. Apart from κ-shape regularity and σ -geometric mesh refinement the proof did not assume further structural assumptions on the triangulations. In particular, simplicial partitions which are obtained by successive bisection tree refinement in the course of adaptive subdivisions are admissible. The approximation results imply the exponential convergence rate exp(−b 3 √ N ) for second order, elliptic PDEs in polygons D ⊂ R 2 (where S denotes the set of corners of D) where solutions belong to the analytic class (i.e., where δ = 1) which are considered, for example, in [2,14,21,31,32]. Theorem 1 also implies the exponential convergence rate exp(−b 4 √ N ) for hp-approximations of electron densities in DFT, due to their analyticity [17] and due to quasioptimality of Galerkin approximations shown, for example, in [3,6] and the references there. In this application, S denotes the set of nuclei, whose centers c ∈ S are assumed known. The extension of [31,32] to Gevrey-regular solutions is essential in this case, as analyticity of electron densities can not be expected, generally, in the presence of empirical (pseudo-)potential functions constructed, for example, from smooth partitions of unity.
Also, unlike other approaches such as plane waves, hp-approximations do not, a priori, impose any specific functional form of the electron densities. Due to the locality of approximation and the separation (2) of the points c ∈ S , we may apply Theorem 1 in each neighborhood ω c , c ∈ S , implying that the total number of degrees of freedom to achieve accuracy ε > 0 in the norm H 1 (D) scales as O(#(S )| log ε| 4 ), i.e. linear scaling in the number #(S ) of nuclei and polylogarithmic scaling in the target accuracy ε. This is analogous to what is reported recently for discontinuous Galerkin discretizations in [24], where Proposition 4 can be used a starting point of proof of an exponential convergence result on tetrahedral meshes; for geometric meshes of hexahedra, analogous results can be found in [30,Sec. 5.2.2]. Exponentially convergent quadrature algorithms for the (singular) electron-pair integrals are available in [11]. The results in the present note are confined to space dimension d ≤ 3. The approach generalizes, however, directly to hp-approximations of point singularities in any dimension d with exponential rate; we remark that I hp N in the terminal layers T ( ) σ of the geometric meshes M σ introduced in Sect. 4.2.8 were built from loworder quasi-interpolants of Scott-Zhang type, which do not require continuity of u near S . Likewise, the exponential convergence rate bound (1) will remain true for linear polynomial degree vectors and, more generally, for degree vectors of bounded variation as introduced in [30]. Also, our construction of I hp N was based on a-priori knowledge of the singular support S . In case S is not known a-priori, adaptive hp-approximations as those in [8] are a method of choice. For these algorithms, the present results establish convergence rate benchmarks.