{'fraction_non_alphanumeric': 0.03662407197858306, 'fraction_numerical': 0.011453647227665387, 'mean_word_length': 4.28886293971683, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FDTD schemes for Maxwell's equations with embedded perfect electric conductors based on the correction function method

In this work, we propose staggered FDTD schemes based on the correction function method (CFM) to discretize Maxwell's equations with embedded perfect electric conductor (PEC) boundary conditions. The CFM uses a minimization procedure to compute a correction to a given FD scheme in the vicinity of the embedded boundary to retain its order. The minimization problem associated with CFM approaches is analyzed in the context of Maxwell's equations with embedded boundaries. In order to obtain a well-posed problem, we propose fictitious interface conditions to fulfill the lack of information, namely the surface current and charge density, on the embedded boundary. Fictitious interfaces can induce some issues for long time simulations and therefore the penalization coefficient associated with fictitious interface conditions must be chosen small enough. We introduce CFM-FDTD schemes based on the well-known Yee scheme and a fourth-order staggered FDTD scheme. Long time simulations and convergence studies are performed in 2-D for various geometries of the embedded boundary. CFM-FDTD schemes have shown high-order convergence.


Introduction
In electromagnetic dynamics, the perfect electric conductor is an important idealized material, that allows surface charges and currents. Since the surface charge and current density are often unknown, PEC walls are modeled through the imposition of boundary conditions that require the continuity of the normal component of the magnetic field and the tangential component of the electric field on the boundary of a domain.
In free-space simulations involving PECs, the interface between a PEC and its surrounding medium is treated as an embedded boundary. Embedded boundary conditions can be difficult to treat, particularly in a finite-difference context. Indeed, challenges include the development of numerical methods that can handle various complex geometries of the PEC without increasing the complexity of a numerical method while retaining high-order accuracy. It is also worth mentioning that high-order schemes are important to diminish the phase error for long time simulations [1].
Many numerical strategies are proposed to achieve high-order accuracy for problems involving embedded PEC walls, such as discontinuous Galerkin (DG) approaches [2], pseudospectral time-domain methods (PSTD) [3,4,5] or finitedifference time-domain (FDTD) schemes [6,7,8,9,10]. A discontinuous Galerkin approach can treat complex geometries of embedded PEC walls by an appropriate mesh. However, a large number of unknowns for high-order accuracy is needed for these approaches. This is due to the use of piecewise polynomial spaces that do not require continuity between two elements of a mesh. Various strategies have been proposed to reduce the computational cost of DG based methods, such as parallel computing strategies or particular choices of basis functions [11,12]. It is also worth mentioning that finite-element approaches with non-body-fitted grids have been developed for electromagnetic problems, but low-order basis functions have been used [13,14,15]. On the other hand, finite-difference time-domain approaches use a simple Cartesian grid and have low computational costs. However, the imposition of embedded boundary con-ditions is far from trivial. A naive approach is to use the Yee scheme [16] with a staircased approximation of the embedded boundary. Unfortunately, this approach leads to a first-order scheme at best and sometimes to non-convergent approximations [8]. To overcome this issue, many FDTD approaches have been proposed, such as overlapping grids [6], contour path methods [7,10] and locally modifications of FD schemes [8,9]. A staircase-free FDTD scheme has been proposed in [8] to recover a second-order scheme without significantly compromising the simplicity of the used FD scheme. They explicitly impose PEC boundary conditions by locally modifying a finite-difference scheme in the vicinity of the boundary. Following the same idea, a fourth-order finite-difference scheme based on the Matched Interface and Boundary (MIB) method [17] has been proposed to handle embedded PEC walls using the vector Helmholtz equation [9]. This FD scheme is obtained by deriving and explicitly imposing jump conditions for PEC walls on the embedded boundary. However, the complexity of a MIB based scheme increases with its order or the complexity of the geometry of the embedded boundary because of the imposition of high-order jump conditions [18,19].
Finally, pseudospectral methods have the advantage to require less grid points per wavelength than FDTD approaches. The Fourier and Chebyshev collocation methods with a multidomain strategy have been used to achieve high-order accuracy [3,4]. These approaches need a multidomain decomposition with an appropriate mesh grid for embedded boundary conditions, and therefore an additional treatment of interfaces between subdomains is needed. An alternative approach is to use a Fourier penalty method [5]. Although this approach does not need a multidomain strategy, there is some stability issues that limit the order of the method in two and three dimensions.
Another avenue to handle embedded boundary conditions is FD schemes based on the Correction Function Method (CFM). The CFM, which was inspired by the Ghost Fluid Method (GFM), has been originally developed to treat Poisson's equations with interface jump conditions with arbitrarily complex interfaces. FD schemes based on the CFM achieve high-order by means of a minimization problem. This numerical strategy has been successfully applied to Poisson problems with constant and piecewise constant coefficients, and interface jumps [20,21]. It is also worth mentioning that a CFM based strategy have also been used to treat 3-D Poisson equation with interface jump conditions [22]. Afterward, extensions of this method have been used to handle the wave equation and Maxwell's equations with constant coefficients and interface jump conditions [23,24]. Briefly, the underline assumption of the CFM is that jumps on the interface can be smoothly extended in the vicinity of the interface. Therefore, a system of partial differential equations (PDEs) coming from the original system of PDEs is derived to model jumps around the interface. The solution of this system of PDEs is called the correction function. A square measure of the error associated with the correction function's system of PDEs is then minimized to compute approximations of the correction function.
Afterward, these approximations are used to correct a given FD scheme that involves nodes in different subdomains. It is worth mentioning that the additional computational cost associated with minimization problems of the CFM is not negligible. Hence, this makes unavoidable the use of parallel computing strategies [25].
As mentioned before, a FDTD strategy based on the CFM has been developed to handle Maxwell's equations with constant coefficients and interface jump conditions [24]. Even though this numerical scheme could also be used to enforce embedded boundary conditions, one needs to impose all information on the boundary, that is both normal and tangential components of each electromagnetic field. This is a major drawback when embedded PEC wall boundary conditions are considered. As shown in this work, a direct application of the numerical scheme proposed in [24] for PEC wall boundary conditions leads to an ill-posed minimization problem because of the lack of information on the embedded boundary. Hence, the main goal of this work is to overcome this issue and proposes high-order FDTD schemes based on the CFM, which are referred as CFM-FDTD schemes, to handle embedded PEC wall boundary conditions. The algorithm presented in this paper also provides an important stepping stone towards a CFM-FDTD approach to handle Maxwell's interface problems.
We first extend CFM-FDTD schemes for PEC problems for which the surface current and charge density are unknown. Afterward, we describe a construction of appropriate local patches that are needed for minimization problems. This construction reduces the computation cost of the CFM while guarantees the uniqueness of the correction function for a given node to be corrected. We then introduce CFM-FDTD schemes based on the well-known Yee scheme and a fourth-order staggered FDTD scheme. Finally, numerical examples based on the transverse magnetic (TM z ) mode are performed to verify the proposed CFM-FDTD schemes.
The paper is organized as follows. In Section 2, we introduce Maxwell's equations with embedded PEC wall conditions. The CFM applied to Maxwell's equations is presented in Section 3. In this section, we propose an extension of the CFM for problems involving a PEC for which the surface current and charge density are unknown. Minimization problems coming from the CFM are analyzed. The impact of the CFM on a FDTD scheme is also investigated and the construction of local patches is described. We then introduce CFM-FDTD schemes based on the Yee scheme in Section 4 and on a fourth-order staggered FDTD scheme in Section 5. Finally, we perform numerical examples to verify the proposed CFM-FDTD schemes in Section 6.

Definition of the Problem
Assume a domain Ω subdivided into two subdomains Ω + and Ω − . The interface Γ between subdomains is independent of time and allows solutions, that is the magnetic field H and the electric field E in this work, to be discontinuous along it. We define H + and E + as the solutions in Ω + , and H − and E − as the solutions in Ω − . The jumps are denoted as We consider the boundary ∂Ω of Ω and a given time interval I = [0, T ]. Assuming linear media and a periodic domain, Maxwell's equations with interface conditions are given by where µ is the magnetic permeability, is the electric permittivity, ρ is the electric charge density, J s is the surface current density, ρ s is the surface charge density andn is the unit normal to the interface Γ pointing toward Ω + . Fig. 1 illustrates a typical geometry of a domain Ω. Without loss of generality, we assume constant coefficients that are such that , µ > 0 and ρ = 0.
For problems involving a PEC, we often do not know the surface current density J s and the surface charge density ρ s . Interface conditions (1e)-(1h) are then reduced ton × E = 0 on Γ × I, Let us assume that subdomain Ω − is a PEC, we then have E − (x, t) = 0, ∀(x, t) ∈ Ω − × I. Considering that the initial condition of the magnetic field is given by conditions on Γ can be considered as embedded boundary conditions, given bŷ for Ω + . In this work, we focus on problems involving PECs and therefore assume Γ to be an embedded boundary of Ω + .

Correction Function Method
The smoothness of solutions is important when one wants to use FD schemes.
Realizing that problem (1) can have discontinuous solutions, standard FD schemes cannot a priori be used around the embedded boundary Γ . The correction function method allows one to circumvent this issue. The purpose of the CFM is to find a correction for a finite difference approximation that involves nodes that belong to different subdomains. To find such a correction, the CFM assumes that solutions in Ω + × I and Ω − × I can be smoothly extended in a small region where Ω Γ ⊂ Ω encloses the embedded boundary Γ , in such a way that the original PDE is still satisfied. A functional that is a square measure of the error of a PDE that describes the behaviour of jumps or correction functions in the vicinity of the embedded boundary is derived. This functional is then minimized in a discrete functional space to obtain approximations of the correction function in Ω Γ × I. In practice, we define a local patch Ω h Γ ⊂ Ω Γ for which the correction function needs to be computed at a node x c ∈ Ω h Γ and a time interval I h Γ = [t n − ∆t Γ , t n ]. The additional computational cost associated with the CFM is not negligible when compared with the original FDTD scheme. In fact, the CFM consumes most of the computational time. However, a parallel implementation of the computation of correction functions can be performed since minimization problems associated with local patches are independent. We refer to [25] for more details about the benefits of a parallel implementation of the CFM.
In the following, we summarize the procedure for the CFM applied on Maxwell's equations when electromagnetic fields are known on the embedded boundary [24]. Afterward, we present an analysis of the minimization problem that is needed for the CFM. The functional to be minimized is then modified and analyzed for embedded PEC walls for which the surface current and charge density are unknown. We also investigate the impact of such a modification on a FDTD scheme. We then describe a construction of local patches that reduces the computation cost of the CFM and ensures an appropriate representation of the embedded boundary within the local patch.
Let us first introduce some notations. The inner product in and we also use the notation for legibility. The correction functions are defined as Let us first assume that the surface current density and charge density are known. Following the same procedure as in [24], one can obtain the following quadratic functional to minimize where c p > 0 is a penalization coefficient and h is the length in space of the patch. We scale the integral over the domain by h to ensure that all terms in the functional J behave in a similar way when the computational grid is refined [24]. To guarantee the divergence-free constraint (1c) and (1d), we minimize the functional J in a divergence-free space-time polynomial spaces, namely where P k denotes the space of polynomials of degree k. It is worth mentioning that basis functions of V are based on high-degree divergence-free basis functions proposed in [11] for discontinuous Galerkin approaches. The problem statement is then where W = V . The following proposition shows that the quadratic functional J has a global minimizer when finite-dimensional functional spaces are used and for an appropriate choice of the penalization coefficient c p . Proof. Let us first notice that problem (2) is an unconstrained minimization problem. In this case, if the hessian matrix coming from a quadratic functional is positive definite, then the critical point is a global minimum. Let us consider v and w to be basis functions of V and W . The hessian matrix coming from problem (2) is given by Requiring the hessian matrix to be positive definite, one finds the following conditions: Since v and w are basis functions of V and W , we have v = 0 and w = 0.
Noticing that v and w cannot be both orthogonal and collinear ton except for the zero element, we also have b = 0 and d = 0. Hence, condition (3) is always satisfied. However, condition (4) leads to the following criterioñ For a sufficiently largec, criterion (5) is satisfied.
Let us now investigate the impact of the CFM based on minimization problem (2) on a given FD scheme. Assume a spatial finite difference operator noted where U is a vector containing m unknowns that estimate electromagnetic fields on the grid points in space. We consider that a correction is needed at r nodes in the vicinity of the embedded boundary, system (6) then becomes where D is a vector containing r values of the correction function coming from problem (2) and A is a rectangular matrix of dimension m × r with either 0 or ±1 as components depending where the correction is needed. Hence, the correction function can therefore be considered as a time-dependent force term.
From Proposition 3.1, we have the existence and unicity of the coefficients of polynomial approximations of the correction function for an appropriate c p .
Since D H (x, t) and D E (x, t) are polynomial functions, and Ω h Γ ×I h Γ is a compact domain, these approximations and their derivatives are bounded in the infinity norm. It is therefore sufficient to investigate the stability of the original FDTD scheme (6) to identify any time-step criteria of the corresponding CFM-FDTD scheme (c.f. Theorem 5.1.1 in [26]). As for the consistency, it can be shown that the order in space of a CFM-FDTD scheme is min{k, n}, where k is the degree of polynomial approximations of the correction function and n is the order of the FD scheme in space (c.f. Proposition Appendix A.1).

PEC wall boundary conditions
Let us now focus on PEC wall boundary conditions. By Proposition 3.1, one cannot just neglect interface conditions (1f) and (1g) for PEC problems for which J s and ρ s are unknown because this leads to an ill-posed minimization problem.
To circumvent this issue, we propose to use fictitious interface conditions, given where N k is the number of fictitious interfaces Γ k,i ⊂ Ω + ∩ Ω h Γ for k = 1, . . . , 4, n k,i is the unit normal associated with fictitious interface Γ k,i , and H * and E * are respectively finite difference approximations of H + and E + in Ω + . In subsection 3.2, we provide more details on the implementation of these fictitious interface conditions. The functional to minimize is then given bỹ where c f > 0 and c p > 0 are penalization coefficient, and N H = N 2 + N 4 and N E = N 1 + N 3 are the total number of fictitious interfaces for each electromagnetic field. The problem statement is then where W = V . The following proposition guarantees that there is a global minimizer for an appropriate choice of the penalization coefficient c f and fictitious interfaces. Proof. The demonstration is similar to the proof presented in Proposition 3.1.
The Hessian matrix coming from problem (10) is given by where a, c and e are the same as in Proposition 3.1, but One can notice that b ≥ 0 and d ≥ 0. However, since there is fictitious interfaces that are collinear and orthogonal to each plane defined by the axis of the coordinate system and for each fictitious interface condition, we haveb = 0 andd = 0. Requiring the Hessian matrix to be positive definite, one finds the following conditions: Condition (11) is always satisfied and condition (12) leads to the following cri- wherec = c f h . For a sufficiently largec, criterion (13) is satisfied.
Let us now investigate the impact of the CFM based on minimization problem (10) on the stability of the original FD scheme. As shown previously, we have a system of the form (7) with D coming from problem (10). However, this is not completely accurate since fictitious interface conditions (8) depend on FD solutions in Ω + . Computing Gateaux derivatives and using a necessary condition to obtain a minimum, we obtain where c contains coefficients of polynomial approximations of the correction function, and b f and b Γ are associated with terms using respectively fictitious interfaces and embedded boundaries. Moreover, we can define a linear operator (10) is well-posed for appropriate fictitious interfaces, c p and c f , which leads to where L is a finite difference operator, I is the identity operator and A is a linear operator that computes polynomial approximations of the correction function at nodes where it is needed. Since problem (10) is well-posed, we assume that the term L A M −1 b Γ can be bounded. It is then sufficient to investigate the stability of to identify any time-step criteria of the corresponding CFM-FDTD scheme [26].
We remark that we recover the original FDTD scheme in the limit when c f → 0 and therefore its properties. We therefore assume that we should be close to the stability condition of the original FDTD scheme for a sufficient small c f . This assumption is supported by numerical examples performed in Section 6.
In order to choose an appropriate penalization coefficient c f , one should consider the following constraints. First, the priority should be given to embedded boundary conditions in problem (10), that is c p > c f . Second, the weight associated with fictitious interface conditions in the minimization problem should diminish as the length of local patches h goes to zero, that is when the mesh grid size diminishes. This again enforces embedded boundary conditions and avoids any stability issues of CFM-FDTD schemes as the time-step size is refined. However, a too small value of c f could lead to poorly conditioned matrices coming from the minimization problem (10).

Implementation of Fictitious Interface Conditions
This short subsection focus on technical details concerning the implementation of fictitious interface conditions (8).
The functional (9)   and Section 5, we propose an initialization strategy for the Yee scheme and a fourth-order FDTD scheme.

Computation of local patches
The computation of an appropriate local patch Ω h Γ is essential for the CFM. The well-posedness of problem (2) and (10) highly depends on the representation of the embedded boundary within the local patch Ω h Γ . Hence, an appropriate local patch is of foremost importance to obtain an accurate approximation of a correction function. In previous CFM-FDTD schemes, a "Node Centered" approach is used to compute local patches [20,24]. This approach consists to define a local patch and solve a minimization problem for each node to be corrected.
Even though it is more expensive than other constructions of local patches, "Node Centered" approaches have the benefit to guarantee the uniqueness of the correction function at a given node. Hence, a common discrete measure of the divergence for staggered FDTD schemes is conserved for some nodes close to the embedded boundary [24].
In this work, as it is done for other immersed boundary methods [27,28], we directly discretize the embedded boundary. Let us now summarize this approach. For simplicity, let us consider Ω ⊂ R 2 . Assume a given embedded boundary Γ that can be parametrized with respect to the parameter s ∈ [s a , s b ].
The number of local patches is given by where L Γ is the estimated arc length of Γ , h is mesh grid size and α is positive constant. In this work, we use α = 2. Hence, centre points of local patches are given by x c,i = (x(s i ), y(s i )), where s i = s a + i ∆s for i = 0, . . . , N s − 1 and ∆s = s b −sa Ns−1 . For a given node to be corrected at x D , we find the closest x c,i and associate the corresponding local patch to x D . We therefore guarantee the uniqueness of the correction function to each node to be corrected while reducing the computational cost, particularly for large stencils, when compared with "Node Centered" approaches. The local patches are square and aligned with the computational grid. The length in space of local patches is h = β h, where h is the mesh grid size and β is a positive constant. It is worth mentioning that the parameter β is chosen in such a way that enough fictitious interfaces can be generated within Ω h Γ and that all nodes to be corrected are associated with a local patch.
As for the time interval I h Γ = [t n −∆t Γ , t n ], we choose ∆t Γ in such a way that I h Γ include the number of time steps needed to construct space-time interpolants associated with fictitious interface conditions.

Application of the CFM to the Yee Scheme
In this section, we apply the CFM to the well-known Yee scheme [16], which is a popular FDTD scheme in computational electromagnetics, with a particular attention on its initialization. Let us recall that the Yee scheme uses a staggered grid in both space and time. We then need to adapt the functional (9), and more precisely the interval of integration in time of fictitious interface conditions, in order to consider a staggered grid in time. Finally, we conclude with pros and cons of such an approach. In the following, we assume that the parameter β has been chosen in such a way that enough fictitious interfaces have been generated within Ω h Γ (see subsection 3.3) and we therefore focus on the time component. Let us first define a staggered grid in time. We consider a time interval I = [0, T ] subdivided into N t subintervals of length ∆t. We then have t n := n ∆t for n = 0, . . . , N t and t n+1/2 := (n + 1/2) ∆t for n = −1, . . . , N t − 1. The magnetic and electric field are respectively defined at t n+1/2 and t n .
According to the Yee scheme, we first compute H 1/2 using initial conditions H −1/2 and E 0 , as illustrated in Fig. 3. In this case, the CFM-Yee scheme needs to provide corrections for the electric field at t 0 , that is D 0 E . The time interval of local patches is then I h Γ = [t −1/2 , t 0 ]. At first sight, we do not have enough information in time within local patches to build accurate enough space-time interpolants for fictitious interface conditions. However, by Faraday's  (9) is also For the computation of E 1 , one needs correction functions for the magnetic field at t 1/2 , that is D 1/2 H , as illustrated in Fig. 3. The time interval of local patches is then I h Γ = [t −1/2 , t 1/2 ]. In this situation, we construct a first-degree polynomial in time for the electric field using again E 0 and ∂ t E −1/2 . As for the magnetic field, we use H −1/2 and H 1/2 to compute a first-degree polynomial in time. The interval of integration in time associated with fictitious interface conditions involving the magnetic field is then [t −1/2 , t 1/2 ] while the one associated with the electric field is [t −1/2 , t 0 ].
Once the initialization of the proposed CFM-Yee scheme is done, we only have two cases to consider, as illustrated in Fig. 4. The first case involves the Another avenue to initialize the proposed CFM-Yee scheme, although it is very specific to some applications, is to consider that H + and E + in the vicinity of the embedded boundary remain unchanged for t ≤ t 0 . Hence, the numerical strategy described previously for Fig. 4 can be directly used. As an example, this approach could be useful for scattering problems.
Finally, the main disadvantage of this approach is the computation cost associated with minimization problems of the CFM at each update of electromagnetic fields. In fact, the CFM consumes most of the computational time when compared with the finite-difference part. However, a parallel implementation of the computation of approximations of the correction function can be performed since minimization problems associated with local patches are independent at a given time step. We refer to [25] for more details about the benefits of a parallel implementation of the CFM. It is also worth mentioning that we do not have to keep whole previous solutions but only approximations associated with fictitious interfaces. Despite this drawback, the proposed CFM-Yee scheme could achieve second-order convergence for appropriate approximations of the correction function (see Proposition Appendix A.1) while treating various complex geometries of the embedded boundary without significantly increasing the complexity of the numerical approach. Moreover, it can also be implemented as a black-box for existing softwares that use the Yee scheme.

Application of the CFM on a Fourth-Order Staggered FDTD Scheme
In this section, we introduce a CFM-FDTD scheme based on a fourth-order staggered FDTD scheme. The staggered space and time grid are defined as in the Yee scheme. Spatial derivatives are estimated using the fourth-order centered approximation. As for time derivatives, many avenues can be chosen, such as staggered Adams-Bashforth or staggered backward differentiation methods [29]. In this work, we choose a fourth-order staggered free-parameter multistep method introduced in [29], which has a maximum imaginary stability boundary close to the leapfrog method used in the Yee scheme.
In the following, we first describe a fourth-order staggered free-parameter multistep method that is used to discretize time derivatives. We assume that previous solutions needed for the initialization of the multistep method are given. Afterward, we introduce the associated CFM-FDTD scheme with a particular attention on the time component. As in Section 4, we assume that the parameter β has been chosen in such a way that enough fictitious interfaces have been generated within Ω h Γ . Let us consider ∂ t H = f H (E) and ∂ t E = f E (H). The considered fourthorder free-parameter method is given by where β 1 = t, β 2 = s, β 3 = 1 22 s + 12 11 , with s = −1 and t = 1.045.
Let us now introduce the corresponding CFM-FDTD scheme. The time scheme (15) needs to compute f H (E) at t n , t n−1 and t n−2 , and f E (H) at t n+1/2 , t n−1/2 and t n−3/2 . Hence, we need to keep previous corresponding correction functions. For the computation of H n+1/2 , we set I h Γ = [t n−7/2 , t n ] and compute D n E . It is worth mentioning that, at the first update of the magnetic field to estimate H 1/2 , we compute the correction function D E at t 0 , t −1 and t −2 , and D H at t −1/2 and t −3/2 . For the update of the electric field, that is E n+1 , we set I h Γ = [t n−3 , t n+1/2 ] and compute D n+1/2 H . As for the Yee scheme, we need to adapt the interval of integration in time of fictitious interface conditions using a similar procedure as in Section 4 for both cases.

Numerical Examples
In the following, we name the CFM-FDTD scheme based on the Yee scheme as CFM-Yee scheme while the one based on a fourth-order staggered FDTD scheme is named CFM-4 th scheme. The error of U = [H x , H y , E z ] T at t n is computed using approximations and analytic solutions of the magnetic field and the electric field at respectively t n − ∆t 2 and t n because of the staggered grid in time.

Transverse Magnetic Mode
Let us consider the transverse magnetic (TM z ) mode. The unknowns are H x (x, y, t), H y (x, y, t) and E z (x, y, t). For a domain Ω ⊂ R 2 and constant physical parameters, Maxwell's equations are then simplified to with the associated boundary, embedded boundary and initial conditions.

A 2-D Staggered Space Discretization
In this subsection, we present a staggered grid in space. Let us consider for respectively the second and fourth order centered finite-difference. Finally, as it is commonly used, we impose E z = 0 and H x = H y = 0 in PEC subdomains.

Problems with an Analytic Solution
In this subsection, we perform numerical examples with analytic solutions to assess the impact of the penalization coefficient c f and to verify the proposed numerical approach. The domain Ω is divided in two subdomains Ω + and Ω − .
Periodic boundary conditions are used on all ∂Ω. We set c p = 1. The mesh grid size is such that h = ∆x = ∆y. We use second and third degree polynomial approximations of the correction function for respectively the CFM-Yee scheme and the CFM-4 th scheme. Hence, this should lead to a second and third order convergence in L 2 -norm (see Proposition Appendix A.1). We set h = 7 h, and we construct E * and H * in Ω + using at least a second degree interpolating polynomial in space for both schemes. Since Ω − is a PEC subdomain, the embedded boundary Γ then encloses subdomain Ω + . The embedded boundary is a circle centered at (0, 0) with unit radius. The physical parameters are = 1 and µ = 1. The time-step size is ∆t = h 2 for both CFM-FDTD schemes. In subdomain Ω + , the solution in cylindrical coordinates is given by where α i,j is the j-th positive real root of the i-order Bessel function of first kind J i . In this numerical example, we choose i = 6 and j = 2.
Let us assess the impact of the penalization coefficient c f on the proposed CFM-FDTD schemes.    Γ is a square of unit length centered at (0, 0) . The physical parameters are = 1 and µ = 1. The time-step size is ∆t = h 2 for both CFM-FDTD schemes.
In Ω + , the solution is given by H + x (x, y, t) = − π n ω sin(m π x) cos(n π y) sin(ω t), H + y (x, y, t) = π m ω cos(m π x) sin(n π y) sin(ω t), E + z (x, y, t) = sin(m π x) cos(n π y) cos(ω t), where ω = π √ m 2 + n 2 with m, n > 0 [30]. In this numerical example, we choose m = n = 4. Convergence plots are illustrated in Fig. 7 for both proposed CFM-FDTD schemes. A second and fourth order convergence are observed for respectively the CFM-Yee scheme and the fourth-order CFM-FDTD scheme in L 2 -norm. The obtained convergence orders are in agreement with the theory.

Two Concentric PEC Cylinders Problem
This problem considers a holed PEC containing a PEC cylinder as illustrated in Figure 8. It is recalled that Ω − is a PEC subdomain. We therefore have subdomain Ω + enclosed by two PEC walls on Γ 1 and Γ 2 . There is two circular The time-step size is ∆t = h 4 for both CFM-FDTD schemes. In subdomain Ω + , the solution in cylindrical coordinates is given by where α ≈ 1.76368380110927, ω ≈ 9.813695999428405, and J i and Y i are the iorder Bessel function of respectively first and second kind [8]. Convergence plots are illustrated in Fig. 9 for both proposed CFM-FDTD schemes. As expected, we observe a second and third order convergence for respectively the CFM-Yee scheme and the CFM-4 th scheme in L 2 -norm.

Problems with a Manufactured Solution
Let us now consider more complex embedded boundaries. To our knowledge, there is no analytical solution for Maxwell's equations with arbitrary embedded boundaries. We therefore consider that embedded boundary conditions are given byn × E + = a(x, t) on Γ × I, where a(x, t) and b(x, t) are known functions. The physical parameters are = 1 and µ = 2. The time-step size is ∆t = h 2 for both CFM-FDTD schemes. We also consider the 5-star and 3-star embedded boundary illustrated in Fig. 10.
Nevertheless, this allows us to assess the performance of the proposed CFM-FDTD schemes in a more general framework. The time interval is I = [0, 1].
We set h = 7 h for both embedded boundaries. The other parameters are the same as problems with analytical solutions in subsection 6.3. Convergence plots for each geometry of the embedded boundary are illustrated in Fig. 11 using both proposed CFM-FDTD schemes. For the 5-star boundary, a second-order convergence is observed for the CFM-Yee scheme. We also observe a fourth order convergence for the CFM-4 th scheme. As for the 3-star boundary, we observe a second-order convergence for the CFM-Yee scheme while a third-order  convergence is obtained for the CFM-4 th scheme.
where σ = 0.1 and γ = −0.3. It is worth mentioning that we use electromagnetic fields given in (17) to compute all previous solutions needed to initialize the timestepping method presented in Section 5. It is recalled that we set H − x = H − y = 0 and E − z = 0 in the PEC subdomain, that is Ω − . The reference solution U is computed using h = 1 1620 and the CFM-4 th scheme. All nodes used for H x , H y and E z in coarser grids with h ∈ { 1 20 , 1 60 , 1 180 , 1 540 } are also part of the finer grid. Fig. 12 illustrates convergence plots for each geometry of the embedded PEC using both CFM-FDTD schemes. We observe a second and fourth order convergence for respectively the CFM-Yee scheme and the CFM-4 th scheme. Fig. 13, Fig. 14        to avoid any stability issues as the time-step size diminishes. We have applied the CFM to the well-known Yee scheme and a fourth-order staggered FDTD scheme. It is worth mentioning that the proposed numerical strategy can be implemented as a black-box to existing softwares. Based on numerical examples, it has been shown that the proposed CFM-FDTD schemes can handle embedded PEC problems with various geometries of the boundary while retaining high-order convergence and without significantly increasing the complexity of the proposed numerical approach. whereÛ is the vector of true solution values, A is a rectangular matrix with either 0 or ±1 as components,D is the vector of true correction function values, L is a spatial finite difference operator of order n that approximates q-order derivatives and F is a source term. A (k + 1)-order approximation of the correction function leads to a corrected FD scheme of order min{n, k − q + 1} in space.

Acknowledgments
Proof. A (k + 1)-order approximation ofD leads to where h = β h is the length of the space-time patch, β is a positive constant and h is the mesh grid size. The discrete operator L, that approximates q-