# High-frequency homogenization for periodic dispersive media

Marie Touboul Department of Mathematics, Imperial College London, London SW7
2AZ, UK<EMAIL_ADDRESS>). Benjamin Vial Department of Mathematics,
Imperial College London, London SW7 2AZ, UK<EMAIL_ADDRESS>). Raphaël
Assier Department of Mathematics, University of Manchester, Manchester, UK
<EMAIL_ADDRESS>). Sébastien Guenneau The Blackett
Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, UK
<EMAIL_ADDRESS>). Richard V. Craster Department of Mathematics,
UMI 2004 Abraham de Moivre-CNRS, Department of Mechanical Engineering,
Imperial College London, London SW7 2AZ, UK<EMAIL_ADDRESS>).

###### Abstract

High-frequency homogenization is used to study dispersive media, containing
inclusions placed periodically, for which the properties of the material
depend on the frequency (Lorentz or Drude model with damping, for example).
Effective properties are obtained near a given point of the dispersion diagram
in frequency-wavenumber space. The asymptotic approximations of the dispersion
diagrams, and the wavefields, so obtained are then cross-validated via
detailed comparison with finite element method simulations in both one and two
dimensions.

##  1 Introduction

Dispersive media have material properties that are frequency dependent and
form an important and commonly occurring class of materials in physics,
material science and engineering, they are of increasing importance in the
design of advanced materials particularly in plasmonics [1, 2]. In addition to
practical relevance, dispersive media have interesting mathematical features
for instance arising from resonances that lead to unusual material properties
(e.g. artificial magnetism and negative refraction in metamaterials [3], zero-
index band gaps [4], or folded bands [5]).
In electromagnetism, e.g. in optics, the electrodynamic properties of metals
differ from those of dielectrics and, in particular, the relative permittivity
is frequency dependent. The underlying physics is that free electrons form a
plasma, are free to move, support a current, and create forces. At low
frequencies in, say, microwaves, the effect is minimal, but as one approaches
the optical regime, at TeraHertz (THz) frequencies, the plasma frequency
becomes commensurate with the optical frequencies and the frequency dependence
of the properties are important. Indeed the field has a long and distinguished
history, the colours of metal glasses were investigated using the Drude model
by Maxwell-Garnet [6], and the reflectance properties of metals and entire
areas of physics such as plasmonics rely upon the properties captured within
dispersive media. Historically, Drude, Lorentz and Sommerfeld pioneered the
material models that are now commonly used [7], and used in our numerical
examples in section 2.5. There is also broader usage of the Drude-Lorenz-
Sommerfeld model in other fields as they are also well suited to describe
resonant media, arising due to high material contrast of inclusions either in
elasticity [8, 9] or in electromagnetism [10, 11], or due to the presence of
geometric resonances in acoustics e.g. Helmholtz resonators [12, 13].

Structuring dispersive media to contain periodically arranged inclusions, or
alternating layers [4] creates additional difficulties due to the introduction
of geometric structure and the resulting scattering and reflection from
inclusions and surfaces; there are applications for photonic crystals studied
with Drude-like models [14], creating complete electromagnetic band gaps, or
even for the more numerically challenging metallo-dielectric photonic crystals
[15]. Critical to the study of wave propagation through periodic media is the
concept of Bloch waves and of dispersion curves that connect frequency to the
phase shift across a single cell of the periodic medium [16, 17]; the
resulting band diagrams then neatly encapsulate the essential physics of the
structure, i.e. band-gaps of forbidden frequencies can arise where waves will
not propagate in the structure. Our aim here is to provide homogenization
models that allow us to side-step numerical issues, generate insight, and
create effective media for wave propagation through dispersive media that
contain periodically arranged inclusions, or consist of repeating layers.

Classically, dynamic homogenization is understood as a low-frequency
approximation to wave propagation in heterogeneous media. A particularly
successful method, for periodic media, is the two-scale asymptotic expansion
method and the notion of slow and fast variables [18, 19, 20, 21]. In the case
of dispersive media, it has been extended to surfaces recently [22]. The idea
of high-frequency homogenization (HFH), introduced in [23], is to use similar
asymptotic methods to approximate how the dispersion relation, and hence the
media behave, near a given point in wavenumber-frequency space that satisfies
the dispersion relation; recent work by [24] pushes this asymptotic analysis
one order further. A uniform approximation is also derived, considering the
fact that some branches do not intersect at the edges of the Brillouin zone
but are close to each other. Other works concerned the inclusion of a source
term [25] and the derivation of the process in the time domain [26]. The
methodology has also been applied to several other configurations such as
discrete lattice media [27, 28], frame structures [29], optics [30], elastic
plates [31], full vector wave systems [32], elastic composites [33],
reticulated structures [34], or imperfect interfaces [35, 36]. The two-scale
approach is also connected to homogenization near a neighbourhood of an edge
gap in the context of approximations of operator resolvents [37] and there are
connections into spectral theory.
Here we extend the HFH method to the case of dispersive media where the
properties of the material depend on the frequency; this is not a routine
extension as the dispersion curves are now complex and additional
complications due to the frequency dependence, including resonances, now
occur. In Section 2 we consider a one-dimensional (1D) setting of waves though
a laminate of alternating layers, and high-frequency homogenization is applied
for different cases: single eigenvalues at the edges, double eigenvalues at
the edges, nearby eigenvalues at the edges, or single eigenvalues outside the
edges (when no damping is considered); numerical examples are then presented
to cross-validate the asymptotic approximations developed. In Section 3, the
asymptotic results are extended to two-dimensions (2D) and then cross-
validated via comparisons with finite element simulations for metallic rods in
a vacuum. The effective parameter obtained with HFH is also used to
investigate properties of the material.

##  2 One-dimensional (1D) case

###  2.1 Setting

We begin with the 1D case and consider linear waves propagating at a given
angular frequency $\omega$ through a dispersive periodic medium of periodicity
$h(>0)$ and with a macroscopic characteristic length $L(>0)$; physically this
would correspond to a laminate of dispersive medium layers where the layers
alternate with different material properties or a bimaterial string
constructed from alternating dispersive media, see Fig. 1 and for clarity of
exposition and notation we will fix one of the media to be non-dispersive. We
denote the physical space variable as $X$; the material parameters
$\hat{a}_{h}(X,\omega)$ and $\hat{b}_{h}(X,\omega)$ are assumed to be
$h$-periodic in $X$ and frequency dependent. The governing equation for the
field $U_{h}$ is:

$\frac{d}{dX}\left(\hat{a}_{h}(X,\omega)\frac{dU_{h}}{dX})\right)+\hat{b}_{h}(X,\omega)\omega^{2}U_{h}=0.$
(1)

This equation is very general and the field $U_{h}$ describes the transverse
electric field for s-polarisation, the transverse magnetic field for
p-polarisation in electromagnetism, the displacement in elasticity, or the
pressure in acoustics; the parameter $\hat{a}_{h}$ is then the inverse of the
permeability, the inverse of the permittivity, the shear modulus, or the
inverse of mass density, respectively, while $\hat{b}_{h}$ denotes the
permittivity, the permeability, the mass density, or the compressibility,
respectively. Henceforth we will assume the elastic setting in terms of
notation.
The unit cell is divided into two parts distinguished by a “volume fraction”
$\phi$. The left part is characterized by constant positive physical
parameters, while in the right part they are frequency dependent and
dispersive:

$\hat{a}_{h}(X,\omega)=\left\\{\begin{aligned} &{a}_{0}&\text{ for
}X\in(0,\phi h)\\\ &\hat{\mathbb{a}}(\omega)&\text{ for }X\in(\phi
h,h)\end{aligned}\right.$ (2) $\hat{b}_{h}(X,\omega)=\left\\{\begin{aligned}
&{b}_{0}&\text{ for }X\in(0,\phi h)\\\ &\hat{\mathbb{b}}(\omega)&\text{ for
}X\in(\phi h,h).\end{aligned}\right.$ (3)

For the frequency dependence the example of a Lorentz type dependence is given
in Appendix A and used in the numerical examples in section 2.5 .
The edges of the periodic cell are assumed, without loss of generality, to be
located at $X_{n}=nh$ for $n\in\mathbb{Z}$, as illustrated in Figure 1 (left).
We further assume that the interfaces across the edges of the periodic cells
are perfect, implying continuity for the displacement $U_{h}$ and the stress
$\hat{a}_{h}\frac{dU}{dX}$ there; the same is assumed within the unit cells at
$n\phi h$.

###  2.2 Non dimensionalization

To non-dimensionalize the physical problem, we introduce a reference wavespeed
$c_{0}=\sqrt{a_{0}/b_{0}}$ and the following non-dimensional quantities

$x=\frac{X}{L},\quad\delta=\frac{h}{L},\quad\Omega=\frac{\omega
h}{c_{0}},\quad\kappa=Lk,\quad u_{\delta}(x)=\frac{U_{h}(X)}{L}.$ (4)

Moreover, by periodicity

$\displaystyle\hat{b}_{h}(X,\omega)=\hat{b}\left(\frac{X}{h},\omega\right)$
and
$\displaystyle\hat{a}_{h}(X,\omega)=\hat{a}\left(\frac{X}{h},\omega\right),$

where $\hat{b}$ and $\hat{a}$ are 1-periodic in their first argument. These
physical quantities are non-dimensionalized by introducing

$\displaystyle
b\left(\frac{x}{\delta},\Omega\right)=\frac{\hat{b}\left(\frac{X}{h},\omega\right)}{b_{0}}$
and $\displaystyle
a\left(\frac{x}{\delta},\Omega\right)=\frac{\hat{a}\left(\frac{X}{h},\omega\right)}{a_{0}}.$

Using these quantities, (1) is rewritten as the non-dimensional governing
equation

$\displaystyle\delta^{2}\frac{\mathrm{d}}{\mathrm{d}x}\left(a\left(\frac{x}{\delta},\Omega\right)\frac{\mathrm{d}u_{\delta}}{\mathrm{d}x}(x)\right)+\Omega^{2}b\left(\frac{x}{\delta},\Omega\right)u_{\delta}(x)$
$\displaystyle=$ $\displaystyle 0.$ (5)

Upon introducing

$\mathbb{a}(\Omega)=\hat{\mathbb{a}}\left(\frac{c_{0}\Omega}{h}\right)\text{
and }\mathbb{b}(\Omega)=\hat{\mathbb{b}}\left(\frac{c_{0}\Omega}{h}\right),$

the adimensionalized physical parameters depend only on the short scale
$\xi=x/\delta$ and not on the long scale $x$ and become

$a(\xi,\Omega)=\left\\{\begin{aligned} &1&\text{ for }\xi\in(0,\phi)\\\
&\mathbb{a}(\Omega)&\text{ for }\xi\in(\phi,1)\end{aligned}\right.$ (6)
$b(\xi,\Omega)=\left\\{\begin{aligned} &1&\text{ for }\xi\in(0,\phi)\\\
&\mathbb{b}(\Omega)&\text{ for }\xi\in(\phi,1).\end{aligned}\right.$ (7)

We still have continuity for $u_{\delta}$ and
$a\frac{\mathrm{d}u_{\delta}}{\mathrm{d}x}$ at the points $n\delta$ and
$n\phi\delta$ for $n\in\mathbb{Z}$ in the geometry setting of Figure 1
(centre).

Figure 1: The three coordinate systems: (left) Physical coordinates, (middle)
non-dimensional long-scale, (right) non-dimensional short-scale.

###  2.3 Floquet-Bloch analysis

The periodicity of the parameters $a(\xi,\Omega)$ and $b(\xi,\Omega)$ defined
in (103)-(104), allows us to write the solution of (5), as
$u_{\delta}(x)=\mathfrak{u}_{\delta}(x)e^{\mathrm{i}\kappa x}$, for a
$\delta$-periodic function $\mathfrak{u}_{\delta}$ and Bloch wavenumber
$\kappa\in[0,\pi/\delta]$. For any Bloch wavenumber $\kappa$, this implies
that
$u_{\delta}\left({\delta^{+}}\right)=e^{\mathrm{i}\kappa\delta}u_{\delta}(0^{+})$
and
$u_{\delta}^{\prime}(\delta^{+})=u_{\delta}^{\prime}(0^{+})e^{\mathrm{i}\kappa\delta}$,
where we use the prime symbol for differentiation. Using perfect contact
conditions at $x=\delta$ and $x=\delta\phi$, the whole problem is reduced to
the unit cell $x\in(0,\delta)$ where (5) must be satisfied together with

$\displaystyle
u_{\delta}(\delta^{-})=e^{\mathrm{i}\kappa\delta}u_{\delta}(0^{+})$ and
$\displaystyle\mathbb{a}(\Omega)u_{\delta}^{\prime}(\delta^{-})=e^{\mathrm{i}\kappa\delta}u_{\delta}^{\prime}(0^{+}),$
(8)

as well as

$\displaystyle u_{\delta}(\delta\phi^{-})=u_{\delta}(\delta\phi^{+})$ and
$\displaystyle
u_{\delta}^{\prime}(\delta\phi^{-})=\mathbb{a}(\Omega)u_{\delta}^{\prime}(\delta\phi^{+}).$
(9)

The coefficients of (5) are piecewise constant (with respect to $\xi$) and so
we have two different equations on $(0,\delta\phi)$ and $(\delta\phi,\delta)$:

$\displaystyle\left\\{\begin{array}[]{rcc}u_{\delta}^{\prime\prime}(x)+\left(\frac{\Omega}{\delta}\right)^{2}u_{\delta}(x)=0&\text{for}&x\in(0,\delta\phi)\\\
u_{\delta}^{\prime\prime}(x)+\left(\frac{\Omega}{\delta}\right)^{2}\mathbb{g}(\Omega)u_{\delta}(x)=0&\text{for}&x\in(\delta\phi,\delta)\end{array}\right.,$
where
$\displaystyle\mathbb{g}(\Omega):=\frac{\mathbb{b}(\Omega)}{\mathbb{a}(\Omega)}\cdot$
(12)

Equations (12) have solution

$\displaystyle\left\\{\begin{array}[]{lc}u_{\delta}(x)=\mathfrak{C}_{1}\cos\left(\frac{\Omega
x}{\delta}\right)+\mathfrak{C}_{2}\sin\left(\frac{\Omega
x}{\delta}\right)&\text{ for $x\in(0,\delta\phi)$}\\\
u_{\delta}(x)=\mathfrak{C}_{3}\cos\left(\frac{\Omega\sqrt{\mathbb{g}(\Omega)}x}{\delta}\right)+\mathfrak{C}_{4}\sin\left(\frac{\Omega\sqrt{\mathbb{g}(\Omega)}x}{\delta}\right)&\text{
for $x\in(\delta\phi,\delta)$}\end{array}\right.\cdot$ (15)

Using the periodicity and interface conditions (8) and (9), the integration
constants $\mathfrak{C}_{i}$ satisfy the following linear system:

$\displaystyle\bm{\mathcal{M}}(\mathfrak{C}_{1},\mathfrak{C}_{2},\mathfrak{C}_{3},\mathfrak{C}_{4})^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}$
$\displaystyle=$
$\displaystyle(0,0,0,0)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},$
(16)

where the $4\times 4$ matrix $\bm{\mathcal{M}}$ is given by

$\displaystyle\left(\begin{array}[]{cccc}e^{\mathrm{i}\kappa\delta}&0&-\cos\left(\Omega\sqrt{\mathbb{g}(\Omega)}\right)&-\sin\left(\Omega\sqrt{\mathbb{g}(\Omega)}\right)\\\
0&e^{\mathrm{i}\kappa\delta}&\mathbb{a}(\Omega)\sqrt{\mathbb{g}(\Omega)}\sin\left(\Omega\sqrt{\mathbb{g}(\Omega)}\right)&-\mathbb{a}(\Omega)\sqrt{\mathbb{g}(\Omega)}\cos\left(\Omega\sqrt{\mathbb{g}(\Omega)}\right)\\\
\cos(\Omega\phi)&\sin(\Omega\phi)&-\cos\left(\Omega\sqrt{\mathbb{g}(\Omega)}\phi\right)&-\sin\left(\Omega\sqrt{\mathbb{g}(\Omega)}\phi\right)\\\
-\sin(\Omega\phi)&\cos(\Omega\phi)&\mathbb{a}(\Omega)\sqrt{\mathbb{g}(\Omega)}\sin\left(\Omega\sqrt{\mathbb{g}(\Omega)}\phi\right)&-\mathbb{a}(\Omega)\sqrt{\mathbb{g}(\Omega)}\cos\left(\Omega\sqrt{\mathbb{g}(\Omega)}\phi\right)\end{array}\right).$

Note that to get the second and fourth lines of $\bm{\mathcal{M}}$, we divided
through by $\Omega/\delta$. The system (16) has non-trivial solutions only
when $\bm{\mathcal{M}}$ is singular. Upon dividing through by
$e^{\mathrm{i}\kappa\delta}2\mathbb{a}(\Omega)\sqrt{\mathbb{g}(\Omega)}$, the
equation $\det(\bm{\mathcal{M}})=0$ reduces to the dispersion relation
$\text{Disp}(\Omega,\kappa)=0$, where

$\displaystyle\text{Disp}(\Omega,\kappa)$
$\displaystyle=\cos(\kappa\delta)-\cos(\Omega\phi)\cos\left(\Omega\sqrt{\mathbb{g}(\Omega)}(\phi-1)\right)$
(17)
$\displaystyle-\mathbb{d}(\Omega)\sin(\Omega\phi)\sin\left(\Omega\sqrt{\mathbb{g}(\Omega)}(\phi-1)\right),$

where the function $\mathbb{d}$ has been defined by

$\displaystyle\mathbb{d}(\Omega)$ $\displaystyle:=$
$\displaystyle\frac{1}{2}\left(\frac{1}{\mathbb{a}(\Omega)\sqrt{\mathbb{g}(\Omega)}}+\mathbb{a}(\Omega)\sqrt{\mathbb{g}(\Omega)}\right).$

Dispersion relations can usefully be thought of as nonlinear eigenvalue
problems and, in this context, it is known that if we take any open connected
domain in the $\Omega$ complex plane on which the entries of
$\bm{\mathcal{M}}$ are holomorphic, then there will be a finite (possibly
zero) number of isolated zeros of $\det(\bm{\mathcal{M}})$ within this domain
(see e.g. Theorem 2.1 in [38]). Therefore the same is true for the solution of
the dispersion relation. One should however be careful about domains that
contain points for which the entries of $\bm{\mathcal{M}}$ are singular (e.g.
branch points or poles). Given the form of $\bm{\mathcal{M}}$, these
potentially problematic points are values of $\Omega$ for which
$\mathbb{g}(\Omega)=0$ or $\mathbb{g}(\Omega)=\infty$. For spectral properties
of absorptive and dispersive photonic crystals, we refer the reader to [39]
and [40], respectively. More details are given in the next two paragraphs for
the case of the Lorentz model (see Appendix A for the expression of the
physical parameters in this case) as it is representative of issues that
arise.

##### Points for which $\mathbb{g}(\Omega)=0$ for the Lorentz model

Given that $\mathbb{g}(\Omega)=\mathbb{b}(\Omega)/\mathbb{a}(\Omega)$, these
are the points for which either $\mathbb{b}(\Omega)=0$ or
$\mathbb{a}(\Omega)=\infty$. So they are points that are solutions of (cf
equations (103) and (104))

$\displaystyle\left(1-\sum_{p\geq
0}\frac{\Omega_{p,i}^{2}}{\Omega(\Omega+\mathrm{i}\gamma_{p,i})-\Omega^{2}_{D,p,i}}\right)=0$
$\displaystyle\quad\text{ for $i=1$ or $i=2$}.$

If there is only one term in each sum, these points can be written down easily
explicitly, but otherwise for several terms in the sums they have to be found
numerically; finding these points is straightforward.

##### Points for which $\mathbb{g}(\Omega)=\infty$ for the Lorentz model

Given that $\mathbb{g}(\Omega)=\mathbb{b}(\Omega)/\mathbb{a}(\Omega)$, these
are the points for which either $\mathbb{b}(\Omega)=\infty$ or
$\mathbb{a}(\Omega)=0$ and are solutions of

$\displaystyle\left(1-\sum_{p\geq
0}\frac{\Omega_{p,i}^{2}}{\Omega(\Omega+i\gamma_{p,i})-\Omega^{2}_{D,p,i}}\right)=\infty$
$\displaystyle\quad\text{ for $i=1$ or $i=2$}.$

The points are found explicitly by nullifying the denominators of each term in
the sums and are given for $i=1,2$ and $p$ by
$\Omega=\mathcal{R}_{p,i}^{\pm}$, where

$\displaystyle\mathcal{R}_{p,i}^{\pm}$ $\displaystyle=$
$\displaystyle\frac{-\mathrm{i}\gamma_{p,i}\pm\sqrt{4\Omega^{2}_{D,p,i}-\gamma_{p,i}^{2}}}{2}\cdot$

In a neighbourhood of these points the theorem mentioned above does not apply,
and some of these will be accumulation points. In other words, if we take any
open connected set containing one of these points, it will contain infinitely
many zeros of the dispersion relation; this phenomenom is illustrated in
Figure 8b. It is interesting to note that these points are independent of the
choice of Bloch wavenumber $\kappa$. In the remaining parts of this paper, we
will aim to provide an asymptotic homogenised approximation to the dispersion
diagram and the corresponding wave field in the vicinity of an exact solution
$(\Omega_{0},\kappa)$ of the dispersion relation; our method works for points
$\Omega_{0}$ that are not too close to an accumulation point. In the vicinity
of accumulation points, another approach is required as some form of resonance
is expected [41, 42, 43, 44].

###  2.4 High-frequency homogenization

We assume that $\delta\ll 1$ and we recall that $\xi=x/\delta$. To start with,
we pick a frequency-wavenumber pair
$(\Omega_{0},\kappa)\in\mathbb{C}\times[0,\pi/\delta]$ that satisfies
$\text{Disp}(\Omega_{0},\kappa)=0$ and is such that we are not too close to an
accumulation point. Following the two-scale expansion technique, we further
assume the usual HFH ansatz for the wave field $u_{\delta}$ and the reduced
frequency $\Omega$:

$\displaystyle u_{\delta}(x)=\sum_{j\geqslant 0}\delta^{j}u_{j}(x,\xi)$ and
$\displaystyle\Omega^{2}=\sum_{\ell\geqslant
0}\delta^{\ell}\Omega_{\ell}^{2},$ (18)

where we treat $x$ and $\xi$ as two independent variables. The latter implies
that $\frac{\mathrm{d}}{\mathrm{d}x}\leftrightarrow\frac{\partial}{\partial
x}+\frac{1}{\delta}\frac{\partial}{\partial\xi}$. We will assume that

$u_{j}(x,\xi+1)=\mathrm{e}^{\mathrm{i}\kappa\delta}u_{j}(x,\xi),$ (19)

so we will restrict the analysis to (0,1) (see Figure 1 right).
Using this ansatz, and considering that the physical parameters are piecewise
constant, the governing equation (5) becomes

$\displaystyle\sum_{j\geqslant
0}\left[a(\xi,\Omega)\left\\{\delta^{j}\frac{\partial^{2}u_{j}}{\partial\xi^{2}}+2\delta^{j+1}\frac{\partial^{2}u_{j}}{\partial
x\partial\xi}+\delta^{j+2}\frac{\partial^{2}u_{j}}{\partial
x^{2}}\right\\}+\sum_{\ell\geqslant
0}\delta^{\ell+j}\Omega_{\ell}^{2}b(\xi,\Omega)u_{j}\right]=0.$ (20)

Importantly, $a$ and $b$ depend implicitly on $\delta$ through $\Omega$ in the
above expression. Therefore we have to write their expansion in powers of
$\delta$. Up to the second order, we find that for $q=\\{a,b\\}$:

$q(\xi,\Omega)=\mathcal{D}_{0}(q)+\delta\mathcal{D}_{1}(q)+\delta^{2}\mathcal{D}_{2}(q)$
(21)

with

$\displaystyle\mathcal{D}_{0}(q)=q(\xi,\Omega_{0}),$ (22)
$\displaystyle\mathcal{D}_{1}(q)=\frac{\Omega_{1}^{2}}{2\Omega_{0}}\
\frac{\partial q}{\partial\Omega}(\xi,\Omega_{0}),$ (23)
$\displaystyle\mathcal{D}_{2}(q)=\frac{\Omega_{1}^{4}}{8\Omega_{0}^{2}}\
\frac{\partial^{2}q}{\partial\Omega^{2}}(\xi,\Omega_{0})+\frac{1}{2\Omega_{0}}\left(\Omega_{2}^{2}-\frac{\Omega_{1}^{4}}{4\Omega_{0}}\right)\
\frac{\partial q}{\partial\Omega}(\xi,\Omega_{0}).$ (24)

We define $A_{i}=\mathcal{D}_{i}(a)$ and $B_{i}=\mathcal{D}_{i}(b)$ for
$i=\\{0,1,2\\}$, where we have to keep in mind that $A_{i}$ and $B_{i}$ depend
on $\xi$ and $\left\\{\Omega_{0},\cdots,\Omega_{i}\right\\}$. For the
dispersion model we chose their expressions are detailed in Appendix A.
We also need to introduce the average operator $\langle\cdot\rangle$ defined
by

$\langle g\rangle=\int_{0}^{1}g(\xi)\mathrm{d}\xi$

for any function $g$.
In the next sections, we will apply high-frequency homogenization to get
asymptotic approximations of both the wavefields and the dispersion diagrams
for all the relevant cases which are represented in Figure 2.

Figure 2: Different cases approximated with high-frequency homogenization.

####  2.4.1 Single eigenvalues at the edges of the Brillouin zone

We first take points at the edges of the Brillouin zone, that is at $\kappa=0$
or $\kappa=\pi/\delta$ which is associated to periodic and antiperiodic
conditions, respectively, for the fields. In these cases, we expect the
dispersion relation to be locally quadratic i.e. we set $\Omega_{1}=0$.
Indeed, the mapping $\Omega_{n}(k)$ denoting the dispersion diagram along a
branch is holomorphic [38] except around accumulation points and at singular
points (when eigenvalues are no longer single ones); combined with
reciprocity, i.e. $\Omega_{n}(k)=\Omega_{n}(-k)$, this gives that
$\mathrm{d}\Omega_{n}/\mathrm{d}k=0$ at edges.

##### Zeroth-order field

Collecting the terms of order $\delta^{0}$, we get in $(0,1)$:

$A_{0}(\xi,\Omega_{0})\frac{\partial^{2}u_{0}}{\partial\xi^{2}}(x,\xi)+\Omega_{0}^{2}B_{0}(\xi,\Omega_{0})u_{0}(x,\xi)=0$
(25)

where $A_{0}$ and $B_{0}$ are piecewise constants defined in (22). We also
have continuity for $u_{0}$ and $A_{0}\frac{\partial u_{0}}{\partial\xi}$ at
$\phi$ and $0$ together with the 1-periodicity/antiperiodicity for $u_{0}$:

$\left\\{\begin{aligned} u_{0}(\phi^{-})&=u_{0}(\phi^{+})\\\ \frac{\partial
u_{0}}{\partial\xi}(\phi^{-})&=\mathbb{a}(\Omega_{0})\frac{\partial
u_{0}}{\partial\xi}(\phi^{+})\\\ u_{0}(1^{-})&=\pm u_{0}(0^{+})\\\
\mathbb{a}(\Omega_{0})\frac{\partial
u_{0}}{\partial\xi}(1^{-})&=\pm\frac{\partial
u_{0}}{\partial\xi}(0^{+}).\end{aligned}\right.$ (26)

As discussed in Section 2.3, we will build asymptotic approximations
sufficiently far away from the accumulation points. We therefore know that
there is a discrete set of eigenvalues and then choose $\Omega_{0}$ which is
assumed to be a simple eigenvalue associated to the eigenfunction
$\mathcal{U}_{0}(\xi,\Omega_{0})$. The zeroth-order field is therefore

$u_{0}(x,\xi)=f_{0}(x)\mathcal{U}_{0}(\xi,\Omega_{0}),$ (27)

where the slowly varying amplitude $f_{0}(x)$ has to be determined.

##### First-order field

Collecting the terms of order $\delta$, and keeping in mind that
$\Omega_{1}=0$, we get on $(0,1)$:

$A_{0}\frac{\partial}{\partial\xi}\left(\frac{\partial
u_{1}}{\partial\xi}+2\frac{\partial u_{0}}{\partial
x}\right)+\Omega_{0}^{2}B_{0}u_{1}=0$ (28)

together with periodicity/antiperiodicity for $u_{1}$, continuity for $u_{1}$
at $\phi$ and $0$, and continuity for $A_{0}\left(\frac{\partial
u_{0}}{\partial x}+\frac{\partial u_{1}}{\partial\xi}\right)$ at $\phi$ and
$0$. Using (27), Equation(28) reduces to:

$A_{0}\frac{\partial^{2}u_{1}}{\partial\xi^{2}}+\Omega_{0}^{2}B_{0}u_{1}=-2A_{0}\mathcal{U}^{\prime}_{0}(\xi)f^{\prime}_{0}(x).$
(29)

Then, we write $u_{1}$ as:

$u_{1}(x,\xi)=f_{1}(x)\mathcal{U}_{0}(\xi)+f_{0}^{\prime}(x)V(\xi),$ (30)

where

$\left\\{\begin{aligned}
&A_{0}V^{\prime\prime}(\xi)+\Omega_{0}^{2}B_{0}V(\xi)=-2A_{0}\mathcal{U}^{\prime}_{0}(\xi)\\\
&V\text{ is periodic/antiperiodic }\\\ &\text{$V$ and
$A_{0}V^{\prime}+A_{0}\mathcal{U}_{0}$ are continuous at $\phi$ and 0
}.\end{aligned}\right.$ (31)

##### Second-order field

Collecting terms of order $\delta^{2}$, we get on $(0,1)$:

$A_{0}\left(\frac{\partial^{2}u_{0}}{\partial
x^{2}}+2\frac{\partial^{2}u_{1}}{\partial
x\partial\xi}+\frac{\partial^{2}u_{2}}{\partial\xi^{2}}\right)+A_{2}\frac{\partial^{2}u_{0}}{\partial\xi^{2}}+\Omega_{0}^{2}(B_{0}u_{2}+B_{2}u_{0})+\Omega_{2}^{2}B_{0}u_{0}=0,$
(32)

together with periodicity/antiperiodicity for $u_{2}$, continuity for $u_{2}$
at $\phi$ and $0$, and continuity for $A_{0}\left(\frac{\partial
u_{2}}{\partial\xi}+\frac{\partial u_{1}}{\partial
x}\right)+A_{2}\frac{\partial u_{0}}{\partial\xi}$ at $\phi$ and $0$. Consider
now the equation

$\langle{u_{2}}\times\eqref{syst_order0_main}-u_{0}\times{\eqref{syst_order2_main}}\rangle=0.$

After integration by part, some algebra and dividing through by $f_{0}$ we
get:

$\langle-
A_{0}\mathcal{U}_{0}^{2}-A_{0}V^{\prime}\mathcal{U}_{0}+A_{0}V\mathcal{U}^{\prime}_{0}\rangle
f_{0}^{\prime\prime}+\langle-\Omega_{2}^{2}B_{0}\mathcal{U}_{0}^{2}-\Omega_{0}^{2}B_{2}\mathcal{U}_{0}^{2}+A_{2}(\mathcal{U}_{0}^{\prime})^{2}\rangle
f_{0}=0.$ (33)

Furthermore, note that when $\Omega_{1}=0$, $A_{2}$ and $B_{2}$ defined in
(24), regardless of the dispersive model chosen, is:

$B_{2}(\Omega_{0},0,\Omega_{2})=\Omega_{2}^{2}\tilde{B}_{2}(\Omega_{0})\text{
and }A_{2}(\Omega_{0},0,\Omega_{2})=\Omega_{2}^{2}\tilde{A}_{2}(\Omega_{0}).$
(34)

Therefore, we get the sought-after effective equation for $f_{0}$ and
$\Omega_{2}$:

$Tf_{0}^{\prime\prime}+\Omega_{2}^{2}f_{0}=0,$ (35)

where

$T=\frac{\langle
A_{0}\mathcal{U}_{0}^{2}+A_{0}V^{\prime}\mathcal{U}_{0}-A_{0}V\mathcal{U}^{\prime}_{0}\rangle}{\langle
B_{0}\mathcal{U}_{0}^{2}+\Omega_{0}^{2}\tilde{B}_{2}\mathcal{U}_{0}^{2}-\tilde{A}_{2}(\mathcal{U}_{0}^{\prime})^{2}\rangle}=\frac{\left\langle
a(\cdot,\Omega_{0})(\mathcal{U}_{0}^{2}+V^{\prime}\mathcal{U}_{0}-V\mathcal{U}^{\prime}_{0})\right\rangle}{\mathcal{S}(\Omega_{0},\mathcal{U}_{0},\mathcal{U}_{0})},$
(36)

for any frequency-dependent functions $a$ and $b$, and where for any two
functions $f$ and $g$, and any reduced frequency $\Omega_{0}$, $\mathcal{S}$
is defined by:

$\mathcal{S}(\Omega_{0},f,g)=\left\langle\left(b(\cdot,\Omega_{0})+\frac{\Omega_{0}}{2}\frac{\partial
b}{\partial\Omega}(\cdot,\Omega_{0})\right)fg-\frac{1}{2\Omega_{0}}\frac{\partial
a}{\partial\Omega}(\cdot,\Omega_{0})f^{\prime}g^{\prime}\right\rangle.$ (37)

Therefore, we get the effective string described by (35) on the long-scale
where the complex material properties are now solely concentrated in a single
effective parameter $T$ (36).
Applying the Bloch-Floquet conditions (19) gives the final expression for the
quadratic term in the dispersion relation

$\Omega_{2}^{2}=T\tilde{\kappa}^{2},$ (38)

with $\tilde{\kappa}=\kappa$ near 0 and $\tilde{\kappa}=\pi/\delta-\kappa$
near $\pi/\delta$ so that the dispersion relation is approximated by

$\Omega\approx\Omega_{0}+\frac{T}{2\Omega_{0}}(\tilde{\kappa}\delta)^{2}.$
(39)

####  2.4.2 Double eigenvalues at the edges of the Brillouin zone

We next consider the case of multiplicity two for the eigenvalue $\Omega_{0}$.
In that case, $\Omega_{1}$ is no longer 0, and the zeroth-order wavefield is
now written as

$u_{0}(x,\xi)=f_{0}^{(1)}(x)\mathcal{U}_{0}^{(1)}(\xi)+f_{0}^{(2)}(x)\mathcal{U}_{0}^{(2)}(\xi)$
(40)

where $\mathcal{U}_{0}^{(1)}(\xi)$ and $\mathcal{U}_{0}^{(2)}(\xi)$ are two
independent eigenfunctions associated to $\Omega_{0}$, while $f_{0}^{(1)}(x)$
and $f_{0}^{(2)}(x)$ are the slow modulation functions to be found.
Consequently, both eigenfunctions satisfy (25) and we denote
$\eqref{syst_order0_main}^{(i)}$ the equation for the $i$th eigenfunction.
Furthermore, because $\Omega_{1}$ is non-zero, the system satisfied by the
first-order field becomes

$A_{0}\frac{\partial}{\partial\xi}\left(\frac{\partial
u_{1}}{\partial\xi}+2\frac{\partial u_{0}}{\partial
x}\right)+A_{1}\frac{\partial^{2}u_{0}}{\partial\xi^{2}}+\Omega_{0}^{2}(B_{0}u_{1}+B_{1}u_{0})+\Omega_{1}^{2}B_{0}u_{0}=0,$
(41)

together with periodicity/antiperiodicity for $u_{1}$, continuity for $u_{1}$
at $\phi$ and $0$, and continuity for $A_{0}\left(\frac{\partial
u_{0}}{\partial x}+\frac{\partial
u_{1}}{\partial\xi}\right)+A_{1}\frac{\partial u_{0}}{\partial\xi}$ at $\phi$
and $0$.

We introduce $\tilde{B}_{1}$ and $\tilde{A}_{1}$ so that

$B_{1}(\Omega_{0},\Omega_{1})=\frac{\Omega_{1}^{2}}{\Omega_{0}^{2}}\tilde{B}_{1}(\Omega_{0})\text{
and
}A_{1}(\Omega_{0},\Omega_{1})=\frac{\Omega_{1}^{2}}{\Omega_{0}^{2}}\tilde{A}_{1}(\Omega_{0}).$
(42)

Then considering the equations $\langle\eqref{syst_order0_main}^{(i)}\times
u_{1}-\eqref{syst_order1_main_Omega1_non_zero}\times\mathcal{U}_{0}^{(i)}\rangle=0$
for $i=1,2$, we get the effective equation for
$\bm{F}=(f_{0}^{(1)},f_{0}^{(2)})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}$:

$\bm{F}^{\prime}(x)=\frac{\Omega_{1}^{2}}{\langle
A_{0}\mathrm{w}_{0}\rangle}\bm{\mathcal{N}}\bm{F},$ (43)

where $\mathrm{w}_{0}$ is the Wronskian defined by

$\mathrm{w}_{0}(\xi)=\mathcal{U}_{0}^{(1)}(\xi){\mathcal{U}_{0}^{(2)}}^{\prime}(\xi)-{\mathcal{U}_{0}^{(1)}}^{\prime}(\xi)\mathcal{U}_{0}^{(2)}(\xi),$
(44)

$\bm{\mathcal{N}}$ is the matrix defined by

$\bm{\mathcal{N}}=\left(\begin{array}[]{cc}\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)})&\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(2)},\mathcal{U}_{0}^{(2)})\\\
-\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(1)})&-\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)})\end{array}\right),$
(45)

and $\mathcal{S}$ is defined in (37).

Regarding the dispersion diagram, using Bloch-Floquet conditions (19) we get
the two opposite slopes (here the upper and lower notation does not stand for
left or right edge of the Brillouin zone but for the upper and lower branch
starting from $\Omega_{0}$)

$\Omega_{1}^{2}=\pm T_{D}\tilde{\kappa}$ (46)

so that

$\Omega\approx\Omega_{0}\pm\frac{T_{D}}{2\Omega_{0}}\tilde{\kappa}\delta$ (47)

with $\tilde{\kappa}=\kappa$ near 0 and $\tilde{\kappa}=\pi/\delta-\kappa$
near $\pi/\delta$, and $T_{D}$ defined by

$\displaystyle T_{D}=\frac{\langle
a(\cdot,\Omega_{0})\mathrm{w}_{0}\rangle}{\left(\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(1)})\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(2)},\mathcal{U}_{0}^{(2)})-\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)})^{2}\right)^{1/2}}.$
(48)

####  2.4.3 Nearby eigenvalues at the edges of the Brillouin zone

Let us assume now that we have two simple eigenvalues close to each other
following [35, 24]. The two nearby simple eigenvalues are denoted by
$\Omega_{0}^{(1)}$ and $\Omega_{0}^{(2)}$. Their proximity is quantified by
writing

$(\Omega_{0}^{(2)})^{2}-(\Omega_{0}^{(1)})^{2}=\alpha\delta,$ (49)

for some constant $\alpha>0$. We denote by $\mathcal{U}_{0}^{(1)}$ and
$\mathcal{U}_{0}^{(2)}$ the eigenfunctions associated to $\Omega_{0}^{(1)}$
and $\Omega_{0}^{(2)}$, respectively.
The ansatz is considered around the eigenvalue $\Omega_{0}^{(1)}$:

$\displaystyle u_{\delta}(x)=\sum_{j\geqslant 0}\delta^{j}u_{j}(x,\xi)$ and
$\displaystyle\Omega^{2}=(\Omega_{0}^{(1)})^{2}+\sum_{\ell\geqslant
1}\delta^{\ell}\Omega_{\ell}^{2},$ (50)

By similarity with the double eigenvalue case, and to take into account the
coupling between both eigenvalues we look for the zeroth-order wavefield as

$u_{0}(x,\xi)=f_{0}^{(1)}(x)\mathcal{U}_{0}^{(1)}(\xi)+f_{0}^{(2)}(x)\mathcal{U}_{0}^{(2)}(\xi).$
(51)

We will use the notation $D_{j}^{(i)}$, for $D=A,\,B$ and $i=1,\,2$, and
$j=0,\,1$ so that we get from Taylor expansions

$A_{0}^{(i)}\frac{\partial^{2}\mathcal{U}_{0}^{(i)}}{\partial\xi^{2}}+(\Omega_{0}^{(i)})^{2}B_{0}^{(i)}\mathcal{U}_{0}^{(i)}=0$
(52)

and

$D_{0}^{(2)}=D_{0}^{(1)}+\frac{\alpha\delta}{2\Omega_{0}^{(1)}}\frac{\partial
d}{\partial\Omega}(\Omega_{0}^{(1)})+\mathcal{O}(\delta^{2})$ (53)

where $d=a,\,b$ for $D=A,\,B$, respectively. This leads to

$\begin{multlined}A_{0}^{(1)}\frac{\partial^{2}u_{0}}{\partial\xi^{2}}+(\Omega_{0}^{(1)})^{2}B_{0}^{(1)}u_{0}=\delta\left[\frac{\alpha}{2A_{0}^{(1)}}\Omega_{0}^{(1)}B_{0}^{(1)}\frac{\partial
a}{\partial\Omega}(\xi,\Omega_{0}^{(1)})\right.\\\ \left.-\alpha
B_{0}^{(1)}-\frac{\alpha}{2}\Omega_{0}^{(1)}\frac{\partial
b}{\partial\Omega}(\xi,\Omega_{0}^{(1)})\right]u_{0}^{(2)}.\end{multlined}A_{0}^{(1)}\frac{\partial^{2}u_{0}}{\partial\xi^{2}}+(\Omega_{0}^{(1)})^{2}B_{0}^{(1)}u_{0}=\delta\left[\frac{\alpha}{2A_{0}^{(1)}}\Omega_{0}^{(1)}B_{0}^{(1)}\frac{\partial
a}{\partial\Omega}(\xi,\Omega_{0}^{(1)})\right.\\\ \left.-\alpha
B_{0}^{(1)}-\frac{\alpha}{2}\Omega_{0}^{(1)}\frac{\partial
b}{\partial\Omega}(\xi,\Omega_{0}^{(1)})\right]u_{0}^{(2)}.$ (54)

Consequently, (25) is satisfied by the zeroth-order wavefield up to right-hand
side residual term of (54) that modifies the equation for the first order
field that now becomes:

$\displaystyle A_{0}^{(1)}\frac{\partial}{\partial\xi}\left(\frac{\partial
u_{1}}{\partial\xi}+2\frac{\partial u_{0}}{\partial
x}\right)+A_{1}^{(1)}\frac{\partial^{2}u_{0}}{\partial\xi^{2}}+(\Omega_{0}^{(1)})^{2}(B_{0}^{(1)}u_{1}+B_{1}^{(1)}u_{0})+\Omega_{1}^{2}B_{0}^{(1)}u_{0}$
(55) $\displaystyle+\alpha
u_{0}^{(2)}\left[\frac{1}{2A_{0}^{(1)}}\Omega_{0}^{(1)}B_{0}^{(1)}\frac{\partial
a}{\partial\Omega}(\xi,\Omega_{0}^{(1)})-B_{0}^{(1)}-\frac{1}{2}\Omega_{0}^{(1)}\frac{\partial
b}{\partial\Omega}(\xi,\Omega_{0}^{(1)})\right]=0,$

together with periodicity/antiperiodicity for $u_{1}$, continuity for $u_{1}$
at $\phi$ and $0$, and continuity for $A_{0}^{(1)}\left(\frac{\partial
u_{0}}{\partial x}+\frac{\partial
u_{1}}{\partial\xi}\right)+A_{1}^{(1)}\frac{\partial u_{0}}{\partial\xi}$ at
$\phi$ and $0$. As in the double eigenvalue case, considering
$\langle\eqref{zero_eigenfunction_nearby}^{(i)}\times
u_{1}-\eqref{syst_order1_nearby}\times\mathcal{U}_{0}^{(i)}\rangle$ for
$i=1,2$ allows to obtain the effective equation for
$\bm{F}=(f_{0}^{(1)},f_{0}^{(2)})^{\mathrm{T}}$:

$\bm{F}^{\prime}(x)=\frac{\Omega_{1}^{2}}{\langle
A_{0}^{(1)}\mathrm{w}_{0}\rangle}\bm{\mathcal{N}}_{\alpha}\bm{F}(x)$ (56)

with $\mathrm{w}_{0}$ still given by (44). The matrix
$\bm{\mathcal{N}}_{\alpha}$ is defined by

$\bm{\mathcal{N}}_{\alpha}=\left(\begin{array}[]{cc}\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)})&\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(2)},\mathcal{U}_{0}^{(2)})\\\
-\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(1)})&-\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)})\end{array}\right)+\frac{\alpha}{\Omega_{1}^{2}}\left(\begin{array}[]{cc}0&-\mathcal{G}_{2}\\\
0&\mathcal{G}_{1}\end{array}\right)$ (57)

with

$\displaystyle\mathcal{G}_{i}$
$\displaystyle=\left\langle(\tilde{B}_{1}^{(1)}+B_{0}^{(1)})\mathcal{U}_{0}^{(i)}\mathcal{U}_{0}^{(2)}-\frac{B_{0}^{(1)}}{A_{0}^{(1)}}\tilde{A}_{1}^{(1)}\mathcal{U}_{0}^{(i)}\mathcal{U}_{0}^{(2)}\right\rangle$
(58) $\displaystyle=\left\langle
b(\cdot,\Omega_{0}^{(1)})\left(1+\frac{\Omega_{0}^{(1)}}{2}\frac{\partial}{\partial\Omega}[\log(b/a)](\cdot,\Omega_{0}^{(1)})\right)\mathcal{U}_{0}^{(i)}\mathcal{U}_{0}^{(2)}\right\rangle.$

The dispersion relation is then obtained by solving

$\mathrm{det}\left(\frac{\Omega_{1}^{2}}{\langle
A_{0}^{(1)}\mathrm{w}_{0}\rangle}\bm{\mathcal{N}}_{\alpha}-\mathrm{i}\kappa\mathbb{I}_{2}\right)=0.$
(59)

One notes that $\bm{\mathcal{N}}_{\alpha}=\bm{\mathcal{N}}$ when $\alpha=0$
and we then recover the double case.

####  2.4.4 Simple eigenvalues inside the Brillouin zone (no damping)

In this section, we get a linear approximation for an arbitrary point inside
(strictly) the Brillouin zone. However, this is possible only if the physical
parameters are real. Consequently, we consider all the damping terms
$\gamma_{p,1}$ and $\gamma_{p,2}$ equal to 0, in the framework of this
subsection only, so that the coefficients $A_{i}$ and $B_{i}$ ($i=1,2$) are
real.
Let us pick a point $(\kappa^{\star},\Omega_{0})$ with
$\kappa^{\star}\notin\left\\{0,\pi/\delta\right\\}$ and $\Omega_{0}$ solution
of the eigenvalue problem satisfied by the zeroth-order field:

$A_{0}\frac{\partial^{2}u_{0}}{\partial\xi^{2}}+\Omega_{0}^{2}B_{0}u_{0}=0\text{
in }(0,\delta)\times(0,1),$ (60)

together with
$u_{0}(x,\xi+1)=\mathrm{e}^{\mathrm{i}\kappa^{\star}\delta}u_{0}(x,\xi)$ in
$(0,\delta)\times(0,1)$, continuity for $u_{0}$ at $\xi=\phi$ and $\xi=0$, and
continuity for $A_{0}\frac{\partial u_{0}}{\partial\xi}$ at $\xi=\phi$ and
$\xi=0$.
For the first order, we get:

$A_{0}\frac{\partial}{\partial\xi}\left(\frac{\partial
u_{1}}{\partial\xi}+2\frac{\partial u_{0}}{\partial
x}\right)+A_{1}\frac{\partial^{2}u_{0}}{\partial\xi^{2}}+\Omega_{0}^{2}(B_{0}u_{1}+B_{1}u_{0})+\Omega_{1}^{2}B_{0}u_{0}=0\text{
in }(0,\delta)\times(0,1),$ (61)

together with
$u_{1}(x,\xi+1)=\mathrm{e}^{\mathrm{i}\kappa^{\star}\delta}u_{1}(x,\xi)$ in
$(0,\delta)\times(0,1)$, continuity for $u_{1}$ at $\xi=\phi$ and $\xi=0$, and
continuity for $A_{0}\left(\frac{\partial u_{0}}{\partial x}+\frac{\partial
u_{1}}{\partial\xi}\right)+A_{1}\frac{\partial u_{0}}{\partial\xi}$ at
$\xi=\phi$ and $\xi=0$.
We consider
$I=\langle\bar{u}_{1}\times\eqref{inside_1D_zero}-\bar{\eqref{inside_1D_first}}\times
u_{0}\rangle=0$. We still have

$\left\langle\frac{\partial}{\partial\xi}\left[A_{0}\left(\frac{\partial\bar{u}_{0}}{\partial
x}+\frac{\partial\bar{u}_{1}}{\partial\xi}\right)u_{0}+A_{1}\frac{\partial\bar{u}_{0}}{\partial\xi}u_{0}-A_{0}\frac{\partial
u_{0}}{\partial\xi}\bar{u}_{1}\right]\right\rangle=0$ (62)

because the quantities above are continuous and 1-periodic. Therefore, $I=0$
reduces to

$\left\langle A_{0}\left(\frac{\partial^{2}\bar{u}_{0}}{\partial\xi\partial
x}u_{0}-\frac{\partial\bar{u}_{0}}{\partial x}\frac{\partial
u_{0}}{\partial\xi}\right)-A_{1}\left|\frac{\partial
u_{0}}{\partial\xi}\right|^{2}+\Omega_{0}^{2}B_{1}|u_{0}|^{2}+\Omega_{1}^{2}B_{0}|u_{0}|^{2}\right\rangle=0.$
(63)

Dividing through by $f_{0}$ and using (60), we end up with a first-order ODE
for $f_{0}$:

$T_{i}f_{0}^{\prime}(x)-\mathrm{i}\Omega_{1}^{2}f_{0}(x)=0$ (64)

with $T_{i}$ defined by

$T_{i}=\frac{2\langle
a(\cdot,\Omega_{0})\Im(\mathcal{U}_{0}^{\prime}\bar{\mathcal{U}}_{0})\rangle}{\mathcal{S}(\Omega_{0},\mathcal{U}_{0},\bar{\mathcal{U}}_{0})},$
(65)

where we remind that $S$ is defined in (37). Regarding the dispersion
relation, applying the Bloch-Floquet conditions (19) gives

$\Omega_{1}^{2}=T_{i}(\kappa-\kappa^{\star}).$ (66)

####  2.4.5 The low-frequency case

We also obtain the classical low-frequency homogenization by considering
$\Omega_{0}=0$ in (25), which leads to the fact that $\mathcal{U}_{0}$ is
uniform, and we, without loss of generality, choose $\mathcal{U}_{0}=1$. Then,
we write $u_{1}(x,\xi)=f_{0}^{\prime}(x)\mathcal{U}_{1}(\xi)$ where $u_{1}$
satisfies (28) with $\Omega_{0}=0$ and $u_{0}=f_{0}$. Therefore $\langle
A_{0}(1+\mathcal{U}_{1}^{\prime})\rangle=\langle 1/A_{0}\rangle^{-1}$.
Integration on a unit cell of (32) for $\Omega_{0}=0$ then leads to the usual
homogenized equation

$\langle 1/A_{0}\rangle^{-1}f^{\prime\prime}_{0}(x)+\Omega_{2}^{2}\langle
B_{0}\rangle f_{0}(x)=0$ (67)

and dispersion relation

$\Omega_{2}^{2}=\langle 1/A_{0}\rangle^{-1}\langle
B_{0}\rangle^{-1}\kappa^{2}.$ (68)

###  2.5 Numerical investigation

We now use two different methods to compute the whole dispersion diagrams: we
either track the zeros of the dispersion function (17) in the complex plane
along a branch, or we use the finite element method (FEM) to directly solve
for (12). The details are given in Appendix B for the latter. Hereafter, the
dispersion diagrams computed either by zero tracking or by FEM from the exact
dispersion function will be denoted as the exact dispersion diagrams in
contrast to the asymptotic approximations obtained by HFH with which they will
be compared.

####  2.5.1 Dielectric and metallic layers (Drude with damping)

Motivated by the configuration of [45], we consider wave propagation through
alternated layers of silver (Ag) and Titanium dioxide (TiO2). Here only the
permittivity is frequency-dependent, following a Drude law in the metal layer
made of silver (Ag). More precisely, we have $a=1$ in both materials, and in
silver

$\mathbb{b}(\Omega)=\frac{b_{\infty}}{b_{0}}\left(1-\frac{\Omega_{1,2}^{2}}{\Omega(\Omega+\mathrm{i}\gamma_{1,2})}\right)$

with $b_{0}=6.2$, $b_{\infty}=1$, and one resonance for $b$ with
$\Omega_{1,2}=5.01$, and $\gamma_{1,2}=0.01$ in (104). The filling ratio of
the dielectric layer is $\phi=10/11$ for a periodicity of $h=110$ nm.
Since there is damping, the frequency solutions of the dispersion relation are
complex and to visualize the dispersion function and its zeros without having
to consider the real and imaginary parts separately, we plot at a given
frequency its phase portrait in the complex plane, see Figure 3a for
$\kappa=0$. Alternatively, we track the zeros of the function along a branch
of the dispersion diagram in the complex plane, see Figure 3b.

(a) Phase portrait in complex plane for the dispersion (i.e. determinant)
function (17) taken at $\kappa=0$

(b) Dispersion diagram in the complex plane
(here $\kappa\in(0\pi/\delta)$)

Figure 3: Representation of the dispersion function and its zeros.

##### Simple eigenvalue approximations at the edges

Firstly, we use the asymptotic approximations obtained for simple eigenvalues
near the edges of the Brillouin zone, see Section 2.4.1. The resulting
asymptotic approximations of the dispersion diagrams for both the real part
and the imaginary part are displayed in Figure 4, where we used the zero
tracking method to compute the diagrams for the exact dispersion relation. The
absolute errors for each of the branches and for both $\kappa=0$ and
$\kappa=\pi/\delta$ are then shown in Figure 5, where we recover that the
quadratic term $\Omega_{2}$ is well taken into account asymptotically and that
in fact the next term is zero so that the error is
$\mathcal{O}((\tilde{\kappa}\delta)^{4})$.

(a) Real part (b) Imaginary part

Figure 4: Superposition of the dispersion diagram from zero tracking (solid
lines) and the asymptotic approximations obtained by HFH in the simple
eigenvalue case (39) around $\tilde{\kappa}=0$ (dashed lines).

(a) Near $\kappa\delta=0$ (b) Near $\kappa\delta=\pi$

Figure 5: Absolute error between the exact dispersion diagram and the HFH
approximation in the simple eigenvalue case (39) in a log-log scale (solid
lines). Dashed lines are reference orders of convergence.

##### Nearby approximations at the edges

Given these numerical solutions, we now compare with the asymptotics and use
the nearby approximations developed in Section 2.4.3 for single eigenvalues
near the edges of the Brillouin zone. Solutions from the FEM method are shown
in Figure 6 for real and imaginary parts for four modes, and in Figure 7 for
the dispersion diagram in the complex plane. It is clearly seen that the
agreement with the dispersion diagram is much longer lived than that of the
simple eigenvalue approximation for both real and imaginary parts. The nearby
approximation leads notably to a better fit of the imaginary parts, which are
quite small due to the fact that the damping coefficients $\gamma_{i}$ are
also small in practice.

(a) Real part for the first and second mode around $\kappa\delta=\pi$. (b)
Imaginary part for the first and second mode around $\kappa\delta=\pi$.

(c) Real part for the second and third mode around $\kappa\delta=0$. (d)
Imaginary part for the second and third mode around $\kappa\delta=0$.

Figure 6: Comparison of the exact dispersion relation (solid lines) and the
effective one obtained by HFH in the simple eigenvalue case (39) (dashed
lines) and with nearby approximations (59) (dotted lines). Figure 7:
Comparison in the complex plane of the exact dispersion relation (solid lines)
and the effective one obtained by HFH for the single eigenvalue case (39)
(dashed lines) and with nearby approximations (59) (dotted lines). The circle
markers correspond to $\kappa\delta=0$.

####  2.5.2 Stack of positive and negative index materials (Lorentz with no
damping)

In a second more challenging example, we now reproduce the results of Li et
al. [4], see Figure 2 of the latter and compare with the asymptotic results.
This consists of a 1D system of periodicity $h=18$ mm with alternate layers of
air (12 mm thick) and of an effective material which is dispersive (6.0 mm
thick). Both parameters in the dispersive medium follow a Lorentz law without
damping in the effective layers, see (103) and (104). More precisely, we set
$\phi=2/3$, $a_{\infty}=a_{0}=b_{\infty}=b_{0}=1$, $\Omega_{1,1}=1.131$,
$\Omega_{D,1,1}=0.34$, $\Omega_{1,2}=1.885$, $\Omega_{D,1,2}=0.3393$,
$\Omega_{2,2}=3.7699$ and $\Omega_{D,2,2}=4.3354$ resulting in

$\displaystyle\mathbb{a}(\Omega)=\left(1+\frac{\Omega_{1,1}^{2}}{\Omega^{2}-\Omega_{D,1,1}^{2}}\right)^{-1}$
(69)
$\displaystyle\mathbb{b}(\Omega)=1+\frac{\Omega_{1,2}^{2}}{\Omega^{2}-\Omega_{D,1,2}^{2}}+\frac{\Omega_{2,2}^{2}}{\Omega^{2}-\Omega_{D,2,2}^{2}}.$

We first compute the dispersion relation using Bloch-Floquet analysis (see
Section 2.3). The logarithm of the dispersion function (17), is plotted in
Figure 8a, the dispersion curve therefore corresponds to the dark lines in the
map. The main features of Figure 2 in [4] are recovered, together with the
appearance of accumulation points, see Figure 8b for phase portraits zoomed-in
around one of these points. However, we will consider the same range of
frequencies as in [4], for which we are away from any of these points and able
to propose high-frequency homogenized approximations.

(a) Map in $(\kappa,\Omega)$ space of the logarithm

(b) Phase portrait around an accumulation point at $\kappa=0$

Figure 8: Dispersion function. (Left) Map of the logarithm of the dispersion
function in the wavenumber-frequency space. The zeros are the solutions of the
dispersion relation, therefore represented by the dark lines. (Right) Phase
portrait in the complex plan zoomed around the accumulation point which occurs
at $\mathcal{R}^{+}_{1,1}=\Omega_{D,1,1}=4.3354$, as predicted in Section 2.3

##### Simple eigenvalue approximations at the edges

We then recover this band diagram near the edges using the quadratic
approximations obtained with HFH in Section 2.4.1. The comparison is given in
Figure 9a where the numerical curves are obtained using zero tracking; the
branches near the edges are well approximated. For a quantitative validation,
we plot the difference between the exact dispersion relation and the quadratic
approximation on a log-log scale. It validates the approximation of the
quadratic term and again underlines that there is no third-order term because
we get an error of order $\mathcal{O}((\delta\tilde{\kappa})^{4})$, see Figure
9b for the case $\kappa=0$.

Obtaining accurate asymptotics for the dispersion curves is a useful
application of the theory, by validating it and by encapsulating the physics
into a coefficient $T$ that allows us to tune or design features. An equally
important application of the theory is to model forcing, that is, to apply a
source in a structured medium and then use the effective equations to model
the response; we now proceed to demonstrate the efficiency of that approach.
We introduce a source term and choose a frequency of excitation close to an
eigenfrequency at $\kappa=0$, and then compare the wavefields for the
microstructured medium using both numerical simulation and the high-frequency
approximations in Figure 10. We first use a point source spatially located at
$\xi=0.8$ and repeated periodically, with frequency
$\Omega=1.01\Omega_{0}^{(4)}=1.027$, close to the fourth eigenfrequency
studied in Fig. (9a); numerically this is modelled by finite elements studying
one unit cell and applying periodic boundary conditions. As a comparison we
solve the effective equation (35) obtained by HFH to get the envelope function
$f_{0}$ and then recover the first order field using Eq. (27). As displayed on
Fig. (10a), an excellent agreement is obtained with the simulations for the
microstructured medium (solid lines) and the homogenized one (dashed lines).
Next, we consider a finite stack consisting of 20 periods of the
microstructured medium and compare it with the effective medium. The point
source is located in the center at $\xi=10+\phi/2$ with frequency
$\Omega=1.01\Omega_{0}^{(3)}=0.752$, and numerically we use Perfectly Matched
Layers [46] on either side to truncate the simulation domain and damp
propagating waves to avoid reflections at the computational domain boundary.
The dashed lines on Fig. (10b) show the field for the long-scale envelope
function $f_{0}$ which is in good agreement with the results from the finite
multilayer stack, albeit with some minor discrepancies likely due to the
finite extent of the stack and boundary effects not taken into account in our
model.

(a) Superposition of the exact dispersion diagram (plain lines) and the
asymptotic approximations obtained by HFH around the edges (dotted lines)

(b) Absolute error near $\kappa\delta=0$ between both dispersion relations in
a log-log scale (plain lines). Dotted lines are reference orders of
convergence.

Figure 9: Comparison of the exact dispersion relation and the effective one
obtained by HFH for the single eigenvalue case (39).

(a) Periodic source located at $\xi=0.8$ with frequency close to the fourth
eigenfrequency in Fig.(9) for $\kappa\delta=0$
($\Omega=1.01\Omega_{0}^{(4)}=1.027$). The effective coefficient is $T=1.336$.
(b) Finite stack of 20 periods with a source located at $\xi=10+\phi/2$ with
frequency close to the third eigenfrequency in Fig.(9) for $\kappa\delta=0$
($\Omega=1.01\Omega_{0}^{(3)}=0.752$). The dashed line shows the envelope
$f_{0}$ obtained by HFH. The effective coefficient is $T=-0.573$.

Figure 10: Comparison of the wavefields for a point source forcing. HFH
approximation is the one of the simple eigenvalue case (35).

##### Inside the Brillouin zone

Finally, we make use of the linear asymptotic approximations inside the
Brillouin zone, i.e. equation (66) of Section 2.4.4, on the same example. We
see in Figure 11a that using this approximation for only three points inside
the Brillouin zone and combining it with the quadratic asymptotic
approximations at the edges, we almost recover the entire dispersion diagram
(obtained with FEM). The effective coefficient $T_{i}$ also gives an insight
on the group velocity in Figure 11b since they are proportional.

(a) Exact dispersion diagram (left) and the linear asymptotic approximations
(right) obtained by HFH at three points inside the Brillouin zone (66) (dashed
lines). We also plot the quadratic HFH approximations around $\kappa\delta=0$
and $\kappa\delta=\pi$ (39) (plain lines)  (b) Evolution of the coefficient
$T_{i}$ (65) as a function of $\kappa\delta$ for the five modes studied in
Figure 11a.

Figure 11: Effective properties obtained by HFH inside the Brillouin zone
(Section 2.4.4).

####  2.5.3 Double eigenvalue case

We now investigate a double eigenvalue case (asymptotic approximations
developed in Section 2.4.2). Choosing $\phi=0.5$, $b=1$ and a Drude model with
no damping for $a$ with the parameters $a_{0}=a_{\infty}=1$,
$\gamma_{1,1}=\Omega_{D,1,1}=0$ and $\Omega_{1,1}=24.3347$ in (103). This
leads to a double eigenvalue $\Omega_{0}=25.1322$ at $\kappa\delta=\pi$, with
$\mathbb{a}(\Omega_{0})\simeq m^{2}$ where $m$ is an integer ($m=4$ here, cf.
[23]). We note that the value of $\Omega_{1,1}$ is close to $\Omega_{0}$ and
corresponds to a pole of $\mathbb{a}$, meaning the behaviour of the material
properties around those frequencies is highly dispersive. Even in this case,
our method recovers the expected linear asymptotics with opposite slopes
characteristic of a degenerate root, as shown in Figure 12. This is confirmed
quantitatively by the curves in Figure 13 representing the errors between the
exact dispersion relation and the linear asymptotic approximations on a log-
log scale, showing an $\mathcal{O}((\delta\tilde{\kappa})^{2})$ convergence
for both branches.

Figure 12: A double eigenvalue case. Comparison of the exact dispersion (green
solid lines) with the linear asymptotics obtained by HFH (47) (black dashed
lines). The colormap on the left panel shows the logarithm of the determinant
whose zeros give the dispersion relation, and the dotted horizontal line
indicates the position of a pole of $\mathbb{a}$ and of the determinant at
$\Omega_{1,1}$. The right panel is a zoom-magnified view close to
$\kappa\delta=\pi$. The computed value of the slope is $T_{D}=17.240$ Figure
13: Absolute error near $\kappa\delta=\pi$ between both dispersion relations
(the exact dispersion relation solved by FEM and the HFH approximation (47))
in a log-log scale (plain lines, blue for the lowest branch “-” and red for
the upper branch “+”). The dotted line is the reference order of convergence.

##  3 Two-dimensional (2D) case

We now extend the results to the 2D case; as the method is very similar to the
one used in 1D, we will highlight only the differences due to the higher
dimensions together with the final asymptotic approximations obtained.

###  3.1 Setting

We consider the Helmholtz equation in a doubly periodic structure on a square
lattice of size $h$, see Figure 14a,

$\bm{\nabla}_{\bm{X}}\cdot\left[\hat{a}_{h}(\bm{X},\omega)\bm{\nabla}_{\bm{X}}U_{h}(\bm{X})\right]+\omega^{2}\hat{b}_{h}(\bm{X},\omega)U_{h}(\bm{X})=0$
(70)

(a) 2D configuration with zoom on a unit cell

(b) The irreducible Brillouin zone, in wavenumber space

Figure 14: 2D configuration in physical space with zoom on an inclusion in the
different coordinate systems (a) and in wavenumber space (b)

The parameters are frequency-dependent inside $\hat{\mathcal{Y}}_{i}$ in the
periodic unit cell, and are simply taken to be constants outside it. As in the
1D case the methodology is developed for any frequency-dependent function and
the Drude-Lorentz model is used for numerical examples (see Appendix A for
details on this model).
As in 1D, we introduce the two-scales $\bm{x}=\bm{X}/L$ and
$\bm{\xi}=\bm{X}/h=\bm{x}/\delta$; we call $\mathcal{Y}$ the unit cell in
$\bm{\xi}$-coordinate, with $\mathcal{Y}_{i}$ the inclusion where the
parameters are frequency-dependent. Except for this geometry difference, the
non-dimensionalization step is the same as in 1D and a Bloch-Floquet analysis
similar to Section 2.3 allows us to get the 2D dispersion relation.

####  3.1.1 Ansatz

We pick a frequency-wavenumber pair
$(\Omega_{0},\bm{\kappa})\in\mathbb{C}\times[0,\pi/\delta]^{2}$ that satisfies
the dispersion relation in the irreducible Brillouin zone, see Figure 14b. The
ansatz for the non-dimensionalized field $u$ and frequency $\Omega$ (3.6),
together with the expansions for both $a$ and $b$ (21) are the same in 2D so
we get the following non-dimensionalized equation:

$\displaystyle\sum_{j\geqslant 0}$
$\displaystyle\left[a(\bm{\xi},\Omega)\left\\{\delta^{j}\Delta_{\bm{\xi}}u_{j}+\delta^{j+1}\left[\bm{\nabla}_{\bm{\xi}}\cdot\bm{\nabla}_{\bm{x}}u_{j}+\bm{\nabla}_{\bm{x}}\cdot\bm{\nabla}_{\bm{\xi}}u_{j}\right]+\delta^{j+2}\Delta_{\bm{x}}u_{j}\right\\}\right.$
(71) $\displaystyle+\sum_{\ell\geqslant
0}\delta^{\ell+j}\Omega_{\ell}^{2}b(\bm{\xi},\Omega)u_{j}]=0.$

####  3.1.2 Zeroth-order field

Collecting the terms of order $\delta^{0}$, using continuity and periodicity,
we get in $\mathcal{Y}$:

$\left\\{\begin{aligned}
&A_{0}\Delta_{\bm{\xi}}u_{0}+\Omega_{0}^{2}B_{0}u_{0}=0\text{ in
}\mathcal{Y}\\\ &u_{0}\text{ continuous at }\partial\mathcal{Y}_{i}\text{ and
}\partial\mathcal{Y}\\\ &A_{0}\bm{\nabla}_{\bm{\xi}}u_{0}\cdot\bm{n}\text{
continuous at }\partial\mathcal{Y}_{i}\text{ and }\partial\mathcal{Y}\\\
&u_{0}(\bm{x},\xi_{1}+1,\xi_{2})=\mathrm{e}^{\mathrm{i}\kappa_{1}\delta}u_{0}(\bm{x},\xi_{1},\xi_{2})\text{
and
}u_{0}(\bm{x},\xi_{1},\xi_{2}+1)=\mathrm{e}^{\mathrm{i}\kappa_{2}\delta}u_{0}(\bm{x},\xi_{1},\xi_{2})\end{aligned}\right.$
(72)

###  3.2 High-frequency homogenization for single eigenvalues at the edges of
the Brillouin zone

We start with the case of the edges of the Brillouin zone $\mathsf{\Gamma}$,
$\mathsf{X}$ and $\mathsf{M}$, for which
$\bm{\kappa}=(0,0),\,(\pi/\delta,0),\,(\pi/\delta,\pi/\delta)$, respectively.
We choose $\Omega_{0}$ which is assumed to be a simple eigenvalue associated
to the eigenfunction $\mathcal{U}_{0}(\bm{\xi},\Omega_{0})$. The zeroth-order
field then writes

$u_{0}(\bm{x},\bm{\xi})=f_{0}(\bm{x})\mathcal{U}_{0}(\bm{\xi},\Omega_{0}),$

where $f_{0}(\bm{x})$ has to be determined.

####  3.2.1 First-order field

For this single eigenvalue case, we assume that $\Omega_{1}=0$ and we are
looking for the quadratic term $\Omega_{2}$. Collecting the terms of order
$\delta$, we get in $\mathcal{Y}$:

$A_{0}\left(\Delta_{\bm{\xi}}u_{1}+2\bm{\nabla}_{\bm{\xi}}\cdot\bm{\nabla}_{\bm{x}}u_{0}\right)+\Omega_{0}^{2}B_{0}u_{1}=0$
(73)

together with:

  * •

$u_{1}(\bm{x},1,\xi_{2})=\pm u_{1}(\bm{x},0,\xi_{2})$ and
$u_{1}(\bm{x},\xi_{1},1)=\pm u_{1}(\bm{x},\xi_{1},0)$,

  * •

continuity for $u_{1}$ at $\partial\mathcal{Y}_{i}$ and $\partial\mathcal{Y}$,

  * •

continuity for
$\left[A_{0}\left(\bm{\nabla}_{\bm{x}}u_{0}+\bm{\nabla}_{\bm{\xi}}u_{1}\right)\right]\cdot\bm{n}$
at $\partial\mathcal{Y}_{i}$ and $\partial\mathcal{Y}$.

Then, we write $u_{1}$ as:

$u_{1}({\bm{x}},\bm{\xi})=f_{1}({\bm{x}})\mathcal{U}_{0}({\bm{\xi}})+\bm{\nabla}_{\bm{x}}f_{0}(\bm{x})\cdot\bm{V}(\bm{\xi})$
(74)

where $\bm{V}=(V_{1},V_{2})$ with for $j=1,2$:

$\left\\{\begin{aligned} &A_{0}\Delta
V_{j}+\Omega_{0}^{2}B_{0}V_{j}=-2A_{0}\frac{\partial\mathcal{U}_{0}}{\partial{{\xi}}_{j}}\text{
in $\mathcal{Y}$}\\\ &V_{j}(1,\xi_{2})=\pm V_{j}(0,\xi_{2})\text{ and
}V_{j}(\xi_{1},1)=\pm V_{j}(\xi_{1},0)\\\ &\text{ continuity for $V_{j}$ at
$\partial\mathcal{Y}_{i}$ and $\partial\mathcal{Y}$}\\\ &\text{ continuity for
$A_{0}\left(\bm{\nabla}V_{j}+\mathcal{U}_{0}\bm{e}_{j}\right)\cdot\bm{n}$ at
$\partial\mathcal{Y}_{i}$ and $\partial\mathcal{Y}$. }\end{aligned}\right.$
(75)

####  3.2.2 Second-order field

Collecting terms of order $\delta^{2}$, we get in $\mathcal{Y}$:

$A_{0}\left(\Delta_{\bm{x}}u_{0}+2\bm{\nabla}_{\bm{\xi}}\cdot\bm{\nabla}_{\bm{x}}u_{1}+\Delta_{\bm{\xi}}u_{2}\right)+A_{2}\Delta_{\bm{\xi}}u_{0}+\Omega_{0}^{2}(B_{0}u_{2}+B_{2}u_{0})+\Omega_{2}^{2}B_{0}u_{0}=0$
(76)

together with:

  * •

$u_{2}(\bm{x},1,\xi_{2})=\pm u_{2}(\bm{x},0,\xi_{2})$ and
$u_{2}(\bm{x},\xi_{1},1)=\pm u_{2}(\bm{x},\xi_{1},0)$

  * •

continuity for $u_{2}$ at $\partial\mathcal{Y}$ and $\partial\mathcal{Y}_{i}$

  * •

continuity for
$\left[A_{0}\left(\bm{\nabla}_{\bm{\xi}}u_{2}+\bm{\nabla}_{\bm{x}}u_{1}\right)+A_{2}\bm{\nabla}_{\bm{\xi}}u_{0}\right]\cdot\bm{n}$
at $\partial\mathcal{Y}$ and $\partial\mathcal{Y}_{i}$.

We introduce the average operator $\langle\cdot\rangle$ in 2D

$\langle f\rangle=\int_{\mathcal{Y}}f(\bm{\xi})\mathrm{d}\bm{\xi}$

and then consider the expression

$\langle{u_{2}}\times\eqref{syst_order0_2D}-u_{0}\times{\eqref{syst_order2_main_2D}}\rangle)=0.$

After integration by parts, some algebra and dividing through by $f_{0}$ we
get the final effective equation:

$T_{ij}\frac{\partial^{2}f_{0}}{\partial x_{i}\partial
x_{j}}+\Omega_{2}^{2}f_{0}=0$ (77)

with

$T_{ij}=\displaystyle\frac{\left\langle
A_{0}\left(V_{i}\frac{\partial\mathcal{U}_{0}}{\partial\xi_{j}}-\frac{\partial
V_{i}}{\partial\xi_{j}}\mathcal{U}_{0}-\mathcal{U}_{0}^{2}\delta_{ij}\right)\right\rangle}{\mathcal{S}(\Omega_{0},\mathcal{U}_{0},\mathcal{U}_{0})}$
(78)

where we defined in 2D

$\mathcal{S}(\Omega_{0},f,g)=\left\langle\left(b(\cdot,\Omega_{0})+\frac{\Omega_{0}}{2}\frac{\partial
b}{\partial\Omega}(\cdot,\Omega_{0})\right)fg-\frac{1}{2\Omega_{0}}\frac{\partial
a}{\partial\Omega}(\cdot,\Omega_{0})\bm{\nabla}f\cdot\bm{\nabla}g\right\rangle$
(79)

and where we sum over the repeated subscript indexes.
From this effective equation, we also get the quadratic term for the
dispersion relation

$\Omega_{2}^{2}=\tilde{\kappa_{i}}T_{ij}\tilde{\kappa_{j}}$ (80)

with $\tilde{\kappa}_{i}=\kappa_{i}-d_{i}$ where $d_{i}=0$ or $\pi/\delta$
depending on the high-symmetry point we choose.
The tensor $\bm{T}$ (78) encapsulates the effective properties of the periodic
structure beyond the quasi-static, classical homogenization, regime.
Typically, this tensor may have eigenvalues of markedly different magnitude or
of opposite sign, which leads to a change of character of the underlying
effective equation, from elliptic to parabolic and from elliptic to
hyperbolic, respectively. The former appears at a frequency near band edges
and the latter near a frequency at which a saddle point occurs in the
corresponding dispersion curves. This has been used notably to design
dielectric photonic crystals with spectacular directive emission in the form
of + and x wave patterns for a source placed inside in the microwave regime
[47]. The present high-frequency algorithm makes possible the extension of
such experiments to the optical wavelengths wherein the periodic assembly of
dielectric rods has a frequency dependent refractive index.

###  3.3 High-frequency homogenization for repeated eigenvalues at the edges
of the Brillouin zone

In this section, we still consider approximation around the edges of the
Brillouin zone, but for the case of repeated eigenvalues that gives rise to a
linear approximation of the dispersion diagram.

####  3.3.1 Zeroth-order field

System (72) still holds but now we assume repeated eigenvalues of multiplicity
$N$. We introduce the associated eigenfunctions $\mathcal{U}^{(j)}$
($j=1,...,N$). The solution for the leading-order problem is now

$u_{0}(\bm{x},\bm{\xi})=\sum_{j=1}^{N}f_{0}^{(j)}(\bm{x})\mathcal{U}^{(j)}_{0}(\bm{\xi}).$
(81)

We will denote $\eqref{syst_order0_2D}^{(l)}$ the system
$\eqref{syst_order0_2D}$ satisfied by
$u_{0}^{(l)}=f_{0}^{(l)}\mathcal{U}^{(l)}_{0}$ for $l\in\\{1,...,N\\}$.

####  3.3.2 First-order field

The main difference is that now $\Omega_{1}\neq 0$, therefore the system for
the first-order field is modified and collecting the terms of order $\delta$
we now get in $\mathcal{Y}$:

$A_{0}\left(\Delta_{\xi}u_{1}+2\bm{\nabla}_{\xi}\cdot\bm{\nabla}_{x}u_{0}\right)+\Omega_{0}^{2}(B_{0}u_{1}+B_{1}u_{0})+\Omega_{1}^{2}B_{0}u_{0}+A_{1}\Delta_{\xi}u_{0}=0$
(82)

together with:

  * •

$u_{1}(\bm{x},1,\xi_{2})=\pm u_{1}(\bm{x},0,\xi_{2})$ and
$u_{1}(\bm{x},\xi_{1},1)=\pm u_{1}(\bm{x},\xi_{1},0)$,

  * •

continuity for $u_{1}$ at $\partial\mathcal{Y}_{i}$ and $\partial\mathcal{Y}$,

  * •

continuity for
$\left[A_{0}\left(\bm{\nabla}_{x}u_{0}+\bm{\nabla}_{\xi}u_{1}\right)\right]\cdot\bm{n}+A_{1}\bm{\nabla}_{\xi}u_{0}\cdot\bm{n}$
at $\partial\mathcal{Y}_{i}$ and $\partial\mathcal{Y}$.

Let us pick one $l\in\\{1,...,N\\}$ and compute
$\langle\eqref{syst_order0_2D}^{(l)}\times
u_{1}-\eqref{syst_order1_main_2D_double}\times u_{0}^{(l)}\rangle$. By
integrating by parts, using different continuity conditions and dividing
through by $f_{0}^{(l)}$ we get the effective equation:

$\sum_{j\neq l}\langle
A_{0}\bm{W}_{0}^{jl}\rangle\cdot\bm{\nabla}_{\bm{x}}f_{0}^{(j)}=-\Omega_{1}^{2}\sum_{j}\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(l)},\mathcal{U}_{0}^{(j)})f_{0}^{(j)}$
(83)

with

$\bm{W}_{0}^{jl}=\mathcal{U}_{0}^{(l)}\bm{\nabla}_{\bm{\xi}}\mathcal{U}_{0}^{(j)}-\mathcal{U}_{0}^{(j)}\bm{\nabla}_{\bm{\xi}}\mathcal{U}_{0}^{(l)}.$
(84)

We set $f_{0}^{(l)}=\hat{f}_{0}^{(l)}\mathrm{exp}(i\tilde{\kappa_{j}}x_{j})$
and get the following system of equations:

$\bm{C}\hat{\bm{F}}_{0}=\bm{0}$ (85)

with the $2\times 2$ matrix $\bm{C}$ defined by

$C_{lj}=\mathrm{i}\langle
a(\cdot,\Omega_{0})\bm{W}_{0}^{jl}\rangle\cdot\tilde{\bm{\kappa}}+\Omega_{1}^{2}\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(l)},\mathcal{U}_{0}^{(j)}).$
(86)

The value of $\Omega_{1}^{2}$ is then obtained by solving
$\mathrm{det}(\bm{C})=0$.

###### Remark 1

For double eigenvalues we get the following expression for the linear term
(opposite slopes)

$\Omega_{1}^{2}=\pm\frac{\tilde{\bm{\kappa}}\cdot\langle
a(\cdot,\Omega_{0})\bm{W_{0}}^{12}\rangle}{\left(\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(1)})\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(2)},\mathcal{U}_{0}^{(2)})-(\mathcal{S}(\Omega_{0},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)}))^{2}\right)^{1/2}}.$
(87)

###  3.4 High-frequency homogenization for nearby eigenvalues at the
edges of the Brillouin zone

In this section, we again consider asymptotic approximations around the edges
of the Brillouin zone, and we assume that the eigenvalues are single but close
to each other. More precisely, we consider $N$ eigenvalues close to each other
so that the distances between them scale into the small parameter $\delta$ and
write

$\delta\alpha_{l}=(\Omega_{0}^{(l)})^{2}-(\Omega_{0}^{(1)})^{2},$ (88)

for $l\in\\{1,...,N\\}$. To take into account their competitive nature, we
assume that the leading order field is

$u_{0}(\bm{x},\bm{\xi})=\sum_{j=1}^{N}f_{0}^{(j)}(\bm{x})\mathcal{U}_{0}^{(j)}(\bm{\xi}).$
(89)

As in 1D, the ansatz is considered around the eigenvalue $\Omega_{0}^{(1)}.$
In that case, the residual term for the zeroth order equation is

$\delta\sum_{q\neq
1}\alpha_{q}u_{0}^{(q)}\left[\frac{1}{2A_{0}^{(1)}}\Omega_{0}^{(1)}B_{0}^{(1)}\frac{\partial
a}{\partial\Omega}(\cdot,\Omega_{0}^{(1)})-B_{0}^{(1)}-\frac{1}{2}\Omega_{0}^{(1)}\frac{\partial
b}{\partial\Omega}(\cdot,\Omega_{0}^{(1)})\right]$ (90)

that will in turn modify the equation for the first order field in
$\mathcal{Y}$ to be:

$\displaystyle
A_{0}^{(1)}\left(\Delta_{\xi}u_{1}+2\bm{\nabla}_{\xi}\cdot\bm{\nabla}_{x}u_{0}\right)+(\Omega_{0}^{(1)})^{2}(B_{0}^{(1)}u_{1}+B_{1}u_{0})+\Omega_{1}^{2}B_{0}^{(1)}u_{0}+A_{1}\Delta_{\xi}u_{0}$
(91) $\displaystyle+\sum_{q\neq
1}\alpha_{q}u_{0}^{(q)}\left[\frac{1}{2A_{0}^{(1)}}\Omega_{0}^{(1)}B_{0}^{(1)}\frac{\partial
a}{\partial\Omega}(\cdot,\Omega_{0}^{(1)})-B_{0}^{(1)}-\frac{1}{2}\Omega_{0}^{(1)}\frac{\partial
b}{\partial\Omega}(\cdot,\Omega_{0}^{(1)})\right]=0,$

together with $u_{1}(\bm{x},1,\xi_{2})=\pm u_{1}(\bm{x},0,\xi_{2})$ and
$u_{1}(\bm{x},\xi_{1},1)=\pm u_{1}(\bm{x},\xi_{1},0)$, continuity for $u_{1}$
at $\partial\mathcal{Y}_{i}$ and $\partial\mathcal{Y}$, and continuity for
$\left[A_{0}\left(\bm{\nabla}_{x}u_{0}+\bm{\nabla}_{\xi}u_{1}\right)\right]\cdot\bm{n}+A_{1}\bm{\nabla}_{\xi}u_{0}\cdot\bm{n}$
at $\partial\mathcal{Y}_{i}$ and $\partial\mathcal{Y}$.

Let us pick one $l\in\\{1,...,N\\}$ and consider
$\langle\eqref{syst_order0_2D}^{(l)}\times
u_{1}-\eqref{syst_order1_main_2D_nearby}\times u_{0}^{(l)}\rangle=0$. By
integrating by parts, using the different continuity conditions, dividing
through by $f_{0}^{(l)}$, and neglecting the higher-order terms we get the
effective equation

$\sum_{j=1}^{N}\langle
A_{0}\bm{W}_{\mathrm{nby}}^{jl}\rangle\cdot\bm{\nabla}_{x}f_{0}^{(j)}=\sum_{j=1}^{N}\\{-\Omega_{1}^{2}\langle\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(l)},\mathcal{U}_{0}^{(j)})\rangle+\alpha_{j}\mathcal{N}^{jl}\\}f_{0}^{(j)}$
(92)

with

$\left\\{\begin{aligned}
&\bm{W}_{\mathrm{nby}}^{jl}=\mathcal{U}_{0}^{(l)}\bm{\nabla}_{\xi}\mathcal{U}_{0}^{(j)}-\mathcal{U}_{0}^{(j)}\bm{\nabla}_{\xi}\mathcal{U}_{0}^{(l)},\\\
&\mathcal{N}^{jl}=\left\langle
b(\cdot,\Omega_{0}^{(1)})\left(1+\frac{\Omega_{0}^{(1)}}{2}\frac{\partial}{\partial\Omega}[\log(b/a)](\cdot,\Omega_{0}^{(1)})\right)\mathcal{U}_{0}^{(j)}\mathcal{U}_{0}^{(l)}\right\rangle.\end{aligned}\right.$
(93)

The linear term $\Omega_{1}$ of the dispersion relation is given by solving
$\mathrm{det}(\bm{C}_{\mathrm{nby}})=0$ where the $N\times N$ matrix
$\bm{C}_{\mathrm{nby}}$ is defined by

$(\bm{C}_{\mathrm{nby}})_{lj}=\mathrm{i}\left\langle
a(\cdot,\Omega_{0}^{(1)})\bm{W}_{\mathrm{nby}}^{jl}\right\rangle\cdot\tilde{\bm{\kappa}}+\Omega_{1}^{2}\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(l)},\mathcal{U}_{0}^{(j)})-\alpha_{j}\mathcal{N}^{jl}.$
(94)

We can notice that we recover the repeated eigenvalues case when the distances
$\alpha_{j}$ tend to 0.

###### Remark 2

In the case of two nearby eigenvalues the dispersion relation is given by

$\displaystyle\Omega_{1}^{4}\left[\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(1)})\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(2)},\mathcal{U}_{0}^{(2)})-(\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)}))^{2}\right]$
(95)
$\displaystyle+\Omega_{1}^{2}\alpha_{2}\left[-\mathcal{N}^{22}\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(1)})+\mathcal{N}^{21}\mathcal{S}(\Omega_{0}^{(1)},\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)})\right]$
$\displaystyle-\left(\langle
a(\cdot,\Omega_{0}^{(1)})\bm{W}_{\mathrm{nby}}^{21}\rangle\cdot\tilde{\bm{\kappa}}\right)^{2}-\mathrm{i}\alpha_{2}\mathcal{N}^{21}\langle
a(\cdot,\Omega_{0}^{(1)})\bm{W}_{\mathrm{nby}}^{21}\rangle\cdot\tilde{\bm{\kappa}}=0.$

###  3.5 High-frequency homogenization for simple eigenvalues inside the
Brillouin zone (no damping)

Then, when no damping is considered, we are able to get linear asymptotic
approximations near $\bm{\kappa}^{\star}$ which is not one of the high-
symmetry points. The method being very similar to the 1D case, we only give
the effective equation

$\bm{T}^{\mathrm{int}}\cdot\bm{\nabla}f_{0}-\mathrm{i}\Omega_{1}^{2}f_{0}=0$
(96)

where $\bm{T}^{\mathrm{int}}=(T_{1}^{\mathrm{int}},T_{2}^{\mathrm{int}})$ with
$T^{\mathrm{int}}_{i}$ given, for $i=1,2$ by:

$T^{\mathrm{int}}_{i}=\frac{2\langle
A_{0}\Im(\frac{\partial\mathcal{U}_{0}}{\partial\xi_{i}}\bar{\mathcal{U}}_{0})\rangle}{\mathcal{S}(\Omega_{0},\mathcal{U}_{0},\bar{\mathcal{U}_{0}})}.$
(97)

And consequently, the dispersion relation reads

$\Omega_{1}^{2}=\bm{T}^{\mathrm{int}}\cdot(\bm{\kappa}-\bm{\kappa}^{\star}).$
(98)

###  3.6 Low-frequency case

Finally, we also obtain the classical low-frequency homogenized equation. In
that case, we consider $\Omega_{0}=0$ in (72), which leads to the fact that
$\mathcal{U}_{0}$ is uniform, say $\mathcal{U}_{0}=1$. Then, we write
$u_{1}(x,\xi)=\bm{\nabla}_{\bm{x}}f_{0}(\bm{x})\cdot\bm{V}(\bm{\xi})$ where
$u_{1}$ satisfies (LABEL:system_h_2D) with $\Omega_{0}=0$ and $u_{0}=f_{0}$.
Therefore $\bm{V}$ satisfies (LABEL:system_h_2D) with $\Omega_{0}=0$ and
$\mathcal{U}_{0}=1$. Integrating (76) on a unit cell for $\Omega_{0}=0$ then
leads to the usual homogenized equation

$\left\langle
A_{0}\left(\delta_{ij}+\frac{\partial}{\partial\xi_{i}}V_{j}\right)\right\rangle\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}f_{0}(x)+\Omega_{2}^{2}\langle B_{0}\rangle f_{0}(x)=0$
(99)

together with the dispersion relation

$\Omega_{2}^{2}=\left\langle
A_{0}\left(\delta_{ij}+\frac{\partial}{\partial\xi_{i}}V_{j}\right)\right\rangle\langle
B_{0}\rangle^{-1}\kappa_{i}\kappa_{j}.$ (100)

It is well-known one can identify the homogenized matrix in (99) making use of
the acoustic band near the origin [18, 21], albeit for non-dispersive media
(see section 3.2 in [21] for a summary of results published back in 1978 in
the first edition of [18]).

###  3.7 Numerical example

We consider here a two-dimensional lattice of dispersive rods as studied by
Brûlé et al. [41]. The material parameters are the permittivity $\varepsilon$
given by a single resonance Drude model
$\varepsilon(\Omega)=1-\Omega_{p}^{2}/(\Omega(\Omega+{\rm i}\gamma))$, with
$\Omega_{p}/2\pi=1.1$, $\gamma/2\pi=0.05$, and non magnetic material with
permeability $\mu=1$. In the TM polarization (s-polarization) case this
corresponds in our notations to $\mathbb{a}=1/\mu$ and
$\mathbb{b}=\varepsilon$, whereas in the TE polarization (p-polarization) case
$\mathbb{a}=1/\varepsilon$ and $\mathbb{b}=\mu$. The dispersion diagrams are
computed by FEM from the exact dispersion function and here again will be
denoted as the exact dispersion diagrams by opposition to the asymptotic
approximations obtained by HFH.

##### Single eigenvalues at symmetry points

The first example is a square array of period $h$ of square rods of size
$L=0.806h$ made of this Drude permittivity in vacuum. Assuming first TM
polarization, we plot the dispersion along the edge of the first Brillouin
zone on Figure 15 for the first three modes. The results obtained by HFH in
Section 3.2 approximate well the dispersion behaviour locally around the
symmetry points. The spectral features showing deformed triangles in the
complex plane and including an intertwining of the second and third band is
well recovered by the HFH asymptotic approximations (see panel (c) in Figure
15).

(a) Real part (b) Imaginary part (c) Complex plane

Figure 15: Simple eigenvalue HFH approximations (80) (dashed lines) for the
first three bands of the square photonic crystal in TM polarization.

##### Nearby eigenvalues

We now make use of the linear asymptotic approximations of Section 3.4 for
single eigenvalues near edges of the Brillouin zone. Our focus is the cluster
of four eigenfrequencies along $\mathsf{\Gamma}\mathsf{X}$ close to the point
$\mathsf{\Gamma}$: the results of FEM computations are displayed on Fig. (16).
The exact dispersion curves (solid lines) are correctly approximated by the
single eigenvalue HFH model (dashed curves), and even better so by the nearby
case (dotted lines). Indeed the nearby approximation leads to a better
prediction of the local behaviour of bands, and as in the 1D numerical example
this is particularly striking for the imaginary parts. The dispersion diagram
in the complex plane and the corresponding asymptotic approximations are
reported in Fig. (17).

(a) Real part for a cluster of four eigenvalues around $\mathsf{\Gamma}$
point. (b) Imaginary part for a cluster of four eigenvalues around
$\mathsf{\Gamma}$ point.

Figure 16: Comparison of the exact dispersion relation approximated by FEM
(solid lines) and the effective one obtained by HFH for the single eigenvalue
case (80) (dashed lines) and with nearby approximations (94) (dotted lines).
Figure 17: Comparison in the complex plane of the exact dispersion relation
computed by FEM (solid lines) and the effective one obtained by HFH for the
single eigenvalue case (80) (dashed lines) and with nearby approximations (94)
(dotted lines) around $\mathsf{\Gamma}$ point. The circle markers correspond
to $\kappa_{x}\delta=0$.

##### Simple eigenvalues inside the Brillouin zone

For eigenvalues not located at symmetry points, we study the same structure
but setting $\gamma=0$. Making use of the linear approximation, (98), obtained
in Section 3.5, we recover locally the behaviour of the bands (see Figure 18).

(a) TE polarization (b) TM polarization

Figure 18: Superposition of the exact dispersion diagram (plain lines) and the
asymptotic approximations obtained by HFH (dashed lines) for points inside the
Brillouin zone (98).

##### Repeated eigenvalues

For certain choices of the dispersive behaviour and material distribution,
there exist accidental degeneracies. In this section we choose circular
inclusions with $R=0.364h$ and use the same Drude model for the permittivity
with $\Omega_{p}/2\pi=1.434$ and $\gamma=0$. The band structure in TM
polarization is represented on Figure 19, where we can see around
$\mathsf{\Gamma}$ point the coalescence of four bands around $\Omega=13.4$.
The repeated eigenvalue approximation obtained by HFH (see Section 3.3) around
this point shows two linear terms with opposite slopes $\Omega_{1}^{2}=\pm
24.64\kappa_{x}$ and two flat bands with slope close to zero that approximate
well the exact curves.

Figure 19: Comparison of the exact dispersion relation (red solid lines) and
the effective one obtained by HFH for the repeated eigenvalue case (86) (black
dashed lines) around $\mathsf{\Gamma}$ point.

##### Field approximation

HFH theory has been developed in this paper for periodic infinite systems.
However, in this paragraph, we use the asymptotic approximations obtained to
describe finite-size systems, neglecting therefore the boundary effects: we
now consider a 14 by 14 square array of rods with the same Drude permittivity
as for the approximation of Figure 15 but with a smaller damping term
$\gamma/2\pi=0.005$. The finite photonic crystal is excited by a line source
at the center with frequency close to the real part of an eigenfrequency of
the periodic system near symmetry point $\mathsf{X}$. The first example is an
array of square rods of size $L=0.806h$ in TM polarization (see Figure 20). At
the chosen frequency, the real parts of the coefficients $T_{ii}$ in the
effective tensor of the HFH theory (78) of Section 3.2 have opposite sign,
$\Re(T_{11})=-2.108$ and $\Re(T_{22})=0.876$ and are of the same order of
magnitude, leading to an effective hyperbolic behavior. The predicted
theoretical wave field distribution is shown in Figure 20b: the wave
propagation is highly directive and is aligned along the diagonals of the
system. The resulting X-shape effect comes from the superposition of two
effective media, for points $\mathsf{X}(\pi,0)$ and $\mathsf{M}(0,\pi)$ of the
Brillouin zone. The full numerical simulations shown in Figure 20a share this
same qualitative feature.
Next we study the TE polarization case for cylindrical rods of diameter
$D=0.91h$. In this case, HFH predicts a distinct anisotropy aligned along the
lattice axis, since $|\Re(T_{11})|\ll|\Re(T_{22})|$. This is indeed what we
observe on full wave simulations shown in Figure 21a with the same qualitative
agreement for the solution of the effective parabolic equation shown in Figure
21b. This directive emission due to the excitation of a surface plasmon-like
mode, where the field is mostly confined at the interface between the
dielectric background and the Drude metal, is well captured by the dispersive
HFH theory, including the field decay as a result of material losses.

(a) Photonic crystal (b) HFH

Figure 20: Field reconstruction for a finite square array of 196 square rods
in TM polarization excited by a line source at the center with frequency
$\Omega=5.496$, close to the real part of an eigenfrequency of the periodic
system $\Omega_{0}=5.442-0.008{\rm i}$ at point $\mathsf{X}$. The dynamic
effective tensor (78) is $\mathbf{T}=\operatorname{diag}(-2.108-0.018{\rm
i},\,0.876+0.003{\rm i})$, indicating hyperbolic behaviour. In both figures
the norm of the field has been normalized to its maximum value.

(a) Photonic crystal (b) HFH

Figure 21: Field reconstruction for a finite square array of 196 circular rods
in TE polarization excited by a line source at the center with frequency
$\Omega=3.509$, close to the real part of an eigenfrequency of the periodic
system $\Omega_{0}=3.506-0.010{\rm i}$ at point $\mathsf{X}$. The dynamic
effective tensor (78) is $\mathbf{T}=\operatorname{diag}(0.005+3\times
10^{-5}{\rm i},\,0.950-0.012{\rm i})$, indicating parabolic behaviour. In both
figures the norm of the field has been normalized to its maximum value.

##  4 Conclusion

In this paper, we have extended the technique of high-frequency homogenization
to dispersive periodic media, within which the physical properties do depend
on the frequency. The work has been performed in both 1D and 2D and we
considered different cases depending on the nature of the point around which
we want to build an effective approximation. Near the edges of the Brillouin
zone, we performed high-frequency homogenization for a frequency which is a
simple eigenvalue, or a repeated one, and we also considered the case where
several single eigenvalues are close to each other. Far from the edges, we
developed an approximation in the case where no damping is considered.
In each of these cases, we were able to develop an approximation for both the
dispersion diagram and the envelope function which defines the wavefield at
the zeroth order. These asymptotic approximations come with an effective
parameter or tensor that encapsulates the dispersive properties of the
considered material. The results have been validated using comparisons with
Finite Element Simulations for different configurations. We also discussed the
interpretation of the effective parameter with respect to the nature of the
wavefield.
Potential extensions of this work include pushing the asymptotics presented in
this paper to higher orders, or extending it to the case of waves in other
periodic dispersive media: for example the full vector Maxwell’s equations
(such as in photonic crystal fibres within which TE and TM waves are usually
fully coupled in oblique incidence), the Navier equations for fully coupled
in-plane pressure and shear waves in phononic crystals, the Kirchhoff-Love
equations for flexural waves in thin plates, or arrays of resonators (such as
Helmholtz resonators, high-contrasted inclusions [48, 10], or bubbles hosting
Minnaert resonant frequencies [49]).

##  Appendix A Expressions of $a$ and $b$ for the Lorentz model

##### Physical parameters

In the case of a Lorentz model, the physical parameters (2) and (3) are

$\hat{a}_{h}(X,\omega)=\left\\{\begin{aligned} &{a}_{0}&\text{ for
}X\in(0,\phi h)\\\ &\hat{\mathbb{a}}(\omega)={a}_{\infty}\left(1-\sum_{p\geq
0}\frac{\omega_{p,1}^{2}}{\omega(\omega+\mathrm{i}\hat{\gamma}_{p,1})-\omega_{D,p,1}^{2}}\right)^{-1}&\text{
for }X\in(\phi h,h)\end{aligned}\right.$ (101)

and

$\hat{b}_{h}(X,\omega)=\left\\{\begin{aligned} &{b}_{0}&\text{ for
}X\in(0,\phi h)\\\ &\hat{\mathbb{b}}(\omega)={b}_{\infty}\left(1-\sum_{p\geq
0}\frac{\omega_{p,2}^{2}}{\omega(\omega+\mathrm{i}\hat{\gamma}_{p,2})-\omega_{D,p,2}^{2}}\right)&\text{
for }X\in(\phi h,h)\end{aligned}\right.$ (102)

respectively, where $\omega_{p,i}$ stands for a plasmon frequency,
$\hat{\gamma}_{p_{i}}$ for a damping coefficient, and $\omega_{D,p,i}$ for a
Lorentz resonant frequency. When all the $\omega_{D,p,i}$ are zero, the
Lorentz model becomes known as the Drude model.

##### Adimensionalized parameters

The adimensionalized physical parameters (6) and (7) for this model then read

$a(\xi,\Omega)=\left\\{\begin{aligned} &1&\text{ for }\xi\in(0,\phi)\\\
&\mathbb{a}(\Omega)=\frac{{a}_{\infty}}{a_{0}}\left(1-\sum_{p\geq
0}\frac{\Omega_{p,1}^{2}}{\Omega(\Omega+\mathrm{i}{\gamma}_{p,1})-\Omega_{D,p,1}^{2}}\right)^{-1}&\text{
for }\xi\in(\phi,1)\end{aligned}\right.$ (103)

and

$b(\xi,\Omega)=\left\\{\begin{aligned} &1&\text{ for }\xi\in(0,\phi)\\\
&\mathbb{b}(\Omega)=\frac{{b}_{\infty}}{b_{0}}\left(1-\sum_{p\geq
0}\frac{\Omega_{p,2}^{2}}{\Omega(\Omega+\mathrm{i}{\gamma}_{p,2})-\Omega_{D,p,2}^{2}}\right)&\text{
for }\xi\in(\phi,1)\end{aligned}\right.$ (104)

respectively, with for $i=1,\,2$:

$\Omega_{p,i}=\frac{\omega_{p,i}h}{c_{0}}\text{,
}\Omega_{D,p,i}=\frac{\omega_{D,p,i}h}{c_{0}}\text{ and
}\gamma_{p,i}=\frac{\hat{\gamma}_{p,i}h}{c_{0}}.$

##  Appendix B Finite element formulation

Using a Lorentz model, we write
$\mathbb{a}(\Omega)=\frac{N_{a}(\Omega)}{D_{a}(\Omega)}$ and
$\mathbb{b}(\Omega)=\frac{N_{b}(\Omega)}{D_{b}(\Omega)}$, where $N_{a}$,
$D_{a}$, $N_{b}$, $D_{b}$ are polynomials of $\Omega$. We have to solve the
following eigenproblem

$\nabla\cdotp\left[a(\xi,\Omega)\nabla
u(\xi)\right]+\Omega^{2}b(\xi,\Omega)u(\xi)=0.$ (105)

The weak formulation of the problem is derived by multiplying Eq. (105) by the
complex conjugate of a test function $v$ and integrating the first term by
part on the unit cell $\mathcal{Y}$:

$-\int_{\mathcal{Y}}a(\Omega)\nabla u(\xi)\cdotp\nabla v^{*}(\xi){\rm
d}\xi+\Omega^{2}\int_{\mathcal{Y}}b(\Omega)u(\xi)v^{*}(\xi){\rm d}\xi=0$ (106)

where the boundary term vanishes because of the quasi-periodic boundary
conditions. We define for:

$\mathbb{A}_{0}(u,v)=-\int_{\mathcal{Y}\setminus\mathcal{Y}_{i}}\nabla
u(\xi)\cdotp\nabla v^{*}(\xi){\rm
d}\xi,\qquad\mathbb{A}_{1}(u,v)=-\int_{\mathcal{Y}_{i}}\nabla
u(\xi)\cdotp\nabla v^{*}(\xi){\rm d}\xi$

and

$\mathbb{B}_{0}(u,v)=\int_{\mathcal{Y}\setminus\mathcal{Y}_{i}}u(\xi)v^{*}(\xi){\rm
d}\xi,\qquad\mathbb{B}_{1}(u,v)=\int_{\mathcal{Y}_{i}}u(\xi)v^{*}(\xi){\rm
d}\xi.$

Plugging the expression for $a$ and $b$ in (106) and rearranging we get:

$\displaystyle\mathcal{T}(\Omega,u,v)=$ $\displaystyle
D_{a}(\Omega)D_{b}(\Omega)\left[\mathbb{A}_{0}(u,v)+\Omega^{2}\mathbb{B}_{0}(u,v)\right]$
$\displaystyle+N_{a}(\Omega)D_{b}(\Omega)\mathbb{A}_{1}(u,v)$
$\displaystyle+\Omega^{2}D_{a}(\Omega)N_{b}(\Omega)\mathbb{B}_{1}(u,v)=0$

which is a polynomial eigenvalue problem solved using the open source FEniCS
finite element library [50] interfaced with the SLEPc eigensolver [51, 52].

## Acknowledgments

MT and RA would like to thank the Isaac Newton Institute for Mathematical
Sciences, Cambridge, for support and hospitality during the programme
Mathematical theory and applications of multiple wave scattering where work on
this paper was undertaken. This work was supported by EPSRC grant no
EP/R014604/. BV and RVC are supported by the H2020 FET-proactive Metamaterial
Enabled Vibration Energy Harvesting (MetaVEH) project under Grant Agreement
No. 952039. SG and RVC were funded by UK Research and Innovation (UKRI) under
the UK government’s Horizon Europe funding guarantee [grant number 10033143].
The authors also thank Ping Sheng for fruitful discussions about the
accumulation points.

## References

  * [1] Vasily Klimov.  Nanoplasmonics.  Taylor Francis, New York, 2014.
  * [2] Stefan A. Maier.  Plasmonics: Fundamentals and Applications.  Springer Netherlands, 2007.
  * [3] David R Smith, John B Pendry, and Mike CK Wiltshire.  Metamaterials and negative refractive index.  science, 305(5685):788–792, 2004.
  * [4] Jensen Li, Lei Zhou, C. T. Chan, and P. Sheng.  Photonic band gap from a stack of positive and negative index materials.  Physical Review Letters, 90(8):083901, feb 2003.
  * [5] P Y Chen, C G Poulton, A A Asatryan, M J Steel, L C Botten, C Martijn de Sterke, and R C McPhedran.  Folded bands in metamaterial photonic crystals.  New Journal of Physics, 13(5):053007, may 2011.
  * [6] J. C. Maxwell Garnet.  Colours in metal glasses and in metallic films.  Phil. Trans. R. Soc. Lond., 203:385–420, 1904.
  * [7] A F J Levi.  The Drude model.  In Essential Classical Mechanics for Device Physics, 2053-2571, pages 6–1 to 6–20. Morgan & Claypool Publishers, 2016.
  * [8] J.-L. Auriault and G. Bonnet.  Dynamique des composites élastiques périodiques.  Arch Mech., 37(4-5):269–284, 1985.
  * [9] J.-L. Auriault and C. Boutin.  Long wavelength inner-resonance cut-off frequencies in elastic composite materials.  International Journal of Solids and Structures, 49(23-24):3269–3281, nov 2012.
  * [10] Guy Bouchitté, Christophe Bourel, and Didier Felbacq.  Homogenization near resonances and artificial magnetism in 3d dielectric metamaterials.
  * [11] Didier Felbacq and Guy Bouchitté.  Theory of mesoscopic magnetism in photonic crystals.  Physical Review Letters, 94(18), may 2005.
  * [12] Nicholas Fang, Dongjuan Xi, Jianyi Xu, Muralidhar Ambati, Werayut Srituravanich, Cheng Sun, and Xiang Zhang.  Ultrasonic metamaterials with negative modulus.  Nature Materials, 5(6):452–456, apr 2006.
  * [13] Michael R. Haberman and Matthew D. Guild.  Acoustic metamaterials.  Physics Today, 69(6):42–48, jun 2016.
  * [14] Alexander Moroz.  Three-dimensional complete photonic-band-gap structures in the visible.  Physical Review Letters, 83(25):5274, 1999.
  * [15] Alexander Moroz.  Metallo-dielectric diamond and zinc-blende photonic crystals.  Physical Review B, 66(11):115109, 2002.
  * [16] L. Brillouin.  Wave Propagation in Periodic Structures: Electric Filters and Crystal Lattices.  Dover books on engineering and engineering physics. McGraw-Hill Book Company, Incorporated, 1946.
  * [17] Calvin H. Wilcox.  Theory of bloch waves.  Journal d'Analyse Mathématique, 33(1):146–167, dec 1978.
  * [18] Alain Bensoussan, Jacques-Louis Lions, and George Papanicolaou.  Asymptotic analysis for periodic structures, volume 374.  American Mathematical Soc., 2011.
  * [19] N. Bakhvalov and G. Panasenko.  Homogenisation: Averaging Processes in Periodic Media.  Springer Netherlands, 1989.
  * [20] D. Cioranescu and P. Donato.  An Introduction to Homogenization.  An Introduction to Homogenization. Oxford University Press, 1999.
  * [21] Carlos Conca and Muthusamy Vanninathan.  Homogenization of periodic structures via bloch decomposition.  SIAM Journal on Applied Mathematics, 57(6):1639–1659, dec 1997\.
  * [22] Nicolas Lebbe, Agnès Maurel, and Kim Pham.  Homogenized transition conditions for plasmonic metasurfaces.  Physical Review B, 107(8):085124, feb 2023.
  * [23] R. V. Craster, J. Kaplunov, and A. V. Pichugin.  High-frequency homogenization for periodic media.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 466(2120):2341–2362, August 2010.
  * [24] Bojan B. Guzina, Shixu Meng, and Othman Oudghiri-Idrissi.  A rational framework for dynamic homogenization at finite wavelengths and frequencies.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 475(2223):20180547, March 2019.
  * [25] Shixu Meng, Othman Oudghiri-Idrissi, and Bojan B. Guzina.  A convergent low-wavenumber, high-frequency homogenization of the wave equation in periodic media with a source term.
  * [26] Davit Harutyunyan, Graeme W. Milton, and Richard V. Craster.  High-frequency homogenization for travelling waves in periodic media.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(2191):20160066, jul 2016.
  * [27] R. V. Craster, J. Kaplunov, and J. Postnova.  High-frequency asymptotics, homogenisation and localisation for lattices.  The Quarterly Journal of Mechanics and Applied Mathematics, 63(4):497–519, jul 2010.
  * [28] D. J. Colquitt, R. V. Craster, and M. Makwana.  High frequency homogenisation for elastic lattices.  The Quarterly Journal of Mechanics and Applied Mathematics, 68(2):203–230, mar 2015.
  * [29] E. Nolde, R.V. Craster, and J. Kaplunov.  High frequency homogenization for structural mechanics.  Journal of the Mechanics and Physics of Solids, 59(3):651–671, mar 2011.
  * [30] Richard V. Craster, Julius Kaplunov, Evgeniya Nolde, and Sebastien Guenneau.  High-frequency homogenization for checkerboard structures: defect modes, ultrarefraction, and all-angle negative refraction.  Journal of the Optical Society of America A, 28(6):1032, may 2011\.
  * [31] T. Antonakakis and R. V. Craster.  High-frequency asymptotics for microstructured thin elastic plates and platonics.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 468(2141):1408–1427, feb 2012.
  * [32] T. Antonakakis, R.V. Craster, and S. Guenneau.  Homogenisation for elastic photonic crystals and dynamic anisotropy.  Journal of the Mechanics and Physics of Solids, 71:84–96, nov 2014\.
  * [33] Claude Boutin, Antoine Rallu, and Stephane Hans.  Large scale modulation of high frequency waves in periodic elastic composites.  Journal of the Mechanics and Physics of Solids, 70:362–381, oct 2014.
  * [34] Antoine Rallu, Stéphane Hans, and Claude Boutin.  Asymptotic analysis of high-frequency modulation in periodic systems. analytical study of discrete and continuous structures.  Journal of the Mechanics and Physics of Solids, 117:123–156, aug 2018.
  * [35] Raphaël C. Assier, Marie Touboul, Bruno Lombard, and Cédric Bellis.  High-frequency homogenization in periodic media with imperfect interfaces.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476(2244):20200402, December 2020.
  * [36] Bojan B. Guzina and Marc Bonnet.  Effective wave motion in periodic discontinua near spectral singularities at finite frequencies and wavenumbers.  Wave Motion, 103:102729, jun 2021.
  * [37] M Sh Birman and TA Suslina.  Homogenization of a multidimensional periodic elliptic operator in a neighborhood of the edge of an internal gap.  Journal of Mathematical Sciences, 136:3682–3690, 2006.
  * [38] Stefan Güttel and Françoise Tisseur.  The nonlinear eigenvalue problem.  Acta Numerica, 26:1–94, 2017.
  * [39] Jean-Michel Combes, Boris Gralak, and Adriaan Tip.  Spectral properties of absorptive photonic crystals.  In Waves in Periodic and Random Media, page 1. Contemporary Mathematics 339, 2002.
  * [40] Christian Engström and Markus Richter.  On the spectrum of an operator pencil with applications to wave propagation in periodic and frequency dependent materials.  SIAM Journal on Applied Mathematics, 70(1):231–247, jan 2009\.
  * [41] Yoann Brûlé, Boris Gralak, and Guillaume Demésy.  Calculation and analysis of the complex band structure of dispersive and dissipative two-dimensional photonic crystals.  JOSA B, 33(4):691–702, April 2016.
  * [42] Masud Mansuripur, Miroslav Kolesik, and Per Jakobsen.  Leaky modes of solid dielectric spheres.  Physical Review A, 96(1):013846, jul 2017.
  * [43] P. Jakobsen, M. Mansuripur, and M. Kolesik.  Leaky-mode expansion of the electromagnetic field inside dispersive spherical cavity.  Journal of Mathematical Physics, 59(3), mar 2018.
  * [44] Mondher Besbes and Christophe Sauvan.  Role of static modes in quasinormal modes expansions: When and how to take them into account?  Mathematics, 10(19):3542, sep 2022.
  * [45] Carlos J. Zapata-Rodríguez, David Pastor, Luis E. Martínez, María T. Caballero, and Juan J. Miret.  Single-polarization double refraction in plasmonic crystals: Considerations on energy flow.  Applied Mechanics and Materials, 472:729–733, jan 2014.
  * [46] Jean-Pierre Berenger.  A perfectly matched layer for the absorption of electromagnetic waves.  Journal of computational physics, 114(2):185–200, 1994.
  * [47] Lauris Ceresoli, Redha Abdeddaim, Tryfon Antonakakis, Ben Maling, Mohammed Chmiaa, Pierre Sabouroux, Gérard Tayeb, Stefan Enoch, Richard V Craster, and Sébastien Guenneau.  Dynamic effective anisotropy: Asymptotics, simulations, and microwave experiments with dielectric fibers.  Physical Review B, 92(17):174307, 2015.
  * [48] Marc Briane and Muthusamy Vanninathan.  First bloch eigenvalue in high contrast media.  Journal of Mathematical Physics, 55(1), 2014.
  * [49] Habib Ammari, Hyundae Lee, and Hai Zhang.  Bloch waves in bubbly crystal near the first band gap: a high-frequency homogenization approach.  SIAM Journal on Mathematical Analysis, 51(1):45–59, 2019.
  * [50] Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells.  The FEniCS Project Version 1.5.  Archive of Numerical Software, 3(100), December 2015.
  * [51] Vicente Hernandez, Jose E. Roman, and Vicente Vidal.  SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems.  ACM Transactions on Mathematical Software, 31(3):351–362, September 2005.
  * [52] Carmen Campos and Jose E. Roman.  Parallel Krylov Solvers for the Polynomial Eigenvalue Problem in SLEPc.  SIAM Journal on Scientific Computing, 38(5):S385–S411, January 2016\.