# Stable GSTC formulation for Maxwell's equations

N. Lebbe, Kim Pham, Agnès Maurel

► To cite this version:

N. Lebbe, Kim Pham, Agnès Maurel. Stable GSTC formulation for Maxwell's equations. IEEE Transactions on Antennas and Propagation, 2022, 70 (8), pp.6825-6840. 10.1109/TAP.2022.3161436. hal-03203013

HAL Id: hal-03203013

https://inria.hal.science/hal-03203013

Submitted on 20 Apr 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Stable GSTC formulation for Maxwell's equations

Nicolas Lebbe*, Kim Pham†, and Agnès Maurel‡, * Université Côte d'Azur, Inria, CNRS, LJAD, 06902 Sophia Antipolis Cedex †IMSIA, CNRS, EDF, CEA, ENSTA Paris, Institut Polytechnique de Paris, 828 Bd des Maréchaux, 91732 Palaiseau, France, ‡Institut Langevin, ESPCI Paris, PSL University, CNRS, 1 rue Jussieu, 75005 Paris, France,

***Abstract***—We revisit the classical zero-thickness Generalized Sheet Transition Conditions (GSTCs) which are a key tool for efficiently designing metafilms able to control the flow of light in a desired way. It is shown that it is more convenient to use an enlarged formulation of the GSTC in which the original metafilm is replaced by GSTCs that exclude the layer from the physical or computational domain. These new "layer" transition conditions have the same form as their "sheet" analogues hence they do not necessitate additional complications in their use; their advantage is that they provide a well-posed problem hence guaranty the stability of numerical schemes in the time-domain. These assessments are demonstrated for an all-dielectric structure; the effective susceptibility tensors are derived thanks to asymptotic analysis combined with homogenization technique and bounds for the susceptibilities entering the balance of energy are provided. While negative constant susceptibilities appear in the classical zero-thickness GSTCs, their values in the enlarged formulation are always positive which ensure the stability of the effective problem. Validation of the effective model is provided by means of comparison with direct numerics in two and three dimensions.

***Index Terms***—Metasurface, metafilms, generalized sheet transition conditions (GSTC), two-scale homogenization technique, asymptotic analysis, numerical implementation of the enlarged GSTCs.

## I. INTRODUCTION

Metasurfaces are smartly engineered two-dimensional structures that offer an attractive alternative to their bulk three-dimensional analogue. Because of their subwavelength thickness, they are less lossy, more compact and of relatively ease of fabrication. Besides these thin composite layers exhibit great abilities for controlling the flow of light (see review papers [1], [2], [3], [4], [5], [6]) and new perspectives are thought using metasurfaces with spatial gradients [7], [8], temporal variations [9], [10] and non-linearities [11]. To describe metasurfaces, effective models have to encapsulate the microscopic distributions of currents and fields resulting in transition conditions relating the field differences at the two opposite sides of the film to their mean values. Through this homogenization process, the effects of the microscopic distributions are encoded in macroscopic parameters which can be interpreted in terms of effective surface currents and surface polarisations. In their most general form, these conditions are called General Sheet Transmission Conditions (GSTCs) [12].

The interest for effective transition conditions predates the emergence of metafilms and the derivation of reliable models for thin homogenous layers has given rise to a vast literature, see e.g. [13]. In a general perspective, that is leaving

aside the actual composition of the film, Idemen and co-workers postulate the validity of the Maxwell's equations in the sense of distributions; doing so, they derive universal transition conditions by introducing Dirac delta-distribution of fields concentrated on a zero-thickness surface and provide the first theoretical justification of the GSTCs [14], [15], [16]. Since these pioneering works, different methods have been developed to account for actual geometries and arrangements of scattering particles in metafilms, including two-scale asymptotic homogenization [17], [18], [19], homogenization of local fields within dipole approximation [20], [21], and a plethora of retrieval methods see e.g. [22], [23], [24], [25] (a historical overview is provided in [20]). In most cases, for both homogeneous films and metafilms, the transition conditions are expressed across a zero-thickness interface, or sheet. The underlying idea is that the actual thickness can be neglected, besides it may be not well defined. Hence the determination of an effective thickness would be submitted to some arbitrary choice. In addition, models with zero-thickness transition conditions clearly depart from models assimilating a metafilm to a layer with effective bulk properties (permittivity, permeability) and effective thickness. Discussions on the drawbacks of such attempts are provided in [26], [27] and in [28] with the perspective of classical bulk- versus interface-homogenizations.

The notion of an arbitrary choice for the thickness of the film boundaries involved in the field jumps and mean values is perfectly legitimate. Indeed, approximate models result from an asymptotic process with some parameter $\delta \ll 1$ measuring the subwavelength dimensions. Being conducted up to order $O(\delta)$ (in general), an infinite family of models exists and these models are all equivalent up to terms in $O(\delta^2)$. We shall see that shifting the film thickness from zero to any $O(\delta)$ provides such family of equivalent models. This could be an incidental remark but it turns out that these approximate problems can be discriminated in terms of their well-posedness or stability. Up to now, this problem has been addressed in the context of metafilms in elastodynamics [29], [30] and in acoustics [31], [32] and it has been shown that unstable formulations may foster numerical instabilities in the time domain. Such instabilities have been reported in [33] for homogeneous thin films, and in [34], [35] for structured thin films. In all cases, instabilities appear when a negative constant coefficient entering the interfacial energy is involved. In [34], the authors conclude that negative and constant surface susceptibilities represent a non-physical system, and thus are not allowed (where constant means without dispersion). This
conclusion is in close relation with the present study since we shall see that devices as simple as dielectric homogeneous and structured films may lead to negative constant susceptibilities when expressed in zero-thickness GSTCs; in contrast we shall prove that constant negative susceptibilities are not possible when using enlarged GSTC formulation in which at least the film thickness has been restored.

The paper is organized as follow. The main results of the analysis are presented in §II for an all-dielectric metafilm composed of scattering particles embedded at the interface between two different substrates (figure 1). The form of the effective susceptibility tensors entering in the GSTCs and the resulting balance of energy are given when considering enlarged transition conditions across an excluded region that contains the particles. Explicit positive bounds for the effective susceptibilities entering the energy balance are provided. The derivation of the model is detailed in §III. To derive the effective susceptibilities we use a matched asymptotic analysis combined with two-scale homogenization. As a warm-up and to set the main ingredients of the asymptotic procedure, we first envision a homogeneous thin film sandwiched between the two substrates. Next we incorporate the two-scale homogenization tools to deal with a film composed by scattering particles. This procedure allows us to recover the expected symmetry of the susceptibility tensors and to establish the bounds thanks to minimization principles. We provide in §IV validations of the effective problem by means of comparisons with direct numerical calculations of the Maxwell equations in two- and three-dimensions (the details of the numerical implementations of the direct and effective problems are freely available https://www.comsol.fr/community/exchange/841/). Some technical calculations are collected in the appendices.

## II. SUMMARY OF THE MAIN RESULTS

### A. The actual problem

We consider the Maxwell's equations in three-dimensions $x = (x, y, z)$ for dielectric materials with permittivities denoted $\epsilon_0\epsilon_r^{\pm}$ in the two substrates, $\epsilon_0\epsilon_r^{sc}$ in the scattering particles and a uniform permeability $\mu = \mu_0$ ($\epsilon_0$ and $\mu_0$ are the free space parameters). The particles are evenly distributed at the interface $x = 0$ between the two substrates, with spacing $d_y$ along $y$ and $d_z$ along $z$. We define $S = d_yd_z$ the area of the periodic cell and $d = \sqrt{S}$. In the time domain, the electric $\mathbf{E}$, displacement $\mathbf{D}$ and magnetizing $\mathbf{H}$ fields satisfy the Maxwell's equations in the actual problem (figure 1), namely

$$ 
\begin{aligned}
\mathbf{rotH} &= \frac{\partial \mathbf{D}}{\partial t}, & \mathbf{D} &= \epsilon_0 \epsilon_r (\mathbf{x}) \mathbf{E}, \\
\mathbf{rotE} &= -\mu_0 \frac{\partial \mathbf{H}}{\partial t}, & \operatorname{div}\mathbf{H} &= 0,
\end{aligned}
\quad (1) 
$$

<figure><img src="image_2.png" /><figcaption>Fig. 1: Actual problem of a metafilm composed by dielectric scattering particles evenly distributed at the interface between two dielectric substrates.</figcaption></figure>

We recall that the above equations imply an energy-conservation law expressed in a bounded domain $\Omega$ as

$$ 
\begin{aligned}
& \frac{d\mathcal{E}}{dt} + \int_{\Sigma} \mathbf{S} \cdot \mathbf{n} \, dS = 0, \\
& \text{with } \mathcal{E} = \frac{1}{2} \int_{\Omega} (\epsilon_0 \epsilon_r^c |\mathbf{E}|^2 + \mu_0 |\mathbf{H}|^2) \, d\mathbf{x}, \quad \mathbf{S} = \mathbf{E} \times \mathbf{H},
\end{aligned}
\quad (2) 
$$

where $\mathcal{E}$ is the electromagnetic energy and $\mathbf{S}$ is the Poynting vector ($\mathbf{n}$ is the unit vector normal to $\Sigma = \partial\Omega$). In the absence of energy flux through $\Sigma$ that is to say if $\Sigma$ is associated with boundary conditions $\mathbf{E} \times \mathbf{n} = 0$ or $\mathbf{H} \times \mathbf{n} = 0$, the electromagnetic energy is conserved in $\Omega$.

### B. The effective problem

The asymptotic analysis combined with homogenization tools provides an effective model in which the array is replaced by transition conditions of the GSTC type (this will be detailed in the forthcoming §III). The effective problem applies outside the thin layer of thickness $e$ bounded by $\Gamma^-$ and $\Gamma^+$, which has been excluded from the physical space (figure 2). As previously said, the excluded region contains entirely the array of particles which means that $e$ is equal or larger than the extend of the particles along $x$. Outside this region, for $x < -e^-$ and $x > e^+$, the Maxwell equations apply

$$ \mathbf{rotH} = \epsilon_0 \epsilon_r^{\pm} \frac{\partial \mathbf{E}}{\partial t}, \quad \mathbf{rotE} = -\mu_0 \frac{\partial \mathbf{H}}{\partial t}, \quad \operatorname{div}\mathbf{H} = 0. \quad (3) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2: Effective problem in which the GSTCs (5) apply across $\Gamma^-$ ($x = -e^-$) and $\Gamma^+$ ($x = e^+$) with $e = e^- + e^+$.</figcaption></figure>
Next, across the excluded region, transition conditions have to be specified, hence we define for any field **A** = **D**, **E**, **H**,

$$ \llbracket \mathbf{A} \rrbracket = \mathbf{A}|_{e^+} - \mathbf{A}|_{-e^+}, \quad \mathbf{A}_{\text{av}} = \frac{1}{2} (\mathbf{A}|_{e^+} + \mathbf{A}|_{-e^+}), \quad (4) $$

and with the above definitions we obtain conditions of the GSTC type, namely

$$ \llbracket \mathbf{H} \times \hat{\mathbf{x}} \rrbracket = -\epsilon_0 \frac{\partial \mathbf{P}_{\parallel}}{\partial t} - \nabla_{\parallel} M_x \times \hat{\mathbf{x}}, \quad (5) $$

$$ \llbracket \mathbf{E} \times \hat{\mathbf{x}} \rrbracket = \mu_0 \frac{\partial \mathbf{M}_{\parallel}}{\partial t} - \nabla_{\parallel} P_x \times \hat{\mathbf{x}}, $$

where $\hat{\mathbf{x}}$ is the unit vector along $Ox$, $\mathbf{P} = P_x \hat{\mathbf{x}} + \mathbf{P}_{\parallel}$, $\mathbf{M} = M_x \hat{\mathbf{x}} + \mathbf{M}_{\parallel}$ (with $\mathbf{P}_{\parallel} \cdot \hat{\mathbf{x}} = \mathbf{M}_{\parallel} \cdot \hat{\mathbf{x}} = 0$), and the operator $\nabla_{\parallel} M_x \times \hat{\mathbf{x}} = \partial_z M_x \hat{\mathbf{y}} - \partial_y M_x \hat{\mathbf{z}}$ $^1$. The vectors **P** and **M** refer to the electric and magnetic polarization densities, respectively. Expectedly, we obtain zero cross-susceptibility tensors hence we have

$$ \mathbf{P} = \bar{\bar{\chi}}_{ee} \tilde{\mathbf{E}}_{av}, \quad \mathbf{M} = \bar{\bar{\chi}}_{mm} \mathbf{H}_{av}, \quad (6) $$

where the subscript “av” refers to the average values of the fields at both sides of the interface, and where $\bar{\bar{\chi}}_{ee}$, and $\bar{\bar{\chi}}_{mm}$ are the electric-to-electric and magnetic-to-magnetic surface susceptibility tensors whose forms are

$$ \bar{\bar{\chi}}_{ee} = \begin{pmatrix} \chi_{ee}^{xx} & \chi_{ee}^{xy} & \chi_{ee}^{xz} \\ \chi_{ee}^{xy} & \chi_{ee}^{yy} & \chi_{ee}^{yz} \\ \chi_{ee}^{xz} & \chi_{ee}^{yz} & \chi_{ee}^{zz} \end{pmatrix}, \quad \bar{\bar{\chi}}_{mm} = e \begin{pmatrix} -1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad (7) $$

and where the 6 effective susceptibilities are defined in terms of the solutions of so-called elementary problems; they are specified in the forthcoming Eqs. (47)-(48). Next, if the microstructure is invariant when $z \to -z$ we have $\chi_{ee}^{xz} = \chi_{ee}^{yz} = 0$, if it is invariant when $y \to -y$ then $\chi_{ee}^{xy} = \chi_{ee}^{yz} = 0$ (this is shown in appendix C). We recover that for both symmetries are realized, the susceptibility tensor is diagonal. Eventually, the average electric field has to be understood as

$$ \tilde{\mathbf{E}}_{av} = \frac{1}{\epsilon_0} D_{x,av} \hat{\mathbf{x}} + \mathbf{E}_{\parallel,av}, \quad (8) $$

as in [19], since it involves the fields being continuous at the dominant order corresponding to vanishing effect of the thin metafilm.

As previously commented, the excluded domain must satisfy two conditions. On the one hand it has to be thin which means of the same order of magnitude than $e$. On the other hand it has to be large enough in order to contain entirely the particles. Doing so, $x \in (-\infty, -e^-)$ and $x \in (e^+, +\infty)$ are composed of the substrates only both in the actual and in the effective problems.

Two explicit parameters will be used in the following which are the volume averages of the permittivity and of the inverse of the permittivity within the excluded region. Denoting the

$$ \llbracket D_x \rrbracket = -\epsilon_0 \text{div} \mathbf{P}_{\parallel}, \quad \llbracket H_x \rrbracket = -\text{div} \mathbf{M}_{\parallel}. $$

<figure><img src="image_4.png" /><figcaption>Fig. 3: Volume $\mathbf{x} \in (-e^-, e^+) \times (0, d_y) \times (0, d_z)$ comprising a single scattering particle, in which $\varphi^-$ and $\varphi^+$ are the filling fraction of the particle for $x < 0$ and $x > 0$.</figcaption></figure>

volume fractions $\varphi^-$ and $\varphi^+$ of the particles for $x < 0$ and $x > 0$ respectively (figure 3), these averages read

$$ \langle \epsilon_r \rangle^{\pm} = (1 - \varphi^{\pm}) \epsilon_r^{\pm} + \varphi^{\pm} \epsilon_r^{\text{sc}}, \quad \langle \epsilon_r^{-1} \rangle^{\pm} = \frac{(1 - \varphi^{\pm})}{\epsilon_r^{\pm}} + \frac{\varphi^{\pm}}{\epsilon_r^{\text{sc}}}, \quad (9) $$

and the averages over the whole excluded layer are given by

$$ \langle \epsilon_r \rangle = \frac{e^{-\langle \epsilon_r \rangle^{-}} + e^{+\langle \epsilon_r \rangle^{+}}}{e}, \quad \langle \epsilon_r^{-1} \rangle = \frac{e^{-\langle \epsilon_r^{-1} \rangle^{-}} + e^{+\langle \epsilon_r^{-1} \rangle^{+}}}{e}. \quad (10) $$

We now move on to the balance of energy in the effective problem, that must be the counterpart of (2). Denoting $\Omega_e$ the domain $\Omega$ punctured by the excluded domain $x \in (-e^-, e^+)$, the boundary of $\Omega_e$ is $\partial\Omega_e = \Sigma_e \cup \Gamma^+ \cup \Gamma^-$. Next, the fluxes of the Poynting vector through $\Gamma^+$ and $\Gamma^-$ make the GSTCs to appear. Specifically, we have

$$ \frac{d}{dt}(\mathcal{E} + \mathcal{E}_r) + \int_{\Sigma_e} \mathbf{S} \cdot \mathbf{n} \, dS = 0, \\ \mathcal{E} = \frac{1}{2} \int_{\Omega_e} (\epsilon_0 \epsilon_r |\mathbf{E}|^2 + \mu_0 |\mathbf{H}|^2) \, d\mathbf{x}, \\ \mathcal{E}_r = \frac{1}{2} \int_{\Gamma} \left( \mu_0 e \mathbf{H}_{av}^2 - \chi_{ee}^{xx} \frac{D_{x,av}^2}{\epsilon_0} + \epsilon_0 (\chi_{ee}^{yy} E_{y,av}^2 + \chi_{ee}^{zz} E_{z,av}^2 + 2\chi_{ee}^{yz} E_{y,av} E_{z,av}) \right) d\mathbf{x}_{\parallel}, \quad (11) $$

where $\mathbf{x}_{\parallel} = (y, z)$. In the absence of fluxes through $\Sigma_e$, the total energy $(\mathcal{E} + \mathcal{E}_r)$ is conserved. This is why, when the effective problem is implemented numerically, the interfacial energy $\mathcal{E}_r$ has to be positive. Indeed, a negative $\mathcal{E}_r$ fosters numerical instability as $\mathcal{E}_r$ diverging to $-\infty$ can be compensated by $\mathcal{E}$ diverging to $+\infty$ (an example of such behaviour is given in [35]). It is worth noticing that the above balance of energy holds for constant effective susceptibilities. When resonances are involved resulting in frequency dependent susceptibilities in the harmonic regime, $\mathcal{E}_r$ (being the opposite of the energy delivered by the metasurface) can be negative as reported in [10] (figure 6 in this reference). It is worth noticing that for resonant metasurfaces the derivation of the energy-conservation law has to be conducted differently (this has been done in [36], [37] for acoustic resonances).

<footnote n="1">¹It is worth noticing that (5) along with (3) and $D_x = \epsilon_0 \epsilon_r E_x$ provide</footnote>
Coming back to (11), we shall prove that with

$$
\begin{align}
-\chi_{\parallel}^{M} &= \frac{e^{-}}{\langle \epsilon_r \rangle^{-}} + \frac{e^{+}}{\langle \epsilon_r \rangle^{+}}, & -\chi_{\parallel}^{m} &= e \langle \epsilon_r^{-1} \rangle, \tag{12} \\
\chi_{\perp}^{m} &= \frac{e^{-}}{\langle \epsilon_r^{-1} \rangle^{-}} + \frac{e^{+}}{\langle \epsilon_r^{-1} \rangle^{+}}, & \chi_{\perp}^{M} &= e \langle \epsilon_r \rangle \nonumber
\end{align}
$$

(and $\langle \epsilon_r \rangle$ and $\langle \epsilon_r^{-1} \rangle$ defined in (10)), the diagonal term of the susceptibility tensor has the following bounds

$$
\left\{
\begin{array}{l}
0 < -\chi_{\parallel}^{M} \le -\chi_{cc}^{xx} \le -\chi_{\parallel}^{m}, \\
0 < \chi_{\perp}^{m} \le \chi_{cc}^{\alpha\alpha} \le \chi_{\perp}^{M},
\end{array}
\right.
\quad \alpha = y, z \tag{13}
$$

and for non-zero off-diagonal term $\chi_{cc}^{yz}$, we also have

$$
\chi_{cc}^{yz} < \min \left( \sqrt{(\chi_{cc}^{yy} - \chi_{\perp}^{m})(\chi_{cc}^{zz} - \chi_{\perp}^{m})}, \sqrt{(\chi_{cc}^{yy} - \chi_{\perp}^{M})(\chi_{cc}^{zz} - \chi_{\perp}^{M})} \right). \quad (14)
$$

It is worth noticing that (13) provides absolute bounds for the diagonal terms of the susceptibility tensor. In contrast, (14) does not provide absolute bounds for $\chi_{cc}^{yz}$. Rather, the inequalities tell us that once $(\chi_{cc}^{yy}, \chi_{cc}^{zz})$ are known, then $\chi_{cc}^{yz}$ is bounded²

### III. DERIVATION OF THE EFFECTIVE PROBLEM

#### A. A warm-up: asymptotic analysis for a homogeneous film

The case of a homogeneous film with permittivity $\epsilon_r^x$ and thickness $e^x$ (figure 4) is interesting because it allows us to introduce the tools of the asymptotic analysis that will be used for a structured metafilm. Besides, it provides explicit expressions of the susceptibilities. The asymptotic analysis is conducted owing to the small parameter $\delta = e$ much smaller than $\lambda$ the minimum wavelength imposed by the source. Note that often, non-dimensional form of the equations are used with $x \to x/\lambda$ and setting $\delta = e/\lambda \ll 1$; without loss of generality, we set $\lambda = 1$ to avoid multiple coordinate notations. Far from the film, the fields have macroscopic variations (at the wavelength-scale), and we define expansions for the macroscopic fields $\mathbf{A} = \mathbf{H}, \mathbf{D}, \mathbf{E}$ as

$$
\mathbf{A} = \sum_{i=0}^{\infty} \delta^i \mathbf{A}^i(\mathbf{x}, t), \quad (15)
$$

In contrast, within the film, rapid variations occur which are accounted for by introducing the microscopic coordinate $\xi_x = \frac{x}{e}$. Accordingly, we define the expansions for the microscopic fields

$$
\mathbf{A} = \sum_{i=0}^{\infty} \delta^i \mathbf{a}^i(\xi_x, \mathbf{x}_{\parallel}, t), \quad (16)
$$

where $\mathbf{x}_{\parallel} = (y, z)$ is kept to account for slow variations along the film. The equations satisfied by the macroscopic fields $(\mathbf{H}^i, \mathbf{E}^i, \mathbf{D}^i)$ are determined by inserting (15) in (1) and by

²Equivalently, we shall establish that for any $(E_y, E_z)$, denoting the quadratic form Q.F. = $(\chi_{cc}^{yy} E_{y,av}^2 + \chi_{cc}^{zz} E_{z,av}^2 + 2\chi_{cc}^{yz} E_{y,av} E_{z,av})$, we have

$$
e \langle \epsilon_r \rangle |\mathbf{E}_{\parallel av}|^2 \ge \text{Q.F.} \ge \left( \frac{e^{-}}{\langle \epsilon_r^{-1} \rangle^{-}} + \frac{e^{+}}{\langle \epsilon_r^{-1} \rangle^{+}} \right) |\mathbf{E}_{\parallel av}|^2 > 0.
$$

with $|\mathbf{E}_{\parallel av}|^2 = (E_{y,av}^2 + E_{z,av}^2)$.

<figure><img src="image_5.png" /><figcaption>Fig. 4: A thin film with thickness $e^x$ and permittivity $\epsilon_r^x$ at the interface between two different substrate is the simplest case of layer.</figcaption></figure>

identifying the terms with same powers in $\delta$. As $\epsilon_r(\mathbf{x}) = \epsilon_r^\pm$ outside the film, $(\mathbf{H}^i, \mathbf{E}^i)$ simply satisfy

$$
\begin{align}
\operatorname{rot}_x \mathbf{H}^i &= \epsilon_0 \epsilon_r^\pm \partial_t \mathbf{E}^i, && \text{in } \mathcal{X}^\pm, \tag{17} \\
\operatorname{rot}_x \mathbf{E}^i &= -\mu_0 \partial_t \mathbf{H}^i, && \operatorname{div}_x \mathbf{H}^i = 0, \nonumber
\end{align}
$$

and there are missing boundary conditions when $x \to 0^\pm$. Alike, the equations satisfied by the microscopic fields $(\mathbf{H}^i, \mathbf{E}^i, \mathbf{d}^i)$ are determined by inserting (16) in (1) with the differential operator

$$
\nabla \to \frac{1}{\delta} \frac{\partial}{\partial \xi_x} \hat{\mathbf{x}} + \nabla_{\parallel}, \quad (18)
$$

and identifying the terms with same powers in $\delta$. The resulting problems have to be complemented with boundary conditions when $\xi_x \to \pm\infty$. The missing conditions on the macroscopic and microscopic problems are provided simultaneously by so-called matching conditions which tell us that the two expansions have to match in some intermediate regions when $x \to 0^\pm$ and $\xi_x \to \pm\infty$, specifically up to terms $O(\delta^2)$

$$
\mathbf{A}^0(\mathbf{x}, t) + \delta \mathbf{A}^1(\mathbf{x}, t) \underset{\xi_x \to \pm\infty}{\sim} \mathbf{a}^0(\xi_x, \mathbf{x}_{\parallel}, t) + \delta \mathbf{a}^1(\xi_x, \mathbf{x}_{\parallel}, t).
$$

To deduce the resulting conditions at each order in $\delta$, we use Taylor expansions $\mathbf{A}^0(\mathbf{x}, t) = \mathbf{A}^0(\mathbf{x}, t) + \delta\xi_x\partial_x\mathbf{A}^0(0^\pm, t) + \dots$. It follows that the matching conditions read

$$
\lim_{\xi_x \to \pm\infty} \mathbf{a}^0 = \mathbf{A}^0|_{0^{\pm}}, \quad (19)
$$

at the zero-order (with $\mathbf{A}|_{0^\pm} = \mathbf{A}(\mathbf{x} = 0^\pm, \mathbf{x}_{\parallel}, t)$) and at the first-order

$$
\lim_{\xi_x \to \pm\infty} \left( \mathbf{a}^1 - \xi_x \left. \frac{\partial \mathbf{A}^0}{\partial x} \right|_{0^\pm} \right) = \mathbf{A}^1|_{0^\pm} . \quad (20)
$$

1) Zero-order problem: Let us start the asymptotic analysis for the homogeneous film. Because of the form of the differential operator in (18), inserting (16) in (1) provides equations at the dominant order $O(\delta^{-1})$. Specifically, we have $\partial_{\xi_x} \mathbf{h}^0 = \partial_{\xi_x} \mathbf{e}^0 = \mathbf{0}$ and $\partial_{\xi_x} \mathbf{d}^0 = \mathbf{0}$. In virtue of the matching conditions (19), we thus have

$$
\mathbf{h}^0 = \mathbf{H}^0|_0, \quad \mathbf{d}^0 = \mathbf{D}^0|_0, \quad \mathbf{e}^0 = \mathbf{E}^0|_0. \quad (21)
$$

Because of its thinness in $O(\delta)$, the film is not seen at the dominant order $O(1)$. Hence we have to go to the next
order and we shall consider the fields ($\mathbf{A}^0 + \delta\mathbf{A}^1$) which approximate the actual fields up to $O(\delta^2)$. If we do so, it makes sense to restore the actual thickness of the film which means to express the transition conditions between $x = -e^-$ and $x = e^+$ rather than between $x = 0^-$ and $x = 0^+$. As the macroscopic fields are continuous in $\mathcal{X}^{\pm}$, they can be Taylor expanded between $x = 0^-$ and $-e^-$ and between $x = 0^+$ and $e^+$. It follows that up to $O(\delta^2)$ we have

$$ 
\begin{aligned}
\mathbf{A}|_{-e^+} &= \mathbf{A}|_{0^-} - e^- \partial_x \mathbf{A}|_{0^-}, \\
\mathbf{A}|_{e^+} &= \mathbf{A}|_{0^+} + e^+ \partial_x \mathbf{A}|_{0^+},
\end{aligned}
\qquad (22) $$

resulting in

$$ 
\begin{aligned}
\llbracket \mathbf{H}^0 \times \hat{\mathbf{x}} \rrbracket &= (e^+ \partial_x \mathbf{H}^0|_{0^+} + e^- \partial_x \mathbf{H}^0|_{0^-}) \times \hat{\mathbf{x}}, \\
\llbracket \mathbf{E}^0 \times \hat{\mathbf{x}} \rrbracket &= (e^+ \partial_x \mathbf{E}^0|_{0^+} + e^- \partial_x \mathbf{E}^0|_{0^-}) \times \hat{\mathbf{x}}.
\end{aligned}
\qquad (23) $$

Expressing the above conditions in terms of the fields being continuous at $x = 0$ owing to (17) we obtain

$$ 
\left\{
\begin{aligned}
\llbracket \mathbf{H}^0 \times \hat{\mathbf{x}} \rrbracket &= -\varepsilon_0 (e^+ \varepsilon_r^+ + e^- \varepsilon_r^-) \partial_t \mathbf{E}_{\parallel}^0 + e \nabla_{\parallel} H_x^0 \times \hat{\mathbf{x}}, \\
\llbracket \mathbf{E}^0 \times \hat{\mathbf{x}} \rrbracket &= \mu_0 e \partial_t \mathbf{H}_{\parallel}^0|_0 + \left( \frac{e^-}{\varepsilon_r^-} + \frac{e^+}{\varepsilon_r^+} \right) \frac{\nabla_{\parallel} D_x^0|_0}{\varepsilon_0} \times \hat{\mathbf{x}}.
\end{aligned}
\right.
\qquad (24) $$

The right hand-side terms are as small as $e = \delta$ and have to be complemented by the contributions of $\delta \llbracket \mathbf{A}^1 \rrbracket$.

2) First-order problem: The useful problems at the order $O(1)$ are set on $(\mathbf{h}^1 \times \hat{\mathbf{x}})$ and $(\mathbf{e}^1 \times \hat{\mathbf{x}})$. They read

$$ 
\begin{aligned}
\frac{\partial}{\partial \xi_x} (\mathbf{h}^1 \times \hat{\mathbf{x}}) &= -\varepsilon_0 \varepsilon_r \partial_t \mathbf{E}_{\parallel}^0|_0 + \nabla_{\parallel} H_x^0|_0 \times \hat{\mathbf{x}}, \\
\frac{\partial}{\partial \xi_x} (\mathbf{e}^1 \times \hat{\mathbf{x}}) &= \mu_0 \partial_t \mathbf{H}_{\parallel}^0|_0 + \frac{\nabla_{\parallel} D_x^0|_0}{\varepsilon_0 \varepsilon_r} \times \hat{\mathbf{x}}.
\end{aligned}
\qquad (25) $$

As the right hand side terms have a dependence in $\xi_x$ through $\varepsilon_r(\xi_x)$ only, the solutions read

$$ 
\begin{aligned}
\mathbf{h}^1 \times \hat{\mathbf{x}} &= (\mathbf{h}^1 \times \hat{\mathbf{x}})|_{\xi_x=0} \\
&\quad -\varepsilon_0 \int_0^{\xi_x} \varepsilon_r(\xi) d\xi \partial_t \mathbf{E}_{\parallel}^0|_0 + \varepsilon_x \nabla_{\parallel} H_x^0|_0 \times \hat{\mathbf{x}}, \\
\mathbf{e}^1 \times \hat{\mathbf{x}} &= (\mathbf{e}^1 \times \hat{\mathbf{x}})|_{\xi_x=0} \\
&\quad + \varepsilon_x \mu_0 \partial_t \mathbf{H}_{\parallel}^0|_0 + \int_0^{\xi_x} \frac{d\xi}{\varepsilon_r(\xi)} \frac{\nabla_{\parallel} D_x^0|_0}{\varepsilon_0} \times \hat{\mathbf{x}}.
\end{aligned}
$$

Accounting for the matching conditions (20) and passing to the limit $\xi_x \to \pm\infty$, the (diverging) terms linear in $\xi_x$ cancel and we obtain

$$ 
\left\{
\begin{aligned}
\delta \llbracket \mathbf{H}^1 \times \hat{\mathbf{x}} \rrbracket &= -\varepsilon_0 (e^{sc} \varepsilon_r^{sc} - e^{sc^-} \varepsilon_r^{sc^-} - e^{sc^+} \varepsilon_r^{sc^+}) \partial_t \mathbf{E}_{\parallel}^0|_0, \\
\delta \llbracket \mathbf{E}^1 \times \hat{\mathbf{x}} \rrbracket &= \left( \frac{e^{sc}}{\varepsilon_r^{sc}} - \frac{e^{sc^-}}{\varepsilon_r^{sc^-}} - \frac{e^{sc^+}}{\varepsilon_r^{sc^+}} \right) \frac{\nabla_{\parallel} D_x^0|_0}{\varepsilon_0} \times \hat{\mathbf{x}}.
\end{aligned}
\right.
\qquad (26) $$

The final transition conditions link $\llbracket \mathbf{A} \rrbracket = \llbracket \mathbf{A}^0 + \delta \mathbf{A}^1 \rrbracket$ to the average fields as defined in (4). Gathering the contributions (24) and (26) we eventually obtain

$$ 
\left\{
\begin{aligned}
\llbracket \mathbf{H} \times \hat{\mathbf{x}} \rrbracket &= -\varepsilon_0 e \langle \varepsilon_r \rangle \partial_t \mathbf{E}_{\parallel,av} + e \nabla_{\parallel} H_{x,av} \times \hat{\mathbf{x}}, \\
\llbracket \mathbf{E} \times \hat{\mathbf{x}} \rrbracket &= \mu_0 e \partial_t \mathbf{H}_{\parallel,av} + e \langle \varepsilon_r^{-1} \rangle \frac{\nabla_{\parallel} D_{x,av}}{\varepsilon_0} \times \hat{\mathbf{x}},
\end{aligned}
\right.
\qquad (27) $$

since in the present case

$$ e \langle \varepsilon_r \rangle = e^{sc} \varepsilon_r^{sc} + (e^- - e^{sc^-}) \varepsilon_r^{-} + (e^+ - e^{sc^+}) \varepsilon_r^{+}. $$

By identifying the above transition conditions with the GSTC in (5), we obtain

$$ \bar{\chi}_{\text{mm}} = e \begin{pmatrix} -1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad \bar{\chi}_{\text{ee}} = e \begin{pmatrix} -\langle \varepsilon_r^{-1} \rangle & 0 & 0 \\ 0 & \langle \varepsilon_r^- \rangle & 0 \\ 0 & 0 & \langle \varepsilon_r^+ \rangle \end{pmatrix}. \qquad (28) $$

It results that the interfacial energy in (11) reads

$$ \mathcal{E}_{\parallel} = \frac{1}{2} \int_{\Gamma} e \left( \mu_0 \mathbf{H}_{av}^2 + \varepsilon_0 \langle \varepsilon_r \rangle \mathbf{E}_{av}^2 \right) d\mathbf{x}_{\parallel}. $$

The conditions (27) are those obtained in [38] in which the authors directly implement the transition conditions in a FDTD scheme from the Ampere's and Faraday's law (Eq. (4)) is our $\llbracket H_z \rrbracket$ and Eq. (7) is our $\llbracket E_z \rrbracket$ up to a second gradient term, see also [13] for a comparison of several schemes among which the unstable zero-thickness condition of [33]. The identification with (27) is made easy identifying $e$ to the grid step $\Delta x$.

## B. Asymptotic analysis for a metafilms

We now move on to the case of a structured metafilm. Most of the work has been done in the previous section and we shall follow basically the same procedure. The additional difficulty comes from the necessity to keep microscopic coordinates along the three dimensions, as the evanescent field excited in the vicinity of each particle have rapid (or microscopic) variations in the three directions.

According to what has been said, we introduce the microscopic coordinate $\xi = \frac{\xi}{\delta}$ and we now choose $\delta = d$ where $d = \sqrt{d_y d_z}$; this provides a rescaled unit cell $Y$ with $|Y| = 1$ (figure 5). Next, far from the metafilm, the macroscopic fields are expended as in (15). In the vicinity of the film, expansions being the equivalent of (16) are used, with

$$ \mathbf{A} = \sum_i \delta^i \mathbf{a}^i(\xi, \mathbf{x}_{\parallel}, t). \qquad (29) $$

As in the classical two-scale homogenization of bulk-metamaterials, the tangential in-plane microscopic variations fields are handled by $\xi_{\parallel}$ while the macroscopic ones are handled by $x_{\parallel}$. Hence, we add the classical hypothesis that the terms $\mathbf{a}^i$ are $\xi_{\parallel}$-periodic. It follows that analogues problems of the cell problems in the classical homogenization will appear. We call them elementary problems and they are now set on $\xi \in \mathcal{Y}$ a strip infinite along $\xi_x$ with $\xi_{\parallel} \in Y$ (figure 5).

The macroscopic fields $(\mathbf{H}^i, \mathbf{E}^i, \mathbf{D}^i)$ satisfy (17) as in the previous section. The equations satisfied by the microscopic fields $(\mathbf{h}^i, \mathbf{e}^i, \mathbf{d}^i)$ are determined by inserting (29) in (1) with now the differential operator being

$$ \nabla \to \frac{1}{\delta} \nabla \xi + \nabla_{\parallel}. \qquad (30) $$

Eventually, the matching conditions between the macroscopic and the microscopic fields read as in (19)-(20). Note
<figure><img src="image_6.png" /><figcaption>Fig. 5: Elementary cell $\mathcal{Y} = \{(\xi_x, \xi_\|) \in \mathbb{R} \times \mathcal{Y}\} (\xi_\| = (\xi_y, \xi_z), |\mathcal{Y}| = 1)$ in which the static problems (35)-(36) are set.</figcaption></figure>

that these conditions tell us that the microscopic fields $(\mathbf{h}^0, \mathbf{e}^0, \mathbf{d}^0)$ do not depend on $\xi_\|$ when $\xi_x \to \pm\infty$, a result that we shall recover. This result is quite intuitive since the rapid variations are associated to evanescent fields excited at the scatterers hence vanishing when moving away from the array (these evanescent fields were absent for homogeneous films which are translational invariant). The complexity due to the presence of evanescent fields makes that explicit microscopic solutions are not available.

1) **Zero-order problem:** As for the film, inserting (29) in (1) along with (30) provides equations at the dominant order $O(\delta^{-1})$. In particular, we obtain for $\mathbf{h}^0$ the set of equations

$$ egin{cases} 
\mathbf{rot}_{\xi} \mathbf{h}^0 = \mathbf{0}, & \text{div}_{\xi} \mathbf{h}^0 = \mathbf{0}, & \text{in } \mathcal{Y}, \\ 
\lim_{\xi_x \to \pm\infty} \mathbf{h}^0 = \mathbf{H}^0|_{0^{\pm}}, & \\ 
\mathbf{h}^0 \text{ continuous and } \xi_{\|} - \text{periodic}. & 
\end{cases} $$

We conclude that $\mathbf{h}^0$ is independent of $\xi$, and using the matching conditions (19), we obtain

$$ \mathbf{h}^0 = \mathbf{H}^0|_0, \quad D_x^0|_0 = D_x^0|_{0^{\pm}}, \qquad (31) $$

and we have used that the continuity of $\mathbf{H}^0$ implies the continuity of $D_x^0$. We also have that $\mathbf{rot}_{\xi} \mathbf{e}^0 = \mathbf{0}$ that we integrate over $\mathcal{Y}$ to get

$$ \lim_{\xi_x \to +\infty} \int_{\mathcal{Y}} (\mathbf{e}^0 \times \hat{\mathbf{x}}) d\xi_{\|} - \lim_{\xi_x \to -\infty} \int_{\mathcal{Y}} (\mathbf{e}^0 \times \hat{\mathbf{x}}) d\xi_{\|} = \mathbf{0}, \quad (32) $$

as the other contributions on $\partial\mathcal{Y}_m$ cancel by periodicity, hence

$$ \mathbf{E}_{\|}^{0}|_{0} = \mathbf{E}_{\|}^{0}|_{0^{\pm}}. $$

Expectedly in the limit of vanishing thickness, the structured film has a vanishing effect on the electromagnetic propagation and we have to go to the next order. Restoring the small but finite thickness of the metafilm is done as in the previous section using (22)-(23). It results the same jumps as in (24).

Note a noticeable difference with the case of a homogeneous films where we had $\mathbf{e}_{\|}^{0} = \mathbf{E}_{\|}^{0}|_{0}$ completely known. Here we have obtained that the macroscopic field $\mathbf{E}_{\|}^{0}$ is continuous but the microscopic field $\mathbf{e}^{0}(\xi, \mathbf{x}_{\|}, t)$ is still unknown.

2) **First-order problem:** The first, non trivial, problem is thus set on $(\mathbf{e}^0, \mathbf{h}^1)$ being the $\xi_{\|}$-periodic solutions of

$$ egin{cases} 
\partial_t \mathbf{e}^0 = \frac{1}{\varepsilon_0 \varepsilon_r(\xi)} (\mathbf{rot}_{\xi} \mathbf{h}^1 + \mathbf{rot}_{\|} \mathbf{H}^0|_0), \\ 
\mathbf{rot}_{\xi} \mathbf{e}^0 = \mathbf{0}, \quad \text{div}_{\xi} \mathbf{h}^1 + \text{div}_{\|} \mathbf{H}^0|_0 = 0 \quad \text{in } \mathcal{Y}, \\ 
\lim_{\xi_x \to \pm\infty} \partial_t \mathbf{e}^0 = \frac{1}{\varepsilon_0 \varepsilon_r^{\pm}} \partial_t D_x^0|_0 \hat{\mathbf{x}} + \partial_t \mathbf{E}_{\|}|_0, \\ 
\mathbf{h}^1 \text{ and } \mathbf{e}^0 \times \mathbf{n} \text{ continuous,} 
\end{cases} \qquad (33) $$

We have used that $\mathbf{h}^0 = \mathbf{H}^0|_0$ from (31) and the matching condition (19) on $\mathbf{e}^0$. The above problem is the analogue of the problem on $\mathbf{h}^1$ in (25) but without the simplification of a one-dimensional problem (in particular $\mathbf{e}^0$ is unknown in (33)).

Using further from (3) that (i) $\mathbf{rot}_{\|} \mathbf{H}^0 = \partial_t D_x^0 \hat{\mathbf{x}} - \mathbf{rot}_{\xi} (\xi_x \nabla_x H_x^0)$ and (ii) $\text{div}_{\|} \mathbf{H}^0 = -\text{div}_{\xi} (\xi_x \nabla_x H_x^0)$, we obtain that (33) is linear with respect to $\partial_t D_x^0|_0$, $\partial_t E_{\alpha}|_0$, $\alpha = y, z$, and the three components of $\nabla_x H_x^0|_0$. For these latter contributions, the dependences are explicit and identical to that in (25). Thus, we set

$$ egin{aligned} & \mathbf{h}^1 = \mathbf{h}^1|_{\xi=0} + \xi_x \nabla_x H_x^0|_0 + \partial_t D_x^0|_0 \mathbf{q}^{(x)} + \varepsilon_0 \partial_t E_{\alpha}^0|_0 \mathbf{q}^{(\alpha)}, \\ & \mathbf{e}^0 = \frac{D_x^0|_0}{\varepsilon_0 \varepsilon_r} (\mathbf{rot}_{\xi} \mathbf{q}^{(x)} + \hat{\mathbf{x}}) + \frac{E_{\alpha}^0|_0}{\varepsilon_r} \mathbf{rot}_{\xi} \mathbf{q}^{(\alpha)}, \end{aligned} \qquad (34) $$

where $\alpha = y, z$, repeated indices mean summation, and where the elementary fields $\mathbf{q}^{(x)}(\xi)$ and $\mathbf{q}^{(\alpha)}(\xi)$ are solutions to

$$ egin{cases} 
\mathbf{rot}_{\xi} \left( \frac{1}{\varepsilon_r} (\mathbf{rot}_{\xi} \mathbf{q}^{(x)} + \hat{\mathbf{x}}) \right) = \mathbf{0}, & \text{div}_{\xi} \mathbf{q}^{(x)} = \mathbf{0}, \\ 
\lim_{\xi_x \to \pm\infty} \mathbf{rot}_{\xi} \mathbf{q}^{(x)} = \mathbf{0}, & \\ 
\mathbf{q}^{(x)} \text{ and } \frac{1}{\varepsilon_r} (\mathbf{rot}_{\xi} \mathbf{q}^{(x)} + \hat{\mathbf{x}}) \times \mathbf{n} \text{ continuous,} & \\ 
\mathbf{q}^{(x)} \text{ and } \frac{1}{\varepsilon_r} (\mathbf{rot}_{\xi} \mathbf{q}^{(x)} + \hat{\mathbf{x}}) \xi_{\|} - \text{periodic}, & 
\end{cases} \qquad (35) $$

and

$$ egin{cases} 
\mathbf{rot}_{\xi} \left( \frac{1}{\varepsilon_r} \mathbf{rot}_{\xi} \mathbf{q}^{(\alpha)} \right) = \mathbf{0}, & \text{div}_{\xi} \mathbf{q}^{(\alpha)} = \mathbf{0}, \\ 
\lim_{\xi_x \to \pm\infty} \mathbf{rot}_{\xi} \mathbf{q}^{(\alpha)} = \varepsilon_r^{\pm} \hat{\alpha}, & \\ 
\mathbf{q}^{(\alpha)} \text{ and } \frac{1}{\varepsilon_r} \mathbf{rot}_{\xi} \mathbf{q}^{(\alpha)} \times \mathbf{n}, \text{ continuous,} & \\ 
\mathbf{q}^{(\alpha)} \text{ and } \frac{1}{\varepsilon_r} \mathbf{rot}_{\xi} \mathbf{q}^{(\alpha)} \xi_{\|} - \text{periodic}. & 
\end{cases} \qquad (36) $$

In (34), the $\mathbf{q}^{(i)}$, $i = x, y, z$, have been set equal to zero at $\xi = 0$ (as they are defined up to a constant). Next, $\mathbf{q}^{(x)}$ tends to constants when $\xi_x \to \pm\infty$ and the two constants differ in
general; the same applies to $(\mathbf{q}^{(y)} + \varepsilon_r^\pm \xi_x \hat{z})$ and $(\mathbf{q}^{(z)} - \varepsilon_r^\pm \xi_x \hat{z})$.
Accordingly, we define

$$
\begin{align}
\mathbf{b}^{(x)} &= \lim_{\xi_x \to +\infty} \mathbf{q}^{(x)} - \lim_{\xi_x \to -\infty} \mathbf{q}^{(x)}, \\
\mathbf{b}^{(y)} &= \lim_{\xi_x \to +\infty} (\mathbf{q}^{(y)} + \varepsilon_r^\pm \xi_x \hat{z}) - \lim_{\xi_x \to -\infty} (\mathbf{q}^{(y)} + \varepsilon_r^\pm \xi_x \hat{z}), \\
\mathbf{b}^{(z)} &= \lim_{\xi_x \to +\infty} (\mathbf{q}^{(z)} - \varepsilon_r^\pm \xi_x \hat{y}) - \lim_{\xi_x \to -\infty} (\mathbf{q}^{(z)} - \varepsilon_r^\pm \xi_x \hat{y}). \tag{37}
\end{align}
$$

Note that the $\mathbf{b}^{(i)}$ are in-plane vectors with

$$b_x^{(i)} = 0, \quad (38)$$

as $\mathbf{q}^{(i)}$ are divergentless hence $\int_y \mathrm{div}\xi \mathbf{q}^{(i)} = \mathbf{b}^{(i)} \cdot \hat{\mathbf{x}} = 0$.

We determine $\mathbf{H}^1|_{0\pm}$ starting with (34) and using that
$\nabla_x \mathbf{H}_x^0 - \partial_t \mathbf{D}_y^0 \hat{\mathbf{z}} + \partial_t \mathbf{D}_z^0 \hat{\mathbf{y}} = \partial_x \mathbf{H}^0$. Doing so, we get

$$
\lim_{\xi_x \to \pm\infty} \mathbf{h}^1 = \mathbf{h}^1|_{\xi=0} + \xi_x \partial_x \mathbf{H}^0|_{0\pm} + \partial_t D_x^0|_0 \lim_{\xi_x \to \pm\infty} \mathbf{q}^{(x)} \\
+ \varepsilon_0 \partial_t E_y^0|_0 \lim_{\xi_x \to \pm\infty} (\mathbf{q}^{(y)} + \varepsilon_r^\pm \xi_x \hat{z}) \\
+ \varepsilon_0 \partial_t E_z^0|_0 \lim_{\xi_x \to \pm\infty} (\mathbf{q}^{(z)} - \varepsilon_r^\pm \xi_x \hat{y}).
$$

In virtue of the matching conditions (20) on $\mathbf{H}^1$, we obtain
the jump in $\mathbf{H}^1$ of the form

$$[\delta \mathbf{H}^1] = \varepsilon_0 d \frac{\partial}{\partial t} \left( \mathbf{b}^{(x)} \frac{D_x^0}{\varepsilon_0} + \mathbf{b}^{(y)} E_y^0 + \mathbf{b}^{(z)} E_z^0 \right) \Big|_0. \quad (39)$$

We have yet to determine the jump of $(\mathbf{E}^1 \times \hat{\mathbf{x}})$. At the dominant order we have integrated $\mathbf{rot}_\xi \mathbf{e}^0 = 0$ in (32) to get the jump of $(\mathbf{E}^0 \times \hat{\mathbf{x}})$. We do the same but now, we have to integrate the form

$$\mathbf{rot}_\xi \mathbf{e}^1 + \mathbf{rot}_\| \mathbf{e}^0 + \mu_0 \partial_t \mathbf{H}^0 \Big|_0 = 0.$$

We shall use integrals over the domain $\mathcal{Y}_m = \{(\xi_x, \xi_\|) \in (-\xi_x^m, \xi_x^m) \times \mathcal{Y}\}$ as we shall manipulate integrals diverging when $\xi_x^m \to \infty$. The first integral is similar to (32). Using the matching conditions (20) written as

$$\mathbf{E}^1|_{0\pm} = \lim_{\xi_x^m \to +\infty} \left( \mathbf{e}^1_{\pm|\pm\xi_x^m} \mp \xi_x^m \partial_x \mathbf{E}^0|_{0\pm} \right),$$

and passing to the limit when $\xi_x^m \to \infty$, we obtain

$$
\lim_{\xi_x^m \to +\infty} \int_{\mathcal{Y}_m} \mathbf{rot}_\xi \mathbf{e}^1 = -(\mathbf{E}^1|_{0+} - \mathbf{E}^1|_{0-}) \times \hat{\mathbf{x}} \\
- \lim_{\xi_x^m \to +\infty} \xi_x^m (\partial_x \mathbf{E}^0|_{0+} + \partial_x \mathbf{E}^0|_{0-}) \times \hat{\mathbf{x}}. \tag{40}
$$

The second integral can be evaluated thanks to the form of $\mathbf{e}^0$ in (34) resulting in

$$
\int_{\mathcal{Y}_m} \mathbf{rot}_\| \mathbf{e}^0 d\xi = \xi_x^m \mathbf{rot}_\| (\mathbf{E}^0|_{0-} + \mathbf{E}^0|_{0+}) \\
+ \mathbf{rot}_\| \left( (\mathbf{c}^{(x)} + C_x \hat{\mathbf{x}}) \frac{D_x^0}{\varepsilon_0} |_0 + (\mathbf{c}^{(\alpha)} + C_y \hat{\alpha}) \mathbf{E}_\alpha^0 |_0 \right). \quad (41)
$$

We have used that $\int_{\mathcal{Y}_m} \frac{d\xi}{\varepsilon_r} = \xi_x^m \left(\frac{1}{\varepsilon_r^+} + \frac{1}{\varepsilon_r^-}\right) + C_x$, and $\int_{\mathcal{Y}_m^-} \frac{\varepsilon_r^j d\xi}{\varepsilon_r} = 2\xi_x^m + C_y$, with

$$
\begin{align}
C_x &= e_\xi \langle \varepsilon_r^{-1} \rangle - e_\xi^\dagger \varepsilon_r - e_\xi^\dagger \xi_x^\dagger \varepsilon_r^\dagger, \\
C_y &= e_\xi^\dagger \varphi^- \left( \varepsilon_r^- \varepsilon_r^\dagger - 1 \right) + e_\xi^\dagger \varphi^+ \left( \varepsilon_r^\dagger \varepsilon_r^- - 1 \right). \tag{42}
\end{align}
$$

where $e_\xi^\dagger = \frac{e_\xi}{d}$, $e_\xi = \frac{e}{d}$ and we have defined

$$
\begin{align}
\mathbf{c}^{(x)} &= \int_{\mathcal{Y}} \frac{1}{\varepsilon_r} \mathbf{rot}_\xi \mathbf{q}^{(x)} d\xi. \\
\mathbf{c}^{(\alpha)} &= \int_{\mathcal{Y}^+} \frac{1}{\varepsilon_r} (\mathbf{rot}_\xi \mathbf{q}^{(\alpha)} - \varepsilon_r^\dagger \hat{\alpha}) d\xi + \int_{\mathcal{Y}^+} \frac{1}{\varepsilon_r} (\mathbf{rot}_\xi \mathbf{q}^{(\alpha)} - \varepsilon_r^\dagger \hat{\alpha}) d\xi. \tag{43}
\end{align}
$$

Eventually the third integral simply reads

$$\mu_0 \int_{\mathcal{Y}_m} \partial_t \mathbf{H}^0|_0 d\xi = 2\xi_x^m \mu_0 \partial_t \mathbf{H}^0|_0. \quad (44)$$

Summing the contributions (40), (41) and (44) and using (17), the terms linear in $\xi_x^m$ cancel and we obtain

$$ (\mathbf{E}^1|_{0+} - \mathbf{E}^1|_{0-}) \times \hat{\mathbf{x}} = \frac{1}{\varepsilon_0} \mathbf{rot}_\| ((\mathbf{c}^{(x)} + C_x \hat{\mathbf{x}}) D_x^0|_0) + \mathbf{rot}_\| ((\mathbf{c}^{(\alpha)} + C_y \hat{\alpha}) \mathbf{E}_\alpha^0|_0). $$

It is now sufficient to use the properties (61) and (66) (the proofs are given in the appendix A) which tell us that $\mathbf{c}^{(x)} = c_x^{(x)} \hat{\mathbf{x}}$ and $\mathbf{c}^{(\alpha)} + C_y \hat{\alpha} = c_x^{(\alpha)} \hat{\mathbf{x}}$, $\alpha = y, z$ with $c_x^{(y)} = b_x^{(x)}$ and $c_x^{(z)} = -b_y^{(x)}$. As $\mathbf{rot}_\| (D_x^0 \hat{\mathbf{x}}) = \nabla_\| D_x^0 \times \hat{\mathbf{x}}$, the above jumps read

$$[\delta \mathbf{E}^1 \times \hat{\mathbf{x}}] = d \nabla_\| \left( (c_x^{(x)} + C_x) \frac{D_x^0}{\varepsilon_0} + b_z^{(x)} E_y^0|_0 - b_y^{(x)} E_z^0|_0 \right) \times \hat{\mathbf{x}}. \quad (45)$$

Gathering the contributions at the zero-order in (24) and that at the first-order in (39) and (45), we eventually obtain the transition conditions which takes the form of (5) along with (6) and (8) namely

$$
\left\{
\begin{array}{l}
[\mathbf{H} \times \hat{\mathbf{x}}] = -\varepsilon_0 \frac{\partial}{\partial t} \left( \chi_{ee}^{\alpha x} \frac{D_{x,av}}{\varepsilon_0} \hat{\alpha} + \chi_{ee}^{\alpha \beta} E_{\beta,av} \hat{\alpha} \right) \\
\qquad + e \nabla_\| H_{x,av} \times \hat{\mathbf{x}}, \\[1em]
[\mathbf{E} \times \hat{\mathbf{x}}] = \mu_0 e \partial_t \mathbf{H}_\| - \nabla_\| \left( \chi_{ee}^{xx} \frac{D_{x,av}}{\varepsilon_0} + \chi_{ee}^{x\alpha} E_{\alpha,av} \right) \times \hat{\mathbf{x}}
\end{array}
\right.
\tag{46}
$$

where the effective susceptibility $\overline{\chi}_{ee}$ is symmetric with diagonal terms

$$
\begin{align}
\chi_{ee}^{xx} &= -(e \langle \varepsilon_r^{-1} \rangle + d c_x^{(x)}), \nonumber \\
\chi_{ee}^{yy} &= (e^{-\varepsilon_r^{-}} + e^{+\varepsilon_r^{+}} - db_z^{(y)}), \quad \chi_{ee}^{zz} = (e^{-\varepsilon_r^{-}} + e^{+\varepsilon_r^{+}} + db_y^{(z)}). \tag{47}
\end{align}
$$

and off-diagonal terms

$$\chi_{ee}^{xy} = -db_z^{(x)}, \quad \chi_{ee}^{xz} = db_y^{(x)}, \quad \chi_{ee}^{yz} = db_y^{(y)}. \quad (48)$$
### C. Energetic properties of the effective model

The conservation of the energy is written starting from (3) and multiplying the first equation by $\mathbf{E}$, the second by $\mathbf{H}$ and summing. We obtain

$$ \frac{d}{dt} \mathcal{E} + \int_{\partial \Omega_e} \mathbf{S} \cdot \mathbf{n} \, dS = 0. $$

In the effective problem, $\partial\Omega_e = \Sigma_e \cup \Gamma$ contains the boundary $\Gamma$ across which $\mathbf{H} \times \hat{\mathbf{x}}$ and $\mathbf{E} \times \hat{\mathbf{x}}$ are not continuous, hence they produce a contribution to the flux of the Poynting vector resulting in

$$ \frac{d}{dt} \mathcal{E} - \int_{\Gamma} \left[ \left[ \mathbf{S} \cdot \hat{\mathbf{x}} \right] \, d\xi_{\parallel} + \int_{\Sigma_e} \mathbf{S} \cdot \mathbf{n} \, dS \right] = 0. $$

The flux through $\Gamma = \{\Gamma^{-}, \Gamma^{+}\}$ reads

$$ \Phi_{\Gamma} = \int_{\Gamma} \left( \left[ \mathbf{E} \times \hat{\mathbf{x}} \right] \cdot \mathbf{H}_{av} - \left[ \mathbf{H} \times \hat{\mathbf{x}} \right] \cdot \tilde{\mathbf{E}}_{av} \right) \, d\mathbf{x}_{\parallel}. $$

The first contribution $I_1 = \int_{\Gamma} \left[ \mathbf{E} \times \hat{\mathbf{x}} \right] \cdot \mathbf{H}_{av} \, d\mathbf{x}_{\parallel}$ reads

$$ I_1 = \frac{1}{2} \frac{d}{dt} \int_{\Gamma} \left( \mu_0 e \left| \mathbf{H}_{\parallel,av} \right|^2 - \chi_{ee}^{xx} \frac{D_{x,av}^2}{\varepsilon_0} \right) d\mathbf{x}_{\parallel} \\ - \int_{\Gamma} \partial_t D_{x,av} (\chi_{ee}^{yy} E_{y,av} + \chi_{ee}^{xz} E_{z,av}) + \text{b.t.} $$

We have used (46)-(48) and we have integrated by parts the terms $(\nabla_{\parallel} E_{\alpha,av} \times \hat{\mathbf{x}}) \cdot \mathbf{H}_{av}$ along with (3). It is worth noting that the integrations by parts make boundary terms (b.t.) appear which depend on the conditions at the ending points of $\Gamma$ along $x$; they are disregarded in the present study. Doing the same with the second contribution $I_2 = - \int_{\Gamma} \left[ \mathbf{H} \times \hat{\mathbf{x}} \right] \cdot \mathbf{E}_{av} d\mathbf{x}_{\parallel}$, we obtain

$$ I_2 = +\frac{1}{2} \frac{d}{dt} \int_{\Gamma} \mu_0 e H_{x,av}^2 \\ +\frac{1}{2} \frac{d}{dt} \int_{\Gamma} \varepsilon_0 (\chi_{ee}^{yy} E_{y,av}^2 + \chi_{ee}^{zz} E_{z,av}^2 + 2\chi_{ee}^{yz} E_{y,av} E_{z,av}) d\mathbf{x}_{\parallel} \\ + \int_{\Gamma} \partial_t D_{x,av} (\chi_{ee}^{xy} E_{y,av} + \chi_{ee}^{xz} E_{z,av}) + \text{b.t.} $$

Eventually, summing the two contributions $I_1, I_2$, we obtain the interfacial energy (11).

1) Positiveness of the interfacial energy: In this section, we shall prove that $\mathcal{E}_{\Gamma}$ in (11) is definite positive, namely we shall prove that

(i) $-\chi_{ee}^{xx} > 0,$

(ii) $\forall (E_y, E_z) \in (\mathbb{R}^*)^2,$

$$ \chi_{ee}^{yy} E_y^2 + \chi_{ee}^{zz} E_z^2 + 2\chi_{ee}^{yz} E_y E_z > 0, $$

with $\bar{\chi}_{ee}$ given in (47)-(48). More specifically, we shall establish bounds for $-\chi_{ee}^{xx}$ and bounds for the quadratic form $(\chi_{ee}^{yy} E_y^2 + \chi_{ee}^{zz} E_z^2 + 2\chi_{ee}^{yz} E_y E_z)$ from which (13)-(14) are deduced. To do so, we use minimization principles [39] which allow to establish bounds on the effective parameters entering the energy.

• To show (i), we use a minimization principle set on $\mathbf{T}^{(x)} = \frac{1}{\varepsilon_r} (\mathbf{rot}_{\xi} \mathbf{q}^{(x)} + \hat{\mathbf{x}})$ which is associated with the elementary problem (35). The minimization principle reads

$$ \begin{array}{l} \forall \tilde{\mathbf{T}} \in \mathcal{T}, \quad \mathcal{E}^*(\tilde{\mathbf{T}}) \ge \mathcal{E}^*(\mathbf{T}^{(x)}), \\ \text{with } \mathcal{E}^*(\tilde{\mathbf{T}}) = \frac{1}{2} \int_{\mathcal{Y}} \varepsilon_r \left| \tilde{\mathbf{T}} - \frac{\hat{\mathbf{x}}}{\varepsilon_r} \right|^2 d\xi, \end{array} \quad (49) $$

and where $\mathcal{T}$ is the set of admissible fields $\tilde{\mathbf{T}}$ being $\xi_{\parallel}$-periodic with $\tilde{\mathbf{T}} \times \mathbf{n}$ continuous and satisfying $\mathbf{rot}_{\xi} \tilde{\mathbf{T}} = \mathbf{0}$ and $\lim_{\xi_x \to \pm\infty} \tilde{\mathbf{T}} = \frac{\hat{\mathbf{x}}}{\varepsilon_r^\pm}$. By definition, we have

$$ \mathcal{E}^*(\mathbf{T}^{(x)}) = \frac{1}{2} \int_{\mathcal{Y}} \frac{1}{\varepsilon_r} |\mathbf{rot}_{\xi} \mathbf{q}^{(x)}|^2 d\xi = -\frac{c_x^{(x)}}{2}. \quad (50) $$

where we have used (65)) for $i = j = x$. This already tells us that $c_x^{(x)} \le 0$ (as $\mathcal{E}^*(\mathbf{T}^{(x)}) \ge 0$). Now, we choose an admissible field $\tilde{\mathbf{T}} \in \mathcal{T}$ of the form

$$ \tilde{\mathbf{T}} = \hat{\mathbf{x}} \begin{cases} \frac{1}{\varepsilon_r^+}, & \xi_x < -e_{\xi}^-, \\ \tau^-, & -e_{\xi}^- < \xi_x < 0, \\ \tau^+, & 0 < \xi_x < e_{\xi}^+, \\ \frac{1}{\varepsilon_r^+}, & \xi_x > e_{\xi}^+ \end{cases} \quad (51) $$

and at this stage, $\tau^+$ and $\tau^-$ are free, real parameters. Inserting (51) in (49), we obtain

$$ \mathcal{E}^*(\tilde{\mathbf{T}}) = \frac{1}{2} \left( e_{\xi}^{-} (\langle \varepsilon_{\tau} \rangle^{-} \tau^{-2} - 2\tau^{-} + \langle \varepsilon_{\tau}^{-1} \rangle^{-}) + e_{\xi}^{+} (\langle \varepsilon_{\tau} \rangle^{+} \tau^{+2} - 2\tau^{+} + \langle \varepsilon_{\tau}^{-1} \rangle^{+}) \right). $$

The minimum of $\mathcal{E}^*(\tilde{\mathbf{T}})$ with respect to $\tau^+$ and $\tau^-$ is reached for $\tau^\pm = \frac{1}{\langle \varepsilon_\tau \rangle^\pm}$ resulting in

$$ c_x^{(x)} \ge -2\mathcal{E}^*(\tilde{\mathbf{T}}) = \frac{e_{\xi}^{-}}{\langle \varepsilon_{\tau} \rangle^{-}} + \frac{e_{\xi}^{+}}{\langle \varepsilon_{\tau} \rangle^{+}} - e_{\xi} \langle \varepsilon_{\tau}^{-1} \rangle. $$

With $-\chi_{ee}^{xx} = dc_x^{(x)} + e\langle\varepsilon_r^{-1}\rangle$ from (47), and $e^\pm = de_\xi^\pm$, $e = de_\xi$ by definition, we get the two bounds on $-\chi_{ee}^{xx}$ announced in (13).

• To show (ii), we use a minimization principle set on

$$ \mathbf{q}(\xi) = E_y (\mathbf{q}^{(y)}(\xi) + Q(\xi_x)\xi_x \hat{\mathbf{z}}) + E_z (\mathbf{q}^{(z)}(\xi) - Q(\xi_x)\xi_x \hat{\mathbf{y}}), $$

for $(E_y, E_z) \in \mathbb{R}^2$ (see $^3$), and where $Q(\xi_x < 0) = \varepsilon_r^-$ and $Q(\xi_x > 0) = \varepsilon_r^+$. With $\mathbf{q}^{(\alpha)}$, $\alpha = y, z$ solutions to (36), $\mathbf{q}$ satisfies

$$ \left\{ \begin{array}{l} \mathbf{rot}_{\xi} \left( \frac{1}{\varepsilon_r} (\mathbf{rot}_{\xi} \mathbf{q} + Q\mathbf{E}_{\parallel}) \right) = \mathbf{0}, \quad \mathrm{div}_{\xi} \mathbf{q} = \mathbf{0}, \\ \lim_{\xi_x \to \pm\infty} \mathbf{rot}_{\xi} \mathbf{q} = \mathbf{0}, \\ \mathbf{q} \text{ and } \frac{1}{\varepsilon_r} (\mathbf{rot}_{\xi} \mathbf{q} + Q\mathbf{E}_{\parallel}) \times \mathbf{n}, \quad \text{continuous,} \\ \mathbf{q} \text{ and } \frac{1}{\varepsilon_r} \mathbf{rot}_{\xi} \mathbf{q} \quad \mathbf{\xi}_{\parallel}\text{-periodic.} \end{array} \right. \quad (52) $$

<footnote>$^3$Note that $(E_y, E_z)$ can be considered as real parameters as they do not depend on the microscopic coordinates.</footnote>
The minimization principle reads

$$ \forall \tilde{q} \in \mathcal{Q}, \quad \mathcal{E}(\tilde{q}) \ge \mathcal{E}(\mathbf{q}), $$

$$ \mathcal{E}(\tilde{q}) = \int_{\mathcal{Q}} \left( \frac{|\operatorname{rot}\xi \tilde{q}|^2}{2\varepsilon_r} + \frac{Q}{\varepsilon_r} \mathbf{E}_{\parallel} \cdot \operatorname{rot}\xi \tilde{q} \right) d\xi + E_y \tilde{b}_z - E_z \tilde{b}_y, \quad (53) $$

with $\mathcal{Q}$ the set of admissible fields $\tilde{q}$ being continuous and $\xi_{\parallel}$-periodic such that $\operatorname{div}\xi\tilde{q} = 0$; as in (37), we have defined $\tilde{b} = \lim_{\xi_x \to +\infty} \tilde{q} - \lim_{\xi_x \to -\infty} \tilde{q}$. To determine $\mathcal{E}(\mathbf{q})$ we use (52) along with (47)-(48). We obtain that $E(\mathbf{q})$ involves the quadratic form Q.F. is $\chi_{ee}^{yy} E_y^2 + \chi_{ee}^{zz} E_z^2 + 2\chi_{ee}^{yz} E_y E_z$ as it reads

$$ \mathcal{E}(\mathbf{q}) = -\frac{\mathrm{Q.F.}}{d} + \frac{|\mathbf{E}_{\parallel}|^2}{2} e_{\xi}^{\tau} \varepsilon_{\tau}^{-} (2 - \varepsilon_{\tau}^{-} \langle \varepsilon_{\tau}^{-1} \rangle^{-}) + \frac{|\mathbf{E}_{\parallel}|^2}{2} e_{\xi}^{\tau} \varepsilon_{\tau}^{+} (2 - \varepsilon_{\tau}^{+} \langle \varepsilon_{\tau}^{-1} \rangle^{+}). \quad (54) $$

We now choose an admissible field of the form

$$ \tilde{q} = (E_y \hat{z} - E_z \hat{y}) \begin{cases} e_{\xi} \tau^{-}, & \xi_x < -e_{\xi}^{-}, \\ \xi_x \tau^{-}, & -e_{\xi}^{-} < \xi_x < 0, \\ \xi_x \tau^{+}, & 0 < \xi_x < e_{\xi}^{+}, \\ e_{\xi}^{\tau}, & \xi_x > e_{\xi}^{+}, \end{cases} $$

which is such that $\operatorname{rot}\xi\tilde{q} = -\mathbf{E}_{\parallel}\tau^{\pm}$ in $\xi_x \in (-e_{\xi}^{-}, e_{\xi}^{+})$ and zero otherwise. As previously, $\tau^{\pm}$ are free parameters at this stage which will be chosen in order to minimize the resulting energy $\mathcal{E}(\tilde{q})$. Doing so, we obtain $\tau^{\pm} = \varepsilon_{\tau}^{\pm} - 1/(\varepsilon_{\tau}^{-1})^{\pm}$ and the resulting energy reads as

$$ \mathcal{E}(\tilde{q}) = -\frac{|\mathbf{E}_{\parallel}|^2}{2} \left( \frac{e_{\xi}^{-}}{\langle \varepsilon_{\tau}^{-1} \rangle^{-}} + \frac{e_{\xi}^{+}}{\langle \varepsilon_{\tau}^{-1} \rangle^{+}} \right) + \frac{|\mathbf{E}_{\parallel}|^2}{2} \left( e_{\xi}^{-} \varepsilon_{\tau}^{-} (2 - \varepsilon_{\tau}^{-} \langle \varepsilon_{\tau}^{-1} \rangle^{-}) + e_{\xi}^{+} \varepsilon_{\tau}^{+} (2 - \varepsilon_{\tau}^{+} \langle \varepsilon_{\tau}^{-1} \rangle^{+}) \right). \quad (55) $$

Reporting (54) and (55) in (53), we obtain a lower bound for the quadratic form

$$ (\chi_{ee}^{yy} E_{y,av}^2 + \chi_{ee}^{zz} E_{z,av}^2 + 2\chi_{ee}^{yz} E_{y,av} E_{z,av}) \ge \chi_{\perp}^{m}(E_x^2 + E_y)^2 > 0, \quad (56) $$

which already tells us that the quadratic form is positive. (13)-(14). The upper bound is obtained by using the dual problem as used to show (i) and we give below the main steps of the calculations. The minimization principle reads

$$ \forall \tilde{T} \in \mathcal{T}, \quad \mathcal{E}^*(\tilde{T}) \ge \mathcal{E}^*(T), $$
with $\mathcal{E}^*(\tilde{T}) = \frac{1}{2} \int_{\mathcal{Q}} \left| \tilde{T} - \frac{Q}{\varepsilon_r} \mathbf{E}_{\parallel} \right|^2, \quad (57) $$

$\mathcal{T}$ being the set of admissible fields $\tilde{T}$ being $\xi_{\parallel}$-periodic, $\tilde{T} \times n$ continuous and satisfying $\operatorname{rot}\xi\tilde{T} = 0$ and $\lim_{\xi_x \to \pm\infty} \tilde{T} = \mathbf{E}_{\parallel}$ The field $\mathbf{T}$ is defined by $\mathbf{T} = \frac{1}{\varepsilon_r}(\operatorname{rot}\xi\mathbf{q} + Q\mathbf{E}_{\parallel})$ and $\mathcal{E}^*(\mathbf{T}) = -\mathcal{E}(\mathbf{q})$ with $\mathcal{E}(\mathbf{q})$ given in (54). We choose $\tilde{T} = \mathbf{E}_{\parallel}$ which is an admissible field resulting in $\mathcal{E}(\tilde{T}) = \frac{|\mathbf{E}_{\parallel}|^2}{2} \left( e_{\xi}^{-} \varphi^{-} \varepsilon_{\tau}^{sc} \left(1 - \frac{e_{\xi}^{-}}{\varepsilon_{\tau}^{sc}}\right)^2 + e_{\xi}^{+} \varphi^{+} \varepsilon_{\tau}^{sc} \left(1 - \frac{e_{\xi}^{+}}{\varepsilon_{\tau}^{sc}}\right)^2 \right)$ which along (54) provides the upper bound

$$ (\chi_{ee}^{yy} E_{y,av}^2 + \chi_{ee}^{zz} E_{z,av}^2 + 2\chi_{ee}^{yz} E_{y,av} E_{z,av}) \ge \chi_{\perp}^{m}(E_x^2 + E_y)^2 > 0. \quad (58) $$

The bounds for the quadratic form in (56) and (58) are equivalent to those announced in (13)-(14).

## IV. VALIDATION OF THE EFFECTIVE MODEL

In this section, we provide illustrations of the efficiency of the effective model by means of comparisons with direct numerics in the harmonic regime in two- and three-dimensions. For the actual and the effective problems, the Finite Element Method (FEM) has been used in Comsol Multiphysics; the GSTCs along with the resolution of the elementary problems have been implemented (see supplementary material). In the numerics, we have used a cylindrical electric source of current density $\mathbf{J}_e = (0, -1, 1)$ and we considered a normalized value of the wavelength $\lambda = 1$. Eventually, we imposed PML at the ends of the computation domain along $x$ and periodic boundary conditions in the $y$ and $z$ coordinates.

### A. Two-dimensional case

We consider two-dimensional geometries that is geometries invariant along the $z$-direction, with substrates of silica ($\varepsilon_r^{-} = 1.45^2$) and air ($\varepsilon_r^{+} = 1$). The reference case is a simple interface

<figure><img src="image_7.png" /><figcaption>Fig. 6: Geometries of the layers across which the GSTCs apply in the effective problem (see figure 3).</figcaption></figure>

between the silica and the air (figure 6(a)). To avoid a trivial effective problem with $e = 0$ (which would be identical to the actual problem), we choose an excluded region of thickness $e = 0.1\lambda$ (with $1/4$ of silica and $3/4$ of air) across which the GSCTs apply. The susceptibility tensor has been calculated numerically, using (47)-(48) with the unknown coefficients ($c_x^{(x)}, b_z^{(y)}, b_y^{(z)}, b_y^{(x)}, b_y^{(y)}$) given by (37) and (43) (in practice we use the equivalent forms (69)-(71) with scalar potentials, see also the supplementary material). We obtain

$$ \text{interface silica/air, } \quad \bar{\bar{\chi}}_{ee} = \begin{pmatrix} -0.087 & 0 & 0 \\ 0 & 0.127 & 0 \\ 0 & 0 & 0.127 \end{pmatrix}. \quad (59) $$

This is conform with (28) (in particular the tensor is diagonal) and with (13). Indeed, as $\langle \varepsilon_{\tau} \rangle^{\pm} = \varepsilon_{\tau}^{\pm}$ and $\langle \varepsilon_{\tau}^{-1} \rangle^{\pm} = 1/\varepsilon_{\tau}^{\pm}$, the bounds are equal with $\chi_{ee}^{xx} = \chi_{\parallel}^{m} = \chi_{\parallel}^{M} \simeq -0.087 < 0$ and $\chi_{ee}^{\alpha\alpha} = \chi_{\perp}^{m} = \chi_{\perp}^{M} \simeq 0.127 > 0$, $\alpha = y, z$.

As a small step further, we consider a uniform film of silicon ($\varepsilon_{\tau}^{sc} = 3.5^2$) of thickness $\lambda/20$ between the silica and the air. The minimum thickness of the excluded layer is $\lambda/20$ but we choose the same thickness $e = \lambda/10$ as for the interface with
<figure><img src="image_8.png" /><figcaption>Fig. 7: Real part of $E_y$ for (a) a single interface between silica and air and (b) a silicon film between silica and air. The field calculated in the actual problem is shown for $y < 0$ and the field calculated with the GSTCs is shown for $y > 0$.</figcaption></figure>

1/4 of silica, 1/4 of air and 1/2 of silicon (figure 6(b)). The susceptibility tensor calculated numerically now reads

$$ \text{silicon layer, } \bar{\bar{\chi}}_{ee} = \begin{pmatrix} -0.041 & 0 & 0 \\ 0 & 0.690 & 0 \\ 0 & 0 & 0.690 \end{pmatrix}, \quad (60) $$

which is again conform with (28) and (13) since $-\chi_{\parallel}^{M} = -\frac{e}{4}(\frac{1}{\varepsilon_r^+} + \frac{1}{\varepsilon_r^-} + \frac{2}{\varepsilon_r^{sc}}) \approx -0.041$ and $\chi_{\perp}^{M} = -\chi_{\perp}^{M} = -\frac{e}{4}(\varepsilon_r^+ + \varepsilon_r^- + 2\varepsilon_r^{sc}) \approx 0.690$.

We report in figure 7 the real part of the electric field $E_y(x, y)$ for $(x, y) \in (-3.5, 2.5)\lambda \times (-3, 3)\lambda$ (the source is sketched by the black circle). The comparison between the actual and the effective solutions is done by representing for $y < 0$ the electric field corresponding to the simulation of the actual problem (1) and for $y > 0$ the field corresponding to the simulation of the effective problem (3) with the GSTC (5)-(7) and using (59) for the interface ans (60) for the film. In both cases, we observe a good agreement between the results of the direct and effective problems.

Finally, we consider a microstructured array made of a periodic arrangement of clover-like silicon particles (figure 6(c)). The actual spacing is $d_y = \lambda/10$ and the particle is centred in the cell resulting in $e^- = e^+ = \lambda/20$. Due to the symmetry of the particle, the tensor is found to be diagonal (see appendix C) and it reads

$$ \text{clover-like particle } \bar{\bar{\chi}}_{ee} = \begin{pmatrix} -0.046 & 0 & 0 \\ 0 & 0.312 & 0 \\ 0 & 0 & 0.604 \end{pmatrix}. $$

Again, and as it should be, the tensor satisfies the bounds in (13). Indeed, with $\varphi^- = \varphi^+ = 0.42$ the fractions of particle, we have $\langle\varepsilon_r^-\rangle = 6.364$, $\langle\varepsilon_r^+\rangle = 5.725$, $\langle\varepsilon_r\rangle = 6.044$ ($\langle\varepsilon_r^{-1}\rangle = 0.310$, $\langle\varepsilon_r^{+1}\rangle = 0.614$ and $\langle\varepsilon_r^{-1}\rangle = 0.462$, from which

$$ -0.046 \le \chi_{ee}^{xx} \le -0.017, \quad 0.242 \le \chi_{ee}^{\alpha\alpha} \le 0.604. $$

The result for $E_y$ is reported in figure 8. Expectedly the scattering pattern is intermediate between the interface with no silicon particles and the silicon film, and its features are well reproduced by the effective model.

<figure><img src="image_9.png" /><figcaption>Fig. 8: Real part of $E_y$ for a layer comprising an array of clover-like silicon scattering particles. Same representation as in figure 7.</figcaption></figure>

## B. Three-dimensional case

We now move on to three-dimensional configurations with an array of cylindrical- and cuboid-shaped silicon particles with the air as unique substrate (figure 9). As the computational cost is higher in three-dimensions we have reduced the size of the domain which now contains a layer of $20 \times 20$ inclusions. With $d_y = d_z = e = \lambda/10$ this corresponds to $(x, y, z) \in (-2, 2)\lambda \times ((-1, 1)\lambda)^2$.

<figure><img src="image_10.png" /><figcaption>Fig. 9: Geometries of the excluded layer for silicon particles with (a) cylindrical shape and (b) cuboid shape rotated of 30° in the $(y, z)$.</figcaption></figure>

The cylindrical particles have dimensions $h = 0.5e$ and $R = 0.3e$, hence a filling fraction $\varphi = 0.14$. The susceptibility tensor calculated numerically reads

$$ \text{cylindrical particles, } \bar{\bar{\chi}}_{ee} = \begin{pmatrix} -0.072 & 0 & 0 \\ 0 & 0.142 & 0 \\ 0 & 0 & 0.142 \end{pmatrix}, $$

which satisfies the bounds (13), namely

$$ -0.087 \le \chi_{ee}^{xx} \le -0.039, \quad 0.115 \le \chi_{ee}^{\alpha\alpha} \le 0.259, $$

(with $\langle\varepsilon_r\rangle = \varphi\varepsilon_r^{sc} + (1-\varphi) \approx 2.59$, $\langle\varepsilon_r^{-1}\rangle = \varphi/\varepsilon_r^{sc} + (1-\varphi) \approx 0.87$). The cuboid particles have dimensions $h = 0.5e$ and $h' = 0.8e$ hence a filling fraction $\varphi = 0.20$. As the
particles are rotated by an angle of 30° in the (y, z)-plane, the susceptibility tensor now involves non-zero $\chi_{ee}^{yz}$ and it reads

$$ \text{cuboid particles, } \bar{\chi}_{ee} = \begin{pmatrix} -0.066 & 0 & 0 \\ 0 & 0.196 & 0.011 \\ 0 & 0.011 & 0.165 \end{pmatrix}. $$

For the diagonal terms, we have from (13)

$$ -0.082 \le \chi_{ee}^{xx} \le -0.031, \quad 0.123 \le \chi_{ee}^{\alpha\alpha} \le 0.325, $$

(with $\langle \epsilon_r \rangle \approx 3.25$, $\langle \epsilon_r^{-1} \rangle \approx 0.82$). Next, the off-diagonal term $\chi_{ee}^{yz}$ satisfies (14) which for given $\chi_{ee}^{xx}$ and $\chi_{ee}^{yy}$ reads

$$ \chi_{ee}^{yz} < \min(0.056, 0.144). $$

<figure><img src="image_11.png" /><figcaption>Fig. 10: Real part of $E_y$ in the effective problem (top) and in the actual problem (bottom) for the cylindrical particles.</figcaption></figure>

The real parts of $E_y$ calculated with the actual microstructure and with the GSTCs are reported in figures 10 and 11. The array of cylindrical particles has a slightly weaker scattering strength than that of cuboid particles, which is attributable to a slightly lower filling fraction. The difference is however small which results in similar scattered fields, and it is worth noticing that this difference is well captured by the effective model.

## V. CONCLUSION

We have derived an effective, homogenized, model leading to GSTCs formulated across an excluded (enlarged) region. The model has two advantages. Being associated to a positive energy, it is unconditionally stable in time with diagonal terms of the permittivity tensor for which bounds have been established. These theoretical bounds are of prime interest as they give the ranges of the susceptibilities that can be obtained with real structures. Next the whole tensor can be calculated without using retrieval method by means of the

<figure><img src="image_12.png" /><figcaption>Fig. 11: Real part of $E_y$ for the cuboid particles. Same representation as in figure 10.</figcaption></figure>

resolutions of simple scalar static problems. If optimization of a metafilm for a given objective is sought, such a model without any adjustable terms is useful as it offers a one-to-one correspondence between the actual microstructure and the effective susceptibility tensor. The enlarged GSTCs have been implemented numerically which allows us to exemplify the accuracy of the model in two and three dimensions. To our knowledge, such a validation has not been proposed in the previous literature.

Extensions of the present study include the case of magnetodielectric materials which is straightforward and will produce a non-trivial magnetic susceptibility tensor being the counterpart of the electric susceptibility tensor. In the same spirit, accounting for metals with small skin depths should be straightforward. More involved extensions concern arrays with local resonances of the Mie type which requires to adapt the asymptotic analysis.

## ACKNOWLEDGMENT

N.L. acknowledge support from French defence procurement agency under the ANR ASTRID Maturation program, grant agreement number ANR-18-ASMA-0006-01. K.P. acknowledge support from the Agence de l'Innovation de Défense (AID) from the Direction Générale de l'Armement (DGA), under grant no. 2019 65 0070 and the Agence Nationale de la Recherche under grant ANR-19-CE08-0006.

## APPENDIX

### A. Properties of the some effective coefficients

We shall use in this section that $\text{rot}_\xi(\xi_x \hat{y}) = \hat{z}$ and $\text{rot}_\xi(\xi_x \hat{z}) = -\hat{y}$; also the identity $\text{div}(\mathbf{A} \times \mathbf{B}) = \text{rot}\mathbf{A}$.
**B** – **A** · **rotB**. When needed, we shall use the domain $\mathcal{Y}_m = \{(\xi_x, \xi_\parallel) \in (-\xi_x^m, \xi_x^m) \times \mathcal{Y}\}$ (as we shall manipulate integrals diverging when $\xi_x^m \to \infty$).

**1) Explicit expressions:** We start by establishing the explicit form of some effective parameters ($\mathbf{b}^{(i)}$, $\mathbf{c}^{(i)}$), $i = x, y, z$, as well as relations between some of them. Specifically we show below that

$$c_{\alpha}^{(x)} = 0, \quad c_{z}^{(\alpha)} = -C_{y}\delta_{\alpha z}, \quad c_{y}^{(\alpha)} = -C_{y}\delta_{\alpha y}, \quad \alpha = y, z, \qquad (61)$$

where $C_y$ is defined in (42).

We consider the vectors $(\xi_x\hat{\beta})$, $\beta = y, z$ being continuous in $\mathcal{Y}$, $\xi_\parallel$-periodic, and such that $\mathrm{rot}(\xi_x\hat{\beta}) \cdot \hat{x} = 0$. We use from (35) that

$$0 = \int_{\mathcal{Y}} (\xi_x \hat{\beta}) \cdot \mathrm{rot}_{\xi} \left( \frac{1}{\varepsilon_r} (\mathrm{rot}_{\xi} \mathbf{q}^{(x)} + \hat{x}) \right) d\xi \\ = \int_{\mathcal{Y}} \frac{1}{\varepsilon_r} \mathrm{rot}_{\xi} (\xi_x \hat{\beta}) \cdot \mathrm{rot}_{\xi} \mathbf{q}^{(x)} d\xi.$$

We have used that the contributions of the terms $\frac{1}{\varepsilon_r}(\xi_x\hat{\beta} \times \mathrm{rot}_{\xi} \mathbf{q}^{(x)}) \cdot \mathbf{n}$ on the boundary $\partial\mathcal{Y}$ vanish by periodicity or because $\mathrm{rotq}^{(x)}$ vanishes for $\xi_x \to \pm\infty$. For $\beta = y$, we obtain $c_z^{(x)} = 0$, and for $\beta = z$, we obtain $c_y^{(x)} = 0$.

We now consider the problem on $\mathbf{q}^{(\alpha)}$ in (36) and we consider the integral

$$0 = \int_{\mathcal{Y}_m} (\xi_x \hat{y}) \cdot \mathrm{rot}_{\xi} \left( \frac{1}{\varepsilon_r} \mathrm{rot}_{\xi} \mathbf{q}^{(\alpha)} \right) = \hat{z} \cdot \int_{\mathcal{Y}_m} \frac{1}{\varepsilon_r} \mathrm{rot}_{\xi} \mathbf{q}^{(\alpha)} d\xi \\ - \xi_x^m \hat{x} \cdot \int_{\pm \xi_x^m} \frac{1}{\varepsilon_r} (\hat{y} \times \mathrm{rot}_{\xi} \mathbf{q}^{(\alpha)}) ds \qquad (62)$$

We have used that the contributions on $\partial\mathcal{Y}_m$ cancel by periodicity except at $\xi_x = \pm\xi_x^m$. We now consider the limit $\xi_x^m \to \infty$. From (43), we have

$$\lim_{\xi_x^m \to +\infty} \int_{\mathcal{Y}_m} \frac{1}{\varepsilon_r} \mathrm{rot}_{\xi} \mathbf{q}^{(\alpha)} d\xi = \mathbf{c}^{(\alpha)} + C_y \hat{\alpha} + \lim_{\xi_x^m \to +\infty} 2\xi_x^m \hat{\alpha}. \qquad (63)$$

Next as $\mathrm{rotq}^{(\alpha)}$ tends to $\hat{\alpha}$ as $\xi_x$ goes to $\pm\infty$, we also have

$$\lim_{\xi_x^m \to +\infty} \int_{\pm \xi_x^m} \frac{1}{\varepsilon_r} (\hat{y} \times \mathrm{rot}_{\xi} \mathbf{q}^{(\alpha)}) ds = (\hat{y} \times \hat{\alpha}). \qquad (64)$$

Using (63)-(64) in (62) in the limit $\xi_x^m \to \infty$, with $(\hat{y} \times \hat{\alpha}) \cdot \hat{x} = \delta_{\alpha z}$, the terms in $\xi_x^m$ cancel and we obtain that $c_z^{(\alpha)} = -C_y\delta_{\alpha z}$. Doing the same using $\mathbf{A} = (\xi_x \hat{z})$, we obtain $c_y^{(\alpha)} = -C_y\delta_{\alpha y}$.

**2) Relations between some coefficients:** We now use $\mathbf{A} = \mathbf{q}^{(j)}$ and proceed as above to get

$$\mathbf{R}_{ij} = c_i^{(j)} + (\hat{x} \times \hat{i}) \cdot \mathbf{b}^{(j)}, \qquad (65)$$

with $\mathbf{R}_{ij} = -\int_{\mathcal{Y}} \frac{1}{\varepsilon_r} \mathrm{rot}_{\xi} \mathbf{q}^{(i)} \cdot \mathrm{rot}_{\xi} \mathbf{q}^{(j)} d\xi$. As $\mathbf{R}_{ij} = \mathbf{R}_{ji}$, we deduce that

$$c_x^{(y)} = b_z^{(x)}, \quad c_x^{(z)} = -b_y^{(x)}, \quad b_z^{(z)} = -b_y^{(y)}. \qquad (66)$$

It follows that the parameters which remain to be determined reduce to 6 effective parameters $(b_y^{(x)}, b_z^{(x)}, b_y^{(y)}, b_z^{(y)}, b_y^{(z)})$ and $c_x^{(x)}$, according to (47)-(48).

## B. Elementary problems in terms of a scalar potential

**1) Scalar potentials:** The elementary fields $\mathbf{q}^{(i)}$ are defined in (35)-(36) as the solutions of two PDEs similar to electrostatic problems. Yet in practice, since these equations exhibits irrotational fields, one can define the following scalar potentials $p^{(x)}$ and $p^{(\alpha)}$, $\alpha = y, z$ on each constant $\varepsilon_r$ region (where the fields are smooth):

$$\mathrm{rot}_{\xi} \mathbf{q}^{(x)} + \hat{x} = \varepsilon_r \nabla_{\xi} p^{(x)}, \quad \mathrm{rot}_{\xi} \mathbf{q}^{(\alpha)} - \varepsilon_r \hat{\alpha} = \varepsilon_r \nabla_{\xi} p^{(\alpha)}.$$

Taking the divergence of the previous formulae we obtain:

$$\mathrm{div}_{\xi} (\varepsilon_r \nabla_{\xi} p^{(x)}) = 0, \quad \mathrm{div}_{\xi} (\varepsilon_r \nabla_{\xi} (p^{(\alpha)} + \xi_{\alpha})) = 0. \qquad (67)$$

We now move on to the boundary conditions of (35)-(36), resulting in the following limits for the scalar potentials

$$\lim_{\xi_x \to \pm\infty} \varepsilon_r^{\pm} \nabla_{\xi} p^{(x)} = \hat{x}, \quad \lim_{\xi_x \to \pm\infty} \varepsilon_r^{\pm} \nabla_{\xi} p^{(\alpha)} = \mathbf{0}. \qquad (68)$$

And in the same way for the continuity and periodicity:

$$\nabla_{\xi} p^{(i)} \times \mathbf{n} \text{ continuous and } p^{(i)}, \nabla_{\xi} p^{(i)} \xi_{\parallel} \text{ periodic.}$$

From the first condition we also infer that $p^{(i)}$ is continuous across the interfaces. Note that the scalar fields $p^{(i)}$ (as for the elementary fields $\mathbf{q}^{(i)}$) are well-defined up to a constant which could be fixed to an arbitrary value anywhere in the simulation domain.

**2) Susceptibilities in terms of the scalar potentials:** We can also simplify the $b_j^{(i)}$ and $c_j^{(i)}$ coefficients and thus the susceptibility calculation. The diagonal susceptibilities read

$$\chi_{cc}^{xx} = -(e^{\langle \varepsilon_r^{-1} \rangle} + dc_x^{(x)}), \\ \chi_{cc}^{yy} = e^{-\varepsilon_r^{-}} + e^{+\varepsilon_r^{+}} - db_z^{(y)}, \quad \chi_{cc}^{zz} = e^{-\varepsilon_r^{-}} + e^{+\varepsilon_r^{+}} + db_y^{(z)}, \qquad (69)$$

with

$$c_x^{(x)} = \int_{\mathcal{Y}} \left( \nabla p^{(x)} - \frac{\hat{x}}{\varepsilon_r} \right) \cdot \hat{x} \, d\xi, \\ b_z^{(y)} = - \int_{\mathcal{Y}} \varepsilon_r \nabla p^{(y)} \cdot \hat{y} \, d\xi + \int_{\mathcal{Y}_+} (\varepsilon_r^{-} - \varepsilon_r^{+}) \, d\xi + \int_{\mathcal{Y}_-} (\varepsilon_r^{+} - \varepsilon_r^{-}) \, d\xi, \\ b_y^{(z)} = \int_{\mathcal{Y}} \varepsilon_r \nabla p^{(z)} \cdot \hat{z} \, d\xi + \int_{\mathcal{Y}_+} (\varepsilon_r^{-} - \varepsilon_r^{+}) \, d\xi + \int_{\mathcal{Y}_-} (\varepsilon_r^{+} - \varepsilon_r^{-}) \, d\xi,$$

and the off-diagonal terms read

$$\chi_{cc}^{xy} = -db_z^{(x)}, \quad \chi_{cc}^{xz} = db_y^{(z)}, \quad \chi_{cc}^{yz} = db_y^{(y)}, \qquad (70)$$

with

$$b_z^{(x)} = - \int_{\mathcal{Y}} \varepsilon_r \nabla p^{(x)} \cdot \hat{y} \, d\xi, \quad b_y^{(x)} = \int_{\mathcal{Y}} \varepsilon_r \nabla p^{(x)} \cdot \hat{z} \, d\xi, \\ b_y^{(y)} = \int_{\mathcal{Y}} \varepsilon_r \nabla p^{(y)} \cdot \hat{z} \, d\xi. \qquad (71)$$

## C. Symmetries of the microstructures

A symmetry of the inclusion geometry with respect to the plane $(O, \hat{x}, \hat{y})$ ($z \to -z$) or with respect to the plane $(O, \hat{x}, \hat{z})$ ($y \to -y$) produces vanishing coefficients and vanishing susceptibilities due to symmetries of the elementary functions. Although it is possible to look at the symmetries of the $\mathbf{q}^{(i)}$,
$i = x, y, z$ satisfying (35)-(36), it is simpler to look at those of the scalar functions $p^{(i)}$ satisfying (67)-(68).

Let us consider a symmetry of the geometry with respect to the plane $(O, \hat{x}, \hat{y})$, that is the symmetry $\xi_z \to -\xi_z$. As the external loadings for the elementary problems (67)-(68) are along the principal directions $(\hat{x}$ for $p^{(x)}$, $\hat{\alpha}$ for $p^{(\alpha)}$), we deduce that $p^{(x)}$ and $p^{(y)}$ are even functions of $\xi_z$ (while $p^{(z)}$ is an odd function of $\xi_z$. From (71), it follows that $b_y^{(x)} = 0$ since $\partial_{\xi_z} p^{(x)}$ is an odd function of $\xi_z$ hence of vanishing integral (the same for $b_y^{(y)} = 0$).

Doing the same for a symmetry with respect to $O\hat{z}$ that is the symmetry $\xi_y \to -\xi_y$, we obtain that $p^{(x)}$ and $p^{(z)}$ are even functions of $\xi_y$ while $p^{(y)}$ is an odd function of $\xi_y$. It results that $b_{\hat{z}}^{(x)} = b_y^{(y)} = 0$. Using these results in (70), we obtain

$$ \text{Symmetry } z \to -z, \quad \chi_{ee}^{xz} = \chi_{ee}^{yz} = 0, $$

$$ \text{Symmetry } y \to -y, \quad \chi_{ee}^{xy} = \chi_{ee}^{yz} = 0. $$

Hence, inclusions with both symmetries $y \to -y$ and $z \to -z$ produce a diagonal susceptibility tensor.

## REFERENCES

[1] A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, "Planar photonics with metasurfaces," *Science*, vol. 339, no. 6125, 2013.

[2] P. Lalanne and P. Chavel, "Metalenses at visible wavelengths: past, present, perspectives," *Laser Photonics Rev.*, vol. 11, no. 3, p. 1600295, 2017.

[3] Y. Kivshar and A. Miroshnichenko, "Meta-optics with mie resonances," *Opt. Photonics News*, vol. 28, no. 1, pp. 24–31, 2017.

[4] A. Li, S. Singh, and D. Sievenpiper, "Metasurfaces and their applications," *Nanophotonics*, vol. 7, no. 6, pp. 989–1011, 2018.

[5] I. Kim, G. Yoon, J. Jang, P. Genevet, K. T. Nam, and J. Rho, "Outfitting next generation displays with optical metasurfaces," *ACS Photonics*, vol. 5, no. 10, pp. 3876–3895, 2018.

[6] C. L. Holloway and E. F. Kuester, "Electromagnetic metamaterials and metasurfaces: historical overview, characterization, and the effect of length scales," in *Dielectric Metamaterials*. Elsevier, 2020, pp. 1–38.

[7] N. M. Estakhri and A. Alù, "Recent progress in gradient metasurfaces," *JOSA B*, vol. 33, no. 2, pp. A21–A30, 2016.

[8] N. Schmitt, N. Georg, G. Brière, D. Loukrezis, S. Héron, S. Lanteri, C. Klitis, M. Sorel, U. Römer, H. De Gersem et al., "Optimization and uncertainty quantification of gradient index metasurfaces," *Optical Materials Express*, vol. 9, no. 2, pp. 892–910, 2019.

[9] H. Pinson and V. Ginis, "General framework for the frequency shifting of electromagnetic pulses using time-dependent surfaces," *Phys. Rev. B*, vol. 99, no. 20, p. 201407, 2019.

[10] N. Chamanara, Y. Vahabzadeh, and C. Caloz, "Simultaneous control of the spatial and temporal spectra of light with space-time varying metasurfaces," *IEEE Trans. Antennas Propag.*, vol. 67, no. 4, pp. 2430–2441, 2019.

[11] K. Achouri, G. D. Bernasconi, J. Butet, and O. J. Martin, "Homogenization and scattering analysis of second-harmonic generation in nonlinear metasurfaces," *IEEE Transactions on Antennas and Propagation*, vol. 66, no. 11, pp. 6061–6075, 2018.

[12] K. Achouri and O. J. Martin, "Fundamental properties and classification of polarization converting bianisotropic metasurfaces," *IEEE Transactions on Antennas and Propagation*, 2021.

[13] J. G. Maloney and G. S. Smith, "A comparison of methods for modeling electrically thin dielectric and conducting sheets in the finite-difference time-domain (fdtd) method," *IEEE Trans. Antennas Propag.*, vol. 41, no. 5, pp. 690–694, 1993.

[14] M. Idemen, "The maxwell's equations in the sense of distributions," *IEEE Trans. Antennas Propag.*, vol. 21, no. 5, pp. 736–738, 1973.

[15] M. Idemen and A. H. Serbest, "Boundary conditions of the electromagnetic field," *Electron. Lett.*, vol. 23, no. 13, pp. 704–705, 1987.

[16] B. Polat, "Approximate boundary relations on anisotropic sheets," *Progress In Electromag. Res.*, vol. 29, pp. 355–392, 2011.

[17] B. Delourme, H. Haddar, and P. Joly, "On the well-posedness, stability and accuracy of an asymptotic model for thin periodic interfaces in electromagnetic scattering problems," *Math. Models Methods Appl. Sci.*, vol. 23, no. 13, pp. 2433–2464, 2013.

[18] C. L. Holloway, E. F. Kuester, and A. Dienstfrey, "A homogenization technique for obtaining generalized sheet transition conditions for an arbitrarily shaped coated wire grating," *Radio Science*, vol. 49, no. 10, pp. 813–850, 2014.

[19] C. L. Holloway and E. F. Kuester, "A homogenization technique for obtaining generalized sheet-transition conditions for a metal film embedded in a magnetodielectric interface," *IEEE Trans. Antennas Propag.*, vol. 64, no. 11, pp. 4671–4686, 2016.

[20] M. Albooyeh, S. Tretyakov, and C. Simovski, "Electromagnetic characterization of bianisotropic metasurfaces on refractive substrates: General theoretical framework," *Annalen der Physik*, vol. 528, no. 9-10, pp. 721–737, 2016.

[21] S. Tretyakov, V. Asadchy, and A. Diaz-Rubio, "Metasurfaces for general control of reflection and transmission," in *Handbook of Metamaterials and Plasmonics*, vol. 29, no. 1. World Scientific, 2017, pp. 249–293.

[22] D. Morits and C. Simovski, "Electromagnetic characterization of planar and bulk metamaterials: A theoretical study," *Phys. Rev. B*, vol. 82, no. 16, p. 165114, 2010.

[23] T. D. Karamanos, A. I. Dimitriadis, and N. V. Kantartzis, "Robust technique for the polarisability matrix retrieval of bianisotropic scatterers via their reflection and transmission coefficients," *IET Microwaves, Antennas & Propagation*, vol. 8, no. 15, pp. 1398–1407, 2014.

[24] C. L. Holloway, E. F. Kuester, and A. H. Haddab, "Using reflection and transmission coefficients to retrieve surface parameters for an anisotropic metascreen: With a discussion on conversion between te and tm polarizations," *J. Appl. Phys.*, vol. 125, no. 9, p. 095102, 2019.

[25] M. Goodarzi and T. Pakizeh, "Retrieving effective surface susceptibilities of high-index metasurfaces based on dipole approximation," *Optics Comm.*, vol. 483, p. 126659, 2020.

[26] C. L. Holloway, E. F. Kuester, J. A. Gordon, J. O'Hara, J. Booth, and D. R. Smith, "An overview of the theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials," *IEEE Antennas Propag. Mag.*, vol. 54, no. 2, pp. 10–35, 2012.

[27] M. Dehmollaian, Y. Vahabzadeh, K. Achouri, and C. Caloz, "Metasurface modeling by a thin slab," *arXiv preprint arXiv:2005.00901*, 2020.

[28] J.-J. Marigo and A. Maurel, "Homogenization models for thin rigid structured surfaces and films," *J. Acoust. Soc. Am.*, vol. 140, no. 1, pp. 260–273, 2016.

[29] M. David, J.-J. Marigo, and C. Pideri, "Homogenized interface model describing inhomogeneities located on a surface," *J. Elas.*, vol. 109, no. 2, pp. 153–187, 2012.

[30] K. Pham, A. Maurel, and J.-J. Marigo, "Revisiting imperfect interface laws for two-dimensional elastodynamics," *Proc. R. Soc. A*, p. to appear, 2021.

[31] B. Delourme, H. Haddar, and P. Joly, "Approximate models for wave propagation across thin periodic interfaces," *J. Math. Pures Appl.*, vol. 98, no. 1, pp. 28–71, 2012.

[32] A. Maurel, K. Pham, and J.-J. Marigo, "Homogenization of thin 3d periodic structures in the time domain–effective boundary and jump conditions," *Fundamentals and Applications of Acoustic Metamaterials: From Seismic to Radio Frequency*, vol. 1, pp. 73–105, 2019.

[33] L.-K. Wu, L.-T. Han et al., "Implementation and application of resistive sheet boundary-condition in the finite-difference time-domain method," *IEEE Trans. Antennas Propag.*, vol. 40, no. 6, pp. 628–633, 1992.

[34] T. J. Smy, S. A. Stewart, J. G. Rahmeier, and S. Gupta, "Fdtd simulation of dispersive metasurfaces with lorentzian surface susceptibilities," *IEEE Access*, vol. 8, pp. 83027–83040, 2020.

[35] B. Delourme, E. Lunéville, J.-J. Marigo, A. Maurel, J.-F. Mercier, and K. Pham, "A stable, unified, model for resonant faraday cages," *Proc. R. Soc. A*, no. to appear, 2021.

[36] M. Touboul, K. Pham, A. Maurel, J.-J. Marigo, B. Lombard, and C. Bellis, "Effective resonant model and simulations in the time-domain of wave scattering from a periodic row of highly-contrasted inclusions," *J. Elas.*, 2020.

[37] K. Pham, J.-F. Mercier, D. Fuster, J.-J. Marigo, and A. Maurel, "Scattering of acoustic waves by a nonlinear resonant bubbly screen," *J. Fluid Mech.*, vol. 906, 2021.

[38] P. A. Tirkas and K. R. Demarest, "Modeling of thin dielectric structures using the finite-difference time-domain technique," *IEEE Trans. Antennas Propag.*, vol. 39, no. 9, pp. 1338–1344, 1991.

[39] M. E. Gurtin, "Variational principles for linear elastodynamics," Brown Univ. Providence, div. Appl. Math., Tech. Rep., 1963.
# Supplementary material: Stable GSTC formulation for Maxwell's equations

Nicolas Lebbe*, Kim Pham†, and Agnès Maurel‡, * Université Côte d'Azur, Inria, CNRS, LJAD, 06902 Sophia Antipolis Cedex †IMSIA, CNRS, EDF, CEA, ENSTA Paris, Institut Polytechnique de Paris, 828 Bd des Maréchaux, 91732 Palaiseau, France, ‡Institut Langevin, ESPCI Paris, PSL University, CNRS, 1 rue Jussieu, 75005 Paris, France,

**Abstract**—In this supplementary material we provide additional informations on the numerical implementation of the GSTCs by means of the associated variational formulation. We also provide informations on the elementary problems solved for the different geometries, in two- and three-dimensions, reported in the main document.

**Index Terms**—numerical implementation of the enlarged GSTCs.

## I. VARIATIONAL FORMULATION OF THE EFFECTIVE PROBLEM

We will also assume that the microstructure is such that the electric susceptibility is of the following form:

$$ \overline{\overline{\chi}}_{ee} = \begin{pmatrix} \chi_{ee}^{xx} & 0 & 0 \\ 0 & \chi_{ee}^{yy} & \chi_{ee}^{yz} \\ 0 & \chi_{ee}^{zy} & \chi_{ee}^{zz} \end{pmatrix}. \quad (1) $$

The validity of the above form is discussed below. We then provide the derivation of the variational formulation of the effective problem, which has been used in the implementation of the effective problem.

### A. Symmetry $x \to -x$ for $\epsilon_r^+ = \epsilon_r^-$

When the array is placed in identical substrates $\epsilon_r = \epsilon_r^+ = \epsilon_r^-$ with particles satisfying the symmetry $x \to -x$, we have that

$$ \chi_{ee}^{xy} = \chi_{ee}^{xz} = 0. \quad (2) $$

Indeed, from (16), the external loading on $p^{(x)}$ is along $\hat{x}$, hence $p^{(x)} \sim \frac{\xi_x}{\epsilon_r}$ as $\xi_x \to \pm\infty$ is an odd function of $\xi_x$ (everywhere). It results that $b_y^{(x)} = b_y^{(x)} = 0$ in (18), hence $\chi_{ee}^{xy} = \chi_{ee}^{xz} = 0$ from (19).

It is worth noticing that when the susceptibility tensor is plain, that is (i) in the absence of such symmetry $x \to -x$ for a single substrate or (ii) in the absence of the two symmetries $y \to -y$ and $z \to -z$ for a single or two substrates (these two symmetries provide a diagonal susceptibility tensor), the variational formulation presented in I has to be adapted. To the best of our knowledge, this case has not been considered in the literature previously and its implementation remains an open problem.

Corresponding author: N. Lebbe (email: lebbe.nicolas@gmail.com).

### B. Derivation of the variational formulation

We assume a harmonic time-dependence $e^{i\omega t}$, hence the formulation of the Maxwell equations used in the FEM reads

$$ \mathrm{rot} \, \mathrm{rot} \, \mathbf{E} - k_0^2 \varepsilon_r \mathbf{E} = 0, \quad \mathbf{H} = -\frac{1}{i\omega\mu_0} \mathrm{rot} \, \mathbf{E}, \quad (3) $$

with $k_0 = \omega\sqrt{\varepsilon_0\mu_0} = \frac{2\pi}{\lambda}$ and $\lambda$ the wavelength (note that considering the $\mathbf{E}$ field instead of $\mathbf{H}$ is not mandatory but will ease the calculation afterward). To obtain the variational formulation associated to the GSTCs we will consider independently the top $\mathcal{D}^- = (-\infty, e^-) \times \mathbb{R}^2$ and bottom $\mathcal{D}^+ = (e^+, \infty) \times \mathbb{R}^2$ domains. After multiplication by a test function $\phi$ of the first equation in (3) and integration, we obtain

$$ \int_{\mathcal{D}^+ \cup \mathcal{D}^-} (\mathbf{rot} \, \mathbf{E} \cdot \mathbf{rot} \, \phi^* - k_0^2 \varepsilon_r \mathbf{E} \cdot \phi^*) \, dv \\ + \int_{\Gamma^-} \hat{\mathbf{x}} \times \mathbf{rot} \, \mathbf{E} \cdot \phi_{\parallel}^* \, ds - \int_{\Gamma^+} \hat{\mathbf{x}} \times \mathbf{rot} \, \mathbf{E} \cdot \phi_{\parallel}^* \, ds = 0. \quad (4) $$

As we have done to derive the energy, the last two integrals are expressed in terms of the jump and average operators (4) (main text) resulting in

$$ \begin{aligned} & - \int_{\Gamma} \left[ \left\| \hat{\mathbf{x}} \times \mathbf{rot} \, \mathbf{E} \cdot \phi_{\parallel}^{*} \right\| \right] \, ds \\ & = \int_{\Gamma} \left( \left\| \mathbf{rot} \, \mathbf{E} \right\| \times \hat{\mathbf{x}} \cdot \phi_{\text{av}, \parallel}^{*} + \mathbf{rot} \, \mathbf{E}_{\text{av}} \times \hat{\mathbf{x}} \cdot \left\| \phi \right\|_{\parallel}^{*} \right) \, ds. \end{aligned} \quad (5) $$

Now, we express the GSTCs (5) (main text) in term of the $\mathbf{E}$ field only (and in the harmonic regime), namely

$$ \begin{aligned} \llbracket \mathbf{rot} \, \mathbf{E} \rrbracket \times \hat{\mathbf{x}} &= -k_0^2 (\overline{\overline{\chi}}_{ee} \mathbf{E}_{av})_{\parallel} + e \nabla_{\parallel} (\mathbf{rot} \, \mathbf{E}_{av})_{x} \times \hat{\mathbf{x}}, \\ \llbracket \mathbf{E} \rrbracket \times \hat{\mathbf{x}} &= -e (\mathbf{rot} \, \mathbf{E}_{av})_{\parallel} - \nabla_{\parallel} (\overline{\overline{\chi}}_{ee} \mathbf{E}_{av})_{x} \times \hat{\mathbf{x}}, \end{aligned} $$

We shall now work on the above forms of the jumps to deduce a suitable form of the variational formulation. To do so, we shall express (i) $\llbracket \mathbf{rot} \mathbf{E} \rrbracket \times \hat{\mathbf{x}}$ as a function of $\mathbf{E}_{av}$ and (ii) $\llbracket \mathbf{E} \rrbracket \times \hat{\mathbf{x}}$ as a function of

$$ \mathbf{U} = (\mathbf{rot} \, \mathbf{E}_{av})_{\parallel}. \quad (6) $$

* For (i), we use that, for an electric susceptibility tensor of the form (1), we have $(\overline{\overline{\chi}}_{ee} \mathbf{E}_{av})_x = \frac{\chi_{ee}^{xx}}{\varepsilon_0} D_{x,av} = \chi_{ee}^{xx} \varepsilon_r E_{x,av}$ (see (8), main text), hence

$$ \llbracket \mathbf{rot} \, \mathbf{E} \rrbracket \times \hat{\mathbf{x}} = -k_0^2 (\overline{\overline{\chi}}_{ee} \mathbf{E}_{av})_{\parallel} + e \nabla_{\parallel} (\mathbf{rot} \, \mathbf{E}_{av})_{x} \times \hat{\mathbf{x}}. \quad (7) $$
Using the above form in the first integral of (5) and integrating by part <fnref n="1" /> we then find that

$$ \int_{\Gamma} [[\mathbf{rot} \mathbf{E}] \times \hat{\mathbf{x}} \cdot \mathbf{\phi}_{\text{av},\parallel}^{*}] \, \mathrm{d}s = -k_{0}^{2} \int_{\Gamma} [\overline{(\bar{\chi}_{\mathrm{ec}} \mathbf{E}_{\text{av}})}_{\parallel} \cdot \mathbf{\phi}_{\text{av},\parallel}^{*}] \, \mathrm{d}s \\ +e \int_{\Gamma} \nabla_{\parallel} \times \mathbf{E}_{\text{av},\parallel} \cdot \nabla_{\parallel} \times \mathbf{\phi}_{\text{av},\parallel}^{*} \, \mathrm{d}s. \qquad (8) $$

• For (ii), we use that, by definition (6), the second integral of (5) reads

$$ \int_{\Gamma} \mathbf{rot} \mathbf{E}_{\text{av}} \times \hat{\mathbf{x}} \cdot [[\mathbf{\phi}]_{\parallel}^{*}] \, \mathrm{d}s = \int_{\Gamma} \mathbf{U} \times \hat{\mathbf{x}} \cdot [[\mathbf{\phi}]_{\parallel}^{*}] \, \mathrm{d}s, \quad (9) $$

and it remains us to establish a variational formulation associated with **U**. To do so, we come back to $[[\mathbf{E}] \times \hat{\mathbf{x}}]$ and remark that with the form of the susceptibility tensor (1) we have $(\overline{\chi}_{\mathrm{ec}} \tilde{\mathbf{E}})_{\mathrm{x}} = \chi_{\mathrm{ec}}^{xx} \frac{1}{\varepsilon_0} D_{x,\mathrm{av}} = \chi_{\mathrm{ec}}^{xx} \varepsilon_t E_{x,\mathrm{av}}$ (see (8), main text). Next, and following [1], we express $\mathbf{E}_{\mathrm{av}}$ as a function of **rot** $\mathbf{rot} \mathbf{E}_{\mathrm{av}}$ using (3) to get that **U** is the (periodic) solution of the partial differential equation on $\Gamma$

$$ \frac{\chi_{\mathrm{ec}}^{xx}}{k_0^2} \nabla_{\parallel} \times (\nabla_{\parallel} \times \mathbf{U}) + e\mathbf{U} = \hat{\mathbf{x}} \times [[\mathbf{E}]], $$

where we have used that $\nabla_{\parallel}(\mathbf{rot} \mathbf{rot} \mathbf{E}_{\mathrm{av}})_x \times \hat{\mathbf{x}} = \nabla_{\parallel} \times (\nabla_{\parallel} \mathbf{U})$. The associated variational formulation reads

$$ \int_{\Gamma} \left( \frac{\chi_{\mathrm{ec}}^{xx}}{k_0^2} \nabla_{\parallel} \times \mathbf{U} \cdot \nabla_{\parallel} \times \psi^* + (e\mathbf{U} - \hat{\mathbf{x}} \times [[\mathbf{E}]]) \cdot \psi^* \right) \, \mathrm{d}s = 0. \quad (10) $$

Using (8), (9) and (10) (5) in (4), we obtain the final variational formulation set on $\mathbf{E}_{\mathrm{av}}$ and $\mathbf{U}$ given with test functions $\phi$ and $\psi$ as

$$ \begin{aligned} 0 = & \int_{\mathcal{D}^+ \cup \mathcal{D}^-} (\mathbf{rot} \mathbf{E} \cdot \mathbf{rot} \phi^* - k_0^2 \varepsilon_t \mathbf{E} \cdot \phi^*) \, \mathrm{d}v \\ & -k_0^2 \int_{\Gamma} (\overline{\chi}_{\mathrm{ec}} \mathbf{E}_{\mathrm{av}})_{\parallel} \cdot \phi_{\mathrm{av},\parallel}^* \, \mathrm{d}s + e \int_{\Gamma} \nabla_{\parallel} \times \mathbf{E}_{\mathrm{av},\parallel} \cdot \nabla_{\parallel} \times \phi_{\mathrm{av},\parallel}^* \, \mathrm{d}s \\ & + \int_{\Gamma} \mathbf{U} \times \hat{\mathbf{x}} \cdot [[\mathbf{\phi}]_{\parallel}^*] \, \mathrm{d}s \\ & + \frac{\chi_{\mathrm{ec}}^{xx}}{k_0^2} \int_{\Gamma} \nabla_{\parallel} \times \mathbf{U} \cdot \nabla_{\parallel} \times \psi^* \, \mathrm{d}s + \int_{\Gamma} (e\mathbf{U} - \hat{\mathbf{x}} \times [[\mathbf{E}]]) \cdot \psi^* \, \mathrm{d}s. \end{aligned} \quad (11) $$

## II. ELEMENTARY PROBLEMS

As previously said, the implementation of the GSTCs and of the associated elementary problems were done using Comsol Multiphysics freely available at https://www.comsol.fr/community/exchange/841/. We detail them below in two- and three-dimensions.

## A. Elementary problem in two-dimensions

In two dimensions, we have $\frac{\partial}{\partial z} = 0$ and the elementary functions (35)-(36) (main text) $\mathbf{q}^{(i)}(\xi_x, \xi_y)$, $i = x, y, z$, no not depend on $z$. Deep simplifications occur and the elementary problems are reduced to the determination of the two scalar functions $q_z^{(x)}$ and $q_z^{(y)}$. The function $q_z^{(x)}$ is solution to

$$ \left\{ \begin{array}{l} \mathrm{div} \, \boldsymbol{\xi} \left( \frac{1}{\varepsilon_t} \nabla \boldsymbol{\xi} \left( q_z^{(x)} + y \right) \right) = 0, \\ \lim_{\xi_x \to \pm\infty} \nabla \boldsymbol{\xi} q_z^{(x)} = 0, \end{array} \right. $$

and, defining $\mathbf{f} = (\frac{1}{\varepsilon_t} \nabla \boldsymbol{\xi} (q_z^{(x)} + y)), (q_z^{(x)}, \mathbf{f} \cdot \mathbf{n})$ continuous and $(q_z^{(y)}, \mathbf{f})$ periodic.

The function $q_z^{(y)}$ is solution to

$$ \left\{ \begin{array}{l} \mathrm{div} \, \boldsymbol{\xi} \left( \frac{1}{\varepsilon_t} \nabla \boldsymbol{\xi} q_z^{(y)} \right) = 0, \\ \lim_{\xi_x \to \pm\infty} \frac{1}{\varepsilon_t} \nabla \boldsymbol{\xi} q_z^{(y)} = -\hat{\mathbf{x}}, \end{array} \right. $$

and, defining $\mathbf{g} = (\frac{1}{\varepsilon_t} \nabla \boldsymbol{\xi} q_z^{(y)}), (q_z^{(y)}, \mathbf{g} \cdot \mathbf{n})$ continuous and $(q_z^{(y)}, \mathbf{g})$ periodic. We recover the two elementary problems involved in [2] in which polarized TM waves in two-dimensions are considered. Next, it is straightforward to obtain that

$$ q_x^{(x)} = q_y^{(x)} = q_x^{(y)} = q_y^{(y)} = q_z^{(z)} = 0, \quad (12) $$

Among the effective parameters (37) and (43) (main text), those involved in the terms of the susceptibility tensor, (47)-(48) (main text), are given by

$$ \begin{aligned} \chi_{\mathrm{ec}}^{xx} &= - (e \langle \varepsilon_t^{-1} \rangle + d c_x^{(x)}), \\ \chi_{\mathrm{ec}}^{yy} &= (e^{-\varepsilon_t} + e^+ \varepsilon_t^+ - d b_z^{(y)}), \quad \chi_{\mathrm{ec}}^{zz} = e \langle \varepsilon_t \rangle. \end{aligned} \quad (13) $$

(see <fnref n="2" />), and vanishing off-diagonal terms

$$ \chi_{\mathrm{ec}}^{xy} = -db_z^{(x)}, \quad \chi_{\mathrm{ec}}^{xz} = db_y^{(x)} = 0, \quad \chi_{\mathrm{ec}}^{yz} = db_y^{(y)} = 0, \quad (14) $$

and we have considered configurations for which $\chi_{\mathrm{ec}}^{xy} = 0$

<figure><img src="image_13.png" /><figcaption>Fig. 1: (a) Geometry of the elementary cell with a clover-like shape particle in two-dimensions. Elementary problems (b) $q_z^{(x)}(\xi)$ (left) and (c) $q_z^{(y)}(\xi)$ (right), $\xi = (\xi_x, \xi_y)$.</figcaption></figure>

<footnote n="1">¹denoting **a** = (rot**E**av) hence $a_x = (\partial_y E_z - \partial_z E_y)$, we have $\int_{\Gamma} \nabla_{\parallel} a_x \times \hat{\mathbf{x}} \cdot \mathbf{\phi}_{\parallel} \, \mathrm{d}s = \int_{\Gamma} (\partial_z a_x \phi_y - \partial_y a_x \phi_z) \, \mathrm{d}s = \int_{\Gamma} a_x (\partial_y \phi_z - \partial_z \phi_y) \, \mathrm{d}s$ up to local boundaries terms (at the extremities of $\Gamma$) which cancel if periodic boundary conditions are assumed in the computational domain.</footnote>

<footnote n="2">²Note that $(q_x^{(z)}, q_y^{(z)})$ are not vanishing, but it is straightforward to show that $p^{(z)}$ satisfying (67) (main text) is constant, which along with (18), produces the announced value of $\chi_{\mathrm{ec}}^{zz}$.</footnote>
(for instance, with symmetry $y \to -y$ resulting in $q_z^{(x)}$ odd with respect to $y$).

$$ 
\begin{aligned}
c_x^{(x)} &= \int_y \frac{1}{\epsilon_t} \nabla q_z^{(x)} \, d\xi, \\
b_z^{(y)} &= \lim_{\xi_x \to +\infty} (q_z^{(y)} + \epsilon_t^+ \xi_x) - \lim_{\xi_x \to -\infty} (q_z^{(y)} + \epsilon_t^- \xi_x), \\
b_z^{(x)} &= \lim_{\xi_x \to +\infty} q_z^{(x)} - \lim_{\xi_x \to -\infty} q_z^{(x)}.
\end{aligned} 
\quad (15) 
$$

For the clover-shaped silicon particle, the effective parameters obtained numerically are collected in the Table I and we report the shapes of the fields $q_z^{(x)}$ and $q_z^{(y)}$ in figure 1.

<table><caption>TABLE I: Numerical computation of the non zero effective coefficients for the clover-shaped silicon inclusion between silica and air in two-dimensions.</caption><thead><tr><th></th><th>c<sub>x</sub><sup>(x)</sup></th><th>b<sub>z</sub><sup>(y)</sup></th><th>b<sub>z</sub><sup>(x)</sup></th></tr></thead><tbody><tr><td>clover</td><td>-0.0288</td><td>1.57</td><td>1.04E-7</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Fig. 2: (a) Geometry of the elementary cell for the cylindrical-shaped microstructure. (b) Elementary solutions $p^{(i)}(\xi)$, $i = x, y, z$.</figcaption></figure>

## B. Elementary problems in three-dimensions

In three-dimensions, as detailed in the appendix B of the main document, we have solved in practice elementary problems set on scalar potentials $p^{(i)}(\xi)$, $i = x, y, z$ from which the susceptibility tensor is deduced. For the sake of completeness, the procedure is recalled below. The scalar potential $p^{(x)}$ is solution to the elementary problem

$$ 
\left\{ 
\begin{array}{l}
\operatorname{div}\xi (\epsilon_t \nabla\xi p^{(x)}) = 0, \\
\lim_{\xi_x \to \pm\infty} \epsilon_t^{\pm} \nabla\xi p^{(x)} = \hat{x}, \\
p^{(x)}, \nabla\xi p^{(x)} \times \mathbf{n} \text{ continuous, } p^{(x)}, \nabla\xi p^{(x)} \xi_{\parallel}\text{-periodic.}
\end{array}
\right. 
\quad (16) 
$$

and $p^{(\alpha)}$, $\alpha = y, z$, are solutions to

$$ 
\left\{ 
\begin{array}{l}
\operatorname{div}\xi (\epsilon_t \nabla\xi (p^{(\alpha)} + \xi_{\alpha})) = 0, \\
\lim_{\xi_x \to \pm\infty} \epsilon_t^{\pm} \nabla\xi p^{(\alpha)} = \mathbf{0}, \\
p^{(\alpha)}, \nabla\xi p^{(\alpha)} \times \mathbf{n} \text{ continuous } p^{(\alpha)}, \nabla\xi p^{(\alpha)} \xi_{\parallel}\text{-periodic.}
\end{array}
\right. 
\quad (17) 
$$

Note that the scalar fields $p^{(i)}$ are well-defined up to a constant which could be fixed to an arbitrary value anywhere in the simulation domain. The effective parameters entering

<figure><img src="image_15.png" /><figcaption>Fig. 3: Case of the rotated cuboid microstructure, same representation as in 2.</figcaption></figure>
in the GSTCs are obtained from the solutions $p^{(i)}$ computed numerically and they are defined by

$$
\begin{align}
c_x^{(x)} &= \int_y \left( \nabla p^{(x)} - \frac{\hat{x}}{\epsilon_r} \right) \cdot \hat{x} \, d\xi, \nonumber \\
b_z^{(y)} &= -\int_y \epsilon_r \nabla p^{(y)} \cdot \hat{y} \, d\xi + e_\xi^\tau \varphi^\tau (\epsilon_r^\tau - \epsilon_r^{sc}) + e_\xi^+ \varphi^+ (\epsilon_r^+ - \epsilon_r^{sc}), \nonumber \\
b_y^{(z)} &= \int_y \epsilon_r \nabla p^{(z)} \cdot \hat{z} \, d\xi - e_\xi^\tau \varphi^\tau (\epsilon_r^\tau - \epsilon_r^{sc}) - e_\xi^+ \varphi^+ (\epsilon_r^+ - \epsilon_r^{sc}), \nonumber \\
b_z^{(x)} &= -\int_y \epsilon_r \nabla p^{(x)} \cdot \hat{y} \, d\xi, \nonumber \\
b_y^{(x)} &= \int_y \epsilon_r \nabla p^{(x)} \cdot \hat{z} \, d\xi, \nonumber \\
b_y^{(y)} &= \int_y \epsilon_r \nabla p^{(y)} \cdot \hat{z} \, d\xi. \tag{18}
\end{align}
$$

From the above coefficients, we eventually obtain the sus-
ceptibility tensor with

$$
\begin{align}
\chi_{ee}^{xx} &= -(e\langle\epsilon_r^{-1}\rangle + dc_x^{(x)}), \nonumber \\
\chi_{ee}^{yy} &= e^\tau\epsilon_r^\tau + e^+\epsilon_r^+ - db_z^{(y)}, \quad \chi_{ee}^{zz} = e^\tau\epsilon_r^\tau + e^+\epsilon_r^+ + db_y^{(z)}, \nonumber \\
\chi_{ee}^{xy} &= -db_z^{(x)}, \quad \chi_{ee}^{xz} = db_y^{(x)}, \quad \chi_{ee}^{yz} = db_y^{(y)}, \tag{19}
\end{align}
$$

We report the same information for the three-dimensional
cases where we have considered silicon particles (cylinders
and cuboids rotated of 30°).

<table><caption>TABLE II: Numerical computation of the effective coefficients for the cylindrical and cuboid particle in the air in three dimensions.</caption>
  <thead>
    <tr>
      <th></th>
      <th>c<sub>x</sub><sup>(x)</sup></th>
      <th>b<sub>z</sub><sup>(y)</sup></th>
      <th>b<sub>y</sub><sup>(z)</sup></th>
      <th>b<sub>z</sub><sup>(x)</sup></th>
      <th>b<sub>y</sub><sup>(x)</sup></th>
      <th>b<sub>y</sub><sup>(y)</sup></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>cylinder</th>
      <td>-0.153</td>
      <td>-0.417</td>
      <td>0.417</td>
      <td>1.65E-5</td>
      <td>1.05E-5</td>
      <td>-7.83E-7</td>
    </tr>
    <tr>
      <th>rotated<br>cuboid</th>
      <td>-0.158</td>
      <td>-0.956</td>
      <td>0.648</td>
      <td>-2.13E-5</td>
      <td>1.29E-4</td>
      <td>0.113</td>
    </tr>
  </tbody>
</table>

REFERENCES

[1] B. Delourme, H. Haddar, and P. Joly, “On the well-posedness, stability and accuracy of an asymptotic model for thin periodic interfaces in electromagnetic scattering problems,” *Mathematical Models and Methods in Applied Sciences*, vol. 23, no. 13, pp. 2433–2464, 2013.

[2] B. Gallas, A. Maurel, J.-J. Marigo, and A. Ourir, “Light scattering by periodic rough surfaces: equivalent jump conditions,” *JOSA A*, vol. 34, no. 12, pp. 2181–2188, 2017.