H Barucq

R Djellouli

A Saint-Guirons

TWO-DIMENSIONAL APPROXIMATE LOCAL DtN BOUNDARY CONDITIONS FOR ELLIPTICAL-SHAPED BOUNDARIES

We propose a new class of approximate local DtN boundary conditions to be applied on elliptical-shaped exterior boundaries when solving acoustic scattering problems by elongated obstacles. These conditions are : (a) exact for the first modes, (b) easy to implement and to parallelize, (c) compatible with the local structure of the computational finite element scheme, and (d) applicable to exterior elliptical-shaped boundaries that are more suitable in terms of cost-effectiveness for surrounding elongated scatterers. We investigate analytically and numerically the effect of the frequency regime and the slenderness of the boundary on the accuracy of these conditions. We also compare their performance to the second order absorbing boundary condition (BGT2) designed by Bayliss, Gunzburger and Turkel when expressed in elliptical coordinates. The analysis reveals that, in the low frequency regime, the new second order DtN condition (DtN2) retains a good level of accuracy regardless of the slenderness of the boundary. In addition, the DtN2 boundary condition outperforms the BGT2 condition. Such superiority is clearly noticeable for large eccentricity values.

Introduction

Exterior Helmholtz problems are classical mathematical models for studying scattering problems arising in many applications such as sonar, radar, geophysical exploration, nondestructive testing, etc. Despite their simplicity, this class of problems is not completely solved particularly from a numerical point of view. For example, the computation of the solutions of these problems requires first to limit it to a finite domain. This is often achieved by surrounding the given scatterer(s) (or radiator) by an artificial boundary that is located at some distance (measured in multiples of wavelength of interest) from its surface. A so-called "nonreflecting" boundary condition is then prescribed on the artificial boundary to represent the "far-field" behavior of the scattered field. The challenge here is the development of a simple but reliable as well as cost-effective computational procedure for representing the far-field behavior of the scattered field. The quest for such conditions is ongoing (see, e.g., the recent review by Turkel in the book [START_REF] Turkel | Iterative methods for the exterior Helmholtz equation including absorbing boundary, conditions[END_REF]).

We propose in this work new two-dimensional approximate local DtN boundary conditions to be employed on elliptical-shaped boundaries that are primary candidates for surrounding elongated scatterers. The idea for constructing such conditions is driven by several considerations, chief among them the following two reasons. First, the widely-used second order absorbing boundary condition (BGT2) designed by Bayliss, Gunzburger and Turkel for circular-shaped boundaries [START_REF] Bayliss | Boundary conditions for the numerical solution of elliptic equations in exterior regions[END_REF] performs poorly when it is expressed in elliptical coordinates and applied to elliptical-shaped boundaries in the low frequency regime [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]. The accuracy deteriorates significantly for large eccentricity values of the boundaries as observed in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]. The damping effect introduced to this condition [START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF] improves the performance for small eccentricity values. However, the modified BGT2 still performs poorly for eccentricity values larger than 0.6 in the (relatively) low frequency regime (see Figure 6 in [START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF]). Hence, there is a need for constructing local absorbing boundary conditions (ABC) that extend the range of satisfactory performance. Second, the two-dimensional approximate local DtN conditions designed for circular-shaped boundaries [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF] outperforms significantly BGT conditions, particularly for low wavenumber values as reported in [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]. However, using these conditions on circular-shaped exterior boundaries when solving scattering problems by elongated scatterers often leads to larger than needed computational domains, which hampers computational efficiency. This suggests that approximate local DtN boundary conditions designed for elliptical-shaped boundaries is an attractive alternative for improving the computational performance.

Given that, this work is devoted to the construction of theses conditions and to the assessment of their performance when employed for solving two-dimensional acoustic scattering problems by elliptical-shaped obstacles. The idea of constructing two-dimensional approximate local DtN boundary conditions is not new. Indeed, as stated earlier, such conditions have been already derived for circularshaped boundaries [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF]. The construction procedure adopted in [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF] is based on the localization of the truncated global DtN boundary condition [START_REF] Keller | Exact nonreflecting boundary conditions[END_REF]. The key ingredient of this procedure is the trigonometric identities that express high order derivatives of sine and cosine functions (see for example Eq. (A4) p. 276 in [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF]). However, this property is not satisfied by the periodic Mathieu functions (see p. 376 in [START_REF] Stratton | Electromagnetic Theory[END_REF]). Consequently, the procedure used in [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF] is no longer applicable to the truncated global DtN boundary operator when expressed in elliptical coordinates [START_REF] Grote | On nonreflecting boundary conditions[END_REF][START_REF] Givoli | Exact representations on artificial interfaces and applications in mechanics[END_REF]. Hence, the construction methodology we propose for deriving the class of approximate local DtN boundary conditions in elliptical coordinates can be viewed as an inverse-type approach. More specifically, we start from a Robin-type boundary condition with unknown coefficients. Unlike the case of polar coordinates, these coefficients depend on the angle θ of the elliptical coordinates. Such dependence is necessary to preserve the symmetry and local nature of the resulting boundary conditions. Then, we require that the considered condition to be an exact representation of the first modes. Consequently, the coefficients are the unique solution of a linear algebraic system.

We assess mathematically and numerically the performance of the constructed approximate local DtN boundary conditions. More specifically, we analyze the effect of low wavenumber and the eccentricity on the performance of these conditions in the case of two-dimensional scattering problems. We adopt the on-surface radiation condition formulation (OSRC) [START_REF] Kriegsmann | A new formulation of electromagnetic wave scattering usin an on-surface radiation boundary condition approach[END_REF] in order to perform analytically this investigation. We note that such formulation is not appropriate for the high frequency regime as observed previously in [START_REF] Antoine | Fast approximate computation of a time-harmonic scattered field using the On-Surface Radiation method[END_REF]. The main interest in the following analyses is to evaluate the performance of the proposed approximate local DtN conditions at low wavenumber to see if relatively small computational domains can be employed in order to avoid excessive computational cost. The OSRC formulation must be viewed as an extreme case while an exterior elliptical-shaped boundary surrounding an elongated scatterer would be less "demanding" on the boundary condition. The analysis herein shows that the constructed second-order local DtN condition retains a good level of accuracy in the low frequency regime for all eccentricity values of the elliptical-shaped boundaries. We must point out that we have also performed a similar investigation when these conditions are employed for solving exterior two-dimensional radiator problems. However, because of space limitation, we do not report in this paper on the results obtained in this case. These results can be found however in [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF].

Preliminaries

The elliptic cylindrical coordinates (ξ, θ, z) are related to the cartesian coordinates (x, y, z) by the transformation: x = a cos θ, y = b sin θ, and z = z where θ ∈ [0, 2π), z ∈ R, and a and b are the semi-major and semi-minor axes, respectively. They are given by a = f cosh ξ, b = f sinh ξ, where f represents the inter-focal distance and ξ is a nonnegative real number. The eccentricity, e on an elliptical cylinder at ξ = ξ 0 is given by e = 1 cosh

ξ 0 = 1 -b 2 a 2 . Note that 0 < e < 1.
Furthermore, any incident plane wave u inc can be expressed in this coordinate system as follows u inc = e ikf cosh ξ(cos θ cos θ 0 +tanh ξ sin θ sin θ 0 ) ,

where the wavenumber k is a positive number and θ 0 is the incident angle. We consider in this work two-dimensional acoustic scattering problems by elliptical-shaped scatterers. Such restriction allows us however to adopt the OSRC formulation and set the elliptical-shaped artificial boundary on top of the boundary of the obstacle. In addition, we assume, for simplicity, that the scatterers are sound-soft [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. Consequently, the acoustic scattered field u scat solution of a sound-soft elliptical-shaped scatterer can then be expressed in terms of Mathieu functions as [START_REF] Bowman | Electromagnetic and Acoustic Scatttering by Simple Shapes[END_REF]

         u scat = - √ 8π ∞ m=0 i m Ce m Re (3) m (kf, cosh ξ)Se m (kf, cos θ) + ∞ m=1 i m Co m Ro (3) m (kf, cosh ξ)So m (kf, cos θ) , (2)
where Ce m and Co m are the Fourier coefficients [START_REF] Bowman | Electromagnetic and Acoustic Scatttering by Simple Shapes[END_REF][START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF][START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF], Re m ) represents the even (resp. odd) normalization factor of the periodic Mathieu functions (see p. 378 in Reference [START_REF] Stratton | Electromagnetic Theory[END_REF]). We recall (for example) that the mth elliptical even mode is given by [START_REF] Senior | Scalar diffraction by a prolate spheroid at low frequencies[END_REF] u m = Re (3) m (kf, cosh ξ)Se m (kf, cos θ) ; m ≥ 0.

(3)

Two-dimensional approximate local DtN boundary conditions in elliptical coordinates

The two-dimensional first-(DtN1) and second-order (DtN2) local Dirichletto-Neumann boundary conditions, defined on the elliptical-shaped boundary ξ = ξ 0 , are given by DtN1 :

∂u ∂ξ = √ 1 -e 2 e α 0 u, (4)
DtN2 :

∂u ∂ξ = √ 1 -e 2 e(a 0 -a 1 ) [(a 0 α 1 -a 1 α 0 )u + (α 1 -α 0 )D θ u] , (5)
where the constants a m are the characteristic numbers of the mth even periodic Mathieu functions Se m [START_REF] Abramovitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF]. The coefficients α m are given by

α m = Re (3) m ′ (eka, e -1 )
Re (3) m (eka, e -1 ) ; m ≥ 0,

and the operator D θ is the following second order differential operator

D θ = ∂ 2 ∂θ 2 - (eka) 2 2 (cos 2θ) I.
The following four remarks are noteworthy:

• First, the boundary conditions ( 4), ( 5) are symmetric and local boundary conditions. The local feature of these conditions is of a great interest from a numerical view point. Indeed, the incorporation of these conditions in any finite element code introduces only mass-and stiffness-type matrices defined on the exterior boundary. The characteristic numbers a n of the nth even periodic Mathieu functions and the coefficients α n can be computed once for all at the preprocessing level. • Second, it must be pointed out that when e = 0 (the ellipse degenerates to a circle), conditions ( 4) and ( 5) are identical to the two-dimensional DtN conditions designed for circular-shaped boundaries [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF][START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]. This property can be easily established [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF] using the asymptotic behavior of the even radial Mathieu functions of the third type Re 21) and ( 22), p. 171 and Eq. (37), p. 200 in Ref. [START_REF] Meixner | Matieusche funktionen and spharoidfunktionen[END_REF]).

• Last, approximate local DtN conditions of order higher than two are inappropriate for conventional finite element implementations since they require regularity higher than C 0 . Consequently, we consider in the following only the first-(DtN1) and second-order (DtN2) boundary conditions. Moreover, the boundary conditions ( 4) and ( 5) satisfy, by construction, the following property that one can easily verify: Lemma 3.1 The first-order DtN boundary condition (DtN1) given by ( 4) is an exact representation of the first mode u 0 given by Eq.( 3). The second-order DtN boundary condition (DtN2) given by ( 5) is an exact representation of the first even mode u 0 and the second even mode u 1 given by Eq.( 3).

Remark 3.2

We have also constructed a second-order boundary condition that is an exact representation of the first even mode and the first odd mode. The expression of the obtained condition is similar to the condition [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF] [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF]. The only difference is the replacement of (a) the even radial Mathieu function of the third type Re [START_REF] Abramovitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF]. Similarly, we have also constructed a second-order boundary condition that is an exact representation, in the least-squares sense, of the first two even mode and the first odd mode [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF]. We have analyzed the performance of these conditions and found that it is identical to the second-order boundary condition (DtN2) given by Eq. ( 5) as reported in [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF]. Given that, along with space limitation, we only consider in this paper the conditions given by Eq. ( 4) and Eq. ( 5) to illustrate the performance of this class of ABCs.

In the following, we assess analytically and numerically the performance of the approximate local DtN1 and DtN2 boundary conditions given by ( 4) and [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF]. More specifically, we analyze the effect of low wavenumber ka and the eccentricity e on the performance of DtN1 and DtN2 in the case of sound-soft scattering problems. We adopt the on-surface radiation condition formulation (OSRC) [START_REF] Kriegsmann | A new formulation of electromagnetic wave scattering usin an on-surface radiation boundary condition approach[END_REF] in order to perform this investigation analytically. As in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF][START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF][START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF], we assess the performance of the ABCs DtN1 and DtN2 using the specific impedance introduced in [START_REF] Geers | Doubly asymptotic approximations for transient motions of submerged structures[END_REF][START_REF] Geers | Third-order doubly asymptotic approximations for computational acoustics[END_REF] as a practical tool for measuring the efficiency of ABCs in the context of the OSRC formulation. This non-dimensional quantity measures the effect of the truncated medium in physical terms. It provides a convenient indicator of the performance of a given approximate representation. In the elliptical coordinates system, the specific impedance can be expressed as

z = i √ 1 -e 2 ka u scat ∂ ∂ξ (u scat )| ξ=ξ 0 . (7)
Therefore, the specific impedance z exact corresponding to the exact solution for two-dimensional sound-soft acoustic scattering problems can be computed analytically using u scat = -u inc (see Eq. ( 1)), and the Fourier series given by Eq. ( 2) to evaluate ∂ ∂ξ (u scat ) | ξ=ξ 0 (see, for more details, [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF][START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF][START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF]).

Mathematical results

The following lemma states the expression of the approximate specific impedances on the boundary of a sound-soft elliptical-shaped scatterer at ξ = ξ 0 . These expressions can be easily derived from Eq. ( 7) by using u scat = -u inc (see Eq. ( 1)), along with substituting u = -u inc into the DtN conditions (4) and [START_REF] Givoli | Nonreflecting boundary conditions for elastic waves[END_REF]. Lemma 4.1 The approximate specific impedance (z DtN1 ) corresponding to the first order DtN boundary condition (DtN1) is given by

z DtN1 = ika e α 0 . (8)
The approximate specific impedance (z DtN2 ) corresponding to the second order DtN boundary condition (DtN2) is given by

z DtN2 = ika e (a 0 -a 1 ) (a 0 α 1 -a 1 α 0 ) -(α 1 -α 0 ) ika λ + (ka) 2 ∂λ ∂θ 2 + e 2 (ka) 2 2 cos(2θ) , (9
) where the constants a m are the characteristic numbers of the mth even periodic Mathieu functions Se m [START_REF] Abramovitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF], the coefficients α m are given by Eq. ( 6), and the parameter λ is given by

λ = cos θ cos θ 0 + √ 1 -e 2 sin θ sin θ 0 . ( 10
)
Remark 4.2 First, we note that when e = 0, that is the ellipse degenerates to a circle, the approximate DtN specific impedances given by Eqs. ( 8) and ( 9) are identical to the ones obtained in the case of circular-shaped scatterers and an angle of incident plane wave θ = θ 0 [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]. Second, we observe that, unlike the exact specific impedance z exact , the DtN1 specific impedance z DtN1 does not depend on the observation angle θ. This result suggests that the DtN1 boundary condition is less accurate than the DtN2.

The next result states the behavior of the specific impedances as ka → 0. Note that the behavior of z exact has been established in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF][START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF], while the behavior of z DtN2 is a direct consequence, as shown in [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF], of using the asymptotic behavior of the even radial Mathieu functions of the third type Re 21) and ( 22), p. 171 and Eq. (37), p. 200 in Ref. [START_REF] Meixner | Matieusche funktionen and spharoidfunktionen[END_REF]).

Proposition 4.3

The asymptotic behavior of the specific impedances of the scattered field on the boundary of an elliptical-shaped obstacle as ka → 0 is given by

z DtN2 ∼ z exact ∼ π ka 2 + ika ln(ka). (11)
Remark 4.4 Unlike the case of BGT2 boundary condition in elliptical coordinates (see Eq. ( 40), p. 3637 in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF], Proposition 4.3 reveals that the asymptotic behaviors of z exact and z DtN2 are identical. Note that this behavior is identical to the case of circular-shaped scatterers [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF].

Illustrative numerical results

We have performed several numerical experiments to assess the performance of the constructed approximate local DtN boundary conditions. However, due the space limitation, we report here the results of one experiment for illustration. These results are obtained for an incident angle θ 0 = 0, two wavenumber values ka = 0.1, 1, and three eccentricity values e = 0.1, 0.6, 0.9. Note that we have computed the specific impedances using the integral representations of the Mathieu functions given by Eqs. (20.7.25)-(20.7.27) p. 737 in Reference [START_REF] Abramovitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF] , and the Mathematica software package. We have used up to 20 terms in the Fourier series (see Eq. ( 2)). The results depicted in Figure 1 clearly illustrates the superiority of DtN2 over BGT2. Indeed, these results indicate that the performance of DtN2 at low wavenumber values is accurate, as expected from Proposition 4.3. In addition, they show that, unlike BGT2 , the accuracy is not sensitive to the eccentricity e. Indeed, one can observe that the accuracy of BGT2 deteriorates significantly when e ≥ 0.6, as already reported in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]. (crossed) and the BGT2 (dashed) for the incident angle θ 0 = 0, and for ka = 0.1 (left) and ka = 1 (right).

Conclusion

We have designed a new class of approximate local ABCs to be applied on elliptical-shaped exterior boundaries when solving acoustic scattering problems by elongated obstacles. These conditions are exact for the first radiation modes, they are easy to implement and to parallelize, and they preserve the local structure of the computational finite element scheme. The analysis reveals that the new second-order boundary condition DtN2 outperforms the BGT2 condition, and extends the range of satisfactory performance to all eccentricity values.

  the even (resp. odd) radial Mathieu functions of the j th type (see p. 378 in Reference[START_REF] Stratton | Electromagnetic Theory[END_REF]), Se m (resp. So m ) represents the even (resp. odd) periodic Mathieu functions (see p. 377 in Reference[START_REF] Stratton | Electromagnetic Theory[END_REF]), and N (e) m (resp. N (o)

( 3 )

 3 m (see Eqs. (

  (b) the characteristic number a 1 by the one corresponding to the odd periodic Mathieu function So 1

Figure 1 :

 1 Figure 1: Absolute value of the specific impedance for the exact (solid), the DtN2

Acknowledgments

The authors acknowledge the support by INRIA/CSUN Associate Team Program. Any opinions, findings, conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of INRIA or CSUN.


Florent Chave
email: florent.chave@inria.fr

Daniele A Di

Pietro †2

Simon Lemaire
email: simon.lemaire@inria.fr

A discrete Weber inequality on three-dimensional hybrid spaces with application to the HHO approximation of magnetostatics

Keywords: Weber inequalities, Hybrid spaces, Polyhedral meshes, Hybrid High-Order methods, Magnetostatics. AMS Subject Classification: 65N08, 65N12, 65N30

We prove a discrete version of the first Weber inequality on three-dimensional hybrid spaces spanned by vectors of polynomials attached to the elements and faces of a polyhedral mesh. We then introduce two Hybrid High-Order methods for the approximation of the magnetostatics model, in both its (first-order) field and (second-order) vector potential formulations. These methods are applicable on general polyhedral meshes, and allow for arbitrary orders of approximation. Leveraging the previously established discrete Weber inequality, we perform a comprehensive analysis of the two methods. We finally validate them on a set of test-cases.

Introduction

Let Ω ⊂ R 3 denote an open, bounded, and connected polyhedral domain. In the study of problems in electromagnetism, Weber inequalities [START_REF] Weber | A local compactness theorem for Maxwell's equations[END_REF] constitute a very powerful tool. They can be viewed as a generalization of the celebrated Poincaré inequality to the case of vector fields belonging to (curl; Ω) ∩ (div; Ω), and featuring either vanishing tangential component (first Weber inequality), or vanishing normal component (second Weber inequality) on the boundary Ω of the domain. We refer the reader to [2, Theorems 3.4.3 and 3.5.3] for a general (from a topological viewpoint) statement of Weber inequalities.

We next state the (continuous) first Weber inequality. For the sake of simplicity, we assume from now on that Ω is simply-connected and that Ω is connected. Under these assumptions, the first and second Betti numbers of Ω are both zero, i.e., Ω does not have tunnels and does not enclose any void. For a deeper insight into the role of the different topological assumptions we make on the domain, we refer to Remark 7. For any measurable set ⊂ Ω, we irrespectively denote by (•, •) and || • || the usual inner products and norms on the scalar-valued space 2 ( ) and on the vector-valued spaces 2 ( ; R ℓ ) for ℓ ∈ {2, 3}. We also set, letting denote the unit normal vector field on Ω pointing out of Ω, 0 (curl; Ω) := { ∈ (curl; Ω) : ×( × ) = 0 on Ω}, and (div 0 ; Ω) := { ∈ (div; Ω) : div = 0 in Ω} = ∈ 2 (Ω; R 3 ) : ( , grad ) Ω = 0 ∀ ∈ 1 0 (Ω) .

The following 2 (Ω; R 3 )-orthogonal decomposition holds true (cf. [2, Proposition 3.7.2]):

0 (curl; Ω) = grad 1 0 (Ω) ⊥ ⊕ 0 (curl; Ω) ∩ (div 0 ; Ω) . (1)
With the assumptions we have made on the topology of the domain Ω, the first Weber inequality reads: For any ∈ 0 (curl; Ω) ∩ (div 0 ; Ω),

Ω ≤ W curl Ω , (2)
for some constant W > 0 only depending on the domain Ω. In this work, we derive a discrete version of the first Weber inequality (2) on (three-dimensional) hybrid spaces spanned by vectors of polynomials attached to the elements and faces of a (polyhedral) mesh, as they can be encountered in Hybridizable Discontinuous Galerkin [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second-order elliptic problems[END_REF] and related [START_REF] Wang | A weak Galerkin finite element method for second-order elliptic problems[END_REF] methods, or in Hybrid High-Order (HHO) [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF][START_REF] Di Pietro | A Hybrid High-Order locking-free method for linear elasticity on general meshes[END_REF] methods; see [START_REF] Cockburn | Bridging the Hybrid High-Order and Hybridizable Discontinuous Galerkin methods[END_REF] for a discussion highlighting the analogies and differences between these two families of methods in the context of scalar variable diffusion. The corresponding result is stated in Theorem 3. The proof extends the general ideas used in [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF]Lemma 2.15] to derive a discrete Poincaré inequality on hybrid spaces. A related (yet weaker) result for (non-hybrid) broken polynomial spaces on tetrahedral meshes has been derived in [START_REF] Chen | A priori and computable a posteriori error estimates for an HDG method for the coercive Maxwell equations[END_REF]Lemma 4.1]. Poincaré-Friedrichs inequalities for complexes of discrete distributional differential forms (introduced in [START_REF] Licht | Complexes of discrete distributional differential forms and their homology theory[END_REF] as a generalization of the distributional finite element de Rham sequences of [START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF]Section 3]) are also proved in the recent work [START_REF] Christiansen | Poincaré-Friedrichs inequalities of complexes of discrete distributional differential forms[END_REF]. Complexes of discrete distributional differential forms express a notion of compatibility that is relieved from global conformity requirements through the presence of "jump terms" in the distributions (with conforming finite elements corresponding to the special case where such terms vanish). A major difference with respect to the present contribution is that we work in a non-compatible setting, i.e., there is no underlying exact discrete complex. Specifically, the sole result related to (discrete) compatibility leveraged in the proof of the discrete Weber inequality is a decomposition of vector-valued polynomial functions inside each mesh element (see Lemma 2 below). Another important difference between the two works is that our results naturally apply to general polyhedral meshes.

In the second part of this paper, we tackle the HHO approximation of magnetostatics, in both its (first-order) field formulation and (second-order, generalized) vector potential formulation. Various discretization methods have been studied in the literature to approximate the magnetostatics equations (or, more generally, Maxwell equations). Conforming finite element discretizations were originally proposed (on tetrahedra, essentially) in the seminal work of Nédélec [START_REF] Nédélec | Mixed finite elements in R 3[END_REF][START_REF] Nédélec | A new family of mixed finite elements in R 3[END_REF] (cf. also [START_REF] Monk | Finite element methods for Maxwell's equations[END_REF]). We also mention [START_REF] Arnold | Finite element exterior calculus[END_REF], in which a unified presentation of conforming finite element methods based on notions from algebraic topology is provided. Nonconforming discretizations include, on standard meshes, the Discontinuous Galerkin method of [START_REF] Perugia | Stabilized interior penalty methods for the time-harmonic Maxwell equations[END_REF] as well as the Hybridizable Discontinuous Galerkin method of [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF][START_REF] Chen | Analysis of a hybridizable discontinuous Galerkin method for the Maxwell operator[END_REF] and, on general polyhedral meshes, the variant [START_REF] Mu | A weak Galerkin finite element method for the Maxwell equations[END_REF] of [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF] and the method of [START_REF] Chen | A superconvergent HDG method for the Maxwell equations[END_REF]; see also [START_REF] Du | A unified error analysis of Hybridizable Discontinuous Galerkin methods for the static Maxwell equations[END_REF]. Methods that support general polyhedral meshes and are built upon discrete spaces that mimick the continuity properties of the spaces appearing in the continuous weak formulation include the Virtual Element methods of [START_REF] Beirão Da Veiga | Lowest-order virtual element approximation of magnetostatic problems[END_REF][START_REF] Beirão Da Veiga | A family of three-dimensional virtual elements with applications to magnetostatics[END_REF][START_REF] Beirão Da Veiga | Virtual element approximations of the vector potential formulation of magnetostatic problems[END_REF], and the fully discrete method of [START_REF] Di Pietro | An arbitrary-order method for magnetostatics on polyhedral meshes based on a discrete de Rham sequence[END_REF] based on the discrete de Rham sequence of [START_REF] Di Pietro | Fully discrete polynomial de Rham sequences of arbitrary degree on polygons and polyhedra[END_REF] (see also [START_REF] Di Pietro | An arbitrary-order discrete de Rham complex on polyhedral meshes: Exactness, Poincaré inequalities, and consistency[END_REF] for recent developments including error estimates). All the hybridized nonconforming methods cited above deal with the approximation of magnetostatics under its (generalized) vector potential formulation. In this paper, we first study an HHO method (which has been briefly introduced in [START_REF] Chave | A three-dimensional Hybrid High-Order method for magnetostatics[END_REF]) for magnetostatics under its field formulation. We take advantage of the fact that the corresponding problem is first-order to avoid locally reconstructing a discrete curl operator as it has to be done for second-order problems (cf. Remark 11). Doing so, we propose a computationally inexpensive and easy-to-implement method. Second, we study an HHO method for magnetostatics under its (generalized) vector potential formulation, that can be seen as a computationally cheaper variant of the method introduced in [START_REF] Chen | A superconvergent HDG method for the Maxwell equations[END_REF] (cf. Remark 17). Our two HHO methods are applicable on general polyhedral meshes, and allow for an arbitrary order of approximation ≥ 0 with proved energy-error of order + 1 (cf. Theorems 14 and 23). Leveraging the previously established discrete Weber inequality, we carry out a comprehensive analysis of the methods, and validate them on a set of test-cases.

The article is organized as follows. In Section 2, we prove the discrete Weber inequality. Then, in Sections 3.1 and 3.2, we tackle the HHO approximation of magnetostatics under its field and vector potential formulations, respectively.

A discrete Weber inequality on hybrid spaces

Discrete setting

We consider a polyhedral mesh M ℎ = (T ℎ , F ℎ ), that is assumed to belong to a regular mesh sequence in the sense of [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF]Definition 1.9]. The set T ℎ is a finite collection of nonempty, disjoint, open polyhedra (called mesh elements) such that Ω = ∈ T ℎ . For convenience reasons that will be made clear in Remark 8 below, we henceforth assume that any mesh element ∈ T ℎ is star-shaped with respect to some interior point . Note that this assumption is not necessary for the results of this article to hold true. The subscript ℎ refers to the meshsize, defined by ℎ := max ∈ T ℎ ℎ , where ℎ denotes the diameter of the element . The set F ℎ collects the planar mesh faces and, for all ∈ T ℎ , we denote by F the set of faces that lie on the boundary of . Boundary faces lying on Ω are collected in the set F b ℎ , and we denote by F i ℎ ≔ F ℎ \ F b ℎ the set of interfaces. For all ∈ F ℎ , we let ℎ denote its diameter and, for all ∈ T ℎ and all ∈ F , denote the unit normal vector to pointing out of . We recall that, since M ℎ belongs to a regular mesh sequence, for all ∈ T ℎ , the quantity card(F ) is bounded from above uniformly in ℎ and, for all ∈ T ℎ and all ∈ F , ℎ is uniformly comparable to ℎ (cf. [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF]Lemma 1.12]). In what follows, we will use the symbol to indicate that an estimate is valid up to a multiplicative constant that may depend on the mesh regularity parameter, the ambient dimension, and (if need be) the polynomial degree, but that is independent of ℎ.

Hybrid spaces

For ∈ N and ( , ℓ) ∈ {(F ℎ , {2}), (T ℎ , {3})}, we respectively denote by P ( ) and P ( ) the vector spaces of ℓ-variate, scalar-valued and R ℓ -valued polynomial functions on of total degree at most . For future use, for any ∈ T ℎ , we let G ( ) := grad P +1 ( ) and R ( ) := curl P +1 ( ) (the notation R standing for "rot"), and we recall that the following (nonorthogonal) decomposition holds true:

P ( ) = G ( ) ⊕ ( -)× R -1 ( ), (3)
with the convention that R -1 ( ) := {0}. For any ∈ F ℎ , we also let G ( ) := grad P +1 ( ) ⊆ P ( ) denote the space of (tangential) gradients of polynomials of total degree up to + 1 on . Finally, we define the broken spaces

P (T ℎ ) := ∈ 2 (Ω) : | ∈ P ( ) ∀ ∈ T ℎ , P (T ℎ ) := ∈ 2 (Ω; R 3 ) : | ∈ P ( ) ∀ ∈ T ℎ ,
as well as the broken subspaces G (T ℎ ) := grad ℎ P +1 (T ℎ ) and R (T ℎ ) := curl ℎ P +1 (T ℎ ) , where grad ℎ (resp. curl ℎ ) denotes the usual broken grad (resp. curl) operator on 1 (T ℎ ) (resp. (curl; T ℎ )).

Let an integer polynomial degree ≥ 0 be given. We introduce the following (global) hybrid spaces:

X +1 ℎ := v ℎ = (v ) ∈ T ℎ , (v ) ∈ F ℎ : v ∈ P +1 ( ) ∀ ∈ T ℎ v ∈ G +1 ( ) ∀ ∈ F ℎ , (4a)
Y +1 ℎ := q ℎ = (q ) ∈ T ℎ , (q ) ∈ F ℎ : q ∈ P ( ) ∀ ∈ T ℎ q ∈ P +1 ( ) ∀ ∈ F ℎ , (4b)
as well as their subspaces incorporating homogeneous essential boundary conditions:

X +1 ℎ,0 := v ℎ ∈ X +1 ℎ : v = 0 ∀ ∈ F b ℎ , Y +1 ℎ,0 := q ℎ ∈ Y +1 ℎ : q = 0 ∀ ∈ F b ℎ . (5)
(Semi)norms on the above spaces are defined in ( 16) and (40) (or (68b)) below. In (4a), the face unknowns for the vectorial variable lie in a strict subspace of P +1 ( ), that is G +1 ( ). This choice is driven by stability purposes. Since the curl operator has a kernel on P +1 ( ) that is composed of all polynomials in G +1 ( ) (this will be justified by Lemma 2 below), we must be able to control the (jump of the) tangential traces of all functions in that space at interfaces, which is the reason why the face unknowns must (at least) lie in G +1 ( ). As far as (4b) is concerned, the necessity to consider face unknowns in P +1 ( ) for the scalar variable is driven by the necessity to reconstruct a grad operator in P +1 ( ) (in turn driven by the necessity to test this discrete grad operator against element unknowns for the vectorial variable themselves in P +1 ( )); cf. Section 2.3 below. Given a mesh element ∈ T ℎ , we respectively denote by X +1 and Y +1 the restrictions of X +1 ℎ and Y +1 ℎ to , and by v := v , (v ) ∈ F ∈ X +1 and q := q , (q ) ∈ F ∈ Y +1 the respective restrictions of generic vectors of polynomials v ℎ ∈ X +1 ℎ and q ℎ ∈ Y +1 ℎ . Also, we let v ℎ and q ℎ (not underlined) be the broken polynomial functions in P +1 (T ℎ ) and in P (T ℎ ) such that

(v ℎ ) | := v
and (q ℎ ) | := q for all ∈ T ℎ .

We next define the interpolators on X +1 ℎ and Y +1 ℎ . To this purpose, we need some preliminary definitions. For any open set ⊆ Ω and any ∈ F ℎ such that ⊂ , , ( ) ∈ 2 ( ; R 2 ) denotes the tangential trace on of ∈ 1 ( ; R 3 ). For ∈ N and ( , ℓ) ∈ {(F ℎ , {2}), (T ℎ , {3})}, we respectively denote by P, and P, the 2 ( )-orthogonal projector onto P ( ) and the 2 ( ; R ℓ )-orthogonal projector onto P ( ). Similarly, G, stands for the 2 ( ; R 2 )-orthogonal projector onto G ( ) and, for S ∈ {G, R}, S, for the 2 ( ; R 3 )-orthogonal projector onto S ( ).

We finally let I +1 X,ℎ : 1 (Ω; R 3 ) → X +1 ℎ and I +1 Y,ℎ : 1 (Ω) → Y +1 ℎ be such that, for all ∈ 1 (Ω; R 3 ) and all ∈ 1 (Ω),

I +1 X,ℎ := +1 P, ( | ) ∈ T ℎ , +1 G, , ( ) ∈ F ℎ , (6a)
I +1 Y,ℎ := P, ( | ) ∈ T ℎ , +1 P, ( | ) ∈ F ℎ . (6b)
For future use, we close this section by introducing the global 2 (Ω)-orthogonal projector P,ℎ , and 2 (Ω; R 3 )-orthogonal projectors P,ℎ and S,ℎ onto, respectively, P (T ℎ ), P (T ℎ ) and S (T ℎ ).

Gradient reconstruction in Y +1

ℎ

We define the global discrete gradient reconstruction operator +1 ℎ : Y +1 ℎ → P +1 (T ℎ ) such that its local restriction +1 : Y +1 → P +1 ( ) to any ∈ T ℎ solves the following problem: For all q ∈ Y +1 ,

+1 q , = -(q , div ) + ∈ F (q , | • ) ∀ ∈ P +1 ( ). (7)
By the Riesz representation theorem in P +1 ( ) equipped with the 2 ( ; R 3 )-inner product, +1 q is uniquely defined. We note the following commutation property (see, e.g., [17, Section 4.2.1]): For all ∈ 1 (Ω), we have

+1 ℎ I +1 Y,ℎ = +1 P,ℎ (grad ). (8)
We also have the following result in the tetrahedral case, which can be exploited to simplify the HHO schemes of Sections 3.1 and 3.2 on matching tetrahedral meshes; see Remarks 9 and 24 below.

Lemma 1 (Norm +1 ℎ

• Ω ). Let T ℎ be a matching tetrahedral mesh. Then, the map

+1 ℎ • Ω defines a norm on Y +1 ℎ,0 .
Proof. Let T ℎ be a matching tetrahedral mesh, and q ℎ ∈ Y +1 ℎ,0 be such that +1 ℎ q ℎ Ω = 0. Then, for all ∈ T ℎ , enforcing +1 q = 0 in the definition [START_REF] Brezzi | Two families of mixed finite elements for second order elliptic problems[END_REF] of this quantity and integrating by parts yields

grad q , + ∈ F q -q | , | • = 0 ∀ ∈ P +1 ( ). (9)
Let N ( ) := G -1 ( ) ⊕ ( -)× R -1 ( ) (with the convention that G -1 ( ) := {0}) denote the Nédélec space of the first kind of degree on (cf. [START_REF] Nédélec | Mixed finite elements in R 3[END_REF]), and let ∈ P +1 ( ) be such that

( , ) = (grad q , ) ∀ ∈ N ( ), (10a)
( | • , ) = (q -q | , ) ∀ ∈ P +1 ( ), ∀ ∈ F . (10b)
The system [START_REF] Chen | HDG and CG methods for the indefinite time-harmonic Maxwell's equations under minimal regularity[END_REF] uniquely defines as a function of the Nédélec space of the second kind of degree +1 on , that is P +1 ( ) (cf. [START_REF] Nédélec | A new family of mixed finite elements in R 3[END_REF][START_REF] Brezzi | Two families of mixed finite elements for second order elliptic problems[END_REF]). Testing (9) with , and using that grad q ∈ G -1 ( ) ⊆ N ( ) and that q -q | ∈ P +1 ( ) for all ∈ F , we infer from (10) that grad q 2 + ∈ F q -q | 2 = 0. Reproducing the same reasoning on all ∈ T ℎ , and using that q ℎ belongs to the space Y +1 ℎ,0 with strongly enforced boundary conditions, finally yields q ℎ = 0 ℎ .

Discrete Weber inequality

Let us begin with a preliminary technical result.

Lemma 2 (Polynomial decomposition). Let ∈ N. Let ∈ T ℎ and ∈ P ( ). Then, there exist ∈ P +1 ( ) and ∈ P ( ) such that = grad + ( -) × curl and

-grad ≤ 2ℎ curl . ( 11
)
Proof. The decomposition = grad + ( -) × curl directly follows from [START_REF] Beirão Da Veiga | A family of three-dimensional virtual elements with applications to magnetostatics[END_REF]. Using the fact that |( -

)× curl | ≤ | -| | curl | ≤ ℎ | curl |, we infer -grad = ( -)× curl ≤ ℎ curl . ( 12
)
We now focus on the right-hand side of [START_REF] Chen | A superconvergent HDG method for the Maxwell equations[END_REF]. Since curl grad = 0, we have curl = curl ( -)× curl .

Using the vector calculus identity curl( × ) = (div ) -[

( -) div(curl ) -[( -)• grad] (curl ) -3 curl + [(curl )• grad] ( -),
where we have used the fact that the div of the curl is zero in the cancellation. Hence, observing that

[(curl )• grad] ( -) = 3 =1 (curl ) 1≤ ≤3 = 3 =1 (curl ) 1≤ ≤3
= curl (with • denoting the -th component of the vector • and the Kronecker symbol),

curl = -2 curl -[( -)• grad] (curl ). (13)
Since

[( -)• grad] (curl ) = 3 =1 ( -) (curl ) 1≤ ≤3
, multiplying (13) bycurl and integrating over , we get

-curl , curl = 2 curl 2 + -, grad | curl | 2 2 = 1 2 curl 2 + ∈ F 1 2 | curl | 2 | , ( -) | • ≥ 1 2 curl 2 ≥ 0, (14)
where we have used an integration by parts formula to pass to the second line, and the fact that is star-shaped with respect to to conclude. From ( 14) and a Cauchy-Schwarz inequality, we infer

curl ≤ 2 curl . (15)
Plugging ( 15) into ( 12), [START_REF] Chen | A priori and computable a posteriori error estimates for an HDG method for the coercive Maxwell equations[END_REF] follows.

We equip the space X +1 ℎ with the seminorm • X,ℎ defined by

v ℎ 2 X,ℎ := curl ℎ v ℎ 2 Ω + ∈ T ℎ ∈ F ℎ -1 +1 G, , (v ) -v 2 . ( 16
)
This seminorm is the discrete counterpart of the 2 (Ω; R 3 )-norm of the curl and is composed of two contributions: the first one is the 2 (Ω; R 3 )-norm of the broken curl of the function obtained patching the polynomials attached to mesh elements; the second one accounts for the difference between the face unknowns and the (gradient part of the) tangential trace of the element unknowns. The kernel of this seminorm is made clear by Remark 4. On the discrete counterpart of the space 0 (curl; Ω) ∩ (div 0 ; Ω) defined by

• X +1 ℎ,0 := w ℎ ∈ X +1 ℎ,0 : w ℎ , +1 ℎ q ℎ Ω = 0 ∀q ℎ ∈ Y +1 ℎ,0 , (17)
the following discrete Weber inequality holds true. Here and in what follows, the circle overset will be reserved to those hybrid spaces that incorporate a discrete divergence-free property (see, in particular, (39) and (67) below).

Theorem 3 (Discrete Weber inequality).

There exists a constant W > 0 independent of ℎ such that,

for all v ℎ ∈ • X +1 ℎ,0 , one has v ℎ Ω ≤ W v ℎ X,ℎ , (18)
where we remind the reader that the term in the left-hand side is the 2 (Ω; R 3 )-norm of the broken polynomial vector v ℎ such that

(v ℎ ) | = v for all ∈ T ℎ . Proof. Let v ℎ ∈ • X +1 ℎ,0 .
Recall the following standard 2 (Ω; R 3 )-orthogonal Helmholtz decomposition (cf., e.g., [2, Proposition 3.7.1]):

2 (Ω; R 3 ) = grad 1 0 (Ω) ⊥ ⊕ (div 0 ; Ω).
By the characterization of divergence-free functions from [2, Theorem 3.4.1], and since v ℎ ∈ P +1 (T ℎ ) ⊂ 2 (Ω; R 3 ) and Ω is connected, we can write

v ℎ = grad + curl , (19)
for some ∈ 1 0 (Ω), and some ∈ 1 (Ω; R 3 ) such that

∫ Ω = 0 and div = 0. Furthermore, 1 (Ω;R 3 ) curl Ω . ( 20
)
Using the 2 (Ω; R 3 )-orthogonal decomposition ( 19) of v ℎ , we have that

v ℎ 2 Ω = (v ℎ , grad ) Ω + v ℎ , curl Ω =: I 1 + I 2 . (21)
For the first term in [START_REF] Di Pietro | A Hybrid High-Order locking-free method for linear elasticity on general meshes[END_REF], setting φ

ℎ := I +1 Y,ℎ ∈ Y +1
ℎ,0 , we infer that

I 1 = v ℎ , +1 P,ℎ (grad ) Ω = v ℎ , +1 ℎ φ ℎ Ω = 0, (22)
where we have used that v ℎ ∈ P +1 (T ℎ ) to insert +1 P,ℎ , followed by the commutation property (8) of

+1

ℎ and the fact that v ℎ ∈ • X +1 ℎ,0 . Let us now estimate the second term in [START_REF] Di Pietro | A Hybrid High-Order locking-free method for linear elasticity on general meshes[END_REF]. We have

I 2 = ∈ T ℎ v , curl = ∈ T ℎ curl v , - ∈ F | × , v | = ∈ T ℎ curl v , - ∈ F , ( × ), , (v ) -v ,
where we have used an integration by parts formula on each mesh element ∈ T ℎ in the first line, and the fact that the jumps of ∈ 1 (Ω; R 3 ) vanish on interfaces along with v = 0 for all ∈ F b ℎ to insert v into the second term in the second line. Applying Cauchy-Schwarz inequalities to the right-hand side, we obtain

I 2 ≤ ∈ T ℎ curlv 2 + ∈ F ℎ -1 , (v ) -v 2 1 /2 × ∈ T ℎ 2 + ∈ F ℎ | × 2 1 /2 . ( 23
)
We focus on the first factor on the right-hand side of [START_REF] Du | A unified error analysis of Hybridizable Discontinuous Galerkin methods for the static Maxwell equations[END_REF]. For ∈ T ℎ and ∈ F , decomposing v ∈ P +1 ( ) along (3) as v = grad + ( -) × curl with ∈ P +2 ( ) and ∈ P +1 ( ), inserting into the norm ± , (grad ) + +1 G, , (v ) , and using the triangle inequality, we infer, since v ∈ G +1 ( ) and , (grad

) = grad ( | ) ∈ G +1 ( ), ∈ T ℎ ∈ F ℎ -1 , (v ) -v 2 ∈ T ℎ ∈ F ℎ -1 +1 G, , (v ) -v 2 + ∈ T ℎ ∈ F ℎ -1 , v -grad 2 + ∈ T ℎ ∈ F ℎ -1 +1 G, , (v -grad ) 2 . (24) Using the 2 ( ; R 2 )-boundedness of +1
G, , a discrete trace inequality (cf., e.g., [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF]Lemma 1.32]), and Lemma 2 with = v , we infer

∈ T ℎ ∈ F ℎ -1 , (v ) -v 2 ∈ T ℎ curl v 2 + ∈ F ℎ -1 +1 G, , (v ) -v 2 . ( 25
)
Now, for the second factor on the right-hand side of ( 23), using that

| | × | ≤ | | |
, a continuous trace inequality (cf., e.g., [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF]Lemma 1.31]), the fact that ℎ ≤ diam(Ω) for all ∈ T ℎ , and concluding with [START_REF] Di Pietro | Fully discrete polynomial de Rham sequences of arbitrary degree on polygons and polyhedra[END_REF], one has

∈ T ℎ 2 + ∈ F ℎ | × 2 1 /2 1 (Ω;R 3 ) curl Ω . (26)
Plugging ( 25) and ( 26) into [START_REF] Du | A unified error analysis of Hybridizable Discontinuous Galerkin methods for the static Maxwell equations[END_REF], and recalling the definition (16) of the • X,ℎ -seminorm yields

I 2 v ℎ X,ℎ curl Ω ≤ v ℎ X,ℎ v ℎ Ω ,
where we have used the 2 (Ω; R 3 )-orthogonality of the decomposition [START_REF] Di Pietro | An arbitrary-order method for magnetostatics on polyhedral meshes based on a discrete de Rham sequence[END_REF] in the last bound. We conclude by combining ( 21), [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF], and this last estimate.

Remark 4 (Control of element unknowns). A direct proof of the fact that, for all

v ℎ ∈ • X +1 ℎ,0 , v ℎ X,ℎ = 0 implies v ℎ Ω =
0 can be obtained as follows. The volumetric term in [START_REF] Di Pietro | A third Strang lemma and an Aubin-Nitsche trick for schemes in fully discrete formulation[END_REF] first yields that, for all ∈ T ℎ , curl(v ℎ | ) = 0 in , meaning that v ℎ | = grad for some ∈ P +2 ( ) by Lemma 2. The boundary term in [START_REF] Di Pietro | A third Strang lemma and an Aubin-Nitsche trick for schemes in fully discrete formulation[END_REF] then yields the continuity of the tangential component of v ℎ at interfaces, as well as

×(v ℎ × ) = 0 on Ω. Hence, v ℎ ∈ 0 (curl; Ω) and curl v ℎ = 0 in Ω. Since v ℎ ∈ • X +1
ℎ,0 , we also have v ℎ ∈ (div 0 ; Ω) by the commutation property [START_REF] Chave | A three-dimensional Hybrid High-Order method for magnetostatics[END_REF]. Finally, by the continuous first Weber inequality (2), v ℎ Ω = 0. Notice that this result is weaker than the quantitative estimate of Theorem 3, as it does not give any information on how the constant W depends on the mesh at hand. Notice also that imposing v ℎ ∈ (div 0 ; Ω) as we do is actually not necessary. In view of the above analysis, it is sufficient and necessary to impose that v ℎ be orthogonal to the gradient of any function in P +2 (T ℎ ) ∩ 0 0 (Ω). This is the approach pursued in [START_REF] Chen | A priori and computable a posteriori error estimates for an HDG method for the coercive Maxwell equations[END_REF] and [START_REF] Chen | HDG and CG methods for the indefinite time-harmonic Maxwell's equations under minimal regularity[END_REF] on tetrahedral meshes.

Corollary 5 (Norm • X,ℎ ). The map • X,ℎ defines a norm on • X +1 ℎ,0 defined by [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF].

Proof. This is a direct consequence of Theorem 3 and of the definition [START_REF] Di Pietro | A third Strang lemma and an Aubin-Nitsche trick for schemes in fully discrete formulation[END_REF]. For all

v ℎ ∈ • X +1 ℎ,0 , if v ℎ X,ℎ = 0, then v ℎ = 0, i.e. v = 0 for all ∈ T ℎ . Then, for all ∈ F ℎ , +1 G, (v ) = v = 0, i.e. v = 0, whence v ℎ = 0 ℎ . Corollary 6 (Generalized discrete Weber inequality). Let d ℎ : Y +1 ℎ × Y +1 ℎ → R be a symmetric positive semi-definite bilinear form such that, for all ∈ 1 0 (Ω), letting φ ℎ := I +1 Y,ℎ ∈ Y +1 ℎ,0 , d ℎ φ ℎ , φ ℎ 1 /2 grad Ω . ( 27
)
Then, there is W > 0 independent of ℎ such that, for all

(v ℎ , r ℎ ) ∈ X +1 ℎ,0 ×Y +1 ℎ,0 satisfying -v ℎ , +1 ℎ q ℎ Ω + d ℎ r ℎ , q ℎ = 0 ∀q ℎ ∈ Y +1 ℎ,0 , (28)
one has

v ℎ Ω ≤ W v ℎ 2 X,ℎ + d ℎ (r ℎ , r ℎ ) 1 /2 . ( 29
)
Proof. We follow the steps of the proof of Theorem 3. If (v ℎ , r ℎ ) ∈ X +1 ℎ,0 ×Y +1 ℎ,0 satisfies (28), then [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF] becomes

I 1 = (v ℎ , grad ) Ω = v ℎ , +1 ℎ φ ℎ Ω = d ℎ r ℎ , φ ℎ .
By the Cauchy-Schwarz inequality and ( 27), we infer

I 1 d ℎ (r ℎ , r ℎ ) 1 /2 grad Ω ≤ d ℎ (r ℎ , r ℎ ) 1 /2 v ℎ Ω , (30)
where we have used the 2 (Ω; R 3 )-orthogonality of the decomposition [START_REF] Di Pietro | An arbitrary-order method for magnetostatics on polyhedral meshes based on a discrete de Rham sequence[END_REF] in the last bound. The rest of the proof is unchanged provided we substitute [START_REF] Nédélec | A new family of mixed finite elements in R 3[END_REF] to [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF].

This last corollary will be instrumental in the analysis of the HHO method of Section 3.2. Before closing this section, two additional remarks are in order.

Remark 7 (Topological assumptions on the domain). The first Weber inequality (2) is actually valid under the sole topological assumption that the boundary of Ω is connected, so that its second Betti number is zero. The same holds true for the discrete Weber inequalities of Theorem 3 and Corollary 6 (and, incidentally, this is also the case for discrete Weber inequalities on spaces with conforming unknowns, see [START_REF] Di Pietro | An arbitrary-order method for magnetostatics on polyhedral meshes based on a discrete de Rham sequence[END_REF]Theorem 19]). In other words, one does not need to assume, as we do, that Ω is simply-connected to prove these results. This last assumption is however necessary in the applicative Section 3 to have equivalence (i) between Problems (31) and (32) in field formulation, and (ii) between Problems (54) and (55) (in the class of potentials satisfying the Coulomb gauge) in vector potential formulation when div = 0.

Remark 8 (Star-shapedness assumption). We have assumed in Section 2.1 that the mesh elements are star-shaped. This assumption has been instrumental to prove Lemma 2 (and, in turn, Theorem 3 and Corollary 6), where it is used to infer a sign for the rightmost term in the second line of [START_REF] Cockburn | Bridging the Hybrid High-Order and Hybridizable Discontinuous Galerkin methods[END_REF].

According to [START_REF] Di Pietro | An arbitrary-order discrete de Rham complex on polyhedral meshes: Exactness, Poincaré inequalities, and consistency[END_REF]Lemma 46], this assumption is in fact not necessary for Lemma 2 to hold true, and all the results of this article actually (seamlessly) extend to the general case of meshes featuring non-necessarily star-shaped mesh elements (in that case, denotes any given interior point of such that contains a ball centered at of radius comparable to ℎ ). Nonetheless, we have preferred to keep this assumption in our analysis because it enables, as opposed to [START_REF] Di Pietro | An arbitrary-order discrete de Rham complex on polyhedral meshes: Exactness, Poincaré inequalities, and consistency[END_REF]Lemma 46], whose proof hinges on a transport argument, to obtain an explicit multiplicative constant in front of ℎ in [START_REF] Chen | A priori and computable a posteriori error estimates for an HDG method for the coercive Maxwell equations[END_REF].

Application to magnetostatics

In this section, we design and analyze HHO methods for the discretization of the magnetostatics equations. Their analysis leverages the discrete Weber inequality of Theorem 3 and its generalization pointed out in Corollary 6. For the sake of simplicity, we focus on homogeneous boundary conditions: the extension of HHO methods to non-homogeneous boundary conditions is standard (see, e.g., [17, Section 2.4]) and will be considered in the numerical example of Section 3.1.5. We recall that we work on regular (polyhedral) mesh sequences (M ℎ ) ℎ>0 in the sense of [17, Definition 1.9], which are characterized by the fact that the sequence of mesh regularity parameters is bounded from below by a strictly positive real number.

Field formulation

The model

The (first-order) field formulation of the magnetostatics problem consists in finding the magnetic field :

Ω → R 3 such that curl = in Ω, ( 31a
) div = 0 in Ω, (31b) ×( × ) = 0 on Ω, (31c)
where the current density : Ω → R 3 is such that div = 0 in Ω and • = 0 on Ω. We consider the following equivalent (cf. Remark 7) weak formulation of Problem [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF], originally introduced in [25, Eq. ( 52)] (see also [START_REF] Kanayama | Three-dimensional magnetostatic analysis using Nédélec's elements[END_REF]): Find ( ,

) ∈ 0 (curl; Ω) × 1 0 (Ω) such that ( , ) + ( , ) = ( , curl ) Ω ∀ ∈ 0 (curl; Ω), (32a)
( , ) = 0 ∀ ∈ 1 0 (Ω), (32b)
where the bilinear forms : (curl; Ω) × (curl; Ω) → R and :

(curl; Ω) × 1 (Ω) → R are given by ( , ) := (curl , curl ) Ω , ( ,
) := ( , grad ) Ω . (33)
The function : Ω → R is the Lagrange multiplier of the divergence-free constraint on the magnetic induction. Testing (32a) with = grad ∈ 0 (curl; Ω), it is readily inferred that = 0 in Ω. By the decomposition (1), Problem [START_REF] Perugia | Stabilized interior penalty methods for the time-harmonic Maxwell equations[END_REF] can then be equivalently rewritten: Find ∈ 0 (curl; Ω) ∩ (div 0 ; Ω) such that

( , ) = ( , curl ) Ω ∀ ∈ 0 (curl; Ω) ∩ (div 0 ; Ω),
whose well-posedness is a direct consequence of the first Weber inequality (2) and of the Lax-Milgram lemma.

The HHO method

We analyze in this section the HHO method for Problem [START_REF] Perugia | Stabilized interior penalty methods for the time-harmonic Maxwell equations[END_REF] we have briefly introduced in [START_REF] Chave | A three-dimensional Hybrid High-Order method for magnetostatics[END_REF]. This HHO method is based on the hybrid spaces introduced in Section 2.2 (the space X +1 ℎ for the magnetic field, and Y +1 ℎ for the Lagrange multiplier). We define the discrete bilinear forms

a ℎ : X +1 ℎ × X +1 ℎ → R, b ℎ : X +1 ℎ × Y +1 ℎ → R, and c ℎ : Y +1 ℎ × Y +1 ℎ → R such that a ℎ (w ℎ , v ℎ ) := curl ℎ w ℎ , curl ℎ v ℎ Ω + s ℎ (w ℎ , v ℎ ), (34a) b ℎ (w ℎ , q ℎ ) := w ℎ , +1 ℎ q ℎ Ω , (34b) c ℎ (r ℎ , q ℎ ) := (r ℎ , q ℎ ) Ω + ∈ T ℎ ∈ F ℎ (r , q ) , (34c)
where

+1 ℎ : Y +1 ℎ → P +1 (T ℎ )
is the gradient reconstruction operator introduced in Section 2.3, and s ℎ :

X +1 ℎ × X +1 ℎ → R is the stabilization bilinear form such that s ℎ (w ℎ , v ℎ ) := ∈ T ℎ ∈ F ℎ -1 +1 G, , (w ) -w , +1 G, , (v ) -v . ( 35
)
The HHO method for Problem (32) then reads:

Find (u ℎ , p ℎ ) ∈ X +1 ℎ,0 × Y +1 ℎ,0 such that a ℎ (u ℎ , v ℎ ) + b ℎ (v ℎ , p ℎ ) = ( , curl ℎ v ℎ ) Ω ∀v ℎ ∈ X +1 ℎ,0 , (36a)
-b ℎ (u ℎ , q ℎ ) + c ℎ (p ℎ , q ℎ ) = 0 ∀q ℎ ∈ Y +1 ℎ,0 . (36b)
Notice that, contrary to , the discrete Lagrange multiplier p ℎ is in general nonzero, as a consequence of the fact that the global discrete gradient +1 ℎ p ℎ is not irrotational. Some remarks are in order.

Remark 9 (The tetrahedral case). On matching tetrahedral meshes, according to Lemma 1,

+1 ℎ • Ω defines a norm on Y +1
ℎ,0 . Hence, in this case, one can consider a modified version of Problem (36) in which c ℎ is removed and for which stability (hence well-posedness) is preserved. This will be justified rigorously by Lemma 21 below (which essentially states an inf-sup condition for b ℎ ). Removing c ℎ ,

u ℎ ∈ • X +1
ℎ,0 holds true and, as a by-product of the commutation property [START_REF] Chave | A three-dimensional Hybrid High-Order method for magnetostatics[END_REF], u ℎ ∈ (div 0 ; Ω). Remark 10 (Improved stability). On general mesh families, and as opposed to the tetrahedral case, one has to add in (36b) a (consistent) positive semi-definite stabilization of the Lagrange multiplier to ensure the stability of the method. An example of such a stabilization is given by d ℎ defined by (59c).

Here, taking advantage of the fact that the continuous Lagrange multiplier is identically zero (as a consequence of (32a)), we choose to add the (consistent) positive-definite contribution c ℎ (p ℎ , q ℎ ) (remark that c ℎ defines a norm on Y +1 ℎ ). This strategy enables to improve the stability of the method without deteriorating its convergence properties. At the opposite, in the model of Section 3.2 below, the Lagrange multiplier may be nonzero and one cannot add the same contribution at the discrete level because it is not consistent anymore. One uses instead the semi-definite (consistent) bilinear form d ℎ of (59c) and can then only prove a weaker stability result (compare Lemma 13 and Lemma 22). To assess the effect of using c ℎ as a stabilization, we compare in Section 3.1.5 (see Figure 2) the numerical results obtained on a matching tetrahedral mesh family for Problem (36) with c ℎ and without c ℎ (which is possible according to Remark 9). Always better results (in absolute value) are obtained when using c ℎ .

Remark 11 (Curl reconstruction).

As opposed to what is done in HHO methods for second-order problems (see Section 3.2 below), we here take advantage of the fact that the problem is first-order to avoid (locally) reconstructing a discrete curl operator. Doing so, (i) it is possible to consider a smaller local space of face unknowns (that does not need to contain P ( )) for ≥ 1 (cf. [START_REF] Chave | A three-dimensional Hybrid High-Order method for magnetostatics[END_REF]Table 1]), and (ii) there is no need to solve a local problem on each mesh element (which may become, for a sequential implementation, rather costly in 3D for large polynomial degrees).

Letting Z +1 ℎ := X +1
ℎ ×Y +1 ℎ and Z +1 ℎ,0 := X +1 ℎ,0 ×Y +1 ℎ,0 , Problem (36) can be equivalently rewritten:

Find (u ℎ , p ℎ ) ∈ Z +1 ℎ,0 such that A ℎ (u ℎ , p ℎ ), (v ℎ , q ℎ ) = ( , curl ℎ v ℎ ) Ω ∀ (v ℎ , q ℎ ) ∈ Z +1 ℎ,0 , (37)
where the bilinear form A ℎ :

Z +1 ℎ × Z +1 ℎ → R is defined by A ℎ (w ℎ , r ℎ ), (v ℎ , q ℎ ) := a ℎ (w ℎ , v ℎ ) + b ℎ (v ℎ , r ℎ ) -b ℎ (w ℎ , q ℎ ) + c ℎ (r ℎ , q ℎ ). (38)
For future use, we also let

• Z +1 ℎ,0 := (w ℎ , r ℎ ) ∈ Z +1 ℎ,0 : -b ℎ (w ℎ , q ℎ ) + c ℎ (r ℎ , q ℎ ) = 0 ∀q ℎ ∈ Y +1 ℎ,0 = (w ℎ , r ℎ ) ∈ Z +1 ℎ,0 : A ℎ (w ℎ , r ℎ ), (0 ℎ , q ℎ ) = 0 ∀q ℎ ∈ Y +1 ℎ,0 .
(39)

Stability analysis

We recall that X +1 ℎ is equipped with the seminorm • X,ℎ defined by [START_REF] Di Pietro | A third Strang lemma and an Aubin-Nitsche trick for schemes in fully discrete formulation[END_REF], which is such that

• X,ℎ = a ℎ (•, •) 1 /2 . We equip Y +1 ℎ with the norm • Y,ℎ := c ℎ (•, •) 1 /2 , (40)
and Z +1 ℎ with the seminorm

(w ℎ , r ℎ ) Z,ℎ := w ℎ 2 X,ℎ + r ℎ 2 Y,ℎ 1 /2 . ( 41
) Lemma 12 (Norm • Z,ℎ ). The map • Z,ℎ defines a norm on • Z +1 ℎ,0 defined by (39).
Proof. The seminorm property being straightforward, it suffices to prove that, for all (w ℎ , r ℎ ) ∈

• Z +1 ℎ,0 , (w ℎ , r ℎ ) Z,ℎ = 0 implies (w ℎ , r ℎ ) = (0 ℎ , 0 ℎ ). Let then (w ℎ , r ℎ ) ∈ • Z +1
ℎ,0 satisfying (w ℎ , r ℎ ) Z,ℎ = 0. We infer that w ℎ X,ℎ = 0 and r ℎ Y,ℎ = 0. Since • Y,ℎ is a norm on Y +1 ℎ , the second relation directly implies that r ℎ = 0 ℎ . Now, owing to the definitions (39), (34b), and [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF], since (w

ℎ , r ℎ ) ∈ • Z +1 ℎ,0
and r ℎ = 0 ℎ , we have w ℎ ∈

• X +1 ℎ,0 . By Corollary 5, this implies w ℎ = 0 ℎ .

Lemma 13 (Well-posedness). For all

z ℎ ∈ Z +1 ℎ , A ℎ z ℎ , z ℎ = z ℎ 2 Z,ℎ . (42)
Hence, Problem (36) is well-posed, and the following a priori bound holds true:

(u ℎ , p ℎ ) Z,ℎ ≤ Ω . ( 43
)
Proof. The identity (42) is a direct consequence of (38) and (41) along with the definitions of • X,ℎ and • Y,ℎ . To prove well-posedness, since the linear system associated to Problem (36) is square, it is sufficient to prove injectivity. Assume that A ℎ (u

ℎ , p ℎ ), (v ℎ , q ℎ ) = 0 for all (v ℎ , q ℎ ) ∈ Z +1 ℎ,0 . Choosing (v ℎ , q ℎ ) = (0 ℎ , q ℎ
) and using (39), we first infer that (u

ℎ , p ℎ ) ∈ • Z +1 ℎ,0 . Choosing (v ℎ , q ℎ ) = (u ℎ , p ℎ
) and using (42), we then get

(u ℎ , p ℎ ) 2 Z,ℎ = A ℎ (u ℎ , p ℎ ), (u ℎ , p ℎ ) = 0,
which, by Lemma 12, eventually yields (u ℎ , p ℎ ) = (0 ℎ , 0 ℎ ). The a priori bound (43) directly follows from (42) with z ℎ = (u ℎ , p ℎ ), (37), the Cauchy-Schwarz inequality, and u ℎ X,ℎ ≤ (u ℎ , p ℎ ) Z,ℎ .

Error analysis

We recall that ( , ) ∈ 0 (curl; Ω) × 1 0 (Ω) denotes the unique solution to Problem [START_REF] Perugia | Stabilized interior penalty methods for the time-harmonic Maxwell equations[END_REF]. We assume from now on that possesses the additional regularity ∈ 1 (Ω; R 3 ), and we let u ℎ := I +1 X,ℎ ∈ X +1 ℎ,0 and p

ℎ := I +1 Y,ℎ ∈ Y +1 ℎ,0 .
In the spirit of [START_REF] Di Pietro | A third Strang lemma and an Aubin-Nitsche trick for schemes in fully discrete formulation[END_REF] (see also [17, Appendix A]), we estimate the errors

X +1 ℎ,0 e ℎ := u ℎ -u ℎ , Y +1 ℎ,0 ϵ ℎ := p ℎ -p ℎ , (44)
where

(u ℎ , p ℎ ) ∈ X +1 ℎ,0 × Y +1 ℎ,
) ∈ Z +1 ℎ,0 solve A ℎ (e ℎ , ϵ ℎ ), (v ℎ , q ℎ ) = l ℎ (v ℎ ) + m ℎ q ℎ ∀(v ℎ , q ℎ ) ∈ Z +1 ℎ,0 , (45)
where we have defined the consistency error linear forms

l ℎ (v ℎ ) := ( , curl ℎ v ℎ ) Ω -a ℎ ( u ℎ , v ℎ ), (46a)
m ℎ q ℎ := b ℎ ( u ℎ , q ℎ ). (46b)
Theorem 14 (Energy-error estimate). Assume that

∈ 0 (curl; Ω) ∩ 1 (Ω; R 3 ) ∩ +2 (T ℎ ; R 3 ).
Then, the following holds true, with (e ℎ , ϵ ℎ ) ∈ Z +1 ℎ,0 defined by (44):

(e ℎ , ϵ ℎ ) Z,ℎ ∈ T ℎ ℎ 2( +1) | | 2 +2 ( ;R 3 ) 1 /2 . ( 47
)
Proof. Since (e ℎ , ϵ ℎ ) ∈ Z +1 ℎ,0 , by (42) with z ℎ = (e ℎ , ϵ ℎ ) and (45), we infer

(e ℎ , ϵ ℎ ) Z,ℎ ≤ max (v ℎ ,q ℎ ) ∈Z +1 ℎ,0 , (v ℎ ,q ℎ ) Z,ℎ =1 l ℎ (v ℎ ) + m ℎ q ℎ . ( 48
)
Let us first focus on l ℎ (v ℎ ) for v ℎ ∈ X +1 ℎ,0 . Combining its definition (46a) with the fact that = curl in Ω, and the definition (6a) of I +1 X,ℎ , we infer

|l ℎ (v ℎ )| = curl -curl ℎ ( +1 P,ℎ ), curl ℎ v ℎ Ω -s ℎ ( u ℎ , v ℎ ) ≤ curl -curl ℎ ( +1 P,ℎ ) 2 Ω + s ℎ ( u ℎ , u ℎ ) 1 /2 v ℎ X,ℎ ,
where we have used the triangle/Cauchy-Schwarz inequalities and the definition (16) of v ℎ X,ℎ to pass to the second line. The quantity curlcurl ℎ ( +1 P,ℎ ) 2 Ω is estimated using the approximation properties of +1 P,ℎ (see, e.g., [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF]Theorem 1.45]). For the quantity s ℎ ( u ℎ , u ℎ ), recalling the definition (35) of s ℎ and using the 2 ( ; R 2 )-boundedness of +1 G, , we have

s ℎ ( u ℎ , u ℎ ) = ∈ T ℎ ∈ F ℎ -1 +1 G, , ( +1 P, ( | ) -) 2 ≤ ∈ T ℎ ∈ F ℎ -1 +1 P, ( | ) -2 ∈ T ℎ ℎ 2( +1) | | 2 +2 ( ;R 3 ) , (49)
where, for all ∈ T ℎ , we have used the approximation properties of +1 P, on the faces of . Gathering the different estimates, we get

|l ℎ (v ℎ )| ∈ T ℎ ℎ 2( +1) | | 2 +2 ( ;R 3 ) 1 /2 v ℎ X,ℎ . (50)
Let us now focus on m ℎ q ℎ for q ℎ ∈ Y +1 ℎ,0 . Starting from (46b), performing an element-byelement integration by parts in [START_REF] Brezzi | Two families of mixed finite elements for second order elliptic problems[END_REF], and using that grad q ∈ G -1 ( ) ⊂ P +1 ( ), we infer

m ℎ q ℎ = ∈ T ℎ (grad q , ) + ∈ F +1 P, ( | ) | • , q -q | = ∈ T ℎ ∈ F ( +1 P, ( | ) -) | • , q -q | ,
where the last identity follows from another element-by-element integration by parts, and from the fact that div = 0 in Ω, and that ∈ 1 (Ω; R 3 ) along with q = 0 for all ∈ F b ℎ . By the triangle and Cauchy-Schwarz inequalities, one then gets

|m ℎ q ℎ | ≤ ∈ T ℎ ∈ F ℎ -1 +1 P, ( | ) -2 1 /2 × ∈ T ℎ ∈ F ℎ q -q | 2 1 /2 . (51)
Using, for all ∈ T ℎ , the approximation properties of +1 P, on the faces of for the first factor on the right-hand side, and the triangle inequality along with a discrete trace inequality (see, e.g., [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF]Lemma 1.32]) for the second factor, we infer

|m ℎ q ℎ | ∈ T ℎ ℎ 2( +1) | | 2 +2 ( ;R 3 ) 1 /2 q ℎ Y,ℎ . (52)
Plugging (50) and ( 52) into (48) for (v ℎ , q ℎ

) such that (v ℎ , q ℎ ) Z,ℎ = 1 finally yields (47).

Numerical results

Let the domain Ω be the unit cube (0, 1) 3 . We consider Problem [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF] with exact solution

( 1 , 2 , 3 ) := cos( 2 ) cos( 3 ) cos( 1 ) cos( 3 ) cos( 1 ) cos( 2 ) . ( 53
)
Clearly, the function satisfies (31b). The source is set according to (31a), and the zero tangential boundary condition (31c) is replaced by the non-homogeneous boundary condition stemming from (53). We solve the discrete Problem (36) with amended right-hand side accounting for the nonhomogeneous boundary condition on two refined mesh sequences, of respectively cubic and regular tetrahedral meshes. For each problem, the element unknowns for both the magnetic field and the Lagrange multiplier are locally eliminated using a Schur complement technique. This step is fully parallelizable. The resulting (condensed) global linear system is solved using the SparseLU direct solver of the Eigen library, on an Intel Xeon E-2176M 2.70GHz×12 with 16GB of RAM (and up to 150GB of swap). For ∈ {0, 1, 2}, we depict on Figures 1 and2, respectively for the cubic and (regular) tetrahedral mesh families, the relative energy-error u ℎ -I +1 X,ℎ X,ℎ / I +1 X,ℎ X,ℎ (top row) and 2 -error u ℎ -+1 P,ℎ Ω / +1 P,ℎ Ω (bottom row) as functions of (i) the meshsize (left column), (ii) the solution time in seconds, i.e. the time needed to solve the (condensed) global linear system (center column), and (iii) the number of (interface) degrees of freedom (DoF) (right column). For the two mesh families, we obtain, as predicted by Theorem 14, an energy-error convergence rate of order + 1. We also observe a convergence rate of order + 2 for the 2 -error. We remark that, whenever the solution is smooth enough (at least locally), raising the polynomial degree is computationally much more efficient than refining the mesh to increase the accuracy. Following Remark 9, we also solve on the (matching) tetrahedral mesh family a modifed version of Problem (36) in which c ℎ is removed, and we display on Figure 2 the results in dashed lines. Also in this case, we obtain the predicted energy-error convergence rate of order + 1, and observe a convergence rate of order + 2 for the 2 -error. We remark that the results using c ℎ are always better (in absolute value) than those obtained without using it.

= 0 = 1 = 2 10 -1.

Vector potential formulation

The model

The (second-order) vector potential formulation of the magnetostatics problem consists, in its generalized form, in finding the magnetic vector potential : Ω → R 3 and the Lagrange multiplier : Ω → R such that

curl curl + grad = in Ω, ( 54a
) div = 0 in Ω, (54b) ×( × ) = 0 on Ω, (54c) = 0 on Ω, (54d)
where the current density : Ω → R 3 is no longer assumed to be divergence-free, whence the introduction of the Lagrange multiplier term in (54a). When div = 0 in Ω, = 0 in Ω and, letting := curl , Problem (54) is then equivalent, in the class of vector potentials satisfying the Coulomb gauge, to the following problem (cf. Remark 7):

curl = in Ω, div = 0 in Ω, • = 0 on Ω. ( 55
)
The function is then the vector potential associated to the magnetic induction . Assuming that ∈ 2 (Ω; R 3 ), we consider the following equivalent weak formulation of Problem (54

): Find ( , ) ∈ 0 (curl; Ω) × 1 0 (Ω) such that ( , ) + ( , ) = ( , ) Ω ∀ ∈ 0 (curl; Ω), (56a)
( , ) = 0 ∀ ∈ 1 0 (Ω), (56b)
where the bilinear forms : (curl; Ω) × (curl; Ω) → R and : (curl; Ω) × 1 (Ω) → R are defined in [START_REF] Wang | A weak Galerkin finite element method for second-order elliptic problems[END_REF]. Using the decomposition (1), Problem (56) can be equivalently rewritten under the following fully decoupled form: Find ∈ 0 (curl; Ω) ∩ (div 0 ; Ω) and ∈ 1 0 (Ω) such that ( , ) = ( , ) Ω ∀ ∈ 0 (curl; Ω) ∩ (div 0 ; Ω), (grad , ) = ( , grad ) Ω ∀ ∈ 1 0 (Ω), whose well-posedness directly follows from the first Weber inequality (2) and from the Lax-Milgram lemma.

The HHO method

We consider the hybrid spaces introduced in Section 2.2 (the space X +1 ℎ for the magnetic vector potential, and Y +1 ℎ for the Lagrange multiplier), up to a slight modification of the space X +1 ℎ defined in (4a). To this end, we introduce, for any ∈ N and any ∈ F ℎ , the space

P ♭ ( ) := P -1 ( ) ⊕ grad P +1 ( ) , (57)
with P +1 ( ) denoting the space of homogeneous scalar-valued polynomials of total degree + 1 on , and the convention that P -1 ( ) := {0}. Consistently with our notation so far, we let P ♭ , denote the 2 ( ; R 2 )-orthogonal projector onto P ♭ ( ). With this new space at hand, we define

X +1 ♭,ℎ := v ℎ = (v ) ∈ T ℎ , (v ) ∈ F ℎ : v ∈ P +1 ( ) ∀ ∈ T ℎ v ∈ P +1 ♭ ( ) ∀ ∈ F ℎ , ( 58
)
which is from now on meant to replace the space X +1 ℎ . The space Y +1 ℎ keeps the same definition (4b). We also introduce the spaces X +1 ♭,ℎ,0 and • X +1 ♭,ℎ,0 , that are respectively obtained through definitions ( 5) and ( 17), up to the replacement therein of X +1 ℎ by X +1 ♭,ℎ , and of X +1 ℎ,0 by X +1 ♭,ℎ,0 . In turn, the interpolator I +1 X,♭,ℎ : 1 (Ω; R 3 ) → X +1 ♭,ℎ is defined as in formula (6a), up to the replacement of the projector +1 G, by +1 P ♭ , . Here, and as opposed to Section 3.1 (cf. Remark 11), because of the fact that we will have to reconstruct a discrete curl operator, we need to consider a space for the vectorial face unknowns that contains P ( ) (this is indeed needed to prove optimal approximation properties for the curl reconstruction operator). Since, for stability reasons, the space for face unknowns must also contain G +1 ( ), we consider the sum of these two spaces, which rewrites as the direct sum (57). Notice that P +1 ♭ ( ) is strictly sandwiched between P ( ) and P +1 ( ). Remark 15 (Validity of the results of Section 2.4). For all ∈ F ℎ , we have G +1 ( ) ⊆ P +1 ♭ ( ), and it can be checked that, up to the replacement of the 2 ( ; R 2 )-orthogonal projector +1 G, onto G +1 ( ) by the projector +1 P ♭ , onto P +1 ♭ ( ), all the results in Section 2.4 remain valid when X +1 ℎ is replaced by X +1 ♭,ℎ as defined in (58), including the discrete Weber inequality of Theorem 3 and its generalization of Corollary 6 (observe, in particular, that the crucial estimates (24)-( 25) still hold true under these changes).

We define the discrete bilinear forms a ℎ :

X +1 ♭,ℎ × X +1 ♭,ℎ → R, b ℎ : X +1 ♭,ℎ × Y +1 ℎ → R, and d ℎ : Y +1 ℎ × Y +1 ℎ → R such that a ℎ (w ℎ , v ℎ ) := ℎ w ℎ , ℎ v ℎ Ω + s ℎ (w ℎ , v ℎ ), (59a) b
ℎ (w ℎ , q ℎ ) := w ℎ , +1 ℎ q ℎ Ω , (59b)
d ℎ (r ℎ , q ℎ ) := ∈ T ℎ ∈ F ℎ r -r | , q -q | , (59c)
where +1 ℎ : Y +1 ℎ → P +1 (T ℎ ) is the gradient reconstruction operator introduced in Section 2.3, and s ℎ :

X +1 ♭,ℎ × X +1 ♭,ℎ → R is the stabilization bilinear form such that s ℎ (w ℎ , v ℎ ) := ∈ T ℎ ∈ F ℎ -1 +1 P ♭ , , (w )
-w , +1 P ♭ , , (v )
-v . (60)
In (59a), ℎ :

X +1 ♭,ℎ → R (T ℎ ) (with R (T ℎ ) defined in Section 2.
2) is the global discrete curl reconstruction operator such that its local restriction : X +1 ♭, → R ( ) to any ∈ T ℎ solves the following well-posed problem: For all v ∈ X +1 ♭, ,

v , = (v , curl ) + ∈ F v , , ( × ) ∀ ∈ R ( ). (61)
With this definition at hand, one can prove the following commutation property.

Lemma 16 (Commutation property). For all ∈ 1 (Ω; R 3 ), we have

ℎ (I +1 X,♭,ℎ ) = R,ℎ (curl ), (62)
where we remind the reader that R,ℎ is the 2 (Ω; R 3 )-orthogonal projector onto R (T ℎ ).

Proof. For ∈ 1 (Ω; R 3 ), let v ℎ := I +1 X,♭,ℎ . Then, for any ∈ T ℎ , v = +1 P, ( | ), +1 P ♭ , , ( ) ∈ F .
Plugging v into (61), and using that curl ∈ P -1 ( ) ⊂ P +1 ( ) and that , ( × ) ∈ P ( ) ⊂ P +1 ♭ ( ) for all ∈ F to remove the projectors, one gets, for all ∈ R ( ),

v , = ( , curl ) + ∈ F | , | × . ( 63
)
Integrating by parts the right-hand side of (63), we readily infer (62).

Remark 17 (Variant on ℎ

). An alternative choice consists in reconstructing the discrete curl in P (T ℎ ), which requires to solve larger local problems for ≥ 1 (for example, the space R 2 ( ) has dimension 26, whereas P 2 ( ) has dimension 30). In this case, the commutation property (62) reads ℎ (I +1 X,♭,ℎ ) = P,ℎ (curl ). This is the approach pursued in [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF] and [START_REF] Chen | A superconvergent HDG method for the Maxwell equations[END_REF]. The numerical tests we have performed (not reported here) indicate that, interestingly, reconstructing the discrete curl in P (T ℎ ) instead of R (T ℎ ), besides being computationally more expensive, sometimes deteriorates the accuracy of the approximation.

The HHO method for Problem (56) reads:

Find (u ℎ , p ℎ ) ∈ X +1 ♭,ℎ,0 × Y +1 ℎ,0 such that a ℎ (u ℎ , v ℎ ) + b ℎ (v ℎ , p ℎ ) = ( , v ℎ ) Ω ∀v ℎ ∈ X +1 ♭,ℎ,0 , (64a)
-b ℎ (u ℎ , q ℎ ) + d ℎ (p ℎ , q ℎ ) = 0 ∀q ℎ ∈ Y +1 ℎ,0 . (64b)
Notice that, as opposed to and in Problem (56), one cannot efficiently solve for u ℎ and p ℎ independently in Problem (64) as the curl-grad orthogonality is lost at the discrete level.

Remark 18 (Divergence-free current density). At the continuous level, remark that div = 0 in Ω implies = 0 in Ω. At the discrete level, when is divergence-free, one can hence mimick the strategy advocated in Section 3.1 and replace in (64b) the bilinear form d ℎ by the (consistent) bilinear form c ℎ given by (34c). Doing so improves the stability of the method without deteriorating its convergence properties (cf. Remark 10).

Letting Z +1 ♭,ℎ := X +1 ♭,ℎ × Y +1 ℎ and Z +1 ♭,ℎ,0 := X +1 ♭,ℎ,0 × Y +1
ℎ,0 , we notice that Problem (64) can be equivalently rewritten:

Find (u ℎ , p ℎ ) ∈ Z +1 ♭,ℎ,0 such that A ℎ (u ℎ , p ℎ ), (v ℎ , q ℎ ) = ( , v ℎ ) Ω ∀ (v ℎ , q ℎ ) ∈ Z +1 ♭,ℎ,0 , (65)
where the bilinear form A ℎ :

Z +1 ♭,ℎ × Z +1 ♭,ℎ → R is defined by A ℎ (w ℎ , r ℎ ), (v ℎ , q ℎ ) := a ℎ (w ℎ , v ℎ ) + b ℎ (v ℎ , r ℎ ) -b ℎ (w ℎ , q ℎ ) + d ℎ (r ℎ , q ℎ ). (66)
We also let, in analogy with (39),

• Z +1 ♭,ℎ,0 := (w ℎ , r ℎ ) ∈ Z +1 ♭,ℎ,0 : -b ℎ (w ℎ , q ℎ ) + d ℎ (r ℎ , q ℎ ) = 0 ∀q ℎ ∈ Y +1 ℎ,0 = (w ℎ , r ℎ ) ∈ Z +1 ♭,ℎ,0 : A ℎ (w ℎ , r ℎ ), (0 ℎ , q ℎ ) = 0 ∀q ℎ ∈ Y +1 ℎ,0 . (67)

Stability analysis

We equip the spaces X +1 ♭,ℎ and Y +1 ℎ with the seminorms

w ℎ X,♭,ℎ := curl ℎ w ℎ 2 Ω + s ℎ (w ℎ , w ℎ ) 1 /2 , (68a)
r ℎ Y,♭,ℎ := ∈ T ℎ ℎ 2 grad r 2 + d ℎ (r ℎ , r ℎ ) 1 /2 . ( 68b
)
One can easily verify that • Y,♭,ℎ defines a norm on Y +1 ℎ,0 . We now equip Z +1 ♭,ℎ with the seminorm

(w ℎ , r ℎ ) Z,♭,ℎ := w ℎ 2 X,♭,ℎ + r ℎ 2 Y,♭,ℎ 1 /2 . ( 69
) Lemma 19 (Norm • Z,♭,ℎ ). The map • Z,♭,ℎ defines a norm on • Z +1
♭,ℎ,0 defined by (67). Proof. The seminorm property being straightforward, we only need to prove that, for all couples

(w ℎ , r ℎ ) ∈ • Z +1 ♭,ℎ,0 , (w ℎ , r ℎ ) Z,♭,ℎ = 0 implies (w ℎ , r ℎ ) = (0 ℎ , 0 ℎ ). Let then (w ℎ , r ℎ ) ∈ • Z +1 ♭,ℎ,0 be such that (w ℎ , r ℎ ) Z,♭,ℎ = 0. We infer that w ℎ X,♭,ℎ = 0 and r ℎ Y,♭,ℎ = 0. Since • Y,♭,ℎ is a norm on Y +1
ℎ,0 , we directly get from the second relation that r ℎ = 0 ℎ . Now, owing to the definitions (67) of

• Z +1
♭,ℎ,0 , (59b) of b ℎ , and to the fact that • X +1 ♭,ℎ,0 is defined as in [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF] with X +1 ℎ,0 replaced by X +1 ♭,ℎ,0 , we infer from (w ℎ , r ℎ ) ∈ ♭,ℎ,0 , hence w ℎ = 0 ℎ , which concludes the proof. We now state some preliminary results for the stability analysis of Problem (64).

Lemma 20 (Equivalences of seminorms). The following holds true:

w ℎ 2 X,♭,ℎ a ℎ (w ℎ , w ℎ ) w ℎ 2 X,♭,ℎ ∀w ℎ ∈ X +1 ♭,ℎ , (70a)
r ℎ 2 Y,♭,ℎ ∈ T ℎ ℎ 2 +1 r 2 + d ℎ (r ℎ , r ℎ ) r ℎ 2 Y,♭,ℎ ∀r ℎ ∈ Y +1 ℎ . (70b)
Proof. Let us prove (70a). Let w ℎ ∈ X +1 ♭,ℎ , and ∈ T ℎ . By the definition (61) of , testing with = curl w ∈ R ( ), integrating by parts, and using the fact that , curl w × ∈ P ( ) ⊂ P +1

♭ ( ), we infer curl w 2 = w , curl w

+ ∈ F +1 P ♭ , , (w )
w 2 = curl w , w - ∈ F +1 P ♭ , , (w )
w , , ( w × ) ,

and we conclude by the same arguments. The proof of (70b) is similar and is omitted for brevity.

Lemma 21 (Control of +1 ℎ ). For all (w ℎ , r ℎ ) ∈ Z +1 ♭,ℎ , there exists v ★ ℎ ∈ X +1 ♭,ℎ,0 satisfying

v ★ ℎ 2 Ω + v ★ ℎ 2 X,♭,ℎ ∈ T ℎ ℎ 2 +1 r 2 ,
such that the following holds true:

∈ T ℎ ℎ 2 +1 r 2 A ℎ (w ℎ , r ℎ ), (v ★ ℎ , 0 ℎ ) + w ℎ 2 X,♭,ℎ . (71)
Proof. Let (w ℎ , r ℎ ) ∈ Z +1 ♭,ℎ . We define v ★ ℎ ∈ X +1 ♭,ℎ,0 such that

v ★ := ℎ 2 +1 r ∀ ∈ T ℎ , v ★ := 0 ∀ ∈ F ℎ . (72)
We immediately verify, since ℎ ≤ diam(Ω) for all ∈ T ℎ , that the following holds true:

v ★ ℎ 2 Ω ∈ T ℎ ℎ 2 +1 r 2 .
From the definitions (72) of v ★ ℎ , (59b) of b ℎ , and (66) of A ℎ , we infer

∈ T ℎ ℎ 2 +1 r 2 = b ℎ (v ★ ℎ , r ℎ ) = A ℎ (w ℎ , r ℎ ), (v ★ ℎ , 0 ℎ ) -a ℎ (w ℎ , v ★ ℎ ).
The Cauchy-Schwarz inequality followed by the second inequality in (70a) then yields

∈ T ℎ ℎ 2 +1 r 2 A ℎ (w ℎ , r ℎ ), (v ★ ℎ , 0 ℎ ) + w ℎ X,♭,ℎ v ★ ℎ X,♭,ℎ . (73)
Using the definitions (72) of v ★ ℎ and (68a) of • X,♭,ℎ , we get

v ★ ℎ 2 X,♭,ℎ = ∈ T ℎ curl v ★ 2 + ∈ F ℎ -1 +1 P ♭ , , (v ★ ) 2 ∈ T ℎ ℎ -2 v ★ 2 = ∈ T ℎ ℎ 2 +1 r 2 , (74)
where we have used the 2 ( ; R 2 )-boundedness of +1 P ♭ , , as well as an inverse inequality together with a discrete trace inequality (see, e.g., [START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF]Lemmas 1.28 and 1.32]). Starting from (73), and using (74) combined with a Young inequality for the last term in the right-hand side eventually yields the expected result (71).

We are now in position to show well-posedness for Problem (64).

Lemma 22 (Well-posedness). For all z

ℎ ∈ Z +1 ♭,ℎ , there exists v ★ ℎ ∈ X +1 ♭,ℎ,0 satisfying v ★ ℎ Ω + v ★ ℎ X,♭,ℎ z ℎ Z,♭,ℎ and such that A ℎ z ℎ , z ℎ + A ℎ z ℎ , (v ★ ℎ , 0 ℎ ) z ℎ 2 Z,♭,ℎ . (75)
Hence, Problem (64) is well-posed, and the following a priori bound holds true:

(u ℎ , p ℎ ) Z,♭,ℎ Ω . (76)
Proof. Let z ℎ = (w ℎ , r ℎ ) ∈ Z +1 ♭,ℎ . By the first inequality in (70a) and ( 66), one has

w ℎ 2 X,♭,ℎ + d ℎ (r ℎ , r ℎ ) a ℎ (w ℎ , w ℎ ) + d ℎ (r ℎ , r ℎ ) = A ℎ z ℎ , z ℎ . (77)
By Lemma 21 combined with (77), one also has

∈ T ℎ ℎ 2 +1 r 2 A ℎ z ℎ , (v ★ ℎ , 0 ℎ ) + A ℎ z ℎ , z ℎ , (78)
for some v ★ ℎ ∈ X +1 ♭,ℎ,0 such that v ★ ℎ Ω + v ★ ℎ X,♭,ℎ r ℎ Y,♭,ℎ ≤ z ℎ Z,♭,ℎ
, where we have used the second inequality in (70b) and the definition (69) of • Z,♭,ℎ . Summing (77) and (78), and using the first inequality in (70b), we infer (75).

To prove well-posedness, since the linear system associated to Problem (64) is square, it is sufficient to prove injectivity. Assume that A ℎ (u

ℎ , p ℎ ), (v ℎ , q ℎ ) = 0 for all (v ℎ , q ℎ ) ∈ Z +1 ♭,ℎ,0 . Taking (v ℎ , q ℎ ) = (0 ℎ , q ℎ
) and using (67), we first infer that (u

ℎ , p ℎ ) ∈ • Z +1 ♭,ℎ,0 . Taking (v ℎ , q ℎ ) = (u ℎ , p ℎ ) and (v ℎ , q ℎ ) = (v ★ ℎ , 0 ℎ
) and using (75), we then get

(u ℎ , p ℎ ) 2 Z,♭,ℎ A ℎ (u ℎ , p ℎ ), (u ℎ , p ℎ ) + A ℎ (u ℎ , p ℎ ), (v ★ ℎ , 0 ℎ ) = 0,
which, by Lemma 19, eventually yields (u ℎ , p ℎ ) = (0 ℎ , 0 ℎ ). To prove the a priori bound (76), we take z ℎ = (u ℎ , p ℎ ) in (75) and we use (65). We get, by the Cauchy-Schwarz inequality,

(u ℎ , p ℎ ) 2 Z,♭,ℎ ( , u ℎ ) Ω + ( , v ★ ℎ ) Ω ≤ Ω u ℎ Ω + v ★ ℎ Ω .
The conclusion follows from the fact that v ★

Error analysis

We recall that ( , ) ∈ 0 (curl; Ω) × 1 0 (Ω) denotes the unique solution to Problem (56). We assume from now on that possesses the additional regularity ∈ 1 (Ω; R 3 ), and we let u ℎ := I +1 X,♭,ℎ ∈ X +1 ♭,ℎ,0 and p ℎ := I +1 Y,ℎ ∈ Y +1 ℎ,0 . We define the errors

X +1 ♭,ℎ,0 e ℎ := u ℎ -u ℎ , Y +1 ℎ,0 ϵ ℎ := p ℎ -p ℎ , (79)
where

(u ℎ , p ℎ ) ∈ X +1 ♭,ℎ,0 × Y +1 ℎ,
) ∈ Z +1 ♭,ℎ,0 solve A ℎ (e ℎ , ϵ ℎ ), (v ℎ , q ℎ ) = l ℎ (v ℎ ) + m ℎ q ℎ ∀ (v ℎ , q ℎ ) ∈ Z +1 ♭,ℎ,0 , (80)
where we have defined the consistency error linear forms

l ℎ (v ℎ ) := ( , v ℎ ) Ω -a ℎ ( u ℎ , v ℎ ) -b ℎ (v ℎ , p ℎ ), (81a)
m ℎ q ℎ := b ℎ ( u ℎ , q ℎ ) -d ℎ ( p ℎ , q ℎ ). (81b)
Theorem 23 (Energy-error estimate). Assume that

∈ 0 (curl; Ω) ∩ 1 (Ω; R 3 ) ∩ +2 (T ℎ ; R 3 ), ∈ 1 0 (Ω) ∩ +1 (T ℎ ).
Then, the following holds true, with (e ℎ , ϵ ℎ ) ∈ Z +1 ♭,ℎ,0 defined by (79):

(e ℎ , ϵ ℎ ) Z,♭,ℎ ∈ T ℎ ℎ 2( +1) | | 2 +2 ( ;R 3 ) + | | 2 +1 ( ) 1 /2 . ( 82
)
Proof. Since (e ℎ , ϵ ℎ ) ∈ Z +1 ♭,ℎ,0 , by (75) with

z ℎ = (e ℎ , ϵ ℎ ) (notice that (v ★ ℎ , 0 ℎ ) ∈ Z +1 ♭,ℎ,0 and (v ★ ℎ , 0 ℎ ) Z,♭,ℎ = v ★ ℎ X,♭,ℎ z ℎ Z,♭,ℎ
) and (80), we infer

(e ℎ , ϵ ℎ ) Z,♭,ℎ max (v ℎ ,q ℎ ) ∈Z +1 ♭,ℎ,0 , (v ℎ ,q ℎ ) Z,♭,ℎ =1 l ℎ (v ℎ ) + m ℎ q ℎ . (83)
Let us first focus on l ℎ (v ℎ ) for v ℎ ∈ X +1 ♭,ℎ,0 . By (81a), the fact that = curl(curl ) + grad in Ω, and element-by-element integration by parts, we infer

l ℎ (v ℎ ) = ∈ T ℎ curl , curl v - ∈ F curl | × , v | + grad , v ℎ Ω -a ℎ ( u ℎ , v ℎ ) -b ℎ (v ℎ , p ℎ ) = ∈ T ℎ curl , curl v - ∈ F , curl × , , (v ) -v -a ℎ ( u ℎ , v ℎ ), (84)
where we have used the fact that the tangential component of curl is continuous across interfaces (as a consequence of the fact that curl ∈ (curl; Ω) ∩ 1 (T ℎ ; R 3 )) along with v = 0 for all ∈ F b ℎ to insert v into the boundary term, together with the fact that grad ,

v ℎ Ω = b ℎ (v ℎ , p ℎ )
as a consequence of the commutation property [START_REF] Chave | A three-dimensional Hybrid High-Order method for magnetostatics[END_REF]. Using the definitions (59a) of a ℎ and (61) of for ∈ T ℎ (testing with = u ∈ R ( )), and integrating by parts, we have

a ℎ ( u ℎ , v ℎ ) = ∈ T ℎ u , curl v - ∈ F , ( u × ), , (v ) -v + s ℎ ( u ℎ , v ℎ ). (85)
Since, by Lemma 16, u = R, curl | for all ∈ T ℎ , a combination of (84) and ( 85) yields (recall that v ∈ P +1 ( ))

l ℎ (v ℎ ) = ∈ T ℎ ∈ F , ( R, (curl | ) -curl )× , , (v ) -v -s ℎ ( u ℎ , v ℎ ).
Applying the triangle and Cauchy-Schwarz inequalities, we get

|l ℎ (v ℎ )| ≤ ∈ T ℎ ∈ F R, (curl | ) -curl , (v ) -v + s ℎ ( u ℎ , u ℎ ) 1 /2 s ℎ (v ℎ , v ℎ ) 1 /2 . (
ℎ 1 /2 R, (curl | ) -curl R, (curl | ) -curl + curl -curl +1 P, ( | ) + ℎ +1 | | +2 ( ;R 3 ) ,
where we have used yet another triangle inequality to insert curl . The second term on the right-hand side is readily estimated using again the approximation properties of +1 P, . As far as the first term is concerned, we have

R, (curl | ) -curl = min ∈P +1 ( ) curl -curl , which finally yields ℎ 1 /2 R, (curl | ) -curl ℎ +1 | | +2 ( ;R 3 ) .
Plugging this last estimate into (86), applying a discrete Cauchy-Schwarz inequality, and using [START_REF] Kikuchi | Mixed formulations for finite element analysis of magnetostatic and electrostatic problems[END_REF] as well as (49) (with +1 P ♭ , instead of +1 G, ), we infer

|l ℎ (v ℎ )| ∈ T ℎ ℎ 2( +1) | | 2 +2 ( ;R 3 ) 1 /2 v ℎ X,♭,ℎ . (87)
Let us now focus on m ℎ q ℎ for q ℎ ∈ Y +1 ℎ,0 . Recalling (51) and (46b), one can readily infer

|b ℎ ( u ℎ , q ℎ )| ∈ T ℎ ℎ 2( +1) | | 2 +2 ( ;R 3 ) 1 /2 d ℎ (q ℎ , q ℎ ) 1 /2 . (
)| ≤ d ℎ ( p ℎ , p ℎ ) 1 /2 d ℎ (q ℎ , q ℎ ) 1 /2 ≤ ∈ T ℎ ∈ F ℎ -P, ( | ) 2 1 /2 d ℎ (q ℎ , q ℎ ) 1 /2 .
By the approximation properties of P, on mesh faces (see, e.g., [17, Theorem 1.45]), we get

|d ℎ ( p ℎ , q ℎ )| ∈ T ℎ ℎ 2( +1) | | 2 +1 ( ) 1 /2 d ℎ (q ℎ , q ℎ ) 1 /2 . (89)
Gathering ( 88)-( 89), and recalling the definition (81b) of m ℎ q ℎ , finally yields

|m ℎ q ℎ | ∈ T ℎ ℎ 2( +1) | | 2 +2 ( ;R 3 ) + | | 2 +1 ( ) 1 /2 q ℎ Y,♭,ℎ . (90)
Plugging ( 87) and ( 90) into (83) for (v ℎ , q ℎ

) such that (v ℎ , q ℎ ) Z,♭,ℎ = 1 finally yields (82).

Remark 24 (The tetrahedral case). On matching tetrahedral meshes, according to Lemma 1,

+1 ℎ • Ω defines a norm on Y +1
ℎ,0 . Hence, in this case, and as already pointed out in [START_REF] Chen | A superconvergent HDG method for the Maxwell equations[END_REF], one can consider a modified version of Problem (64) in which d ℎ is removed and for which stability (hence wellposedness) is preserved. This is a direct consequence of Lemma 21, which essentially states an inf-sup condition for b ℎ provided

+1 ℎ • Ω is a norm on Y +1 ℎ,0 . Removing d ℎ , u ℎ ∈ • X +1
♭,ℎ,0 holds true and, as a by-product of the commutation property [START_REF] Chave | A three-dimensional Hybrid High-Order method for magnetostatics[END_REF], u ℎ ∈ (div 0 ; Ω). Furthermore, a close inspection of the proof of Theorem 23 shows that, in this case, one ends up with an energy-error estimate that is free of Lagrange multiplier contribution. This allows one to reproduce at the discrete level the following structure of Problem (54). When the current density is given by the gradient of some function ∈ 1 0 (Ω), then = 0 and = . At the discrete level, when = grad , one then gets u ℎ = 0 ℎ and p ℎ = I +1 Y,ℎ (note that this can be observed in practice up to machine precision only if the computation of the right-hand side is also performed up to machine precision).

Numerical results

Let the domain Ω be the unit cube (0, 1) 3 . We consider Problem (54) with exact solution ( 1 , 2 , 3 ) := sin( 2 ) sin( 3) sin( 1) sin( 3) sin( 1) sin( 2) , ( 1 , 2 , 3 ) := sin( 1 ) sin( 2) sin( 3 ).

Clearly, the magnetic vector potential and the Lagrange multiplier satisfy (54b), (54c), and (54d). The current density is set according to (54a).

As in Section 3.1.5, we solve the discrete Problem (64) on two refined mesh sequences, of respectively cubic and regular tetrahedral meshes. For each problem, the element unknowns for both the magnetic vector potential and the Lagrange multiplier are locally eliminated using a Schur complement technique. This step is fully parallelizable. The resulting (condensed) global linear system is solved using the SparseLU direct solver of the Eigen library on the same architecture as in Section 3.1.5. For ∈ {0, 1, 2}, we depict on Figures 3 and4, respectively for the cubic and (regular) tetrahedral mesh families, the relative energy-error u ℎ -I +1 (bottom row) as functions of (i) the meshsize (left column), (ii) the solution time in seconds, i.e. the time needed to solve the (condensed) global linear system (center column), and (iii) the number of (interface) DoF (right column). For the two mesh families, we obtain, as predicted by Theorem 23, a convergence rate of the energy-error on the magnetic vector potential and of the 2 -like-error on the Lagrange multiplier of order + 1. We also observe a convergence rate of order + 2 for the 2 -error on the magnetic vector potential. Following Remark 24, we also solve on the (matching) tetrahedral mesh family a modified version of Problem (64) in which d ℎ is removed, and we display on , which remains meaningful for = 0. Also in this case, we obtain the predicted energy-error on the magnetic vector potential and 2 -like-error on the Lagrange multiplier of order + 1, and we observe a convergence rate of order + 2 for the 2 -error on the magnetic vector potential. The numerical tests show that removing d ℎ leads to a slightly more accurate approximation of the magnetic vector potential (and a slightly less accurate approximation of the Lagrange multiplier, but this latter result is not meaningful since two different error measures are used).

Ω

  Figure 1: Relative energy-error u ℎ -I +1 X,ℎ X,ℎ / I +1 X,ℎ X,ℎ (top row) and 2 -error

Ω

  Figure 2: Relative energy-error u ℎ -I +1 X,ℎ X,ℎ / I +1 X,ℎ X,ℎ (top row) and 2 -error

)

  Z,♭,ℎ , and from the combination of Remark 15 with the generalized discrete Weber inequality (29) of Corollary 6 applied to (u ℎ , p ℎ ) satisfying (28) (one can easily check that d ℎ satisfies (27)) to bound u ℎ Ω .

Ω

  Figure 4 the results in dashed lines. We plot the relative energy-error u ℎ -(middle row), and since we remove the contribution d ℎ , we replace the 2 -like-error measure p ℎ -I +1 Y,ℎ Y,♭,ℎ (bottom row) by the measure ∈ T ℎ ℎ 2 +1 p -I +1 Y, ( | ) 2 1 /2

10 - 1 . 4

 14 Figure 3: Relative energy-error u ℎ -I +1 X,♭,ℎ X,♭,ℎ / I +1 X,♭,ℎ X,♭,ℎ (top row), 2 -error

  0 is the unique solution to Problem (36). Notice that, since = 0 in Ω, we actually have p

	ℎ	= 0	ℎ and ϵ	ℎ = p	ℎ	. Recalling (37) and (38), the errors (e	ℎ , ϵ	ℎ

  w , , (curl w × ) .

	By the Cauchy-Schwarz inequality, a discrete trace inequality (see, e.g., [17, Lemma 1.32]), and
	recalling the definition (59a) of a ℎ , we get curl ℎ w ℎ	2 Ω	a ℎ (w	ℎ , w	ℎ ), and the first inequality
	in (70a) follows by adding s ℎ (w	ℎ , w	ℎ ) to both sides. To prove the second inequality, we test (61) with
	=	w ∈ R ( ) to infer

  , using a triangle inequality, a discrete trace inequality (see, e.g.,[START_REF] Di Pietro | The Hybrid High-Order method for polytopal meshes. Design, analysis, and applications[END_REF] Lemma 1.32]) on R, (curl | )curl

	+1 P, ( | ) , and the approximation properties of +1 P, on mesh faces (see, e.g., [17, Theorem
	1.45]) for curl	+1 P, (

86) Let us focus on R, (curl | )curl for ∈ F . Adding/subtracting curl +1 P, ( | ) | )curl , we infer

  88)Now, applying the Cauchy-Schwarz inequality, recalling the definition (59c) of d ℎ , noticing that ) | for all ∈ T ℎ and ∈ F , and using the 2 ( )-boundedness of

	P, ( | ) | = +1 P, P, ( +1 P, , we infer
	|d ℎ ( p

| ℎ , q ℎ

Acknowledgments

The authors thank Lorenzo Botti (University of Bergamo) for giving them access to his 3D C++ code SpaFEDTe (https://github.com/SpaFEDTe/spafedte.github.com). Florent Chave and Simon Lemaire also acknowledge support from the LabEx CEMPI (ANR-11-LABX-0007-01).