Bulk-surface virtual element method for systems of PDEs in two-space dimension

In this paper we consider a coupled bulk-surface PDE in two space dimensions. The model consists of a PDE in the bulk that is coupled to another PDE on the surface through general nonlinear boundary conditions. For such a system we propose a novel method, based on coupling a virtual element method [Beir\~ao da Veiga et al., 2013] in the bulk domain to a surface finite element method [Dziuk&Elliott, 2013] on the surface. The proposed method, which we coin the Bulk-Surface Virtual Element Method (BSVEM) includes, as a special case, the bulk-surface finite element method (BSFEM) on triangular meshes [Madzvamuse&Chung, 2016]. The method exhibits second-order convergence in space, provided the exact solution is $H^{2+1/4}$ in the bulk and $H^2$ on the surface, where the additional $\frac{1}{4}$ is required only in the simultaneous presence of surface curvature and non-triangular elements. Two novel techniques introduced in our analysis are (i) an $L^2$-preserving inverse trace operator for the analysis of boundary conditions and (ii) the Sobolev extension as a replacement of the lifting operator [Elliott&Ranner, 2013] for sufficiently smooth exact solutions. The generality of the polygonal mesh can be exploited to optimize the computational time of matrix assembly. The method takes an optimised matrix-vector form that also simplifies the known special case of BSFEM on triangular meshes [Madzvamuse&Chung, 2016]. Three numerical examples illustrate our findings.

PDEs posed on the surface Γ := ∂Ω through either linear or non-linear coupling, see for instance [37]. For ease of presentation, we confine the exposition to the case m = n = 1, even if the whole study applies to the case of arbitrary m and n. In the time-independent case, let u(x) and v(x) be the bulk and surface variables obeying the following elliptic bulk-surface (BS) linear problem: x ∈ Ω; where α, β > 0, while ∆ and ∆ Γ denote the Laplace and Laplace-Beltrami operators respectively, and ν denotes the outward unit normal vector field on Γ (see Appendix A for notations and definitions). The above model (1) was considered in [27]. For a nonlinear time-dependent generalisation, let u(x, t) and v(x, t) be the bulk and surface variables obeying the following bulk-surface reaction-diffusion system (BSRDS): where T is the final time and q(u), r(u, v), s(u, v) are possibly nonlinear functions. The model comprises several time-dependent BSPDE models existing in the literature, see for example [22,29,37,41]. The quickly growing interest toward BSPDEs arises from the numerous applications of such PDE problems in different areas, such as cellular biological systems [22,29,36,40] or fluid dynamics [15,17,34], among many other applications. Among the various state-of-the art numerical methods for the spatial discretisation of BSPDEs existing in the literature we mention finite elements [27,33,37,38], trace finite elements [32], cut finite elements [17] and discontinuous Galerkin methods [20].
On the other hand, the Virtual Element Method (VEM) for bulk-only PDEs is a recent extension of the well-known finite element method (FEM) for the numerical approximation of several classes of PDEs on flat domains [5] or surfaces [31]. The key feature of VEM is that of being a polygonal method, i.e. it handles elements of a quite general polygonal shape, rather than just of triangular shape [5]. The success of virtual elements is due to several advantages arising from non-polygonal mesh generality. A non-exhaustive list of such advantages includes: (i) computationally cheap mesh pasting [12,19,31], (ii) efficient adaptive algorithms [18], flexible approximation of the domain and in particular of its boundary [23], and the possibility of enforcing higher regularity to the numerical solution [4,9,16], just to mention a few. Thanks to these advantages, several extensions of the original VEM for the Poisson equation [5] were developed for numerous PDE problems, such as heat [46] and wave equations [45], reaction-diffusion systems [1], Cahn-Hilliard equation [4], Stokes equation [8], linear elasticity [6], plate bending [16], fracture problems [11], eigenvalue problems [39] and many more.
The purpose of the present paper is to introduce a bulk-surface virtual element method (BSVEM) for the spatial discretisation of a coupled system of BSPDEs in two space dimensions. The proposed method combines the VEM for the bulk equation(s) with the surface finite element method (SFEM) for the surface equation(s). We apply the proposed method to (i) the linear elliptic BS Poisson problem (1) and (ii) the BSRDS (2).
The main novelty in the present study is devoted to error analysis. In fact, the simultaneous presence of non-triangular elements and boundary approximation error (which cannot be neglected in the context of BSPDEs, because the boundary is itself the domain of a surface PDE) provide new numerical analysis challenges. We prove that the proposed method possesses optimal second-order convergence in space, provided the exact solution is H 2+1/4 in the bulk and H 2 on the surface. H denotes the Hilbert space. This is slightly more than the usual requirement of H 2 both in the bulk and on the surface [33]. However, our analysis requires this slightly higher regularity assumption only in the simultaneous presence of a curved boundary Γ and non-triangular elements close to the boundary, which is a novel case. Otherwise, our results fall back to the known cases in the literature, see for instance [33] for the case of triangular BSFEM and [1] for the case of polygonal VEM in the absence of curvature in the boundary Γ. The extension to higher order cases would require the usage of curved elements, because otherwise the geometric error arising from the boundary approximation would dominate over the numerical error, see [24]. This challenge will be addressed in future studies. The proposed analysis has three by-products: 1. The bulk-VEM of lowest polynomial order k = 1 possesses optimal convergence also in the presence of curved boundaries. A first work in this direction is found in [10], in which the authors consider polygonal elements with a curved boundary that match the exact domain in order to avoid errors arising from boundary approximation. Subsequently, in [14], the need of matching the exact domain was removed by introducing suitable corrections in the method.
In the present work, instead, we obtain similar results in the low polynomial case k = 1 by harnessing in a novel way the geometric error estimates of BS polyhedral domains in [27].

2.
A potential alternative approach to the lifting operator used in the analysis of SPDEs [26] and BSPDEs [27] is the Sobolev extension operator. In fact, we prove that, if a function is H 2+1/4 instead of H 2 in the two-dimensional bulk, the Sobolev extension retains optimal approximation properties of lifting. Moreover, the Sobolev extension of a function has the property of preserving its W m,p class, while its lift does not because it is not C k for any positive integer k. This property, which is crucial in our analysis, is potentially beneficial to the error analysis of bulk-only or BSVEM approximation of more general PDEs, where the boundary curvature was not accounted for, see for instance [5,8,31,39,45,46]. 3. We construct a special inverse trace operator that we use for accounting for general boundary conditions. Like the standard inverse trace operator, our inverse trace maps a function v ∈ H 1 (Γ) to a function v B ∈ H 1 (Ω) such that v B H 1 (Ω) ≤ C v H 1 (Γ) , with C depending on the domain Ω. Our inverse trace has the stronger property that v B L 2 (Ω) ≤ C v L 2 (Γ) and |v B | H 1 (Ω) ≤ C v H 1 (Γ) . This means that the proposed operator preserves both L 2 and H 1 norms up to the same multiplicative constants, i.e. it is an L 2 -preserving inverse trace.
The proposed method has all the benefits of polygonal meshes, two of which will be illustrated in the present work. First, the usage of suitable polygons drastically reduces the computational complexity of matrix assembly on equal meshsize in comparison to the triangular BSFEM. Similar results are obtained in the literature through other methods, such as trace FEMs [32] or cut FEMs [17]. Second, a curved portion of the boundary can be approximated with a single element with many edges. The BSVEM lends itself to other advantages due to its polygonal nature, such as efficient algorithms for adaptivity or mesh pasting, see for instance [18]. These aspects will be addressed in future studies. Hence, the structure of our paper is as follows. In Section 2, we elaborate on the weak formulations, existence and regularity for problems (1) and (2). In Section 3, we introduce polygonal BS meshes, analyse geometric error, define suitable function spaces, analyse their approximation properties and present the spatial discretisation of the considered BSPDE problems. In Section 4, we carry out the convergence error analysis for the parabolic case, the main result being optimal second-order spatial convergence in the L 2 norm, both in the bulk and on the surface. In Section 5 we present an IMEX-Euler time discretisation of the parabolic problem. In Section 6 we show that polygonal meshes can be exploited to significantly reduce the computational time of the matrix assembly. In Section 7 we illustrate our findings through three numerical examples. The first example shows (i) optimised matrix assembly and (ii) optimal convergence for the elliptic problem (1). The second example shows optimal convergence in space and time for the parabolic problem (2) for the linear case. The third example compares the BSFEM-and BSVEM-IMEX Euler approximations of the wave pinning RDS considered in [22]. In Appendix A we provide some basic definitions and results required in our analysis. In Appendix B we report some lengthy proofs involved in Section 4.
The problem of existence, uniqueness and regularity for the parabolic problem (5) is much more complicated and strictly depends on the nature of the kinetics q(·), r(·) and of the coupling kinetics s(·). For the remainder of this work we will adopt the following set of assumptions.
Assumption 1 (Existence, uniqueness and regularity for problem (5)). We assume that: • Γ is a C 3 surface, q, r, s are C 2 and globally Lipschitz functions.
• There exists a unique solution (u, v) that fulfils where T > 0 is the final time and C > 0 depends on Ω, q C 2 (R) , r C 2 (R 2 ) and s C 2 (R 2 ) .
In many applications, assuming globally Lipschitz kinetics is too restrictive and an ad-hoc analysis is required, see for instance [29]. However, there are notable examples of BSRDS with globally Lipschitz kinetics, such as the wave pinning model studied in [22] and considered in the numerical example in Section 7.3. From here onwards, we shall assume that the weak parabolic problem (5) has a unique and sufficiently regular solution.

The Bulk-Surface Virtual Element Method
In this section we will introduce the Bulk-Surface Virtual Element Method (BSVEM) through the following steps: • describe the polygonal BS meshes that will be used in the method (Subsection 3.1); • quantify the geometric error arising from polygonal approximation of BS domains (Subsection 3.2); • introduce the discrete function spaces and bilinear forms used in the method and their approximation properties (Subsections 3.3-3.4); • present the spatially discrete formulations of the elliptic-and parabolic problems (1) and (2), respectively (Subsection 3.5).

Polygonal bulk-surface meshes
Let h > 0 be a positive number called meshsize and let Ω h = ∪ E∈E h E be a polygonal approximation of the bulk Ω, where E h is a set of non-degenerate polygons. The polygonal bulk Ω h automatically induces a piecewise linear approximation Γ h of Γ, defined by Γ h = ∂Ω h , exactly as in the case of triangular meshes, see for example [27]. Notice that we can write Γ h = ∪ F ∈F h F , where F h is the set of the edges of Ω h that lie on Γ h . We assume that: (F1) the diameter of each element E ∈ E h does not exceed h; (F2) for any two distinct elements E 1 , E 2 ∈ E h , their intersection E 1 ∩ E 2 is either empty, or a common vertex, or a common edge.
(F3) all nodes of Γ h lie on Γ; (F4) every edge F ∈ F h is contained in the Fermi stripe U of Γ (see Fig. 1 for an illustration).
(V1) there exists γ 1 > 0 such that every E ∈ E h is star-shaped with respect to a ball of radius where h E is the diameter of E; (V2) there exists γ 2 > 0 such that for all E ∈ E h , the distance between any two nodes of E is at least γ 2 h E .
Assumptions (F1)-(F4) are standard in the SFEM literature, see for instance [26], while assumptions (V1)-(V2) are standard in the VEM literature, see for instance [5]. The combined assumptions (F1)-(V2) will prove sufficient in our BS setting. In the following definitions and results we provide the necessary theory for estimating the geometric error arising from the boundary approximation.
Definition 1 (Essentials of polygonal BS meshes). An edgeē of any element E ∈ E h is called a boundary edge ifē ⊂ Γ h , otherwiseē is called an inner edge. Let BE(E) and IE(E) be the sets of boundary and inner edges of E, respectively. An element E ∈ E h is called an external element if it has at least one boundary edge, otherwise E is called an internal element. Let Ω B be the discrete narrow band defined as the union of the external elements of Ω h as illustrated in Fig. 1(b). From Assumption (F4), for any boundary edgeē, we have that a(ē) ⊂ Γ, where a is the normal projection defined in Lemma A1. Hence, for sufficiently small h and for all E ∈ E h , it is possible to define the exact elementȆ as the compact set enclosed by the edges IE(E) ∪ {a(ē)|ē ∈ BE(E)}, see Fig. 2 for an illustration.
Remark 1 (Properties of polygonal BS meshes). For any BS mesh (Ω h , Γ h ) of meshsize h > 0 it holds that: • for sufficiently small h > 0, the discrete narrow band Ω B is contained in the Fermi stripe U as shown in Fig. 1(b) (blue colour); • the collection of all exact elements is a coverage of the exact bulk Ω, that is ∪ E∈E hȆ = Ω.
Let N ∈ N and let x i , i = 1, . . . , N , be the nodes of Ω h , which can be ordered in an arbitrary way. However, if Ω has a rectangular shape and the nodes are ordered along a Cartesian grid, the matrices associated with the method will have a block-tridiagonal structure. Let M ∈ N, M < N and assume that the nodes of Γ h are x k , k = 1, . . . , M , i.e. the first M nodes of Ω h . Throughout the paper we need the following reduction matrix R ∈ R N ×M defined as R : for all k = 1, . . . , M , where δ ik is the Kronecker symbol. The reduction matrix R fulfils the following two properties: In what follows, we will use the matrix R for an optimised implementation of the BSVEM.         Fig. 3(c) according to Def. 1. Figure 3: On the mapping G E between a triangular element E and its curved counterpartȆ. In the case of one boundary edge (subplots (a)-(b)), a C 2 homeomorphism G E : E →Ȇ is known to exist, see [27]. With two adjacent boundary edges (subplots (c)-(d)) the mapping G E cannot be smooth, because G E maps the non-smooth curveē 1 ∪ē 2 onto the smooth curve a(ē 1 ) ∪ a(ē 2 ). The approach proposed in Lemma 1 solves this issue.

Variational crime
We now consider the geometric error due to the boundary approximation. Since the surface variational crime in surface finite elements is well-understood [26], we will mainly focus on the variational crime in the bulk. To this end, it is useful to analyse the relation between any element E ∈ E h and its exact counterpartȆ, see Def. 1. For the special case of triangular meshes with at most one boundary edge per element ( Fig. 3(a)-(b)), there exists a C 2 homeomorphism G E : E →Ȇ that is quadratically close to the identity with respect to the meshsize, see [27]. Instead, when an element E has adjacent, non-collinear boundary edges (which can occur even in triangular meshes, see Fig.  3(c)-(d)), such a smooth homeomorphism does not exist. In fact, a smooth mapping cannot map the adjacent boundary edges of E -a curve with corners -onto a portion of the smooth curve Γ. However, in the following result we show the existence of a homeomorphism between E andȆ with slightly weaker regularity, which is sufficient for our purposes.
Lemma 1 (Parameterisation of the exact geometry). Let h be sufficiently small (depending on the curvature of Γ and the mesh regularity parameter γ 1 ) and let E h fulfil assumptions (F1)-(V2). There exists a homeomorphism G : Ω h → Ω such that G ∈ W 1,∞ (Ω h ) and G| Ω h \Ω B = Id; where a is the normal projection defined in Lemma A1, JG is the Jacobian of G and C is a constant that depends on Γ and the constants γ 1 , γ 2 are those considered in Assumptions (V1)-(V2).
Proof. Let E ∈ E h . By Assumption (V2), E is star-shaped with respect to a ball B E of diameter R E ≥ γ 2 h E , let x E be the center of such ball ( Fig. 4(a)). Forē ∈ BE(E) let Tē be the triangle spanned byē and x E ( Fig. 4(b)). Let us now consider the collection of all the Tē's defined as follows: We need to prove that BT h is quasi-uniform, i.e we need to prove that for all E ∈ E h andē ∈ BE(E), the triangle Tē has an inscribed ball of diameter greater or equal tō γh E , where the constantγ > 0 depends on γ 1 and γ 2 , only. To this end, if hē is the height of Tē relative to the basisē ( Fig. 4(c)), then we have In addition, since no edge of Tē is longer than h E , our claim follows. Let nowTē be the curved triangle corresponding to Tē ( Fig. 4(d)), which again is well-defined by assuming h sufficiently small depending on Γ andγ, which in turn depends on γ 1 and γ 2 . Since the triangulation BT h is quasi-uniform, then, from [27, Proposition 4.7 and its proof] there exists a diffeomorphism Gē : Tē →Tē such that where C depends only onγ, which in turn depends only on γ 1 and γ 2 . We are ready to construct the mapping G : Ω h → Ω as follows: Property (18) can now be rephrased by saying that G, restricted to any inner edge of Ω h , is the identity. In addition, since all the Gē's are homeomorphisms, we obtain that G is a homeomorphism between Ω h and Ω. Finally, from (17), (19), (20) and (22) we obtain the desired estimates (10)- (13).
Remark 2 (Virtual elements for bulk-only PDEs). Lemma 1 has an important consequence in the analysis of boundary approximation for VEMs for the case of bulk PDEs. For triangular finite elements, boundary approximation is a well-understood topic, see for instance [21]. For VEMs, the first work in this direction is [10], in which a VEM on exact curved polygons is considered, in order to take out the geometric error. However, in the lowest-order VEM on polygonal domains, it is empirically known that the geometric error does not prevent optimality, as discussed in [14]. To the best of the authors' knowledge, the present study provides, as a by-product, the first rigorous proof of this fact.
Thanks to Lemma 1 it is possible to define bulk-and surface-lifting operators. Lemma 1 also enables us to show the equivalence of Sobolev norms under lifting as illustrated next.  Lemma 2 (Equivalence of norms under lifting). There exists two constants c 2 > c 1 > 0 depending on Γ and γ 2 such that, for all V : Ω h → R and for all W : Γ h → R, Proof. Estimates (23) We are ready to estimate the effect of lifting on bulk and surface integrals.
Lemma 3 allows to rephrase integrals on the exact domain as integrals on the discrete domain, and vice-versa, up to a small error that is O(h) in the bulk and O(h 2 ) on the surface.

Remark 3 (Preservation of regularity under lifting).
For E ∈ E h the inverse lift of an H 2 (Ȇ) function is not, in general, H 2 (E), cf. Remark A1 and Lemma 1. This problem does not arise in the context of triangular BSFEM in the presence of curved boundaries, see for instance [33], or bulk virtual elements in the absence of curved boundaries, see for instance [1]. In the context of bulkor bulk-surface VEM in the presence of curved boundaries, our analysis requires full H 2 -regularity of the exact solution mapped on the polygonal domain. Because u −ℓ does not retain the regularity class of u, we need an alternative mapping of u on Γ h . To this end, we recall the following Theorem.
We are now able to approximate the exact solution u throughũ |Ω h , that is the restriction to Ω h of the Sobolev extensionũ of the exact solution u. We present the following variant of Lemma 3 in order to quantify the geometric error of the Sobolev extension.

Virtual element space and operators in the bulk
In this section, we define virtual element spaces on polygons and polygonal domains by following [7]. Let E be a polygon in R 2 . A preliminary virtual element space on E is given bỹ where P 1 (E) is the space of linear polynomials on the polygon E. The functions in V(E) are not known in closed form, but we are able to use them in a spatially discrete method, hence the name virtual. Let us consider the elliptic projection Using Green's formula, it is easy to see that the operator Π ∇ E is computable, see [3] for the details. The so-called enhanced virtual element space in two dimensions is now defined as follows: The practical usability of the space V(E) stems from the following result.
Hence, the nodal values constitute a set of degrees of freedom.
For s = 1, 2 we define the broken bulk Sobolev seminorms as |u| s,Ω,h := E∈E h |u |E | H s (E) . The approximation properties of the space V(E) are given by the following result.
where C is a constant that depends only on γ 1 .
. If E is a convex polygon, from the properties of the Poisson problem on convex Lipschitz domains, it holds that V(E) ⊂ H 2 (E), see for instance [42]. [42]. In either case, V(E) ⊂ C 0 (E), see Theorem A4. We will account for this regularity issue in devising a numerical method with optimal convergence.
The discontinuous and continuous bulk virtual element spaces are defined by pasting local spaces: Thanks to Remark 4, the only source of discontinuity in V Ω,h are jumps across edges. In V Ω we consider the Lagrange basis {ϕ i , i = 1, . . . , N } where, for each i = 1, . . . , N , ϕ i is the unique V Ω function such that ϕ i (x j ) = δ ij , for all j = 1, . . . , N , with δ ij being the Kronecker symbol. The set {ϕ i , i = 1, . . . , N } is a basis of V Ω thanks to Proposition 1.
In the remainder of this section, let E be an element of Ω h . The stabilizing form The As shown in [3], is not a new projector, the presentation and the analysis of the method benefit from the usage of the equivalent definition (43). Moreover, since Π 0 E = Π ∇ E , the boundedness property of projection operators in Hilbert spaces translates to We are now ready to introduce the approximate The following result easily follows from the construction of the bilinear forms of a E and m E .
Proposition 3 (Stability and consistency). The bilinear forms a h and m h are consistent, i.e. for all v ∈ V(E) and p ∈ P 1 (E) The bilinear forms a h and m h are stable, meaning that there exists two constants 0 < α * < α * depending on γ 2 such that, for all v ∈ V(E) Proof. See [5].
We observe from (48) and (49) that the error in the approximate bilinear forms a E and m E is not a function of the meshsize h, see also [5]. Nevertheless, we will show that the method retains optimal convergence thanks to the consistency properties (47). The global bilinear forms are defined by pasting the corresponding local bilinear forms as follows: The bilinear form m h is not sufficient to discretise load terms like Ω f ϕ, because f is not in the space V Ω . We resolve this issue by combining the approaches in [1] and [30,31]. From [1] we take the usage of the projection operator Π ∇ , while from [30,31] we take the usage of the Lagrange interpolant. In the context of virtual elements, the Lagrange interpolant is defined as follows.
Definition 3 (Virtual Lagrange Interpolant). Given an element-wise continuous function f : The following result provides an estimate of the interpolation error in the bulk.
Proposition 4 (Interpolation error in the bulk).
Unlike projection operators, we may not assert that the interpolant I Ω is bounded in L 2 (Ω h ). However from (51) we have the quasi-boundedness of I Ω in L 2 (Ω h ): where C > 0 depends only on γ 1 .

Finite element space and operators on the surface
Let F ∈ F h be an edge of the approximated curve Γ h . The local finite element space on F is the space P 1 (F ) of linear polynomials on F . The global finite element space on Γ h is defined by Before introducing the spatially discrete formulations, we are left to treat terms like Γ gϕ, since g is not in the boundary finite element space V Γ .
Definition 4 (Surface Lagrange interpolant). If g : Γ h → R is a continuous function, the Lagrange interpolant I Γ (g) of g is the unique V Γ function such that I Γ (g)(x i ) = g(x i ) for all i = 1, . . . , M .
We consider the broken surface Sobolev norm |v| 2,Γ,h := F ∈F h |V |F | H 2 (F ) . The following are basic properties of Lagrange interpolation.

The spatially discrete formulations
We are now ready to introduce the BSVEM discretisation of the weak problems (4) and (5). The discrete counterpart of the elliptic problem (4) is: find U ∈ V Ω and V ∈ V Γ such that for all ϕ ∈ V Ω and ψ ∈ V Γ . We express the spatially discrete solution (U, V ) in the Lagrange bases as follows: Hence, problem (56) is equivalent to: find ξ := (ξ i , . . . , ξ N ) T ∈ R N and η := (η 1 , . . . , for all j = 1, . . . , N and l = 1, . . . , M . We consider the matrices a Ω i,j := a h (ϕ i , ϕ j ); and m Ω i,j := m h (ϕ i , ϕ j ), i, j = 1, . . . , N ; Moreover, we consider the column vectors f ∈ R N and g ∈ R M defined by By using (53) we can now rewrite the discrete formulation (58) in matrix-vector form as a block (N + M ) × (N + M ) linear algebraic system: In compact form, (62) reads which is uniquely solvable, as a consequence of the positive semi-definiteness of a h and the positive definiteness of m h . The spatial discretisation of the parabolic problem (5) Remark 5 (Special cases). If every element E ∈ E h is convex or f is linear, optimal convergence is retained by replacing the term m h (Π 0 I Ω (q(Π 0 U )), ϕ) with m h (I Ω (q(U )), ϕ), i.e by removing the projection operator Π 0 . By expressing the time-dependent semi-discrete solution (U, V ) in the Lagrange bases as follows the fully discrete problem can be written in matrix-vector form as an (M + N ) × (M + N ) nonlinear ODE system of the form: Remark 6 (Implementation). Thanks to the reduction matrix R, we are able to implement the spatially discrete problems (62) and (67) by using only two kinds of mass matrix (M Ω and M Γ ), two kinds of stiffness matrix (A Ω and A Γ ) and R itself. In previous works on BSFEM as illustrated in [37], five kinds of mass matrix were used to evaluate the integrals appearing in the spatially discrete formulation. We stress once again that, since the pre-existing BSFEM is a special case of the proposed BSVEM, this work provides, as a by-product, an optimised matrix implementation of the BSFEM.

Stability and convergence analysis
The spatially discrete parabolic problem (64) fulfils the following stability estimates.
Lemma 6 (Stability estimates for the spatially discrete parabolic problem (64)). There exists C > 0 depending on q, r, s and Ω such that Proof. The proof relies on standard energy techniques. By choosing ϕ = U and ψ = V in (64), using the Lipschitz continuity of q, r, s, the Young's inequality and summing over the equations we have where c > 0 is arbitrarily small, thanks to the Young's inequality. By applying (24), (A.11) and (48) to the last term in (70) we can choose c such that we have By using (49) into (71) we have An application of Grönwall's lemma to (72) yields which yields (68) after an application of (24). Similarly, by choosing ϕ =U and ψ =V in (64) and summing over the equations we have where we have exploited the Lipschitz continuity of q, r, s and the Young's inequality and (49). From (74) we immediately get By applying Grönwall's lemma and then using (48)-(49) we obtain (69).
To derive error estimates for the spatially discrete solution we need suitable bulk and surface Ritz projections. The surface Ritz projection is taken from [28], while the bulk Ritz projection is tailor-made.

Definition 5 (Surface Ritz projection). The surface Ritz projection of a function
Theorem 3 (Error bounds for the surface Ritz projection). The surface Ritz projection fulfils the optimal a priori error bound where C > 0 depends only on Γ and ℓ is the lifting.
We now define a suitable Ritz projection in the bulk.
Definition 6 (Bulk Ritz projection). The bulk Ritz projection of a function u ∈ H 1 (Ω) is the unique function Ru ∈ V Ω such that where −ℓ is the inverse lifting.
In the following theorems we show that the bulk Ritz projection of a sufficiently regular function fulfils optimal a priori error bounds in H 1 (Ω), H 1 (Γ), L 2 (Γ) and L 2 (Ω) norms.

Theorem 4 (H 1 (Ω) a priori error bound for the bulk Ritz projection). For any
where ℓ is the lifting and C depends on Ω and the constants γ 1 and γ 2 in (V1)-(V2).
Proof. We set e h := Ru −ũ. From (27), (33), (40), (47), (48) and (78) we have By using (24), (33) and (80) we get In order to prove the L 2 convergence, care must be taken about inverse trace operators. A consequence of Theorem A3 in the Appendix A is that, given v ∈ H 1 (Γ), there exists v B ∈ H 1 (Ω) such that Tr(v B ) = v and v B H 1 (Ω) ≤ C v H 1 (Γ) . However, for our purposes, we need the existence of a constant C > 0 such that the bounds v B L 2 (Ω) ≤ C v L 2 (Γ) and |v B | H 1 (Ω) ≤ C v H 1 (Γ) are simultaneously fulfilled, namely we need a L 2 -preserving inverse trace operator. In the following result, we prove the existence of such an operator under special assumptions on the regularity of Γ.
Lemma 7 (L 2 -preserving inverse trace). Assume the boundary Γ is C 3 . Then, for any Proof. With δ as defined in Theorem A1, we take 0 ≤ δ 0 ≤ δ and we define By using (A.9) we have Since the decomposition (d(x), a(x)) is unique (see Lemma A1) and all the points x ∈ Γ s share the same distance d(x), then the mapping a s := a |Γs : Γ s → Γ is invertible. Moreover, since a 0 = Id |Γ (which implies ∇ Γ a = 1) and a ∈ C 2 (U ) (see Remark A1), we can choose δ 0 small enough such that 0 < c ≤ ∇ Γs a s ≤ C. Hence, a −1 s is C 2 as well, which implies that ∇ Γ a −1 s ≤ C. Hence by setting y = a(x) we have where · is the Euclidean norm. By combining (83) and (84) we obtain the first inequality in (81), since δ 0 depends only on Γ. An application of the chain rule and Leibniz's rule yields Since d(x) and a(x) are both C 2 on the compact set U δ ⊃ U δ 0 and −δ 0 ≤ d(x) ≤ 0, we obtain where | · | is the absolute value. Thanks to the continuous pasting in (82), it is easy to show that the bound (85) for the distributional gradient ∇v B still holds on the junction Γ δ 0 between U δ 0 and Ω \ U δ 0 . From (A.9), (85) and the Young's inequality we have which proves the second inequality in (81).
Proof. See Appendix B.
Theorem 6 (Convergence of the BSVEM for the parabolic case). Assume that the kinetics q, r, s are C 2 and globally Lipschitz continuous (or at least Lipschitz on the range of the discrete solution). Assume that the exact solution (u, v) of the parabolic problem (2) is such that u, u t ∈ L ∞ ([0, T ]; H 2+1/4 (Ω)) and v, v t , Tr(u), Tr(u t ) ∈ L ∞ ([0, T ]; H 2 (Γ)). Let (U, V ) be the solution of (64). Then it holds that where the constant C depends on the diffusion coefficient d u , the final time T and on the following norms: Proof. See Appendix B.

Time discretisation of the parabolic problem
For the time discretisation of the semi-discrete formulation (67) of the parabolic problem (2) we use the IMEX (IMplicit-EXplicit) Euler method, which approximates diffusion terms implicitly and reaction-and boundary terms explicitly, see for instance [30]. This choice is meant to make the implementation as simple as possible. However, if the boundary conditions and/or the reaction terms are linear, the aforementioned terms can be approximated implicitly without involving any additional nonlinear solver. We choose a timestep τ > 0 and we consider the equally spaced discrete times t n := nτ , for n = 0, . . . , N T , with N T := T τ . For n = 0, . . . , N T we denote by ξ n and η n the numerical solution at time t n . The IMEX Euler time discretisation of (67) reads for n = 0, . . . , N T − 1, where R is the reduction matrix introduced in (9). By solving for ξ n+1 and η n+1 we obtain the following time stepping scheme: for n = 0, . . . , N T − 1. Hence, two linear systems must be solved at each time step, where the matrix coefficients are the same for all n.

Construction of meshes optimised for matrix assembly
One of the benefits of using suitable polygonal elements is the reduction in computational complexity of the matrix assembly. This is useful especially for (i) time-independent problems and (ii) timedependent problems on evolving domains, where matrix assembly might take the vast majority of the overall computational time.
As an example, given an arbitrarily shaped domain Ω with C 1 boundary Γ, we construct a polygonal mesh designed for fast matrix assembly, by proceeding as follows. Suppose that the bulk Ω is contained in a square Q. We discretise Q with a Cartesian grid made up of square mesh elements and assume that at least one of such squares is fully contained in the interior of Ω (Fig. 5(a)). Then we discard the elements that are outside Q (Fig. 5(b)). Finally, we project on Γ the nodes that are outside Ω, thereby producing a discrete narrow band of irregular quadrilaterals (highlighted in purple in Fig. 5(c)). The mesh is then post-processed in two steps as illustrated in Fig. 6 to comply with Assumptions (V1)-(V2). Let ε > 0 be a threshold. First, all hǫ-close nodes are merged, see Figs. 6(a)-(b). After that, the only chance for a polygon with at most four vertices fulfilling (V1) to not fulfil (V2) is the presence of small angles. In this case, all ǫ-small angles are reduced to the null angle, see Figs. 6(c)-(d). It is worth remarking that Virtual Elements were proven to be robust with respect to distorted elements (see [7]), so ǫ can be chosen small.
The resulting mesh Ω h is made up of equal square elements, except for the elements that are close to Γ, as we can see in It is worth remarking that the advantage described in this section becomes even more striking in higher space dimension: if the embedding space has dimension D ∈ N, D ≥ 2, then the procedure described in this section reduces the computational complexity of matrix assembly from O(1/h D ) to O(1/h D−1 ). Matrix assembly optimization can be also achieved through alternative approaches, such as cut FEM [17] or trace FEM [32]. However, in these works, the authors adopt a level set representation of the boundary Γ, which we do not need in this study, as we exploit the usage of arbitrary polygons to approximate Γ.

Numerical simulations
In this section we provide three numerical examples to compare the performances of BSFEM and BSVEM. In the first example, we show (i) optimal convergence in the case of the elliptic problem (1) and (ii) the optimised matrix assembly introduced in Section 6. In the second example, we show optimal convergence in space and time in the case of the parabolic problem (2). In the third (a) Step 1. Before: some hǫ-close nodes are present.
(d) Step 2. After: all ǫ-small angles are reduced to the null angle and the element is eliminated.

Experiment 1: The elliptic problem
We start by constructing an exact solution for the elliptic problem (1) on the unit circle Ω := {(x, y) ∈ R 2 |x 2 + y 2 = 1} by using the fact that −∆ Γ xy = 4xy, i.e. the function w(x, y) := xy is an eigenfunction of the Laplace-Beltrami operator on the unit circumference Γ = ∂Ω. Specifically, in (1) we choose the following load terms f (x, y) := xy, for (x, y) ∈ Ω and g(x, y) := 2 + 5(α+2) β xy for (x, y) ∈ Γ. Here α and β are the parameters that appear in the model, such that the exact solution is given by u(x, y) := xy for (x, y) ∈ Ω and v(x, y) := α+2 β xy for (x, y) ∈ Γ. For illustrative purposes, we choose α = 1 and β = 2. We solve the problem with BSFEM and BSVEM on two respective sequences of five meshes with decreasing meshsizes. For BSVEM, we use optimised meshes as in Section 6, while for BSFEM we use quasi-uniform Delaunay triangular meshes. See Fig. 7 for an illustration. For each mesh, we measure: • the L 2 (Ω) × L 2 (Γ) and L ∞ (Ω) × L ∞ (Γ) errors of the numerical solution and the respective experimental orders of convergence (EOCs), that are quadratic in the meshsize. In this case BSVEM is almost three times as accurate as BSFEM on similar meshsizes, probably because the mesh reflects the symmetry of the problem; • the condition number cond ell of the linear system (63). The BSVEM is approximately four times as ill conditioned as BSFEM on similar meshsizes. An approach to reduce this gap is described in [13]; • only for BSVEM, the number N Ω of all elements and the number N Γ of elements that intersect the boundary. The first is proportional to 1 h 2 , while the latter is proportional to 1 h , as predicted in Section 6. This illustrates that only O( 1 h ) local matrices must be computed during matrix assembly, even if the mesh has O( 1 h 2 ) elements. In Fig. 8 we show the BSVEM numerical solution on the finest mesh. In Tabs. 1-2 we show the errors, the respective EOCs, and the numbers N Ω , N Γ and cond ell .
(a) Quasi-uniform Delaunay triangular mesh of the kind used for BSFEM.    Table 1, with meshsize h = 4.7611e-2. Left plot shows the solution u in the bulk and the right plot shows the solution v on the surface (curve) Γ.  (1). For BSVEM, EOCs in both L 2 (Ω) × L 2 (Γ) and L ∞ (Ω)× L ∞ (Γ) norms are approximately two. Moreover, matrix assembly requires the computation of   (1), solved with BSFEM on a sequence of triangular meshes with meshsizes as similar as possible as those in Tab. 1. The EOCs in both L 2 (Ω) × L 2 (Γ) and L ∞ (Ω) × L ∞ (Γ) norms are approximately two and the errors almost three times as large w.r.t. the corresponding errors of BSVEM.

Experiment 3: Nonlinear parabolic problem
In this section we compare the BSFEM-and BSVEM-IMEX Euler numerical solutions of the bulksurface wave pinning (BSWP) model considered in [22,Fig. 4]. The non-dimensionalised BSWP model seeks to find a bulk concentration b : Ω×[0, T ] → R and a surface concentration a :  where the kinetic f is of the form f (a, b) := k 0 + γ(a 2 )/(1 + a 2 ) b − a. We solve the BSWP model (92) on the unit circle Ω := {(x, y) ∈ R 2 |x 2 + y 2 ≤ 1}. We choose the parameters of the model as ε 2 = 0.001, k 0 = 0.05 and γ = 0.79. The initial condition, plotted in Fig. 9(a), is b(x, y, 0) = 2.487, for (x, y) ∈ Ω and a(x, y, 0) = 0.309 + 0.35(1 + sign(x)) exp(−20y 2 ) for (x, y) ∈ Γ. The final time is T = 4.5. For the BSFEM we take a quasi-uniform Delaunay triangular mesh, while for the BSVEM we take a polygonal mesh designed for optimised matrix assembly, as explained in Section 6. An illustrative coarser representation of the meshes typically used is shown in Fig. 7. The details of these meshes are reported in Table 5. As we can see, on almost equal meshsizes, the BSVEM generates (i) a significantly large number of boundary nodes, which translates into better boundary approximation and (ii) less number of elements in the bulk, which implies faster matrix assembly. For both spatial discretisations, the time discretisation is computed with timestep τ = 2e-3. The condition numbers cond Ω := cond(M Ω + τ d u K Ω ) and cond Γ := cond(M Γ + τ d v K Γ ) of BSVEM on our polygonal mesh are only about 2.5 times as large as those of BSFEM on a quasi-uniform Delaunay mesh. This result, which could be further improved through an orthogonal polynomial basis approach [13] or through a different time integrator, confirms the well-known robustness of VEM with respect to general polygonal meshes [7]. The BSVEM solution at various times is shown in Fig. 9.

Conclusions
In this study, we have considered a bulk-surface virtual element method (BSVEM) for the numerical approximation of linear elliptic and semilinear parabolic coupled BSPDE problems on smooth BS domains. The proposed method simultaneously extends the BSFEM for BSRDSs [37] and the VEM for linear elliptic [5] and semilinear parabolic [1] bulk PDEs. The method has a simplified vectormatrix form that can be exploited also in the special case of the BSFEM considered in [37]. We have introduced polygonal BS meshes in two space dimensions and we have shown that the geometric error arising from domain approximation is O(h) in the bulk and O(h 2 ) on the surface, where h is the meshsize, exactly as in the special case of the triangular BS meshes used in the BSFEM [27]. Suitable polygonal meshes reduce the asymptotic computational complexity of matrix assembly from O(1/h 2 ) to O(1/h). In future studies we will show that polygonal meshes can unlock new efficient adaptive refinement strategies.
We have introduced novel theory to address the challenges of the error analysis. First, if the exact solution (u, v) is H 2+1/4 (Ω) × H 2 (Γ), the lifting operator can be replaced with the Sobolev extension operator. Second, if Γ is sufficiently smooth, there exists an L 2 -preserving inverse trace operator that preserves the L 2 norm and the H 1 seminorm up to the same scale factor. Third, we have introduced a tailor-made Ritz projection in the bulk for the VEM that accounts for boundary conditions and that fulfils optimal error bounds. By using our bulk-Ritz projection we have drawn two consequences. First, the lowest order bulk-VEM [5] retains optimal convergence even in the simultaneous presence of curved boundaries and non-zero boundary conditions, a result that lacked a rigorous proof in the literature. Second, the proposed BSVEM possesses optimal spatial convergence, that is Numerical examples validate our findings in terms of (i) convergence rate in space and time for both the elliptic and the parabolic case and (ii) the computational advantages of polygonal meshes. Given the encouraging results, the extension of BSVEM to higher space dimension and different BSPDE problems is under development. Moreover, the generalisation of problem (2) to evolving BS domains, which would comprise additional models addressed in the literature such as [35], is left for future work.
Lemma A1 (Fermi coordinates). If Γ is a C k , k ≥ 2 surface, there exists an open neighbourhood U ⊂ R 3 of Γ such that every x ∈ U admits a unique decomposition of the form The maximal open set U with this property is called the Fermi stripe of Γ (see Fig. 1(a)), a(x) is called the normal projection onto Γ and (a(x), d(x)) are called the Fermi coordinates of x. The oriented distance function fulfils d ∈ C k (U ).
Definition A1 (C 1 (Γ) functions). A function u : Γ → R is said to be C 1 (Γ) if there exist an open neighbourhood U of Γ and a C 1 functionû : U → R such thatû |Γ = u, i.e.û is a C 1 extension of u off Γ.
Definition A2 (Tangential gradient and tangential derivatives). The tangential gradient ∇ Γ u of a function u ∈ C 1 (Γ) is defined by The result of the computation in (A.2) is independent of the choice of the extensionû. The components D x u and D y u of the tangential gradient ∇ Γ u are called the tangential derivatives of u.
Definition A3 (C k (Γ) functions). For k ∈ N, k > 1, a function u : Γ → R is said to be C k (Γ) if it is C 1 (Γ) and its tangential derivatives are C k−1 (Γ).
Since d ∈ C k (U ) from Lemma A1 and ν ∈ C k−1 (Γ) from (A.1), then F ∈ C k−1 (Ω × Γ). Since the normal projection a(x) of x ∈ U can be regarded as the solution of the implicit equation F (x, a(x)) = 0, then a ∈ C k−1 (U ) as well.
Definition A4 (Laplace-Beltrami operator). The Laplace-Beltrami operator ∆ Γ u of a function u ∈ C 2 (Γ) is defined by A2. Bulk-and surface function spaces  The following definition can be found in [25].
Definition A7 (Fractional Sobolev function spaces). Let s ∈ (0, 1) and p ∈ [1, +∞). For u : Ω → R, the bulk-fractional Sobolev norm is defined by with · being the Euclidean norm. If s > 1, s / ∈ N, by decomposing s as s = m + σ, where m ∈ N and σ ∈ (0, 1), the fractional Sobolev norm is defined as For any (integer or non-integer) s ∈ [0, +∞), the Sobolev-Slobodeckij space W s,p (Ω) is defined as where the norm · W s,p (Ω) is understood as (A.3) or (A.6) according to whether s is an integer or a non-integer number.
Theorem A2 (Narrow band trace inequality). With the notations of the previous theorem, there exists C > 0 depending on Ω such that any u ∈ H 1 (Ω) fulfils Proof. See [27].
Theorem A3 (Trace theorem and inverse trace theorem). Let k ∈ N, 1 2 < s ≤ k and assume that the boundary Γ is a C k curve. 1 Then there exists a bounded operator Tr : H s (Ω) → H s− 1 2 (Γ), called the trace operator, such that Tr(u) = u |Γ . The trace operator fulfils The trace operator has a continuous inverse operator Tr −1 : Proof. See [43] or [44].
A simple consequence of Theorem A3 is the following Corollary A1 (Normal trace theorem). Let k ∈ N \ {1}, 3 2 < s ≤ k and assume that the boundary Γ is a C k curve. There exists a bounded operator T n : H s (Ω) → L 2 (Γ), called the normal trace operator, such that T n(u) = ∇u · n, with n being the outward unit normal vector field on Γ. The normal trace fulfils T n(u) Proof. Let u ∈ H s (Ω), then ∇u ∈ H s−1 (Ω). From Theorem A3 we have that Tr(∇u) ∈ H s−3/2 (Γ).
If · denotes the Euclidean norm in R 2 , we have that Theorem A4 (Sobolev embeddings). Assume that Ω has a Lipschitz boundary.
Proof. See [2] for the case of integer-order Sobolev spaces and [25] for the fractional case.

Appendix B: proofs of main theorems
Proof of Theorem 5 By assumption we have that Tr(u) ∈ H 2 (Γ). Since Γ is a C 3 curve, then for all F ∈ F h , the normal projection a : F → a(F ) is a C 2 function, see Remark A1. Hence, (u −ℓ ) |F ∈ H 2 (F ) and From (25), (54), the second equation in (78), and (B.1) we have that with the terms in H 2+1/4 (Ω) norm arising only in the simultaneous presence of curved boundaries and non-triangular boundary elements, which proves (87).

(B.25)
To estimate the fourth term we proceed as in (B.24). By using (30), (54)-(55), (B.1), (B.18), the C 2 regularity and the Lipschitz continuity of r, we have that (B.26) We estimate the fifth term as in (B.26) by also using (A.11):  By applying Grönwall's lemma and accounting for the h 2 -accuracy of the initial conditions (65), we get the desired estimate.