Éphane St

Lanteri

Diego Paredes
email: diego.paredes@pucv.cl

Claire Scheid
email: claire.scheid@inria.fr

F R Éd

Éric Valentin
email: valentin@lncc.br

THE MULTISCALE HYBRID-MIXED METHOD FOR THE MAXWELL EQUATIONS IN HETEROGENEOUS MEDIA *

Keywords: Maxwell equations, wave model, hybrid method, finite element, multiscale AMS subject classifications. 65N12, 65N15, 65N30

In this work, we address time dependent wave propagation problems with strong multiscale features (in space and time). Our goal is to design a family of innovative high performance numerical methods suitable to the simulation of such multiscale problems. Particularly, we extend the multiscale hybrid-mixed (MHM) finite element method for the two-and three-dimensional time-dependent Maxwell equations with heterogeneous coefficients. The MHM method arises from the decomposition of the exact electric and magnetic fields in terms of the solutions of locally independent Maxwell problems tied together with a one-field formulation on top of a coarse-mesh skeleton. The multiscale basis functions, which are responsible for upscaling, are driven by local Maxwell problems with tangential component of the magnetic field prescribed on faces. A high-order discontinuous Galerkin method in space combined with a second-order explicit leap-frog scheme in time discretizes the local problems. This makes the MHM method effective and yields a staggered algorithm within a "divide-and-conquer" framework. Several two-dimensional numerical tests assess the optimal convergence of the MHM method and its capacity to preserve the energy principle, as well as its accuracy to solve heterogeneous media problems on coarse meshes.

1. Introduction. Partial differential equations embedding multiscale features occur in a wide range of scientific and technological applications involving wave propagation in heterogeneous media. We consider here the case of electromagnetic wave propagation, more precisely, light propagation (e.g., optical wave) in interaction with nanometer scale structures, i.e., nanophotonics. Nanophotonics is the field of science and technology aiming at establishing and using the peculiar properties of light and light-matter interaction in various nanostructures. Nanostructured photonic crystal are very attractive as these structures serve as excellent waveguides. In fact, light paths are created inside these structures removing background electromagnetic modes and keeping only the desired band of frequencies. Generally the medium is a periodic structure, but interesting phenomena arise when band gaps are included to increase (or decrease) local dielectric function for a certain frequency range. Such a gap breaks the periodicity and creates simple cavity modes. These modes can be combined to generate more complex structures where the fields deviate even further from symmetry, pushing them toward the photonic band gap. Overall, photonic-crystal waveguiding structures have led to the design of optical components such as sensors or/and logic gates that are orders of magnitude smaller than currently available devices (see [START_REF] Busch | Periodic nanostructures for photonics[END_REF] for an interesting review).

Such a scenario represents a big challenge for numerical methods in view of computational simulations. In a broad sense, very fine meshes must be adopted to approximate the high-frequency field distributions as they are greatly impacted by the multiscale structures of the medium and the anisotropic nature of the cavity regions. The idea of using high-order polynomial interpolations on coarse meshes, as proposed in [START_REF] Hesthaven | High-order accurate methods in time-domain computational electromagnetics[END_REF], for instance, is not a feasible option since interfaces between different material media and faces of the partition must coincide to avoid spurious numerical modes. Consequently, effective numerical methods for long-time simulations require a fine mesh to fit the complex geometry with high-order approximation to minimize dispersion. This leads to cost-prohibitive large-scale computations, placing realistic three-dimensional problems out of reach. In the quest to overcome such a shortcoming, high-order numerical methods have been devised on building-block structures so as to better match massive parallel facilities. Indeed, it is natural to strive for alternatives to standard techniques such as the finite difference FDTD scheme [START_REF] Collino | Conservative space-time mesh refinement methods for the FDTD solution of Maxwells equations[END_REF] or the finite volume methods for which parallelization was not a goal-driven feature. An interesting alternative is the discontinuous Galerkin method for time domain (DGTD) (see [START_REF] Catella | An implicit DGTD method for solving the two-dimensional Maxwell equations on unstructured triangular meshes[END_REF][START_REF] Descombes | Recent advances on a DGTD method for time-domain electromagnetics[END_REF][START_REF] Leger | A parallel non-conforming multi-element DGTD method for the simulation of electromagnetic wave interaction with metallic nanoparticles[END_REF]) which encompass some of these requirements at the price of increasing the total number of degrees of freedom. A bias is that the DGTD method still relies on (global) fine meshes so as to incorporate the influence of the heterogeneous geometries into numerical solutions.

Multiscale finite element methods have attracted great attention by their capacity to be accurate on coarse meshes (see [START_REF] Abdulle | The heterogeneous multiscale method[END_REF][START_REF] Arbogast | Analysis of a two-scale locally conservative subgrid upscaling for elliptic problems[END_REF][START_REF] Arbogast | A multiscale mortar mixed finite element method[END_REF][START_REF] Efendiev | Multiscale finite element problems with highly oscillatory coefficients[END_REF] just to cite a few). Starting with the seminal work by Babuska and Osborn [START_REF] Babuska | Generalized finite element methods: Their performance and their relation to mixed methods[END_REF], the multiscale methods were further extended to higher dimension in [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF]. The approach carries the concept of multiscale basis functions, which upscale to the coarse mesh. They are driven by independent local problems which make the multiscale methods particularly attractive to implement in parallel environments. Recently, the classical hybridization procedure [START_REF] Raviart | Primal hybrid finite element methods for 2nd order elliptic equations[END_REF] has been used to devise a new family of high-order multiscale finite element methods called the multiscale hybrid-mixed (MHM) method. Particularly adapted to handle multiscale and/or high-contrast coefficients on coarse meshes, the MHM method permits face-crossing interfaces endowed in local boundary conditions. Originally proposed for the Laplace problem in [START_REF] Harder | A family of multiscale hybrid-mixed finite element methods for the Darcy equation with rough coefficients[END_REF], the MHM method has been analyzed in [START_REF] Araya | Multiscale hybrid-mixed methods[END_REF] and further extended to other elliptic problems in [START_REF] Harder | A hybrid-mixed method for elasticity[END_REF][START_REF] Harder | On a multiscale hybrid-mixed method for advective-reactive dominated problems with heterogenous coefficients[END_REF] and mixed problems in [START_REF] Araya | Multiscale hybrid-mixed method for the Stokes and Brinkman equations-the method[END_REF][START_REF] Pereira | A locking-free MHM method for elasticity[END_REF]. Also, its robustness with respect to (small) physical parameters has been established in [START_REF] Paredes | On the robustness of multiscale hybrid-mixed methods[END_REF], and an abstract general setting to develop and analyze MHM methods has been proposed in [START_REF] Harder | Foundations of the MHM method[END_REF].

This work extends the MHM method to the first-order Maxwell system in time domain, with the upshot that the physical coefficients may be highly heterogeneous as found in wave propagation in nano-structured problems. In the vast literature on numerical methods for wave propagation models, some methods share similarities (and/or the same goals) with the MHM methods. One is the DG method in [START_REF] Cockburn | Locally divergence-free discontinuous galerkin methods for the Maxwell equations[END_REF] proposed for the two-dimensional first-order Maxwell equation with constant coefficients. The hybrid discontinuous Galerkin methods first presented in [START_REF] Cockburn | A characterization of hybridized mixed methods for second order elliptic problems[END_REF] and extended to the frequency and time domain Maxwell equation in homogeneous media in [START_REF] Li | Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2D time-harmonic Maxwell's equations[END_REF][START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF] and [START_REF] Christophe | An implicit hybridized discontinuous Galerkin method for the 3D time-domain Maxwell equations[END_REF] share with the MHM method a similar global-local algorithm. Others include the heterogeneous multiscale method for the wave equation [START_REF] Abdulle | Finite element heterogeneous multiscale method for the wave equation: Long time effects[END_REF], the discontinuous Petrov-Galerkin method [START_REF] Gopalakrishnan | Dispersive and dissipative errors in the DPG method with scaled norms for Helmholtz equation[END_REF], and the localizable orthogonal decomposition method [START_REF] Peterseim | Eliminating the pollution effect in Helmholtz problems by local subscale correction[END_REF] for the Helmholtz equation with constant coefficient or also numerical methods built up on top of homogenization techniques or asymptotic expansions as the one presented in [START_REF] Gopalakrishnan | Asymptotic and numerical techniques for resonances of thin photonic structures[END_REF] for the Helmholtz and in [START_REF] Zhanga | A multiscale approach and a hybrid Fe-Be algorithm for heterogeneous scattering of Maxwell's equations[END_REF] for the frequency domain Maxwell equations with multiscale coefficients. However, the primal hybridization of the Maxwell model selected as the starting point in this work as well as the nature of the solution decomposition leads to a different global-local family of numerical methods compared with the ones proposed in the aforementioned papers. Also, the MHM method may be seen as a genuine multiscale method, in principle, as it embeds built-in upscaling through multilevel submeshes if necessary. Thereby the proposed method is, to the best of our knowledge, the first in the literature of multiscale finite element method to be devised for the first-order Maxwell equations in time domain with heterogeneous coefficients.

1.1. The model and main results. To highlight the main points involved in the construction of the MHM method, let Ω ⊂ R 3 be an open bounded simply connected domain with Lipschitz boundary ∂Ω, where ∂Ω is a perfect electrically conducting boundary and (0, T ) is the time interval, with T > 0 fixed. We consider the Maxwell problem to find the electric field e : (0, T ) × Ω → R 3 and the magnetic field h : (0,

T ) × Ω → R 3 such that                    ε ∂ t e -∇ × h = f in (0, T ) × Ω, µ ∂ t h + ∇ × e = 0 in (0, T ) × Ω, ∇ • (ε e) = ρ in (0, T ) × Ω, ∇ • (µ h) = 0 in (0, T ) × Ω, e × n = 0 on (0, T ) × ∂Ω, e = e 0 , h = h 0 at t = 0 on Ω, (1.1)
where n is the unit outward normal vector on ∂Ω, f is the electric current density, and ρ is the charge density. They verify the compatibility relation

∇ • f -∂ t ρ = 0 in (0, T ) × Ω, (1.2)
and e 0 and h 0 are given regular functions with values in R 3 . Here µ and ε denote the magnetic permeability and the electric permittivity, respectively, which are in general third-order symmetric tensors possibly depending on x := (x 1 , x 2 , x 3 ) T ∈ Ω and embedding multiple geometrical scales. They are both measurable, uniformly positive definite, and bounded, i.e., there exist positive constants c µ min and c µ max , c ε min and c ε max such that

c µ min |ζ| 2 ≤ ζ T µ(x) ζ ≤ c µ max |ζ| 2 for all ζ ∈ R 3 , for all x ∈ Ω, (1.3) and c ε min |ζ| 2 ≤ ζ T ε(x) ζ ≤ c ε max |ζ| 2 for all ζ ∈ R 3 , for all x ∈ Ω, (1.4)
where |.| stands for the Euclidean norm. We recall that, for a given vector function

v := (v 1 , v 2 , v 3 ) T with values in R 3 , ∇ × v is a vector function with values in R 3 defined by ∇ × v := (∂ x2 v 3 -∂ x3 v 2 , ∂ x3 v 1 -∂ x1 v 3 , ∂ x1 v 2 -∂ x2 v 1 ) T and ∇ • v is a scalar function with value in R defined by ∇ • v := ∂ x1 v 1 + ∂ x2 v 2 + ∂ x3 v 3 .
The MHM approach is constructive. It starts from a hybrid formulation of the Maxwell model (1.1) defined on a (coarse) partition T H of Ω, formed by elements K ∈ T H , and on a partition T ∆t of the time interval [0, T ] decomposed in elements I n := [t n-1 , t n ] ∈ T ∆t , where n ∈ {1, . . . , N }, N ∈ N * . Next, the continuity of the variables is relaxed on the set of boundaries ∂K for all K ∈ T H , hereafter denoted by ∂T H . A weak continuity of the tangential component of the electric field e is imposed on the boundary elements in ∂T H , whereas the continuity in time is imposed strongly on each time-step t n . More specifically, we propose the following hybrid weak formulation: For each

I n ∈ T ∆t , find (e n ; h n , λ n ) ∈ C 0 (I n ; V)×C 0 (I n ; V)×C 0 (I n ; Λ) such that                          K∈T H K ε ∂ t e n • v - K h n • ∇ × v + λ n , v ∂K = K∈T H K f • v for all v ∈ V, K∈T H K µ ∂ t h n • w + K∈T H K ∇ × e n • w = 0 for all w ∈ V, K∈T H ν, e n ∂K = 0 for all ν ∈ Λ, e n (t n-1 , •) = e n-1 (t n-1 , •) and h n (t n-1 , •) = h n-1 (t n-1 , •) , (1.5)
with e 0 (0, •) = e 0 and h 0 (0,

•) = h 0 in Ω. (1.6)
Here, we denote by ., . ∂K the paring product between H -1/2 (∂K) and H 1/2 (∂K), for all K ∈ T H , and the space V stands for the broken H 1 space, i.e.,

V := v ∈ L 2 (Ω) : v | K ∈ H 1 (K) for all K ∈ T H , (1.7)
and Λ the space of the restriction of the tangential component of functions in H(curl; Ω) to the boundaries ∂K, i.e.,

Λ := ¶ v × n K | ∂K ∈ H -1/2 (∂K) for all K ∈ T H : v ∈ H(curl; Ω) © , (1.8)
with n K being the outward normal vector on ∂K. The spaces have their usual meaning, and they will be detailed in the next sections. The global solution in time (e, h, λ) is defined such that (e, h, λ)| In = (e n , h n , λ n ) for all I n ∈ T ∆t . Hereafter, it is said that (e, h, λ) satisfy the hybrid formulation (1.5)- (1.6). Importantly, the solution of the classical weak form of (1.1) and the solution (e, h, λ) of (1.5)-(1.6) coincide. This result will be precisely addressed in the next section.

Next, before defining the discretization in space and time, formulation (1.5) is reduced to the statement of a system of locally and globally defined problems. To this end, we use the linearity of the Maxwell equations and decompose e n = e n,λ + e n,f and h n = h n,λ + h n,f , (1.9)

where, from the two first equations in (1.5), the functions (e n,λ , h n,λ ) are defined as the solution of the local Maxwell problem (in a distributional sense)

                 ε ∂ t e n,λ -∇ × h n,λ = 0 in I n × K, µ ∂ t h n,λ + ∇ × e n,λ = 0 in I n × K, n K × h n,λ = -λ n on I n × ∂K,
e n,λ (t n-1 , •) = 0 and h n,λ (t n-1 , •) = 0, e 0,λ = 0 and h 0,λ = 0 in K at t = 0, (1.10) and (e n,f , h n,f ) solves

                 ε ∂ t e n,f -∇ × h n,f = f in I n × K, µ ∂ t h n,f + ∇ × e n,f = 0 in I n × K, n K × h n,f = 0 on I n × ∂K, e n,f (t n-1 , •) = e n-1 (t n-1 , •) and h n,f (t n-1 , •) = h n-1 (t n-1 , •),
e 0,f = e 0 and h 0,f = h 0 in K at t = 0.

(1.11)

Owing to characterization (1.9) the third equation in (1.5) becomes a global one-field face formulation for the Lagrange multiplier, i.e., for all t ∈ I n ∈ T ∆t , find λ n (t) ∈ Λ such that

K∈T H ν, e n,λ ∂K = - K∈T H ν, e n,f ∂K for all ν ∈ Λ. (1.12)
The MHM method is a byproduct of the discretization of the coupled problems (1.10)- (1.12).

The one-level MHM arises from the choice of a finite dimensional space of C 0 (0, T ; Λ) (possibly nonconforming) and looking for the approximate solution λ H (t) of λ(t) in such a space satisfying (1.12) at discrete times. In this work, the Lagrange multipliers λ H (t) are assumed to be constant in each I n ∈ T ∆t and polynomial of high degree on ∂T H . Interestingly, we verify that such a simple choice leads to a second-order time scheme. Also, by going discrete, problems (1.10)-(1.12) may be decoupled as follows:

• the local problem (1.10) is responsible for computing the multiscale basis, with its boundary condition λ n replaced by the space-time basis of λ H ; • the global problem (1.12) is responsible for the calculation of the degrees of freedom of λ H , which depend on I n ∈ T ∆t . As a result, the one-level approximation (e H , h H ) of (e, h) stems from the combination of the multiscale basis functions, the degrees of freedom of λ H , and the solution of (1.11). A detailed description of (e H , h H ) is the subject of section 4.

The aforementioned one-level MHM method relies on closed formulas for the multiscale basis functions and for the solution of (1.11) which are not available in general. As such, one has to resort to a second level of discretization to make the method (1.12) effective. Interestingly, such choices can be quite general in the sense that any numerical method with approximate properties could be locally adopted. However, in this setting, the global problem is naturally impacted by the choice made at the second level. We propose to associate the discontinuous Galerkin (DG) method proposed in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] to discretize the local problems in space with the explicit secondorder leap-frog scheme for the time marching process. These options are particularly attractive in practice as they lead to a fully explicit scheme (see [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] for details) which fulfills the discrete energy principle. Also, the discontinuous interpolation at the local level turns out to be of paramount importance to preserve the quality of the numerical solution in the case of interface crossing faces. A detailed description of the two-level MHM method and its properties will be presented in section 5.

The theoretical foundations of the MHM strategy for the Maxwell equations is one of the main contributions of this work. Particularly, we demonstrate the equivalence between the hybrid formulation (1.5) and the original weak form of the Maxwell equations and we establish the well-posedness of local problems (1.10) and (1.11). In regard to the numerical aspects of the one-and two-level methods, we prove that they preserve the energy conservation properties under a certain CFL condition for the latter. We also dedicate this work to a complete numerical validation of the new MHM method. Specifically, we verify that the method yields optimal convergence for all variable and superconvergence for the electric field in their natural norms, and it outperforms standard methods. Also, a wave propagation problem defined on a heterogeneous nanoscale device shows that the method is accurate on coarse meshes when the faces are not aligned with interfaces. In view of such novelties, the error analysis of the semi-and fully discrete MHM methods are left to be addressed in a future work.

The paper is organized as follows. The next section states the problem and addresses some preliminary theoretical results. The characterization of the exact fields in terms of global-local problems is the subject of section 3. Section 4 is devoted to the one-level MHM method. The fully discrete two-level MHM method is presented in section 5, and section 6 is dedicated to the numerical validation of the MHM method. Appendix A details the proofs of the theoretical result anticipated in the previous sections, and Appendix B presents the MHM method for the transverse magnetic (TM) two-dimensional model. Conclusions follow in section 7.

Statement and preliminaries.

As usual, we assume that the initial electric and magnetic fields e 0 and h 0 satisfy

∇ • (ε e 0 ) = ρ(0, •) and ∇ • (µ h 0 ) = 0 in Ω. (2.1)
Hence, problem (1.1) simplifies observing that the third equation (Gauss law) and the fourth equation (absence of magnetic monopoles) in (1.1) stem from the first (Faraday law) and second (Ampere law) equations in (1.1). As a result, problem (1.1) can be rewritten in an equivalent way as follows: Find e : (0, T ) × Ω → R 3 and h : (0,

T ) × Ω → R 3 such that          ε ∂ t e -∇ × h = f in (0, T ) × Ω, µ ∂ t h + ∇ × e = 0 in (0, T ) × Ω, e × n = 0 on (0, T ) × ∂Ω, e = e 0 , h = h 0 at t = 0 on Ω, (2.2)
constrained by (1.2). The following spaces, equipped with their usual inner-product and induced norms, will be used in what follows:

H(curl; Ω) := v ∈ L 2 (Ω) : ∇ × v ∈ L 2 (Ω) , H 0 (curl; Ω) := {v ∈ H(curl; Ω) : v × n | ∂Ω = 0} , H(div; Ω) := v ∈ L 2 (Ω) : ∇ • v ∈ L 2 (Ω) , H(div, µ; Ω) := v ∈ L 2 (Ω) : ∇ • (µ v) ∈ L 2 (Ω) , H 0 (div, µ; Ω) := {v ∈ H(div, µ; Ω) : v • n | ∂Ω = 0} , H(div, ε; Ω) := v ∈ L 2 (Ω) : ∇ • (ε v) ∈ L 2 (Ω) ,
The standard weak formulation of problem (2.2) consists of finding (e, h) ∈ C 0 (0, T ;

H 0 (curl; Ω)) × C 0 (0, T ; H(curl; Ω)) such that ® (ε ∂ t e, v) Ω -(∇ × h, v) Ω = (f , v) Ω for all v ∈ H(curl; Ω), (µ ∂ t h, w) Ω + (e, ∇ × w) Ω = 0 for all w ∈ H(curl; Ω), (2.3)
where (•, •) D stands for the usual L 2 (D) inner-product, and D is an open set, and we denote by • 0,D the induced norm. We recall that problem (2.3) is well-posed (cf. [START_REF] Ciarlet | Modéles et méthodes pour les équations de Maxwell[END_REF]) as follows.

Theorem 1. Assume that Ω is an open, bounded, simply connected domain with Lipschitz boundary such that the neighborhood of each point of ∂Ω is located in one side of ∂Ω. If the regularity conditions

(i) ρ ∈ C 1 (0, T ; L 2 (Ω)), f ∈ C 0 (0, T ; H(div; Ω)), and ∂ t f ∈ L 2 (0, T ; L 2 (Ω)), (ii) e 0 ∈ H 0 (curl; Ω) ∩ H(div, ε; Ω) and h 0 ∈ H(curl; Ω) ∩ H 0 (div, µ; Ω), hold, then there exists a unique solution (e, h) to (2.3) such that (e, ∂ t e) ∈ C 0 (0, T ; H 0 (curl; Ω) ∩ H(div, ε; Ω)) × C 0 (0, T ; L 2 (Ω)), (h, ∂ t h) ∈ C 0 (0, T ; H(curl; Ω) ∩ H 0 (div, µ; Ω)) × C 0 (0, T ; L 2 (Ω)).
Remark 1. In the case ε and µ are smooth functions in Ω, and Ω is of class C 1,1 or convex, the solution (e, h) of (2.3) belongs to the more regular space C 0 (0, T ; H 1 (Ω))× C 0 (0, T ; H 1 (Ω)) (cf. [START_REF] Amroucche | Vector Potential in Three Dimensional Nonsmooth Domains[END_REF]). When these coefficients are piecewise smooth functions (piecewise constant, for instance) the fields turn out to be piecewise H 1 in the subdomains where the coefficients are smooth (see, e.g., [START_REF] Costabel | Singularities of Maxwell interface problems[END_REF]). Hereafter, we suppose the solutions of (2.3) are H 1 (Ω) regular. To be complete, we shall make some comments on what one should change if this assumption is relaxed every time it is used. Now, instead of working directly with problem (2.3), we adopt the following perspective: we seek (e, h) as the solution of the Maxwell model (2.3) in a weaker, broken space which relaxes continuity and localizes computations. Ultimately, we find that the approach allows for the construction of e and h using local problems. From this perspective, we first partition the space-time domain Ω × [0, T ]. We discretize the time interval [0, T ] in N + 1 points 0 = t 0 < t 1 . . . < t N -1 < t N = T with N ∈ N * (not necessarily) equally distributed, define for n ∈ {1, . . . , N }, I n := [t n-1 , t n ], and collect them in the set T ∆t , where ∆t := max n=1,...,N |I n |. The partition of the space domain Ω can be general (as the partitions used in the virtual finite element method [START_REF] Veiga | Basic principles of virtual element methods[END_REF], for instance). However, we focus here on partitions based on a family of regular meshes {T H } H>0 , where H is the characteristic length of T H . The mesh can be very general, composed of heterogeneous element geometries. Each element K has a boundary ∂K consisting of faces F , and we collect in ∂T H the boundaries ∂K and in F H the faces associated with T H . We denote by F ∂ the set of faces on ∂Ω and by F 0 := F H \ F ∂ the set of internal faces. To each F ∈ F H , we associate a normal n, taking care to ensure this is facing outward on ∂Ω for each F ∈ F ∂ . For each K ∈ T H , we further denote by n K the outward normal on ∂K and let n K F := n K | F for each F ⊂ ∂K. We recall that the broken space V given in (1.7) and the space

W := v ∈ L 2 (Ω) : v | K ∈ H(curl; K) for all K ∈ T H
are Hilbert spaces with the respective scalar products

(u, v) V := K∈T H (u, v) K + K∈T H (∇u, ∇v) K , (u, v) W := K∈T H (u, v) K + K∈T H (∇ × u, ∇ × v) K ,
and we denote their induced norms by • V and • W , respectively. Also, the space Λ given in (1.8) is a Banach space when equipped with the following quotient norm:

ν Λ := inf v∈H(curl;Ω) v×n K | ∂K =ν,∀K∈T H v W .
We denote by (•, •) T H and (•, •) ∂T H the summation of the respective inner (or dual) products over the sets K and ∂K, respectively, given by

(w, v) T H := K∈T H K w • v dx and (ν, v) ∂T H := K∈T H ν, v ∂K ,
where w, v ∈ V and ν ∈ Λ. Here, •, • ∂K is uniquely defined via the dual product involving H -1/2 (∂K) and

H 1/2 (∂K) through ν, v ∂K := K ∇ × w • v dx - K w • ∇ × v dx,
where w × n K = ν on ∂K, and we denote it by (•,

•) ∂K if ν ∈ L 2 (∂K).
Owing to this notation and setting for n ∈ {1, . . . , N }

D n := C 0 (I n ; V) ∩ C 1 (I n ; L 2 (Ω)) × C 0 (I n ; V) ∩ C 1 (I n ; L 2 (Ω)) × C 0 (I n ; Λ), (2.4)
we propose the following hybrid version of problem (2.3). Let n ∈ {1, . . . , N }. For all

I n ∈ T ∆t , find (e n , h n , λ n ) ∈ D n such that                (ε ∂ t e n , v) T H -(h n , ∇ × v) T H + (λ n , v) ∂T H = (f , v) T H for all v ∈ V, (µ ∂ t h n , w) T H + (∇ × e n , w) T H = 0 for all w ∈ V, (ν, e n ) ∂T H = 0 for all ν ∈ Λ, e n (t n-1 , •) = e n-1 (t n-1 , •) and h n (t n-1 , •) = h n-1 (t n-1 , •) on Ω, e 0 (0, •) = e 0 and h 0 (0, •) = h 0 at t = 0 on Ω. (2.5)
We emphasize that the first three equations in (2.5) should be read as being valid for all t ∈ I n (fixed although arbitrary), and the global solution (e, h, λ) is defined such that (e, h, λ) | In = (e n , h n , λ n ), for all I n ∈ T ∆t , and hereafter we will say that (e, h, λ) is a solution of (2.5).

Remark 2. When the V-regularity assumption for solutions is relaxed, the hybrid formulation (2.5) still makes sense by replacing the space V by W and interpreting the operator (•, •) ∂T H as the unique density extension to Λ × W of (•, •) ∂T H defined on Λ × V (see Appendix A for details).

Observe that the variables in (2.5) belong a priori to a larger space than the solution of the original problem (2.3). However, the space of Lagrange multipliers imposes H(curl; Ω)-conformity on the solution and the respect of the boundary condition e × n | ∂Ω = 0. Also, problem (2.5) has a unique solution (in the sense of Remark 2) where (e, h) coincides with the solution of problem (2.3). These results are summarized next and the proof detailed in Appendix A.

Lemma 2. The function (e, h, λ) is the solution of (2.5) if and only if (e, h) solves (2.3) in the sense of Theorem 1. Furthermore, it holds that λ = -n K × h on (0, T ) × ∂K.

(2.6) Remark 3. Absorbing boundary conditions can be prescribed in (2.2). They are taken into account as a natural condition within the hybrid formulation (2.5). Indeed, assume that the following boundary condition is set on a subset ∂Ω A of the boundary ∂Ω:

n × e -n × (α h × n) = n × g on (0, T ) × ∂Ω A , (2.7)
where the given right-hand-side g is defined by

g := e i -α h i × n and α 2 := ε -1 µ.
The functions (e i , h i ) are the incident electric and magnetic fields, respectively. Hence, using equivalence (2.6), condition (2.7) is prescribed by replacing the third equation in (2.5) by

(ν, e) ∂T H = (ν, g) ∂Ω A -(ν, α λ) ∂Ω A for all ν ∈ Λ.
We are ready to propose an equivalent global-local formulation to (2.5) which turns out to be more convenient to build "divide and conquer" numerical algorithms.

A global-local formulation.

Rather than selecting a pair of finite subspaces of V × V × Λ at this point, we rewrite (2.5) in an equivalent form which is suitable to reduce the statement to a system of locally and globally defined problems. Such an approach guides the definition of stable finite subspaces composed of functions which incorporate multiple scales into the basis functions.

First, observe that problem (2.5) is equivalent, for all I n ∈ T ∆t , to find (e n , h n , λ n ) ∈ D n such that (ν, e n ) ∂T H = 0 for all ν ∈ Λ, (3.1) and for all K ∈ T H ,

     (ε ∂ t e n , v) K -(h n , ∇ × v) K = (f , v) K -λ n , v ∂K , (µ ∂ t h n , w) K + (∇ × e n , w) K = 0, e n (t n-1 , •) = e n-1 (t n-1 , •) and h n (t n-1 , •) = h n-1 (t n-1 , •) on K, (3.2) for all (v, w) ∈ V(K) × V(K)
, where we adopt the convention (e 0 (0, •), h 0 (0, •)) = (e 0 , h 0 ) and V(K) stands for the space of functions in V restricted to K. The corresponding electric and magnetic fields in (3.1)-(3.2) respect the classical energy principle. This is the subject of the next lemma. Lemma 3. Let E : (0, T ) → R be defined as follows:

E(t) := 1 2 (ε e, e) T H + 1 2 (µ h, h) T H , (3.3)
where (e, h) solve (3.1)-(3.2). Then, it holds that

E(t) = E(0) + t 0 (f , e) T H for all t ∈ (0, T ). (3.4)
Proof. Let I n ∈ T ∆t , and choose e n and h n as test functions in the first and second equations of (3.2), respectively. This leads us to

     (ε ∂ t e n , e n ) K -(h n , ∇ × e n ) K = (f , e n ) K -λ n , e n ∂K , (µ ∂ t h n , h n ) K + (∇ × e n , h n ) K = 0, e n (t n-1 , •) = e n-1 (t n-1 , •) and h n (t n-1 , •) = h n-1 (t n-1 , •) on K. (3.5)
Summing up the resulting equations over all the elements K of T H , we arrive at

(ε ∂ t e n , e n ) T H + (µ ∂ t h n , h n ) T H = (f , e n ) T H -(λ n , e n ) ∂T H . (3.6) Next, using (3.1) in (3.6), we get (ε ∂ t e n , e n ) T H + (µ ∂ t h n , h n ) T H = (f , e n ) T H ,
and from the regularity of e n , we conclude that

∂ t E = (f , e) T H . (3.7)
The result follows by integrating (3.7) in the interval [0, t], with t ∈ (0, T ).

Notice that (3.2) implies that e and h can be computed in each element K ∈ T H from f and λ once the latter is known. Owing to such property, we can use the linearity of the Maxwell operator to decompose the fields as e n = e n,λ + e n,f and h n = h n,λ + h n,f for all n ∈ {1, . . . , N }, (3.8) where (e n,λ , h n,λ ) satisfies

       (ε ∂ t e n,λ , v) K -(h n,λ , ∇ × v) K = -λ n , v ∂K for all v ∈ V(K), (µ ∂ t h n,λ , w) K + (∇ × e n,λ , w) K = 0 for all w ∈ V(K), e n,λ (t n-1 , •) = 0 and h n,λ (t n-1 , •) = 0, (3.9)
and (e n,f , h n,f ) solves

       (ε ∂ t e n,f , v) K -(h n,f , ∇ × v) K = (f , v) K for all v ∈ V(K), (µ ∂ t h n,f , w) K + (∇ × e n,f , w) K = 0 for all w ∈ V(K), e n,f (t n-1 , •) = e n-1 (t n-1 , •) and h n,f (t n-1 , •) = h n-1 (t n-1 , •), (3.10)
with the same convention as above (e 0 , h 0 )(0, •) := (e 0 , h 0 ). The corresponding strong versions of (3.9) and (3.10) (meant in a distributional sense) were anticipated in (1.10) and (1.11), respectively. Also, the initial boundary value problems (3.9) and (3.10) are well-posed as proved in Lemma 6 (see Appendix A).

Owing to characterization (3.8) and replacing it in (3.1), we can rewrite the global problem (3.1) in an equivalent form in which the Lagrange multipliers are the leading variable, i.e., for each I n ∈ T ∆t and t ∈ I n , we look for λ n (t) ∈ Λ such that (ν, e n,λ ) ∂T H = -(ν, e n,f ) ∂T H for all ν ∈ Λ.

) ∂T H + (ν, α λ n ) ∂Ω A = -(ν, e n,f ) ∂T H + (ν, g) ∂Ω A for all ν ∈ Λ.
4. Space-time discretization: The first level. Since λ uniquely determines e λ , it is enough to select a finite dimensional space Q ∆t H of C 0 (0, T ; Λ) (possibly nonconforming) in order to discretize problem (3.11) in space and time.

To this end, we define P 0 (I n ) the space of piecewise constant functions in the time intervals I n ∈ T ∆t and set Λ H ⊂ Λ as the following polynomial space:

Λ H = Λ m l := ν H ∈ Λ : ν H | F ∈ P m l (F ) 3 for all F ∈ F H . (4.1)
Here, the space P m l (F ) is the space of discontinuous polynomial functions on F of degree less than or equal to l ≥ 0 and defined on an equally spaced partition of F , composed of m elements (m ≥ 1). In what follows, the space Λ l stands for Λ 1 l . Owing to these choices, we set

Q ∆t H := In∈T∆t P 0 (I n ) ⊗ Λ H , (4.2)
where Λ H is given in (4.1), and observe that such a space is nonconforming in C 0 (0, T, Λ) with respect to the time variable. This corresponds to the simplest spacetime element. More involved time interpolation choices may be used with the upside that the conformity is recovered. It is worth mentioning that the simplicity of (4.2) associated with its second-order accuracy in time (see section 6) makes the space (4.2) the best compromise between precision and computational cost.

Next, mimicking (3.8), the one-level approximated fields are decomposed as

e n H = e n,λ H + e n,f and h n H = h n,λ H + h n,f , (4.3)
and the one-level MHM method results from the standard Galerkin method applied to (3.11), i.e., for all

I n ∈ T ∆t , find λ n H ∈ Λ H such that (ν H , e n,λ H ) ∂T H = -(ν H , e n,f ) ∂T H for all ν H ∈ Λ H , (4.4)
where λ H | In = λ n H . The function (e n,f , h n,f ) solves (3.10) and, using (3.9), the function (e n,λ H , h n,λ H ) satisfies

       ε ∂ t e n,λ H , v K -h n,λ H , ∇ × v K = -(λ n H , v) ∂K for all v ∈ V(K), µ ∂ t h n,λ H , w K + ∇ × e n,λ H , w K = 0 for all w ∈ V(K), e n,λ H (t n-1 , •) = 0 and h n,λ H (t n-1 , •) = 0 on Ω. (4.5)
Well-posedness of (4.5) follows from Lemma 6 (see Appendix A) and the corresponding strong problem associated to (4.5) is given in (1.10) (with λ n replaced by λ n H ). We analogously define the one-level global electric and magnetic fields e H and h H such that their restrictions to I n ∈ T ∆t coincide with e n H and h n H , respectively. Hence, the local conservation property of their exact counterpart (under the same assumption (2.1)) holds following the proof in Lemma 3. Specifically, the discrete energy E H given by

E H (t) := 1 2 (ε e H , e H ) T H + 1 2 (µ h H , h H ) T H (4.6)
is preserved as follows:

E H (t) = E H (0) + t 0 (f , e H ) T H . (4.7)
It is instructive to consider (e n,λ H , h n,λ H ) in more detail as it plays a central role in (4.4). Let {ψ i }, with i ∈ {1, . . . , dim Λ H }, be a basis for Λ H and set {ψ i 1 In }, with (i, n) ∈ {1, . . . , dim Λ H } × {1, . . . , N } a basis for Q ∆t H , where 1 In stands for the characteristic function of I n ∈ T ∆t . If β n i ∈ R are the corresponding degrees of freedom associated to the basis of Q ∆t H , then we write

λ n H = dim Λ H i=1 β n i ψ i for all n ∈ {1, . . . , N }. (4.8) Now, define the set (η n,e i , η n,h i ) ⊂ D n × D n , with i ∈ {1, . . . , dim Λ H }, such that (η n,e i , η n,h i ) | K×In satisfies the following Maxwell problem for all K × I n ∈ T H × T ∆t : ε ∂ t η n,e i , v K -η n,h i , ∇ × v K = -(ψ i , v) ∂K for all v ∈ V(K), µ ∂ t η n,h i , w K + ∇ × η n,e i , w K = 0 for all w ∈ V(K), (4.9)
where (η n,e i (t n-1 , •), η n,h i (t n-1 , •)) = (0, 0). Here ψ i changes its sign according to the sign of n K F • n, for each F ⊂ ∂K, as a consequence of Λ H in (4.1). As a result of (4.8) and (4.9), we can write (e n,λ H , h n,λ H ) the solution of (4.5) in the form

e n,λ H = dim Λ H i=1 β n i η n,e i and h n,λ H = dim Λ H i=1 β n i η n,h i on I n ∈ T ∆t , (4.10)
and then from (4.3), on each I n ∈ T ∆t we get

e n H = dim Λ H i=1 β n i η n,e i + e n,f and h n H = dim Λ H i=1 β n i η n,h i + h n,f . (4.11)
The one-level approach assumes that close formulas are available for (η n,e i , η n,h i ) and (e n,f , h n,f ). As such, to compute (4.11) effectively, it remains to obtain the degrees of freedom β n i . To this end, we use the face-based global formulation (4.4) evaluated on N discrete time-steps (one per time slab I n ). Indeed, replacing (4.10) in (4.4) it results that β n i , for i ∈ {1, . . . , dim(Λ H )} and n ∈ {1, . . . , N }, satisfy the following discrete system:

dim Λ H i=1 K∈T H β n i ∂K ψ j η n,e i ds = - K∈T H ∂K
ψ j e n,f ds for all j ∈ {1, . . . , dim(Λ H )}, where η n,e i and e n,f are evaluated at one time per interval I n ∈ T ∆t (for instance, at t n+1/2 := tn-1+tn 2 ).

Remark 5. Observe that problem (4.9) is equivalent (in a distributional sense) to the following strong problem:

             ε ∂ t η n,e i -∇ × η n,h i = 0 in I n × K, µ ∂ t η n,h i + ∇ × η n,e i = 0 in I n × K, n K × η n,h i = -ψ i on I n × ∂K, η n,e i (t n-1 , •) = 0 and η n,h i (t n-1 , •) = 0 in K. (4.12)
Thereby, the functions (η n,e i , η n,h i ) are in one-to-one correspondence with the basis ψ i and are responsible for the upscaling since they are designed to naturally incorporate the influence of physical coefficients and geometry. This multiscale basis spans a conforming local subspace of D n × D n . As usual, the global basis functions (η e i , η h i ) and (e f , h f ) are defined such that their restrictions coincide with (η n,e i , η n,h i ) and (e n,f , h n,f ) on I n ∈ T ∆t , respectively, and we notice that the support of the multiscale basis functions is restricted to two elements K ∈ T H with clear (positive) impact on the sparsity of the underlying matrices.

Remark 6. The one-level MHM method is nonconforming in H 1 (Ω) as the normal component of the variables may be discontinuous across faces for all t ∈ (0, T ). On the other hand, it is an H(curl; Ω)-conforming method for the magnetic field since its tangential component is continuous across faces (see (1.10))

As mentioned, the one-level MHM method assumes that the basis functions (η e i , η h i ) and (e f , h f ) are known exactly. However, this is generally not feasible, meaning a second level of discretization is necessary to make the method practical. This is the subject of the next section.

The two-level MHM method.

5.1. The semidiscrete MHM method. We discretize local problems (4.9) (or (4.12)) and (3.10) with respect to the space variable by selecting a local finite dimensional space V h (K) whose functions are defined over a regular partition of K, denoted by T K h h>0 . Such a partition may differ in each K ∈ T H and h denotes the characteristic length of T K h . Also, let F K h be the set of faces on T K h , and let F K 0 be the set of internal faces. To each face γ ∈ F K h we associate a normal vector n τ γ , taking care to ensure this is facing outward on ∂τ , for all τ ∈ T K h . The global finite dimensional space is then given by

V h := K∈T H V h (K), (5.1)
and on top of (5.1) a second-level numerical method is proposed. First observe that the MHM algorithm is quite general, and standard face-based elements [START_REF] Nédélec | Mixed finite elements in R 3[END_REF] or nodal continuous stabilized methods [START_REF] Badia | A nodal-based finite element approximation of the maxwell problem suitable for singular solutions[END_REF], for instance, can be incorporated into the MHM algorithm at this point.

In this work, we focus on a simple nodal, nonconforming finite element space defined as a piecewise discontinuous polynomial space in each element K ∈ T H , i.e., V h (K) ⊂ V(K) is given by

V h (K) := v h ∈ L 2 (K) : v h | τ ∈ [P k (τ )] 3 for all τ ∈ T K h , (5.2)
where P k (τ ) is the space of polynomial functions on τ of degree less than or equal to k ≥ 1. Notice that, whatever choice of numerical method we make at the second level, the format of the global problem remains unchanged and reads as follows: For all

I n ∈ T ∆t , find λ n H ∈ Λ H such that (ν H , e n,λ H h ) ∂T H = -(ν H , e n,f h ) ∂T H for all ν H ∈ Λ H , (5.3)
where e n,λ H h and e n,f h belong to V h (as well as the corresponding h n,λ H h and h n,f h ). They correspond to the second-level approximation of e n,λ H and e n,f . As a result, the semidiscrete approximation of the exact solution is given by (e, h) ≈ (e H,h , h H,h ) := (e λ H h + e f h , h λ H h + h f h ), (5.4) where, as usual, the global functions in time (e λ H h , h λ H h ) and (e f h , h f h ) are defined through their local (in time) counterpart (e n,λ H h , h n,λ H h ) and (e n,f h , h n,f h ). Now, we detail the nonconforming method [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] to compute the semidiscrete basis functions of (e λ H h , h λ H h ) and (e f h , e f h ) needed to solve (5.3). Seen as an alternative to the standard Galerkin method on continuous polynomial interpolations, the DG method (within the MHM framework) consists of relaxing the continuity of the interpolator on the faces γ of a submesh T K h for all K ∈ T H . Such an option can be particularly attractive when high-contrast heterogeneities still persist within local problems, making their resolution computationally expensive. In these cases a second level of parallelization may be needed. Another important application of the DG method lies in when material interfaces are not aligned with faces. In such a case, the space Λ m l in (4.1) with m ≥ 1 must be adopted, and then it is important to take into account such discontinuities with precision within the second-level approximations. This attractive feature of the MHM method will be explored in the numerical section (see section 6).

To introduce the DG method proposed in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] in the context of the MHM methodology, we need first some extra notation. Let us denote by

{v h } := 1 2 ((v h ) | τ + (v h ) | τ )
the mean value of a vector function v h on faces γ ∈ F K h (with τ and τ being the two elements sharing the face γ), which coincides with the value of v h if γ ⊂ ∂K. We define the jump on a face γ ∈ F K h as

v h := n τ γ × v h | τ + n τ γ × v h | τ on γ ⊂ ∂τ ∩ ∂τ , where τ, τ ∈ T K h , (5.5) with v h = n τ γ × v h | τ if γ ⊂ ∂K,
and the inner-products on the submesh T K h are given by

(•, •) T K h := τ ∈T K h (•, •) τ and (•, •) F K h := γ∈F K h (•, •) γ .
Owing to these definitions, and for each K ∈ T H , the DG method applied to (4.5) consists of finding (e n,λ

H h , h n,λ H h ) ∈ C 0 (I n , V h (K)) × C 0 (I n , V h (K)), for all I n ∈ T ∆t , such that            ε ∂ t e n,λ H h , v h T K h -h n,λ H h , ∇ × v h T K h + Ä ¶ h n,λ H h © , v h ä F K 0 = -(λ n H , v h ) ∂K , µ ∂ t h n,λ H h , w h T K h + e n,λ H h , ∇ × w h T K h - Ä ¶ e n,λ H h © , w h ä F K h = 0, e n,λ H h (t n-1 , •) = 0 and h n,λ H h (t n-1 , •) = 0 (5.6) for all (v h , w h ) ∈ V h (K) × V h (K).
Similarly, the DG method applied to (3.10) consists of finding (e n,f h ,

h n,f h ) ∈ V h (K) × V h (K) such that            ε ∂ t e n,f h , v h T K h -h n,f h , ∇ × v h T K h - Ä ¶ h n,f h © , v h ä F K 0 = (f , v h ) T K h , µ ∂ t h n,f h , w h T H,h + e n,f h , ∇ × w h T K h + Ä ¶ e n,f h © , w h ä F K h = 0, e n,f h (t n-1 , •) = e n-1 h (t n-1 , .) and h n,f h (t n-1 , •) = h n-1 h (t n-1 , •) (5.7) for all (v h , w h ) ∈ V h (K) × V h (K)
, with the convention e 0 h (0, •) = e 0 and h 0 h (0, •) = h 0 . Here e 0 and h 0 must be understood as well-chosen projections of the exact initial conditions in the sense that such approximations do not undermine convergence rates.

Remark 7. The second-level discrete energy given by

E H,h (t) := 1 2 (ε e H,h , e H,h ) T H + 1 2 (µ h H,h , h H,h ) T H ,
where (e H,h , h H,h ) are given in (5.4), fulfills the following energy principle:

E H,h (t) = E H,h (0) + t 0 (f , e H,h ) T H .
The proof of this result follows closely the one in Lemma 3, using the central flux formulation presented in the DG method [START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF], and the global problem (5.3). We refer to [START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] for the further details.

The fully discrete MHM method.

In this section, we devise the fully discrete two-level MHM method by setting up a time-marching scheme for the local problems (5.6) and (5.7). An explicit local finite difference time scheme is proposed which yields an effective MHM method. It turns out that, after some manipulations, the scheme is equivalent to the second-order leap-frog time scheme.

First, we need some additional notation associated to the (local) time discretization. We decompose I n in a (nonnecessary) uniform partition, i.e.,

t n = t 0 n < t 1 n < • • • < t Nn n = t n+1 , with N n ∈ N * and ∆t k n = t k+1 n -t k n
, and k = 0, . . . , N n -1, and set ∆t n := max k∈[0,Nn-1] ∆t k n . Observe that such a partition can be set differently in each element K ∈ T h , which naturally induces a local time-stepping scheme. Indeed, it is well known that when combined with an explicit time integration method and in the presence of an unstructured locally refined mesh, high-order DGTD methods suffer from a severe time-step size restriction. A possible alternative to overcome this limitation is to use smaller time-steps, given by a local stability criterion, where the smallest elements are. The local character of DG formulations is a desirable feature for the development of explicit local time-stepping schemes. Such techniques were developed for the second-order wave equation in [START_REF] Diaz | Energy conserving explicit local time-stepping for second-order wave equations[END_REF][START_REF] Diaz | Multi-level explicit local time-stepping methods for second-order wave equations[END_REF]. In [START_REF] Piperno | Symplectic local time stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF], a second-order symplectic local time-stepping DGTD method is proposed for the Maxwell equations in a nonconducting medium, based on the Störmer-Verlet method. Grote and Mitkova derived local time-stepping methods of arbitrarily high accuracy for the Maxwell equations from the standard leap-frog scheme in [START_REF] Grote | Explicit local time-stepping methods for Maxwel's equations[END_REF]. More recently, Grote, Mehlin, and Mitkova [START_REF] Grote | Runge-Kutta-based explicit local time-stepping methods for wave propagation[END_REF] designed Runge-Kutta-based explicit local time-stepping methods for time-dependent wave propagation problems.

The previously cited local time-stepping strategies straightforwardly apply to the local DGTD solvers in the MHM framework proposed in this paper. We adopt here a different viewpoint for illustrating the use of local time-steps. For the sake of clarity, we drop such a dependency in the notation. Next, we explore the intrinsic local timestepping of the MHM method within a simple framework by assuming that T ∆t is uniform with characteristic length ∆t := t n+1 -t n . Moreover, for each I n ∈ T ∆t , we consider a three-point subdivision t 0 n = t n , t 1 n = tn+tn+1 2 and t 2 n = t n+1 (observe that t 1 n = t n+ 1 2 ) for the electric field and a two-point subdivision t 0 n = t n and t 1 n = t n+1 for the magnetic field.

We denote e n,i H,h (resp., h n,i H,h ) the approximated field of e n H,h (resp., h n H,h ) at the discrete local time t i , i ∈ {n -1, n -1 2 , n} (resp., i ∈ {n -1, n}). On each I n , we propose the following scheme to approximate (e n H,h , h n H,h ) at time t i , adopting the DG method in space presented in the last section:

                                                 Ñ ε e n,n-1 2 H,h -e n,n-1 H,h ∆t 2 , v h é T K h - Ä h n,n-1 H,h , ∇ × v h ä T K h + Ä ¶ h n,n-1 H,h © , v h ä F K 0 = -(λ n-1 H , v h ) ∂K + f n-1 , v h T K h , µ h n,n H,h -h n,n-1 H,h ∆t , w h T K h + e n,n-1 2 H,h , ∇ × w h T K h -e n,n-1 2 H,h , w h F K h = 0, Ñ ε e n,n H,h -e n,n-1 2 H,h ∆t 2 , v h é T K h - Ä h n,n H,h , ∇ × v h ä T K h + Ä ¶ h n,n H,h © , v h ä F K 0 = -(λ n H , v h ) ∂K + (f n , v h ) T K h , e n,n-1 H,h = e n-1,n-1 H,h and h n,n-1 H,h = h n-1,n-1 H,h , (5.8)
with the initial conditions e 0,0 H,h = e 0 and h 0,0 H,h = h 0 . Interestingly, by rearranging the contribution terms from each I n ∈ T ∆t , we recognize the explicit second-order leap-frog time scheme. In what follows, we simplify the notation and use e

n+ 1 2 H,h and h n H,h to mean e n,n+ 1 2 H,h
and h n,n H,h . As a result, the fully discrete MHM method reads as follows: For all n ∈ {1, . . . , N }, find (e

n+ 1 2 H,h , h n H,h , λ n H ) ∈ V h × V h × Λ H such that ν H , e n+ 1 2

H,h ∂T H

= 0 for all ν H ∈ Λ H , (5.9)

                     Ç µ h n H,h -h n-1 H,h ∆t , w h å T K h = -e n-1 2 H,h , ∇ × w h T K h + e n-1 2 H,h , w h F K h , Ñ ε e n+ 1 2 H,h -e n-1 2 H,h ∆t , v h é T K h = (f n , v h ) T K h + h n H,h , ∇ × v h T K h -h n H,h , v h F K 0 -(λ n H , v h ) ∂K (5.10) for all (v h , w h ) ∈ V h (K) × V h (K).
Here we set e

1 2 H,h = e 0 + ε -1 ∆t ∇ × h 0 and h 0 H,h = h 0 and f n = f (t n ).
Remark 8. The DG method naturally increases the number of degrees of freedom needed to approximate the basis functions. However, we recall that this is an attractive alternative since the two-level MHM method would now enjoy the flexibility to be implemented using a multilevel parallel strategy.

It is important to verify the conditions under which the discrete energy principle associated with the solution of the fully discrete MHM method (5.9)-(5.10) is preserved. Such a stability result is closely related to a CFL condition as shown next.

Lemma 4. Let E n+ 1 2
H,h be defined as follows:

E n+ 1 2 H,h := 1 2 ε e n+ 1 2 H,h , e n+ 1 2 H,h T H + 1 2 Ä µ h n H,h , h n+1 H,h ä T H
for n = 1, . . . , N, (5.11) where (e

n+ 1 2 H,h , h n H,h
) is the solution of (5.9)-(5.10). Then, E

n+ 1 2
H,h is a quadratic form under the CFL condition ∆t ≤ C h (where C is a generic constant that is independant of ∆t and H, and h), and it holds that

E n+ 1 2 H,h = E n-1 2 H,h + 1 2 ∆t(f n , e n+ 1 2 H,h + e n-1 2 H,h ) T H for n = 1, . . . , N. (5.12) Proof. Take e n-1 2 H,h + e n+ 1 2
H,h as a test function in the second equation of (5.10) and h n H,h in the first equation of (5.10) evaluated at time-steps n and n + 1, to get

                                           Ñ ε e n+ 1 2 H,h -e n-1 2 H,h ∆t , e n+ 1 2 H,h + e n-1 2 H,h é T K h = f n , e n+ 1 2 H,h + e n-1 2 H,h T K h + h n H,h , ∇ × (e n+ 1 2 H,h + e n-1 2 H,h ) T K h -h n H,h , e n+ 1 2 H,h + e n-1 2 H,h F K 0 -(λ n H , e n+ 1 2 H,h + e n-1 2 H
,h ) ∂K , Ç µ h n H,h -h n-1 H,h ∆t , h n H,h å T K h = -e n-1 2 H,h , ∇ × h n H,h T K h + e n- 1 2 H,h , h n H,h F K h
, Ç µ h n+1 H,h -h n H,h ∆t , h n H,h å T K h = -e n+ 1 2 H,h , ∇ × h n H,h T K h + e n+ 1 2 H,h , h n H,h F K h . (5.13)
Now, summing up the above equations over K ∈ T H , and using the properties of fluxes (see [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF][START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] for these standard arguments), we arrive at

K∈T H (ε e n+ 1 2 H,h , e n+ 1 2 H,h ) T K h -(ε e n-1 2 H,h , e n-1 2 H,h ) T K h +(µ h n H,h , h n+1 H,h ) T K h -(µ h n-1 H,h , h n H,h ) T K h = ∆t K∈T H ï f n , e n+ 1 2 H,h + e n-1 2 H,h T K h -∆t λ n H , e n+ 1 2 H,h + e n-1 2 H,h F K h ò .
From (5.9) and using definition (5.11), it holds that

E n+ 1 2 H,h -E n-1 2 H,h = 1 2 ∆t(f n , e n+ 1 2 H,h + e n-1 2 H,h ) T H . (5.14)
Furthermore, under the CFL condition ∆t ≤ C h (see [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF]) the fully discrete energy

E n+ 1 2
H,h is a quadratic form in terms of the unknowns. Indeed, the term

1 2 (µ h n H,h , h n+1 H,h ) T H can be rewritten as 1 2 (µ h n H,h , h n H,h ) T H + 1 2 (µ h n H,h , h n+1 H,h -h n H,h ) T H .
Using the third equation in (5.13), and integrating it by parts, it holds that

1 2 (µ h n H,h , h n+1 H,h -h n H,h ) T H = ∆t 2 K∈T H ∇ × e n+ 1 2 H,h , h n H,h T K h + e n+ 1 2 H,h , h n H,h F K h . (5.15)
The result follows using inverse inequalities for both terms in the right-hand side of (5.15) and absorbing the negative contributions of the L 2 -norm of h n H,h and e

n+ 1 2 H,h
into the corresponding positive terms that defines the energy function (under the CFL condition). Now, in preparation for the practical algorithm induced by the MHM method (5.9)-(5.10), we remark that the system (5.9)-(5.10) may be decoupled once we bring out their dependency on the solutions from the previous time-step. Indeed, owing to the linearity of problem (5.10), the electric field e where the functions (e

λ,n+ 1 2 H,h , e f ,n+ 1 2 h ) satisfy              1 ∆t ε e λ,n+ 1 2 H,h , v h K = -(λ n H , v h ) ∂K , 1 ∆t ε e f ,n+ 1 2 h , v h K = (f n , v h ) K + h n H,h , ∇ × v h T K h -h n H,h , v h F K 0 + 1 ∆t ε e n-1 2 H,h , v h K , (5.17) and we chose (e λ, 1 2 H,h , e f , 1 2
h ) as a (convergent) approximation of their continuous counterpart. Regarding h n H,h , it appears (again) from (5.10) using the following postprocessing:

1 ∆t µ h n H,h , w h K = 1 ∆t µ h n-1 H,h , w h K -e n-1 2 H,h , ∇ × w h T K h + e n-1 2 H,h , w h F K h , (5.18)
where h 0 H,h = h 0 . Finally, substituting (5.16)

) ∂T H = -(ν H , e f ,n+ 1 2 h
) ∂T H for all ν H ∈ Λ H . (5.19) 6. Numerical results. This section starts with a description of the algorithm underlying the two-level MHM method presented in section 5. Then, an extensive validation of the MHM method through two numerical tests in two-dimensional domains (see the TM model in Appendix B for details) is presented. Theoretical aspects of the MHM method are highlighted in the first test, for which an analytical solution is available. We focus on convergence aspects of the MHM method, conservation properties, and comparison with the standard DG method proposed in [START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF]. The second test assesses the capacity of the MHM method to deal with heterogeneous media on coarse meshes and with interfaces which are not aligned with the faces of the partition. In what follows, we assume the nondimensional version of the Maxwell equations.

6.1. The staggered MHM algorithm. We outline the steps involved in the practical algorithm behind the MHM method (5.17)- (5.19). To this end, we need some additional notation. Define the integer numbers

m H := dim Λ H and m K h := dim V h (K). Let {ψ 1 , ψ 2 , . . . , ψ m H } and {φ K 1 , φ K 2 , . . . , φ K m K h } be
λ n H = m H i=1 β n i ψ i and η e i,h | K = m K h l=1 c K i,l φ K l and e f ,n+1/2 h | K = m K h l=1 γ n+1/2,K l φ K l .
As a result of (5.17), the functions η e i,h | K satisfy

ε η e i,h , φ K k K = -∆t (ψ i , φ K k ) ∂K for all φ K k ∈ V h (K). (6.1)
The problem (5.19) is equivalent to the following positive definite linear system:

A β n = F n+1/2 (6.2)
for all n ∈ {1, . . . , N }, where (

β n ) T = (β n 1 , β n 2 , . . . , β n m H ) ∈ R m H . Note that A is an m H × m H -matrix with entries given by A = (a ij ) ∈ R m H ×m H with a ij := -(ψ i , η e j,h ) ∂T H , 1 ≤ i, j ≤ m H ,
and F n+1/2 is an m H -vector whose entries are

F n+1/2 = (f i ) ∈ R m H with f i := (ψ i , e f ,n+1/2 h ) ∂T H , 1 ≤ i ≤ m H .
In what follows, we denote by

C j,K , H n-1 K , E n-1 2 K , E f ,n+ 1 2 K
, and F n K the local vectors composed by the degrees of freedom of η e j,h

| K , h n-1 H,h | K , e n-1 2 H,h | K , e f ,n+ 1 2 h
| K , and f n h | K , respectively. They are computed using the following MHM algorithm. Remark 9. Since step 5 in the algorithm is independent of the time-step n, the vector C j,K can be computed once and for all at n = 1. Thanks to the structure of the DG method, the matrices involved in steps 5-7 are block diagonal. As such, they can be factored ahead of time resulting in explicit local solvers. Also interesting is that the local matrices used to build up the global system in step 10 have already been obtained in the previous steps. Thereby, the global matrices involved in step 10 can be assembled "on the fly." Particularly, the "action" of the inverse of the block diagonal matrix M ε has already been taken into account in step 5 through completely independent local problems. This feature makes the MHM method suitable to be implemented in parallel architectures. Finally, we observe that the left-hand side of the linear system in step 10 (Schur complement) corresponds to a symmetric positive definite matrix.

6.

2. An analytical solution. We now assess the theoretical aspects of the MHM method presented in the previous sections for the two-dimensional Maxwell equations (see Appendix B). To this end, we consider the analytical solution given by H,h and h 0 H,h from the initial conditions e 0 and h 0 2: for n = 1 to N do 3:

h x (t, x, y) = - √ 2
for K ∈ T H do 4:

Compute the local matrices, with k, l ∈ {1, . . . , m K h } and j ∈ {1, . . . , m H }

M δ K := (m K lk ), m K lk := δ φ K l , φ K k T K h , δ = 1, ε, µ and
K K := (k K kl ), k K kl = φ K l , ∇ × φ K k T K h -( ¶ φ K l © , φ K k ) F K h and S j,K := (s K jk ), s K jk = (ψ j , φ K k ) ∂K 5: Compute C j,K = {c K j,l } from M ε K C j,K = -∆t S j,K 6
:

Compute H n K = {d n,K l } from M µ K H n K = M µ K H n-1 K -∆t K K E n-1 2 K and update h n H,h h n H,h | K = m K h l=1 d n,K l φ K l 7: Compute E f ,n+ 1 2 K = {γ n+1/2,K l } from M ε K E f ,n+ 1 2 K = M ε K E n-1 2 K + ∆t K K H n K + ∆t M 1 K F n K 8: Assemble M ε K , C j,K , S j,K and E f ,n+ 1 2 K into M ε , C, S and E f ,n+ 1 2 9:
end for 10:

Solve the global system

S T (M ε ) -1 S C β n = - 1 ∆t S T E f ,n+ 1 2 11: Update e n+ 1 2 H,h e n+ 1 2 H,h = m H j=1 β n j K∈T H m K h l=1 c K j,l φ K l + K∈T H m K h l=1 γ n+1/2,K l φ K l 12: end for
From this choice, we prescribe the dielectric boundary condition e z (t, x, y) = 0 on ∂Ω, initial conditions h x (0, x, y) = h y (0, x, y) = 0 and e z (0, x, y) = sin(2 π x) sin(2 π y) in Ω, and the right-hand side j z (t, x, y) = 0 in Ω. The MHM method (5.19) is solved with linear (Λ H = Λ 1 ) and quadratic (Λ H = Λ 2 ) polynomial interpolation on faces. Regarding the second level, we adopt a one-element submesh with a P l+2 (K) interpolation. The convergence is measured using the following discrete norms:

(h -h H , e -e H ) ∞ := max n=0,...,N -1 h x,n -h x,n H,h 2
0,Ω + h y,n -h y,n H,h 2
0,Ω + e z,n -e z,n H,h 2
0,Ω 1/2 , (h -h H , e -e H ) curl,∞ := max n=0,...,N -1 h x,n -h x,n H,h 2
W + h y,n -h y,n H,h 2
W + e z,n -e z,n H,h 2 W 1/2
.

First, we look at the linear interpolation case on faces (Λ H = Λ 1 ). In Figures . Interestingly, we found that the error associated to the electric field is superconvergent (see Figure 6.4 on the right). Next, we compare the MHM method using linear interpolation on faces (Λ H = Λ 1 ) with the DG method [START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] using quadratic interpolation. To this end, we measure the error in the • ∞ and the • curl,∞ norms. We observe in Figure 6.5 that the MHM method outperforms the DG method. Indeed, the MHM method is optimally convergent in both norms, whereas the DG method is suboptimal (e.g., nonconvergent) in the • ∞ norm (e.g., • curl,∞ norm) in some cases. It is worth recalling that the MHM method achieves such an efficiency despite using the DG method as a second-level solver, which validates the robustness of the approach.

Also, we verify in Figure 6.6 (left) that the discrete energy is fully preserved in accordance with the theory. The local magnetic field conservation is investigated in Figure 6.6 (right) although the method has not imposed such a constraint. Surprisingly, we observe that max K∈T H | K ∇ • h H,h dx| stays bounded and close to zero as time increases.

6.3.

A nano-waveguide problem. This is a prototype problem of a photonic crystal structure in the emerging nano-photonics area [START_REF] Joannopoulos | Photonics crystals: Molding the Flow of Light[END_REF]. The idea here is to simulate an idealized waveguide device before it has undertaken an optimization procedure to maximize the performance of directional transmission (see [START_REF] Santos | Designing novel photonic devices by bioinspired computing[END_REF], for example). The periodic distribution and sizes of the photonic crystal structure are set for simplicity as the MHM method relies on neither periodicity nor scale separation assumptions. The photonic-crystal type represents a nano-structuring device encapsulated in a 10 × 10 square which is composed of 15 cylindrical holes of radius 0.3125 (see Figure 6.7 for details). Here µ = 1 and the value of the permittivity ε corresponds to silicium within the holes (ε = 3.14), silica (ε = 1.5) in the device enclosing the holes, and surrounding air (ε = 1). We use absorbing boundary conditions with an incident plane wave of frequency 0.5 (see Remark 3).

We solve the problem for a time period equal to T = 16/ √ 2 and set ∆t = 10 -2 . First, we adopt the DG method [START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] on top of a quadrilateral mesh composed of 65,536 elements using quadratic interpolations to compute a reference solution. This leads to 589,824 degrees of freedom. We set a first coarse mesh with 256 quadrilateral elements (see Figure 6.8 on the right). On top of this mesh, we use the MHM method with Λ H = Λ 8 1 and arrive at 8,704 degrees of freedom total. Submeshes account for 64 elements with quadratic interpolation. As such, the total number of degrees of freedom associated with the approximation of the basis functions is 147,456. We recall that such basis computations are completely independent of one another. A comparison between the reference electric field computed on the fine mesh and the one obtained from the MHM method with its overlying mesh is depicted in Figure 6.9. Observe that both electric fields coincide and that the MHM method is able to precisely account for the impact of crossing-face interfaces on the solutions. Next, using a much coarser mesh of 64 elements (see Figure 6.9 on the left) with Λ H = Λ 16 1 , we observe in Figure 6.10 that the solution remains fairly close to the reference. Although this corresponds to only 4,608 degrees of freedom, the approximate electric field shows good agreement with the solution (see Figure 6.10 on the right). Also in Figure 6.10 (left), we .

present the approximate electric field using the DG method [START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] on a mesh of 1,024 elements with quadratic interpolation (corresponding to 9,216 degrees of freedom).

As expected, the solution from the DG method is less precise since the mesh is notaligned with the material interfaces. In conclusion, the MHM method provides an accurate framework to simulate wave propagation phenomena on coarse meshes when faces are nonaligned with interfaces. The upshot is that the shortcomings found in standard finite element methods are overcome when they are used as two-level solvers within the MHM method.

Using a reference solution computed on a fine mesh with 1,048,576 P 2 (K) elements (9,437,184 degrees of freedom total), we calculate the error in the L 2 norm at the final time-step (denoted by • 0,Ω ) between such a reference and the solutions from the DG and MHM methods used to generate Proof. Denote by W * the space of continuous linear forms on W for the norm induced by the scalar product (•, •) W . Since L belongs to W * , then from the Riesz-Fréchet theorem, it holds that there exists f ∈ W such that L(v) = (f , v) W . This can be rewritten as

L(v) = K∈T H (f , v) K + K∈T H (p, ∇ × v) K (A.2)
with f ∈ L 2 (Ω), and p = ∇ × f ∈ L 2 (K) for all K ∈ T H . Since L vanishes on H 0 (curl; Ω), we get As a result, we have f = -∇×p, which implies that ∇×p ∈ L 2 (Ω) (in a distributional sense) and then p ∈ H(curl; Ω). Also, from (A.2) and the characterization of f , it holds that

L(v) = - K∈T H (∇ × p, v) K + K∈T H (p, ∇ × v) K , (A.3)
and then using Green's formula, and p ∈ H(curl; Ω), we arrive at

L(v) = K∈T H (p × n K , v) ∂K for all v ∈ V. (A.4)
We conclude there exists ν ∈ Λ (e.g., ν := p × n K on ∂K for all K ∈ T H ) such that

L(v) = (ν, v) ∂T H for all v ∈ V. (A.5)
Let us prove that ν is unique. Suppose that (ν, v) ∂T H = 0 for all v ∈ V. (A.6) From the surjectivity of the trace application from H 1 (K) to H

(3. 11 ) 4 .

 114 Remark If an absorbing boundary condition is prescribed on ∂Ω A ⊂ ∂Ω, then from Remark 3, (3.11) becomes (ν, e n,λ

2 Fig. 6 .

 26 Fig. 6.1. Optimal convergence history with respect to H (left) and ∆t (right) in (h -h H , ee H ) ∞ norm. Here Λ H = Λ 1 .

Fig. 6 . 6 .

 66 Fig. 6.6. History of the discrete energy conservation (left) and history of the local magnetic field conservation (right). Here Λ H = Λ 1 .

Fig. 6 . 7 .

 67 Fig. 6.7. The photonic crystal structure.

Fig. 6 . 8 .

 68 Fig. 6.8. Coarse meshes.

Fig. 6 . 9 .

 69 Fig.6.9. The electric field using the DG method with 589,824 degrees of freedom (left) and using the MHM method with 156,160 degrees of freedom (right) at time t =16

√ 2 .Fig. 6 . 10 .

 2610 Fig. 6.10. Comparison between the electric field on a mesh of 9,216 degrees of freedom from the DG method (left) and on a mesh of 4,608 degrees of freedom from the MHM method (right). Here t = 16 √ 2

  Figure 6.10. Notably, we found that the MHM method yields a drastic decrease in error which amounts to a factor of order 10b(ν, v) := (ν, v) ∂T H .Furthermore, b can be extended in a unique way to Λ × W such thatL(v) = b(ν, v) for all v ∈ W,where b is a continuous bilinear form on Λ × W that uniquely extends b. Moreover,H 0 (curl; Ω) = ¶ v ∈ W : b(ν, v) = 0 for all ν ∈ Λ © . (A.1)

  v) K + K∈T H (p, ∇ × v) K = 0 for all v ∈ H 0 (curl; Ω),and then it holds that (f , v) Ω + (p, ∇ × v) Ω = 0 for all v ∈ H 0 (curl; Ω).

1 2 (

 2 ∂K), we deduce that ν = 0 in H -1 2 (∂K) for all K ∈ T H , and then uniqueness follows. Now, since V is dense in W and b(ν, •) is uniformly continuous on V with respect to the norm • W , b(ν, •) can be extended in a unique manner to W by density. Let us denote by b(ν, •) such an extension for all given ν ∈ Λ. By the uniqueness of this

  a basis for Λ H and V h (K), ∈ {1, . . . , N } and i ∈ {1, . . . , m H }, such that

	γ	n+1/2,K 1	, γ 2 n+1/2,K	, . . . , γ	n+1/2,K h m K	, with n
	respectively. There exist real numbers β n 1 , β n 2 , . . . , β n m H and c K i,1 , c K i,2 , . . . , c K i,m K h	, and

  [START_REF] Arbogast | Analysis of a two-scale locally conservative subgrid upscaling for elliptic problems[END_REF].1 and 6.2, we verify that the MHM method achieves error optimality with respect to H and ∆t. Specifically, we find that the numerical solutions converge as O(H 2 ) and O(∆t 2 ) in the • ∞ -norm and O(H) in the • curl,∞ -norm. We perform the same analysis with Λ H = Λ 2 . The expected optimal convergence rates are again recovered (e.g., O(H 3 ) and O(∆t 2 ), and O(H 2 ) in the • ∞ and • curl,∞ , respectively). This is depicted in Figures 6.3and 6.4 (on the left)

  1. Optimal convergence history with respect to H (left) and ∆t (right) in (h -h H , ee H ) ∞ norm. Here Λ H = Λ 1 . Optimal convergence history with respect to H in the (h -h H , ee H ) curl,∞ norm. Here Λ H = Λ 1 . Fig. 6.3. Optimal convergence history with respect to H (left) and ∆t (right) in the (h -h H , ee H ) ∞ norm. Here Λ H = Λ 2 . Optimal convergence history with respect to H in the (h -h H , ee H ) curl,∞ norm with Λ H = Λ 2 (left). Superconvergence with respect to H for the ee H curl,∞ norm (right) using Λ H = Λ 1 and Λ H = Λ 2 . Comparison between the MHM method and the DG method with respect to the (hh H , ee H ) ∞ (left) and (h -h H , ee H ) curl,∞ (right) norms.

		-1 10								4.096
	h , eH -||∞ e) H -||( h	-2 10 10 -3 10 -4 -5 10	||(hH -e)||curl,∞ h, eH -	0.0125 0.025 0.05		DG, k=1		∞ , ||curl , eH -e) H -h ||( h	1.024 0.256 0.064 0.016 0.004		DG, k=1
					0.0063		MHM, l=1			0.001	MHM	MHM, l=1
		10 -6					DG, k=2 MHM, l=2			0.0003	H	DG, k=2 MHM, l=2
	||( h H -h , eH -e) ||∞ ||(hH -h, eH -e)||curl,∞ E n H,h	Fig. 6.2. 0.0156 10 -6 10 -5 10 -4 0.0078 0.00005 0.0001 0.0002 0.0004 0.0008 0.0016 0.0032 Fig. 6.4. 0.0156 0 0.5 0.1249 0.1249 0.1249 0.1250 0.1250 0.125 0.1250 0.1250 0.1251 0.1251 Fig. 6.5. t 0.0312 H 0.0156 0.0312 H 1 1.5 2 0.1251	0.0078 0.0313 0.0625 0.0625 H 2.5	3	0.0156 MHM H 3 0.0625 0.125 MHM H 2 0.125 3.5 4	||( h H -h , eH -e) ||∞ ||eH -e||curl , ∞	0.0313 0.0156 0.0156 0.0003 0 0.5 0 0.2 -6 10 H 10 -5 10 -4 10 -5 10 -4 10 -3 10 -2 0.4 0.6 0.8 1 1.2 maxK K ∇ • h dx 1.4 ×10 -13	0.0625 0.0312 0.0007 0.0312 1 1.5	∆t H H 2 t	0.0014 0.0625 0.0625 2.5	3	0.0028 ∆ t 2 MHM 0.125 MHM, l=1 MHM, l=2 H 2 H 3 0.125 3.5 4

The research leading to these results has received funding from the European Union's Horizon 2020 Programme (2014-2020) and from the Brazilian Ministry of Science, Technology and Innovation through Rede Nacional de Pesquisa (RNP) under the HPC4E project (www.hpc4e.eu), grant agreement 689772. The second author was funded by CONICYT/Chile through FONDECYT project 1181572 and CNPq/Brazil. The fourth author was funded by CNPq/Brazil, CAPES/Brazil, and INRIA/France.

.

compared to the DG method. Observe that such a result is in complete accordance with Figure 6.10. Also, we study the convergence of the numerical solution from the MHM method using Λ H = Λ 1 . Since the exact solution is not regular due to the jumping coefficients, we would expect a decrease in the rate of convergence. Indeed, we observe in Figure 6.11 that the MHM solution achieves an error of order O(H 3/2 ) in the • 0,Ω norm.

Conclusion.

This work introduced the foundation of the MHM method for the Maxwell equations in the time domain with heterogeneous coefficients. It relied on independent local Maxwell problems to devise basis functions which are brought together by a one-field face-based global problem. As a result, the approximate electric and magnetic fields are driven by face-based degrees of freedom as they are obtained from a simple processing of the Lagrange multipliers. A DG method combined with a leap-frog scheme was used as a second-level solver to approximate the local Maxwell problems, which made the MHM method practical. With such a choice, the present method is shown to be optimally convergent and to achieve high-order accuracy in space (and second-order in time), with the upshot of preserving a discrete energy principle. In addition, the approximated electric field was found to be superconvergent. The MHM method naturally embeds an upscaling procedure which turns out to be suitable to simulate wave propagation problems in highly heterogenous media on coarse meshes. We conclude that the MHM method, which is particularly adapted to parallel computing environments, emerges as a realistic option to handle threedimensional wave propagation problems in highly heterogeneous media as found in the nano-photonic field. Nevertheless, some open questions still remain, such as the cost-effectiveness of the MHM method on modern parallel environments when compared to other domain decomposition approaches, and its numerical analysis. These important subjects will be addressed in forthcoming works.

Appendix A. This section outlines the proof of the main lemma used throughout the manuscript. Lemma 5. Let L be a linear continuous form on W. Hence, L vanishes on H 0 (curl; Ω) if and only if there exists a unique ν ∈ Λ such that

where b is the continuous bilinear form on Λ × V defined by extension, b(•, •) is also a bilinear form, which is also continuous on Λ × W as it is an extension (by density) of a continuous bilinear form. Furthermore, since L and b are continuous on Λ × W and coincide on a dense subspace, we get

Let B ∈ L(W, Λ ) be the linear operator associated to b, and let B ∈ L(Λ, W ) be its adjoint. We recall that the polar space of a set X ⊂ W is defined by

We proved above that B is an isomorphism from Λ onto H 0 (curl; Ω) • , and then the range of B coincides with H 0 (curl; Ω) • . Denoting by Ker( B) the null-space of operator B, and using that B is a continuous operator between Banach spaces, we get (see [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF], for instance)

where we used that H 0 (curl; Ω) is closed in W with respect to • W , and then (A.1) follows.

We are ready to prove Lemma 2.

Proof of Lemma 2. Assume that (e, h, λ) ∈ C 0 (0, T ; W)×C 0 (0, T ; W)×C 0 (0, T ; Λ) solves (2.5), and take I n ∈ T ∆t . Using Lemma 5 and the third equation in (2.5), the electric field e n ∈ C 0 (I n ; H 0 (curl; Ω)). Next, we choose v and w in H 0 (curl; Ω) in (2.5), and integrate the second equation in (2.5) by parts. Using e n ∈ C 0 (I n ; H 0 (curl; Ω)), the second equation in (2.3) holds. Now, integrating the first equation in (2.5) by parts yields

and, thus, the first equation in (2.3) is satisfied in the L 2 sense in each K ∈ T H and for all t ∈ I n . Observe that (A.9) implies that λ n = -n K × h n | ∂K for all K in T H and I n ∈ T ∆t , and then h n ∈ C 0 (I n ; H(curl; Ω)). Now, collecting the previous results, we conclude problem (2.3) holds for all I n ∈ T ∆t , and

The fields also belong to C 1 (0, T ; L 2 (Ω)) from (2.5) and the regularity result in (A.12), namely,

We next verify that the divergence constraints in (1.1) are satisfied in [0, T ]. We start proving the result in the time interval I 1 := [t 0 , t 1 ] by testing problem (2.5) with ∇β, ∇η, where (β, η) are infinitely differentiable functions with compact support in Ω. As a result, for all t ∈ I 1 , it holds that

Using that ∇β ∈ H 0 (curl; Ω), it turns out that (λ 1 , ∇β) ∂T H = 0 since λ(t) ∈ Λ for all t ∈ I 1 . It results that (in a distributional sense)

and then from (1.2) and (2.1), we get

Following the same reasoning, we prove that if

and the result follows. Now, assume that (e, ∂ t e) ∈ C 0 (0, T ; H 0 (curl; Ω)∩H(div, ε; Ω))×C 0 (0, T ; L 2 (Ω)) and (h, ∂ t h) ∈ C 0 (0, T ; H(curl; Ω) ∩ H 0 (div, µ; Ω)) × C 0 (0, T ; L 2 (Ω)) solves (2.3). We define for each t ∈ (0, T ), the following continuous linear form:

Noticing that L t vanishes on H 0 (curl; Ω), from Lemma 5 there exists a unique λ(t) ∈ Λ such that

Here (•, •) ∂T H stands for the unique density extension to Λ × W of (•, •) ∂T H defined on Λ × V (corresponding to the operator b(•, •) in Lemma 5). Thereby, (e, h, λ) satisfies the first equation in (2.5) on each I n ∈ T ∆t and respects the continuity condition at t n for each n ∈ {0, N -1}. Also, owing to the characterization of H 0 (curl; Ω) in (A.1), the function e satisfies the third equation in (2.5). Also, the second equation in (2.3) implies immediately that the second equation in (2.5) holds. Next, from the definition of L t in (A. [START_REF] Costabel | Singularities of Maxwell interface problems[END_REF], we arrive at

and then we get

where (•, •) ∂K is the extension to Λ × W of the natural scalar product on ∂K defined on Λ × V. On the other hand, we have that (A.20) holds, and then

which leads to λ = -n K × h on ∂K for all K ∈ T H . As a result, the assumed regularity for h implies the claimed regularity for λ (and thus for λ n ).

Next, we address the well-posedness of the local Maxwell problems.

Lemma 6. Assume ε and µ satisfy (1.4) and (i)

), for all K ∈ T H and n ∈ {1, . . . , N }. Hence, for all I n ∈ T ∆t , the following initial boundary value problems are well-posed: find (e n,λ , h n,λ ) such that

and (e n,f , h n,f ) such that

with (e 0 , h 0 )(0, •) := (e 0 , h 0 ).

Proof. We first address the case of ε and µ are constant functions in K ∈ T H . Let K ∈ T H , n ∈ {1, . . . , N }. Hereafter, we skip the superscripts for the sake of clarity Observe that (A.22) and (A.23) can be written in the following second-order form in D (I n ): [START_REF] Efendiev | Multiscale finite element problems with highly oscillatory coefficients[END_REF] where (u n-1 , v n-1 ) ∈ V(K) × L 2 (K), and

Here, u n-1 and v n-1 stand for the initial conditions of (A.22) or (A.23) and

for system (A.22), or

for system (A.23). Such a derivation is standard and can be found in [8, Chapter 7], for instance. For (A.27), the result follows from [20, Theorem 1, p. 558]. Regarding (A.26), we cannot directly apply such a theorem since its right-hand side is not continuous in L 2 (I n , L 2 (K)) but rather in L 2 (I n , H(curl; K)). We then revisit the proof in [20, Theorem 1] and remark that the well-posedness result for (A.26) holds if ∂ t u and u are bounded in L ∞ (I n , L 2 (K)) and L ∞ (I n , H(curl; K)), respectively. Thereby, we establish next such a priori estimates for the system (A.24) with (A.26).

Consider two positive constants β and α such that

and then integrating (A.29) over [t n-1 , t], t ∈ I n , we get

By integration by parts, we rewrite the first term in the right-hand side of (A.30) as

and then

From a weighted Young's inequality with positive constants δ and δ, it holds that t tn-1

where

Using (A.28) and a weighted Young's inequality again with a positive constant δ, we arrive at

Next, choosing δ = 2 α one obtains that there exists a positive constant M such that

and the a priori estimate follows from the Grönwall inequality. The well-posedness result arises following closely the proof in [20, Theorem 1]. The proof generalizes to the case functions ε and µ satisfy (1.4). It follows by replacing (A.24) by

where

and following closely the proof of the constant case using (1.4).

Appendix B. For sake of completeness, the MHM method is presented for the two-dimensional Maxwell equations (the TM model) in this section. First, we recall that the TM problem consists of finding the electric field e z : (0, T ) × Ω → R and the magnetic fields

where j z is the electric current density, and e 0 , h x 0 and h y 0 are given regular functions with values in R 2 .

Here, the space for the Lagrange multipliers reduces to

where s stands for the tangential unit vector on ∂K. Next, denote by V the space of functions in L 2 (Ω) such that their restriction to K belongs to H 1 (K) for all K ∈ T H . We skip the intermediate steps leading to the fully discrete version of the MHM method as they follow closely the three-dimensional case. The notation adopted below is that used in the three-dimensional case, if not mentioned otherwise. We select Λ H ⊂ Λ as the polynomial space

where we recall that P m l (F ) is the space of discontinuous polynomial functions on F of degree less than or equal to l ≥ 0, and defined on an equally spaced partition of F composed of m elements (m ≥ 1). The two-dimensional version the MHM method using the leap-frog time scheme and the DG method proposed in [START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] adopts the nonconforming finite dimensional spaces

The fully discrete MHM method reads as follows: Find (e

H,h = h x 0 , and h y,0 H,h = h y 0 . The exact electric and magnetic fields obtained from (B.1) and evaluated at time-steps n + 1/2 and n, respectively, are approximated as follows: ) satisfy, in each K ∈ T H , the following local problems, respectively: