<figure><img src="image_1.png" /></figure>

HAL
open science

# A priori error analysis of a virtual element approximation-based boundary element method for the electric field integral equation

Alexis Touzalin, Emanuele Arcese, Sébastien Pernet

► To cite this version:

Alexis Touzalin, Emanuele Arcese, Sébastien Pernet. A priori error analysis of a virtual element approximation-based boundary element method for the electric field integral equation. 2024. hal-04529074

HAL Id: hal-04529074

https://hal.science/hal-04529074

Preprint submitted on 2 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_1.png" /></figure>

HAL
open science

# A priori error analysis of a virtual element approximation-based boundary element method for the electric field integral equation

Alexis Touzalin, Emanuele Arcese, Sébastien Pernet

► To cite this version:

Alexis Touzalin, Emanuele Arcese, Sébastien Pernet. A priori error analysis of a virtual element approximation-based boundary element method for the electric field integral equation. 2024. hal-04529074

HAL Id: hal-04529074

https://hal.science/hal-04529074v1

Preprint submitted on 2 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# A priori error analysis of a virtual element approximation-based boundary element method for the electric field integral equation

ALEXIS TOUZALIN* AND EMANUELE ARCESE

CEA CESTA, 15 Avenue des Sablières, 33114, Le Barp, France

*Corresponding author: alexis.touzalin@cea.fr

AND

SÉBASTIEN PERNET

DTIS, ONERA, Université de Toulouse, 2 Avenue Edouard Belin, 31000, Toulouse, France

We develop a Galerkin discretization of the electric field integral equation based on a lowest-order virtual element approximation for surface meshes composed of polygons. This new boundary element discretization relies on the divergence-conforming virtual element counterparts of the classical Raviart-Thomas finite elements on simplices and is particularly suited to handle hanging nodes. We prove the well-posedness of the resulting stabilization-free numerical scheme by establishing the stable-uniform character of the discrete inf-sup condition in the natural norm for polyhedral surfaces. Moreover, we demonstrate through an *a priori* error analysis the quasi-optimal convergence of the scheme, leading to the same convergence rate as that of the classical Raviart-Thomas boundary element scheme. Finally, numerical experiments involving scattering problems are presented in order to give more insight into the behavior of the virtual boundary element scheme in terms of *h*-convergence and accuracy as a function of the regularity of solutions and meshes.

**Keywords:** Virtual elements; Maxwell's equations; Boundary integral operators; Polygonal meshes.

## 1. Introduction

With significant advances in matrix compression and preconditioning techniques, the boundary element method (BEM), based on boundary integral equations, has become a competitive tool to numerically simulate the scattering or radiation of time-harmonic electromagnetic waves. Applications can be found, for example, in antenna design (Francavilla et al., 2012; Stutzman and Thiele, 2012; Vipiana et al., 2010), in Radar Cross Section (RCS) analysis (Collino et al., 2008; Kong and Sheng, 2018; Peng et al., 2016; Stupfel, 2015) and in ElectroMagnetic Compatibility (EMC) (Li et al., 2014; Solis et al., 2020). An integral equation widely used in engineering to describe such physical phenomena is the so-known Electric Field Integral Equation (EFIE). A common discretization of the EFIE on the surface of the physical body by means of Galerkin methods relies on Raviart-Thomas (RT) divergence-conforming finite elements (Raviart and Thomas, 1977) (or also known as the Rao-Wilton-Glisson finite elements, RWG (Rao et al., 1982) in short). The first comprehensive study, dating back to the 80s, of this conforming approximation for any finite element space order on smooth surfaces was done in Bendali (1984a,b). In the beginning of the 2000s, various efforts were devoted to extending the theoretical understanding of the convergence behavior of these Galerkin schemes (also involving other divergence-conforming spaces, such as the Brezzi-Douglas-Marini ones, BDM (Brezzi et al., 1985) in short) to general surfaces, including the Lipschitz ones with boundary (Buffa and Christiansen, 2003) as well as the polyhedral ones with either planar faces (Buffa et al., 2002a; Hiptmair and Schwab, 2002) or curvilinear faces (Buffa and Hiptmair, 2003). However, all these studies are confined to the standard
# A priori error analysis of a virtual element approximation-based boundary element method for the electric field integral equation

ALEXIS TOUZALIN* AND EMANUELE ARCESE

CEA CESTA, 15 Avenue des Sablières, 33114, Le Barp, France

*Corresponding author: alexis.touzalin@cea.fr

AND

SÉBASTIEN PERNET

DTIS, ONERA, Université de Toulouse, 2 Avenue Edouard Belin, 31000, Toulouse, France

We develop a Galerkin discretization of the electric field integral equation based on a lowest-order virtual element approximation for surface meshes composed of polygons. This new boundary element discretization relies on the divergence-conforming virtual element counterparts of the classical Raviart-Thomas finite elements on simplices and is particularly suited to handle hanging nodes. We prove the well-posedness of the resulting stabilization-free numerical scheme by establishing the stable-uniform character of the discrete inf-sup condition in the natural norm for polyhedral surfaces. Moreover, we demonstrate through an *a priori* error analysis the quasi-optimal convergence of the scheme, leading to the same convergence rate as that of the classical Raviart-Thomas boundary element scheme. Finally, numerical experiments involving scattering problems are presented in order to give more insight into the behavior of the virtual boundary element scheme in terms of *h*-convergence and accuracy as a function of the regularity of solutions and meshes.

**Keywords:** Virtual elements; Maxwell's equations; Boundary integral operators; Polygonal meshes.

## 1. Introduction

With significant advances in matrix compression and preconditioning techniques, the boundary element method (BEM), based on boundary integral equations, has become a competitive tool to numerically simulate the scattering or radiation of time-harmonic electromagnetic waves. Applications can be found, for example, in antenna design (Francavilla et al., 2012; Stutzman and Thiele, 2012; Vipiana et al., 2010), in Radar Cross Section (RCS) analysis (Collino et al., 2008; Kong and Sheng, 2018; Peng et al., 2016; Stupfel, 2015) and in ElectroMagnetic Compatibility (EMC) (Li et al., 2014; Solis et al., 2020). An integral equation widely used in engineering to describe such physical phenomena is the so-known Electric Field Integral Equation (EFIE). A common discretization of the EFIE on the surface of the physical body by means of Galerkin methods relies on Raviart-Thomas (RT) divergence-conforming finite elements (Raviart and Thomas, 1977) (or also known as the Rao-Wilton-Glisson finite elements, RWG (Rao et al., 1982) in short). The first comprehensive study, dating back to the 80s, of this conforming approximation for any finite element space order on smooth surfaces was done in Bendali (1984a,b). In the beginning of the 2000s, various efforts were devoted to extending the theoretical understanding of the convergence behavior of these Galerkin schemes (also involving other divergence-conforming spaces, such as the Brezzi-Douglas-Marini ones, BDM (Brezzi et al., 1985) in short) to general surfaces, including the Lipschitz ones with boundary (Buffa and Christiansen, 2003) as well as the polyhedral ones with either planar faces (Buffa et al., 2002a; Hiptmair and Schwab, 2002) or curvilinear faces (Buffa and Hiptmair, 2003). However, all these studies are confined to the standard
finite element setting which prohibits the presence of hanging nodes within the mesh of the body surface. This means that mesh elements can intersect either in a common vertex or along a common side which, in turn, cannot contain any additional nodes that are not element vertices. These hanging nodes may appear in many situations, as for instance during local refinement/coarsening of the mesh within design and optimization processes. Nevertheless, in real-life electromagnetic problems it can be convenient to handle a nonmatching mesh, which can be built by independently partitioning each part of the domain into elements whose size is adapted to capture the local electromagnetic features. This would avoid generating an oversized conventional mesh (without hanging nodes) with eventually narrow elements, leading to the solution of an overly large and potentially ill-conditioned linear system, originating from the EFIE discretization (the reader can refer to Adrian et al. (2021) for a recent review and discussion on discretization-related ill-conditioning sources for the EFIE).

To overcome this defect, nonconforming approximations such as the Discontinuous Galerkin (DG) methods, which enforce only weakly the solution continuity across the element's boundaries, have been successfully devised in the partial differential equations (PDE) community (Cockburn et al., 2000; Hesthaven and Warburton, 2008; Houston et al., 2005; Melenk et al., 2013). Some works about DG approximations of boundary integral equations for (quasi-)elliptic scalar problems have been proposed in (Heuer and Karkulik, 2017; Heuer and Meddahi, 2013; Heuer and Salmerón, 2017). In the realm of electromagnetism, a first extension of the DG formalism to boundary integral equations, resulting in an interior penalty variational formulation, has recently been studied in Peng et al. (2013). Compared to the conventional divergence-conforming approximation of the EFIE, this DG formulation involves two additional terms, namely a consistency-type term and an interior penalty stabilization term. In a subsequent study (Peng et al., 2016), the authors have further investigated these typical terms coming up with an antisymmetric variant of the previous DG formulation. While good accuracy has been observed through the several numerical experiments performed on both academic and more general geometries (see, e.g., (Echeverri Bautista et al., 2015; Peng et al., 2016, 2013) and references therein), these DG schemes present two main drawbacks. First, the practical interest of employing discontinuous polynomial basis functions within DG spaces results in doubling the unknowns at the boundaries of each mesh element, which leads to an increase in the linear system size. Afterwards, although the stabilization term is necessary for the well-posedness of the discrete problem, this term can negatively impact the conditioning of the linear system and the convergence of the numerical method. And the optimal choice of such a term cannot be obtained through theoretical analysis. For example, in the nonsymmetric formulation (Peng et al., 2016), authors have chosen a stabilization term on the basis of numerical experiments that is proportional to $\log(h)$, with $h$ denoting the mesh-size, in order to ensure the scheme's optimal convergence. Moreover, as far as we known, a complete mathematical analysis of this DG discretization of the EFIE is still missing. A first step toward a theoretical analysis of the interior penalty DG formulations of boundary integral equation for the Helmholtz problem was made in (Heuer and Salmerón, 2017; Messai and Pernet, 2020). Finally, there also are other Galerkin discretization schemes for the boundary integral equations that handle simplicial surface meshes with hanging nodes worth mentioning (Bendali et al., 2012; Martin et al., 2023).

Recently, a new conforming approximation paradigm called the virtual element method (VEM), firstly introduced in Beirão da Veiga et al. (2013), has considerably attracted the engineering community's attention for solving PDE. The VEM owns a great success as, being a generalization of the standard finite element variational setting, the method supports very general meshes with elements having an almost arbitrary shape; in this framework indeed, any hanging node is simply interpreted as a vertex of degenerate elements (say, e.g., a triangle with five nodes on its boundary is considered as a polygon with five edges). Roughly speaking, the main idea behind the VEM is that, on each
finite element setting which prohibits the presence of hanging nodes within the mesh of the body surface. This means that mesh elements can intersect either in a common vertex or along a common side which, in turn, cannot contain any additional nodes that are not element vertices. These hanging nodes may appear in many situations, as for instance during local refinement/coarsening of the mesh within design and optimization processes. Nevertheless, in real-life electromagnetic problems it can be convenient to handle a nonmatching mesh, which can be built by independently partitioning each part of the domain into elements whose size is adapted to capture the local electromagnetic features. This would avoid generating an oversized conventional mesh (without hanging nodes) with eventually narrow elements, leading to the solution of an overly large and potentially ill-conditioned linear system, originating from the EFIE discretization (the reader can refer to Adrian et al. (2021) for a recent review and discussion on discretization-related ill-conditioning sources for the EFIE).

To overcome this defect, nonconforming approximations such as the Discontinuous Galerkin (DG) methods, which enforce only weakly the solution continuity across the element's boundaries, have been successfully devised in the partial differential equations (PDE) community (Cockburn et al., 2000; Hesthaven and Warburton, 2008; Houston et al., 2005; Melenk et al., 2013). Some works about DG approximations of boundary integral equations for (quasi-)elliptic scalar problems have been proposed in (Heuer and Karkulik, 2017; Heuer and Meddahi, 2013; Heuer and Salmerón, 2017). In the realm of electromagnetism, a first extension of the DG formalism to boundary integral equations, resulting in an interior penalty variational formulation, has recently been studied in Peng et al. (2013). Compared to the conventional divergence-conforming approximation of the EFIE, this DG formulation involves two additional terms, namely a consistency-type term and an interior penalty stabilization term. In a subsequent study (Peng et al., 2016), the authors have further investigated these typical terms coming up with an antisymmetric variant of the previous DG formulation. While good accuracy has been observed through the several numerical experiments performed on both academic and more general geometries (see, e.g., (Echeverri Bautista et al., 2015; Peng et al., 2016, 2013) and references therein), these DG schemes present two main drawbacks. First, the practical interest of employing discontinuous polynomial basis functions within DG spaces results in doubling the unknowns at the boundaries of each mesh element, which leads to an increase in the linear system size. Afterwards, although the stabilization term is necessary for the well-posedness of the discrete problem, this term can negatively impact the conditioning of the linear system and the convergence of the numerical method. And the optimal choice of such a term cannot be obtained through theoretical analysis. For example, in the nonsymmetric formulation (Peng et al., 2016), authors have chosen a stabilization term on the basis of numerical experiments that is proportional to $\log(h)$, with $h$ denoting the mesh-size, in order to ensure the scheme's optimal convergence. Moreover, as far as we known, a complete mathematical analysis of this DG discretization of the EFIE is still missing. A first step toward a theoretical analysis of the interior penalty DG formulations of boundary integral equation for the Helmholtz problem was made in (Heuer and Salmerón, 2017; Messai and Pernet, 2020). Finally, there also are other Galerkin discretization schemes for the boundary integral equations that handle simplicial surface meshes with hanging nodes worth mentioning (Bendali et al., 2012; Martin et al., 2023).

Recently, a new conforming approximation paradigm called the virtual element method (VEM), firstly introduced in Beirão da Veiga et al. (2013), has considerably attracted the engineering community's attention for solving PDE. The VEM owns a great success as, being a generalization of the standard finite element variational setting, the method supports very general meshes with elements having an almost arbitrary shape; in this framework indeed, any hanging node is simply interpreted as a vertex of degenerate elements (say, e.g., a triangle with five nodes on its boundary is considered as a polygon with five edges). Roughly speaking, the main idea behind the VEM is that, on each
mesh element, the approximation space is spanned by virtual basis functions that are not explicitly
known. These functions consist of polynomials up to a certain degree ensuring the desired accuracy
plus other general functions being solution to local PDE. The use of some local projection operators,
only defined through the knowledge of the degrees of freedom, associated to the virtual basis functions,
allows to compute the discrete bilinear forms with the required accuracy and hence, to build the whole
algebraic problem up to stabilization. An explicit stabilization term is then generally included in the
weak formulation to make the underlying discrete scalar product equivalent to the continuous one.
Many contributions are available in the literature (Beirão da Veiga et al., 2013, 2018, 2017, 2022;
Beirão da Veiga and Mascotto, 2022) including the extension of VEM to curvilinear mesh elements
(Beirão da Veiga et al., 2019; Dassi et al., 2022). In electromagnetism, few works exist (Barnafi et al.,
2023; Beirão da Veiga et al., 2018; Beirão da Veiga et al., 2018; Beirão da Veiga et al., 2022) (these
lists are not exhaustive).

In this work, we propose for the first time to extend the principle of the VEM to the approximation
of boundary integral operators used in the EFIE. This aims at providing the new boundary element
discretization, which involves the natural virtual element counterparts of the classical lowest-order RT
elements, with robustness in respect of the hanging nodes and hence, with flexibility in handling more
general mesh elements than simplices. A feature that can help address the structural and geometrical
complexities of the scattering or radiating bodies. Our purpose in the present work is to theoretically
study and numerically analyze a novel Galerkin boundary element method based on a virtual element
approximation space (termed here V-BEM, short for virtual boundary element method) for discretizing
the EFIE on nonsmooth surfaces. We particularly focus on the scattering of electromagnetic plane
waves by a perfectly conducting body, which is assumed, in this paper, to be a polyhedron whose
boundary consists of a finite number of planar faces. This analysis may be clearly extended to both
curvilinear piecewise polyhedral surfaces, by employing particular Lipschitz-isomorphisms from the
approximated surface to the physical one (see, e.g., Section 3.5 of Christiansen (2003)), and open
surfaces, by considering the suitable functional framework (Buffa and Christiansen, 2003).

The strongly non-local variational setting of the EFIE (involving, e.g., pseudo-differential operators, non-local norms) does not allow the direct application of the standard VEM approach for discretization, which would require the introduction of a local projector defined with respect to the operators used in the boundary integral formulations. A natural solution we opt for to address this problem is to go a little beyond the VEM framework by simply considering a $L^2$-projector from the literature. This however results in a perturbed discrete weak formulation of the EFIE, whose analysis (i.e. well-posedness and convergence) cannot be carried out on the basis of existing theoretical results in this field (see, e.g., (Buffa and Hiptmair, 2003; Christiansen, 2003; Hiptmair and Schwab, 2002)). We devise the new EFIE discretization by making use of two key components, coming both from the VEM literature: the divergence-conforming lowest-order virtual element space (the serendipity variant) (Beirão da Veiga and Mascotto, 2022), along with the element-wise $L^2$-orthogonal projection (Beirão da Veiga et al., 2018). We point out that the proposed weak formulation gets rid of the typical stabilization term of the VEM. Our choice for this approximation setting is motivated by the fact that, according to Beirão da Veiga et al. (2017), such a virtual space allows for the extension of the lowest-order RT space on polygonal elements. And as a result, the proposed V-BEM can be interpreted as a discretization scheme combining the virtual element approximation on (degenerate) polygons with the classical RT approximation whenever the elements are triangles.

The main contribution of this work is twofold: on the one hand, we prove a uniform discrete *inf-sup*
condition of the new stabilization-free V-BEM formulation and establish an *a priori* error estimate for
the resulting scheme in the natural norm, i.e. the $H^{-1/2}(\operatorname{div}\Gamma, \Gamma)$-norm. On the other hand, we provide
mesh element, the approximation space is spanned by virtual basis functions that are not explicitly known. These functions consist of polynomials up to a certain degree ensuring the desired accuracy plus other general functions being solution to local PDE. The use of some local projection operators, only defined through the knowledge of the degrees of freedom, associated to the virtual basis functions, allows to compute the discrete bilinear forms with the required accuracy and hence, to build the whole algebraic problem up to stabilization. An explicit stabilization term is then generally included in the weak formulation to make the underlying discrete scalar product equivalent to the continuous one. Many contributions are available in the literature (Beirão da Veiga et al., 2013, 2018, 2017, 2022; Beirão da Veiga and Mascotto, 2022) including the extension of VEM to curvilinear mesh elements (Beirão da Veiga et al., 2019; Dassi et al., 2022). In electromagnetism, few works exist (Barnafi et al., 2023; Beirão da Veiga et al., 2018; Beirão da Veiga et al., 2018; Beirão da Veiga et al., 2022) (these lists are not exhaustive).

In this work, we propose for the first time to extend the principle of the VEM to the approximation of boundary integral operators used in the EFIE. This aims at providing the new boundary element discretization, which involves the natural virtual element counterparts of the classical lowest-order RT elements, with robustness in respect of the hanging nodes and hence, with flexibility in handling more general mesh elements than simplices. A feature that can help address the structural and geometrical complexities of the scattering or radiating bodies. Our purpose in the present work is to theoretically study and numerically analyze a novel Galerkin boundary element method based on a virtual element approximation space (termed here V-BEM, short for virtual boundary element method) for discretizing the EFIE on nonsmooth surfaces. We particularly focus on the scattering of electromagnetic plane waves by a perfectly conducting body, which is assumed, in this paper, to be a polyhedron whose boundary consists of a finite number of planar faces. This analysis may be clearly extended to both curvilinear piecewise polyhedral surfaces, by employing particular Lipschitz-isomorphisms from the approximated surface to the physical one (see, e.g., Section 3.5 of Christiansen (2003)), and open surfaces, by considering the suitable functional framework (Buffa and Christiansen, 2003).

The strongly non-local variational setting of the EFIE (involving, e.g., pseudo-differential operators, non-local norms) does not allow the direct application of the standard VEM approach for discretization, which would require the introduction of a local projector defined with respect to the operators used in the boundary integral formulations. A natural solution we opt for to address this problem is to go a little beyond the VEM framework by simply considering a $L^2$-projector from the literature. This however results in a perturbed discrete weak formulation of the EFIE, whose analysis (i.e. well-posedness and convergence) cannot be carried out on the basis of existing theoretical results in this field (see, e.g., (Buffa and Hiptmair, 2003; Christiansen, 2003; Hiptmair and Schwab, 2002)). We devise the new EFIE discretization by making use of two key components, coming both from the VEM literature: the divergence-conforming lowest-order virtual element space (the serendipity variant) (Beirão da Veiga and Mascotto, 2022), along with the element-wise $L^2$-orthogonal projection (Beirão da Veiga et al., 2018). We point out that the proposed weak formulation gets rid of the typical stabilization term of the VEM. Our choice for this approximation setting is motivated by the fact that, according to Beirão da Veiga et al. (2017), such a virtual space allows for the extension of the lowest-order RT space on polygonal elements. And as a result, the proposed V-BEM can be interpreted as a discretization scheme combining the virtual element approximation on (degenerate) polygons with the classical RT approximation whenever the elements are triangles.

The main contribution of this work is twofold: on the one hand, we prove a uniform discrete *inf-sup* condition of the new stabilization-free V-BEM formulation and establish an *a priori* error estimate for the resulting scheme in the natural norm, i.e. the $H^{-1/2}(\operatorname{div}\Gamma, \Gamma)$-norm. On the other hand, we provide
a broader understanding of the scheme behavior through a numerical analysis of the h-convergence on surface currents and RCS for both smooth and nonsmooth scattering bodies.

The outline of the paper is the following. In section 2, we recall some notations and important functional spaces, as well as the mesh regularity assumptions. The EFIE model problem is detailed in section 3, wherein we introduce the new Galerkin discretization by the divergence-conforming virtual elements. In section 4, we summarize the main theoretical results of the paper, such as the h-uniform discrete inf-sup condition and the a priori error estimate. Section 5 gathers all the proofs and technical tools needed to prove the results of the previous section. Section 6 is dedicated to numerical experiments aiming at comparing the h-convergence of the proposed V-BEM with the classical BEM approach based on the RT approximation space. Finally, we draw some conclusions in section 7.

Let us point out, from the very beginning, that to keep the paper understandable we only investigate lowest-order RT and virtual approximation spaces on planar elements, any generalization to higher order spaces being straightforward. Furthermore, in the present paper, we do not describe any issues related to either the implementation or the solution of the induced linear system. All these questions will be the subject of a future article.

Nota Bene: In the rest of the article, we denote $a \lesssim b$ (or $a \gtrsim b$) if $\exists C > 0$, such that $a \le Cb$ ($a \ge Cb$, respectively), where $a$ and $b$ are comparable objects and $C$ is independent of discretization, but can vary with the problem parameters (body geometry, wavenumber, etc.).

## 2. Notation and functional spaces

Let us recall standard definitions and main properties for Hilbert and Sobolev spaces which allow to correctly define trace and integral operators in the context of polyhedral domains. Some basic assumptions on the surface meshes used for partitioning the domain are stated as well. To begin with, let $\Omega$ be a closed polyhedron in $\mathbb{R}^3$ whose boundary $\Gamma := \partial\Omega$ is assumed to be split into $N_\Gamma$ planar faces $(\Gamma_i)_{i=1,...,N_\Gamma}$, and let $n \in L^\infty(\Gamma)$ be the unit normal vector field on the boundary pointing outward of $\Omega$.

### 2.1. Functional spaces

At each of these faces $\Gamma_i$, we associate a 2D local Cartesian representation using a local orthonormal coordinate system $(O_i, \tau_1^i, \tau_2^i)$. A point $x \in \Gamma_i$ can be represented in two ways by using either its classical Cartesian variables $(x_1, x_2, x_3)$ or the local coordinate system, i.e. $x^{2D} = (x_1, x_2)$. We also define the orthogonal variable associated to $x^{2D}$ by $x^{2D,\perp} := (-x_2, x_1)$. In two dimensions, let $v$ be a vector and $q$ a scalar, we denote the divergence operator and the curl operators as

$$ \nabla \cdot v = \frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2}, \quad \mathrm{curl} v = \frac{\partial v_2}{\partial x_1} - \frac{\partial v_1}{\partial x_2}, \quad \mathrm{curl} q = \left( \frac{\partial q}{\partial x_2}, -\frac{\partial q}{\partial x_1} \right). \quad (2.1) $$

On the other hand, if $v$ is a three-dimensional vector field $(v_1, v_2, v_3)$, we consider

$$ \nabla \cdot v = \frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2} + \frac{\partial v_3}{\partial x_3}, \quad \mathrm{curl} v = \left( \frac{\partial v_3}{\partial x_2} - \frac{\partial v_2}{\partial x_3}, \frac{\partial v_1}{\partial x_3} - \frac{\partial v_3}{\partial x_1}, \frac{\partial v_2}{\partial x_1} - \frac{\partial v_1}{\partial x_2} \right). \quad (2.2) $$

In doing so, there will be no ambiguity in the paper when using 2D and 3D operators. It is now possible to recall the standard definitions of some Hilbert and Sobolev spaces. With $s \ge 0$, we make use of
a broader understanding of the scheme behavior through a numerical analysis of the h-convergence on surface currents and RCS for both smooth and nonsmooth scattering bodies.

The outline of the paper is the following. In section 2, we recall some notations and important functional spaces, as well as the mesh regularity assumptions. The EFIE model problem is detailed in section 3, wherein we introduce the new Galerkin discretization by the divergence-conforming virtual elements. In section 4, we summarize the main theoretical results of the paper, such as the h-uniform discrete inf-sup condition and the a priori error estimate. Section 5 gathers all the proofs and technical tools needed to prove the results of the previous section. Section 6 is dedicated to numerical experiments aiming at comparing the h-convergence of the proposed V-BEM with the classical BEM approach based on the RT approximation space. Finally, we draw some conclusions in section 7.

Let us point out, from the very beginning, that to keep the paper understandable we only investigate lowest-order RT and virtual approximation spaces on planar elements, any generalization to higher order spaces being straightforward. Furthermore, in the present paper, we do not describe any issues related to either the implementation or the solution of the induced linear system. All these questions will be the subject of a future article.

Nota Bene: In the rest of the article, we denote a ≲ b (or a ≳ b) if ∃C > 0, such that a ≤ Cb (a ≥ Cb, respectively), where a and b are comparable objects and C is independent of discretization, but can vary with the problem parameters (body geometry, wavenumber, etc.).

## 2. Notation and functional spaces

Let us recall standard definitions and main properties for Hilbert and Sobolev spaces which allow to correctly define trace and integral operators in the context of polyhedral domains. Some basic assumptions on the surface meshes used for partitioning the domain are stated as well. To begin with, let Ω be a closed polyhedron in R³ whose boundary Γ := ∂Ω is assumed to be split into NΓ planar faces (Γᵢ)ᵢ=₁,··,NΓ, and let n ∈ L∞(Γ) be the unit normal vector field on the boundary pointing outward of Ω.

### 2.1. Functional spaces

At each of these faces Γᵢ, we associate a 2D local Cartesian representation using a local orthonormal coordinate system (Oᵢ, τ₁ᵢ, τ₂ᵢ). A point x ∈ Γᵢ can be represented in two ways by using either its classical Cartesian variables (x₁, x₂, x₃) or the local coordinate system, i.e. x²D = (x₁, x₂). We also define the orthogonal variable associated to x²D by x²D,⊥ := (-x₂, x₁). In two dimensions, let v be a vector and q a scalar, we denote the divergence operator and the curl operators as

$$ \nabla \cdot v = \frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2}, \quad \mathrm{curl} v = \frac{\partial v_2}{\partial x_1} - \frac{\partial v_1}{\partial x_2}, \quad \mathrm{curl} q = \left( \frac{\partial q}{\partial x_2}, -\frac{\partial q}{\partial x_1} \right). \quad (2.1) $$

On the other hand, if v is a three-dimensional vector field (v₁, v₂, v₃), we consider

$$ \nabla \cdot v = \frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2} + \frac{\partial v_3}{\partial x_3}, \quad \mathrm{curl} v = \left( \frac{\partial v_3}{\partial x_2} - \frac{\partial v_2}{\partial x_3}, \frac{\partial v_1}{\partial x_3} - \frac{\partial v_3}{\partial x_1}, \frac{\partial v_2}{\partial x_1} - \frac{\partial v_1}{\partial x_2} \right). \quad (2.2) $$

In doing so, there will be no ambiguity in the paper when using 2D and 3D operators. It is now possible to recall the standard definitions of some Hilbert and Sobolev spaces. With s ≥ 0, we make use of
the Sobolev spaces $H_{\text{loc}}^s(\Omega)$ with $L^2 \equiv H^0 \equiv H$, where the loc sub-fix implies that the space lives on any compact of $\Omega$. Of course, Hilbert spaces $H^s(\Omega)$ appear when the sub-fix is dropped. These spaces are equipped with the norms $\|\cdot\|_{H^s(\Omega)}$. It is worth noting that spaces written here in bold letters denote spaces involving vectors, whereas spaces in normal letters involve scalars like the Hilbert space $H^s(\Omega)$. Furthermore, with $\mathbf{A}$ a first-order differential operator, we use the definitions of the spaces given in Buffa and Hiptmair (2003):

$$H_{\text{loc}}^{s}(\mathbf{A}, \Omega) = \{\mathbf{v} \in H_{\text{loc}}^{s}(\Omega) \mid \mathbf{A}\mathbf{v} \in H_{\text{loc}}^{s}(\Omega)\}, \qquad (2.3)$$

$$H^s(\mathbf{A}, \Omega) = \{\mathbf{v} \in H^s(\Omega) \mid \mathbf{A}\mathbf{v} \in H^s(\Omega)\}, \qquad (2.4)$$

$$H_{\text{loc}}^{s}(\mathbf{A}0, \Omega) = \{\mathbf{v} \in H_{\text{loc}}^{s}(\Omega) \mid \mathbf{A}\mathbf{v} = 0\}, \qquad (2.5)$$

$$H^s(\mathbf{A}0, \Omega) = \{\mathbf{v} \in H^s(\Omega) \mid \mathbf{A}\mathbf{v} = 0\}. \qquad (2.6)$$

Some trace operators are used, especially the tangential trace and the tangential component, which are defined for a sufficiently smooth function $\mathbf{v} : \Omega \to \mathbb{C}^3$, by $\gamma_t \mathbf{v} = \mathbf{n} \times \mathbf{v}_{|\Gamma}$ and $\pi_t \mathbf{v} = \gamma_t \mathbf{v} \times \mathbf{n}$, respectively. In particular, these operators are well-defined for all functions in the spaces $H^{s+1/2}(\Omega)$ for $0 < s < 1$ due to the continuity results of the standard trace operator in Lipschitz domain (see theorem 3.38 in McLean and McLean (2000)). Consequently, for $s \in (0, 1)$, applying the tangential trace or the tangential component operators on Hilbert spaces allows us to get the following tangential spaces:

$$H_{\perp}^{s}(\Gamma) := \gamma_t (H^{s+1/2}(\Omega)), \qquad H_{\|}^{s}(\Gamma) := \pi_t (H^{s+1/2}(\Omega)), \qquad (2.7)$$

where $H_{\perp}^s(\Gamma)$, $H_{\|}^s(\Gamma) \subset L_t^2(\Gamma)$, with $L_t^2(\Gamma)$ being the tangential space defined as $\{\mathbf{v} \in (L^2(\Gamma))^3 \mid \mathbf{u} \cdot \mathbf{n} = 0\}$. The dual spaces of $H_{\perp}^s(\Gamma)$ and $H_{\|}^s(\Gamma)$ are $H_{\perp}^{-s}(\Gamma)$ and $H_{\|}^{-s}(\Gamma)$ with respect to the $L_t^2(\Gamma)$ pivot space for $s \in (0, 1)$, respectively, and the related duality pairings are denoted $\langle \cdot, \cdot \rangle_{\perp, s, -s}$ and $\langle \cdot, \cdot \rangle_{\|, s, -s}$ (Hiptmair and Schwab, 2002). In order to make the notation clearer in throughout the paper, we will improperly shorten the sub-fix of these duality pairings for $s = 1/2$, which, for instance, for the vector and scalar cases read

$$\text{if } a \in H^{1/2}(\Gamma), b \in H^{-1/2}(\Gamma) \text{ then } \langle a, b \rangle_{1/2, -1/2} \text{ is written as } \langle a, b \rangle, \qquad (2.8)$$

$$\text{if } \mathbf{a} \in \mathbf{H}_{\|}^{1/2}(\Gamma), \mathbf{b} \in \mathbf{H}_{\|}^{-1/2}(\Gamma) \text{ then } \langle \mathbf{a}, \mathbf{b} \rangle_{\|, 1/2, -1/2} \text{ is written as } \langle \mathbf{a}, \mathbf{b} \rangle_{\|}. \qquad (2.9)$$

Whenever ambiguity arises, the full notation of the duality pairings is used. Since $\Omega$ is assumed to be a polyhedron, the tangential gradient can be defined face by face: let $\mathbf{v} : \Omega \to \mathbb{C}$ be a sufficiently smooth function,

$$\nabla_{\Gamma} \mathbf{v}(\mathbf{x}) = \nabla_{\Gamma_i} \mathbf{v}(\mathbf{x}) = \pi_t (\nabla \mathbf{v}) (\mathbf{x}), \quad \text{if } \mathbf{x} \in \Gamma_i. \qquad (2.10)$$

By using the local orthonormal coordinate system $(O_i, \tau_1^i, \tau_2^i)$, $\nabla_{\Gamma_i}$ can be equivalently defined as 2D gradient applied to the function $v^{2D}(x_1, x_2) := v(O_i + x_1\tau_1^i + x_2\tau_2^i)$, i.e.

$$\nabla_{\Gamma_i} v = \nabla v^{2D} = \left( \frac{\partial v^{2D}}{\partial x_1}, \frac{\partial v^{2D}}{\partial x_2} \right)^T. \qquad (2.11)$$

In the same way, the vectorial tangential curl operator is defined as

$$\mathrm{curl}_{\Gamma} v(\mathbf{x}) = \mathrm{curl}_{\Gamma_i} v(\mathbf{x}) = \gamma_t (\nabla v)(\mathbf{x}) = \mathrm{curl} v^{2D}(\mathbf{x}^{2D}), \quad \text{if } \mathbf{x} \in \Gamma_i. \qquad (2.12)$$
the Sobolev spaces $H_{\text{loc}}^s(\Omega)$ with $L^2 \equiv H^0 \equiv H$, where the loc sub-fix implies that the space lives on any compact of $\Omega$. Of course, Hilbert spaces $H^s(\Omega)$ appear when the sub-fix is dropped. These spaces are equipped with the norms $\|\cdot\|_{H^s(\Omega)}$. It is worth noting that spaces written here in bold letters denote spaces involving vectors, whereas spaces in normal letters involve scalars like the Hilbert space $H^s(\Omega)$. Furthermore, with $\mathbf{A}$ a first-order differential operator, we use the definitions of the spaces given in Buffa and Hiptmair (2003):

$$H_{\text{loc}}^{s}(\mathbf{A}, \Omega) = \{\mathbf{v} \in H_{\text{loc}}^{s}(\Omega) \mid \mathbf{A}\mathbf{v} \in H_{\text{loc}}^{s}(\Omega)\}, \qquad (2.3)$$

$$H^s(\mathbf{A}, \Omega) = \{\mathbf{v} \in H^s(\Omega) \mid \mathbf{A}\mathbf{v} \in H^s(\Omega)\}, \qquad (2.4)$$

$$H_{\text{loc}}^{s}(\mathbf{A}0, \Omega) = \{\mathbf{v} \in H_{\text{loc}}^{s}(\Omega) \mid \mathbf{A}\mathbf{v} = 0\}, \qquad (2.5)$$

$$H^s(\mathbf{A}0, \Omega) = \{\mathbf{v} \in H^s(\Omega) \mid \mathbf{A}\mathbf{v} = 0\}. \qquad (2.6)$$

Some trace operators are used, especially the tangential trace and the tangential component, which are defined for a sufficiently smooth function $\mathbf{v} : \Omega \to \mathbb{C}^3$, by $\gamma_t \mathbf{v} = \mathbf{n} \times \mathbf{v}_{|\Gamma}$ and $\pi_t \mathbf{v} = \gamma_t \mathbf{v} \times \mathbf{n}$, respectively. In particular, these operators are well-defined for all functions in the spaces $H^{s+1/2}(\Omega)$ for $0 < s < 1$ due to the continuity results of the standard trace operator in Lipschitz domain (see theorem 3.38 in McLean and McLean (2000)). Consequently, for $s \in (0, 1)$, applying the tangential trace or the tangential component operators on Hilbert spaces allows us to get the following tangential spaces:

$$H_{\perp}^{s}(\Gamma) := \gamma_t (H^{s+1/2}(\Omega)), \qquad H_{\|}^{s}(\Gamma) := \pi_t (H^{s+1/2}(\Omega)), \qquad (2.7)$$

where $H_{\perp}^s(\Gamma)$, $H_{\|}^s(\Gamma) \subset L_t^2(\Gamma)$, with $L_t^2(\Gamma)$ being the tangential space defined as $\{\mathbf{v} \in (L^2(\Gamma))^3 \mid \mathbf{u} \cdot \mathbf{n} = 0\}$. The dual spaces of $H_{\perp}^s(\Gamma)$ and $H_{\|}^s(\Gamma)$ are $H_{\perp}^{-s}(\Gamma)$ and $H_{\|}^{-s}(\Gamma)$ with respect to the $L_t^2(\Gamma)$ pivot space for $s \in (0, 1)$, respectively, and the related duality pairings are denoted $\langle \cdot, \cdot \rangle_{\perp, s, -s}$ and $\langle \cdot, \cdot \rangle_{\|, s, -s}$ (Hiptmair and Schwab, 2002). In order to make the notation clearer in throughout the paper, we will improperly shorten the sub-fix of these duality pairings for $s = 1/2$, which, for instance, for the vector and scalar cases read

$$\text{if } a \in H^{1/2}(\Gamma), b \in H^{-1/2}(\Gamma) \text{ then } \langle a, b \rangle_{1/2, -1/2} \text{ is written as } \langle a, b \rangle, \qquad (2.8)$$

$$\text{if } \mathbf{a} \in \mathbf{H}_{\|}^{1/2}(\Gamma), \mathbf{b} \in \mathbf{H}_{\|}^{-1/2}(\Gamma) \text{ then } \langle \mathbf{a}, \mathbf{b} \rangle_{\|, 1/2, -1/2} \text{ is written as } \langle \mathbf{a}, \mathbf{b} \rangle_{\|}. \qquad (2.9)$$

Whenever ambiguity arises, the full notation of the duality pairings is used. Since $\Omega$ is assumed to be a polyhedron, the tangential gradient can be defined face by face: let $\mathbf{v} : \Omega \to \mathbb{C}$ be a sufficiently smooth function,

$$\nabla_{\Gamma} \mathbf{v}(\mathbf{x}) = \nabla_{\Gamma_i} \mathbf{v}(\mathbf{x}) = \pi_t (\nabla \mathbf{v}) (\mathbf{x}), \quad \text{if } \mathbf{x} \in \Gamma_i. \qquad (2.10)$$

By using the local orthonormal coordinate system $(O_i, \tau_1^i, \tau_2^i)$, $\nabla_{\Gamma_i}$ can be equivalently defined as 2D gradient applied to the function $v^{2D}(x_1, x_2) := v(O_i + x_1\tau_1^i + x_2\tau_2^i)$, i.e.

$$\nabla_{\Gamma_i} v = \nabla v^{2D} = \left( \frac{\partial v^{2D}}{\partial x_1}, \frac{\partial v^{2D}}{\partial x_2} \right)^T. \qquad (2.11)$$

In the same way, the vectorial tangential curl operator is defined as

$$\mathrm{curl}_{\Gamma} v(\mathbf{x}) = \mathrm{curl}_{\Gamma_i} v(\mathbf{x}) = \gamma_t (\nabla v)(\mathbf{x}) = \mathrm{curl} v^{2D}(\mathbf{x}^{2D}), \quad \text{if } \mathbf{x} \in \Gamma_i. \qquad (2.12)$$
These two tangential operators can be extended as a linear continuous mapping from $H^{3/2}(\Gamma)$ to $H_{\parallel}^{1/2}(\Gamma)$ and $H_{\perp}^{1/2}(\Gamma)$, respectively (Buffa and Ciarlet, 2001a,b). By duality, we define the tangential operators $\nabla_{\Gamma}$ and $\text{curl}_{\Gamma}$ as the adjoint operators of $-\nabla_{\Gamma}$ and $\text{curl}_{\Gamma}$, respectively:

$$ \nabla_{\Gamma} : H_{\parallel}^{-1/2}(\Gamma) \to H^{-3/2}(\Gamma) \text{ defined by } \langle \nabla_{\Gamma} \cdot \mathbf{v}, \varphi \rangle_{-3/2,3/2} := - \langle \mathbf{v}, \nabla_{\Gamma} \varphi \rangle_{\parallel}, \quad (2.13) $$

$$ \text{curl}_{\Gamma} : H_{\perp}^{-1/2}(\Gamma) \to H^{-3/2}(\Gamma) \text{ defined by } \langle \text{curl}_{\Gamma} \mathbf{v}, \varphi \rangle_{-3/2,3/2} := \langle \mathbf{v}, \text{curl}_{\Gamma} \varphi \rangle_{\perp,1/2,-1/2}. \quad (2.14) $$

In the context of sufficiently smooth tangential vector fields, the 2D representation face by face of the operators $\nabla_{\Gamma}$ and $\text{curl}_{\Gamma}$ can be also used:

$$ \nabla_{\Gamma_i} \cdot \mathbf{v} = \nabla \cdot \mathbf{v}^{2D}, \quad \text{curl}_{\Gamma_i} \mathbf{v} = \text{curl} \mathbf{v}^{2D}, \quad (2.15) $$

where $\mathbf{v}^{2D}(x_1, x_2) := (\mathbf{v}(O_i + x_1 \boldsymbol{\tau}_1^i + x_2 \boldsymbol{\tau}_2^i) \cdot \boldsymbol{\tau}_1^i, \mathbf{v}(O_i + x_1 \boldsymbol{\tau}_1^i + x_2 \boldsymbol{\tau}_2^i) \cdot \boldsymbol{\tau}_2^i)^T$. With the tangential operators $\nabla_{\Gamma}$ and $\text{curl}_{\Gamma}$, we recall the following tangential spaces, especially $H^{-1/2}(\text{div}_{\Gamma}, \Gamma)$ that will be used intensively throughout the article:

$$ H^{-1/2}(\text{div}_{\Gamma}, \Gamma) := \{ \mathbf{v} \in H_{\parallel}^{-1/2}(\Gamma) \mid \nabla_{\Gamma} \cdot \mathbf{v} \in H^{-1/2}(\Gamma) \}, \quad (2.16) $$

$$ H^{-1/2}(\text{curl}_{\Gamma}, \Gamma) := \{ \mathbf{v} \in H_{\perp}^{-1/2}(\Gamma) \mid \text{curl}_{\Gamma} \mathbf{v} \in H^{-1/2}(\Gamma) \}. \quad (2.17) $$

Their associated norms are defined as

$$ \|\mathbf{v}\|_{H^{-1/2}(\text{div}_{\Gamma}, \Gamma)}^2 := \|\mathbf{v}\|_{H_{\parallel}^{-1/2}(\Gamma)}^2 + \|\nabla_{\Gamma} \cdot \mathbf{v}\|_{H^{-1/2}(\Gamma)}^2, \quad (2.18) $$

$$ \|\mathbf{v}\|_{H^{-1/2}(\text{curl}_{\Gamma}, \Gamma)}^2 := \|\mathbf{v}\|_{H_{\perp}^{-1/2}(\Gamma)}^2 + \|\text{curl}_{\Gamma} \mathbf{v}\|_{H^{-1/2}(\Gamma)}^2, \quad (2.19) $$

where the $H_{\parallel}^{-1/2}(\Gamma)$- and $H_{\perp}^{1/2}(\Gamma)$-norms are:

$$ \|\mathbf{v}\|_{H_{\parallel}^{-1/2}(\Gamma)} := \sup_{\mathbf{w} \in H_{\parallel}^{1/2}(\Gamma)} \frac{|\langle \mathbf{w}, \mathbf{v} \rangle_{\parallel}|}{\|\mathbf{w}\|_{H_{\parallel}^{1/2}(\Gamma)}}, \quad (2.20) $$

$$ \|\mathbf{v}\|_{H_{\parallel}^{1/2}(\Gamma)} := \inf_{\mathbf{w} \in H_{\parallel}^{1/2}(\Gamma)} \{ \|\mathbf{w}\|_{H^{1/2}(\Gamma)} \mid \pi_t \mathbf{w} = \mathbf{v} \}. \quad (2.21) $$

Now, we can give two fundamental results about $\gamma_t$ and $\pi_t$ (see Buffa and Ciarlet (2001a,b); Buffa et al. (2002b)): these operators can be extended as continuous mapping from $\mathbf{H}$ (curl, $\Omega$) to $\mathbf{H}^{-1/2}(\text{div}_{\Gamma}, \Gamma)$ and $\mathbf{H}^{-1/2}(\text{curl}_{\Gamma}, \Gamma)$, respectively.

Finally, we introduce spaces with local extra-regularity in order to reveal the h-convergence rate of corollary 4.6. The space $\mathbf{H}^{\sigma}(\text{div}_{\Gamma}, \Gamma)$ is defined as follows

$$ \mathbf{H}^{\sigma}(\text{div}_{\Gamma}, \Gamma) = \{ \mathbf{v} \in \mathbf{H}_{-}^{\sigma}(\Gamma) : \text{div}_{\Gamma}(\mathbf{v}) \in \mathbf{H}^{\sigma}(\Gamma) \}, \quad (2.22) $$

with $\mathbf{H}_{-}^{\sigma}(\Gamma) := \{ \mathbf{v} \in L_t^2(\Gamma) : \forall i = 1, \dots, N_{\Gamma}, \mathbf{v}|_{\Gamma_i} \in \mathbf{H}_t^{\sigma}(\Gamma_i) \}$ and

$$ \|\mathbf{v}\|_{\mathbf{H}^{\sigma}(\text{div}_{\Gamma}, \Gamma)}^2 := \|\mathbf{v}\|_{\mathbf{H}_{-}^{\sigma}(\Gamma)}^2 + \|\nabla_{\Gamma} \cdot \mathbf{v}\|_{\mathbf{H}^{\sigma}(\Gamma)}^2, \quad (2.23) $$

with $\|\mathbf{v}\|_{\mathbf{H}^{\sigma}(\Gamma)}^2 = \sum_{j=1}^{N_{\Gamma}} \|\mathbf{v}\|_{\mathbf{H}_{\sigma}^{\sigma}(\Gamma_j)}^2$ and $\mathbf{H}_t^{\sigma}(\Gamma_j)$ being the space of tangential vector fields of regularity $\sigma \ge 0$ on the plane surface $\Gamma_j$.
These two tangential operators can be extended as a linear continuous mapping from $H^{3/2}(\Gamma)$ to $H_{\parallel}^{1/2}(\Gamma)$ and $H_{\perp}^{1/2}(\Gamma)$, respectively (Buffa and Ciarlet, 2001a,b). By duality, we define the tangential operators $\nabla_{\Gamma}$ and $\text{curl}_{\Gamma}$ as the adjoint operators of $-\nabla_{\Gamma}$ and $\text{curl}_{\Gamma}$, respectively:

$$ \nabla_{\Gamma} : H_{\parallel}^{-1/2}(\Gamma) \to H^{-3/2}(\Gamma) \text{ defined by } \langle \nabla_{\Gamma} \cdot \mathbf{v}, \varphi \rangle_{-3/2,3/2} := - \langle \mathbf{v}, \nabla_{\Gamma} \varphi \rangle_{\parallel}, \quad (2.13) $$

$$ \text{curl}_{\Gamma} : H_{\perp}^{-1/2}(\Gamma) \to H^{-3/2}(\Gamma) \text{ defined by } \langle \text{curl}_{\Gamma} \mathbf{v}, \varphi \rangle_{-3/2,3/2} := \langle \mathbf{v}, \text{curl}_{\Gamma} \varphi \rangle_{\perp,1/2,-1/2}. \quad (2.14) $$

In the context of sufficiently smooth tangential vector fields, the 2D representation face by face of the operators $\nabla_{\Gamma}$ and $\text{curl}_{\Gamma}$ can be also used:

$$ \nabla_{\Gamma_i} \cdot \mathbf{v} = \nabla \cdot \mathbf{v}^{2D}, \quad \text{curl}_{\Gamma_i} \mathbf{v} = \text{curl} \mathbf{v}^{2D}, \quad (2.15) $$

where $\mathbf{v}^{2D}(x_1, x_2) := (\mathbf{v}(O_i + x_1 \boldsymbol{\tau}_1^i + x_2 \boldsymbol{\tau}_2^i) \cdot \boldsymbol{\tau}_1^i, \mathbf{v}(O_i + x_1 \boldsymbol{\tau}_1^i + x_2 \boldsymbol{\tau}_2^i) \cdot \boldsymbol{\tau}_2^i)^T$. With the tangential operators $\nabla_{\Gamma}$ and $\text{curl}_{\Gamma}$, we recall the following tangential spaces, especially $H^{-1/2}(\text{div}_{\Gamma}, \Gamma)$ that will be used intensively throughout the article:

$$ H^{-1/2}(\text{div}_{\Gamma}, \Gamma) := \{ \mathbf{v} \in H_{\parallel}^{-1/2}(\Gamma) \mid \nabla_{\Gamma} \cdot \mathbf{v} \in H^{-1/2}(\Gamma) \}, \quad (2.16) $$

$$ H^{-1/2}(\text{curl}_{\Gamma}, \Gamma) := \{ \mathbf{v} \in H_{\perp}^{-1/2}(\Gamma) \mid \text{curl}_{\Gamma} \mathbf{v} \in H^{-1/2}(\Gamma) \}. \quad (2.17) $$

Their associated norms are defined as

$$ \|\mathbf{v}\|_{H^{-1/2}(\text{div}_{\Gamma}, \Gamma)}^2 := \|\mathbf{v}\|_{H_{\parallel}^{-1/2}(\Gamma)}^2 + \|\nabla_{\Gamma} \cdot \mathbf{v}\|_{H^{-1/2}(\Gamma)}^2, \quad (2.18) $$

$$ \|\mathbf{v}\|_{H^{-1/2}(\text{curl}_{\Gamma}, \Gamma)}^2 := \|\mathbf{v}\|_{H_{\perp}^{-1/2}(\Gamma)}^2 + \|\text{curl}_{\Gamma} \mathbf{v}\|_{H^{-1/2}(\Gamma)}^2, \quad (2.19) $$

where the $H_{\parallel}^{-1/2}(\Gamma)$- and $H_{\perp}^{1/2}(\Gamma)$-norms are:

$$ \|\mathbf{v}\|_{H_{\parallel}^{-1/2}(\Gamma)} := \sup_{\mathbf{w} \in H_{\parallel}^{1/2}(\Gamma)} \frac{|\langle \mathbf{w}, \mathbf{v} \rangle_{\parallel}|}{\|\mathbf{w}\|_{H_{\parallel}^{1/2}(\Gamma)}}, \quad (2.20) $$

$$ \|\mathbf{v}\|_{H_{\parallel}^{1/2}(\Gamma)} := \inf_{\mathbf{w} \in H_{\parallel}^{1/2}(\Gamma)} \{ \|\mathbf{w}\|_{H^{1/2}(\Gamma)} \mid \pi_t \mathbf{w} = \mathbf{v} \}. \quad (2.21) $$

Now, we can give two fundamental results about $\gamma_t$ and $\pi_t$ (see Buffa and Ciarlet (2001a,b); Buffa et al. (2002b)): these operators can be extended as continuous mapping from $\mathbf{H}$ (curl, $\Omega$) to $\mathbf{H}^{-1/2}(\text{div}_{\Gamma}, \Gamma)$ and $\mathbf{H}^{-1/2}(\text{curl}_{\Gamma}, \Gamma)$, respectively.

Finally, we introduce spaces with local extra-regularity in order to reveal the h-convergence rate of corollary 4.6. The space $\mathbf{H}^{\sigma}(\text{div}_{\Gamma}, \Gamma)$ is defined as follows

$$ \mathbf{H}^{\sigma}(\text{div}_{\Gamma}, \Gamma) = \{ \mathbf{v} \in \mathbf{H}_{-}^{\sigma}(\Gamma) : \text{div}_{\Gamma}(\mathbf{v}) \in \mathbf{H}^{\sigma}(\Gamma) \}, \quad (2.22) $$

with $\mathbf{H}_{-}^{\sigma}(\Gamma) := \{ \mathbf{v} \in L_t^2(\Gamma) : \forall i = 1, \dots, N_{\Gamma}, \mathbf{v}|_{\Gamma_i} \in \mathbf{H}_t^{\sigma}(\Gamma_i) \}$ and

$$ \|\mathbf{v}\|_{\mathbf{H}^{\sigma}(\text{div}_{\Gamma}, \Gamma)}^2 := \|\mathbf{v}\|_{\mathbf{H}_{-}^{\sigma}(\Gamma)}^2 + \|\nabla_{\Gamma} \cdot \mathbf{v}\|_{\mathbf{H}^{\sigma}(\Gamma)}^2, \quad (2.23) $$

with $\|\mathbf{v}\|_{\mathbf{H}^{\sigma}(\Gamma)}^2 = \sum_{j=1}^{N_{\Gamma}} \|\mathbf{v}\|_{\mathbf{H}_{\sigma}^{\sigma}(\Gamma_j)}^2$ and $\mathbf{H}_t^{\sigma}(\Gamma_j)$ being the space of tangential vector fields of regularity $\sigma \ge 0$ on the plane surface $\Gamma_j$.
## 2.2. Mesh regularity properties

Henceforth, let $(\mathcal{T}_h)_{h>0}$ be a family of meshes constituted of $N_{\mathcal{T}_h}$ polygonal elements that provides an approximate surface corresponding exactly to $\Gamma$. For each element $K$ of a given $\mathcal{T}_h$, its set of edges is $\mathcal{E}_K$, its diameter is $h_K$, being the maximum distance between two vertices of $K$, and at each edge $e \in \mathcal{E}_K$ of length $h_e$ is associated the unit normal vector $n_e^K$ pointing outward the element $K$ and living on its plane. The center of mass of an element $K$ is denoted $x_K$. We also denote by $n_K$ the restriction of $n$ to $K$ and by $h_K^{\min}$ the minimum distance between two vertices of $K$. In the rest of the manuscript, we adopt from the virtual element literature the following standard regularity assumptions for polygonal meshes (for more details see, e.g., (Beirão da Veiga et al., 2022) and section 2 of Beirão da Veiga and Mascotto (2022)): $\forall h > 0, \exists \gamma \in (0, 1)$ independent of $h$ such that

* every polygonal element $K \in \mathcal{T}_h$ is star-shaped with respect to a ball of diameter $\gamma h_K$,

* for every element $K \in \mathcal{T}_h$, for each edge $e \in \mathcal{E}_K$, $\gamma h_K \le h_e$.

# 3. Continuous and discrete problems

In this section, the well-known EFIE used to model the scattering of time-harmonic waves by a perfectly conducting body in homogeneous medium and its classical RT finite element approximation are first recalled. Next, we devise a new discretization of the EFIE by using a divergence-conforming virtual element space of lowest order, which we simply term V-BEM.

## 3.1. From the continuous model to its classical discrete formulation

We consider the first-order Maxwell equations in time-harmonic domain (with time dependence of $e^{-i\omega t}$):

$$ \left\{ \begin{array}{ll} \text{curl}\mathbf{E} - i\kappa Z_0 \mathbf{H} = 0, & \text{in } \mathbb{R}^3 \setminus \bar{\Omega}, \\ \text{curl}\mathbf{H} + i\kappa Z_0^{-1} \mathbf{E} = 0, & \text{in } \mathbb{R}^3 \setminus \bar{\Omega}, \end{array} \right. \qquad (3.1) $$

where $\mathbf{E}$ and $\mathbf{H}$ are complex-valued vectorial functions corresponding to the total electric and magnetic fields, respectively, $Z_0 = \sqrt{\mu_0/\varepsilon_0}$, with $\varepsilon_0$ and $\mu_0$ being the electric permittivity and magnetic permeability in the vacuum respectively, and $\kappa = \omega/c_0 > 0$ is the wavenumber, with $c_0 = 1/\sqrt{\varepsilon_0\mu_0}$ being the speed of light. The perfectly conducting nature of the scattering body occupying $\Omega$ is taken into account by imposing a Perfect Electric Conductor (PEC) boundary condition:

$$ \mathbf{E} \times \mathbf{n} = 0 \quad \text{on } \Gamma. \qquad (3.2) $$

Finally, the problem is closed by using the so-called Silver-Müller radiation condition at infinity which consists in selecting only outgoing waves:

$$ \lim_{|\mathbf{x}| \to +\infty} |\mathbf{x}| \left( Z_0(\mathbf{H} - \mathbf{H}^I)(\mathbf{x}) \times \frac{\mathbf{x}}{|\mathbf{x}|} - (\mathbf{E} - \mathbf{E}^I)(\mathbf{x}) \right) = 0, \qquad (3.3) $$

where $(\mathbf{E}^I, \mathbf{H}^I)$ is an incident plane wave solution of (3.1) in the entire domain $\mathbb{R}^3$.

The solution $(\mathbf{E}, \mathbf{H})$ of this problem can be parametrized by the electric $J = \gamma_t \mathbf{H}$ and magnetic $M = -\gamma_t \mathbf{E}$ currents corresponding to complex-valued tangential vector fields living on 2-manifold
## 2.2. Mesh regularity properties

Henceforth, let $(\mathcal{T}_h)_{h>0}$ be a family of meshes constituted of $N_{\mathcal{T}_h}$ polygonal elements that provides an approximate surface corresponding exactly to $\Gamma$. For each element $K$ of a given $\mathcal{T}_h$, its set of edges is $\mathcal{E}_K$, its diameter is $h_K$, being the maximum distance between two vertices of $K$, and at each edge $e \in \mathcal{E}_K$ of length $h_e$ is associated the unit normal vector $n_e^K$ pointing outward the element $K$ and living on its plane. The center of mass of an element $K$ is denoted $\mathbf{x}_K$. We also denote by $n_K$ the restriction of $\mathbf{n}$ to $K$ and by $h_K^{\min}$ the minimum distance between two vertices of $K$. In the rest of the manuscript, we adopt from the virtual element literature the following standard regularity assumptions for polygonal meshes (for more details see, e.g., (Beirão da Veiga et al., 2022) and section 2 of Beirão da Veiga and Mascotto (2022)): $\forall h > 0, \exists \gamma \in (0, 1)$ independent of $h$ such that

* every polygonal element $K \in \mathcal{T}_h$ is star-shaped with respect to a ball of diameter $\gamma h_K$,

* for every element $K \in \mathcal{T}_h$, for each edge $e \in \mathcal{E}_K$, $\gamma h_K \le h_e$.

# 3. Continuous and discrete problems

In this section, the well-known EFIE used to model the scattering of time-harmonic waves by a perfectly conducting body in homogeneous medium and its classical RT finite element approximation are first recalled. Next, we devise a new discretization of the EFIE by using a divergence-conforming virtual element space of lowest order, which we simply term V-BEM.

## 3.1. From the continuous model to its classical discrete formulation

We consider the first-order Maxwell equations in time-harmonic domain (with time dependence of $e^{-i\omega t}$):

$$ \left\{ \begin{array}{ll} \text{curl}\mathbf{E} - i\kappa Z_0 \mathbf{H} = 0, & \text{in } \mathbb{R}^3 \setminus \bar{\Omega}, \\ \text{curl}\mathbf{H} + i\kappa Z_0^{-1} \mathbf{E} = 0, & \text{in } \mathbb{R}^3 \setminus \bar{\Omega}, \end{array} \right. \qquad (3.1) $$

where $\mathbf{E}$ and $\mathbf{H}$ are complex-valued vectorial functions corresponding to the total electric and magnetic fields, respectively, $Z_0 = \sqrt{\mu_0/\varepsilon_0}$, with $\varepsilon_0$ and $\mu_0$ being the electric permittivity and magnetic permeability in the vacuum respectively, and $\kappa = \omega/c_0 > 0$ is the wavenumber, with $c_0 = 1/\sqrt{\varepsilon_0\mu_0}$ being the speed of light. The perfectly conducting nature of the scattering body occupying $\Omega$ is taken into account by imposing a Perfect Electric Conductor (PEC) boundary condition:

$$ \mathbf{E} \times \mathbf{n} = 0 \quad \text{on } \Gamma. \qquad (3.2) $$

Finally, the problem is closed by using the so-called Silver-Müller radiation condition at infinity which consists in selecting only outgoing waves:

$$ \lim_{|\mathbf{x}| \to +\infty} |\mathbf{x}| \left( Z_0(\mathbf{H} - \mathbf{H}^I)(\mathbf{x}) \times \frac{\mathbf{x}}{|\mathbf{x}|} - (\mathbf{E} - \mathbf{E}^I)(\mathbf{x}) \right) = 0, \qquad (3.3) $$

where $(\mathbf{E}^I, \mathbf{H}^I)$ is an incident plane wave solution of (3.1) in the entire domain $\mathbb{R}^3$.

The solution $(\mathbf{E}, \mathbf{H})$ of this problem can be parametrized by the electric $\mathbf{J} = \gamma_t \mathbf{H}$ and magnetic $\mathbf{M} = -\gamma_t \mathbf{E}$ currents corresponding to complex-valued tangential vector fields living on 2-manifold
Γ. For that, the Stratton-Chu representation formulae (Monk, 2003; Stratton and Chu, 1939) are used:

$$ \left\{ \begin{array}{ll} \mathbf{E}(\mathbf{x}) = \mathbf{E}^I(\mathbf{x}) + \iota\kappa Z_0(\mathcal{T}\mathbf{J})(\mathbf{x}) + (\mathcal{K}\mathbf{M})(\mathbf{x}), & \mathbf{x} \in \mathbb{R}^3 \setminus \bar{\Omega}, \\ \mathbf{H}(\mathbf{x}) = \mathbf{H}^I(\mathbf{x}) + \iota\kappa Z_0^{-1}(\mathcal{T}\mathbf{M})(\mathbf{x}) - (\mathcal{K}\mathbf{J})(\mathbf{x}), & \mathbf{x} \in \mathbb{R}^3 \setminus \bar{\Omega}, \end{array} \right. \qquad (3.4) $$

where the potential operators are defined as

$$ \mathcal{T}: \left| \begin{array}{ll} \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma) & \rightarrow \mathbf{H}_{\mathrm{loc}}(\mathrm{curl}^2, \mathbb{R}^3 \setminus \Gamma) \cap \mathbf{H}_{\mathrm{loc}}(\mathrm{div}0, \mathbb{R}^3 \setminus \Gamma) \\ \mathbf{J}(\mathbf{x}) & \rightarrow \left(\frac{1}{\kappa^2}\nabla\nabla\cdot+1\right) \int_{\Gamma} G_{\kappa}(\mathbf{x}-\mathbf{y})\mathbf{J}(\mathbf{y})d\gamma_y \end{array} \right. , \qquad (3.5) $$

$$ \mathcal{K}: \left| \begin{array}{ll} \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma) & \rightarrow \mathbf{H}_{\mathrm{loc}}(\mathrm{curl}^2, \mathbb{R}^3 \setminus \Gamma) \cap \mathbf{H}_{\mathrm{loc}}(\mathrm{div}0, \mathbb{R}^3 \setminus \Gamma) \\ \mathbf{J}(\mathbf{x}) & \rightarrow -\nabla \times \int_{\Gamma} G_{\kappa}(\mathbf{x}-\mathbf{y})\mathbf{J}(\mathbf{y})d\gamma_y \end{array} \right. , \qquad (3.6) $$

with $G_{\kappa}$ being the Green function associated to the radiating solutions of the 3-D Helmholtz equation:

$$ G_{\kappa}(\mathbf{x}) = \frac{e^{\iota\kappa|\mathbf{x}|}}{4\pi|\mathbf{x}|}, \quad \mathbf{x} \neq 0. \qquad (3.7) $$

In this paper, we only work with the EFIE that is obtained by computing the tangential component of the first equation of (3.4) and by using the PEC boundary condition $M = 0$:

$$ \pi_t \mathbf{E}(\mathbf{x}) = \pi_t \mathbf{E}^I(\mathbf{x}) + \iota\kappa Z_0 \pi_t(\mathcal{T}\mathbf{J})(\mathbf{x}) \Rightarrow \pi_t \mathbf{E}^I(\mathbf{x}) + \iota\kappa Z_0 \pi_t(\mathcal{T}\mathbf{J})(\mathbf{x}) = 0, \text{ for } \mathbf{x} \in \Gamma. \qquad (3.8) $$

The tangential component of $\mathcal{T}$ is a boundary integral operator which induces the following weak formulation of (3.8)

$$ \begin{array}{ll} \text{Find } \mathbf{J} \in \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma) \text{ such that, } & \forall \mathbf{J}' \in \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma), \\ a(\mathbf{J}, \mathbf{J}') = f(\mathbf{J}'). & \end{array} \qquad (3.9) $$

Here, $a$ is the sesquilinear form defined as follows

$$ a: \left| \begin{array}{ll} \left( \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma) \right)^2 & \rightarrow \mathbb{C} \\ (\mathbf{J}, \mathbf{J}') & \rightarrow \langle \mathbf{V}_k \mathbf{J}, \mathbf{J}' \rangle_{\parallel} - \frac{1}{\kappa^2} \langle \mathbf{V}_k \nabla_{\Gamma} \cdot \mathbf{J}, \nabla_{\Gamma} \cdot \mathbf{J}' \rangle, \end{array} \right. \qquad (3.10) $$

where, according to Buffa and Hiptmair (2003):

$$ \begin{array}{ll} \mathbf{V}_k: \mathbf{H}_{\parallel}^{-1/2}(\Gamma) \rightarrow \mathbf{H}_{\parallel}^{1/2}(\Gamma), & \\ \langle \mathbf{V}_k \mathbf{J}, \mathbf{J}' \rangle_{\parallel} = \int_{\Gamma} \int_{\Gamma} G_{\kappa}(\mathbf{x}-\mathbf{y}) \mathbf{J}(\mathbf{y}) \cdot \mathbf{J}'(\mathbf{x}) d\gamma_y d\gamma_x, & (3.11) \end{array} $$

$$ \begin{array}{ll} \mathbf{V}_k: H^{-1/2}(\Gamma) \rightarrow H^{1/2}(\Gamma), & \\ \langle \mathbf{V}_k \nabla_{\Gamma} \cdot \mathbf{J}, \nabla_{\Gamma} \cdot \mathbf{J}' \rangle = \int_{\Gamma} \int_{\Gamma} G_{\kappa}(\mathbf{x}-\mathbf{y}) \nabla_{\Gamma} \cdot \mathbf{J}(\mathbf{y}) \nabla_{\Gamma} \cdot \mathbf{J}'(\mathbf{x}) d\gamma_y d\gamma_x. & (3.12) \end{array} $$
Γ. For that, the Stratton-Chu representation formulae (Monk, 2003; Stratton and Chu, 1939) are used:

$$ \left\{ \begin{array}{ll} \mathbf{E}(\mathbf{x}) = \mathbf{E}^I(\mathbf{x}) + \iota\kappa Z_0(\mathcal{T}\mathbf{J})(\mathbf{x}) + (\mathcal{K}\mathbf{M})(\mathbf{x}), & \mathbf{x} \in \mathbb{R}^3 \setminus \bar{\Omega}, \\ \mathbf{H}(\mathbf{x}) = \mathbf{H}^I(\mathbf{x}) + \iota\kappa Z_0^{-1}(\mathcal{T}\mathbf{M})(\mathbf{x}) - (\mathcal{K}\mathbf{J})(\mathbf{x}), & \mathbf{x} \in \mathbb{R}^3 \setminus \bar{\Omega}, \end{array} \right. \qquad (3.4) $$

where the potential operators are defined as

$$ \mathcal{T}: \left| \begin{array}{ll} \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma) & \rightarrow \mathbf{H}_{\mathrm{loc}}(\mathrm{curl}^2, \mathbb{R}^3 \setminus \Gamma) \cap \mathbf{H}_{\mathrm{loc}}(\mathrm{div}0, \mathbb{R}^3 \setminus \Gamma) \\ \mathbf{J}(\mathbf{x}) & \rightarrow \left(\frac{1}{\kappa^2}\nabla\nabla\cdot+1\right) \int_{\Gamma} G_{\kappa}(\mathbf{x}-\mathbf{y})\mathbf{J}(\mathbf{y})d\gamma_y \end{array} \right. , \qquad (3.5) $$

$$ \mathcal{K}: \left| \begin{array}{ll} \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma) & \rightarrow \mathbf{H}_{\mathrm{loc}}(\mathrm{curl}^2, \mathbb{R}^3 \setminus \Gamma) \cap \mathbf{H}_{\mathrm{loc}}(\mathrm{div}0, \mathbb{R}^3 \setminus \Gamma) \\ \mathbf{J}(\mathbf{x}) & \rightarrow -\nabla \times \int_{\Gamma} G_{\kappa}(\mathbf{x}-\mathbf{y})\mathbf{J}(\mathbf{y})d\gamma_y \end{array} \right. , \qquad (3.6) $$

with $G_κ$ being the Green function associated to the radiating solutions of the 3-D Helmholtz equation:

$$ G_{\kappa}(\mathbf{x}) = \frac{e^{\iota\kappa|\mathbf{x}|}}{4\pi|\mathbf{x}|}, \quad \mathbf{x} \neq 0. \qquad (3.7) $$

In this paper, we only work with the EFIE that is obtained by computing the tangential component of the first equation of (3.4) and by using the PEC boundary condition $M = 0$:

$$ \pi_t \mathbf{E}(\mathbf{x}) = \pi_t \mathbf{E}^I(\mathbf{x}) + \iota\kappa Z_0 \pi_t(\mathcal{T}\mathbf{J})(\mathbf{x}) \Rightarrow \pi_t \mathbf{E}^I(\mathbf{x}) + \iota\kappa Z_0 \pi_t(\mathcal{T}\mathbf{J})(\mathbf{x}) = 0, \text{ for } \mathbf{x} \in \Gamma. \qquad (3.8) $$

The tangential component of $\mathcal{T}$ is a boundary integral operator which induces the following weak formulation of (3.8)

$$ \begin{array}{l} \text{Find } \mathbf{J} \in \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma) \text{ such that, } \forall \mathbf{J}' \in \mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma), \\ a(\mathbf{J}, \mathbf{J}') = f(\mathbf{J}'). \end{array} \qquad (3.9) $$

Here, $a$ is the sesquilinear form defined as follows

$$ a: \left| \begin{array}{ll} \left(\mathbf{H}^{-1/2}(\mathrm{div}_{\Gamma}, \Gamma)\right)^2 & \rightarrow \mathbb{C} \\ (\mathbf{J}, \mathbf{J}') & \rightarrow \langle \mathbf{V}_k \mathbf{J}, \mathbf{J}' \rangle_{\parallel} - \frac{1}{\kappa^2} \langle \mathbf{V}_k \nabla_{\Gamma} \cdot \mathbf{J}, \nabla_{\Gamma} \cdot \mathbf{J}' \rangle, \end{array} \right. \qquad (3.10) $$

where, according to Buffa and Hiptmair (2003):

$$ \begin{array}{l} \mathbf{V}_k: \mathbf{H}_{\parallel}^{-1/2}(\Gamma) \rightarrow \mathbf{H}_{\parallel}^{1/2}(\Gamma), \\ \langle \mathbf{V}_k \mathbf{J}, \mathbf{J}' \rangle_{\parallel} = \int_{\Gamma} \int_{\Gamma} G_{\kappa}(\mathbf{x}-\mathbf{y}) \mathbf{J}(\mathbf{y}) \cdot \mathbf{J}'(\mathbf{x}) d\gamma_y d\gamma_x, \end{array} \qquad (3.11) $$

$$ \begin{array}{l} \mathbf{V}_k: H^{-1/2}(\Gamma) \rightarrow H^{1/2}(\Gamma), \\ \langle \mathbf{V}_k \nabla_{\Gamma} \cdot \mathbf{J}, \nabla_{\Gamma} \cdot \mathbf{J}' \rangle = \int_{\Gamma} \int_{\Gamma} G_{\kappa}(\mathbf{x}-\mathbf{y}) \nabla_{\Gamma} \cdot \mathbf{J}(\mathbf{y}) \nabla_{\Gamma} \cdot \mathbf{J}'(\mathbf{x}) d\gamma_y d\gamma_x. \end{array} \qquad (3.12) $$
The right-hand side of the EFIE (3.9) has the following form

$$f(\mathbf{J}') = -\frac{1}{i\kappa Z_0} \langle \pi_t \mathbf{E}^I, \mathbf{J}' \rangle_{\parallel}. \qquad (3.13)$$

It is well known that this continuous problem is well-posed except for a discrete number of wavenumbers corresponding to resonance frequencies of the enlightened body (Buffa and Hiptmair, 2003). In this paper, we assume to be far from those frequencies in order to discretize the problem. The classical approach of Galerkin discretization for this equation is to use a boundary element method based on the lowest-order RT finite element space (Raviart and Thomas, 1977) defined on a simplicial mesh $T_h$ (i.e. composed of triangles):

$$\mathcal{R}\mathcal{T}_h^0 = \{ \mathbf{v}_h \in \mathbf{H} (\mathrm{div}_{\Gamma}, \Gamma) \mid \forall K \in \mathcal{T}_h, \mathbf{v}_{h|_K}^{2\mathrm{D}} \in \mathcal{R}\mathcal{T}_0(K) \}, \qquad (3.14)$$

$$\mathcal{R}\mathcal{T}_0(K) = \{ (\mathbb{P}_0(K))^2 \oplus \mathbf{x}_G^{2\mathrm{D}}\mathbb{P}_0(K) \}, \qquad (3.15)$$

where $\forall \mathbf{x} \in K$, $\mathbf{x}_G^{2\mathrm{D}} = \mathbf{x}^{2\mathrm{D}} - \mathbf{x}_K^{2\mathrm{D}}$ and $\mathbb{P}_s(X)$ denotes the space of polynomials of total degree $\le s$ defined over $X$. By looking for solutions on this RT finite dimensional space, the discrete weak formulation of (3.9) can be then written as

$$\begin{gathered} \text{Find } \mathbf{J}_h \in \mathcal{R}\mathcal{T}_h^0 \text{ such that, } \forall \mathbf{J}'_h \in \mathcal{R}\mathcal{T}_h^0, \\ a(\mathbf{J}_h, \mathbf{J}'_h) = f(\mathbf{J}'_h). \end{gathered} \qquad (3.16)$$

As an approximation of a Fredholm-type equation, the problem (3.16) is well-posed for sufficiently fine meshes (Buffa and Hiptmair, 2003) and it will be used to numerically evaluate the performance of the new numerical scheme based on the natural virtual element counterpart of the above RT elements devised in the next subsection.

## 3.2. A Virtual Boundary Element Method (V-BEM) for the EFIE

From now on, the mesh $T_h$ is assumed to be composed of polygonal elements and satisfies the regularity properties stated in section 2.2. We emphasize that, here, the term polygon also refers to any triangular- (or quadrangular-) shaped element with more than three (four, respectively) vertices, some of which may be collinear. The standard construction of a VEM is based on two main ingredients:

* A local finite dimensional space on each polygon $K$ whose virtual functions are implicitly defined as solutions to PDE respecting the nature of the global continuity of the physical solution.

* An operator of projection from the local virtual space onto a suited polynomial space that can be explicitly computed from the only knowledge of the degrees of freedom associated to the virtual functions. A stabilization term is generally included as a complement to the projection in order to correct the rank deficiency introduced in the discrete setting.

In the context of the EFIE, we propose a Galerkin discretization that relies on the divergence-conforming serendipity virtual space (Beirão da Veiga et al., 2018; Beirão da Veiga and Mascotto, 2022, (3.14)): for each polygon $K$ the local virtual space reads

$$\begin{gathered} V_0^e(K) = \{ \mathbf{v}^{2\mathrm{D}} : K \to \mathbb{C}^2 \mid \forall e \in \mathcal{E}_K, \mathbf{v}_e^{2\mathrm{D}} \cdot \mathbf{n}_e^K \in \mathbb{P}_0(e), \nabla \cdot \mathbf{v}^{2\mathrm{D}} \in \mathbb{P}_0(K), \\ \mathrm{curl} \mathbf{v}^{2\mathrm{D}} \in \mathbb{P}_0(K), \int_K \mathbf{v}^{2\mathrm{D}} \cdot \mathbf{x}_G^{2\mathrm{D},\perp} d\gamma = 0 \}. \end{gathered} \qquad (3.17)$$
The right-hand side of the EFIE (3.9) has the following form

$$f(\mathbf{J}') = -\frac{1}{i\kappa Z_0} \langle \pi_t \mathbf{E}^I, \mathbf{J}' \rangle_{\parallel} \quad (3.13)$$

It is well known that this continuous problem is well-posed except for a discrete number of wavenumbers corresponding to resonance frequencies of the enlightened body (Buffa and Hiptmair, 2003). In this paper, we assume to be far from those frequencies in order to discretize the problem. The classical approach of Galerkin discretization for this equation is to use a boundary element method based on the lowest-order RT finite element space (Raviart and Thomas, 1977) defined on a simplicial mesh $T_h$ (i.e. composed of triangles):

$$\mathcal{R}\mathcal{T}_h^0 = \{ \mathbf{v}_h \in \mathbf{H} (\mathrm{div}_{\Gamma}, \Gamma) \mid \forall K \in \mathcal{T}_h, \mathbf{v}_{h|_K}^{2\mathrm{D}} \in \mathcal{R}\mathcal{T}_0(K) \}, \quad (3.14)$$

$$\mathcal{R}\mathcal{T}_0(K) = \{(\mathbb{P}_0(K))^2 \oplus \mathbf{x}_G^{2\mathrm{D}}\mathbb{P}_0(K)\}, \quad (3.15)$$

where $\forall \mathbf{x} \in K$, $\mathbf{x}_G^{2\mathrm{D}} = \mathbf{x}^{2\mathrm{D}} - \mathbf{x}_K^{2\mathrm{D}}$ and $\mathbb{P}_s(X)$ denotes the space of polynomials of total degree $\le s$ defined over $X$. By looking for solutions on this RT finite dimensional space, the discrete weak formulation of (3.9) can be then written as

$$\begin{gathered} \text{Find } \mathbf{J}_h \in \mathcal{R}\mathcal{T}_h^0 \text{ such that, } \forall \mathbf{J}'_h \in \mathcal{R}\mathcal{T}_h^0, \\ a(\mathbf{J}_h, \mathbf{J}'_h) = f(\mathbf{J}'_h). \end{gathered} \quad (3.16)$$

As an approximation of a Fredholm-type equation, the problem (3.16) is well-posed for sufficiently fine meshes (Buffa and Hiptmair, 2003) and it will be used to numerically evaluate the performance of the new numerical scheme based on the natural virtual element counterpart of the above RT elements devised in the next subsection.

## 3.2. A Virtual Boundary Element Method (V-BEM) for the EFIE

From now on, the mesh $T_h$ is assumed to be composed of polygonal elements and satisfies the regularity properties stated in section 2.2. We emphasize that, here, the term polygon also refers to any triangular- (or quadrangular-) shaped element with more than three (four, respectively) vertices, some of which may be collinear. The standard construction of a VEM is based on two main ingredients:

* A local finite dimensional space on each polygon $K$ whose virtual functions are implicitly defined as solutions to PDE respecting the nature of the global continuity of the physical solution.

* An operator of projection from the local virtual space onto a suited polynomial space that can be explicitly computed from the only knowledge of the degrees of freedom associated to the virtual functions. A stabilization term is generally included as a complement to the projection in order to correct the rank deficiency introduced in the discrete setting.

In the context of the EFIE, we propose a Galerkin discretization that relies on the divergence-conforming serendipity virtual space (Beirão da Veiga et al., 2018; Beirão da Veiga and Mascotto, 2022, (3.14)): for each polygon $K$ the local virtual space reads

$$\begin{gathered} V_0^e(K) = \{ \mathbf{v}^{2\mathrm{D}} : K \to \mathbb{C}^2 \mid \forall e \in \mathcal{E}_K, \mathbf{v}_e^{2\mathrm{D}} \cdot \mathbf{n}_e^K \in \mathbb{P}_0(e), \nabla \cdot \mathbf{v}^{2\mathrm{D}} \in \mathbb{P}_0(K), \\ \mathrm{curl} \mathbf{v}^{2\mathrm{D}} \in \mathbb{P}_0(K), \int_K \mathbf{v}^{2\mathrm{D}} \cdot \mathbf{x}_G^{2\mathrm{D},\perp} d\gamma = 0 \}. \end{gathered} \quad (3.17)$$
**Remark 3.1** The integral constraint inside $V_0^e$ is used as a fixed degree of freedom to help compute the local projection operator (see (3.24)).

The dimension of this space is $\dim(V_0^e) = \mathrm{card}(\mathcal{E}_K)$ and each edge of $K$ is associated to a degree of freedom that corresponds to an edge flux defined as follows

$$v \mapsto \Lambda_e^K(v) := \int_e (v \cdot n_e^K) p_0 ds, \quad \forall p_0 \in \mathbb{P}_0(e). \qquad (3.18)$$

This virtual element $(\Lambda_e^K, K, V_0^e)$ is thus unisolvant (Beirão da Veiga et al., 2018) in the sense of Ciarlet (Ciarlet, 2002) and the associated local basis functions are: let $e \in \mathcal{E}_K$,

$$\varphi_e^K \in V_0^e(K), \quad \forall \tilde{e} \in \mathcal{E}_K, \quad \Lambda_{\tilde{e}}^K(\varphi_e^K) = \delta_{e\tilde{e}}. \qquad (3.19)$$

We remark that this space verifies an important (polynomial) consistency property: $RT_0(K) \subset V_0^e(K)$ (Beirão da Veiga et al., 2018, Section 3). A particular case is when $K$ is a triangle: $RT_0(K) = V_0^e(K)$; in other words, we obtain exactly the RT approximation previously presented (and the degrees of freedom of both spaces coincide).

For a general polygon, the basis functions are non-polynomial and are only (implicitly) characterized by the four constraints used in (3.17). Nevertheless, some fundamental informations can be exactly computed from the only knowledge of the degrees of freedom (3.18): let $v: K \to \mathbb{C}^3$ such that $v \cdot n_K = 0$ and $v^{2D} \in V_0^e(K)$,

* the tangential divergence of $v$:

$$\nabla_{\Gamma} \cdot v = \frac{1}{|K|} \int_K \nabla \cdot v^{2D} d\gamma = \frac{1}{|K|} \int_{\partial K} (v^{2D} \cdot n_e^K)_{|_{\partial K}} ds = \frac{1}{|K|} \sum_{e \in \mathcal{E}_K} \Lambda_e^K (v^{2D}), \qquad (3.20)$$

where the property $\nabla_{\Gamma} \cdot v \in \mathbb{P}_0(K)$ is used,

* the $L^2$-orthogonal projection on $\mathbb{P}_{1,t}(K) := \{p \in (\mathbb{P}_1(K))^3 : p \cdot n_K = 0\}$: $\Pi_K^0 v \in \mathbb{P}_{1,t}(K)$ is defined by

$$\Pi_K^0 v = (\Pi_K^{0,2D} v^{2D})_1 \tau_1 + (\Pi_K^{0,2D} v^{2D})_2 \tau_2 \quad \text{where } \Pi_K^{0,2D}: V_0^e(K) \to (\mathbb{P}_1(K))^2 \qquad (3.21)$$

$$\int_K \Pi_K^{0,2D} v^{2D} \cdot p d\gamma = \int_K v^{2D} \cdot p d\gamma, \forall p \in (\mathbb{P}_1(K))^2. \qquad (3.22)$$

By using the following well-known decomposition of the polynomial space $(\mathbb{P}_1(K))^2$:

$$(\mathbb{P}_1(K))^2 = \nabla \mathbb{P}_2(K) \oplus x_K^{2D, \perp} \mathbb{P}_0(K), \qquad (3.23)$$

the right-hand side of (3.22) can be exactly computed in a same way as in Beirão da Veiga et al. (2018, (3.15)): let $p \in (\mathbb{P}_1(K))^2$ such that $p = \nabla p_2 + p_0 x_K^{2D, \perp}$,
**Remark 3.1** The integral constraint inside $V_0^e$ is used as a fixed degree of freedom to help compute the local projection operator (see (3.24)).

The dimension of this space is $\dim(V_0^e) = \mathrm{card}(\mathcal{E}_K)$ and each edge of $K$ is associated to a degree of freedom that corresponds to an edge flux defined as follows

$$v \mapsto \Lambda_e^K(v) := \int_e (v \cdot n_e^K) p_0 ds, \quad \forall p_0 \in \mathbb{P}_0(e). \qquad (3.18)$$

This virtual element $(\Lambda_e^K, K, V_0^e)$ is thus unisolvant (Beirão da Veiga et al., 2018) in the sense of Ciarlet (Ciarlet, 2002) and the associated local basis functions are: let $e \in \mathcal{E}_K$,

$$\varphi_e^K \in V_0^e(K), \quad \forall \tilde{e} \in \mathcal{E}_K, \quad \Lambda_{\tilde{e}}^K(\varphi_e^K) = \delta_{e\tilde{e}}. \qquad (3.19)$$

We remark that this space verifies an important (polynomial) consistency property: $RT_0(K) \subset V_0^e(K)$ (Beirão da Veiga et al., 2018, Section 3). A particular case is when $K$ is a triangle: $RT_0(K) = V_0^e(K)$; in other words, we obtain exactly the RT approximation previously presented (and the degrees of freedom of both spaces coincide).

For a general polygon, the basis functions are non-polynomial and are only (implicitly) characterized by the four constraints used in (3.17). Nevertheless, some fundamental informations can be exactly computed from the only knowledge of the degrees of freedom (3.18): let $v: K \to \mathbb{C}^3$ such that $v \cdot n_K = 0$ and $v^{2D} \in V_0^e(K)$,

* the tangential divergence of $v$:

$$\nabla_{\Gamma} \cdot v = \frac{1}{|K|} \int_K \nabla \cdot v^{2D} d\gamma = \frac{1}{|K|} \int_{\partial K} (v^{2D} \cdot n_e^K)_{|_{\partial K}} ds = \frac{1}{|K|} \sum_{e \in \mathcal{E}_K} \Lambda_e^K (v^{2D}), \qquad (3.20)$$

where the property $\nabla_{\Gamma} \cdot v \in \mathbb{P}_0(K)$ is used,

* the $L^2$-orthogonal projection on $\mathbb{P}_{1,t}(K) := \{p \in (\mathbb{P}_1(K))^3 : p \cdot n_K = 0\}$: $\Pi_K^0 v \in \mathbb{P}_{1,t}(K)$ is defined by

$$\Pi_K^0 v = (\Pi_K^{0,2D} v^{2D})_1 \tau_1 + (\Pi_K^{0,2D} v^{2D})_2 \tau_2 \quad \text{where} \quad \Pi_K^{0,2D}: V_0^e(K) \to (\mathbb{P}_1(K))^2 \qquad (3.21)$$

$$\int_K \Pi_K^{0,2D} v^{2D} \cdot p d\gamma = \int_K v^{2D} \cdot p d\gamma, \forall p \in (\mathbb{P}_1(K))^2. \qquad (3.22)$$

By using the following well-known decomposition of the polynomial space $(\mathbb{P}_1(K))^2$:

$$(\mathbb{P}_1(K))^2 = \nabla \mathbb{P}_2(K) \oplus x_K^{2D, \perp} \mathbb{P}_0(K), \qquad (3.23)$$

the right-hand side of (3.22) can be exactly computed in a same way as in Beirão da Veiga et al. (2018, (3.15)): let $p \in (\mathbb{P}_1(K))^2$ such that $p = \nabla p_2 + p_0 x_K^{2D, \perp}$,
$$
\begin{align}
\int_K \mathbf{v}^{2\mathrm{D}} \cdot \mathbf{p} d\gamma &= -\int_K \nabla \cdot \mathbf{v}^{2\mathrm{D}} p_2 d\gamma + \int_{\partial K} (\mathbf{v}^{2\mathrm{D}} \cdot \mathbf{n}_e^K)_{|_{\partial K}} p_2 ds + \int_K \mathbf{v}^{2\mathrm{D}} \cdot \mathbf{x}_K^{2\mathrm{D},\perp} p_0 d\gamma, \nonumber \\
&= -\int_K \nabla \cdot \mathbf{v}^{2\mathrm{D}} p_2 d\gamma + \int_{\partial K} (\mathbf{v}^{2\mathrm{D}} \cdot \mathbf{n}_e^K)_{|_{\partial K}} p_2 ds, \nonumber \\
&= -\left(\sum_{e \in \mathcal{E}_K} \Lambda_e^K (\mathbf{v}^{2\mathrm{D}})\right) \frac{1}{|K|} \int_K p_2 d\gamma + \sum_{e \in \mathcal{E}_K} \Lambda_e^K (\mathbf{v}^{2\mathrm{D}}) \int_e p_2 ds. \tag{3.24}
\end{align}
$$

With all these local ingredients, we can define the divergence-conforming global virtual space as follows

$$
\mathcal{V}_h = \{ \mathbf{v} \in \mathbf{H} (\mathrm{div}_{\Gamma}, \Gamma) \mid \mathbf{v}_{|_K}^{2\mathrm{D}} \in \mathcal{V}_0^e (K), \forall K \in \mathcal{T}_h \}. \quad (3.25)
$$

The L²-projection Π₀<sub>h</sub> of a virtual function v<sub>h</sub> ∈ V<sub>h</sub> onto the polynomial space

$$
\mathbb{P}_{1,t}(\mathcal{T}_h) := \{ \mathbf{p} : \Gamma \to \mathbb{C}^3 : \forall K \in \mathcal{T}_h, \mathbf{p}|_{K \cap \Gamma} \in \mathbb{P}_{1,t}(K) \}, \quad (3.26)
$$

is obviously defined element-wise: ∀K ∈ T_h, (Π_h^0 v_h)|_K = Π_K^0 (v|_K).

We note that the dimension of $\mathcal{V}_h$ is equal to the number of all element edges within the mesh $\mathcal{T}_h$
and the virtual basis functions are defined in the canonical way: for all edges $e = K \cap T$ with $K$ and $T$
two adjacent elements of $\mathcal{T}_h$,

$$
\varphi_e := \begin{cases} \varphi_e^K & \text{in } K, \\ -\varphi_e^T & \text{in } T, \\ 0 & \text{otherwise.} \end{cases} \tag{3.27}
$$

**Remark 3.2** The position of the sign – in (3.27) is purely arbitrary.

From the global virtual space and the associated L²-projection, we introduce two virtual boundary element discretizations of the EFIE. The first one corresponds to the practical numerical scheme (in what follows, we refer to it as V-BEM) and uses the above L²-projection operator in order to compute the term in the sesquilinear and linear forms, which are not computable due to the non-explicit knowledge of the virtual basis functions:

$$
\text{Find } \mathbf{J}_h \in \mathcal{V}_h \text{ such that } \forall \mathbf{J}'_h \in \mathcal{V}_h,
$$
$$
a_h(\mathbf{J}_h, \mathbf{J}'_h) = f_h(\mathbf{J}'_h). \tag{3.28}
$$

with

$$
a_h(\mathbf{J}_h, \mathbf{J}'_h) := \langle \mathbf{V}_K \Pi_h^0 \mathbf{J}_h, \Pi_h^0 \mathbf{J}'_h \rangle_{\|} - \frac{1}{\kappa^2} \langle V_k \nabla_{\Gamma} \cdot \mathbf{J}_h, \nabla_{\Gamma} \cdot \mathbf{J}'_h \rangle, \quad (3.29)
$$

$$
f_h(\mathbf{J}'_h) := - \frac{1}{i\kappa Z_0} \langle \pi_t \mathbf{E}^I, \Pi_h^0 \mathbf{J}'_h \rangle_{\|}. \quad (3.30)
$$

The second one is purely theoretical (we thus refer to it as ideal V-BEM that alludes to discretization
we would derive in an ideal world in which the basis functions are explicitly known) and it will be used
$$
\begin{align}
\int_K \mathbf{v}^{2\mathrm{D}} \cdot \boldsymbol{p} d\gamma &= -\int_K \nabla \cdot \mathbf{v}^{2\mathrm{D}} p_2 d\gamma + \int_{\partial K} (\mathbf{v}^{2\mathrm{D}} \cdot \boldsymbol{n}_e^K)_{|_{\partial K}} p_2 ds + \int_K \mathbf{v}^{2\mathrm{D}} \cdot \boldsymbol{x}_K^{2\mathrm{D},\perp} p_0 d\gamma, \nonumber \\
&= -\int_K \nabla \cdot \mathbf{v}^{2\mathrm{D}} p_2 d\gamma + \int_{\partial K} (\mathbf{v}^{2\mathrm{D}} \cdot \boldsymbol{n}_e^K)_{|_{\partial K}} p_2 ds, \nonumber \\
&= -\left(\sum_{e \in \mathcal{E}_K} \Lambda_e^K (\mathbf{v}^{2\mathrm{D}})\right) \frac{1}{|K|} \int_K p_2 d\gamma + \sum_{e \in \mathcal{E}_K} \Lambda_e^K (\mathbf{v}^{2\mathrm{D}}) \int_e p_2 ds. \tag{3.24}
\end{align}
$$

With all these local ingredients, we can define the divergence-conforming global virtual space as follows

$$
\mathcal{V}_h = \{ \mathbf{v} \in \mathbf{H} (\mathrm{div}_{\Gamma}, \Gamma) \mid \mathbf{v}_{|_K}^{2\mathrm{D}} \in \mathcal{V}_0^e (K), \forall K \in \mathcal{T}_h \}. \quad (3.25)
$$

The $L^2$-projection $\Pi_h^0$ of a virtual function $v_h \in \mathcal{V}_h$ onto the polynomial space

$$
\mathbb{P}_{1,t}(\mathcal{T}_h) := \{ \mathbf{p} : \Gamma \to \mathbb{C}^3 : \forall K \in \mathcal{T}_h, \mathbf{p}|_{K \cap \Gamma} \in \mathbb{P}_{1,t}(K) \}, \quad (3.26)
$$

is obviously defined element-wise: $\forall K \in \mathcal{T}_h$, $(\Pi_h^0 v_h)|_{\mathcal{K}} = \Pi_h^0 (v|_{\mathcal{K}})$.

We note that the dimension of $\mathcal{V}_h$ is equal to the number of all element edges within the mesh $\mathcal{T}_h$
and the virtual basis functions are defined in the canonical way: for all edges $e = K \cap T$ with $K$ and $T$
two adjacent elements of $\mathcal{T}_h$,

$$
\varphi_e := \begin{cases} \varphi_e^K & \text{in } K, \\ -\varphi_e^T & \text{in } T, \\ 0 & \text{otherwise.} \end{cases} \tag{3.27}
$$

**Remark 3.2** The position of the sign – in (3.27) is purely arbitrary.

From the global virtual space and the associated L²-projection, we introduce two virtual boundary element discretizations of the EFIE. The first one corresponds to the practical numerical scheme (in what follows, we refer to it as V-BEM) and uses the above L²-projection operator in order to compute the term in the sesquilinear and linear forms, which are not computable due to the non-explicit knowledge of the virtual basis functions:

$$
\text{Find } \mathbf{J}_h \in \mathcal{V}_h \text{ such that } \forall \mathbf{J}'_h \in \mathcal{V}_h,
$$
$$
a_h(\mathbf{J}_h, \mathbf{J}'_h) = f_h(\mathbf{J}'_h). \tag{3.28}
$$

with

$$
a_h(\mathbf{J}_h, \mathbf{J}'_h) := \langle \mathbf{V}_K \Pi_h^0 \mathbf{J}_h, \Pi_h^0 \mathbf{J}'_h \rangle_{\|} - \frac{1}{\kappa^2} \langle V_k \nabla_{\Gamma} \cdot \mathbf{J}_h, \nabla_{\Gamma} \cdot \mathbf{J}'_h \rangle, \quad (3.29)
$$

$$
f_h(\mathbf{J}'_h) := - \frac{1}{i\kappa Z_0} \langle \pi_t \mathbf{E}^I, \Pi_h^0 \mathbf{J}'_h \rangle_{\|}. \quad (3.30)
$$

The second one is purely theoretical (we thus refer to it as ideal V-BEM that alludes to discretization
we would derive in an ideal world in which the basis functions are explicitly known) and it will be used
to analyze the former formulation:

$$ \text{Find } \mathbf{J}_h \in \mathcal{V}_h \text{ such that } \forall \mathbf{J}'_h \in \mathcal{V}_h, \\ a(\mathbf{J}_h, \mathbf{J}'_h) = f(\mathbf{J}'_h). \qquad (3.31) $$

**Remark 3.3** If the surface mesh is only composed of triangles, these two formulations are identical and correspond to the classical lowest-order RT discretization of the EFIE. Moreover, on hybrid meshes composed of simplices and polygons, the new V-BEM turns out to be an inherent coupling between the classical RT approximation on triangles and the virtual element approximation.

**Remark 3.4** As common in the standard VEM approach, a stabilization term is necessary in order to balance the lost of information due to the use of the projection, whereas in our situation the stabilization term is not mandatory. Indeed, as we will see in section 4, the weak formulation V-BEM (3.28) does not include any typical VEM stabilization while being well-posed and leading to a quasi-optimal numerical scheme.

With all the aforementioned tools and formulations, we are ready to carry out the error analysis of the virtual element approximation of the EFIE.

## 4. Main theoretical results: well-posedness and error analysis

In this section, we present the main results associated to the theoretical analysis of the V-BEM approach (3.28) proposed to discretize the EFIE. This mainly concerns the well-posedness of the new discrete formulation and the a priori error estimate, as well as some necessary intermediate results. To begin with, the following lemma provides the continuity property exhibited by the sesquilinear form $a_h$.

**Lemma 4.1 (Continuity)** The sesquilinear form $a_h$ is uniformly bounded on $\mathcal{V}_h \times \mathcal{V}_h$.

The following proposition establishes the asymptotic consistency of the V-BEM.

**Proposition 4.2 (Asymptotic consistency)** Let $(\mathcal{T}_h)_{h>0}$ be a family of surface polygonal meshes satisfying the regularity properties of section 2.2 and let $\mathcal{P}_h$ be a quasi-optimal projector from $\mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)$ to $\mathcal{V}_h$, i.e. $\mathcal{P}_h : \mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma) \to \mathcal{V}_h$, such that $\forall \mathbf{J} \in \mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)$, the following error estimate holds:

$$ \|\mathbf{J} - \mathcal{P}_h \mathbf{J}\|_{\mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)} \lesssim \inf_{\mathbf{J}'_h \in \mathcal{V}_h} \|\mathbf{J} - \mathbf{J}'_h\|_{\mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)}. \qquad (4.1) $$

Then for $\mathbf{J} \in \mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma) \cap \mathbf{H}_{\|\|}^{-1/2+\varepsilon}(\Gamma)$ for $\varepsilon > 0$, the consistency error $R_h(\mathbf{J})$ defined by

$$ R_h(\mathbf{J}) = \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{|f_h(\mathbf{J}'_h) - a_h(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|}{\|\mathbf{J}'_h\|_{\mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)}} \qquad (4.2) $$

tends to 0 when $h \to 0^+$.

The discrete inf-sup condition, uniform in mesh-size, on $a_h$ is then provided by the following theorem.
to analyze the former formulation:

$$
\begin{equation}
\begin{gathered}
\text{Find } \mathbf{J}_h \in \mathcal{V}_h \text{ such that } \forall \mathbf{J}'_h \in \mathcal{V}_h, \\
a(\mathbf{J}_h, \mathbf{J}'_h) = f(\mathbf{J}'_h).
\end{gathered}
\tag{3.31}
\end{equation}
$$

**Remark 3.3** If the surface mesh is only composed of triangles, these two formulations are identical and correspond to the classical lowest-order RT discretization of the EFIE. Moreover, on hybrid meshes composed of simplices and polygons, the new V-BEM turns out to be an inherent coupling between the classical RT approximation on triangles and the virtual element approximation.

**Remark 3.4** As common in the standard VEM approach, a stabilization term is necessary in order to balance the lost of information due to the use of the projection, whereas in our situation the stabilization term is not mandatory. Indeed, as we will see in section 4, the weak formulation V-BEM (3.28) does not include any typical VEM stabilization while being well-posed and leading to a quasi-optimal numerical scheme.

With all the aforementioned tools and formulations, we are ready to carry out the error analysis of
the virtual element approximation of the EFIE.

**4. Main theoretical results: well-posedness and error analysis**

In this section, we present the main results associated to the theoretical analysis of the V-BEM approach (3.28) proposed to discretize the EFIE. This mainly concerns the well-posedness of the new discrete formulation and the a priori error estimate, as well as some necessary intermediate results. To begin with, the following lemma provides the continuity property exhibited by the sesquilinear form $a_h$.

**Lemma 4.1 (Continuity)** The sesquilinear form $a_h$ is uniformly bounded on $V_h \times V_h$.

The following proposition establishes the asymptotic consistency of the V-BEM.

**Proposition 4.2 (Asymptotic consistency)** Let $(\mathcal{T}_h)_{h>0}$ be a family of surface polygonal meshes satisfying the regularity properties of section 2.2 and let $\mathcal{P}_h$ be a quasi-optimal projector from $\mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)$ to $\mathcal{V}_h$, i.e. $\mathcal{P}_h : \mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma) \to \mathcal{V}_h$, such that $\forall \mathbf{J} \in \mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)$, the following error estimate holds:

$$
\| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)} \lesssim \inf_{\mathbf{J}'_h \in \mathcal{V}_h} \| \mathbf{J} - \mathbf{J}'_h \|_{\mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)}. \quad (4.1)
$$

Then for $\mathbf{J} \in \mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma) \cap \mathbf{H}_\parallel^{-1/2+\varepsilon}(\Gamma)$ for $\varepsilon > 0$, the consistency error $R_h(\mathbf{J})$ defined by

$$
R_h(\mathbf{J}) = \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{|f_h(\mathbf{J}'_h) - a_h(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|}{\|\mathbf{J}'_h\|_{\mathbf{H}^{-1/2}(\text{div}_\Gamma, \Gamma)}} \quad (4.2)
$$

tends to 0 when $h \to 0^+$.

The discrete *inf-sup* condition, uniform in mesh-size, on $a_h$ is then provided by the following theorem.
**Theorem 4.3** (h-Uniform discrete inf-sup condition) Let $(T_h)_{h>0}$ be a family of surface polygonal meshes satisfying the regularity properties of section 2.2. Then, there exist $\tilde{C} > 0$ and $h_0 > 0$ such that $\forall h \in (0, h_0)$, the discrete inf-sup condition holds:

$$ \inf_{J_h \in \mathcal{V}_h} \sup_{J'_h \in \mathcal{V}_h} \frac{|a_h(J_h, J'_h)|}{\|J_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)} \|J'_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)}} \geq \tilde{C}. \quad (4.3) $$

In particular, this inf-sup condition is h-uniform in the interval $(0, h_0)$.

**Remark 4.4** Theorem 4.3 implies that the new V-BEM formulation (3.28) is well-posed.

With all above results, we can state an a priori error estimate theorem.

**Theorem 4.5** (A priori error estimate) Let $(T_h)_{h>0}$ be a family of surface polygonal meshes satisfying the regularity properties of section 2.2. Then the V-BEM scheme (3.28) verifies the following error estimate,

$$ \|\mathbf{J} - \mathbf{J}_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)} \lesssim \frac{1}{\tilde{C}} R_h(\mathbf{J}) + \inf_{\tilde{J}_h \in \mathcal{V}_h} \|\mathbf{J} - \tilde{J}_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)}, \quad (4.4) $$

where $\mathbf{J}$ is the solution of the EFIE formulation (3.9) and $\tilde{C}$ is the uniform inf-sup discrete constant (see theorem 4.3). In particular, if $\mathbf{J} \in H^{-1/2}(\text{div}_\Gamma, \Gamma) \cap H_{\|\|}^{-1/2+\varepsilon}(\Gamma)$ for $\varepsilon > 0$,

$$ \lim_{h \to 0^+} \| \mathbf{J} - \mathbf{J}_h \|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)} = 0. \quad (4.5) $$

In order to explicit the h-convergence rate of the solutions of the discretized problem (3.28) toward the continuous solution of the EFIE (3.9), we need to consider sufficiently smooth solutions and data. More precisely, from the classical regularity results (see theorem 9 of Buffa et al. (2002a)) of the EFIE solution on nonsmooth domains and for a source f corresponding to the tangential component of a regular field like the one related to the incident plane wave considered in this paper, we can consider that $\mathbf{J} \in H^\sigma(\text{div}_\Gamma, \Gamma)$ for a certain regularity parameter $\sigma \in [0, 1/2)$.

**Corollary 4.6** (Convergence rate) Under the hypothesis of theorem 4.5, if $\mathbf{J} \in H^\sigma(\text{div}_\Gamma, \Gamma)$ for $\sigma \in (0, 1/2)$ such that $V_k \mathbf{J} \in H_{\|\|}^{\sigma+1}(\Gamma)$, then the following a priori error estimate holds

$$ \|\mathbf{J} - \mathbf{J}_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)} \lesssim h^{1/2+\sigma} \left( \|\mathbf{f}\|_{H^2(\Gamma)} + \|V_k \mathbf{J}\|_{H^{\sigma+1}(\Gamma)} + \|\mathbf{J}\|_{H^\sigma(\text{div}_\Gamma, \Gamma)} \right). \quad (4.6) $$

**Remark 4.7** The estimate (4.6) shows that the V-BEM exhibits the same convergence rate as the classical lowest-order RT element scheme (see Buffa and Christiansen (2003) or Buffa and Hiptmair (2003)), but in the context of polygonal meshes.

**Remark 4.8** These results can be extended to smooth surfaces by introducing Lipschitz-isomorphisms from the polyhedral approximate surface $\Gamma_h$ to the exact one $\Gamma$. In this case, the estimate (4.6) is expected to lead to an optimal convergence rate, i.e. $h^{3/2}$, like the lowest-order RT discretization (see Christiansen (2003)).
**Theorem 4.3** (h-Uniform discrete inf-sup condition) Let $(T_h)_{h>0}$ be a family of surface polygonal meshes satisfying the regularity properties of section 2.2. Then, there exist $\tilde{C} > 0$ and $h_0 > 0$ such that $\forall h \in (0, h_0)$, the discrete inf-sup condition holds:

$$ \inf_{J_h \in \mathcal{V}_h} \sup_{J'_h \in \mathcal{V}_h} \frac{|a_h(J_h, J'_h)|}{\|J_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)} \|J'_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)}} \geq \tilde{C}. \quad (4.3) $$

In particular, this inf-sup condition is h-uniform in the interval $(0, h_0)$.

**Remark 4.4** Theorem 4.3 implies that the new V-BEM formulation (3.28) is well-posed.

With all above results, we can state an a priori error estimate theorem.

**Theorem 4.5** (A priori error estimate) Let $(T_h)_{h>0}$ be a family of surface polygonal meshes satisfying the regularity properties of section 2.2. Then the V-BEM scheme (3.28) verifies the following error estimate,

$$ \|\mathbf{J} - \mathbf{J}_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)} \lesssim \frac{1}{\tilde{C}} R_h(\mathbf{J}) + \inf_{\tilde{J}_h \in \mathcal{V}_h} \|\mathbf{J} - \tilde{J}_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)}, \quad (4.4) $$

where $\mathbf{J}$ is the solution of the EFIE formulation (3.9) and $\tilde{C}$ is the uniform inf-sup discrete constant (see theorem 4.3). In particular, if $\mathbf{J} \in H^{-1/2}(\text{div}_\Gamma, \Gamma) \cap H_{\|\|}^{-1/2+\varepsilon}(\Gamma)$ for $\varepsilon > 0$,

$$ \lim_{h \to 0^+} \| \mathbf{J} - \mathbf{J}_h \|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)} = 0. \quad (4.5) $$

In order to explicit the h-convergence rate of the solutions of the discretized problem (3.28) toward the continuous solution of the EFIE (3.9), we need to consider sufficiently smooth solutions and data. More precisely, from the classical regularity results (see theorem 9 of Buffa et al. (2002a)) of the EFIE solution on nonsmooth domains and for a source $f$ corresponding to the tangential component of a regular field like the one related to the incident plane wave considered in this paper, we can consider that $\mathbf{J} \in H^\sigma(\text{div}_\Gamma, \Gamma)$ for a certain regularity parameter $\sigma \in [0, 1/2)$.

**Corollary 4.6** (Convergence rate) Under the hypothesis of theorem 4.5, if $\mathbf{J} \in H^\sigma(\text{div}_\Gamma, \Gamma)$ for $\sigma \in (0, 1/2)$ such that $V_k \mathbf{J} \in H_{\|\|}^{\sigma+1}(\Gamma)$, then the following a priori error estimate holds

$$ \|\mathbf{J} - \mathbf{J}_h\|_{H^{-1/2}(\text{div}_\Gamma, \Gamma)} \lesssim h^{1/2+\sigma} \left( \|\mathbf{f}\|_{H^2(\Gamma)} + \|V_k \mathbf{J}\|_{H^{\sigma+1}(\Gamma)} + \|\mathbf{J}\|_{H^\sigma(\text{div}_\Gamma, \Gamma)} \right). \quad (4.6) $$

**Remark 4.7** The estimate (4.6) shows that the V-BEM exhibits the same convergence rate as the classical lowest-order RT element scheme (see Buffa and Christiansen (2003) or Buffa and Hiptmair (2003)), but in the context of polygonal meshes.

**Remark 4.8** These results can be extended to smooth surfaces by introducing Lipschitz-isomorphisms from the polyhedral approximate surface $\Gamma_h$ to the exact one $\Gamma$. In this case, the estimate (4.6) is expected to lead to an optimal convergence rate, i.e. $h^{3/2}$, like the lowest-order RT discretization (see Christiansen (2003)).
## 5. Proofs

In this section, we are going to prove all aforementioned theoretical results concerning the new V-BEM formulation. To this aim, the section is outlined into 4 parts, each of which provides comprehensive proofs and is requisite for the next part. In the rest of this section, for the sake of clarity, we use the notation $W = H^{-1/2}(\text{div}_\Gamma, \Gamma)$ and assume that $(\mathcal{T}_h)_{h>0}$ is a family of surface polygonal meshes satisfying the regularity properties of section 2.2.

### 5.1. Some technical results associated to the virtual element space $\mathcal{V}_h$

The following lemma gives an $h$-inverse estimate in the finite-dimensional virtual space $\mathcal{V}_h$. This allows us to establish a fundamental $h$-absorption property (in the sense that a given quantity $f$ is $h$-absorbable with respect to quantity $g$ if and only if $f = o(g)$ as $h \to 0^+$) associated to the $L^2$-projection operator (3.24) for $-1/2$ Sobolev norm: lemma 5.2 and corollary 5.3. Finally, a partial de Rham diagram is given in lemma 5.4.

**Lemma 5.1** (Inverse estimate) Let $(\mathcal{T}_h)_{h>0}$ be a family of meshes, then the following inverse estimate holds:

$$ \forall \mathbf{J}_h \in \mathcal{V}_h, \quad \|\mathbf{J}_h\|_0 \lesssim h^{-1/2} \|\mathbf{J}_h\|_{\mathbf{H}_{\parallel}^{-1/2}(\Gamma)}. \qquad (5.1) $$

*Proof. Step 1:* We begin by proving the following result: let $\mathbf{J}_h \in \mathcal{V}_h$ and $K \in \mathcal{T}_h$,

$$ \|\mathbf{J}_h\|_{0,K} \le h_K^{-1/2} \|\mathbf{J}_h\|_{\mathbf{H}_{\parallel}^{-1/2}(K)}. \qquad (5.2) $$

To this end, we use the inverse estimate $\|\mathbf{J}_h\|_{0,K} \lesssim h_K^{-1} \|\mathbf{J}_h\|_{\mathbf{H}_{\parallel}^{-1/2}(K)}$ that is proved in lemma 5.3 of Beirão da Veiga and Mascotto (2022) in a 2D context. Let us recall that for a planar element $K$, $\mathbf{H}_{\parallel}^{-1}(K) = \mathbf{H}_{\perp}^{-1}(K) = \mathbf{H}_t^{-1}(K)$. Next, by invoking a Sobolev space interpolation argument (see appendix B of McLean and McLean (2000)) associated to the interpolation property $\mathbf{H}_{\parallel}^{-1/2}(\Gamma) = [\mathbf{H}_{\parallel}^{-1}(\Gamma), \mathcal{L}_t^2(\Gamma)]_{\frac{1}{2}}$, we immediately obtain the result.

**Step 2:** By summing (5.2) over all elements:

$$ \|\mathbf{J}_h\|_0^2 = \sum_{K \in \mathcal{T}_h} \|\mathbf{J}_{h|_K}\|_{0,K}^2 \lesssim h^{-1} \sum_{K \in \mathcal{T}_h} \|\mathbf{J}_{h|_K}\|_{\mathbf{H}_{\parallel}^{-1/2}(K)}^2. \qquad (5.3) $$

Finally, by using the estimate $\sum_{K \in \mathcal{T}_h} \|\mathbf{J}_{h|_K}\|_{\mathbf{H}_{\parallel}^{-1/2}(K)}^2 \le \|\mathbf{J}_h\|_{\mathbf{H}_{\parallel}^{-1/2}(\Gamma)}^2$, the estimate (5.1) is obtained.

Lemma 5.1 □

**Lemma 5.2** $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that

$$ \sup_{\widetilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\| (Id - \Pi_h^0) \widetilde{\mathcal{J}}_h \|_{\mathbf{H}_{\parallel}^{-1/2}(\Gamma)}}{\| \widetilde{\mathcal{J}}_h \|_{\mathbf{H}_{\parallel}^{-1/2}(\Gamma)}} \le \beta_h. \qquad (5.4) $$
## 5. Proofs

In this section, we are going to prove all aforementioned theoretical results concerning the new V-BEM formulation. To this aim, the section is outlined into 4 parts, each of which provides comprehensive proofs and is requisite for the next part. In the rest of this section, for the sake of clarity, we use the notation $W = H^{-1/2}(\text{div}_{\Gamma}, \Gamma)$ and assume that $(\mathcal{T}_h)_{h>0}$ is a family of surface polygonal meshes satisfying the regularity properties of section 2.2.

### 5.1. Some technical results associated to the virtual element space $\mathcal{V}_h$

The following lemma gives an $h$-inverse estimate in the finite-dimensional virtual space $\mathcal{V}_h$. This allows us to establish a fundamental $h$-absorption property (in the sense that a given quantity $f$ is $h$-absorbable with respect to quantity $g$ if and only if $f = o(g)$ as $h \to 0^+$) associated to the $L^2$-projection operator (3.24) for $-1/2$ Sobolev norm: lemma 5.2 and corollary 5.3. Finally, a partial de Rham diagram is given in lemma 5.4.

**Lemma 5.1** (Inverse estimate) Let $(\mathcal{T}_h)_{h>0}$ be a family of meshes, then the following inverse estimate holds:

$$ \forall \mathbf{J}_h \in \mathcal{V}_h, \quad \|\mathbf{J}_h\|_0 \lesssim h^{-1/2} \|\mathbf{J}_h\|_{\mathbf{H}_{\parallel}^{-1/2}(\Gamma)}. \qquad (5.1) $$

*Proof. Step 1:* We begin by proving the following result: let $\mathbf{J}_h \in \mathcal{V}_h$ and $K \in \mathcal{T}_h$,

$$ \|\mathbf{J}_h\|_{0,K} \le h_K^{-1/2} \|\mathbf{J}_h\|_{\mathbf{H}_{\parallel}^{-1/2}(K)}. \qquad (5.2) $$

To this end, we use the inverse estimate $\|\mathbf{J}_h\|_{0,K} \lesssim h_K^{-1} \|\mathbf{J}_h\|_{\mathbf{H}_{\parallel}^{-1/2}(K)}$ that is proved in lemma 5.3 of Beirão da Veiga and Mascotto (2022) in a 2D context. Let us recall that for a planar element $K$, $\mathbf{H}_{\parallel}^{-1}(K) = \mathbf{H}_{\perp}^{-1}(K) = \mathbf{H}_t^{-1}(K)$. Next, by invoking a Sobolev space interpolation argument (see appendix B of McLean and McLean (2000)) associated to the interpolation property $\mathbf{H}_{\parallel}^{-1/2}(\Gamma) = [\mathbf{H}_{\parallel}^{-1}(\Gamma), \mathcal{L}_t^2(\Gamma)]_{\frac{1}{2}}$, we immediately obtain the result.

**Step 2:** By summing (5.2) over all elements:

$$ \|\mathbf{J}_h\|_0^2 = \sum_{K \in \mathcal{T}_h} \|\mathbf{J}_{h|_K}\|_{0,K}^2 \lesssim h^{-1} \sum_{K \in \mathcal{T}_h} \|\mathbf{J}_{h|_K}\|_{\mathbf{H}_{\parallel}^{-1/2}(K)}^2. \qquad (5.3) $$

Finally, by using the estimate $\sum_{K \in \mathcal{T}_h} \|\mathbf{J}_{h|_K}\|_{\mathbf{H}_{\parallel}^{-1/2}(K)}^2 \le \|\mathbf{J}_h\|_{\mathbf{H}_{\parallel}^{-1/2}(\Gamma)}^2$, the estimate (5.1) is obtained.

Lemma 5.1 □

**Lemma 5.2** $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that

$$ \sup_{\widetilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\| (Id - \Pi_h^0) \widetilde{\mathcal{J}}_h \|_{\mathbf{H}_{\parallel}^{-1/2}(\Gamma)}}{\| \widetilde{\mathcal{J}}_h \|_{\mathbf{H}_{\parallel}^{-1/2}(\Gamma)}} \le \beta_h. \qquad (5.4) $$
*Proof.* From the definition of the dual norm, we have:

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} = \sup_{\tilde{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, \hat{\mathcal{J}} \right\rangle_{\parallel, -1/2, 1/2}}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.5) $$

Since $(Id - \Pi_h^0) \tilde{\mathcal{J}}_h \in L_t^2(\Gamma)$ and $\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma) \subset L_t^2(\Gamma)$, the duality pairing of (5.5) can be written as the $L^2$-scalar product:

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} = \sup_{\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, \hat{\mathcal{J}} \right\rangle_0}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.6) $$

Furthermore, from the definition of the projection operator $\Pi_h^0$, we have that

$$ \left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, q \right\rangle_0 = 0 \quad \forall q \in \mathbb{P}_{1,t}(\mathcal{T}_h), \quad (5.7) $$

which implies that $\forall q \in \mathbb{P}_{1,t}(\mathcal{T}_h)$:

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} = \sup_{\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, \hat{\mathcal{J}} - q \right\rangle_0}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.8) $$

In particular, by taking $q = \Pi_h^0 \hat{\mathcal{J}}$:

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} = \sup_{\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, (Id - \Pi_h^0) \hat{\mathcal{J}} \right\rangle_0}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.9) $$

Then, using Cauchy-Schwarz inequality (CS) yields

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} \le \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_0 \sup_{\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\| (Id - \Pi_h^0) \hat{\mathcal{J}} \right\|_0}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.10) $$

Besides, if $\hat{\mathcal{J}} \in H_{\parallel}^1(\Gamma)$, then the projection's error is upper bounded:

$$ \left\| (Id - \Pi_h^0) \hat{\mathcal{J}} \right\|_0 \lesssim h \left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^1(\Gamma)}, \quad (5.11) $$

and by using an interpolation argument associated to $H_{\parallel}^{1/2}(\Gamma) := [L_t^2(\Gamma), H_{\parallel}^1(\Gamma)]_{1/2}$, we obtain the following error estimate for $\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)$:

$$ \left\| (Id - \Pi_h^0) \hat{\mathcal{J}} \right\|_0 \lesssim h^{1/2} \left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}, \quad (5.12) $$
*Proof.* From the definition of the dual norm, we have:

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} = \sup_{\tilde{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, \hat{\mathcal{J}} \right\rangle_{\parallel, -1/2, 1/2}}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.5) $$

Since $(Id - \Pi_h^0) \tilde{\mathcal{J}}_h \in L_t^2(\Gamma)$ and $\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma) \subset L_t^2(\Gamma)$, the duality pairing of (5.5) can be written as the $L^2$-scalar product:

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} = \sup_{\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, \hat{\mathcal{J}} \right\rangle_0}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.6) $$

Furthermore, from the definition of the projection operator $\Pi_h^0$, we have that

$$ \left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, q \right\rangle_0 = 0 \quad \forall q \in \mathbb{P}_{1,t}(\mathcal{T}_h), \quad (5.7) $$

which implies that $\forall q \in \mathbb{P}_{1,t}(\mathcal{T}_h)$:

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} = \sup_{\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, \hat{\mathcal{J}} - q \right\rangle_0}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.8) $$

In particular, by taking $q = \Pi_h^0 \hat{\mathcal{J}}$:

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} = \sup_{\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\langle (Id - \Pi_h^0) \tilde{\mathcal{J}}_h, (Id - \Pi_h^0) \hat{\mathcal{J}} \right\rangle_0}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.9) $$

Then, using Cauchy-Schwarz inequality (CS) yields

$$ \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_{H_{\parallel}^{-1/2}(\Gamma)} \le \left\| (Id - \Pi_h^0) \tilde{\mathcal{J}}_h \right\|_0 \sup_{\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)} \frac{\left\| (Id - \Pi_h^0) \hat{\mathcal{J}} \right\|_0}{\left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}}. \quad (5.10) $$

Besides, if $\hat{\mathcal{J}} \in H_{\parallel}^1(\Gamma)$, then the projection's error is upper bounded:

$$ \left\| (Id - \Pi_h^0) \hat{\mathcal{J}} \right\|_0 \lesssim h \left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^1(\Gamma)}, \quad (5.11) $$

and by using an interpolation argument associated to $H_{\parallel}^{1/2}(\Gamma) := [L_t^2(\Gamma), H_{\parallel}^1(\Gamma)]_{1/2}$, we obtain the following error estimate for $\hat{\mathcal{J}} \in H_{\parallel}^{1/2}(\Gamma)$:

$$ \left\| (Id - \Pi_h^0) \hat{\mathcal{J}} \right\|_0 \lesssim h^{1/2} \left\| \hat{\mathcal{J}} \right\|_{H_{\parallel}^{1/2}(\Gamma)}, \quad (5.12) $$
By inserting (5.12) into (5.10), we get the inequality:

$$ \|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)} \lesssim h^{1/2} \|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_0, \quad (5.13) $$

which means that

$$ \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \lesssim h^{1/2} \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_0}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}. \quad (5.14) $$

By using the inverse estimate from lemma 5.1, we have that

$$ \frac{h^{1/2}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \lesssim \frac{1}{\|\tilde{\mathcal{J}}_h\|_0}, \quad (5.15) $$

thus:

$$ \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \lesssim \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_0}{\|\tilde{\mathcal{J}}_h\|_0} \quad (5.16) \\ \leq \|Id - \Pi_h^0\|_{h \to 0^+} = 0 $$

Then, $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that

$$ \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \leq \beta_h. \quad (5.17) $$

Lemma 5.2 □

Therefore, we can deduce the following corollary.

**Corollary 5.3** $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that $\forall \mathcal{J}_h \in \mathcal{V}_h$,

$$ \|(Id - \Pi_h^0) \mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)} \leq \beta_h \|\mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}. \quad (5.18) $$

*Proof.* With lemma 5.2, $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that

$$ \frac{\|(Id - \Pi_h^0) \mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \leq \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \leq \beta_h, \quad (5.19) $$

$$ \Rightarrow \|(Id - \Pi_h^0) \mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)} \leq \beta_h \|\mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}. \quad (5.20) $$

Corollary 5.3 □
By inserting (5.12) into (5.10), we get the inequality:

$$ \|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)} \lesssim h^{1/2} \|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_0, \quad (5.13) $$

which means that

$$ \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \lesssim h^{1/2} \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_0}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}. \quad (5.14) $$

By using the inverse estimate from lemma 5.1, we have that

$$ \frac{h^{1/2}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \lesssim \frac{1}{\|\tilde{\mathcal{J}}_h\|_0}, \quad (5.15) $$

thus:

$$ \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \lesssim \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_0}{\|\tilde{\mathcal{J}}_h\|_0} \quad (5.16) \\ \leq \|Id - \Pi_h^0\|_{h \to 0^+} = 0 $$

Then, $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that

$$ \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \leq \beta_h. \quad (5.17) $$

Lemma 5.2 □

Therefore, we can deduce the following corollary.

**Corollary 5.3** $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that $\forall \mathcal{J}_h \in \mathcal{V}_h$,

$$ \|(Id - \Pi_h^0) \mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)} \leq \beta_h \|\mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}. \quad (5.18) $$

*Proof.* With lemma 5.2, $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that

$$ \frac{\|(Id - \Pi_h^0) \mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \leq \sup_{\tilde{\mathcal{J}}_h \in \mathcal{V}_h} \frac{\|(Id - \Pi_h^0) \tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}}{\|\tilde{\mathcal{J}}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}} \leq \beta_h, \quad (5.19) $$

$$ \Rightarrow \|(Id - \Pi_h^0) \mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)} \leq \beta_h \|\mathcal{J}_h\|_{\mathcal{H}_{\|}^{-1/2}(\Gamma)}. \quad (5.20) $$

Corollary 5.3 □
The next lemma provides the commuting diagram property of the VEM interpolation operator.

**Lemma 5.4** Let $P_h^X$ and $Q_h$ be the VEM interpolation operator on $V_h$ and the $L^2$-projection into $Q_h = \mathbb{P}_0(\mathcal{T}_h)$, respectively. Then the diagram 1 is commuting, i.e. $\forall \mathbf{v} \in \mathbf{H}(\mathrm{div}_{\Gamma}, \Gamma)$, $\nabla_{\Gamma} \cdot (P_h^X \mathbf{v}) = Q_h(\nabla_{\Gamma} \cdot \mathbf{v})$.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Commuting diagram.</figcaption></figure>

*Proof.* It is sufficient to locally prove the commuting property on each element $K \in \mathcal{T}_h$. Let $\mathbf{v} \in \mathbf{H}(\mathrm{div}_{\Gamma}, \Gamma)$, the VEM interpolation is defined on the basis of the degrees of freedom:

$$P_h^X \mathbf{v}|_K = \sum_{e \in \mathcal{E}_K} \left( \int_e \mathbf{v}|_K \cdot \mathbf{n}_e^K ds \right) \boldsymbol{\varphi}_e^K. \qquad (5.21)$$

The local characterization of the basis functions immediately leads to

$$\nabla_{\Gamma} \cdot (P_h^X \mathbf{v}|_K) = \sum_{e \in \mathcal{E}_K} \left( \int_e \mathbf{v}|_K \cdot \mathbf{n}_e^K ds \right) \nabla_{\Gamma} \cdot \boldsymbol{\varphi}_e^K, \qquad (5.22)$$

$$= \frac{1}{|K|} \sum_{e \in \mathcal{E}_K} \left( \int_e \mathbf{v}|_K \cdot \mathbf{n}_e^K ds \right). \qquad (5.23)$$

Thus, we have that $\nabla_{\Gamma} \cdot (P_h^X \mathbf{v}|_K) \in \mathbb{P}_0(\mathcal{T}_h)$. Now, from the definition of the $L^2$-projection $Q_h$, we have:

$$\int_K Q_h (\nabla_{\Gamma} \cdot \mathbf{v}) d\gamma = \int_K \nabla_{\Gamma} \cdot \mathbf{v} d\gamma, \qquad (5.24)$$

$$\Rightarrow |K| Q_h (\nabla_{\Gamma} \cdot \mathbf{v})|_K = \int_{\partial K} \mathbf{v}|_K \cdot \mathbf{n}_{\partial K} ds, \qquad (5.25)$$

$$\Rightarrow Q_h (\nabla_{\Gamma} \cdot \mathbf{v})|_K = \frac{1}{|K|} \sum_{e \in \mathcal{E}_K} \left( \int_e \mathbf{v}|_K \cdot \mathbf{n}_e^K ds \right). \qquad (5.26)$$

Therefore,

$$\nabla_{\Gamma} \cdot (P_h^X \mathbf{v}|_K) = Q_h (\nabla_{\Gamma} \cdot \mathbf{v})|_K. \qquad (5.27)$$

Lemma 5.4 □
The next lemma provides the commuting diagram property of the VEM interpolation operator.

**Lemma 5.4** Let $P_h^X$ and $Q_h$ be the VEM interpolation operator on $V_h$ and the $L^2$-projection into $Q_h = \mathbb{P}_0(T_h)$, respectively. Then the diagram 1 is commuting, i.e. $\forall \mathbf{v} \in \mathbf{H}(\mathrm{div}_\Gamma, \Gamma)$, $\nabla_\Gamma \cdot (P_h^X \mathbf{v}) = Q_h(\nabla_\Gamma \cdot \mathbf{v})$.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Commuting diagram.</figcaption></figure>

*Proof.* It is sufficient to locally prove the commuting property on each element $K \in T_h$. Let $\mathbf{v} \in \mathbf{H}(\mathrm{div}_\Gamma, \Gamma)$, the VEM interpolation is defined on the basis of the degrees of freedom:

$$P_h^X \mathbf{v}|_K = \sum_{e \in \mathcal{E}_K} \left( \int_e \mathbf{v}|_K \cdot \mathbf{n}_e^K \, ds \right) \boldsymbol{\varphi}_e^K. \qquad (5.21)$$

The local characterization of the basis functions immediately leads to

$$\nabla_\Gamma \cdot (P_h^X \mathbf{v}|_K) = \sum_{e \in \mathcal{E}_K} \left( \int_e \mathbf{v}|_K \cdot \mathbf{n}_e^K \, ds \right) \nabla_\Gamma \cdot \boldsymbol{\varphi}_e^K, \qquad (5.22)$$

$$= \frac{1}{|K|} \sum_{e \in \mathcal{E}_K} \left( \int_e \mathbf{v}|_K \cdot \mathbf{n}_e^K \, ds \right). \qquad (5.23)$$

Thus, we have that $\nabla_\Gamma \cdot (P_h^X \mathbf{v}|_K) \in \mathbb{P}_0(T_h)$. Now, from the definition of the $L^2$-projection $Q_h$, we have:

$$\int_K Q_h (\nabla_\Gamma \cdot \mathbf{v}) d\gamma = \int_K \nabla_\Gamma \cdot \mathbf{v} d\gamma, \qquad (5.24)$$

$$\Rightarrow |K| Q_h (\nabla_\Gamma \cdot \mathbf{v})|_K = \int_{\partial K} \mathbf{v}|_K \cdot \mathbf{n}_{\partial K} ds, \qquad (5.25)$$

$$\Rightarrow Q_h (\nabla_\Gamma \cdot \mathbf{v})|_K = \frac{1}{|K|} \sum_{e \in \mathcal{E}_K} \left( \int_e \mathbf{v}|_K \cdot \mathbf{n}_e^K ds \right). \qquad (5.26)$$

Therefore,

$$\nabla_\Gamma \cdot (P_h^X \mathbf{v}|_K) = Q_h (\nabla_\Gamma \cdot \mathbf{v})|_K. \qquad (5.27)$$

Lemma 5.4 □
## 5.2. Discrete inf-sup condition of the ideal V-BEM

In this part, we establish the discrete inf-sup condition of the ideal V-BEM (3.31), which does not involve any VEM projection operator. To this aim, we follow the methodology proposed in Buffa and Hiptmair (2003) in the context of boundary element discretization of the EFIE by RT finite elements. Let us recall that the weak formulation of the EFIE (3.9) can be split into a coercive form and a compact one. In particular, from a Helmholtz-type decomposition of any tangential vector field $J \in W$, i.e. $J = \nabla_{\Gamma}\varphi + n \times \nabla_{\Gamma}\psi =: \mathcal{R}_{\Gamma}(J) + \mathcal{Z}_{\Gamma}(J)$, we can define a test-function $J' = (\chi_{\Gamma} + T)\bar{J}$, where $\chi_{\Gamma}J := \mathcal{R}_{\Gamma}(J) - \mathcal{Z}_{\Gamma}(J)$ and $T: W \to W$ is a compact operator, with $\bar{J}$ being the complex conjugate of $J$, such that, $\exists \alpha > 0$,

$$ \forall J \in W, \quad \frac{|a(J, J')|}{\|J\|_{W}} \geq \alpha \|J'\|_{W}. \qquad (5.28) $$

The analysis of the discrete formulation involving RT discretization (Buffa and Hiptmair, 2003) is carried out by mimicking the proof of the above continuous inf-sup condition, i.e. making use of a discrete Helmholtz-type decomposition, $J'_h = (P_h^X \circ \chi_\Gamma + P_h^T \circ T) \bar{J}_h$ where $P_h^X$ is the local RT interpolation operator and $P_h^T$ is the W-orthogonal projector. Now, when it comes to applying this methodology to the ideal V-BEM scheme, $P_h^X$ is instead intended to be the VEM interpolation operator. The next two lemmas give h-absorption fundamental results associated to both operators T and $\chi_\Gamma$ in the virtual element setting.

**Lemma 5.5 (h-Absorption for T)** $\exists \varepsilon_h > 0$ such that, $\lim_{h \to 0^+} \varepsilon_h = 0, \forall J \in W$,

$$ \|(Id - P_h^T) \circ T(J)\|_W \leq \varepsilon_h \|J\|_W. \qquad (5.29) $$

*Proof.* The proof is the same as for the equation (66) of Buffa and Hiptmair (2003), which comes from corollary 10.5 of Kress (2014). Lemma 5.5 □

**Lemma 5.6 (h-Absorption for $\chi_\Gamma$)** For all $J_h \in V_h$, the following estimate holds:

$$ |a(J_h, (Id - P_h^X) \circ \chi_\Gamma (\bar{J}_h))| \lesssim h^{1/2} \|J_h\|_W^2. \qquad (5.30) $$

*Proof.* Since $(Id - P_h^X) \circ \chi_\Gamma (\bar{J}_h) = \overline{(Id - P_h^X) \circ \chi_\Gamma (J_h)}$ and from the continuity of the sesquilinear form a, we have that

$$ |a(J_h, (Id - P_h^X) \circ \chi_\Gamma (\bar{J}_h))| \lesssim \|J_h\|_W \|(Id - P_h^X) \circ \chi_\Gamma (J_h)\|_W. \qquad (5.31) $$

From the commuting diagram 5.4 and the property $\nabla_{\Gamma} \cdot \chi_{\Gamma}(J_h) = \nabla_{\Gamma} \cdot J_h$, we can derive:

$$ \nabla_{\Gamma} \cdot (Id - P_h^X) \circ \mathcal{R}^{\Gamma}(J_h) = \nabla_{\Gamma} \cdot \mathcal{R}^{\Gamma}(J_h) - Q_h \circ \nabla_{\Gamma} \cdot (\mathcal{R}^{\Gamma}(J_h)), \qquad (5.32) $$

$$ = \nabla_{\Gamma} \cdot J_h - Q_h \circ \nabla_{\Gamma} \cdot (J_h) = 0. \qquad (5.33) $$

That implies:

$$ \|(Id - P_h^X) \circ \chi_\Gamma (J_h)\|_W = \|(Id - P_h^X) \circ \mathcal{R}^{\Gamma}(J_h)\|_{H_{\parallel}^{-1/2}(\Gamma)}. \qquad (5.34) $$
## 5.2. Discrete inf-sup condition of the ideal V-BEM

In this part, we establish the discrete inf-sup condition of the ideal V-BEM (3.31), which does not involve any VEM projection operator. To this aim, we follow the methodology proposed in Buffa and Hiptmair (2003) in the context of boundary element discretization of the EFIE by RT finite elements. Let us recall that the weak formulation of the EFIE (3.9) can be split into a coercive form and a compact one. In particular, from a Helmholtz-type decomposition of any tangential vector field $J \in W$, i.e. $J = \nabla_{\Gamma}\varphi + n \times \nabla_{\Gamma}\psi =: \mathcal{R}_{\Gamma}(J) + \mathcal{Z}_{\Gamma}(J)$, we can define a test-function $J' = (\chi_{\Gamma} + T)\bar{J}$, where $\chi_{\Gamma}J := \mathcal{R}_{\Gamma}(J) - \mathcal{Z}_{\Gamma}(J)$ and $T: W \to W$ is a compact operator, with $\bar{J}$ being the complex conjugate of $J$, such that, $\exists \alpha > 0$,

$$ \forall J \in W, \quad \frac{|a(J, J')|}{\|J\|_{W}} \ge \alpha \|J'\|_{W}. \qquad (5.28) $$

The analysis of the discrete formulation involving RT discretization (Buffa and Hiptmair, 2003) is carried out by mimicking the proof of the above continuous inf-sup condition, i.e. making use of a discrete Helmholtz-type decomposition, $J'_h = (P_h^X \circ \chi_\Gamma + P_h^T \circ T) \bar{J}_h$ where $P_h^X$ is the local RT interpolation operator and $P_h^T$ is the W-orthogonal projector. Now, when it comes to applying this methodology to the ideal V-BEM scheme, $P_h^X$ is instead intended to be the VEM interpolation operator. The next two lemmas give h-absorption fundamental results associated to both operators $T$ and $\chi_\Gamma$ in the virtual element setting.

**Lemma 5.5** (h-Absorption for $T$) $\exists \varepsilon_h > 0$ such that, $\lim_{h \to 0^+} \varepsilon_h = 0, \forall J \in W$,

$$ \|(Id - P_h^T) \circ T(J)\|_{W} \le \varepsilon_h \|J\|_{W}. \qquad (5.29) $$

*Proof.* The proof is the same as for the equation (66) of Buffa and Hiptmair (2003), which comes from corollary 10.5 of Kress (2014). Lemma 5.5 □

**Lemma 5.6** (h-Absorption for $\chi_\Gamma$) For all $J_h \in V_h$, the following estimate holds:

$$ |a(J_h, (Id - P_h^X) \circ \chi_\Gamma (\bar{J}_h))| \lesssim h^{1/2} \|J_h\|_W^2. \qquad (5.30) $$

*Proof.* Since $(Id - P_h^X) \circ \chi_\Gamma (\bar{J}_h) = (Id - P_h^X) \circ \chi_\Gamma (J_h)$ and from the continuity of the sesquilinear form $a$, we have that

$$ |a(J_h, (Id - P_h^X) \circ \chi_\Gamma (\bar{J}_h))| \lesssim \|J_h\|_W \|(Id - P_h^X) \circ \chi_\Gamma (J_h)\|_W. \qquad (5.31) $$

From the commuting diagram 5.4 and the property $\nabla_\Gamma \cdot \chi_\Gamma(J_h) = \nabla_\Gamma \cdot J_h$, we can derive:

$$ \nabla_{\Gamma} \cdot (Id - P_{h}^{X}) \circ \mathcal{R}^{\Gamma}(J_{h}) = \nabla_{\Gamma} \cdot \mathcal{R}^{\Gamma}(J_{h}) - Q_{h} \circ \nabla_{\Gamma} \cdot (\mathcal{R}^{\Gamma}(J_{h})), \qquad (5.32) $$

$$ = \nabla_{\Gamma} \cdot J_{h} - Q_{h} \circ \nabla_{\Gamma} \cdot (J_{h}) = 0. \qquad (5.33) $$

That implies:

$$ \|(Id - P_h^X) \circ \chi_\Gamma (J_h)\|_W = \|(Id - P_h^X) \circ \mathcal{R}^\Gamma (J_h)\|_{H_{\|}^{-1/2}(\Gamma)}. \qquad (5.34) $$
Since $\mathcal{R}^{\Gamma}(\mathbf{J}_h) \in \mathbf{H}_{\parallel}^{1/2}(\Gamma)$ and $\nabla_{\Gamma} \cdot \mathcal{R}^{\Gamma}(\mathbf{J}_h) \in \mathcal{Q}_h$, from a straightforward extension of lemma 16 and a direct application of lemma 2 of Buffa and Hiptmair (2003), we have:

$$ \|(Id - P_h^X) \circ \mathcal{R}^{\Gamma}(\mathbf{J}_h)\|_{\mathbf{L}_{\parallel}^2(\Gamma)} \lesssim h^{1/2} \|\mathcal{R}^{\Gamma}(\mathbf{J}_h)\|_{\mathbf{H}_{\parallel}^{1/2}(\Gamma)}, \quad (5.35) $$

$$ \lesssim h^{1/2} \|\nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{-1/2}(\Gamma)}. \quad (5.36) $$

By using (5.36) in (5.34) and the $\mathbf{W}$-norm, (5.34) leads to

$$ \|(Id - P_h^X) \circ \chi_{\Gamma}(\mathbf{J}_h)\|_{\mathbf{W}} \lesssim h^{1/2} \|\nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{-1/2}(\Gamma)}, \quad (5.37) $$

$$ \lesssim h^{1/2} \|\mathbf{J}_h\|_{\mathbf{W}}. \quad (5.38) $$

Finally, combining (5.31) and (5.38), we get the result.

Lemma 5.6 □

Now, we can prove the discrete inf-sup condition on the sesquilinear form $a$ in (3.31).

**Theorem 5.7** ($h$-Uniform discrete inf-sup condition of the ideal V-BEM) Let $(\mathcal{T}_h)_{h>0}$ be a family of surface polygonal meshes. Then, it exists $C > 0$ and $h_0 > 0$ such that $\forall h \in (0, h_0)$,

$$ \inf_{\mathbf{J}_h \in \mathcal{V}_h} \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{|a(\mathbf{J}_h, \mathbf{J}'_h)|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}'_h\|_{\mathbf{W}}} \geq C. \quad (5.39) $$

In particular, this inf-sup condition is $h$-uniform in the interval $(0, h_0)$.

*Proof.* On the basis of the methodology of Buffa and Hiptmair (2003), we choose as a test-function the following candidate:

$$ \mathbf{J}'_h = P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}). \quad (5.40) $$

Therefore, we can write that

$$ \begin{aligned} |a(\mathbf{J}_h, \mathbf{J}'_h)| &= |a(\mathbf{J}_h, P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}))|, \\ &= |a(\mathbf{J}_h, (\chi_{\Gamma} + \mathbf{T})(\overline{\mathbf{J}_h}))| \\ &\quad - |a(\mathbf{J}_h, ((Id - P_h^X) \circ \chi_{\Gamma} + (Id - P_h^T) \circ \mathbf{T})(\overline{\mathbf{J}_h}))|. \end{aligned} \quad (5.41) $$

By using lemma 10 from Buffa and Hiptmair (2003) and the continuity of $a$, $\exists C_G > 0$ and $\exists \tilde{C} > 0$ such that

$$ \begin{aligned} |a(\mathbf{J}_h, P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}))| &\geq C_G \|\mathbf{J}_h\|_{\mathbf{W}}^2 - |a(\mathbf{J}_h, (Id - P_h^X) \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}))| \\ &\quad - \tilde{C} \|\mathbf{J}_h\|_{\mathbf{W}} \|(Id - P_h^T) \circ \mathbf{T}(\overline{\mathbf{J}_h})\|_{\mathbf{W}}. \end{aligned} \quad (5.42) $$

Combining lemmas 5.5 and 5.6 in (5.42) yields: $\exists C > 0$ (independent of $h$),

$$ |a(\mathbf{J}_h, P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}))| \geq C_G \|\mathbf{J}_h\|_{\mathbf{W}}^2 - C h^{1/2} \|\mathbf{J}_h\|_{\mathbf{W}}^2 - \tilde{C} \varepsilon_h \|\mathbf{J}_h\|_{\mathbf{W}}^2. \quad (5.43) $$

Thus, $\exists h_0 > 0$ such that $\forall h < h_0$, $-\frac{1}{C_G}(\tilde{C}\varepsilon_h + Ch^{1/2}) > 1/2$, and the following estimate holds:

$$ |a(\mathbf{J}_h, P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}))| \geq \frac{C_G}{2} \|\mathbf{J}_h\|_{\mathbf{W}}^2. \quad (5.44) $$

Finally, by using the estimate $\|\mathbf{J}'_h\|_{\mathbf{W}} \lesssim \|\mathbf{J}_h\|_{\mathbf{W}}$, (5.44) leads to the result.

Theorem 5.7 □
Since $\mathcal{R}^{\Gamma}(\mathbf{J}_h) \in \mathbf{H}_{\parallel}^{1/2}(\Gamma)$ and $\nabla_{\Gamma} \cdot \mathcal{R}^{\Gamma}(\mathbf{J}_h) \in \mathcal{Q}_h$, from a straightforward extension of lemma 16 and a direct application of lemma 2 of Buffa and Hiptmair (2003), we have:

$$ \|(Id - P_h^X) \circ \mathcal{R}^{\Gamma}(\mathbf{J}_h)\|_{\mathbf{L}_{\parallel}^2(\Gamma)} \lesssim h^{1/2} \|\mathcal{R}^{\Gamma}(\mathbf{J}_h)\|_{\mathbf{H}_{\parallel}^{1/2}(\Gamma)}, \quad (5.35) $$

$$ \lesssim h^{1/2} \|\nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{-1/2}(\Gamma)}. \quad (5.36) $$

By using (5.36) in (5.34) and the $\mathbf{W}$-norm, (5.34) leads to

$$ \|(Id - P_h^X) \circ \chi_{\Gamma}(\mathbf{J}_h)\|_{\mathbf{W}} \lesssim h^{1/2} \|\nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{-1/2}(\Gamma)}, \quad (5.37) $$

$$ \lesssim h^{1/2} \|\mathbf{J}_h\|_{\mathbf{W}}. \quad (5.38) $$

Finally, combining (5.31) and (5.38), we get the result.

Lemma 5.6 □

Now, we can prove the discrete inf-sup condition on the sesquilinear form $a$ in (3.31).

**Theorem 5.7** (h-Uniform discrete inf-sup condition of the ideal V-BEM) Let $(\mathcal{T}_h)_{h>0}$ be a family of surface polygonal meshes. Then, it exists $C > 0$ and $h_0 > 0$ such that $\forall h \in (0, h_0)$,

$$ \inf_{\mathbf{J}_h \in \mathcal{V}_h} \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{|a(\mathbf{J}_h, \mathbf{J}'_h)|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}'_h\|_{\mathbf{W}}} \geq C. \quad (5.39) $$

In particular, this inf-sup condition is h-uniform in the interval $(0, h_0)$.

*Proof.* On the basis of the methodology of Buffa and Hiptmair (2003), we choose as a test-function the following candidate:

$$ \mathbf{J}'_h = P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}). \quad (5.40) $$

Therefore, we can write that

$$ \begin{aligned} |a(\mathbf{J}_h, \mathbf{J}'_h)| &= |a(\mathbf{J}_h, P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}))|, \\ &= |a(\mathbf{J}_h, (\chi_{\Gamma} + \mathbf{T})(\overline{\mathbf{J}_h}))| \\ &\quad - |a(\mathbf{J}_h, ((Id - P_h^X) \circ \chi_{\Gamma} + (Id - P_h^T) \circ \mathbf{T})(\overline{\mathbf{J}_h}))|. \end{aligned} \quad (5.41) $$

By using lemma 10 from Buffa and Hiptmair (2003) and the continuity of $a$, $\exists C_G > 0$ and $\exists \tilde{C} > 0$ such that

$$ \begin{aligned} |a(\mathbf{J}_h, P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}))| &\geq C_G \|\mathbf{J}_h\|_{\mathbf{W}}^2 - |a(\mathbf{J}_h, (Id - P_h^X) \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}))| \\ &\quad - \tilde{C} \|\mathbf{J}_h\|_{\mathbf{W}} \|(Id - P_h^T) \circ \mathbf{T}(\overline{\mathbf{J}_h})\|_{\mathbf{W}}. \end{aligned} \quad (5.42) $$

Combining lemmas 5.5 and 5.6 in (5.42) yields: $\exists C > 0$ (independent of $h$),

$$ |a(\mathbf{J}_h, P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}))| \geq C_G \|\mathbf{J}_h\|_{\mathbf{W}}^2 - C h^{1/2} \|\mathbf{J}_h\|_{\mathbf{W}}^2 - \tilde{C} \varepsilon_h \|\mathbf{J}_h\|_{\mathbf{W}}^2. \quad (5.43) $$

Thus, $\exists h_0 > 0$ such that $\forall h < h_0$, $-\frac{1}{C_G} (\tilde{C} \varepsilon_h + Ch^{1/2}) > 1/2$, and the following estimate holds:

$$ |a(\mathbf{J}_h, P_h^X \circ \chi_{\Gamma}(\overline{\mathbf{J}_h}) + P_h^T \circ T(\overline{\mathbf{J}_h}))| \geq \frac{C_G}{2} \|\mathbf{J}_h\|_{\mathbf{W}}^2. \quad (5.44) $$

Finally, by using the estimate $\|\mathbf{J}'_h\|_{\mathbf{W}} \lesssim \|\mathbf{J}_h\|_{\mathbf{W}}$, (5.44) leads to the result.

Theorem 5.7 □
### 5.3. Discrete inf-sup condition of the V-BEM

In order to prove a h-uniform discrete inf-sup condition on the perturbed sesquilinear form $a_h$ related to the V-BEM (3.28), we rely on both the methodology from Buffa and Hiptmair (2003), (i.e. (5.40)) and the stability result for the ideal V-BEM (see theorem 5.7).

To begin with, we need to derive a technical lemma by using the tools from section 5.1 that provides a h-absorption result on the difference of the sesquilinear form a (3.31) and its perturbation $a_h$ (3.28) with respect to the W-norm.

**Lemma 5.8** $\exists \varepsilon_h \underset{h \to 0^+}{\longrightarrow} 0$ such that, $\forall \mathbf{J}_h \in \mathcal{V}_h$,

$$ \mathcal{E}_h = \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|a(\mathbf{J}_h, \widetilde{\mathbf{J}}_h) - a_h(\mathbf{J}_h, \widetilde{\mathbf{J}}_h)|}{\|\widetilde{\mathbf{J}}_h\|_W} \le \varepsilon_h \|\mathbf{J}_h\|_W, \quad (5.45) $$

*Proof.* Since the difference of the divergence terms vanishes and by the characterization of W-norm, the following estimate holds:

$$ \mathcal{E}_h = \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|a(\mathbf{J}_h, \widetilde{\mathbf{J}}_h) - a_h(\mathbf{J}_h, \widetilde{\mathbf{J}}_h)|}{\|\widetilde{\mathbf{J}}_h\|_W} \quad (5.46) $$

$$ = \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|\langle \mathbf{V}_k \mathbf{J}_h, \widetilde{\mathbf{J}}_h \rangle_{\|} - \langle \mathbf{V}_k \Pi_h^0 \mathbf{J}_h, \Pi_h^0 \widetilde{\mathbf{J}}_h \rangle_{\|}|}{\|\widetilde{\mathbf{J}}_h\|_W}, $$

$$ = \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|\langle \mathbf{V}_k ((Id - \Pi_h^0) \mathbf{J}_h), \widetilde{\mathbf{J}}_h \rangle_{\|} + \langle \mathbf{V}_k \Pi_h^0 \mathbf{J}_h, ((Id - \Pi_h^0) \widetilde{\mathbf{J}}_h) \rangle_{\|}|}{\|\widetilde{\mathbf{J}}_h\|_W}, $$

$$ \le \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{\overbrace{\langle \mathbf{V}_k ((Id - \Pi_h^0) \mathbf{J}_h), \widetilde{\mathbf{J}}_h \rangle_{\|}}^{A_1} + \overbrace{\langle \mathbf{V}_k \Pi_h^0 \mathbf{J}_h, ((Id - \Pi_h^0) \widetilde{\mathbf{J}}_h) \rangle_{\|}}^{A_2}}{\|\widetilde{\mathbf{J}}_h\|_{H_{\|}^{-1/2}(\Gamma)}}. \quad (5.47) $$

Since $\mathbf{V}_k$ is continuous from $H_{\|}^{-1/2}(\Gamma)$ to $H_{\|}^{1/2}(\Gamma)$ (corollary 3 of Buffa and Hiptmair (2003)), the terms $A_1$ and $A_2$ of (5.47) can be upper bounded as follows:

$$ \begin{aligned} A_1 &\lesssim \|\mathbf{V}_k((Id - \Pi_h^0)\mathbf{J}_h)\|_{H_{\|}^{1/2}(\Gamma)} \|\widetilde{\mathbf{J}}_h\|_{H_{\|}^{-1/2}(\Gamma)} \\ &\lesssim \|(Id - \Pi_h^0)\mathbf{J}_h\|_{H_{\|}^{-1/2}(\Gamma)} \|\widetilde{\mathbf{J}}_h\|_{H_{\|}^{-1/2}(\Gamma)}, \end{aligned} \quad (5.48) $$
### 5.3. Discrete inf-sup condition of the V-BEM

In order to prove a h-uniform discrete inf-sup condition on the perturbed sesquilinear form $a_h$ related to the V-BEM (3.28), we rely on both the methodology from Buffa and Hiptmair (2003), (i.e. (5.40)) and the stability result for the ideal V-BEM (see theorem 5.7).

To begin with, we need to derive a technical lemma by using the tools from section 5.1 that provides a h-absorption result on the difference of the sesquilinear form a (3.31) and its perturbation $a_h$ (3.28) with respect to the W-norm.

**Lemma 5.8** $\exists \varepsilon_h \underset{h \to 0^+}{\longrightarrow} 0$ such that, $\forall \mathbf{J}_h \in \mathcal{V}_h$,

$$ \mathcal{E}_h = \sup_{\tilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|a(\mathbf{J}_h, \tilde{\mathbf{J}}_h) - a_h(\mathbf{J}_h, \tilde{\mathbf{J}}_h)|}{\|\tilde{\mathbf{J}}_h\|_W} \le \varepsilon_h \|\mathbf{J}_h\|_W, \quad (5.45) $$

*Proof.* Since the difference of the divergence terms vanishes and by the characterization of W-norm, the following estimate holds:

$$ \mathcal{E}_h = \sup_{\tilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|a(\mathbf{J}_h, \tilde{\mathbf{J}}_h) - a_h(\mathbf{J}_h, \tilde{\mathbf{J}}_h)|}{\|\tilde{\mathbf{J}}_h\|_W} \quad (5.46) $$

$$ = \sup_{\tilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|\langle \mathbf{V}_k \mathbf{J}_h, \tilde{\mathbf{J}}_h \rangle_{\|} - \langle \mathbf{V}_k \Pi_h^0 \mathbf{J}_h, \Pi_h^0 \tilde{\mathbf{J}}_h \rangle_{\|}|}{\|\tilde{\mathbf{J}}_h\|_W}, $$

$$ = \sup_{\tilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|\langle \mathbf{V}_k ((Id - \Pi_h^0) \mathbf{J}_h), \tilde{\mathbf{J}}_h \rangle_{\|} + \langle \mathbf{V}_k \Pi_h^0 \mathbf{J}_h, ((Id - \Pi_h^0) \tilde{\mathbf{J}}_h) \rangle_{\|}|}{\|\tilde{\mathbf{J}}_h\|_W}, $$

$$ \le \sup_{\tilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{\overbrace{\langle \mathbf{V}_k ((Id - \Pi_h^0) \mathbf{J}_h), \tilde{\mathbf{J}}_h \rangle_{\|}}^{A_1} + \overbrace{\langle \mathbf{V}_k \Pi_h^0 \mathbf{J}_h, ((Id - \Pi_h^0) \tilde{\mathbf{J}}_h) \rangle_{\|}}^{A_2}}{\|\tilde{\mathbf{J}}_h\|_{H_{\|}^{-1/2}(\Gamma)}}. \quad (5.47) $$

Since $\mathbf{V}_k$ is continuous from $H_{\|}^{-1/2}(\Gamma)$ to $H_{\|}^{1/2}(\Gamma)$ (corollary 3 of Buffa and Hiptmair (2003)), the terms $A_1$ and $A_2$ of (5.47) can be upper bounded as follows:

$$ \begin{aligned} A_1 & \lesssim \|\mathbf{V}_k((Id - \Pi_h^0)\mathbf{J}_h)\|_{H_{\|}^{1/2}(\Gamma)} \|\tilde{\mathbf{J}}_h\|_{H_{\|}^{-1/2}(\Gamma)} \\ & \lesssim \|(Id - \Pi_h^0)\mathbf{J}_h\|_{H_{\|}^{-1/2}(\Gamma)} \|\tilde{\mathbf{J}}_h\|_{H_{\|}^{-1/2}(\Gamma)}, \end{aligned} \quad (5.48) $$
$$
\begin{align}
A_2 &\lesssim \|\mathbf{V}_k (\Pi_h^0 \mathbf{J}_h)\|_{\mathbf{H}_{\|}^{1/2}(\Gamma)} \| (Id - \Pi_h^0) \widetilde{\mathbf{J}}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \nonumber \\
&\lesssim \|\Pi_h^0 \mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \| (Id - \Pi_h^0) \widetilde{\mathbf{J}}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}. \tag{5.49}
\end{align}
$$

Therefore, (5.47) can be written as

$$
\mathcal{E}_h \lesssim \|(Id - \Pi_h^0) \mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} + \|\Pi_h^0 \mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{\| (Id - \Pi_h^0) \widetilde{\mathbf{J}}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}}{\| \widetilde{\mathbf{J}}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}}. \quad (5.50)
$$

The right-hand side of (5.50) can be estimated by using corollary 5.3 for the first term and lemma 5.2 for the second term: $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that

$$
\mathcal{E}_h \lesssim \beta_h \left( \|\mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} + \|\Pi_h^0 \mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \right). \quad (5.51)
$$

Now, by making use again of corollary 5.3, the second term on the right-hand side of (5.51) can be
bounded from above:

$$
\| \Pi_h^0 \mathbf{J}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \le \| \mathbf{J}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} + \|(Id - \Pi_h^0) \mathbf{J}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}, \quad (5.52)
$$

$$
\le \|\mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} + \beta_h \|\mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}, \quad (5.53)
$$

$$
\le (1 + \beta_h) \| \mathbf{J}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} . \qquad (5.54)
$$

By inserting (5.54) into (5.51) and by definition of the W-norm, we get

$$
\mathcal{E}_h \lesssim \varepsilon_h \|\mathbf{J}_h\|_{\mathbf{W}}, \quad (5.55)
$$

$$
\text{where } \varepsilon_h = \beta_h (2 + \beta_h) \underset{h \to 0^+}{\longrightarrow} 0.
$$

Lemma 5.8 □

With all these technical tools, we are in a position to prove theorem 4.3.

*Proof of theorem 4.3 (h-uniform discrete inf-sup condition of V-BEM).* Firstly, by adding and subtracing the sesquilinear form *a* in the *inf-sup* condition’s argument, we get

$$
\frac{|a_h(\mathbf{J}_h, \mathbf{J}_h') - a(\mathbf{J}_h, \mathbf{J}_h') + a(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \ge \frac{|a(\mathbf{J}_h, \mathbf{J}_h')| - |a(\mathbf{J}_h, \mathbf{J}_h') - a_h(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}}. \quad (5.56)
$$

By considering the test-function $\mathbf{J}'_h$ (5.40) as in the proof of theorem 5.7, we immediately obtain (5.39):
$\exists C > 0$ such that, $\forall \mathbf{J}_h \in \mathcal{V}_h$,

$$
\frac{|a(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \geq C. \quad (5.57)
$$

By inserting (5.57) into (5.56), we get

$$
\frac{|a_h(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \geq C - \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|a(\mathbf{J}_h, \widetilde{\mathbf{J}}_h) - a_h(\mathbf{J}_h, \widetilde{\mathbf{J}}_h)|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\widetilde{\mathbf{J}}_h\|_{\mathbf{W}}}. \quad (5.58)
$$
$$
\begin{align}
A_2 &\lesssim \|\mathbf{V}_k (\Pi_h^0 \mathbf{J}_h)\|_{\mathbf{H}_{\|}^{1/2}(\Gamma)} \| (Id - \Pi_h^0) \widetilde{\mathbf{J}}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \nonumber \\
&\lesssim \|\Pi_h^0 \mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \| (Id - \Pi_h^0) \widetilde{\mathbf{J}}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}. \tag{5.49}
\end{align}
$$

Therefore, (5.47) can be written as

$$
\mathcal{E}_h \lesssim \|(Id - \Pi_h^0) \mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} + \|\Pi_h^0 \mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{\| (Id - \Pi_h^0) \widetilde{\mathbf{J}}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}}{\| \widetilde{\mathbf{J}}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}}. \quad (5.50)
$$

The right-hand side of (5.50) can be estimated by using corollary 5.3 for the first term and lemma 5.2 for the second term: $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$ such that

$$
\mathcal{E}_h \lesssim \beta_h \left( \|\mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} + \|\Pi_h^0 \mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \right). \quad (5.51)
$$

Now, by making use again of corollary 5.3, the second term on the right-hand side of (5.51) can be
bounded from above:

$$
\| \Pi_h^0 \mathbf{J}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} \le \| \mathbf{J}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} + \|(Id - \Pi_h^0) \mathbf{J}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}, \quad (5.52)
$$

$$
\le \|\mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} + \beta_h \|\mathbf{J}_h\|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)}, \quad (5.53)
$$

$$
\le (1 + \beta_h) \| \mathbf{J}_h \|_{\mathbf{H}_{\|}^{-1/2}(\Gamma)} . \quad (5.54)
$$

By inserting (5.54) into (5.51) and by definition of the W-norm, we get

$$
\mathcal{E}_h \lesssim \varepsilon_h \|\mathbf{J}_h\|_{\mathbf{W}}, \quad (5.55)
$$

$$
\text{where } \varepsilon_h = \beta_h (2 + \beta_h) \underset{h \to 0^+}{\longrightarrow} 0.
$$

Lemma 5.8 □

With all these technical tools, we are in a position to prove theorem 4.3.

*Proof of theorem 4.3 (h-uniform discrete inf-sup condition of V-BEM).* Firstly, by adding and subtracing the sesquilinear form *a* in the *inf-sup* condition’s argument, we get

$$
\frac{|a_h(\mathbf{J}_h, \mathbf{J}_h') - a(\mathbf{J}_h, \mathbf{J}_h') + a(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \ge \frac{|a(\mathbf{J}_h, \mathbf{J}_h') - |a(\mathbf{J}_h, \mathbf{J}_h') - a_h(\mathbf{J}_h, \mathbf{J}_h')||}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}}. \quad (5.56)
$$

By considering the test-function $\mathbf{J}'_h$ (5.40) as in the proof of theorem 5.7, we immediately obtain (5.39):
$\exists C > 0$ such that, $\forall \mathbf{J}_h \in \mathcal{V}_h$,

$$
\frac{|a(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \geq C. \qquad (5.57)
$$

By inserting (5.57) into (5.56), we get

$$
\frac{|a_h(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \geq C - \sup_{\widetilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|a(\mathbf{J}_h, \widetilde{\mathbf{J}}_h) - a_h(\mathbf{J}_h, \widetilde{\mathbf{J}}_h)|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\widetilde{\mathbf{J}}_h\|_{\mathbf{W}}}. \quad (5.58)
$$
From lemma 5.8, (5.58) becomes

$$ \frac{|a_h(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \geq (C - C_1 \varepsilon_h), \qquad (5.59) $$

with $C_1 > 0$ and $\varepsilon_h \underset{h \to 0^+}{\longrightarrow} 0$.

To conclude, $\exists h_0 > 0$ such that, $\forall h < h_0$, $(C - C_1 \varepsilon_h) > C/2 > 0$ and $\forall \mathbf{J}_h \in \mathcal{V}_h$,

$$ \sup_{\tilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|a_h(\mathbf{J}_h, \tilde{\mathbf{J}}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\tilde{\mathbf{J}}_h\|_{\mathbf{W}}} \geq \frac{|a_h(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \geq \frac{C}{2} > 0. \qquad (5.60) $$

Theorem. 4.3 □

## 5.4. A priori error estimates

We are here left with the a priori error analysis of the V-BEM scheme. To this end, we need to establish two intermediate results in order to prove theorem 4.5: the continuity and the asymptotic consistency of the perturbed sesquilinear form $a_h$.

*Proof of lemma 4.1 (Continuity).* Firstly, by making use of the triangular and the CS inequalities and the continuity of $V_k$ and $V_k$ (corollary 3 of Buffa and Hiptmair (2003)), we develop $a_h$ as follows

$$ |a_h(\mathbf{J}_h, \mathbf{J}_h')| = \left| \langle V_k \Pi_h^0 \mathbf{J}_h, \Pi_h^0 \mathbf{J}_h' \rangle_{\|} - \frac{1}{\kappa^2} \langle V_k \nabla_{\Gamma} \cdot \mathbf{J}_h, \nabla_{\Gamma} \cdot \mathbf{J}_h' \rangle_{\|} \right|, \qquad (5.61) $$

$$ \leq \left| \langle V_k \Pi_h^0 \mathbf{J}_h, \Pi_h^0 \mathbf{J}_h' \rangle_{\|} \right| + \frac{1}{\kappa^2} \left| \langle V_k \nabla_{\Gamma} \cdot \mathbf{J}_h, \nabla_{\Gamma} \cdot \mathbf{J}_h' \rangle_{\|} \right|, \qquad (5.62) $$

$$ \leq \|V_k \Pi_h^0 \mathbf{J}_h\|_{H_{\|}^{1/2}(\Gamma)} \|\Pi_h^0 \mathbf{J}_h'\|_{H_{\|}^{-1/2}(\Gamma)} \\ + \frac{1}{\kappa^2} \|V_k \nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{1/2}(\Gamma)} \|\nabla_{\Gamma} \cdot \mathbf{J}_h'\|_{H^{-1/2}(\Gamma)}, \qquad (5.63) $$

$$ \lesssim \|\Pi_h^0 \mathbf{J}_h\|_{H_{\|}^{-1/2}(\Gamma)} \|\Pi_h^0 \mathbf{J}_h'\|_{H_{\|}^{-1/2}(\Gamma)} \\ + \frac{1}{\kappa^2} \|\nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{-1/2}(\Gamma)} \|\nabla_{\Gamma} \cdot \mathbf{J}_h'\|_{H^{-1/2}(\Gamma)}. \qquad (5.64) $$

By using the estimate (5.54) from the proof of the inf-sup condition, we get the continuity of $a_h$, with $\beta_h \underset{h \to 0^+}{\longrightarrow} 0$,

$$ \begin{aligned} |a_h(\mathbf{J}_h, \mathbf{J}_h')| & \lesssim (1 + \beta_h)^2 \|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}} + \frac{1}{\kappa^2} \|\nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{-1/2}(\Gamma)} \|\nabla_{\Gamma} \cdot \mathbf{J}_h'\|_{H^{-1/2}(\Gamma)}, \\ & \lesssim \left( (1 + \beta_h)^2 + \frac{1}{\kappa^2} \right) \|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}. \end{aligned} \qquad (5.65) $$

Lemma 4.1 □
From lemma 5.8, (5.58) becomes

$$ \frac{|a_h(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \geq (C - C_1\varepsilon_h), \qquad (5.59) $$

with $C_1 > 0$ and $\varepsilon_h \underset{h \to 0^+}{\longrightarrow} 0$.

To conclude, $\exists h_0 > 0$ such that, $\forall h < h_0$, $(C - C_1\varepsilon_h) > C/2 > 0$ and $\forall \mathbf{J}_h \in \mathcal{V}_h$,

$$ \sup_{\tilde{\mathbf{J}}_h \in \mathcal{V}_h} \frac{|a_h(\mathbf{J}_h, \tilde{\mathbf{J}}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\tilde{\mathbf{J}}_h\|_{\mathbf{W}}} \geq \frac{|a_h(\mathbf{J}_h, \mathbf{J}_h')|}{\|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}} \geq \frac{C}{2} > 0. \qquad (5.60) $$

Theorem. 4.3 □

## 5.4. A priori error estimates

We are here left with the a priori error analysis of the V-BEM scheme. To this end, we need to establish two intermediate results in order to prove theorem 4.5: the continuity and the asymptotic consistency of the perturbed sesquilinear form $a_h$.

*Proof of lemma 4.1 (Continuity).* Firstly, by making use of the triangular and the CS inequalities and the continuity of $V_k$ and $V_k$ (corollary 3 of Buffa and Hiptmair (2003)), we develop $a_h$ as follows

$$ |a_h(\mathbf{J}_h, \mathbf{J}_h')| = \left| \langle V_k \Pi_h^0 \mathbf{J}_h, \Pi_h^0 \mathbf{J}_h' \rangle_{\parallel} - \frac{1}{\kappa^2} \langle V_k \nabla_{\Gamma} \cdot \mathbf{J}_h, \nabla_{\Gamma} \cdot \mathbf{J}_h' \rangle \right|, \qquad (5.61) $$

$$ \leq \left| \langle V_k \Pi_h^0 \mathbf{J}_h, \Pi_h^0 \mathbf{J}_h' \rangle_{\parallel} \right| + \frac{1}{\kappa^2} \left| \langle V_k \nabla_{\Gamma} \cdot \mathbf{J}_h, \nabla_{\Gamma} \cdot \mathbf{J}_h' \rangle \right|, \qquad (5.62) $$

$$ \leq \|V_k \Pi_h^0 \mathbf{J}_h\|_{H_{\parallel}^{1/2}(\Gamma)} \|\Pi_h^0 \mathbf{J}_h'\|_{H_{\parallel}^{-1/2}(\Gamma)} \\ + \frac{1}{\kappa^2} \|V_k \nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{1/2}(\Gamma)} \|\nabla_{\Gamma} \cdot \mathbf{J}_h'\|_{H^{-1/2}(\Gamma)}, \qquad (5.63) $$

$$ \lesssim \|\Pi_h^0 \mathbf{J}_h\|_{H_{\parallel}^{-1/2}(\Gamma)} \|\Pi_h^0 \mathbf{J}_h'\|_{H_{\parallel}^{-1/2}(\Gamma)} \\ + \frac{1}{\kappa^2} \|\nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{-1/2}(\Gamma)} \|\nabla_{\Gamma} \cdot \mathbf{J}_h'\|_{H^{-1/2}(\Gamma)}. \qquad (5.64) $$

By using the estimate (5.54) from the proof of the inf-sup condition, we get the continuity of $a_h$, with $\beta_h \underset{h \to 0^+}{\longrightarrow} 0$,

$$ \begin{aligned} |a_h(\mathbf{J}_h, \mathbf{J}_h')| & \lesssim (1 + \beta_h)^2 \|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}} + \frac{1}{\kappa^2} \|\nabla_{\Gamma} \cdot \mathbf{J}_h\|_{H^{-1/2}(\Gamma)} \|\nabla_{\Gamma} \cdot \mathbf{J}_h'\|_{H^{-1/2}(\Gamma)}, \\ & \lesssim \left( (1 + \beta_h)^2 + \frac{1}{\kappa^2} \right) \|\mathbf{J}_h\|_{\mathbf{W}} \|\mathbf{J}_h'\|_{\mathbf{W}}. \end{aligned} \qquad (5.65) $$

Lemma 4.1 □
*Proof of proposition 4.2 (Asymptotic consistency).* Firstly, by adding and subtracting the sesquilinear form $a$, as well as the right-hand side $f$ in the sup's argument in (4.2), we have that

$$
\begin{align}
R_h(\mathbf{J}) &\le \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{|f_h(\mathbf{J}'_h) - f(\mathbf{J}'_h)| + \underbrace{|f(\mathbf{J}'_h) - a(\mathbf{J}, \mathbf{J}'_h)|}_{=0} + |a(\mathbf{J}, \mathbf{J}'_h) - a_h(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|}{\|\mathbf{J}'_h\|_W}, \nonumber \\
&\le \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \left( \frac{|f_h(\mathbf{J}'_h) - f(\mathbf{J}'_h)| + |a(\mathbf{J}, \mathbf{J}'_h) - a(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|}{\|\mathbf{J}'_h\|_W} \right. \nonumber \\
&\qquad \left. + \frac{|a(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h) - a_h(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|}{\|\mathbf{J}'_h\|_W} \right) \tag{5.66}
\end{align}
$$

By using the following property associated to the $L^2$-projector:

$$ \forall \mathbf{f}, \mathbf{g} \in L_t^2(\Gamma), \int_{\Gamma} \mathbf{f} \cdot \Pi_h^0 \mathbf{g} \, d\gamma = \int_{\Gamma} \Pi_h^0 \mathbf{f} \cdot \mathbf{g} \, d\gamma, \quad (5.67) $$

and the continuity of the sesquilinear form $a$, the first two terms in (5.66) can be upper bounded as follows

$$ |f_h(\mathbf{J}'_h) - f(\mathbf{J}'_h)| \le \|(Id - \Pi_h^0) \mathbf{f}\|_0 \|\mathbf{J}'_h\|_0, \quad (5.68) $$

$$ |a(\mathbf{J}, \mathbf{J}'_h) - a(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)| = |a(\mathbf{J} - \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|, \quad (5.69) $$

$$ \lesssim \|\mathbf{J} - \mathcal{P}_h \mathbf{J}\|_W \|\mathbf{J}'_h\|_W. \quad (5.70) $$

We remark that in (5.68), $f$ identifies the Riesz representation of the right-hand side $f$. Moreover, since the difference of the divergence terms vanishes in the last term of (5.66), we can write, by using again (5.67), the following estimate

$$
\begin{align}
R_h(\mathbf{J}) &\le \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\| (Id - \Pi_h^0) \mathbf{f} \|_0 \| \mathbf{J}'_h \|_0 + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_W \| \mathbf{J}'_h \|_W}{\| \mathbf{J}'_h \|_W} \nonumber \\
&\quad + \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\left| \langle \mathbf{V}_k \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h \rangle_{\|} - \langle \mathbf{V}_k \Pi_h^0 \mathcal{P}_h \mathbf{J}, \Pi_h^0 \mathbf{J}'_h \rangle_{\|} \right|}{\| \mathbf{J}'_h \|_W}, \tag{5.71}
\end{align}
$$

$$
\begin{align}
&\le \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\| (Id - \Pi_h^0) \mathbf{f} \|_0 \| \mathbf{J}'_h \|_0 + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_W \| \mathbf{J}'_h \|_W}{\| \mathbf{J}'_h \|_W} \nonumber \\
&\quad + \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\left| \langle \mathbf{V}_k \mathcal{P}_h \mathbf{J} - \Pi_h^0 \mathbf{V}_k \Pi_h^0 \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h \rangle_{\|} \right|}{\| \mathbf{J}'_h \|_W}. \tag{5.72}
\end{align}
$$
*Proof of proposition 4.2 (Asymptotic consistency).* Firstly, by adding and subtracting the sesquilinear form $a$, as well as the right-hand side $f$ in the sup's argument in (4.2), we have that

$$
\begin{align}
R_h(\mathbf{J}) &\le \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{|f_h(\mathbf{J}'_h) - f(\mathbf{J}'_h)| + \underbrace{|f(\mathbf{J}'_h) - a(\mathbf{J}, \mathbf{J}'_h)|}_{=0} + |a(\mathbf{J}, \mathbf{J}'_h) - a_h(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|}{\|\mathbf{J}'_h\|_W}, \nonumber \\
&\le \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \left( \frac{|f_h(\mathbf{J}'_h) - f(\mathbf{J}'_h)| + |a(\mathbf{J}, \mathbf{J}'_h) - a(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|}{\|\mathbf{J}'_h\|_W} \right. \nonumber \\
&\qquad \left. + \frac{|a(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h) - a_h(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|}{\|\mathbf{J}'_h\|_W} \right) \tag{5.66}
\end{align}
$$

By using the following property associated to the $L^2$-projector:

$$ \forall \mathbf{f}, \mathbf{g} \in \mathcal{L}_t^2(\Gamma), \int_{\Gamma} \mathbf{f} \cdot \Pi_h^0 \mathbf{g} \, d\gamma = \int_{\Gamma} \Pi_h^0 \mathbf{f} \cdot \mathbf{g} \, d\gamma, \quad (5.67) $$

and the continuity of the sesquilinear form $a$, the first two terms in (5.66) can be upper bounded as follows

$$ |f_h(\mathbf{J}'_h) - f(\mathbf{J}'_h)| \le \|(Id - \Pi_h^0) \mathbf{f}\|_0 \|\mathbf{J}'_h\|_0, \quad (5.68) $$

$$ |a(\mathbf{J}, \mathbf{J}'_h) - a(\mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)| = |a(\mathbf{J} - \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h)|, \quad (5.69) $$

$$ \lesssim \|\mathbf{J} - \mathcal{P}_h \mathbf{J}\|_{\mathbf{W}} \|\mathbf{J}'_h\|_{\mathbf{W}}. \quad (5.70) $$

We remark that in (5.68), $f$ identifies the Riesz representation of the right-hand side $f$. Moreover, since the difference of the divergence terms vanishes in the last term of (5.66), we can write, by using again (5.67), the following estimate

$$
\begin{align}
R_h(\mathbf{J}) &\le \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\| (Id - \Pi_h^0) \mathbf{f} \|_0 \| \mathbf{J}'_h \|_0 + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}} \| \mathbf{J}'_h \|_{\mathbf{W}}}{\| \mathbf{J}'_h \|_{\mathbf{W}}} \nonumber \\
&\quad + \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\left| \langle \mathbf{V}_k \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h \rangle_{\|} - \langle \mathbf{V}_k \Pi_h^0 \mathcal{P}_h \mathbf{J}, \Pi_h^0 \mathbf{J}'_h \rangle_{\|} \right|}{\| \mathbf{J}'_h \|_{\mathbf{W}}}, \tag{5.71}
\end{align}
$$

$$
\begin{align}
&\le \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\| (Id - \Pi_h^0) \mathbf{f} \|_0 \| \mathbf{J}'_h \|_0 + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}} \| \mathbf{J}'_h \|_{\mathbf{W}}}{\| \mathbf{J}'_h \|_{\mathbf{W}}} \nonumber \\
&\quad + \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\left| \langle \mathbf{V}_k \mathcal{P}_h \mathbf{J} - \Pi_h^0 \mathbf{V}_k \Pi_h^0 \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h \rangle_{\|} \right|}{\| \mathbf{J}'_h \|_{\mathbf{W}}}. \tag{5.72}
\end{align}
$$
By using the inverse estimate 5.1 between $H_{\parallel}^{-1/2}(\Gamma)$ and $L_t^2(\Gamma)$, (5.72) becomes:

$$ R_h(\mathbf{J}) \lesssim h^{-1/2} \|(Id - \Pi_h^0) \mathbf{f}\|_0 + \|\mathbf{J} - \mathcal{P}_h \mathbf{J}\|_{\mathbf{W}} \\ + \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\underbrace{\langle \mathbf{V}_k \mathcal{P}_h \mathbf{J} - \Pi_h^0 \mathbf{V}_k \Pi_h^0 \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h \rangle_{\parallel}}_{A_1}}{\|\mathbf{J}'_h\|_{\mathbf{W}}}. \quad (5.73) $$

By rewriting $A_1$ in terms of some $L^2$-projection errors, we get:

$$ \begin{aligned} A_1 &= |\langle (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} + (Id - \Pi_h^0) \mathbf{V}_k (\mathcal{P}_h \mathbf{J} - \mathbf{J}) + \Pi_h^0 \mathbf{V}_k (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h \rangle_{\parallel}| \\ &\leq \left( \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k (\mathcal{P}_h \mathbf{J} - \mathbf{J}) \|_0 \right. \\ &\qquad \left. + \| \Pi_h^0 \mathbf{V}_k (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_0 \right) \| \mathbf{J}'_h \|_0. \end{aligned} \quad (5.74) $$

So, from the inverse estimate 5.1, we obtain a new upper-bound of $R_h$:

$$ R_h(\mathbf{J}) \lesssim h^{-1/2} \left( \| (Id - \Pi_h^0) \mathbf{f} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 \right. \\ \qquad \left. + \| (Id - \Pi_h^0) \mathbf{V}_k (\mathcal{P}_h \mathbf{J} - \mathbf{J}) \|_0 + \| \Pi_h^0 \mathbf{V}_k (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_0 \right) \\ \qquad + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}}. \quad (5.75) $$

Now, by using the projection error from (5.12), the continuity properties of $V_k$, the $L^2$-stability of $\Pi_h^0$ i.e $\|\Pi_h^0 \mathbf{u}\|_0 \leq \|\mathbf{u}\|_0$ and the projection error estimate

$$ \| (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_{H_{\parallel}^{-1}(\Gamma)} \lesssim h^{1/2} \| (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_{H_{\parallel}^{-1/2}(\Gamma)}, $$

obtained by the standard duality argument (5.72) leads to:

$$ R_h(\mathbf{J}) \lesssim h^{-1/2} \left( \| (Id - \Pi_h^0) \mathbf{f} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 + h^{1/2} \| \mathbf{V}_k (\mathcal{P}_h \mathbf{J} - \mathbf{J}) \|_{H_{\parallel}^{1/2}(\Gamma)} \right. \\ \qquad \left. + \| \mathbf{V}_k (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_0 \right) + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}}, \quad (5.76) $$

$$ \lesssim h^{-1/2} \left( \| (Id - \Pi_h^0) \mathbf{f} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 + \| (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_{H_{\parallel}^{-1}(\Gamma)} \right) \\ + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}} + \|(\mathcal{P}_h \mathbf{J} - \mathbf{J})\|_{H_{\parallel}^{-1/2}(\Gamma)}, \quad (5.77) $$

$$ \lesssim h^{-1/2} \left( \| (Id - \Pi_h^0) \mathbf{f} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 + h^{1/2} \| (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_{H_{\parallel}^{-1/2}(\Gamma)} \right) \\ + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}}. \quad (5.78) $$
By using the inverse estimate 5.1 between $H_{\parallel}^{-1/2}(\Gamma)$ and $L_t^2(\Gamma)$, (5.72) becomes:

$$ R_h(\mathbf{J}) \lesssim h^{-1/2} \|(Id - \Pi_h^0) \mathbf{f}\|_0 + \|\mathbf{J} - \mathcal{P}_h \mathbf{J}\|_{\mathbf{W}} \\ + \sup_{\mathbf{J}'_h \in \mathcal{V}_h} \frac{\underbrace{\langle \mathbf{V}_k \mathcal{P}_h \mathbf{J} - \Pi_h^0 \mathbf{V}_k \Pi_h^0 \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h \rangle_{\parallel}}_{A_1}}{\|\mathbf{J}'_h\|_{\mathbf{W}}}. \quad (5.73) $$

By rewriting $A_1$ in terms of some $L^2$-projection errors, we get:

$$ A_1 = \left| \begin{aligned} & \langle (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} + (Id - \Pi_h^0) \mathbf{V}_k (\mathcal{P}_h \mathbf{J} - \mathbf{J}) + \Pi_h^0 \mathbf{V}_k (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J}, \mathbf{J}'_h \rangle_{\parallel} \\ & \le \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k (\mathcal{P}_h \mathbf{J} - \mathbf{J}) \|_0 \\ & \quad + \| \Pi_h^0 \mathbf{V}_k (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_0 \end{aligned} \right| \|\mathbf{J}'_h\|_0. \quad (5.74) $$

So, from the inverse estimate 5.1, we obtain a new upper-bound of $R_h$:

$$ R_h(\mathbf{J}) \lesssim h^{-1/2} \left( \begin{aligned} & \| (Id - \Pi_h^0) \mathbf{f} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 \\ & + \| (Id - \Pi_h^0) \mathbf{V}_k (\mathcal{P}_h \mathbf{J} - \mathbf{J}) \|_0 + \| \Pi_h^0 \mathbf{V}_k (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_0 \\ & + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}} \end{aligned} \right). \quad (5.75) $$

Now, by using the projection error from (5.12), the continuity properties of $V_k$, the $L^2$-stability of $\Pi_h^0$ i.e $\|\Pi_h^0 \mathbf{u}\|_0 \le \|\mathbf{u}\|_0$ and the projection error estimate

$$ \| (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_{H_{\parallel}^{-1}(\Gamma)} \lesssim h^{1/2} \| (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_{H_{\parallel}^{-1/2}(\Gamma)}, $$

obtained by the standard duality argument (5.72) leads to:

$$ R_h(\mathbf{J}) \lesssim h^{-1/2} \left( \begin{aligned} & \| (Id - \Pi_h^0) \mathbf{f} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 + h^{1/2} \| \mathbf{V}_k (\mathcal{P}_h \mathbf{J} - \mathbf{J}) \|_{H_{\parallel}^{1/2}(\Gamma)} \\ & + \| \mathbf{V}_k (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_0 \end{aligned} \right) + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}}, \quad (5.76) $$

$$ \lesssim h^{-1/2} \left( \begin{aligned} & \| (Id - \Pi_h^0) \mathbf{f} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 + \| (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_{H_{\parallel}^{-1}(\Gamma)} \\ & + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}} + \| (\mathcal{P}_h \mathbf{J} - \mathbf{J}) \|_{H_{\parallel}^{-1/2}(\Gamma)} \end{aligned} \right), \quad (5.77) $$

$$ \lesssim h^{-1/2} \left( \begin{aligned} & \| (Id - \Pi_h^0) \mathbf{f} \|_0 + \| (Id - \Pi_h^0) \mathbf{V}_k \mathbf{J} \|_0 + h^{1/2} \| (Id - \Pi_h^0) \mathcal{P}_h \mathbf{J} \|_{H_{\parallel}^{-1/2}(\Gamma)} \\ & + \| \mathbf{J} - \mathcal{P}_h \mathbf{J} \|_{\mathbf{W}} \end{aligned} \right). \quad (5.78) $$
Finally, from corollary 5.3 and the W-stability of $P_h$, $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$, such that

$$R_h(\mathbf{J}) \lesssim h^{-1/2} \left( \|(Id - \Pi_h^0) \mathbf{f}\|_0 + \|(Id - \Pi_h^0) \mathbf{V}_k \mathbf{J}\|_0 \right) + \beta_h \|\mathbf{J}\|_{\mathbf{W}} + \|\mathbf{J} - P_h \mathbf{J}\|_{\mathbf{W}}, \quad (5.79)$$

and we get the result by using the extra-regularity assumption, i.e. $\mathbf{J} \in \mathbf{W} \cap \mathbf{H}_{\|}^{-1/2+\epsilon}(\Gamma)$, and the regularity of the right-hand side $\mathbf{f}$ that is the (smooth) tangential component of a plane wave.

Prop. 4.2 □

Now, we are able to give some results about the h-convergence of the V-BEM.

**Proof of theorem 4.5 (A priori error estimate).** The estimate (4.4) is obtained from theorem 2.24 of Ern and Guermond (2004) by using the continuity (lemma 4.1) and the asymptotic consistency (proposition 4.2) results as well as the uniform discrete inf-sup condition (theorem 4.3). The convergence of the scheme is ensured by the approximability property of $V_h$, which is induced by the local RT polynomial consistency of the virtual element approximation.

Theorem 4.5 □

As a consequence, we can establish the convergence rate of the V-BEM from some extra-regularity assumptions. To that end, we need two intermediate lemmas for proving corollary 4.6.

**Lemma 5.9 (Convergence rate of consistency error)** Let $\mathbf{J} \in \mathbf{H}^{\sigma}(\mathrm{div}_{\Gamma}, \Gamma)$ for $\sigma \in [0, 1/2)$ such that $\mathbf{V}_k \mathbf{J} \in \mathbf{H}_{\|}^{\sigma+1}(\Gamma)$. The following estimate of the consistency error holds:

$$R_h(\mathbf{J}) \lesssim h^{1/2+\sigma} \left( \|\mathbf{f}\|_{\mathbf{H}_{\|}^2(\Gamma)} + \|\mathbf{V}_k \mathbf{J}\|_{\mathbf{H}_{\|}^{\sigma+1}(\Gamma)} + \|\mathbf{J}\|_{\mathbf{H}_{\|}^{\sigma}(\Gamma)} \right). \quad (5.80)$$

*Proof.* From the regularity assumptions, $\mathbf{J} \in L_t^2(\Gamma)$ and the consistency error can be defined by taking $P_h = \mathrm{Id}$. Consequently, (5.76) gives the simplified estimate:

$$R_h(\mathbf{J}) \lesssim h^{-1/2} \left( \|(Id - \Pi_h^0) \mathbf{f}\|_0 + \|(Id - \Pi_h^0) \mathbf{V}_k \mathbf{J}\|_0 + \|\mathbf{V}_k (Id - \Pi_h^0) \mathbf{J}\|_0 \right). \quad (5.81)$$

We now estimate each term of the right-hand side of (5.81):

* $\mathbf{f} \in \mathbf{H}_{\|}^2(\Gamma)$ because it is a tangential component of a plane wave. So, the error estimate associated to the $L^2$-projection $\Pi_h^0$ immediately gives

$$\|(Id - \Pi_h^0) \mathbf{f}\|_0 \le C h^2 \|\mathbf{f}\|_{\mathbf{H}_{\|}^2(\Gamma)} \Rightarrow h^{-1/2} \|(Id - \Pi_h^0) \mathbf{f}\|_0 \lesssim h^{3/2} \|\mathbf{f}\|_{\mathbf{H}_{\|}^2(\Gamma)}. \quad (5.82)$$

* Again, from the properties of $\Pi_h^0$ and the assumption $\mathbf{V}_k \mathbf{J} \in \mathbf{H}_{\|}^{\sigma+1}(\Gamma)$, we get

$$\|(Id - \Pi_h^0) \mathbf{V}_k \mathbf{J}\|_0 \le C h^{\sigma+1} \|\mathbf{V}_k \mathbf{J}\|_{\mathbf{H}_{\|}^{\sigma+1}(\Gamma)}, \quad (5.83)$$

$$\Rightarrow h^{-1/2} \|(Id - \Pi_h^0) \mathbf{V}_k \mathbf{J}\|_0 \lesssim h^{\sigma+1/2} \|\mathbf{V}_k \mathbf{J}\|_{\mathbf{H}_{\|}^{\sigma+1}(\Gamma)}. \quad (5.84)$$

* Finally, from the continuity of $\mathbf{V}_k$ from $\mathbf{H}_{\|}^{-1}(\Gamma)$ to $L_t^2(\Gamma)$ (Buffa et al. (2003, th. 3.8), Heuer and Meddahi (2013, (4.1))) and a standard duality argument for $\Pi_h^0$, we obtain the last estimate:
Finally, from corollary 5.3 and the W-stability of $P_h$, $\exists \beta_h \underset{h \to 0^+}{\longrightarrow} 0$, such that

$$R_h(J) \lesssim h^{-1/2} \left( \|(Id - \Pi_h^0) f\|_0 + \|(Id - \Pi_h^0) V_k J\|_0 \right) + \beta_h \|J\|_W + \|J - P_h J\|_W, \quad (5.79)$$

and we get the result by using the extra-regularity assumption, i.e. $J \in W \cap H_{\parallel}^{-1/2+\epsilon}(\Gamma)$, and the regularity of the right-hand side $f$ that is the (smooth) tangential component of a plane wave.

Prop. 4.2 □

Now, we are able to give some results about the h-convergence of the V-BEM.

**Proof of theorem 4.5 (A priori error estimate).** The estimate (4.4) is obtained from theorem 2.24 of Ern and Guermond (2004) by using the continuity (lemma 4.1) and the asymptotic consistency (proposition 4.2) results as well as the uniform discrete inf-sup condition (theorem 4.3). The convergence of the scheme is ensured by the approximability property of $V_h$, which is induced by the local RT polynomial consistency of the virtual element approximation.

Theorem 4.5 □

As a consequence, we can establish the convergence rate of the V-BEM from some extra-regularity assumptions. To that end, we need two intermediate lemmas for proving corollary 4.6.

**Lemma 5.9 (Convergence rate of consistency error)** Let $J \in H^\sigma(\mathrm{div}_\Gamma, \Gamma)$ for $\sigma \in [0, 1/2)$ such that $V_k J \in H_{\parallel}^{\sigma+1}(\Gamma)$. The following estimate of the consistency error holds:

$$R_h(J) \lesssim h^{1/2+\sigma} \left( \|f\|_{H^2(\Gamma)} + \|V_k J\|_{H_{\parallel}^{\sigma+1}(\Gamma)} + \|J\|_{H_{\parallel}^{\sigma}(\Gamma)} \right). \quad (5.80)$$

*Proof.* From the regularity assumptions, $J \in L_t^2(\Gamma)$ and the consistency error can be defined by taking $P_h = Id$. Consequently, (5.76) gives the simplified estimate:

$$R_h(J) \lesssim h^{-1/2} \left( \|(Id - \Pi_h^0) f\|_0 + \|(Id - \Pi_h^0) V_k J\|_0 + \|V_k (Id - \Pi_h^0) J\|_0 \right). \quad (5.81)$$

We now estimate each term of the right-hand side of (5.81):

* $f \in H_{\parallel}^2(\Gamma)$ because it is a tangential component of a plane wave. So, the error estimate associated to the $L^2$-projection $\Pi_h^0$ immediately gives

$$\|(Id - \Pi_h^0) f\|_0 \le C h^2 \|f\|_{H^2(\Gamma)} \Rightarrow h^{-1/2} \|(Id - \Pi_h^0) f\|_0 \lesssim h^{3/2} \|f\|_{H^2(\Gamma)}. \quad (5.82)$$

* Again, from the properties of $\Pi_h^0$ and the assumption $V_k J \in H_{\parallel}^{\sigma+1}(\Gamma)$, we get

$$\|(Id - \Pi_h^0) V_k J\|_0 \le C h^{\sigma+1} \|V_k J\|_{H_{\parallel}^{\sigma+1}(\Gamma)}, \quad (5.83)$$

$$\Rightarrow h^{-1/2} \|(Id - \Pi_h^0) V_k J\|_0 \lesssim h^{\sigma+1/2} \|V_k J\|_{H_{\parallel}^{\sigma+1}(\Gamma)}. \quad (5.84)$$

* Finally, from the continuity of $V_k$ from $H_{\parallel}^{-1}(\Gamma)$ to $L_t^2(\Gamma)$ (Buffa et al. (2003, th. 3.8), Heuer and Meddahi (2013, (4.1))) and a standard duality argument for $\Pi_h^0$, we obtain the last estimate:
$$ \|V_k (Id - \Pi_h^0) J\|_0 \lesssim \|(Id - \Pi_h^0) J\|_{H_{\parallel}^{-1}(\Gamma)}, \quad (5.85) $$

$$ \lesssim h \|(Id - \Pi_h^0) J\|_0, $$

$$ \lesssim h^{1+\sigma} \|J\|_{H_{-\sigma}^{\sigma}(\Gamma)}, $$

$$ \Rightarrow h^{-1/2} \|V_k (Id - \Pi_h^0) J\|_0 \lesssim h^{1/2+\sigma} \|J\|_{H^{\sigma}(\Gamma)}. \quad (5.86) $$

Hence, by combining (5.82), (5.84) and (5.86), we obtain the result.

Lemma 5.9 □

**Lemma 5.10** Let $J \in H^{\sigma}(\mathrm{div}_{\Gamma}, \Gamma)$ for $\sigma \in (0, 1/2)$. The optimal approximation error can be estimated by

$$ \inf_{\xi_h \in \mathcal{V}_h} \|J - \xi_h\|_{W} \lesssim h^{1/2+\sigma} \|J\|_{H^{\sigma}(\mathrm{div}_{\Gamma}, \Gamma)}. \quad (5.87) $$

*Proof.* The proof is straightforward by using Proposition 4.6 in Buffa and Christiansen (2003). In fact, the VEM interpolation operator $P_h^X$ satisfies the same properties as its lowest-order RT finite element counterparts:

* the virtual degrees of freedom are the edge fluxes,

* the a priori error estimate : $\forall u \in W$ such that $\forall j \in [1, N_{\Gamma}]$, $u|_{\Gamma_j} \in H^s(\mathrm{div}_{\Gamma}, \Gamma_j)$ with $s > 0$,

$$ \|u - P_h^X u\|_{H^0(\mathrm{div}_{\Gamma}, \Gamma)} \lesssim h^s \sum_{j=1}^{N_{\Gamma}} \|u|_{\Gamma_j}\|_{H^s(\mathrm{div}_{\Gamma}, \Gamma_j)}. \quad (5.88) $$

The estimate (5.88) is directly obtained from both the local VEM error estimates (theorem A.1, proposition 3.2 and corollary 3.3 of Beirão da Veiga and Mascotto (2022)) and the results of Ern and Guermond (2021) (chapters 16 and 17).

Lemma 5.10 □

*Proof of corollary 4.6 (Convergence rate).* The corollary holds by using lemmas 5.9 and 5.10 on the error estimate of theorem 4.5.

Corollary 4.6 □

## 6. Numerical results

In this section, we present numerical experiments in order to analyze the behavior in terms of h-convergence and accuracy of the proposed virtual boundary element scheme through the comparison with the BEM approach based on the conventional lowest-order triangular RT elements. The symmetric indefinite non-Hermitian linear systems originating from the EFIE discretization on the basis of both V-BEM and BEM are solved using a direct method, with fully distribute memory (MPI-type) parallelism, which relies on LU factorization algorithms from LAPACK library. All numerical simulations were run on the CEA's high-performance computing architectures.

For this study, two representative electromagnetic scattering test cases involving solutions with different regularities are considered:

* A smooth geometry corresponding to a metallic sphere, of radius 1 m, illuminated by a linearly-polarized incident plane wave of frequency 500 MHz whose wave vector is in the opposite direction of z-axis (i.e. (0°, 0°) in spherical coordinates, see figure 2). In this case, the analytical solution has $C^{\infty}$-regularity and is available through Mie series (Monk, 2003; Nedelec, 2001).
$$ \|V_k (Id - \Pi_h^0) J\|_0 \lesssim \|(Id - \Pi_h^0) J\|_{H_{\parallel}^{-1}(\Gamma)}, \quad (5.85) $$

$$ \lesssim h \|(Id - \Pi_h^0) J\|_0, $$

$$ \lesssim h^{1+\sigma} \|J\|_{H_{-\parallel}^{\sigma}(\Gamma)}, $$

$$ \Rightarrow h^{-1/2} \|V_k (Id - \Pi_h^0) J\|_0 \lesssim h^{1/2+\sigma} \|J\|_{H^{\sigma}(\Gamma)}. \quad (5.86) $$

Hence, by combining (5.82), (5.84) and (5.86), we obtain the result.

Lemma 5.9 □

**Lemma 5.10** Let $J \in H^{\sigma}(\mathrm{div}_{\Gamma}, \Gamma)$ for $\sigma \in (0, 1/2)$. The optimal approximation error can be estimated by

$$ \inf_{\xi_h \in \mathcal{V}_h} \|J - \xi_h\|_{W} \lesssim h^{1/2+\sigma} \|J\|_{H^{\sigma}(\mathrm{div}_{\Gamma}, \Gamma)}. \quad (5.87) $$

*Proof.* The proof is straightforward by using Proposition 4.6 in Buffa and Christiansen (2003). In fact, the VEM interpolation operator $P_h^X$ satisfies the same properties as its lowest-order RT finite element counterparts:

* the virtual degrees of freedom are the edge fluxes,

* the a priori error estimate : $\forall u \in W$ such that $\forall j \in [1, N_{\Gamma}]$, $u|_{\Gamma_j} \in H^s(\mathrm{div}_{\Gamma}, \Gamma_j)$ with $s > 0$,

$$ \|u - P_h^X u\|_{H^0(\mathrm{div}_{\Gamma}, \Gamma)} \lesssim h^s \sum_{j=1}^{N_{\Gamma}} \|u|_{\Gamma_j}\|_{H^s(\mathrm{div}_{\Gamma}, \Gamma_j)}. \quad (5.88) $$

The estimate (5.88) is directly obtained from both the local VEM error estimates (theorem A.1, proposition 3.2 and corollary 3.3 of Beirão da Veiga and Mascotto (2022)) and the results of Ern and Guermond (2021) (chapters 16 and 17).

Lemma 5.10 □

*Proof of corollary 4.6 (Convergence rate).* The corollary holds by using lemmas 5.9 and 5.10 on the error estimate of theorem 4.5.

Corollary 4.6 □

## 6. Numerical results

In this section, we present numerical experiments in order to analyze the behavior in terms of $h$-convergence and accuracy of the proposed virtual boundary element scheme through the comparison with the BEM approach based on the conventional lowest-order triangular RT elements. The symmetric indefinite non-Hermitian linear systems originating from the EFIE discretization on the basis of both V-BEM and BEM are solved using a direct method, with fully distribute memory (MPI-type) parallelism, which relies on LU factorization algorithms from LAPACK library. All numerical simulations were run on the CEA's high-performance computing architectures.

For this study, two representative electromagnetic scattering test cases involving solutions with different regularities are considered:

* A smooth geometry corresponding to a metallic sphere, of radius 1 m, illuminated by a linearly-polarized incident plane wave of frequency 500 MHz whose wave vector is in the opposite direction of z-axis (i.e. (0°, 0°) in spherical coordinates, see figure 2). In this case, the analytical solution has $C^{\infty}$-regularity and is available through Mie series (Monk, 2003; Nedelec, 2001).
* A nonsmooth geometry corresponding to a metallic cube, of 1 m of side, illuminated by a linearly-polarized incident plane wave of frequency 800 MHz whose wave vector has (45°, 45°) direction in spherical coordinates (or rather, points along the body diagonal from the cube top, see figure 2). In this case, the geometrical singularities (edges and corners) induce low-regularity electric currents and, as the analytical solution is not explicitly known, we replace it by a reference solution (referred to as J, by abusing of notation) computed on a sufficiently fine simplicial mesh.

The quantities of interest for the study are the following:

* The relative error on J in the $H^0(\text{div}_\Gamma, \Gamma)$-type norm defined as

$$e_r(J_h) := \frac{\sqrt{\|J - \Pi_h^0 J_h\|_0^2 + \frac{1}{\kappa^2} \|\nabla_\Gamma \cdot J - \nabla_\Gamma \cdot J_h\|_0^2}}{\sqrt{\|J\|_0^2 + \frac{1}{\kappa^2} \|\nabla_\Gamma \cdot J\|_0^2}}, \quad (6.1)$$

where $J_h$ is replaced by $\Pi_h^0 J_h$ in the first part of the norm due to the implicit definition of virtual functions. Note that we chose to evaluate the solution error in the above stronger norm with respect to the $H^{-1/2}(\text{div}_\Gamma, \Gamma)$-norm, as the latter is difficult to compute mainly due to its non-local character.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Example of three types of meshes of the sphere and the cube surface with an equivalent $h_{mean}$ (see (6.6)): the maximum number of edges within a polygon for the sphere and cube meshes is 6 (ad-hoc) and 11 (random), respectively.</figcaption></figure>
* A nonsmooth geometry corresponding to a metallic cube, of 1 m of side, illuminated by a linearly-polarized incident plane wave of frequency 800 MHz whose wave vector has (45°, 45°) direction in spherical coordinates (or rather, points along the body diagonal from the cube top, see figure 2). In this case, the geometrical singularities (edges and corners) induce low-regularity electric currents and, as the analytical solution is not explicitly known, we replace it by a reference solution (referred to as J, by abusing of notation) computed on a sufficiently fine simplicial mesh.

The quantities of interest for the study are the following:

* The relative error on J in the $H^0(\text{div}_\Gamma, \Gamma)$-type norm defined as

$$e_r(J_h) := \frac{\sqrt{\|J - \Pi_h^0 J_h\|_0^2 + \frac{1}{\kappa^2} \|\nabla_\Gamma \cdot J - \nabla_\Gamma \cdot J_h\|_0^2}}{\sqrt{\|J\|_0^2 + \frac{1}{\kappa^2} \|\nabla_\Gamma \cdot J\|_0^2}}, \quad (6.1)$$

where $J_h$ is replaced by $\Pi_h^0 J_h$ in the first part of the norm due to the implicit definition of virtual functions. Note that we chose to evaluate the solution error in the above stronger norm with respect to the $H^{-1/2}(\text{div}_\Gamma, \Gamma)$-norm, as the latter is difficult to compute mainly due to its non-local character.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Example of three types of meshes of the sphere and the cube surface with an equivalent $h_{mean}$ (see (6.6)): the maximum number of edges within a polygon for the sphere and cube meshes is 6 (ad-hoc) and 11 (random), respectively.</figcaption></figure>
* The monostatic and bistatic RCS errors are defined as, respectively,

$$E_{RCS}^{m} = \left| \sigma(-\hat{\xi}_{I}, \mathbf{J}) - \sigma(-\hat{\xi}_{I}, \Pi_{h}^{0} \mathbf{J}_{h}) \right|, \qquad (6.2)$$

$$E_{RCS}^{b}(\hat{\xi}_{s}) = \left| \sigma(\hat{\xi}_{s}, \mathbf{J}) - \sigma(\hat{\xi}_{s}, \Pi_{h}^{0} \mathbf{J}_{h}) \right|^{2}, \qquad (6.3)$$

where $\hat{\xi}_{I}$ and $\hat{\xi}_{s}$ are the unit vectors associated to the directions of the incident wave and the scattered wave, respectively, and $\sigma$ is related to the scattering amplitude (Nedelec, 2001) and is defined as

$$\sigma(\hat{\xi}_{s}, \mathbf{v}) = \frac{-i\kappa}{\sqrt{4\pi}} \int_{\Gamma_{\text{or}} \Gamma_{h}} e^{i\kappa(x \cdot \hat{\xi}_{s})} (\mathbf{v}(x) - \hat{\xi}_{s}(\mathbf{v}(x) \cdot \hat{\xi}_{s})) \cdot \hat{e}_{r} d\gamma_{x}, \qquad (6.4)$$

with $\hat{e}_{r}$ being the direction of the scattered wave's polarization.

* The bistatic RCS, being a relevant observable in radar applications, characterizes the ability of a body to send back the incident electromagnetic wave energy into different directions and is defined, with $\mathbf{v}$ being $\mathbf{J}$ or $\Pi_{h}^{0} \mathbf{J}_{h}$, as

$$\text{RCS}(\hat{\xi}_{s}) = \left| \sigma(\hat{\xi}_{s}, \mathbf{v}) \right|^{2}. \qquad (6.5)$$

**Remark 6.1** We highlight the fact that $\Pi_{h}^{0} \mathbf{J}_{h} = \mathbf{J}_{h}$ on triangular RT elements, hence the above quantity definitions are modified accordingly.

**Remark 6.2** For the sake of conciseness, only numerical results associated to a given linear polarization for both incident and scattered electromagnetic fields are presented. Nevertheless, we get the same conclusion for the other polarization.

For the experiments, we consider two types of polygonal surface meshes in order to test the V-BEM approach. The first type of mesh consists of triangles and n-sided regular polygons (n > 3), which is constructed from a simple post-treatment of meshes made of triangles or/and quadrangles generated through Gmsh mesher (Geuzaine and Remacle, 2009). We refer to this mesh type as ad-hoc (examples of such a mesh for both sphere and cube are shown in figure 2). The second type of mesh fully consists of convex arbitrary-shape polygons generated through Voro++ library (Rycroft, 2009). As a result, this approach leads to less regular mesh than the previous one, in the sense that an important local heterogeneity between edges lengths can appear. (see, table 1). This type of mesh is instead identified

<table><caption>TABLE 1 Median and maximum of $h_{r} = h_{K}/h_{K}^{min}$ for a representative subset of sphere meshes (left) and all cube surface meshes considered in figure 5 (right).</caption><thead><tr><th colspan="6">Sphere</th><th colspan="6">Cube</th></tr><tr><th colspan="2">triang</th><th colspan="2">ad-hoc</th><th colspan="2">random</th><th colspan="2">triang</th><th colspan="2">ad-hoc</th><th colspan="2">random</th></tr><tr><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th></tr></thead><tbody><tr><td>1.04</td><td>2.44</td><td>1.07</td><td>9.61</td><td>10.13</td><td>560.80</td><td>1.001</td><td>1.431</td><td>1.831</td><td>3.442</td><td>8.434</td><td>19.959</td></tr><tr><td>1.04</td><td>2.72</td><td>1.05</td><td>3.54</td><td>9.92</td><td>4815.31</td><td>1.001</td><td>1.431</td><td>1.327</td><td>3.331</td><td>8.524</td><td>19.995</td></tr><tr><td>1.03</td><td>1.76</td><td>1.05</td><td>3.79</td><td>10.74</td><td>51031.05</td><td>1.001</td><td>1.431</td><td>1.234</td><td>3.875</td><td>8.566</td><td>19.999</td></tr><tr><td>1.02</td><td>2.57</td><td>1.04</td><td>4.28</td><td>10.60</td><td>31520.29</td><td>1.001</td><td>1.431</td><td>1.064</td><td>3.795</td><td>8.722</td><td>20.000</td></tr></tbody></table>
* The monostatic and bistatic RCS errors are defined as, respectively,

$$E_{RCS}^{m} = \left| \sigma(-\hat{\xi}_{I}, \mathbf{J}) - \sigma(-\hat{\xi}_{I}, \Pi_{h}^{0} \mathbf{J}_{h}) \right|, \qquad (6.2)$$

$$E_{RCS}^{b}(\hat{\xi}_{s}) = \left| \sigma(\hat{\xi}_{s}, \mathbf{J}) - \sigma(\hat{\xi}_{s}, \Pi_{h}^{0} \mathbf{J}_{h}) \right|^{2}, \qquad (6.3)$$

where $\hat{\xi}_{I}$ and $\hat{\xi}_{s}$ are the unit vectors associated to the directions of the incident wave and the scattered wave, respectively, and $\sigma$ is related to the scattering amplitude (Nedelec, 2001) and is defined as

$$\sigma(\hat{\xi}_{s}, \mathbf{v}) = \frac{-i\kappa}{\sqrt{4\pi}} \int_{\Gamma_{\text{or}} \Gamma_{h}} e^{i\kappa(x \cdot \hat{\xi}_{s})} (\mathbf{v}(x) - \hat{\xi}_{s}(\mathbf{v}(x) \cdot \hat{\xi}_{s})) \cdot \hat{e}_{r} d\gamma_{x}, \qquad (6.4)$$

with $\hat{e}_{r}$ being the direction of the scattered wave's polarization.

* The bistatic RCS, being a relevant observable in radar applications, characterizes the ability of a body to send back the incident electromagnetic wave energy into different directions and is defined, with $\mathbf{v}$ being $\mathbf{J}$ or $\Pi_{h}^{0} \mathbf{J}_{h}$, as

$$\text{RCS}(\hat{\xi}_{s}) = \left| \sigma(\hat{\xi}_{s}, \mathbf{v}) \right|^{2}. \qquad (6.5)$$

**Remark 6.1** We highlight the fact that $\Pi_{h}^{0} \mathbf{J}_{h} = \mathbf{J}_{h}$ on triangular RT elements, hence the above quantity definitions are modified accordingly.

**Remark 6.2** For the sake of conciseness, only numerical results associated to a given linear polarization for both incident and scattered electromagnetic fields are presented. Nevertheless, we get the same conclusion for the other polarization.

For the experiments, we consider two types of polygonal surface meshes in order to test the V-BEM approach. The first type of mesh consists of triangles and $n$-sided regular polygons ($n > 3$), which is constructed from a simple post-treatment of meshes made of triangles or/and quadrangles generated through Gmsh mesher (Geuzaine and Remacle, 2009). We refer to this mesh type as ad-hoc (examples of such a mesh for both sphere and cube are shown in figure 2). The second type of mesh fully consists of convex arbitrary-shape polygons generated through Voro++ library (Rycroft, 2009). As a result, this approach leads to less regular mesh than the previous one, in the sense that an important local heterogeneity between edges lengths can appear. (see, table 1). This type of mesh is instead identified

<table><caption>TABLE 1 Median and maximum of $h_{r} = h_{K}/h_{K}^{\text{min}}$ for a representative subset of sphere meshes (left) and all cube surface meshes considered in figure 5 (right).</caption><thead><tr><th colspan="6">Sphere</th><th colspan="6">Cube</th></tr><tr><th colspan="2">triang</th><th colspan="2">ad-hoc</th><th colspan="2">random</th><th colspan="2">triang</th><th colspan="2">ad-hoc</th><th colspan="2">random</th></tr><tr><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th><th>h<sub>r</sub><sup>med</sup></th><th>h<sub>r</sub><sup>max</sup></th></tr></thead><tbody><tr><td>1.04</td><td>2.44</td><td>1.07</td><td>9.61</td><td>10.13</td><td>560.80</td><td>1.001</td><td>1.431</td><td>1.831</td><td>3.442</td><td>8.434</td><td>19.959</td></tr><tr><td>1.04</td><td>2.72</td><td>1.05</td><td>3.54</td><td>9.92</td><td>4815.31</td><td>1.001</td><td>1.431</td><td>1.327</td><td>3.331</td><td>8.524</td><td>19.995</td></tr><tr><td>1.03</td><td>1.76</td><td>1.05</td><td>3.79</td><td>10.74</td><td>51031.05</td><td>1.001</td><td>1.431</td><td>1.234</td><td>3.875</td><td>8.566</td><td>19.999</td></tr><tr><td>1.02</td><td>2.57</td><td>1.04</td><td>4.28</td><td>10.60</td><td>31520.29</td><td>1.001</td><td>1.431</td><td>1.064</td><td>3.795</td><td>8.722</td><td>20.000</td></tr></tbody></table>
with the keyword random (examples of such a mesh are shown in figure 2). Let us point out that the techniques of construction employed ensure the planarity of each element within meshes. For the results comparison, we further consider classical triangular meshes, which we refer to them as triang (examples of such a mesh are shown in figure 2) for all BEM-based simulations. Moreover, we can define an average mesh-size $h_{mean}$ for each mesh as follows

$$h_{mean}(T_h) = \frac{1}{N_{T_h}} \sum_{K \in T_h} h_K. \qquad (6.6)$$

**Remark 6.3** We point out that, contrary to the polyhedral geometries considered in the above theoretical sections, the planar surface meshes of the sphere induce approximate surfaces, $\Gamma_h$. As a consequence, a difficulty appears when computing (6.1) since the analytical and discrete solutions have different definition domains. To overcome this issue, in the sphere case, we project the analytical solution onto $\Gamma_h$ as follows:

$$\forall \mathbf{x}_h \in \Gamma_h, \tilde{\mathbf{J}}(\mathbf{x}_h) := \mathbf{J}(\mathbf{x}), \quad \text{where } \mathbf{x} = \arg \min_{z \in \Gamma} |\mathbf{x}_h - \mathbf{z}|, \qquad (6.7)$$

and hence, we substitute $\mathbf{J}$ by the projected analytical solution $\tilde{\mathbf{J}}$ in (6.1).

First, we investigate the convergence rates obtained with V-BEM in both sphere and cube cases. To this aim, we consider a family of successively refined meshes for each triang, ad-hoc and random configuration. The characteristics of a representative subset of such meshes are given in table 1. We can remark in figure 3 and in figure 5 that both V-BEM and BEM approaches convergence with same rates in terms of relative electric current and monostatic RCS errors, regardless the mesh regularity, whether for the sphere or the cube.

Looking at the accuracy, we note that for a given mesh-size, the bistatic RCS results, shown in figure 6, obtained using the V-BEM have the same trends as those of the BEM approach for both the sphere and the cube. We point out that despite the important edge heterogeneity in random meshes, the RCS errors $E_{RCS}^b$ are of the same order. Particularly, for the spherical geometry, the error related

<figure><img src="image_4.png" /><figcaption>FIG. 3. Relative errors $e_r$ on $\mathbf{J}$ as a function of the mesh-size $h_{mean}$ (in m) for the sphere (left) and the cube (right).</figcaption></figure>
with the keyword random (examples of such a mesh are shown in figure 2). Let us point out that the techniques of construction employed ensure the planarity of each element within meshes. For the results comparison, we further consider classical triangular meshes, which we refer to them as triang (examples of such a mesh are shown in figure 2) for all BEM-based simulations. Moreover, we can define an average mesh-size $h_{mean}$ for each mesh as follows

$$h_{mean}(T_h) = \frac{1}{N_{T_h}} \sum_{K \in T_h} h_K. \qquad (6.6)$$

**Remark 6.3** We point out that, contrary to the polyhedral geometries considered in the above theoretical sections, the planar surface meshes of the sphere induce approximate surfaces, $\Gamma_h$. As a consequence, a difficulty appears when computing (6.1) since the analytical and discrete solutions have different definition domains. To overcome this issue, in the sphere case, we project the analytical solution onto $\Gamma_h$ as follows:

$$\forall \mathbf{x}_h \in \Gamma_h, \tilde{\mathbf{J}}(\mathbf{x}_h) := \mathbf{J}(\mathbf{x}), \quad \text{where } \mathbf{x} = \arg \min_{z \in \Gamma} |\mathbf{x}_h - \mathbf{z}|, \qquad (6.7)$$

and hence, we substitute $\mathbf{J}$ by the projected analytical solution $\tilde{\mathbf{J}}$ in (6.1).

First, we investigate the convergence rates obtained with V-BEM in both sphere and cube cases. To this aim, we consider a family of successively refined meshes for each triang, ad-hoc and random configuration. The characteristics of a representative subset of such meshes are given in table 1. We can remark in figure 3 and in figure 5 that both V-BEM and BEM approaches convergence with same rates in terms of relative electric current and monostatic RCS errors, regardless the mesh regularity, whether for the sphere or the cube.

Looking at the accuracy, we note that for a given mesh-size, the bistatic RCS results, shown in figure 6, obtained using the V-BEM have the same trends as those of the BEM approach for both the sphere and the cube. We point out that despite the important edge heterogeneity in random meshes, the RCS errors $E_{RCS}^b$ are of the same order. Particularly, for the spherical geometry, the error related

<figure><img src="image_4.png" /><figcaption>FIG. 3. Relative errors $e_r$ on $\mathbf{J}$ as a function of the mesh-size $h_{mean}$ (in m) for the sphere (left) and the cube (right).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Distribution of the real part of the x-component of $J_h$ (in A·m⁻¹) on triang (BEM) and random (V-BEM) meshes of the sphere (top). Module of the real part of $J_h$ (in A·m⁻¹) on triang (BEM) and ad-hoc (V-BEM) meshes of the cube surface (bottom).</figcaption></figure>

to random meshes is globally better because of the higher number of dofs (162k) in comparison with those associated to triang (72k dofs) and ad-hoc (61k dofs) meshes. On the other hand, we observe that the ad-hoc result is as accurate as the triang one. As for the scattering by the cube, being the

<figure><img src="image_6.png" /><figcaption>FIG. 5. Monostatic RCS errors $E_{RCS}^m$ as a function of the mesh-size $h_{mean}$ (in m) for the sphere (left) and the cube (right).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Distribution of the real part of the *x*-component of $J_h$ (in A·m⁻¹) on triang (BEM) and random (V-BEM) meshes of the sphere (top). Module of the real part of $J_h$ (in A·m⁻¹) on triang (BEM) and ad-hoc (V-BEM) meshes of the cube surface (bottom).</figcaption></figure>

to random meshes is globally better because of the higher number of dofs (162k) in comparison with those associated to triang (72k dofs) and ad-hoc (61k dofs) meshes. On the other hand, we observe that the ad-hoc result is as accurate as the triang one. As for the scattering by the cube, being the

<figure><img src="image_6.png" /><figcaption>FIG. 5. Monostatic RCS errors $E_{RCS}^m$ as a function of the mesh-size $h_{mean}$ (in m) for the sphere (left) and the cube (right).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Bistatic RCS (in dB·m⁻²) and bistatic RCS error $E_{RCS}^b$ as a function of the observation angle (in degree) representing the angular variation of the scattered-field vector $\hat{\xi}_s$ with respect to the incident one $\hat{\xi}_I$, for the sphere (left) and the cube (right) for a given $h_{mean}$ corresponding to: the 5th triang and ad-hoc meshes, the 4th random mesh in figure 5 (sphere), and the finest meshes in figure 5 (cube).</figcaption></figure>

number of dofs similar in the three mesh types (i.e. 395k for triang, 392k for ad-hoc and 372k for random), the three simulations lead to equivalent bistatic RCS errors. Finally, figure 4 shows the same qualitative behavior of the electric current on the sphere and on the cube surface obtained from V-BEM and BEM calculations.

## 7. Conclusions

In this paper, we designed a boundary element method relying on a lowest-order virtual element approximation for the solution of the electric field integral equation on surface meshes composed of polygons or classical elements (e.g., triangles and quadrangles) within which hanging nodes may appear. We established the well-posedness of the new 1stabilization-free discrete weak formulation by proving the *h*-uniform *inf-sup* condition in the natural norm for polyhedral surfaces. Through an *a priori* error analysis, we then demonstrated that the resulting scheme convergences quasi-optimally and, more specifically at the same rate as that of the classical lowest-order RT boundary element scheme on simplicial meshes when considering sufficiently smooth solutions. Finally, numerical simulations of scattering phenomena, involving solutions of different regularities, showed that the proposed scheme behaves, in terms of *h*-convergence and accuracy, as its RT discretization counterpart on simplicial meshes while handling surface meshes that combine a broad range of elements, from, e.g., triangles to convex polygons eventually featuring important edge heterogeneity.

## REFERENCES

S. B. Adrian, A. Dely, D. Consoli, A. Merlini, and F. P. Andriulli. Electromagnetic Integral Equations: Insights in Conditioning and Preconditioning. *IEEE Open Journal of Antennas and Propagation*, 2:1143–1174, 2021.

N. A. Barnafi, F. Dassi, and S. Scacchi. Parallel block preconditioners for virtual element discretizations of the time-dependent Maxwell equations. *Journal of Computational Physics*, 478:111970, Apr. 2023.

L. Beirão da Veiga and L. Mascotto. Interpolation and stability properties of low-order face and edge virtual element spaces. *IMA Journal of Numerical Analysis*, 43(2):828–851, 03 2022.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Bistatic RCS (in dB·m⁻²) and bistatic RCS error $E_{RCS}^b$ as a function of the observation angle (in degree) representing the angular variation of the scattered-field vector $\hat{\xi}_s$ with respect to the incident one $\hat{\xi}_I$, for the sphere (left) and the cube (right) for a given $h_{mean}$ corresponding to: the 5th triang and ad-hoc meshes, the 4th random mesh in figure 5 (sphere), and the finest meshes in figure 5 (cube).</figcaption></figure>

number of dofs similar in the three mesh types (i.e. 395k for triang, 392k for ad-hoc and 372k for random), the three simulations lead to equivalent bistatic RCS errors. Finally, figure 4 shows the same qualitative behavior of the electric current on the sphere and on the cube surface obtained from V-BEM and BEM calculations.

## 7. Conclusions

In this paper, we designed a boundary element method relying on a lowest-order virtual element approximation for the solution of the electric field integral equation on surface meshes composed of polygons or classical elements (e.g., triangles and quadrangles) within which hanging nodes may appear. We established the well-posedness of the new 1stabilization-free discrete weak formulation by proving the *h*-uniform *inf-sup* condition in the natural norm for polyhedral surfaces. Through an *a priori* error analysis, we then demonstrated that the resulting scheme convergences quasi-optimally and, more specifically at the same rate as that of the classical lowest-order RT boundary element scheme on simplicial meshes when considering sufficiently smooth solutions. Finally, numerical simulations of scattering phenomena, involving solutions of different regularities, showed that the proposed scheme behaves, in terms of *h*-convergence and accuracy, as its RT discretization counterpart on simplicial meshes while handling surface meshes that combine a broad range of elements, from, e.g., triangles to convex polygons eventually featuring important edge heterogeneity.

## REFERENCES

S. B. Adrian, A. Dely, D. Consoli, A. Merlini, and F. P. Andriulli. Electromagnetic Integral Equations: Insights in Conditioning and Preconditioning. *IEEE Open Journal of Antennas and Propagation*, 2:1143–1174, 2021.

N. A. Barnafi, F. Dassi, and S. Scacchi. Parallel block preconditioners for virtual element discretizations of the time-dependent Maxwell equations. *Journal of Computational Physics*, 478:111970, Apr. 2023.

L. Beirão da Veiga and L. Mascotto. Interpolation and stability properties of low-order face and edge virtual element spaces. *IMA Journal of Numerical Analysis*, 43(2):828–851, 03 2022.