%~Mouliné par MaN_auto v.0.40.4 (569efd63) 2026-03-25 11:13:35 \documentclass[CRMECA,Unicode,XML,screen,biblatex,published]{cedram} \addbibresource{CRMECA_Lombard_20260076.bib} \usepackage{bm} \usepackage[labelformat=simple]{subcaption} \renewcommand\thesubfigure{(\alph{subfigure})} \newcommand{\clO}{\mathcal{O}} \newcommand{\clI}{\mathcal{I}} \newcommand{\clD}{\mathcal{D}} \newcommand{\clU}{\mathcal{U}} \newcommand{\clB}{\mathcal{B}} \newcommand{\clV}{\mathcal{V}} \newcommand{\clW}{\mathcal{W}} \newcommand{\clF}{\mathcal{F}} \newcommand{\clA}{\mathcal{A}} \newcommand{\clC}{\mathcal{C}} \newcommand{\bbR}{\mathbb{R}} \newcommand{\bbZ}{\mathbb{Z}} \newcommand{\bbN}{\mathbb{N}} \newcommand{\bbM}{\mathbb{M}} \newcommand{\rmd}{\mathrm{d}} \newcommand{\rmi}{\mathrm{i}} \newcommand{\rme}{\mathrm{e}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \graphicspath{{./figures/}} \newcommand*{\mk}{\mkern -1mu} \newcommand*{\Mk}{\mkern -2mu} \newcommand*{\mK}{\mkern 1mu} \newcommand*{\MK}{\mkern 2mu} %\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red} \newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}} \newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel} \newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel} \newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel} \newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %~ This fixes spacing issues around \left( and \right) and other similar symbols \let\originalleft\left \let\originalright\right \renewcommand{\left}{\mathopen{}\mathclose\bgroup\originalleft} \renewcommand{\right}{\aftergroup\egroup\originalright} \newenvironment{alcases}{\left\{\begingroup\arraycolsep=0pt\renewcommand{\arraystretch}{1.2}\begin{array}{r@{\;}l@{\quad}l}}{\end{array}\endgroup\right.} \let\oldmapsto\mapsto \renewcommand*{\mapsto}{\mathchoice{\longmapsto}{\oldmapsto}{\oldmapsto}{\oldmapsto}} \renewcommand*{\to}{\mathchoice{\longrightarrow}{\rightarrow}{\rightarrow}{\rightarrow}} \let\oldleftrightarrow\leftrightarrow \renewcommand*{\leftrightarrow}{\mathchoice{\longleftrightarrow}{\oldleftrightarrow}{\oldleftrightarrow}{\oldleftrightarrow}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \title{Predicting topologically protected interface state with high-frequency homogenization} \alttitle{Calcul par homogénéisation haute-fréquence d'un mode d'interface protégé topologiquement} \keywords{ \kwd{periodic media} \kwd{high-frequency homogenization} \kwd{topologically protected interface modes} \kwd{surface impedance} \kwd{Floquet--Bloch theory}} \altkeywords{ \kwd{milieux périodiques} \kwd{homogénéisation haute-frequence} \kwd{modes d'interface topologiquement protégés} \kwd{impédance de surface} \kwd{théorie de Floquet--Bloch}} \author{\firstname{Marie} \lastname{Touboul}\CDRorcid{0000-0003-4052-0994}\IsCorresp} \address{POEMS, CNRS, INRIA, Institut Polytechnique de Paris, 91120, Palaiseau, France} \email[M.~Touboul]{marie.touboul@ensta.fr} \author{\firstname{Bruno} \lastname{Lombard}\CDRorcid{0000-0002-5648-2737}} \address{Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR 7031, Marseille, France} \email[B.~Lombard]{lombard@lma.cnrs-mrs.fr} \author{\firstname{Antonin} \lastname{Coutant}\CDRorcid{0000-0002-5526-9184}} \addressSameAs{1}{Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR 7031, Marseille, France} \email[A.~Coutant]{coutant@lma.cnrs-mrs.fr} \begin{abstract} When two semi-infinite periodic media are joined together, a localized interface mode may exist, whose frequency belongs to their common band gap. Moreover, if certain spatial symmetries are satisfied, this mode is topologically protected and thus is robust to defects. A method has recently been proposed to identify the existence and the frequency of this mode, based on the computation of surface impedances at all the frequencies in the gap. In this work, we approximate the surface impedances thanks to high-frequency effective models, and therefore get a prediction of topologically protected interface states while only computing the solution of an eigenvalue problem at the edges of the bandgaps. We also show that the nearby eigenvalues high-frequency effective models give rise to a better approximation of the surface impedance. \end{abstract} \begin{altabstract} Quand deux milieux périodiques semi-infinis sont accolés, un mode localisé à l'interface peut exister. Sa fréquence appartient à l'intersection des bandes de fréquences interdites de chacun des deux milieux. De plus, si certaines symétries matérielles sont satisfaites, alors ce mode est protégé topologiquement et est donc robuste aux imperfections. Une méthode a récemment été proposée pour identifier l'existence et la fréquence de ce mode, basée sur le calcul des impédances de surface des deux milieux, et ce à toutes les fréquences à l'intérieur des bandes interdites. Ici, nous approchons les impédances de surface par des modèles effectifs déduits de l'homogénéisation haute-fréquence. La détermination du mode d'interface protégé topologiquement revient alors simplement à calculer un seul problème aux valeurs propres aux frontières des bandes interdites. On montre aussi que le modèle effectif reposant sur l'hypothèse de valeurs propres voisines donne une meilleure approximation des impédances de surface. \end{altabstract} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \dateposted{2026-04-08} \begin{document} %\input{CR-pagedemetas} \maketitle \section{Introduction}\label{SecIntro} It is well known that periodic media present some frequency band gaps, in which wave propagation is forbidden. When two semi-infinite periodic media, named phononic or photonic crystals, are joined together, an interface mode may exist in the band gap common to both media. A~central question is to study not only its existence but also its stability i.e. the robustness against defects~\cite{Ma19,Ozawa19}. If certain spatial symmetries are satisfied, appropriate topological invariants can be obtained from the bulk properties of the crystals~\cite{Asboth16,ProdanShultzBaldes}. The \emph{bulk-edge correspondence} then allows one to relate these invariants with the presence of interface localized modes, which are topologically protected~\cite{Delplace11}. Despite the broad applicability of the bulk-edge correspondence, its application\linebreak requires case-by-case analyses. One approach to prove the correspondence can rely on calculations of the surface impedances of the crystals in contact~\cite{Xiao14}. This approach, which was originally limited to {bilayered} configurations, has been recently generalized to any class of centro-symmetric crystals~\cite{Coutant2024,Thiang2023}. Applications to PCs with dispersion~\cite{Alexopoulos2025a} or damping~\cite{Alexopoulos2025b} have also been~investigated. In any case, the frequency of the topologically-protected mode can then be evaluated very accurately, at the price of force-brute calculations of the Bloch dispersion diagrams. On the other hand, high-frequency homogenization has been developed to provide an effective model around a given point of the dispersion diagram~\cite{Craster2010}. The method relies on identifying this point and its associated eigenfunction. This eigenfunction accounts for the fast-scale variations which are predominant at high-frequencies, a particular case being standing waves at the corners of the Brillouin zone. As one goes further away from this point, this eigenfunction is modulated by a long-scale modulation function. Through high-frequency homogenization (HFH), based on two-scale asymptotic procedures, one gets an effective ordinary differential equation for this modulation function together with an approximation of the dispersion\linebreak relation as we go away from the known point. All the information is then encapsulated in an effe\-ctive parameter or tensor which appears in the ODE for the modulation function. The idea of\linebreak using Bloch expansion to analyze asymptotic behavior dates back to~\cite{Bensoussan1978}, and the HFH methodology has been successfully applied in various contexts: discrete lattice media~\cite{Craster2010a,Colquitt2015}, frame structures~\cite{Nolde2011}, optics~\cite{Craster2011}, elastic plates~\cite{Antonakakis2012}, full vector wave systems~\cite{Antonakakis2014}, elastic composites~\cite{Boutin2014}, reticulated structures~\cite{Rallu2018}, imperfect interfaces~\cite{Assier2020,Guzina2021}, dispersive media~\cite{Touboul2024a}. Other works concerned the derivation of the process in the time domain~\cite{Harutyunyan2016}, higher-order asymptotic analysis and consideration of branches close to one another~\cite{Guzina2019}, and inclusion of a source term~\cite{Meng2021}. The purpose of this work is to bridge the gap between the study of topologically protected states with surface impedance on the one hand and high-frequency homogenization on the other. More precisely, our aim is to use directly the effective models given by high-frequency homogenization at the borders of the gap to approximate the surface impedance and therefore the frequency at which a protected state would occur within a gap. The paper is organized as follows. Firstly, in Section~\ref{SecContinu}, we recall the setting: two-semi infinite Phononic Crystal (PCs) in contact, and how the surface impedance is related to the existence of an interface localized mode at their interface. Secondly, in Section~\ref{Sec:HFH_single} the classical results of high-frequency homogenization around the edges of the dispersion diagram are recalled in the simple case of non-dispersive bulk properties with perfect contact at the interfaces. However, the methodology presented in this paper can easily be extended to more exotic configurations where HFH has already been developed. The HFH effective model is then used to approximate the surface impedance and therefore to predict a topologically-protected mode. A numerical example is given to illustrate the accuracy of this approximation. Thirdly, in Section~\ref{Sec:HFH_nearby}, we deal with another sort of HFH effective model by considering the fact that the edges of the band gap are nearby eigenvalues. Another approximation of the surface impedance is therefore derived based on this second effective model which is compared to the previous one through a numerical example, and is shown to be much more accurate. {All technical details concerning HFH are referred to in the Appendices.} \section{Physical modeling and surface impedance}\label{SecContinu} \subsection{The periodic setting} We consider the 1D Helmholtz equation, which describes linear elastic wave propagation or acoustics in mechanics, or TE or TM polarization in electromagnetism. We adopt here the notations of 1D linear elasticity to follow~\cite{Coutant2024}. The displacement field $u_h$ at a given frequency~$\omega$ therefore satisfies: \begin{equation}\label{Helmholtz} \frac{\rmd }{\rmd x}\left(E_h(x)\frac{\rmd u_h}{\rmd x}\right)+\rho_h(x)\omega^2 u_h =0, \end{equation} with $x$ being the space coordinate and $E_h$ and $\rho_h$ the Young's modulus and the mass density, respectively. The physical parameters $\rho_h$ and $E_h$ are assumed to be real and $h$-periodic i.e. \begin{equation}\label{rho_and_mu} \rho_h(x) = \rho\left(\frac{x}{h}\right) \quad \text{and} \quad E_h(x) = E\left(\frac{x}{h}\right), \end{equation} where $\rho$ and $E$ are 1-periodic functions. We choose the edges of the periodic cells to be located at $nh$ with $n\in\bbZ$ and for any 1-periodic function $g$ we introduce its mean value over a periodic cell \begin{equation}\label{eq:mean_g} \langle g\rangle = \int_0^1 g(\xi)\rmd \xi. \end{equation} The first Brillouin zone is denoted by $\clB=[\frac{-\pi}{h},\frac{\pi}{h}]$. Given a Bloch wavenumber $k\in \clB$, the squared angular frequency $\lambda=\omega^2$ is given by the eigenvalue problem: \begin{equation}\label{eq:eigenvalue} \begin{cases} \frac{\rmd }{\rmd x}\left(E\left(\frac{x}{h}\right)\frac{\rmd u_h}{\rmd x}\right)+\rho\left(\frac{x}{h}\right)\omega^2 u_h =0 \quad\text{on $(0,h)$},\\ u_h(h^+) = e^{\rmi kh}u_h(0^+) \quad\text{and}\quad u'_h(h^+) = e^{\rmi kh}u'_h(0^+). \end{cases} \end{equation} This eigenvalue problem has an infinite set of possibly repeated real positive eigenvalues $\lambda_n(k)=\omega_n^2(k)$ ordered by increasing value. For each integer $n$, the $n^{\rm th}$ branch $k\mapsto\lambda_n(k)$, is continuous with respect to $k\in\clB$ and one notes \begin{equation} \lambda_n^- =\min_{k\in\clB} \lambda_n(k) \quad\text{and}\quad \lambda_n^+= \max_{k\in\clB} \lambda_n(k). \end{equation} Therefore the entire range of possible squared angular frequencies is given by \begin{equation} \bigcup_{n\geq 1}\left[\lambda_n^-,\lambda_n^+\right]. \end{equation} If $\lambda_{n+1}^->\lambda_n^+$, there is a frequency gap $[\lambda_n^+,\lambda_{n+1}^-]$. {The eigenfunctions $u_h$ in~\eqref{eq:eigenvalue} are referred to as Bloch modes.} \subsection{Two semi-infinite PCs and surface impedance}\label{Sec:Results_Imp} We consider two semi-infinite PCs in contact at $x=0$. The interface at $x=0$ is at the boundary of the unit cells for both crystals. Both the left crystal ($x<0$) and right crystal ($x>0$) present a band structure as presented in the previous section. Let us assume that there is a non-empty overlap $\clI$ of the gaps for the left and the right crystal. For each crystal, and for a given frequency $\omega$ in the common bandgap $\clI$, two modes are solution of~\eqref{Helmholtz}: one increasing exponentially and one decreasing exponentially as $x\to \pm\infty$. We select $v(\cdot,\omega)$ the unique solution decreasing when $x\to -\infty$ (resp. $x\to +\infty$) for the left (resp. right) crystal. { An interface localized solution must be proportional with these decreasing evanescent modes on both sides of $x=0$. Existence of an interface mode is thus equivalent to the existence of non-zero $(\alpha,\beta)$ so that the solution $u_h$ reads \begin{equation} u_h(x) = \begin{cases} \alpha v(x,\omega) &\text{for } x<0 \\ \beta v(x,\omega) &\text{for } x>0. \\ \end{cases} \end{equation} Continuity of the displacement and stress leads to \begin{equation} \begin{pmatrix} v(0^-,\omega) & v(0^+,\omega)\\ E(0^-)\,v'(0^-,\omega) & E(0^+)\,v'(0^+,\omega) \end{pmatrix} \begin{pmatrix} \alpha\\ -\beta \end{pmatrix} = \begin{pmatrix} 0\\ 0 \end{pmatrix}. \end{equation} Non-trivial coefficients exist iff the determinant of the matrix is zero: \begin{equation}\label{Det-Impedance} v(0^-,\omega)\,E(0^+)\,v'(0^+,\omega)-v(0^+,\omega)\,E(0^-)\,v'(0^-,\omega)=0. \end{equation} We introduce the surface impedances~\cite{Xiao2014} for the left and right crystals: \begin{equation}\label{eq:def_impedance} Z_L(\omega) = -\frac{v(0^-,\omega)}{E(0^-)v'(0^-,\omega)} \quad \text{and}\quad Z_R(\omega) = +\frac{v(0^+,\omega)}{E(0^+)v'(0^+,\omega)}, \end{equation} where $L$ (resp. $R$ index) denotes the left (resp. right) semi-infinite PC. From~\eqref{Det-Impedance}--\eqref{eq:def_impedance}, one obtains the necessary and sufficient condition \begin{equation}\label{sum_impedances} Z_L(\omega)+Z_R(\omega) = 0 \end{equation} for the existence of an interface localized mode at frequency $\omega\in\clI.$ } In~\cite{Coutant2024}, this notion of surface impedance is used for PCs with mirror symmetry to prove the existence and uniqueness of topologically protected interface modes if there is a symmetry inversion of the Bloch modes at gap edges between the left and right crystal. More precisely the proof is based on the three following lemmas: \begin{lemm}\label{lemma:monotony} In a given gap, the surface impedance decreases with frequency. \end{lemm} From now on, we also assume that the PCs are mirror symmetric, or equivalently reflection symmetric with respect to the centre of each cell: on $[0,h]$, one has $\rho(x)=\rho(h-x)$ and $E(x)=E(h-x)$. \begin{lemm}\label{remark_symmetry} For PCs with mirror symmetry, let us consider a Bloch mode $u_h^{(n)}$ associated to the non-degenerated $n^{\rm th}$ eigenvalue at an edge of the gap ($k=0$ or $\frac{\pi}{h}$). Then: \begin{itemize} \item either $u_h^{(n)}$ is symmetric and $\bigl(u_h^{(n)}\bigr)'$ is antisymmetric, leading to $\bigl(u_h^{(n)}\bigr)'(0)=0$ for $k=0$ and $u_h^{(n)}(0)=0$ for $k=\frac{\pi}{h}$, \item or $u_h^{(n)}$ is antisymmetric and $\bigl(u_h^{(n)}\bigr)'$ is symmetric, leading to $u_h^{(n)}(0)=0$ for $k=0$ and $\bigl(u_h^{(n)}\bigr)'(0)=0$ for $k=\frac{\pi}{h}$. \end{itemize} Moreover, the Bloch modes on each edge of a gap attain different symmetries. Therefore, for PCs with mirror symmetry, the Bloch mode vanishes at one of the two edges of the band gap. \end{lemm} \begin{lemm}\label{Lemma:topo} The symmetries of $u_h^{(n)}$ are maintained by continuous deformations of the $n^{\rm th}$ band, as long as it remains isolated. \end{lemm} Based on the three previous lemmas and on the definition of the surface impedance, one deduces the following key result: \begin{theo}\label{lemma:limits} The surface impedance for a given crystal with mirror symmetry either varies from $+\infty$ to 0 or from $0$ to $-\infty$ in a gap, depending on the symmetry of the Bloch modes. If the Bloch modes of the left and right crystal at the edges of the gap have the same symmetry then $Z_L$ and $Z_R$ have the same sign and $Z_L+Z_R$ never cancels out: there is no solution to~\eqref{sum_impedances} and therefore no interface mode in the common gap. On the contrary, if there is a symmetry inversion then $Z_L+Z_R$ changes sign in $\clI$; by continuity and monotony, it therefore cancels out a unique time, proving the existence and uniqueness of an interface mode in $\clI$. Moreover, Lemma~\ref{Lemma:topo} ensures that it is topologically protected. \end{theo} \section{Approximating the surface impedance based on HFH with single eigenvalues}\label{Sec:HFH_single} From now on, one considers continuity of the displacement $u_h$ and of the stress $E\frac{\rmd u_h}{\rmd x}$, but all the results can be easily extended to imperfect contact conditions~\cite{Assier2020,Guzina2021}. We start by recalling the effective model obtained by high-frequency homogenization, and then we use it to approximate the surface impedance. \subsection{High-frequency effective model} \subsubsection{Two-scale expansion} To recall the results of high-frequency homogenization, we introduce the following non-dimensionalized quantities and variables: \begin{equation}\label{eq:nondimparam} X = \frac{x}{L}, \quad \delta = \frac{h}{L}, \quad \Omega = \frac{\omega h}{c_0}, \quad \kappa = L k, \quad U_{\delta} (X) = \frac{u_h (x)}{L}, \quad \alpha = \frac{\rho}{\rho_0}, \quad \beta=\frac{E}{E_0}, \end{equation} with $\rho_0$, $E_0$ and $c_0=\sqrt{\frac{E_0}{\rho_0}}$ some arbitrary reference physical parameters, and $L$ a macroscopic length much larger than $h$ (typically the slow scale at which the modulation function evolves). In non-dimensionalized coordinates, the eigenvalue problem~\eqref{eq:eigenvalue} becomes \begin{equation}\label{eq:eigenvalue_adim} \begin{cases} {\delta^2}\,\frac{\rmd }{\rmd X}\left(\beta\left(\frac{X}{\delta}\right)\frac{\rmd U_\delta}{\rmd X}\right)+\alpha\left(\frac{X}{\delta}\right)\Omega^2 U_\delta =0 \quad\text{on $(0,\delta)$},\\ U_\delta(\delta^+) = e^{\rmi \kappa\delta}U_\delta(0^+) \quad\text{and}\quad U'_\delta(\delta^+) = e^{\rmi \kappa\delta}U'_\delta(0^+). \end{cases} \end{equation} High-frequency homogenization model relies on the assumption that $\delta$ is a small parameter \begin{equation}\label{small_parameter} \delta \ll 1 \end{equation} and on the introduction of the associated fast spatial scale \begin{equation}\label{fast_variable} \xi = \frac{X}{\delta} \end{equation} which accounts for the fast small-scale variations of the physical parameters. We consider approximations of $U_\delta$ and of $\Omega$ around a non-dimensionalized point $(\kappa^\star,\Omega_0)$ of the dispersion diagram such that $\Omega_0^2$ is a single eigenvalue of~\eqref{eq:eigenvalue_adim}. Following the two-scale expansion technique, one then assumes the following ansatz for the non-dimensionalized field and frequency: \begin{equation}\label{eq:asymptoticrepfield} U_{\delta} (X) = \sum_{j \geqslant 0} \delta^j U_j (X, \xi) \quad\text{and}\quad \Omega^2 = \sum_{\ell \geqslant 0} \delta^{\ell} \Omega_{\ell}^2, \end{equation} where we treat $X$ and $\xi$ as two independent space variables. It implies that \begin{equation}\label{eq:exp_derivative} \frac{\rmd }{\rmd X} \leftrightarrow \frac{\partial}{\partial X} + \frac{1}{\delta} \frac{\partial}{\partial \xi}. \end{equation} Moreover, we assume that the Floquet--Bloch conditions are satisfied by the fast spatial scale: \begin{equation}\label{eq:BF_cond} U_j (X, \xi + 1)=\rme^{\rmi \kappa^\star\delta} U_j (X, \xi). \end{equation} \begin{rema}\label{rk:sg_Omega_sq} In~\eqref{eq:asymptoticrepfield}, $\Omega_\ell$ ($\ell \geq 1$) is not necessarily real, allowing for negative values of $\Omega_\ell^2$. \end{rema} \subsubsection{Leading-order effective equation} One can show that the leading-order field $U_0$ reads: \begin{equation}\label{decompo_u0_single} U_0(X,\xi)=f_0(X)\clU_0(\xi), \end{equation} with $\clU_0$ the eigenfunction associated to $\Omega_0$ for the following eigenvalue problem on a unit cell: \begin{equation}\label{eq:eigenvalue_unit_cell} \begin{cases} \frac{\partial}{\partial\xi }\left(\beta\left(\xi\right)\frac{\partial \clU_0}{\partial\xi}\left(\xi\right)\right)+\alpha\left(\xi\right)\Omega_0^2\clU_0\left(\xi\right) =0 \quad\text{on $(0,1)$},\\ \clU_0(1^+) = e^{\rmi \kappa^\star\delta}\clU_0(0^+) \quad\text{and}\quad \clU'_0(1^+) = e^{\rmi \kappa^\star\delta}\clU'_0(0^+), \end{cases} \end{equation} and $f_0$ the envelop function, or modulation function to be found. At next order, the first-order field can be written as \begin{equation}\label{eq:decompo_u1} U_1(X,\xi)=f_1(X)\clU_0(\xi)+f'_0(X) \clU_1(\xi), \end{equation} with $f_1$ the first-order envelop function and $\clU_1$ solution of the following cell problem \begin{equation}\label{system_U1} \begin{cases} \partial_\xi\left(\beta(\xi) \frac{\partial \clU_1}{\partial_\xi}(\xi)+\beta(\xi)\clU_0(\xi)\right)+\Omega_0^2 \alpha(\xi)\clU_1(\xi)=-\beta(\xi)\frac{\partial \clU_0}{\partial_\xi}(\xi)& \text{on $(0,1)$},\\ \clU_1(1^+)=e^{\rmi \kappa^\star \delta} \clU_1(0^+) \quad\text{and}\quad \clU_1'(1^+)=e^{\rmi \kappa^\star \delta} \clU_1'(0^+), \\ \clU_1 \quad\text{and}\quad \beta \clU_1'+\beta \clU_0 &\text{are continuous}. \end{cases} \end{equation} In practice we write \begin{equation}\label{decompo_U1} \clU_1(\xi) = \clV_1(\xi)-\xi \clU_0(\xi), \end{equation} so that the unknown $\clV_1$ now solves the same ODE than $\clU_0$ with different boundary conditions, namely: \begin{equation}\label{system_V} \begin{cases} \partial_\xi\left(\beta(\xi) \frac{\partial \clV_1}{\partial_\xi}(\xi)\right)+\Omega_0^2 \alpha(\xi)\clV_1(\xi)=0& \text{on $(0,1)$},\\ \clV_1(1^+)-\clU_0(1^+)=e^{\rmi \kappa^\star\delta} \clV_1(0^+) \quad\text{and}\quad \clV_1'(1^+)-\clU'_0(1^+)=e^{\rmi \kappa^\star \delta } \clV_1'(0^+), \\ \clV_1\quad\text{and}\quad\beta \clV_1' &\text{are continuous}. \end{cases} \end{equation} In the case of an approximation around an edge of the dispersion diagram ($\kappa^\star = 0$ or $\frac{\pi}{\delta}$), the sought-after effective equation for $f_0$ reads: \begin{equation}\label{syst_f} T f_0'' +\Omega_2^2f_0=0, \end{equation} with $f_0''$ standing for the second-order derivative of $f_0$ and $\Omega_2$ the quadratic term in~\eqref{eq:asymptoticrepfield} which is still to be found. The effective parameter $T$ is given by \begin{equation}\label{def_Ti} T = \frac{\left\langle \beta\clU_0^2+\beta \clU_1'\clU_0-\beta \clU_1\clU'_0\right\rangle}{\left\langle \alpha\clU_0^2\right\rangle}=\frac{\left\langle \beta \clV_1'\clU_0-\beta \clV_1\clU'_0\right\rangle}{\left\langle \alpha\clU_0^2\right\rangle}=\frac{\left\langle \beta \clW_1\right\rangle}{\left\langle \alpha\clU_0^2\right\rangle}, \end{equation} where one introduces \begin{equation}\label{def:wronskian} \clW_1 = \clV_1'\clU_0-\clV_1\clU'_0 \end{equation} which is the wronskian associated to the second-order ODE $(\beta g')'+\Omega_0^2\alpha g = 0$, satisfying therefore $(\beta \clW_1)'=0$. Together with continuity conditions, it leads to $\beta\clW_1$ constant on (0,1) and to the simplified expression for $T$: \begin{equation}\label{T_simpl} T =\frac{\beta(0^+)\clW_1(0^+)}{C}, \end{equation} where we have introduced \begin{equation} C = \left\langle \alpha\clU_0^2\right\rangle >0. \end{equation} Regarding the dispersion relation, the linear term $\Omega_1$ vanishes i.e. \begin{equation}\label{eq:Omega1vanishes} \Omega_1=0 \end{equation} while the quadratic term at a point $\kappa$ around $\kappa^\star$ is given by: \begin{equation}\label{Omega_2_final_1D} \Omega_2^2= T\left(\kappa-\kappa^\star\right)^2. \end{equation} Using Floquet--Bloch conditions then gives the expression of the envelop function in terms of the wavenumber \begin{equation}\label{expression_f0} f_0(X) = e^{i\left(\kappa-\kappa^\star\right)X}. \end{equation} The computation of these approximations and of the expression of the effective parameter $T$ can be found out in various references~\cite{Assier2020,Craster2010,Meng2021,Touboul2024a}. It is recalled in Appendix~\ref{Sec:App_Single_leading} for the sake of completeness. \subsubsection{First-order effective equation} A large part of the literature on HFH only considers the leading-order approximation. However, higher-order approximations have been investigated in~\cite{Guzina2019,Guzina2021,Meng2021}. In this work we will use the first-order approximation in order to have a satisfying approximation of the surface impedance. Indeed, in the case of centro-symmetric PC, the Bloch modes vanish at the beginning or end of the gap; see Lemma~\ref{remark_symmetry}. Using only the leading-order term $U_0(X,\delta)$ would yield a crude null approximation of the fields in the gap, and hence of the surface impedance. The first-order envelop function $f_1$ introduced in~\eqref{eq:decompo_u1} is solution of the same effective equation as $f_0$ \begin{equation}\label{eq:effective_eq_f1} T f_1'' +\Omega_2^2f_1=0, \end{equation} leading consequently to the same envelop function as at the leading order~\eqref{expression_f0}: \begin{equation}\label{eq:total_first_envelop} f_0(X)+\delta f_1(X) = e^{\rmi \left(\kappa-\kappa^\star\right)X}. \end{equation} Therefore, the first-order approximation of the field is given by \begin{equation}\label{eq:approx_first_order_field} U^{(1)}(X,\xi)= U_0(X,\xi)+\delta U_1(X,\xi) = e^{\rmi \left(\kappa-\kappa^\star\right)X} \clU_0(\xi)+ \delta\,i\left(\kappa-\kappa^\star\right) e^{\rmi \left(\kappa-\kappa^\star\right)X}\clU_1(\xi), \end{equation} with $\clU_0$ the eigenfunction associated to $\Omega_0$ solution of~\eqref{eq:eigenvalue_unit_cell} and $\clU_1$ given by~\eqref{decompo_U1}. Here again, the details of the derivation of the first-order effective equation is recalled in Appendix~\ref{Sec:App_Single_first}. \subsection{Approximation of the surface impedance} In order to predict the occurrence and frequency of a localized interface mode from~\eqref{sum_impedances}, the approach followed in~\cite{Coutant2024} was to compute $Z_L$ and $Z_R$ in~\eqref{eq:def_impedance} for any frequency in the common gap $\clI$. The aim of this section is rather to use the HFH approximation described in the previous section around an edge of the dispersion diagram, to get an approximation of these quantities and therefore of the frequency at which an interface localized mode occurs. {Once the edges of the band gaps are identified, the impedances are evaluated through a single analytical calculation, which is an appealing feature of HFH.} The couple $(\kappa^\star,\Omega_0)$, with $\kappa^\star=\{0,\frac{\pi}{h}\}$, will now stand for one of the four edges of the gaps (beginning or end of the gap, for either the right or left crystal) for which there's a non empty overlap $\clI$. One therefore needs to compute \begin{equation} Z_{R,L} = \pm \frac{v\left(0^\pm,\omega\right)}{E\left(0^\pm\right)v'\left(0^\pm,\omega\right)}, \end{equation} where $v(\cdot,\omega)$ is the unique (up to a multiplicative factor) solution decreasing when $x\to \pm \infty$. Using~\eqref{eq:nondimparam}, the impedances in non-dimensionalized coordinates read \begin{equation} Z_{R,L} = \pm\frac{L}{E_0}\frac{V\left(0^\pm,\Omega\right)}{\beta\left(0^\pm\right)V'\left(0^\pm,\Omega\right)}, \end{equation} where $V(\cdot,\Omega)=\frac{v(\cdot,\omega)}{L}$ the non-dimensionalized mode decreasing at $X\pm \infty$. We start by approximating the numerator, i.e. the field at 0, by $U^{(1)}(0,0)$ by $U^{(1)}$ defined in~\eqref{eq:approx_first_order_field}. The only term which depends on the frequency within the gap is $\kappa$, which satisfies \begin{equation}\label{eq:decompo_kappa_gap} \kappa=\kappa^\star+\rmi \kappa_I,\quad \text{with}\quad \kappa_I\in\bbR. \end{equation} From~\eqref{eq:total_first_envelop} together with the condition of decrease when $x\to \pm \infty$, one gets \begin{equation}\label{expression_f0_kappaI} U^{(1)}(X,\xi)= e^{-\kappa_I X} \clU_0(\xi) -\delta\,\kappa_I e^{-\kappa_I X}\clU_1(\xi), \end{equation} where we choose $\kappa_I>0$ for the right crystal and $\kappa_I<0$ for the left crystal. Therefore, the field at the interface is approximated by \begin{equation}\label{approx_field} V(0,\Omega) \approx U^{(1)}(0,0)=\clU_0(0)-\delta\,\kappa_I \clU_1(0). \end{equation} Using~\eqref{eq:exp_derivative}, we then approximate the denominator, i.e. the stress at the interface, by writing \begin{equation}\label{eq:approx_sigma} \beta \,U_\delta'= \frac{1}{\delta} \beta\,\partial_\xi U_0+\beta\,\partial_X U_0+\beta\,\partial_\xi U_1+\clO(\delta) \end{equation} so that at the interface \begin{equation}\label{eq:approx_sigma1} \beta\left(0^\pm\right) V'\left(0^\pm,\Omega\right) \approx \frac{ \beta\left(0^\pm\right)}{\delta} \left(\clU'_0\left(0^\pm\right)-\delta\,\kappa_I\left(\clU_0\left(0^\pm\right)+\clU_1'\left(0^\pm\right)\right) \right). \end{equation} We also deduce from~\eqref{eq:Omega1vanishes} and~\eqref{Omega_2_final_1D} together with~\eqref{eq:decompo_kappa_gap} inserted in the ansatz on $\Omega^2$~\eqref{eq:asymptoticrepfield} that \begin{equation}\label{Omega_gap} \Omega^2 \approx \Omega_0^2 -\delta^2\kappa_I^2 T \end{equation} so that \begin{equation}\label{kappa_I} \delta\,\kappa_I \approx \pm \sqrt{\frac{\Omega_0^2-\Omega^2}{T}}, \end{equation} with the upper sign for the right crystal, and the lower sign for the left crystal. Finally, one gets \begin{equation} \label{eq:final_impedance_single} Z_{R,L} = \pm {h}\frac{ \clU_0(0)-\delta\,\kappa_I \clU_1(0)}{ E(0^+)\Bigl(\clU'_0(0^+)-\delta\,\kappa_I\bigl(\clU_0(0)+\clU_1'(0^+)\bigr)\Bigr)}. \end{equation} In the above expression, given an edge of a band gap, $\clU_0$, $\clU_1$ and their derivatives are computed once for all (computation is detailed in Appendix~\ref{Sec:Comp_U0}). The only dependence on $\Omega$ is then analytical through $\delta\,\kappa_I$ given by~\eqref{kappa_I}. \begin{rema}\label{sg_T} From~\eqref{Omega_gap}, one deduces that $T$ is negative near the entry of a gap (for which $\Omega^2 \geq \Omega_0^2$) and $T$ is positive near the exit of the gap (for which $\Omega^2 \leq \Omega_0^2$). \end{rema} \begin{rema} The expression of~\eqref{eq:final_impedance_single} allows to recover the result of Lemma~\ref{lemma:monotony}. Indeed, one can compute \begin{equation}\label{eq:deriv_Z} \begin{aligned}[b] &\frac{\rmd Z_{R,L}}{\rmd \Omega}\\ & = \pm {h}\frac{ -\clU_1(0)\frac{\rmd (\delta\kappa_I)}{\rmd \Omega}\Big(\clU'_0(0^+)-\delta\kappa_I\big(\clU_0(0)+\clU_1'(0^+)\big)\Big)+\frac{\rmd (\delta\kappa_I)}{\rmd \Omega}\big(\clU_0(0)+\clU_1'(0^+)\big)\Big(\clU_0(0)-\delta\kappa_I \clU_1(0) \Big)}{ E(0^+)\Big(\clU'_0(0^+)-\delta\kappa_I\big(\clU_0(0)+\clU_1'(0^+)\big)\Big)^2}\\ & = \pm {h}\frac{\rmd (\delta\kappa_I)}{\rmd \Omega}\frac{ \clU_0(0)^2+\clU'_1(0^+)\clU_0(0)-\clU_1(0)\clU'_0(0^+)}{ E(0^+)\Big(\clU'_0(0^+)-\delta\kappa_I\big(\clU_0(0)+\clU_1'(0^+)\big)\Big)^2} \\ & = \pm {h}\frac{\rmd (\delta\kappa_I)}{\rmd \Omega}\frac{\beta(0^+) \clW_1(0^+)}{ \beta(0^+)E(0^+)\Big(\clU'_0(0^+)-\delta\kappa_I\big(\clU_0(0)+\clU_1'(0^+)\big)\Big)^2} \\ & = \pm {h}\frac{\rmd (\delta\kappa_I)}{\rmd \Omega}\frac{T C}{ \beta(0^+)E(0^+)\Big(\clU'_0(0^+)-\delta\kappa_I\big(\clU_0(0)+\clU_1'(0^+)\big)\Big)^2}, \end{aligned} \end{equation} where we have used~\eqref{T_simpl}. Therefore $\frac{\rmd Z_{R,L}}{\rmd \Omega}$ has the same sign as $\pm T \frac{\rmd (\delta\kappa_I)}{\rmd \Omega}$, {where $\kappa_I$ is defined in~\eqref{kappa_I}. The sign of the quantities involved depend on the PC considered (L or R) and on the point where expansion is done (entry or exit of the gap), leading to four cases. Table~\ref{TabSign} sums up the signs involved, using Remark~\ref{sg_T}, and yields to $\pm T \frac{\rmd (\delta\kappa_I)}{\rmd \Omega}<0$ in all cases.} We therefore recover that $Z_{R,L}$ is always decreasing in the gap, as stated in Lemma~\ref{lemma:monotony}. \end{rema} \begin{table}[htbp] \caption{Signs of the quantities involved in $\frac{\rmd Z_{R,L}}{\rmd \Omega}$.} \label{TabSign} \centering \begin{tabular}{|l|llll|} \hline & (L, entry) & (L, exit) & (R, entry) & (R, exit) \\ \hline $\pm$ coefficient in~\eqref{eq:deriv_Z} & - & - & + & + \\ \hline sign of $T$ & - & + & - & + \\ \hline sign of $\frac{\rmd (\delta\kappa_I)}{\rmd \Omega}$& - & + & + & - \\ \hline sign of $\frac{\rmd Z_{R,L}}{\rmd \Omega}$ & - & - & - & -\\ \hline \end{tabular} \end{table} \subsection{Numerical example}\label{Sec:Num_single} \begin{figure}[htbp] \centering \subcaptionbox{\label{fig1a}}{\includegraphics[width=0.48\linewidth]{left_mode_begin_gap}}\hfill \subcaptionbox{\label{fig1b}}{\includegraphics[width=0.48\linewidth]{right_mode_begin_gap}} \\ \subcaptionbox{\label{fig1c}}{\includegraphics[width=0.48\linewidth]{left_mode_end_gap}} \subcaptionbox{\label{fig1d}}{\includegraphics[width=0.48\linewidth]{right_mode_end_gap}} \caption{\label{Fig:Simple_Mode} Mode $\clU_0$ over a unit cell for the left crystal (left row) and right crystal (right row) at the entry of the gap (top row) and the exit of the gap (bottom row) {between bands 7 and 8, at $\kappa^\star=\frac{\pi}{\delta}$}. The dotted vertical line denote the different layers of the unit cell.} \end{figure} Let us consider the bilayered periodic structure described by the following centro-symmetric unit cell on $(0,h)$: \begin{equation}\label{eq:def_bilayer_centro_sym} (\rho,E)(x) = \begin{cases} (\rho_A,E_A) = (1000, E_A) & \text{if}\quad 00. \end{equation} The ansatz is considered around the first eigenvalue $\Omega_0^{(1)}$: \begin{equation}\label{eq:asymptoticnearby} U_{\delta} (X) = \sum_{j \geqslant 0} \delta^j U_j (X, \xi) \quad\text{and}\quad \Omega^2 = \left(\Omega_0^{(1)}\right)^2+\sum_{\ell \geqslant 1} \delta^{\ell} \Omega_{\ell}^2. \end{equation} To take into account the coupling between both eigenvalues, the leading-order wavefield reads \begin{equation}\label{eq:nearby_zero_order} U_0(X,\xi) = f_0^{(1)}(X)\clU_0^{(1)}(\xi)+f_0^{(2)}(X)\clU_0^{(2)}(\xi), \end{equation} where $\clU_0^{(1)}$ and $\clU_0^{(2)}$ are the {eigenvectors} associated to $\Omega_0^{(1)}$ and $\Omega_0^{(2)}$, respectively. This expression is consistent with the decomposition onto the basis of two independent eigenfunctions in the case $\gamma=0$ where we get a Dirac point. In the literature, a similar but not identical approach has been followed~\cite{Fefferman2014}, in the sense that the size of the gap is of order $\delta$. However, the meaning of the small parameter differs: in~\cite{Fefferman2014}, the authors start from a configuration with a Dirac point, and add a perturbation of order $\delta$, while here we work with a fixed configuration without Dirac point ($\gamma\neq 0$). High-frequency homogenization comes with an effective equation for the vector \begin{equation}\label{eq:def_F0} {\bm{F}}_0(X) = \left(f_0^{(1)}(X) \,\, f_0^{(2)}(X)\right)^T \end{equation} containing the two modulation functions. It is also useful to introduce the following coefficients \begin{equation}\label{eq:def_coeff_nearby} C_i = \left\langle\alpha \bigl(\clU_0^{(i)}\bigr)^2\right \rangle \;\text{ for } i=1,2 \quad\text{and}\quad \clW = \left\langle\beta\bigl(\clU_0^{(1)}\partial_\xi\clU_0^{(2)}-\partial_\xi\clU_0^{(1)}\clU_0^{(2)}\bigr) \right \rangle:=\langle \beta w \rangle. \end{equation} Contrary to the single (not close) eigenvalues case, the linear term $\Omega_1^2$ in the expansion \begin{equation}\label{expansion_nearby} \Omega^2 = \bigl(\Omega_0^{(1)}\bigr)^2+\delta\Omega_1^2 +\clO\bigl(\delta^2\bigr) \end{equation} does not cancel out and is involved in the effective equation for $\bm{F}_0$, which is given by \begin{equation}\label{eff_ODE_nearby} \bm{F}_0'(X) = \frac{\Omega_1^2}{\clW} \bbN \bm{F}_0(X), \quad\text{with}\quad \bbN = \begin{pmatrix} 0& \left(1-\frac{\gamma}{\Omega_1^2}\right) C_2 \\ - C_1 & 0 \end{pmatrix}. \end{equation} By developing the above system, one gets for $j=1,2$ \begin{equation} T_n^2\ \bigl(f_0^{(j)}\bigr)''+\Omega_1^4\left(1-\frac{\gamma}{\Omega_1^2}\right)\,f_0^{(j)}=0, \end{equation} where the effective parameter in the nearby case is defined by \begin{equation}\label{eq:effective_parameter_nearby} T_n = \frac{\clW}{\sqrt{C_1 C_2 }}. \end{equation} Regarding the dispersion relation, the linear term at a point $\kappa$ around $\kappa{^\star}$ is given by the following expression \begin{equation}\label{eq:Omega1_nearby} \delta\Omega_1^2 = \frac{\delta}{2}\left(\gamma \mp \sqrt{4T_n^2\left(\kappa^\star-\kappa\right)^2+\gamma^2} \right), \end{equation} where the upper (resp. lower) sign designs the branch emerging from $(\Omega_0^{(1)}$ (resp. $(\Omega_0^{(2)}$). Here again, these results can be found in~\cite{Assier2020,Guzina2019,Touboul2024a} and are rederived in Appendix~\ref{Sec:App_Nearby} for the sake of completeness. \subsection{Approximation of the surface impedance} The effective equation~\eqref{eff_ODE_nearby} combined with Floquet--Bloch conditions gives \begin{equation}\label{eq:expr_F0_nearby} \bm{F}_0(X) = \left( \begin{array}{c} -\rmi \frac{\clW\left(\kappa -\kappa^\star\right)}{C_1\Omega_1^2 } \\ 1 \end{array} \right) e^{\rmi \left(\kappa -\kappa^\star\right)X}, \end{equation} as detailed in Appendix~\ref{Sec:App_Nearby}. Therefore, in the bandgap, i.e. where~\eqref{eq:decompo_kappa_gap} holds, the leading-order field reads \begin{equation}\label{eq:in_gap} U_0(X,\xi) = e^{-\kappa_I X} \left(\frac{\clW \kappa_I\delta}{C_1\delta\Omega_1^2 } \clU_0^{(1)}(\xi)+\clU_0^{(2)}(\xi) \right), \end{equation} where we choose again $\kappa_I>0$ for the right crystal and $\kappa_I<0$ for the left crystal to have the decreasing condition at infinity. We therefore get the following approximations for the impedance: \begin{equation} Z_{R,L} \approx \pm h\frac{U_0(0,0)}{E(0^+)\partial_\xi U_0(0,0^+)}= \pm h\frac{\frac{\clW \kappa_I\delta}{C_1\delta\Omega_1^2 } \clU_0^{(1)}(0)+\clU_0^{(2)}(0) }{E(0^+)\left(\frac{\clW \kappa_I\delta}{C_1\delta\Omega_1^2 } {\clU_0^{(1)}}'(0^+)+{\clU_0^{(2)}}'(0^+) \right)}. \end{equation} This expression can be further simplified by taking the dominant term in $\clW$ defined by~\eqref{eq:def_coeff_nearby}. The first step is to notice, using~\eqref{eq:distance_nearby}, that \begin{equation} (\beta w)'(\xi) = -\delta\gamma\alpha(\xi)\clU_0^{(1)}(\xi)\clU_0^{(2)}(\xi). \end{equation} Therefore \begin{equation} \beta(\xi) w(\xi) = \beta(0^+) w(0^+) -\delta\gamma\int_0^\xi\alpha(s)\clU_0^{(1)}(s)\clU_0^{(2)}(s)\rmd s \end{equation} leading to \begin{equation} \begin{aligned} \clW &= \langle\beta w\rangle\\ &= \beta(0^+) w(0^+)-\delta\gamma \int_0^1(1-s)\alpha(s)\clU_0^{(1)}(s)\clU_0^{(2)}(s)\rmd s\\ &= \beta(0^+) w(0^+)+\clO(\delta). \end{aligned} \end{equation} Consequently, we will use the following final approximation for the impedance: \begin{equation}\label{eq:impedance_nearby} Z_{R,L} = \pm \frac{h}{E(0^+)}\frac{\frac{\beta(0^+) w(0^+) \kappa_I\delta}{C_1\delta\Omega_1^2 } \clU_0^{(1)}(0)+\clU_0^{(2)}(0) }{ \frac{\beta(0^+) w(0^+) \kappa_I\delta}{C_1\delta\Omega_1^2 } {\clU_0^{(1)}}'(0^+)+{\clU_0^{(2)}}'(0^+)}. \end{equation} \begin{rema} Again, this approximation of the impedance allows to recover the result of monotony of Lemma~\ref{lemma:monotony}. Indeed, one can compute \begin{equation}\label{eq:dZ_nearby} \begin{aligned}[b] \frac{\rmd Z_{R,L}}{\rmd \Omega} &= \pm\frac{h}{E(0^+)}\frac{\beta(0^+) w(0^+)}{C_1} \frac{\rmd }{\rmd \Omega}\left(\frac{ \kappa_I\delta}{\delta\Omega_1^2 }\right)\frac{w(0+)}{\left(\frac{\beta(0^+) w(0^+) \kappa_I\delta}{C_1\delta\Omega_1^2 } {\clU_0^{(1)}}'(0^+)+{\clU_0^{(2)}}'(0^+) \right)^2} \\ &= \pm \frac{\rmd }{\rmd \Omega}\left(\frac{ \kappa_I\delta}{\delta\Omega_1^2 }\right)\frac{h w(0+)^2}{E_0 C_1\left(\frac{\beta(0^+) w(0^+) \kappa_I\delta}{C_1\delta\Omega_1^2 } {\clU_0^{(1)}}'(0^+)+{\clU_0^{(2)}}'(0^+) \right)^2}. \end{aligned} \end{equation} By inverting~\eqref{eq:Omega1_nearby} in the band gap, one can show that \begin{equation} \frac{\kappa_I\delta}{\delta\Omega_1^2} = \pm \sqrt{\frac{\frac{\gamma\delta}{\delta\Omega_1^2}-1}{T_n^2}}. \end{equation} Therefore, \begin{equation} \begin{aligned}[b] \frac{\rmd }{\rmd \Omega}\left(\frac{ \kappa_I\delta}{\delta\Omega_1^2 }\right) &= \mp \frac{1}{2} \gamma\delta \frac{\bigl(\delta\Omega_1^2\bigr)'}{\bigl(\delta\Omega_1^2\bigr)^2T_n^2}\left(\frac{\frac{\gamma\delta}{\delta\Omega_1^2}-1}{T_n^2}\right)^{\frac{-1}{2}} \\ &= - \frac{1}{2} \gamma\delta \frac{\bigl(\delta\Omega_1^2\bigr)'}{\delta\Omega_1^2 T_n^2}\frac{1}{\kappa_I\delta}, \end{aligned} \end{equation} where $ (\delta\Omega_1^2)' =\frac{\rmd (\delta\Omega_1^2)}{\rmd \Omega}\geq 0$. We also recall that $\kappa_I\delta$ is positive (resp. negative) for the right\linebreak (resp. left) crystal. Consequently, $\frac{\rmd }{\rmd \Omega}\bigl(\frac{ \kappa_I\delta}{\delta\Omega_1^2 }\bigr)$ is negative (resp. positive) for the right (resp. left) crystal, and by~\eqref{eq:dZ_nearby} $Z_{R,L}$ is decreasing. \end{rema} \begin{figure}[!hb] \centering \subcaptionbox{\label{fig5a}}{ \includegraphics[width=0.48\linewidth]{left_crystal_real_nearby}}\hfill \subcaptionbox{\label{fig5b}}{ \includegraphics[width=0.48\linewidth]{right_crystal_real_nearby}} \\ \subcaptionbox{\label{fig5c}}{ \includegraphics[width=0.48\linewidth]{left_crystal_imagl_nearby}} \hfill \subcaptionbox{\label{fig5d}}{\includegraphics[width=0.48\linewidth]{right_crystal_imagl_nearby}} \caption{\label{Fig:DD_nearby} Dispersion diagrams for the left crystal (left row) and right crystal (right row). Two real branches (first row) and imaginary part in the band gap defined by these two branches (bottom row). Physical parameters are given by~\eqref{eq:def_bilayer_centro_sym} and~\eqref{param_crystal_LR}. The red plain line denotes the exact dispersion diagram obtained by solving the eigenvalue problem~\eqref{eq:eigenvalue_adim} at any frequency, while the dashed black line is the HFH approximation around the edge for nearby eigenvalues given by~\eqref{eq:Omega1_nearby}.} \end{figure} \subsection{Numerical example} We deal with the same example as in Section~\ref{Sec:Num_single} but we now use the nearby eigenvalues HFH effective model instead of the classical one. The approximations of the dispersion diagrams for the same branches as previously are plotted in Figure~\ref{Fig:DD_nearby}. We can see a much longer-lived approximation on the range of wavenumber considered for both the real and imaginary part. We especially get a much better fit within the all range of the band gap. The approximation of the surface impedance for both the left and right crystal is represented together with the exact one in Figure~\ref{Fig:Impedance_nearby}. Contrary to the approximation obtained with the classical HFH model, the approximation of the surface impedance is now continuous and presents an excellent agreement with the exact surface impedance even around the middle of the gap.\linebreak Both the approximations of $Z_L(\Omega)+Z_R(\Omega)$ and their exact values cancel out for $\Omega=15.708$ showing a very precise prediction of the frequency at which the interface protected mode occurs. {For this frequency, the topological protected mode obtained by the nearby HFH approximation is plotted in Figure~\ref{Fig:ModeTopoHFHNearby}}. \begin{figure}[!h] \centering \subcaptionbox{\label{fig6a}}{ \includegraphics[width=0.48\linewidth]{Z_L_nearby}} \hfill \subcaptionbox{\label{fig6b}}{\includegraphics[width=0.48\linewidth]{Z_R_nearby}} \\ \subcaptionbox{\label{fig6c}}{ \includegraphics[width=0.48\linewidth]{ZR_ZL_full_nby}}\hfill \subcaptionbox{\label{fig6d}}{ \includegraphics[width=0.48\linewidth]{ZR_ZL_full_nby_zoom}} \caption{\label{Fig:Impedance_nearby} Top row: surface impedance for the left crystal (a) and right crystal (b). The red plain line denotes the exact surface impedance computed at any frequency, while the dashed black line is the HFH approximation for nearby eigenvalues given by~\eqref{eq:impedance_nearby}. Bottom row: sum of the surface impedances (c) and zoom (d) around the zero i.e. the frequency of the topologically protected mode and its approximation.} \end{figure} \begin{figure}[!hb] \centering \includegraphics[width=0.5\linewidth]{mode_topo_nearby} \caption{Topologically protected mode obtained by HFH in the nearby eigenvalues case. The two semi-infinite PCs are in contact at $x=0$ and we represent in blue (resp. orange) the mode in the left (resp. right) crystal. The dashed vertical lines represent the boundaries of the periodic cells.} \label{Fig:ModeTopoHFHNearby} \end{figure} \begin{figure}[!h] \centering \subcaptionbox{\label{fig8a}}{ \includegraphics[width=0.48\linewidth]{N10-F857-Iter40000}}\hfill \subcaptionbox{\label{fig8b}}{ \includegraphics[width=0.48\linewidth]{TopoN10-Iter59263}} \caption{{Time-domain} simulation of wave propagation across two finite networks of $N=10$ cells (the vertical lines denote the locations of the interfaces between the different cells) glued at the center of the domain. A source emits a monochromatic wave impacting the left PC. In (a), the frequency is $\Omega=13.461$, where no topological interface mode exists. In (b), the frequency is $\Omega=15.708$ where the total surface impedance vanishes. An interface mode is generated at the interface between PC-L and PC-R denoted by the red plain vertical line.} \label{Fig:ModeTopoSimu} \end{figure} Lastly, we illustrate the scattering of waves by this finite-size network of cells~\cite{KalozoumisPR18}. For this purpose, we use time-domain simulations. The integration methods are described, for example, in~\cite{Darche2025}. Transparent boundary conditions are used on the left and right extremities, in order to avoid spurious reflections by the edges of the computational domain. A monochromatic source is placed on the left of the PC-L. In Figure~\ref{fig8a}, the frequency is $\Omega=13.461$, which belongs to the common gap between bands 6 and 7. In this gap, no topologically protected interface mode exists. As a consequence, one observes the generation of an evanescent mode that decreases exponentially from the left boundary of the PC-L. In Figure~\ref{fig8b}, one uses $\Omega=15.708$, identified as the frequency of the topological interface mode thanks to HFH. We then observe the generation of an evanescent mode centred on the interface between PC-L and PC-R. \section{Conclusion}\label{SecConclu} We have showed that it is possible to use effective models got from high-frequency homogenization in order to get an approximation of the surface impedance for a given phononic crystal. Consequently, we can deduce, whether or not, and at which possible frequency, a topologically protected interface state would occur when putting two semi-infinite phononic crystals in contact. We compared the quality of the approximation deduced from classical high-frequency effective models to the one deduced from nearby eigenvalues high-frequency effective models. We showed that the latter gives a much better approximation of the surface impedance within a band gap, and therefore a more precise prediction of the occurrence of a topologically protected state. {A natural extension of this work would be its generalization to higher dimensional systems, such as analogues of the quantum spin Hall effect or the Valley Hall effect~\cite{Yves2022}. High frequency homogenization would allow us not only to obtain approximate dispersion relation for topological edge or interface waves, but also provide powerful approximation to study their transport properties~\cite{Torres2024}.} \section*{Declaration of interests} The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations. %\section*{Conflicts of interest} % %The authors declare no competing financial interest. \section*{Acknowledgments} The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. \appendix %\def\theequation{\Alph{section}\arabic{equation}} %\setcounter{equation}{0} %\def\thetable{\Alph{section}\arabic{table}} %\setcounter{figure}{0} \section{Derivation of the effective model for single eigenvalues}\label{Sec:App_Single} \subsection{Leading-order effective equation}\label{Sec:App_Single_leading} We recall that in the present case we consider points at the edges of the Brillouin zone, that is at $\kappa^\star=0$ or $\kappa^\star=\frac{\pi}{\delta}$ which is associated to periodic and antiperiodic conditions, respectively, for the fields. Using the ansatz~\eqref{eq:asymptoticrepfield} and~\eqref{eq:exp_derivative}, the governing equation in~\eqref{eq:eigenvalue_adim} becomes \begin{align}\label{eq:allinonegoveq} \sum_{j \geqslant 0} \left[\delta^j \frac{\partial}{\partial\xi}\left(\beta \frac{\partial U_j}{\partial \xi}\right) + \delta^{j+1}\left(\beta\frac{\partial^2U_j}{\partial X\partial \xi}+\frac{\partial}{\partial\xi}\left(\beta \frac{\partial U_j}{\partial X}\right)\right)+\delta^{j+2}\beta\frac{\partial^2U_j}{\partial X^2} + \sum_{\ell \geqslant 0} \delta^{\ell + j} \Omega_{\ell}^2 \alpha U_j \right] = 0. \end{align} Collecting the terms of order $\delta^0$, we get in $(0,1)$: \begin{equation}\label{syst_order0_main} \frac{\partial}{\partial\xi}\left(\beta\frac{\partial U_0}{\partial\xi}\right)+\Omega_0^2 \alpha U_0 = 0 \end{equation} together with continuity for $U_0$ and $\beta U'_0$ at possible interfaces, and $ U_0 (X, \xi + 1)=\rme^{\rmi \kappa^\star\delta} U_0 (X, \xi)=\pm U_0 (X, \xi)$. From this equation, one deduces the form of the leading-order field~\eqref{decompo_u0_single}. Going to next order and collecting the term of order $\delta$, one gets on $(0,1)$: \begin{equation}\label{eq:system_order1} \frac{\partial}{\partial\xi}\left(\beta\left(\frac{\partial U_1}{\partial\xi}+\frac{\partial U_0}{\partial X}\right)\right)+\beta\frac{\partial^2 U_0}{\partial X\partial\xi}+\Omega_0^2 \alpha U_1+\Omega_1^2 \alpha U_0 = 0 \end{equation} together with continuity for $U_1$ and $\beta\left(\frac{\partial U_1}{\partial\xi}+\frac{\partial U_0}{\partial X}\right)$, and periodicity/antiperiodicity conditions. One then computes $\langle \eqref{syst_order0_main} \times U_1 - \eqref{eq:system_order1}\times U_0\rangle$ and gets after integration by parts \begin{equation} \left\langle \beta \frac{\partial U_0}{\partial X}\frac{\partial U_0}{\partial \xi} -\beta \frac{\partial^2U_0}{\partial X\partial\xi}U_0-\Omega_1^2\alpha U_0^2\right\rangle = 0. \end{equation} After making use of~\eqref{decompo_u0_single}, this reduces to \begin{equation} \Omega_1^2 \left\langle \alpha \clU_0^2 \right\rangle = 0 \end{equation} leading therefore to $\Omega_1=0.$ Together with~\eqref{eq:system_order1}, this gives the expression~\eqref{eq:decompo_u1} for $U_1$. Collecting terms of order $\delta^2$, we get on $(0,1)$: \begin{equation}\label{eq:syst_u2} \frac{\partial}{\partial\xi}\left(\beta\left(\frac{\partial U_2}{\partial\xi}+\frac{\partial U_1}{\partial X}\right)\right)+\beta\frac{\partial^2 U_1}{\partial X\partial\xi}+\beta\frac{\partial^2 U_0}{\partial X^2}+\Omega_0^2 \alpha U_2+\Omega_2^2 \alpha U_0 = 0, \end{equation} together with continuity for $U_2$ and $\beta\left(\frac{\partial U_2}{\partial\xi}+\frac{\partial U_1}{\partial X}\right)$, and periodicity/antiperiodicity for $U_2$. One then computes $\langle \eqref{syst_order0_main} \times U_2 - \eqref{eq:syst_u2}\times U_0\rangle$ and gets the final effective equation for $f_0$~\eqref{syst_f} after integration by parts and making use of~\eqref{decompo_u0_single} and~\eqref{eq:decompo_u1}. In this effective equation, $T\neq 0$ but it can be either positive or negative. Therefore $f_0$ can be written \begin{equation}\label{eq:f_0_exp} f_0(X) = Ae^{\rmi \sqrt{\frac{\Omega_2^2}{T}}X}+Be^{-\rmi \sqrt{\frac{\Omega_2^2}{T}}X}. \end{equation} \begin{rema}\label{Remark-Omega2T} The value of $T$ and the expression of $f_0$ is valid for both the real and imaginary bands. On a pass-band (real band), $\frac{\Omega_2^2}{T}$ will be positive so that $f_0$ describes propagative waves, while in a gap (imaginary band) $\frac{\Omega_2^2}{T}$ will be negative describing therefore evanescent waves. \end{rema} We now recall that $\xi=\frac{X}{\delta}$ and that the field has to satisfy the Floquet--Bloch conditions, therefore leading to \begin{equation}\label{eq:BF_U0} U_0(X+\delta,\xi+1) = \rme^{\rmi \kappa \delta} U_0(X,\xi). \end{equation} Using in the above equation~\eqref{decompo_u0_single} and the periodicity/antiperiodicity of $\clU_0$ one gets \begin{equation} f_0(X+\delta)\rme^{\rmi \kappa^\star \delta}\clU_0(\xi) =\rme^{\rmi \kappa \delta} f_0(X)\clU_0(\xi). \end{equation} We recall that $\kappa^\star = 0$ or $\frac{\pi}{\delta}$ is an edge of the dispersion diagram, while $\kappa$ is a point around $\kappa^\star$ at which we are using our effective model. Therefore \begin{equation} f_0(X+\delta) = \rme^{\rmi \left(\kappa-\kappa^\star\right) \delta}f_0(X). \end{equation} Using~\eqref{eq:f_0_exp} in the above equation yields \begin{equation} Ae^{\rmi \sqrt{\frac{\Omega_2^2}{T}}\delta} = A\rme^{\rmi \left(\kappa-\kappa^\star\right) \delta} \quad\text{and}\quad Be^{-\rmi \sqrt{\frac{\Omega_2^2}{T}}\delta} = B\rme^{\rmi \left(\kappa-\kappa^\star\right) \delta}, \end{equation} and keeping in mind that $\frac{\Omega_2^2}{T}$ is positive leads to \begin{equation} \begin{cases} \sqrt{\frac{\Omega_2^2}{T}} =\left(\kappa-\kappa^\star\right) & \text{and}\quad B=0\quad \text{for } \kappa^\star=0 \\ \sqrt{\frac{\Omega_2^2}{T}} =\left(\kappa^\star-\kappa\right) & \text{and}\quad A=0\quad \text{for } \kappa^\star=\frac{\pi}{\delta}. \end{cases} \end{equation} In any case, up to a constant ($A$ or $B$) that we take equal to $1$, $f_0$ satisfies~\eqref{expression_f0}. \subsection{First-order effective equation}\label{Sec:App_Single_first} Inserting~\eqref{decompo_u0_single} and~\eqref{eq:decompo_u1} in~\eqref{eq:syst_u2}, and using~\eqref{syst_f} one gets: \begin{equation} \frac{\partial}{\partial\xi}\left(\beta\frac{\partial U_2}{\partial\xi}\right)+\Omega_0^2\alpha U_2 + \left(\frac{\partial}{\partial\xi}\left(\beta \clU_0\right)+\beta\clU_0'\right)f_1'+\left(\frac{\partial}{\partial\xi}\left(\beta \clU_1\right)+\beta \clU_1'+(\beta-\alpha T)\clU_0\right)f_0''. \end{equation} Therefore $U_2$ can be written as \begin{equation}\label{eq:decompo_U2} U_2(X,\xi) = f_2(X)\clU_0(\xi)+f_1'(X)\clU_1(\xi)+f_0''(X)\clU_2(\xi), \end{equation} with $\clU_2$ solution of the following cell problem: \begin{equation}\label{system_W} \begin{cases} \partial_\xi(\beta(\xi) \clU_2'(\xi)+\beta(\xi)\clU_1(\xi))+\Omega_0^2 \alpha(\xi)\clU_2(\xi)\\ \qquad=-\beta(\xi)\clU_1'(\xi)-(\beta(\xi)-\alpha(\xi) T)\clU_0(\xi) &\text{on $(0,1)$},\\ \clU_2(1^+)=e^{\rmi \kappa^\star\delta} \clU_2(0^+) \quad\text{and}\quad \clU_2'(1^+)=e^{\rmi \kappa^\star\delta} \clU_2'(0^+), \\ \clU_2\quad\text{and}\quad\beta \clU_2'+\beta \clU_1&\text{are continuous}. \end{cases} \end{equation} Collecting terms of order $\delta^3$, we get on $(0,1)$: \begin{multline}\label{eq:syst_u3} \frac{\partial}{\partial\xi}\left(\beta(\xi)\left(\frac{\partial U_3}{\partial\xi}(X,\xi)+\frac{\partial U_2}{\partial X}(X,\xi)\right)\right) \\+\beta(\xi)\frac{\partial^2 U_2}{\partial X\partial\xi} +\beta(\xi)\frac{\partial^2 U_1}{\partial X^2}+\Omega_0^2 \alpha(\xi)U_3(X,\xi)+\Omega_2^2 \alpha(\xi)U_1(X,\xi) = 0, \end{multline} together with continuity for $U_3$ and $\beta\bigl(\frac{\partial U_3}{\partial\xi}+\frac{\partial U_2}{\partial X}\bigr)$, and periodicity/antiperiodicity for $U_3$. We notice that there is no $\Omega_3$ in the above expression. It is justified by the following remark. \begin{rema} For a given branch, the mapping that associates to a wavenumber the frequency is holomorphic (except when eigenvalues are no longer single ones, which is excluded here). Consequently, the function $\Omega(\kappa)$ can be represented by a Taylor series expansion. Combined with reciprocity, i.e., $\Omega(\kappa)=\Omega(-\kappa)$, this gives that all odd powers in the ansatz on $\Omega$ vanish, as proved previously for $\Omega_1$. \end{rema} One then computes $\langle\eqref{syst_order0_main} \times U_3 -\eqref{eq:syst_u3}\times U_0\rangle$ and gets, using integration by parts together with~\eqref{decompo_u0_single}, \eqref{eq:decompo_u1}, \eqref{eq:decompo_U2} and~\eqref{syst_f} differentiated once: \begin{equation}\label{eq:integration_u3} T f_1''(X)+\Omega_2^2 f_1(X)+\frac{1}{\left\langle \alpha\clU_0^2\right\rangle}\left\langle \beta\left(\clU_2'\clU_0-\clU_2\clU_0'\right)+(\beta-\alpha T)\clU_1\clU_0\right\rangle f_0'''(X) = 0. \end{equation} Moreover, $\langle \eqref{system_W} \times \clU_1-\eqref{system_U1}\times \clU_2\rangle$ leads to \begin{equation} \left\langle \beta\left(\clU_2'\clU_0-\clU_2\clU_0')+(\beta-\alpha T\right)\clU_1\clU_0\right\rangle=0. \end{equation} We therefore get~\eqref{eq:effective_eq_f1}. \section{Derivation of the effective model for nearby eigenvalues}\label{Sec:App_Nearby} Because $\clU_0^{(i)}$ ($i=1,2$) are eigenfunctions associated to $\Omega_0^{(i)}$, they satisfy: \begin{equation}\label{eq:eigenvalue_nearby} \frac{\partial}{\partial\xi}\left(\beta \frac{\partial \clU_0^{(i)}}{\partial\xi}\right)+\bigl(\Omega_0^{(i)}\bigr)^2\alpha\clU_0^{(i)} =0. \end{equation} We introduce the notation~\eqref{eq:eigenvalue_nearby}$^{(i)}$ for this equation at a given $i$, and we denote by $U_0^{(i)}=f_0^{(i)}\clU_0^{(i)}$ for $i=1,2$. One can then check that \begin{equation} \frac{\partial}{\partial\xi}\left(\beta \frac{\partial U_0}{\partial\xi}\right)+\bigl(\Omega_0^{(1)}\bigr)^2\alpha U_0 = -\delta\alpha\gamma U_0^{(2)}, \end{equation} where we recall that $\gamma$ defines the distance through the two eigenvalues through~\eqref{eq:distance_nearby}. The right-hand side residual term will therefore modify the governing equation for $U_1$, which, collecting the terms of order $\delta$, reads: \begin{equation}\label{eq:system_order1_nearby} \frac{\partial}{\partial\xi}\left(\beta\left(\frac{\partial U_1}{\partial\xi}+\frac{\partial U_0}{\partial X}\right)\right)+\beta\frac{\partial^2 U_0}{\partial X\partial\xi}+\bigl(\Omega_0^{(1)}\bigr)^2 \alpha U_1+\Omega_1^2 \alpha U_0-\alpha\gamma U_0^{(2)} = 0 \end{equation} together with continuity for $U_1$ and $\beta\bigl(\frac{\partial U_1}{\partial\xi}+\frac{\partial U_0}{\partial X}\bigr)$ and periodicity/antiperiodicity for $U_1$. One then considers $\langle\eqref{eq:eigenvalue_nearby}^{(i)}\times U_1-\eqref{eq:system_order1_nearby}\times \clU_0^{(i)}\rangle $ for $i=1,2$ to obtain two equations that can be recast as the first-order ODE~\eqref{eff_ODE_nearby}. One notes that to get~\eqref{eff_ODE_nearby}, we used $\langle\alpha \clU_0^{(1)}\clU_0^{(2)} \rangle=0$ since the eigenfunctions are orthogonal. The two eigenvalues of $\bbN$ are given by \begin{equation}\label{eq:eigenvalues_nearby} \lambda_i = (-1)^i \sqrt{\left(\frac{\gamma}{\Omega_1^2}-1\right)C_1C_2} \end{equation} with $i=1,2$ while the associated eigenvectors are \begin{equation}\label{eq:eigenvectors_nearby} \bm{\clF}_{\lambda_i} = \left(-\frac{\lambda_i}{C_1},\,1\right)^T. \end{equation} Therefore, $\bm{F}_0$ can be written as \begin{equation}\label{eq:F0_eigenvalues_nearby} \bm{F}_0(X) = \clA_1 \bm{\clF}_{\lambda_1} e^{\Omega_1^2/\clW\lambda_1 X}+\clA_2 \bm{\clF}_{\lambda_2} e^{\Omega_1^2/\clW\lambda_2 X}. \end{equation} Moreover, applying the Floquet--Bloch conditions~\eqref{eq:BF_U0} yields \begin{equation} \rme^{\rmi \kappa^\star \delta}\left(f_0^{(1)}(X+\delta)\clU_0^{(1)}(\xi)+ f_0^{(2)}(X+\delta)\clU_0^{(2)}(\xi)\right) =\rme^{\rmi \kappa \delta}\left(f_0^{(1)}(X)\clU_0^{(1)}(\xi) +f_0^{(2)}(X)\clU_0^{(2)}(\xi)\right). \end{equation} Since $\clU_0^{(1)}$ and $\clU_0^{(2)}$ are linearly independent, we get from the above equation \begin{equation}\label{eq:bf_f0_nearby} \bm{F}_0(X+\delta) = \rme^{\rmi \left(\kappa-\kappa^\star\right)\delta }\bm{F}_0(X). \end{equation} Combining~\eqref{eq:bf_f0_nearby} and~\eqref{eq:F0_eigenvalues_nearby}, we obtain \begin{equation} \clA_1 \bm{\clF}_{\lambda_1} e^{\Omega_1^2/\clW\lambda_1 \delta } = \clA_1 \bm{\clF}_{\lambda_1} \rme^{\rmi \left(\kappa-\kappa^\star\right)\delta } \quad\text{and}\quad \clA_2 \bm{\clF}_{\lambda_2} e^{\Omega_1^2/\clW\lambda_2 \delta } = \clA_2 \bm{\clF}_{\lambda_2} \rme^{\rmi \left(\kappa-\kappa^\star\right)\delta }. \end{equation} Since $\lambda_1=-\lambda_2$, we deduce \begin{equation}\label{eq:syst_nearby} \begin{cases} \left(\Omega_1^2/\clW\right)\lambda_1 =\rmi \left(\kappa-\kappa^\star\right) & \text{and}\quad \clA_2=0, \\ \text{or}\\ \left(\Omega_1^2/\clW\right)\lambda_2 =\rmi \left(\kappa-\kappa^\star\right) & \text{and}\quad \clA_1=0. \end{cases} \end{equation} These two cases can be summarized, using~\eqref{eq:eigenvalues_nearby} and~\eqref{eq:effective_parameter_nearby}, as \begin{equation} \pm\frac{\Omega_1^2}{T_n}\sqrt{\frac{\gamma}{\Omega_1^2}-1}=\rmi \left(\kappa-\kappa^\star\right), \end{equation} whose roots are given by~\eqref{eq:Omega1_nearby}. Moreover, using~\eqref{eq:syst_nearby} together with~\eqref{eq:F0_eigenvalues_nearby} and~\eqref{eq:eigenvectors_nearby}, we get the expression of the envelop function in terms of the wavenumber given by~\eqref{eq:expr_F0_nearby}. \section{Computing \texorpdfstring{$\clU_0$}{U\_0} and \texorpdfstring{$\clU_1$}{U\_1}}\label{Sec:Comp_U0} In this section, we detail the computation of $\clU_0$ and $\clU_1$ for the centro-symmetric unit cell described by~\eqref{eq:def_bilayer_centro_sym}. The strategy is exactly the same for any piecewise-constant unit cell. \subsection{Computation of \texorpdfstring{$\clU_0$}{U\_0}} We recall that $\clU_0$ is solution of~\eqref{eq:eigenvalue_unit_cell} where $\alpha$ and $\beta$ are piecewise-constant. $\clU_0$ can therefore be written as \begin{equation} \clU_0(\xi) = \begin{cases} \clC_1 \cos\left(\frac{\Omega_0}{\bar{c}_A}\xi\right)+ \clC_2 \sin\left(\frac{\Omega_0}{\bar{c}_A}\xi\right) & \text{for $0<\xi<\varphi_1$ }\\[0.25em] \clC_3 \cos\left(\frac{\Omega_0}{\bar{c}_B}\xi\right)+ \clC_4 \sin\left(\frac{\Omega_0}{\bar{c}_B}\xi\right) & \text{for $\varphi_1<\xi<\varphi_2$ }\\[0.25em] \clC_5 \cos\left(\frac{\Omega_0}{\bar{c}_A}\xi\right)+ \clC_6 \sin\left(\frac{\Omega_0}{\bar{c}_A}\xi\right) & \text{for $\varphi_2<\xi<1$,} \end{cases} \end{equation} where $\varphi_1 = h_A/(2h)$, $\varphi_2=(h-h_A)/(2h)$, and $\bar{c}_j=\sqrt{\beta_j/\alpha_j}$ for $j=A,B$. Using the Floquet--Bloch conditions $\clU_0(1^+) = e^{\rmi \kappa^\star\delta}\clU_0(0^+)$ and $\clU'_0(1^+) = e^{\rmi \kappa^\star\delta}\clU'_0(0^+)$ and the continuity conditions for $\clU_0$ and $\beta \clU_0$ at 0, $\varphi_1$ and $\varphi_2$, the integration constants $\clC_i$ satisfy the following linear system \begin{equation}\label{linear_syst_U0} \bbM(\Omega_0,\kappa^\star) \bm{C} = \bm{0}, \end{equation} with $\bm{C}$ the vector containing the integration constants, i.e. $\bm{C}[i]=\clC_i$ for $i=1,\dots,6$, and $\bbM(\Omega_0,\kappa^\star)$ the $6\times 6$ matrix encapsulating the continuity and Floquet--Bloch conditions which is given by: \[ \begin{pmatrix} e^{\rmi \kappa^\star\delta} & 0 & 0 & 0 & -\cos\left(\frac{\Omega_0}{\bar{c}_A}\right) & -\sin\left(\frac{\Omega_0}{\bar{c}_A}\right) \\[0.35em] 0 & \bar{Z}_A e^{\rmi \kappa^\star\delta} & 0 & 0 & \bar{Z}_A \sin\left(\frac{\Omega_0}{\bar{c}_A} \right)& - \bar{Z}_A \cos\left(\frac{\Omega_0}{\bar{c}_A}\right) \\[0.35em] \cos\left(\frac{\Omega_0}{\bar{c}_A}\varphi_1 \right) & \sin\left(\frac{\Omega_0}{\bar{c}_A} \varphi_1 \right) & -\cos\left(\frac{\Omega_0}{\bar{c}_B}\varphi_1 \right) & -\sin\left(\frac{\Omega_0}{\bar{c}_B} \varphi_1\right) & 0 & 0 \\[0.35em] - \bar{Z}_A \sin\left(\frac{\Omega_0}{\bar{c}_A} \varphi_1 \right) & \bar{Z}_A \cos\left(\frac{\Omega_0}{\bar{c}_A} \varphi_1 \right) & \bar{Z}_B \sin\left(\frac{\Omega_0}{\bar{c}_B}\varphi_1 \right) & -\bar{Z}_B \cos\left(\frac{\Omega_0}{\bar{c}_B} \varphi_1 \right) & 0 & 0 \\[0.35em] 0 & 0 & \cos\left(\frac{\Omega_0}{\bar{c}_B}\varphi_2 \right) & \sin\left(\frac{\Omega_0}{\bar{c}_B} \varphi_2\right) & -\cos\left(\frac{\Omega_0}{\bar{c}_A}\varphi_2 \right) & -\sin\left(\frac{\Omega_0}{\bar{c}_A }\varphi_2 \right) \\[0.35em] 0 & 0 & -\bar{Z}_B \sin\left(\frac{\Omega_0}{\bar{c}_B} \varphi_2 \right) & \bar{Z}_B \cos\left(\frac{\Omega_0}{\bar{c}_B} \varphi_2 \right) & \bar{Z}_A \sin\left(\frac{\Omega_0}{\bar{c}_A} \varphi_2 \right) & -\bar{Z}_A \cos\left(\frac{\Omega_0}{\bar{c}_A} \varphi_2 \right) \end{pmatrix} \] where $\bar{Z}_j=\sqrt{\alpha_j\beta_j}$ with $j=A,B$. The frequency $\Omega_0$ is then solution of $\det(\bbM)(\Omega_0,\kappa^\star)=0$, and $\bm{C}$ is given by the kernel of $\bbM(\Omega_0,\kappa^\star)$ for the value of $\Omega_0$ previously found. \subsection{Computation of \texorpdfstring{$\clU_1$}{U\_1}} Since $\clU_1$ reads $\clU_1=\clV_1-\xi\clU_0$ with $\clV_1$ solution of~\eqref{system_V} and $\clU_0$ already known, we only have to compute $\clV_1$. Similarly as for $\clU_0$, $\clV_1$ can be written as \begin{equation} \clV_1(\xi) = \begin{cases} \clD_1 \cos\left(\frac{\Omega_0}{\bar{c}_A}\xi\right)+ \clD_2 \sin\left(\frac{\Omega_0}{\bar{c}_A}\xi\right) & \text{for } 0<\xi<\varphi_1 \\[0.35em] \clD_3 \cos\left(\frac{\Omega_0}{\bar{c}_B}\xi\right)+ \clD_4 \sin\left(\frac{\Omega_0}{\bar{c}_B}\xi\right) & \text{for }\varphi_1<\xi<\varphi_2\\[0.35em] \clD_5 \cos\left(\frac{\Omega_0}{\bar{c}_A}\xi\right)+ \clD_6 \sin\left(\frac{\Omega_0}{\bar{c}_A}\xi\right) & \text{for }\varphi_2<\xi<1, \end{cases} \end{equation} where $\Omega_0$ is now known, and with some new unknown coefficients $\clD_i$ which are gathered in $\bm{D}$ solution of \begin{equation}\label{syst_V1} \bbM \bm{D} = \bm{V} \end{equation} with \begin{equation} \bm{V} = \begin{pmatrix} & -\clC_5\cos\left(\frac{\Omega_0}{\bar{c}_A}\right)-\clC_6\sin\left(\frac{\Omega_0}{\bar{c}_A}\right)\\ & -\bar{Z}_A e^{\rmi \kappa^\star\delta}\clC_2 \\ &0\\ &0\\ &0\\ &0 \end{pmatrix}. \end{equation} One therefore only needs to invert numerically~\eqref{syst_V1} to get $\clV_1$ and therefore $\clU_1$. \printbibliography \end{document}