# A homogenized model accounting for dispersion, interfaces and source points for transient waves in 1D periodic media

R Cornaggia, Bruno Lombard

► To cite this version:

R Cornaggia, Bruno Lombard. A homogenized model accounting for dispersion, interfaces and source points for transient waves in 1D periodic media. ESAIM: Mathematical Modelling and Numerical Analysis, 2023, 57 (3), pp.1413-1444. 10.1051/m2an/2023027. hal-03652455v2

HAL Id: hal-03652455

https://hal.science/hal-03652455v2

Submitted on 24 Mar 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A HOMOGENIZED MODEL ACCOUNTING FOR DISPERSION, INTERFACES AND SOURCE POINTS FOR TRANSIENT WAVES IN 1D PERIODIC MEDIA.

R. CORNAGGIA<fnref n="1" /> AND B. LOMBARD<fnref n="2" />

**Abstract.** A homogenized model is proposed for linear waves in 1D microstructured media. It combines second-order asymptotic homogenization (to account for dispersion) and interface correctors (for transmission from or towards homogeneous media). A new bound on a second-order effective coefficient is proven, ensuring well-posedness of the homogenized model whatever the microstructure. Based on an analogy with existing enriched continua, the evolution equations are reformulated as a dispersive hyperbolic system. The efficiency of the model is illustrated via time-domain numerical simulations. An extension to Dirac source terms is also proposed.

**1991 Mathematics Subject Classification.** 35F45,74J10,74Q10.

March 2022.

## 1. INTRODUCTION

Understanding, modelling and controlling wave propagation in heterogeneous media is of major interest in numerous engineering domains, e.g. seismic simulation and protection, sound and vibration attenuation, non-destructive testing, etc. The particular case of periodically varying media, see e.g. [28] for sonic and phononic crystals, gathered much attention. Indeed the interaction between waves and the periodic medium produces pronounced dispersive features, including band-gaps (“forbidden” range of frequencies for which waves decay exponentially) at higher frequencies. On the modelling side, robust methods exploit the periodicity to establish effective models linking the salient macroscopic features of the wave propagation to the microstructure. These effective models enable efficient numerical simulation, and they also provide a solid ground to implement topological optimization algorithms of architected materials [3, 14].

In this paper, the focus is on the long wavelength regime, compared to the characteristic size of the microstructure. In this regime (prior to the first band-gap), higher-order models accounting for dispersion effects are studied since [35], notably thanks to double-scale asymptotic homogenization methods [9], see e.g. [3, 14] for models of scalar waves in two-dimensional media. One-dimensional situations are particularly well documented: [1, 4, 5, 18, 19, 27] among other computed the higher-order homogenized coefficients and designed well-posed models thanks to the “Boussinesq trick” [2] that enables the partial or total permutation of space and time derivatives without loss of the asymptotic order of approximation. This trick was generalized in [41] that established a whole family of models. Introducing an additional modelling tool, [36] recently proposed a

<footnote>Keywords and phrases: waves in periodic media, high-order homogenization, effective transmission conditions.</footnote>

<footnote n="1">¹ Sorbonne Université, CNRS, Institut ∂'Alembert, UMR 7190, Paris, France (remi.cornaggia@sorbonne-universite.fr)</footnote>

<footnote n="2">² Aix-Marseille Univ, CNRS, Centrale Marseille, LMA UMR 7031, Marseille, France (lombard@lma.cnrs-mrs.fr)</footnote>
reformulation of these homogenized wave equations into hyperbolic systems well-suited for time-domain simulation [23], based on a phenomenological stress-gradient model presented in [21].

On the other hand, using these models in realistic bounded domains, i.e. addressing either boundary conditions or transmission conditions towards another homogeneous or microstructured medium, is an issue that is still a major research topic. In two or three dimensions, oscillating boundary layers that appear near these boundaries of interfaces, can be accounted for using corrector functions to complement the inner expansion and to obtain convergence estimates, see e.g. [34] for eigenvalue problems in bounded domains, [26] for elliptic problems with interfaces and [10, 11, 20] and the references therein for time-harmonic transmission problems. However, (i) these boundary or transmission correctors are complex objects whose theoretical study is quite involved [6, 22] and (ii) their direct computation is often as costly as the one of the original problem posed on the microstructured domain.

Fewer works address the “practical” implementation of these correctors, i.e. the design of enriched boundary or transmission conditions suitable for numerical implementation and simulation. A notable exception is the contribution of S. Fliss and collaborators [8, 20, 40]. These authors consider two-dimensional time-harmonic problems in half-planes or strips and address interface neighborhoods thanks to matched asymptotic expansions or “enriched” double-scale expansions accounting explicitly for local correctors. Also using matched expansions, effective transmission conditions are built for scattering problems by layered “slabs” of materials with interfaces either orthogonal [31] or aligned [32] with the layers, and in the context of water waves by [33]. For purely one-dimensional propagation, [15] proposes boundary and transmission conditions for wave problems homogenized up to second order, still for time-harmonic loads and fields. This methodology is also applied to non-uniformly oscillating media in [37].

Finally, the time-domain simulation of these homogenized wave models and associated boundary and transmission conditions is almost inexistant, to the best or our knowledge, although a proposal was made recently in [8, Chap. 5] for a half-plane with a Dirichlet boundary condition. In this context, the main purpose of the present paper is to build a well-posed “total” model to simulate 1D transient wave propagation through microstructured domains, combining (i) the family of second-order models from [41] and their reformulation as hyperbolic systems from [36], and (ii) the effective transmission conditions proposed in [15], adapted to the time-domain context to ensure well-posedness. An extension of the proposed tools to Dirac source points immersed in the microstructure is also given.

In Section 2, dispersive models for one-dimensional waves obtained from second-order homogenization are recalled, and a new result on the second-order homogenized coefficient is given (Proposition 1). Then the proposed hyperbolic system is presented and the conditions for its well-posedness are established (Proposition 2). Section 3 gives relevant corrections for transmission conditions between homogeneous and homogenized media, and states the stability conditions for the total model incorporating the hyperbolic system and these conditions (Proposition 3). Section 4 presents the treatment of source points. The improvements brought by the proposed models are illustrated by a set of time-domain simulations in Section 5, and Section 6 finally summarizes the findings of the papers and proposes several extensions. Auxiliary definitions and proofs are gathered in Appendix A.

## 2. SECOND-ORDER HOMOGENIZED MODEL IN UNBOUNDED SPACE

This section focuses on free waves in an unbounded periodic elastic medium, characterized by density $ρ_ℓ(x) = ρ(x/ℓ)$ and Young's modulus $E_ℓ(x) = E(x/ℓ)$, in terms of 1-periodic functions $(ρ, E)$ and the periodicity length $ℓ$. The material displacement is denoted $u_ℓ(x, t)$, and $v_ℓ = ∂_t u_ℓ$ and $σ_ℓ$ are the velocity and stress fields. Combining the momentum balance equation $∂_t(ρ_ℓ v_ℓ) = ∂_x σ_ℓ$ and the elastic constitutive law $σ_ℓ = E_ℓ ∂_x u_ℓ$, the source-free wave motion is then described equivalently by the wave equation satisfied by $u_ℓ$,
or the first-order system satisfied by $(v_\ell, \sigma_\ell)$:

$$ \rho_{\ell} \partial_{tt} u_{\ell} - \partial_{x}(E_{\ell} \partial_{x} u_{\ell}) = 0 \quad \Leftrightarrow \quad \begin{cases} \partial_{t} v_{\ell} - \frac{1}{\rho_{\ell}} \partial_{x} \sigma_{\ell} = 0, \\ \partial_{t} \sigma_{\ell} - E_{\ell} \partial_{x} v_{\ell} = 0. \end{cases} \qquad (1) $$

The main results obtained by second-order homogenization at low frequencies and wavenumbers, and moderate variations of material parameters $(E, \rho)$ (so that the ratio between period and wavelength is the only small scale of the problem), are summarized below. Since the derivation of these results is now classical, they are not justified: the reader is referred to [41] and the references therein for more details on the formal expansions, to [15] for error estimates in the time-harmonic case, and to [1, 27] for error estimates in the transient case.

**Remark 1.** Rigorous asymptotic analysis should be done in terms of a non-dimensional parameter $\varepsilon = \ell/\lambda$, where $\lambda$ is a "reference" wavelength that serves as a spatial scale. For time-harmonic problems, this wavelength can be defined, see e.g. [15, Remark 1]. In the present work, the focus is on modelling rather than on asymptotic analysis, and transient waves with broad frequency –and therefore wavelength– spectra are considered, hence all the expansions and estimates below are formally given in terms of $\ell$ for simplicity: $o(1)$, $o(\ell)$ and $o(\ell^2)$ remainders indicate leading-, first- and second-order approximations. Discussions on the choice of reference wavelength $\lambda$ and the influence of ratio $\varepsilon$ are postponed to numerical illustrations in Section 5.1.

## 2.1. Formal expansions, cell functions and macroscopic fields

The two-scale formal expansion begins with the following ansatz on the solution:

$$ u_{\ell}(x, t) = \sum_{n} \ell^{n} u_{n}(x, x/\ell, t), \quad y \mapsto u_{n}(x, y, t) \text{ is 1-periodic,} $$

where the spatial dependency is split between the “slow” space variable $x$ and its “fast” counterpart $y = x/\ell$. Injecting this expression into the wave equation (1) and solving iteratively the resulting cascade of equations, see [41], one finds separated variable solutions of the form:

$$ \begin{aligned} u_0(x, y, t) &= \bar{u}_0(x, t) \\ u_1(x, y, t) &= \bar{u}_1(x, t) + P_1(y) \partial_x \bar{u}_0(x, t) \\ u_2(x, y, t) &= \bar{u}_2(x, t) + P_1(y) \partial_x \bar{u}_1(x, t) + P_2(y) \partial_{xx} \bar{u}_0(x, t) \end{aligned} \qquad (2) $$

written in terms of slowly varying macroscopic fields $\bar{u}_j$ and fastly oscillating correctors involving 1-periodic cell functions $P_j$, modulated by the successive space derivatives of the previous macroscopic fields. The cell functions are solutions of static cell problems posed on the scaled periodicity cell $[0, 1]$, as recalled in Appendix A.1.

The formal process also provides wave equations satisfied by the macroscopic fields $\bar{u}_n$, whose source terms depend on the previously determined fields $\{\bar{u}_j\}_{j<n}$. Here we prefer to work with total macroscopic fields $U_n$ instead, that concatenate all the contributions of order up to $n$. In what follows, we accordingly define:

$$ U_0 = \bar{u}_0, \quad U_1 = \bar{u}_0 + \ell \bar{u}_1, \quad U_2 = \bar{u}_0 + \ell \bar{u}_1 + \ell^2 \bar{u}_2. $$

Section 2.2 below provides and discusses the wave equations satisfied by these fields.

Finally, (i) defining the total macroscopic velocities $V_n = \partial_t U_n$, and stresses $S_n = E_0 \partial_x U_n$, where $E_0$ is the effective Young's modulus to be defined later, (ii) summing the contributions of eqs. (2), (iii) using replacements e.g. $\bar{u}_0 = U_1 + o(1)$ in the corrector terms and (iv) neglecting higher-order contributions, the total $n$-th order
approximations of $(u_\ell, v_\ell, \sigma_\ell)$ are given by:

$$
\begin{align}
u_{\ell}(x,t) &= U_n(x,t) + \sum_{j=1}^{n} \ell^j P_j \left(\frac{x}{\ell}\right) \partial_x^j U_n(x,t) + o(\ell^n), \\
v_{\ell}(x,t) &= V_n(x,t) + \sum_{j=1}^{n} \ell^j P_j \left(\frac{x}{\ell}\right) \partial_x^j V_n(x,t) + o(\ell^n), \tag{3} \\
\sigma_{\ell}(x,t) &= S_n(x,t) + \sum_{j=1}^{n} \ell^j Q_j \left(\frac{x}{\ell}\right) \partial_x^j S_n(x,t) + o(\ell^n),
\end{align}
$$

where specific stress cell functions $Q_j$ are introduced [15], see Appendix A.1 and Remark 2 below.

**Remark 2.** The literature often focuses on the displacement $u_\ell$ and on the wave equation it satisfies (left of (1)). The stress field $\sigma_\ell = E_\ell \partial_x u_\ell$, particularly relevant when considering boundary or transmission conditions, is then treated as a byproduct of the analysis, whose approximation accuracy suffers from the differentiation. On the other hand, the system featuring the two fields $(v_\ell, \sigma_\ell)$ in (1) emphasizes the similar roles played by these fields, that should be reflected by equally accurate approximations. This is the main motivation for the functions $Q_j$ associated to the stress field $\sigma_\ell$, as proposed in this 1D waves context by [15]. The two-fields formalism is also used for deriving boundary and transmission correctors [31–33] and providing associated error estimates in higher dimensions, see [10, 11] and the references therein.

## 2.2. Homogenized wave models up to second order

Effective equations satisfied by the macroscopic fields $(U_n, V_n, S_n)$, also provided by the two-scale approach, are now discussed at orders 0, 1 and 2.

**Homogenized model at leading order.** This model is given similarly to (1) by a wave equation on $U_0$ or a system on $(V_0, S_0)$:

$$ \partial_{tt} U_0 - c_0^2 \partial_{xx} U_0 = 0 \quad \Leftrightarrow \quad \left\{ \begin{array}{l} \partial_t V_0 - \frac{1}{\rho_0} \partial_x S_0 = 0, \\ \partial_t S_0 - E_0 \partial_x V_0 = 0. \end{array} \right. \qquad (4) $$

The constant homogenized coefficients at leading order $\rho_0$ and $E_0$ are given in 1D by the usual formula:

$$ \rho_0 = \int_0^1 \rho(y)dy, \quad E_0 = \left( \int_0^1 E^{-1}(y)dy \right)^{-1}, \qquad (5) $$

i.e. the arithmetic mean and the harmonic mean of their oscillating counterparts $(\rho, E)$. The celerity $c_0$ is defined naturally as:

$$ c_0 = \sqrt{E_0/\rho_0}. \qquad (6) $$

**First-order homogenized model.** As shown by [12, Section 2.3] or [15, Lemma 1], the first-order macroscopic fields $(U_1, V_1, S_1)$ also satisfy the homogenized equations (4), i.e. there is no first-order contribution to the wave equation.

**Second-order homogenized model.** The model (4) is valid in the quasistatic limit, i.e. for wavelengths much longer than $\ell$. For shorter wavelengths, dispersive effects (i.e. frequency-dependent wave velocity) cannot be neglected and must be accounted for in the homogenized models. This is done by incorporating second-order terms into the wave equation, i.e. using the equations satisfied by $(U_2, V_2, S_2)$. In [41], it was shown that
second-order expansions lead to a family of enriched wave equations for the macroscopic displacement $U_2$:

$$ \partial_{tt}U_2 - c_0^2\partial_{xx}U_2 - \ell^2 \left( \beta_x c_0^2\partial_{xxxx}U_2 - \beta_m\partial_{xxtt}U_2 - \frac{\beta_t}{c_0^2}\partial_{tttt}U_2 \right) = 0, \quad \text{with} \quad \beta_x - \beta_m - \beta_t = \beta, \quad (7) $$

where the homogenization process provides the values of the homogenized velocity $c_0$ given by (6), and the second-order coefficient $\beta$, for which a new expression is provided here:

**Proposition 1.** The second-order homogenized coefficient $\beta$ can be computed in terms of the material coefficients $(\rho, E)$, the homogenized coefficients $(\rho_0, E_0)$ given by (5) and the first cell functions $(P_1, Q_1)$:

$$ \beta = \frac{1}{2} \int_{0}^{1} [P_{1}(y) - Q_{1}(y)]^{2} \left( \frac{E_{0}}{E(y)} + \frac{\rho(y)}{\rho_{0}} \right) dy. $$

As a consequence, $\beta$ is non-negative.

The proof of this result is given in Appendix A.2 and relies on reciprocity identities on the cell functions $(P_1, Q_1, P_2, Q_2)$ entering the initial definition (60) of $\beta$ given in [41]. This important non-negativity property is emphasized as it plays a key role in the analysis of the ensuing model (Section 2.3). It is trivially satisfied for bilaminates, see eq. (67), and it was postulated in [36, Sect. 4.3] after many numerical tests on multi-laminated composites, but it is proven here for any periodic medium for the first time, to the best of our knowledge, for non-constant $\rho$: see [19, Remark 1] when $\rho = 1$.

**Remark 3.** The family of models (7), parametrized by the coefficients $(\beta_x, \beta_m, \beta_t)$, is obtained thanks to the wave equation $\partial_{tt}U_2 = c_0^2\partial_{xx}U_2 + O(\ell^2)$ at leading-order, which enables to exchange second-order space and time derivatives in higher-order terms. This possibility, sometimes called “Boussinesq trick” [2], was used by [19] to exhibit a well-posed model with only “mixed” derivatives $(\beta_x = \beta_t = 0, \beta_m = -\beta < 0)$, to replace ill-posed models with fourth-order space derivatives $(\beta_x = \beta > 0, \beta_t = \beta_m = 0)$ that arise naturally from the formal homogenization process (when $\rho = 1$). Later, [27] also studied the “mixed” model, derived an equivalent Korteweg-de-Vries equation and proved error estimates for long-time propagation. The study [1] also proved well-posedness and established error estimates for a two-parameter family of models $(\beta_x - \beta_m = \beta, \beta_t = 0)$. The family (7) introduced in [41] is a generalization of these models which provides two degrees of freedom to choose a particular model. The model can therefore be optimized e.g. by fitting the dispersion curve of the periodic medium, see [36], while still satisfying the well-posedness requirements.

In higher dimensions, terms involving fourth-order space derivatives, associated with a high-order tensor-valued coefficient, are in general impossible to remove. However, the Boussinesq trick still can be performed to obtain well-posed models [2].

Finally, in this work as in [15], only the cases where $\beta_x = 0$ are considered because (i) a fourth-order space derivative in the wave equation would come along with additional boundary or transmission conditions that are not easy to define (although a proposal is made by [8, Chap. 5] for a Dirichlet boundary condition) and (ii) the analysis performed in [36] showed that this term leads to unstable models when $\beta_x > 0$, involving an additional modelling constraint. The so-called (mt) model, featuring fourth-order mixed and time derivatives, is therefore used hereafter:

$$ \partial_{tt}U_2 - c_0^2\partial_{xx}U_2 + \ell^2 \left( \beta_m\partial_{xxtt}U_2 + \frac{\beta_t}{c_0^2}\partial_{tttt}U_2 \right) = 0, \quad \text{with} \quad -\beta_m - \beta_t = \beta, \quad (8) $$

and only one degree of freedom is left to choose the couple $(\beta_m, \beta_t)$.

To investigate the properties of this wave equation, an equivalent system is now presented. Using this formalism enables to reuse theoretical results specific to hyperbolic systems (stability, existence of solutions), as well as dedicated numerical methods (integration schemes, discretization of interfaces). An overview may be found in the reference books [23, 29].
## 2.3. Stress-gradient formalism

The upcoming system comes from a stress-gradient model introduced in [21] and studied in [36]. This phenomenological model features four fields: the usual displacement $u$ and stress $\sigma$, but also a microdisplacement $\phi$ and a microstress $r$ to account for microstructural effects. These fields satisfy the constitutive relations and equilibrium equations:

$$ \begin{cases} \sigma = E_{\text{SG}}\partial_x(u + \phi), \\ r = D\phi, \end{cases} \quad \begin{cases} \rho_{\text{SG}}\partial_{tt}u = \partial_x\sigma, \\ \rho_{\text{SG}}J\partial_{tt}\phi = \partial_x\sigma - r, \end{cases} \qquad (9) $$

written in terms of density and Young's modulus ($\rho_{\text{SG}}, E_{\text{SG}}$) and two micro-stiffness and micro-inertia coefficients ($D, J$). All these parameters are strictly positive real numbers, and $D$ should not be confused with a differentiation operator.

**Relationships with the second-order homogenized model.** Combining the relations (9), as done in [21, Section 5.3], provides a fourth-order enriched wave equation for $u$:

$$ \partial_{tt}u - \frac{E_{\text{SG}}}{\rho_{\text{SG}}}\partial_{xx}u - \frac{E_{\text{SG}}}{D}(1+J)\partial_{xxtt}u + \frac{\rho_{\text{SG}}J}{D}\partial_{tttt}u = 0, \qquad (10) $$

which is formally the same equation than the (mt) homogenized equation given by (8). More precisely, the equivalence is obtained by identifying the coefficients ($E_{\text{SG}}, \rho_{\text{SG}}, J, D$) as:

$$ E_{\text{SG}} = E_0, \quad \rho_{\text{SG}} = \rho_0, \quad J = \frac{\beta_t}{\beta} \quad \text{and} \quad D = \frac{E_0}{\ell^2\beta}. \qquad (11) $$

Moreover, identifying the terms in the two wave equations (8) and (10), and using the relations (9) and (11), the fields of the stress-gradient model satisfy the following relations with the macroscopic displacement $U_2$ and stress $S_2 = E_0\partial_x U_2$:

$$ \begin{cases} u = U_2 \\ \sigma = S_2 + E_0\partial_x\phi = S_2 + \ell^2\beta\partial_x r \end{cases} \quad \text{and} \quad \begin{cases} \partial_{xx}\phi = \ell^2 \left[ -\frac{\beta_m}{c_0^2}\partial_{xxtt}U_2 - \frac{\beta_t}{c_0^4}\partial_{tttt}U_2 \right], \\ \partial_{xx}r = \frac{E_0}{\beta} \left[ -\frac{\beta_m}{c_0^2}\partial_{xxtt}U_2 - \frac{\beta_t}{c_0^4}\partial_{tttt}U_2 \right]. \end{cases} \qquad (12) $$

Several remarks may be done from the above relations:

* The microdisplacement $\phi$ is a second-order term, while the microstress $r$ is in fact a leading-order contribution: from (9) one has $r = \partial_x\sigma + o(1)$.

* The parameter $J$ controls the choice of (mt) model: $J = 1$ corresponds to the “time” model ($t$) (with $\beta_m = 0, \beta_t = -\beta$), and $J = 0$ corresponds to the “mixed” model ($m$) (with $\beta_m = -\beta, \beta_t = 0$) in the wave equation (10). However $J = 0$ makes no sense in the original stress-gradient system (9) that therefore cannot describe the ($m$) model.

**Equivalent system.** To pursue towards the derivation of a hyperbolic system equivalent to the wave equation (8), the total macroscopic velocity $w$, incorporating the microscopic fluctuations, and the microscopic velocity $\varphi = o(\ell)$ are introduced:

$$ w = \partial_t(u + \phi), \quad \varphi = \partial_t\phi. \qquad (13) $$
Then, slightly reformulating the equations (9), the four fields $(w, \sigma, \varphi, r)$ are found to satisfy the system:

$$ \left\{ \begin{array}{l} \partial_t w - \frac{a}{\rho_0} \partial_x \sigma = -\frac{a-1}{\rho_0} r, \\ \partial_t \sigma - E_0 \partial_x w = 0, \\ \partial_t \varphi - \frac{a-1}{\rho_0} \partial_x \sigma = -\frac{a-1}{\rho_0} r, \\ \partial_t r = \frac{E_0}{\ell^2 \beta} \varphi, \end{array} \right. \qquad (14) $$

where the convenient parameter $a$ controls the choice of $(mt)$ model:

$$ a = \frac{J+1}{J} = \frac{-\beta_m}{\beta_t}. $$

**Remark 4.** In [36] a similar system was introduced, featuring the classical velocity $v = \partial_t u$ instead of $w$. The present choice slightly simplifies the notation when introducing interfaces in Section 3.

**Total fields approximations.** Finally, approximations of the total fields $(v_\ell, \sigma_\ell)$ are given by the expansions (3) in terms of the macroscopic fields $(V_2, S_2)$. From the relations (12), these macroscopic fields can be recovered a posteriori from the stress-gradient velocity and stress $(w, \sigma)$ and the auxiliary fields $(\varphi, r)$ as:

$$ V_2 = w - \varphi \quad \text{and} \quad S_2 = \sigma - \ell^2 \beta \partial_x r. \qquad (15) $$

However, a simpler way to compute the total fields, without computing the fields $(V_2, S_2)$, is to use the following approximations:

$$ \begin{aligned} v_{\ell}(x,t) &= w(x) + \ell P_1 \left(\frac{x}{\ell}\right) \partial_x w(x,t) + \ell^2 P_2 \left(\frac{x}{\ell}\right) \partial_{xx} w(x,t) - \varphi(x,t) + o(\ell^2), \\ \sigma_{\ell}(x,t) &= \sigma(x,t) + \ell Q_1 \left(\frac{x}{\ell}\right) \partial_x \sigma(x,t) + \ell^2 Q_2 \left(\frac{x}{\ell}\right) \partial_{xx} \sigma(x,t) - \ell^2 \beta \partial_x r(x,t) + o(\ell^2). \end{aligned} \qquad (16) $$

They are formally justified by the second-order amplitude of the fluctuations $V_2 - w$ and $S_2 - \sigma$.

## 2.4. Model properties

The system (14) above can be recast in matricial form:

$$ \partial_t U + A \cdot \partial_x U = S \cdot U, \qquad (17) $$

where $U := (w, \sigma, \varphi, r)^T$ and:

$$ A = \begin{bmatrix} 0 & -a/\rho_0 & 0 & 0 \\ -E_0 & 0 & 0 & 0 \\ 0 & -(a-1)/\rho_0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \quad \text{and} \quad S = \begin{bmatrix} 0 & 0 & 0 & -(a-1)/\rho_0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -(a-1)/\rho_0 \\ 0 & 0 & E_0/(\ell^2\beta) & 0 \end{bmatrix}. $$

**Hyperbolicity.** The matrices $A$ and $S$ have eigenvalues:

$$ \mathrm{Sp}(A) = \{0, 0, \pm c_0 \sqrt{a}\} \quad \text{and} \quad \mathrm{Sp}(S) = \{0, 0, \pm c_0 / (\ell \sqrt{-\beta_t})\}. $$
The eigenspaces of the matrix **A** are:

$$
\begin{align*}
\mathbf{A} \cdot \mathbf{U} &= 0 && \Rightarrow \mathbf{U} \in \operatorname{Vect} \{(0, 0, 1, 0)^\mathrm{T}, (0, 0, 0, 1)^\mathrm{T}\}, \\
\mathbf{A} \cdot \mathbf{U} &= \pm c_0 \sqrt{a} \mathbf{U} && \Rightarrow \mathbf{U} \in \operatorname{Vect} \{(a, \mp \rho_0 c_0 \sqrt{a}, a - 1, 0)^\mathrm{T}\}.
\end{align*}
$$

For $a \neq 0$, all these vectors are linearly independent and therefore they form a basis of $\mathbb{R}^4$. To obtain a hyperbolic system, the eigenvalues of **A** need to be real. To ensure that the null solution is stable, the eigenvalues of **S** need to be imaginary. These two conditions imply that the coefficients $(a, \beta_m, \beta_t)$ should satisfy:

$$
\left\{
\begin{array}{ll}
\text{Hyperbolicity:} & a > 0 \\
\text{Stability:} & \beta_t > 0
\end{array}
\right.
\implies
\left\{
\begin{array}{ll}
\frac{\beta_m}{\beta_t} < 0 & \implies \left\{ \begin{array}{l} \beta_m < 0, \\ \beta_t > 0 \end{array} \right. \\
\beta_t > 0 &
\end{array}
\right.
\tag{18}
$$

As already pointed out by [36] and omitted in previous studies [15], these conditions complement the condition
$-\beta_m - \beta_t = \beta$ coming from the homogenization process.

For latter use, the following exponential matrix, needed in the numerical integration of (17), is introduced:

$$
\exp(\mathbf{S}t) =
\begin{bmatrix}
1 & 0 & \cos \omega_0 t - 1 & -\frac{1}{\nu} \sin \omega_0 t \\
0 & 1 & 0 & 0 \\
0 & 0 & \cos \omega_0 t & -\frac{1}{\nu} \sin \omega_0 t \\
0 & 0 & \nu \sin \omega_0 t & \cos \omega_0 t
\end{bmatrix},
\quad (19)
$$

where:

$$
\omega_0 = \frac{c_0}{\ell} \frac{1}{\sqrt{\beta_t}} \quad \text{and} \quad \nu = \frac{\rho_0 c_0}{\sqrt{\beta(a-1)} \ell}.
$$

**Energy conservation.** Introducing the symmetrizer matrix

$$
\mathbf{M} = \begin{bmatrix} \rho_0 & 0 & -\rho_0 & 0 \\ 0 & 1/E_0 & 0 & 0 \\ -\rho_0 & 0 & a\rho_0/(a-1) & 0 \\ 0 & 0 & 0 & \ell^2\beta/E_0 \end{bmatrix},
$$

one obtains:

$$
\mathbf{M} \cdot \mathbf{A} = \begin{bmatrix} 0 & -1 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \quad \text{and} \quad \mathbf{M} \cdot \mathbf{S} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 \\ 0 & 0 & 1 & 0 \end{bmatrix},
$$

i.e. **M** · **A** is symmetric, and **M** · **S** is skew-symmetric. Multiplying the system (17) by **U**^T · **M** and making use of these (skew)symmetries, one obtains:

$$
\frac{1}{2} \partial_t (\mathbf{U}^\mathrm{T} \cdot \mathbf{M} \cdot \mathbf{U}) + \frac{1}{2} \partial_x (\mathbf{U}^\mathrm{T} \cdot \mathbf{M} \cdot \mathbf{A} \cdot \mathbf{U}) = 0. \quad (20)
$$

Considering compactly supported initial conditions or source terms, and integrating over a space-time domain
Ω × [0, T], where Ω is chosen large enough so that **U**(*·*, *T*) = **0** on ∂Ω, one obtains the energy conservation:

$$
\mathcal{E}^{\text{vol}}(T) = \mathcal{E}^{\text{vol}}(0).
$$
The volume energy $\mathcal{E}^{\text{vol}}$ is the sum of a kinetic and elastic energy given as:

$$ \begin{aligned} \mathcal{E}^{\text{vol}} &= \frac{1}{2} \int_{\Omega} \mathbf{U}^{\text{T}} \cdot \mathbf{M} \cdot \mathbf{U} \, dx = \frac{1}{2} \int_{\Omega} \rho_0 \left( w^2 - 2w\varphi + \frac{a}{a-1}\varphi^2 \right) + \frac{1}{E_0} (\sigma^2 + \ell^2\beta r^2) \, dx \\ &= \frac{1}{2} \int_{\Omega} \rho_0 \left( v^2 + \frac{\beta_t}{\beta}\varphi^2 \right) + \frac{1}{E_0} (\sigma^2 + \ell^2\beta r^2) \, dx. \end{aligned} \quad (21) $$

In the second expression, the macroscopic velocity $v = \partial_t u$ is introduced and $1 - a = \beta/\beta_t$ is used to emphasize that the energy is positive since $\beta > 0$ (Prop. 1) and $\beta_t > 0$ (as required by the stability condition (18)).

Finally, the results of this part are summarized in the following proposition:

**Proposition 2.** If the parameters $(\beta_m, \beta_t)$ satisfying the relation $-\beta_m - \beta_t = \beta$ are chosen so that $\beta_m < 0$ and $\beta_t > 0$, the system (14) is hyperbolic, the null solution is stable, and the associated positive volume energy $\mathcal{E}^{\text{vol}}$ defined by (21) is conserved.

### 3. FIRST-ORDER TRANSMISSION CONDITIONS AND TOTAL MODEL

Now that an effective model is proposed for unbounded microstructured media, this section focuses on how to use this model for bounded domains made of microstructured materials. In the same way as higher-order correctors were introduced in the wave equation, higher-order interface correctors will be introduced. These correctors will be defined as in the previous study [15], that focused on time-harmonic problems and in which rigorous error estimates are proven to support this proposition, following earlier works on boundary correctors, see [10, 11]. In the present work, an extension is made to the transient case without rigorous justification: our purpose is to establish a well-posed and practical model rather than to prove asymptotic accuracy (that is nevertheless observed in numerical examples).

To fix ideas, the transmission of waves from a homogeneous domain $\Omega^- = \{x < 0\}$ to a microstructured domain $\Omega^+ = \{x > 0\}$ is first addressed in detail. Then the extension to a “slab” $\Omega =]0, L[$ surrounded at both extremities by homogeneous media is given.

For future use, let us introduce the jump $\llbracket f \rrbracket_0$ and the mean $\langle f \rangle_0$ of a function $f$ across an interface at $x = 0$ as:

$$ \llbracket f \rrbracket_0 = f^+ - f^- \quad \text{and} \quad \langle f \rangle_0 = \frac{1}{2} (f^+ + f^-), \quad \text{with } f^\pm = \lim_{x \to 0^\pm} f(x), $$

and we give the following relations:

$$ f^{\pm} = \langle f \rangle_0 \pm \llbracket f \rrbracket_0/2 \quad \text{and} \quad \llbracket fg \rrbracket_0 = \llbracket f \rrbracket_0 \langle g \rangle_0 + \langle f \rangle_0 \llbracket g \rrbracket_0. \quad (22) $$

#### 3.1. Microstructured transmission problem

The original problem to be approximated is defined as follows: the fields $(v_\ell, \sigma_\ell)$ satisfy the system (1) with constant coefficients $(\rho_-, E_\ell)$ in $\Omega^- = \{x < 0\}$ and with oscillating coefficients $(\rho_\ell, E_\ell)$ in $\Omega^+ = \{x > 0\}$. A perfect interface is considered at $x = 0$, i.e. the velocity and stress are continuous:

$$ \llbracket v_\ell \rrbracket_0 = 0, \quad \llbracket \sigma_\ell \rrbracket_0 = 0. \quad (23) $$

To complete this problem, initial conditions are finally given:

$$ v_\ell(x, 0) = v_\star(x) \quad \text{and} \quad \sigma_\ell(x, 0) = \sigma_\star(x), \quad (24) $$

where $v_\star$ and $\sigma_\star$ are entirely supported in the homogeneous domain $\Omega^-$ to simplify the analysis. Indeed, initial conditions supported in the microstructured domain $\Omega^+$ should be addressed in the homogenization process, which is outside the scope of the present work: again the reader is referred to [27] for insights on these issues
in the case of an unbounded domain. In Section 4, we will study the case of a Dirac source point immersed in the microstructured medium.

In the following parts, the focus will be on the interface conditions at $x = 0$ when the wave propagation in the right domain $\Omega^+$ is described by an effective model.

## 3.2. The homogenized model at leading order

The model at leading order is first presented. The macroscopic fields are noted $(V_0, S_0)$ and the correctors in the approximation (3) are ignored, i.e. the approximations $(v_\ell, \sigma_\ell) \approx (V_0, S_0)$ are used. Using the effective model at the leading order (4) in $\Omega^+$, these fields satisfy (i) the systems:

$$ \Omega^{-} : \left\{ \begin{array}{l} \partial_t V_0 - \frac{1}{\rho_-} \partial_x S_0 = 0, \\ \partial_t S_0 - E_- \partial_x V_0 = 0, \end{array} \right. \quad \Omega^{+} : \left\{ \begin{array}{l} \partial_t V_0 - \frac{1}{\rho_0} \partial_x S_0 = 0, \\ \partial_t S_0 - E_0 \partial_x V_0 = 0, \end{array} \right. \qquad (25) $$

(ii) the interface conditions at $x = 0$ inherited from the conditions (23) (see [10] and [32, Appendix 1.A.]):

$$ \llbracket V_0 \rrbracket_0 = 0, \quad \text{and} \quad \llbracket S_0 \rrbracket_0 = 0, \qquad (26) $$

and (iii) the initial conditions (24) with $(v_\ell, \sigma_\ell)$ replaced by $(V_0, S_0)$.

To prepare the ensuing analysis for higher-order models, the energy associated with this problem is now studied. Similarly to what is done in Section 2, one obtains from the systems (25):

$$ \partial_t \mathcal{E}_0^{\text{vol}} + \mathcal{D}_0^{\text{int}} = 0, \qquad (27) $$

where the volume energy $\mathcal{E}_0^{\text{vol}}$ is now:

$$ \mathcal{E}_0^{\text{vol}} = \mathcal{E}_0^{\text{vol}^{-}} + \mathcal{E}_0^{\text{vol}^{+}} \quad \text{with} \quad \left\{ \begin{array}{l} \mathcal{E}_0^{\text{vol}^{-}} = \frac{1}{2} \int_{\Omega^{-}} \rho_{-} V_0^2 + \frac{S_0^2}{E_{-}} dx, \\ \mathcal{E}_0^{\text{vol}^{+}} = \frac{1}{2} \int_{\Omega^{+}} \rho_{0} V_0^2 + \frac{S_0^2}{E_{0}} dx, \end{array} \right. \qquad (28) $$

and the interface term is:

$$ \mathcal{D}_0^{\text{int}} = \llbracket V_0 S_0 \rrbracket_0 = \llbracket V_0 \rrbracket_0 \langle S_0 \rangle_0 + \langle V_0 \rangle_0 \llbracket S_0 \rrbracket_0. \qquad (29) $$

Then, a sufficient condition for the stability of the problem is this term to be the time derivative of a positive interface energy: if there exists $\mathcal{E}_0^{\text{int}} > 0$ such that $\mathcal{D}_0^{\text{int}} = \partial_t \mathcal{E}_0^{\text{int}}$, then (i) the total energy is conserved from (27) and (ii) the volume term $\mathcal{E}_0^{\text{vol}}$ is bounded:

$$ \partial_t (\mathcal{E}_0^{\text{vol}} + \mathcal{E}_0^{\text{int}}) = 0 \quad \Rightarrow \quad 0 \le \mathcal{E}_0^{\text{vol}}(t) = (\mathcal{E}_0^{\text{vol}}(0) + \mathcal{E}_0^{\text{int}}(0)) - \mathcal{E}_0^{\text{int}}(t) \le \mathcal{E}_0^{\text{vol}}(0), $$

because the initial conditions on $(V_0, S_0)$ ensure $\mathcal{E}_0^{\text{int}}(0) = 0$ and $\mathcal{E}_0^{\text{int}} \ge 0$ by assumption. This property then ensures that the fields $(V_0, S_0)$ remain bounded in time.

In the leading-order case, the interface term $\mathcal{D}_0^{\text{int}}$ is identically canceled by the transmission conditions (26), and therefore from (27) the energy $\mathcal{E}_0^{\text{vol}}$ is conserved.

## 3.3. First-order homogenized model: correcting the transmission conditions

First-order approximations from (3) are now considered:

$$ v_{\ell}(x, t) = V_{1}(x, t) + \ell P_{1}\left(\frac{x}{\ell}\right) \partial_{x} V_{1}(x, t) + o(\ell), \quad \sigma_{\ell}(x, t) = S_{1}(x, t) + \ell Q_{1}\left(\frac{x}{\ell}\right) \partial_{x} S_{1}(x, t) + o(\ell), \qquad (30) $$
where $(P_1, Q_1) = (0, 0)$ in the homogeneous domain $\Omega^-$. As seen in Section 2.2, the first-order macroscopic fields $(V_1, S_1)$ satisfy the same non-dispersive systems (25) as at leading order in each of the two half-spaces. On the other hand, replacing the fields $(v_\ell, \sigma_\ell)$ by their approximations (30) above in the transmission conditions (23) and neglecting the remainders, as proposed by [15], one obtains imperfect non-symmetrical conditions on the jump of the macroscopic fields:

$$ \begin{cases} \llbracket V_1 \rrbracket_0 = -\ell P_1(0) (\partial_x V_1)^+, \\ \llbracket S_1 \rrbracket_0 = -\ell Q_1(0) (\partial_x S_1)^+, \end{cases} \qquad (31) $$

i.e. a first-order correction compared to the conditions at the leading order (26).

A similar energy analysis than at the leading order results in

$$ \partial_t \mathcal{E}_1^{\text{vol}} + \mathcal{D}_1^{\text{int}} = 0, \qquad (32) $$

where the total volume energy $\mathcal{E}_1^{\text{vol}}$ is the same as $\mathcal{E}_0^{\text{vol}}$ given by (28) with $(V_0, S_0)$ replaced by $(V_1, S_1)$. The same analogy holds for the interface term $\mathcal{D}_1^{\text{int}}$ compared with $\mathcal{D}_0^{\text{int}}$ in (29). Given the conditions (31), one has:

$$ \mathcal{D}_1^{\text{int}} = \llbracket V_1 \rrbracket_0 \langle S_1 \rangle_0 + \langle V_1 \rangle_0 \llbracket S_1 \rrbracket_0 = -\ell [P_1(0)(\partial_x V_1)^+ \langle S_1 \rangle_0 + Q_1(0)(\partial_x S_1)^+ \langle V_1 \rangle_0]. \qquad (33) $$

There is no guarantee that $\mathcal{D}_1^{\text{int}}$ is the derivative of a positive interface energy. The negative case may result in an ill-posed problem. The following paragraphs explain how to deal with this issue, through symmetrization and stabilization of the jump conditions.

**Spring-mass transmission conditions.** To symmetrize the relations (31), the limit values e.g. $\partial_x V_1^+$ need to be written in terms of mean values across the interfaces, e.g. $\langle \partial_t S_1 \rangle_0$. The process introduces additional approximations, up to residuals that should be at least of second order (i.e. in $o(\ell)$) to preserve the overall first-order approximation. The main tools to do so are the bulk equations that link the time and space derivatives of $(V_1, S_1)$ in $\Omega^+$:

$$ \partial_x V_1 = E_0^{-1} \partial_t S_1 \quad \text{and} \quad \partial_x S_1 = \rho_0 \partial_t V_1 \quad \text{in } \Omega^+. \qquad (34) $$

Then, using also the jump-mean relations (22), the following approximations are found at the leading order:

$$ \begin{aligned} (\partial_x V_1)^+ &= E_0^{-1} \partial_t S_1^+, & (\partial_x S_1)^+ &= \rho_0 \partial_t V_1^+, \\ &= E_0^{-1} \partial_t \left[ \langle S_1 \rangle_0 + \frac{1}{2} \llbracket S_1 \rrbracket_0 \right], & &= \rho_0 \partial_t \left[ \langle V_1 \rangle_0 + \frac{1}{2} \llbracket V_1 \rrbracket_0 \right], \quad (35) \\ &= E_0^{-1} \partial_t \langle S_1 \rangle_0 + o(1), & &= \rho_0 \partial_t \langle V_1 \rangle_0 + o(1), \end{aligned} $$

where the jump conditions (31) were used to neglect the first-order jumps $\llbracket S_1 \rrbracket_0$ and $\llbracket V_1 \rrbracket_0$. Combining these approximations with (31) and neglecting second-order terms, one finally obtains:

$$ \begin{cases} \llbracket V_1 \rrbracket_0 = -\ell E_0^{-1} P_1(0) \partial_t \langle S_1 \rangle_0, \\ \llbracket S_1 \rrbracket_0 = -\ell \rho_0 Q_1(0) \partial_t \langle V_1 \rangle_0. \end{cases} \qquad (36) $$

This form is better known as spring-mass transmission conditions, that may also model e.g. a thin layer of elastic material between the considered media [30, App. 1]. A new interface term $\mathcal{D}_{1,\text{sm}}^{\text{int}}$ then replaces $\mathcal{D}_1^{\text{int}}$ in (32)-(33), conveniently written as the derivative of an interface energy $\mathcal{E}_{1,\text{sm}}^{\text{int}}$:

$$ \mathcal{D}_{1,\text{sm}}^{\text{int}} = \partial_t \mathcal{E}_{1,\text{sm}}^{\text{int}} \quad \text{with:} \quad \mathcal{E}_{1,\text{sm}}^{\text{int}} = -\frac{\ell}{2} [E_0^{-1} P_1(0) \langle S_1 \rangle_0^2 + \rho_0 Q_1(0) \langle V_1 \rangle_0^2]. \qquad (37) $$

However, there is still no guarantee that this interface energy remains positive for all time, except in the particular cases where $P_1(0) \le 0$ and $Q_1(0) \le 0$. It is why a last modelling step is necessary, as addressed now.
**Stabilization introducing an enlarged interphase.** To obtain centered transmission conditions that are associated to a positive interface energy, a thick interface $I^d = [-d\ell, d\ell]$ of thickness $2d\ell$, centered on the point $x = 0$, is introduced, and the transmission conditions are written across this interface. This method has been first used for microstructured interface homogenization [17], but also for a 2D transmission problem similar to the present case in [20, 40].

Using Taylor expansions of the macroscopic fields in the homogeneous and homogenized media, e.g. $V_1(0^-) = V_1(-d\ell) + d\ell \partial_x V_1(-d\ell) + o(\ell)$ and keeping only first-order terms in $\ell$, the conditions (31) become:

$$ \left\{ \begin{array}{l} V_1(-d\ell) + d\ell \partial_x V_1(-d\ell) = V_1(d\ell) + \ell(P_1(0) - d) \partial_x V_1(d\ell), \\ S_1(-d\ell) + d\ell \partial_x S_1(-d\ell) = S_1(d\ell) + \ell(Q_1(0) - d) \partial_x S_1(d\ell). \end{array} \right. \qquad (38) $$

Then, using the same first-order approximations (35) than above, one obtains equivalent (up to $o(\ell)$) spring-mass conditions:

$$ \left\{ \begin{array}{l} \llbracket V_1 \rrbracket_d = \ell A_1 \partial_t \langle S_1 \rangle_d, \qquad A_1 = -P_1(0)E_0^{-1} + d(E_0^{-1} + E_0^{-1}), \\ \llbracket S_1 \rrbracket_d = \ell B_1 \partial_t \langle V_1 \rangle_d, \qquad B_1 = -Q_1(0)\rho_0 + d(\rho_0 + \rho_0), \end{array} \right. \qquad (39) $$

where $\llbracket \cdot \rrbracket_d$ and $\langle \cdot \rangle_d$ are the jump and mean values of fields across the interface $I^d$. Similarly to (37), the interface term in the energy is then written:

$$ \mathcal{D}_{1,d}^{\text{int}} = \partial_t \mathcal{E}_{1,d}^{\text{int}} \quad \text{with:} \quad \mathcal{E}_{1,d}^{\text{int}} = \frac{\ell}{2} [A_1 \langle S_1 \rangle_d^2 + B_1 \langle V_1 \rangle_d^2]. $$

A sufficient condition for the interface energy $\mathcal{E}_{1,d}^{\text{int}}$ to be positive is that the two coefficients $(A_1, B_1)$ are positive. This is achieved by choosing the interface parameter $d$ such that:

$$ d \ge d_{\min} := \max \left( \frac{E_{-}P_{1}(0)}{E_{-} + E_{0}}, \frac{\rho_{0}Q_{1}(0)}{\rho_{-} + \rho_{0}}, 0 \right). \qquad (40) $$

As already noticed, when $P_1(0) \le 0$ and $Q_1(0) \le 0$, there is no need for a thick interface: one can choose $d = d_{\min} = 0$ and the conditions (36) are retrieved. In the other cases, choosing $d = d_{\min}$ will cancel one of the factors $(A_1, B_1)$ in the relations (39), i.e. one has either $\llbracket V_1 \rrbracket_d = 0$ or $\llbracket S_1 \rrbracket_d = 0$.

**Remark 5.** Instead of introducing a thick interface, another way to design stable interface conditions from (36) is to slightly modify the homogenized model by choosing non-zero means to the cell functions $(P_1, Q_1)$ to ensure $P_1(0) \le 0$ and $Q_1(0) \le 0$. This method was introduced in [1] to design stable wave models in unbounded media, and applied in [8] to stabilize enriched boundary conditions at the edge of a half-plane. It was implemented by the authors of the present paper, and is indeed a good alternative to treat one interface. But since it requires a global change of the model (in the whole microstructured domain), it cannot be easily extended of the case of a slab treated below in Section 3.5, or any case where more than one interface needs to be addressed. The thick interface method, that introduces a local approximation, was therefore prefered.

**Remark 6.** When $E = E_0$ is constant, i.e. when $\rho$ is the only varying parameter, one has $P_1 = 0$ (see Appendix A.1) and therefore there is no first-order correction on the velocity jump in the spring-mass conditions (36), but still a correction on the stress. This was shown already in [32] for electromagnetic waves (u being the electric field and $\rho$ being the permittivity). However, depending on the sign of $Q_1(0)$, stabilizing the interface condition may still be necessary, and artificial corrections on the velocity jump may appear e.g. if the “thick” interface method is chosen, resulting in conditions (39).

### 3.4. Total “hybrid” model incorporating second-order dispersion

Finally, the second-order model studied in Section 2 is implemented in the microstructured domain to account for dispersion. Using the four fields $(w, \sigma, \varphi, r)$ of the stress-gradient system in $\Omega^+$, while $(v_\ell, \sigma_\ell) \approx (w, \sigma)$ in
$\Omega^-$, these fields satisfy:

$$ \Omega^{-} : \left\{ \begin{array}{l} \partial_t w - \frac{1}{\rho_-} \partial_x \sigma = 0, \\ \partial_t \sigma - E_- \partial_x w = 0, \end{array} \right. \qquad \Omega^{+} : \left\{ \begin{array}{l} \partial_t w - \frac{a}{\rho_0} \partial_x \sigma = -\frac{a-1}{\rho_0} r, \\ \partial_t \sigma - E_0 \partial_x w = 0, \\ \partial_t \varphi - \frac{a-1}{\rho_0} \partial_x \sigma = -\frac{a-1}{\rho_0} r, \\ \partial_t r = \frac{E_0}{\ell^2 \beta} \varphi. \end{array} \right. \qquad (41) $$

At this stage, an important remark is that the two auxiliary fields $(\varphi, r)$ satisfy ordinary differential equations in time, and therefore no boundary condition are needed for these fields (but initial conditions must be added). The systems (41) must therefore be complemented by transmission conditions on $(w, \sigma)$ only, as in the leading-order and first-order cases, and by initial conditions for all fields. Since we choose initial conditions only supported in the left domain $\Omega^-$, we set:

$$ w(x, 0) = v_{*}(x), \quad \sigma(x, 0) = \sigma_{*}(x), \quad \varphi(x, 0) = 0 \quad \text{and} \quad r(x, 0) = 0. \qquad (42) $$

Then, combining the results on energies of the stress-gradient system collected in Section 2.3, and the additional term coming from the interface as specified above, one obtains:

$$ \partial_t \mathcal{E}^{\text{vol}} + \mathcal{D}^{\text{int}} = 0, $$

where the total energy $\mathcal{E}^{\text{vol}}$ is now:

$$ \mathcal{E}^{\text{vol}} = \mathcal{E}^{\text{vol}^{-}} + \mathcal{E}^{\text{vol}^{+}} \quad \text{with} \quad \left\{ \begin{array}{l} \mathcal{E}^{\text{vol}^{-}} = \frac{1}{2} \int_{\Omega^{-}} \rho_{-} w^2 + \frac{\sigma^2}{E_{-}} dx, \\ \mathcal{E}^{\text{vol}^{+}} = \frac{1}{2} \int_{\Omega^{+}} \rho_{0} \left( w^2 - 2w\varphi + \frac{a}{a-1}\varphi^2 \right) + \frac{1}{E_{0}} \left( \sigma^2 + \ell^2 \beta r^2 \right) dx, \end{array} \right. \qquad (43) $$

and the interface term $\mathcal{D}^{\text{int}}$ is again given by:

$$ \mathcal{D}^{\text{int}} = \llbracket w\sigma \rrbracket_0 = \llbracket w \rrbracket_0 \langle \sigma \rangle_0 + \langle w \rangle_0 \llbracket \sigma \rrbracket_0. $$

This term is formally identical to the ones that appear for the first-order model (33).

At this stage, it would be natural to use the second-order approximations in (3) to design transmission conditions for $(w, \sigma)$. However, for technical reasons discussed in Remark 7 below, only *first-order* transmission conditions are designed to obtain a “hybrid” model, using the terminology introduced in [32].

As in the previous section, the chosen total fields approximations are:

$$ v_{\ell}(x, t) = w(x, t) + \ell P_{1}\left(\frac{x}{\ell}\right) \partial_{x} w(x, t) + o(\ell), \quad \sigma_{\ell}(x, t) = \sigma(x, t) + \ell Q_{1}\left(\frac{x}{\ell}\right) \partial_{x} \sigma(x, t) + o(\ell), \qquad (44) $$

and the same analysis can be performed. Indeed, the links between the fields derivatives given by (34) can be replaced by:

$$ \partial_x w = E_0^{-1} \partial_t \sigma \quad \text{and} \quad \partial_x \sigma = \rho_0 \partial_t w + o(\ell) \quad \text{in } \Omega^+. \qquad (45) $$

The second relation in (45) is obtained by taking the difference between the first and third equation of (41) in $\Omega^+$, and by using $\varphi = o(\ell)$. The second-order remainder $o(\ell)$ in (45) has no influence on the approximations at the leading order (35). Consequently, the same first-order spring-mass transmission conditions (39), written
across a thick interface, are applied to $(w, \sigma)$:

$$ \left\{ \begin{array}{ll} \llbracket w \rrbracket_d = \ell A_1 \partial_t \langle \sigma \rangle_d, & A_1 = -P_1(0)E_0^{-1} + d(E_0^{-1} + E_0^{-1}), \\ \llbracket \sigma \rrbracket_d = \ell B_1 \partial_t \langle w \rangle_d, & B_1 = -Q_1(0)\rho_0 + d(\rho_- + \rho_0). \end{array} \right. \qquad (46) $$

The various asymptotic formal expansions are still valid and indicate that the transmission conditions (46) should result in an overall first-order approximation. In this way, the overall problem including (i) the systems (41), (ii) the initial conditions (42) and (iii) the transmission conditions (46), is proven to be stable using the same arguments than in the previous section. It is summarized in the next proposition.

**Proposition 3.** We consider as "total model" the system (41) with the symmetrized jump conditions (46) across a thick interface, where the interface parameter $d$ satisfies (40). Associated with this system is a constant energy over time $\mathcal{E} = \mathcal{E}^{\text{vol}} + \mathcal{E}_d^{\text{int}}$, where $\mathcal{E}^{\text{vol}} \ge 0$ is given by (43), and $\mathcal{E}_d^{\text{int}}$ is given by

$$ \mathcal{E}_d^{\text{int}} = \frac{\ell}{2} [A_1 \langle \sigma \rangle_d^2 + B_1 \langle w \rangle_d^2] \ge 0. $$

**Remark 7.** Imposing transmission conditions to the second-order approximation (16) of the total fields, as done in the time-harmonic case in [15], would lead to

$$ \begin{aligned} w^{-} &= w^{+} + \ell P_{1}(0)(\partial_{x}w)^{+} + \ell^{2} P_{2}(0)(\partial_{xx}w)^{+} - \varphi^{+}, \\ \sigma^{-} &= \sigma^{+} + \ell Q_{1}(0)(\partial_{x}\sigma)^{+} + \ell^{2} Q_{2}(0)(\partial_{xx}\sigma)^{+} - \ell^{2}\beta(\partial_{x}r)^{+}, \end{aligned} $$

where $w^{\pm} = w(0^{\pm}, t)$ and similarly for $(\sigma, \varphi, r)$ and their derivatives. The same ideas than in Section 3.3 can be applied to obtain symmetrized jump conditions on a thick interface of the form:

$$ \left\{ \begin{array}{l} \llbracket w \rrbracket_d = \ell A_1 \partial_t \langle \sigma \rangle_d + \ell^2 A_2 \partial_{tt} \langle w \rangle_d + \varphi(d\ell), \\ \llbracket \sigma \rrbracket_d = \ell B_1 \partial_t \langle w \rangle_d + \ell^2 (B_2 \partial_{tt} \langle \sigma \rangle_d + \beta(\partial_x r)(d\ell)), \end{array} \right. $$

with $(A_1, B_1)$ given by (46) and new coefficients $(A_2, B_2)$, generalizing (46). However, these second-order terms feature (i) (non-classical) second-order time derivatives of $(w, \sigma)$ and (ii) the traces of auxiliary fields $(\varphi, \partial_x r)$, and our attempts to prove the stability of the resulting model were unsuccessful so far. Moreover, the "hybrid" model described in Proposition 3 already provides a good compromise between implementation easiness (classical spring-mass transmission conditions are used) and both qualitative and quantitative performances, as illustrated in Section 5 below.

### 3.5. Microstructured slab bounded by two homogeneous domains

To end this section, one extends the previous model to the case of a microstructured medium bounded at both extremities by homogeneous domains. It is representative of many real experiments of wave propagation. The microstructured domain is $\Omega^L = ]0, L[$, while $\Omega^+ = \{x > L\}$ now denotes the right homogeneous domain, characterized by coefficients $(\rho_+, E_+)$. Additional transmission conditions must be designed for the second interface at $x = L$. The continuity of the first-order approximations (44) is written:

$$ \left\{ \begin{array}{l} w(L^{-}) + \ell P_{1}(y_{L}) \partial_{x} w(L^{-}) = w(L^{+}), \\ \sigma(L^{-}) + \ell Q_{1}(y_{L}) \partial_{x} \sigma(L^{-}) = \sigma(L^{+}), \end{array} \right. $$

where $y_L = L/\ell$. Then, a similar analysis as in Section 3.3 leads to a reformulation of these conditions involving the jumps $\llbracket \cdot \rrbracket_{d'}$ and means $\langle \cdot \rangle_{d'}$ across a thick interface $[L - d'\ell, L + d'\ell]$ as:

$$ \left\{ \begin{array}{ll} \llbracket w \rrbracket_{d'} = \ell A'_1 \partial_t \langle \sigma \rangle_{d'} & A'_1 = P_1(y_L)E_0^{-1} + d'(E_0^{-1} + E_+^{-1}), \\ \llbracket \sigma \rrbracket_{d'} = \ell B'_1 \partial_t \langle w \rangle_{d'} & B'_1 = Q_1(y_L)\rho_0 + d'(\rho_0 + \rho_+), \end{array} \right. \qquad (47) $$
where the two coefficients ($A'_1, B'_1$) must be positive for these conditions to be associated with a positive interface energy. The thickness $d'$ of this second interface must be chosen so that:

$$d' \ge d'_{\min} := \max \left( -\frac{E_+ P_1(y_L)}{E_+ + E_0}, -\frac{\rho_0 Q_1(y_L)}{\rho_+ + \rho_0}, 0 \right), \quad (48)$$

a condition similar to its counterpart (40) except for the signs.

## 4. CORRECTORS FOR A DIRAC SOURCE TERM

So far, no source term was considered in the evolution equations. The initial conditions were supported in the homogeneous domain, thus avoiding to adapt them to the microstructure [27]. To complete the modelling framework, the case of a single force source point immersed in the microstructured medium is now addressed, using the tools previously developed for interfaces.

### 4.1. Source terms and jump conditions

A point force at point $x_s$ is modeled by a Dirac source term in the initial wave equation:

$$\rho_{\ell}(x) \partial_{tt} u_{\ell}(x, t) - \partial_{x} (E_{\ell}(x) \partial_{x} u_{\ell}(x, t)) = \rho_{\ell}(x) g(t) \delta_{s}(x), \quad \delta_{s}(x) := \delta(x - x_{s}), \quad (49)$$

where $\delta$ is the Dirac distribution and $g$ the time-dependent amplitude of the source. To simplify the analysis, and having in mind the microstructures with fixed periodicity length $\ell$ considered in the numerical examples, we furthermore assume that the point source position is fixed at $y_s \in [0, 1]$ inside a unit cell, and we note $\rho_s = \rho(y_s)$ the density at this position, that does not depend on $\ell$.

The above equation (49) is to be understood in the weak sense, i.e. :

$$\int_{\Omega} (\rho_{\ell} \partial_{tt} u_{\ell} v + E_{\ell} \partial_{x} u_{\ell} \partial_{x} v) dx = [E_{\ell}(\partial_{x} u_{\ell})v](x^{+}) - [E_{\ell}(\partial_{x} u_{\ell})v](x^{-}) + \rho_{s} v(x_{s})g, \quad \forall v \in H(\Omega), \quad (50)$$

with a suitable choice of domain $\Omega = ]x^{-}, x^{+}[$ and functional space $H(\Omega)$. As in the previous sections, space and time dependency are omitted in (50) and below for compactness.

Such a source point has been addressed by [12], and even extended to double-couple source (modelling ruptures of seismic faults) and to higher dimensions (for non-periodic homogenization) in [13]. However, the provided expression, relying on energy conservation arguments, only treat first-order homogenization. The present proposal exploits the fact that in 1D, a Dirac distribution yields a jump condition on the solution. As a consequence, the formalism developed in the previous sections can be reused here to complete the two-scale formal expansion.

**Remark 8.** This approach does not hold in higher dimensions, for which an analysis of the correction for a source point should be incorporated into the two-scale asymptotic expansion, see [13]. However, it could be applied to line sources in 2D or area sources in 3D, modelling e.g. a distribution of charges in electromagnetism, arrays of microphones in acoustics etc.

**Replacing the source point by a jump condition.** In 1D, a Dirac source term entails a discontinuity in the stress $\sigma_{\ell} = E_{\ell} \partial_{x} u_{\ell}$ (accordingly with the physical interpretation of a Dirac source term as a point force). Indeed, integrating the wave equation (49) between $x_s - \xi$ and $x_s + \xi$, one obtains:

$$\sigma_{\ell}(x_s + \xi) - \sigma_{\ell}(x_s - \xi) - \int_{x_s - \xi}^{x_s + \xi} \rho_{\ell} \partial_{tt} u_{\ell} \, dx = - \rho_s g,$$
and then letting $\xi \to 0$ gives the announced result. Accordingly, the wave equation (49) may be rewritten as the system:

$$ \left\{ \begin{array}{ll} \partial_t v_\ell - \frac{1}{\rho_\ell} \partial_x \sigma_\ell = 0 & x < x_s \text{ and } x > x_s, \\ \partial_t \sigma_\ell - E_\ell \partial_x v_\ell = 0 & x < x_s \text{ and } x > x_s, \\ [v_\ell]_s = 0 & x = x_s, \\ [\sigma_\ell]_s = -\rho_s g & x = x_s, \end{array} \right. \qquad (51) $$

where $[\cdot]_s$ denotes the jump at source location $x_s$. This problem can be homogenized as the previously seen transmission problem, except that the same periodic (or homogenized) medium occupies both half-spaces, and the successive corrections will account for the discontinuity of $\sigma_\ell$.

## 4.2. Families of models

**Leading-order homogenization.** At the leading order, as in Section 3.2, $E_\ell$ and $\rho_\ell$ are replaced by the constant $E_0$ and $\rho_0$ in the bulk equations of system (51), and the source term is unchanged (still featuring the value $\rho_s$).

**First-order homogenization.** At first order, as in Section 3.3, the bulk equations satisfied by the macroscopic fields denoted ($V_1, S_1$) are the same as at leading order, and the transmission conditions are applied to the first-order approximations of the full fields (30), leading to:

$$ x < x_s \text{ and } x > x_s : \left\{ \begin{array}{ll} \partial_t V_1 - \frac{1}{\rho_0} \partial_x S_1 = 0 & \text{and} \\ \partial_t S_1 - E_0 \partial_x V_1 = 0 & \end{array} \right. \qquad \left\{ \begin{array}{l} [V_1]_s + \ell P_1(y_s) [\partial_x V_1]_s = 0, \\ [S_1]_s + \ell Q_1(y_s) [\partial_x S_1]_s = -\rho_s g, \end{array} \right. \qquad (52) $$

where $y_s = x_s/\ell$ is the source position inside the cell that contains $x_s$. Then, the bulk equations provide approximations at the leading order of the jumps of spatial derivatives:

$$ [\partial_x V_1]_s = \frac{1}{E_0} \partial_t [S_1]_s = -\frac{\rho_s}{E_0} \partial_t g + o(1) \quad \text{and} \quad [\partial_x S_1]_s = \rho_0 \partial_t [V_1]_s = o(1). \qquad (53) $$

Neglecting the $o(\ell)$ terms in (52), one obtains first-order conditions:

$$ \left\{ \begin{array}{l} [V_1]_s = \ell P_1(y_s) \frac{\rho_s}{E_0} \partial_t g, \\ [S_1]_s = -\rho_s g, \end{array} \right. \qquad (54) $$

i.e. the system in (52) features an additional inhomogeneous jump condition on the macroscopic velocity. The resulting system may finally be written in the equivalent form:

$$ \left\{ \begin{array}{l} \partial_t V_1 - \frac{1}{\rho_0} \partial_x S_1 = \frac{\rho_s}{\rho_0} \delta_s g, \\ \partial_t S_1 - E_0 \partial_x V_1 = -\ell \rho_s P_1(y_s) \partial_t g \delta_s. \end{array} \right. \qquad (55) $$

This expression was also obtained in [12, Section 2.4] using an energy conservation argument.

**Application to the stress-gradient system.** Finally, the stress-gradient system (20) satisfied by the macroscopic fields $(w, \sigma, \varphi, r)$ is considered. As in Section 3.4, only first-order corrections are applied to the source term for simplicity. Hence, the analysis performed above for these corrections remains valid: the main difference is that the stress-velocity relation $\partial_x \sigma = \rho_0 \partial_t w + o(\ell)$ is a first-order approximation rather than an equality, but
it does not change the approximations at the leading order (53). Therefore jump conditions (54) are applied to $(w, \sigma)$.

Coming back to Dirac notations, the stress spatial derivative $\partial_x\sigma$ is involved in two equations of the system, and hence the source term must be accordingly distributed:

$$ \left\{ \begin{aligned} \partial_t w - \frac{a}{\rho_0} \partial_x \sigma &= -\frac{a-1}{\rho_0} r + \frac{a \rho_s}{\rho_0} g \delta_s, \\ \partial_t \sigma - E_0 \partial_x w &= -\ell \rho_s P_1(y_s) \partial_t g \delta_s, \\ \partial_t \varphi - \frac{a-1}{\rho_0} \partial_x \sigma &= -\frac{a-1}{\rho_0} r + \frac{(a-1) \rho_s}{\rho_0} g \delta_s, \\ \partial_t r &= \frac{E_0}{\ell^2 \beta} \varphi. \end{aligned} \right. \qquad (56) $$

As in Section 3.4, this is a hybrid model combining first-order corrections for the source term and second-order corrections in the wave equation to account for dispersive effects. If transmission correctors were derived at the second order, see Remark 7, applying them to the problem (51) would certainly enable to obtain a complete second-order model.

# 5. NUMERICAL EXPERIMENTS

The proposed models are now applied to simulate the wave propagation into an example bilaminate material. The material properties, initial conditions and source terms are first given (Section 5.1), the numerical implementation is briefly described (Section 5.2), and two test-cases are addressed: (i) a wave transmitted from a homogeneous medium towards the bilaminate (Section 5.3) and (ii) a wave emitted from a source point embedded in a slab of bilaminate, surrounded by homogeneous media (Section 5.4).

## 5.1. Microstructured medium, initial conditions and source modelling

The numerical illustrations are provided for bilaminate materials, for which analytical expressions of the homogenized coefficients and cell problems are available, see Appendix A.3. Furthermore, an analytical formula for $(\beta_m, \beta_t)$ (or, here, a) was proposed in [15] to achieve a better fit of the exact dispersion curve for low frequencies, see (68). This specific choice of $(mt)$ model was also discussed in [36] and found to compete reasonably with numerically optimized models that achieve an overall fit of this curve on the whole Brillouin zone. Finally, this choice is numerically found to be compatible with the stability conditions $\{\beta_m < 0, \beta_t > 0\}$ established by (18).

<table><caption>**Figure 1.** Unit cell of a bilaminate with two phases (A, B) and phase ratio $\alpha$.</caption><tr><td>E<sub>A</sub></td><td>E<sub>B</sub> = &gamma;<sub>E</sub>E<sub>A</sub></td></tr><tr><td>&rho;<sub>A</sub></td><td>&rho;<sub>B</sub> = &gamma;<sub>&rho;</sub>&rho;<sub>A</sub></td></tr><tr><td>&alpha;</td><td>1 &minus; &alpha;</td></tr></table>

The bilaminate investigated here has a periodicity length $\ell = 20$ m, with a phase ratio $\alpha = 0.25$ (figure 1). The physical parameters are $\rho_A = 1000$ kg/m$^3$, $E_A = 10^9$ Pa, and $\rho_B = 1500$ kg/m$^3$, $E_B = 6 \times 10^9$ Pa. The associated homogenized parameters are $\rho_0 = 1375$ kg/m$^3$, $E_0 = 2.66 \times 10^9$ Pa and $\beta = 0.01959$, see Appendix A.3. The optimized parameters of the $(mt)$ model (68) are $\beta_m = -0.107071$ and $\beta_t = 0.087481$.
<figure><img src="image_2.png" /><figcaption>Figure 2. Time evolution of the source (58) at $f_c = 3$ Hz (a). Spectrum of the source for various central frequencies $f_c$, as a function of the asymptotic parameter $\epsilon_c$ (b).</figcaption></figure>

In the following examples, two types of forcing are considered: (i) a Cauchy problem (24) with an initial right-going wave

$$ \begin{pmatrix} v_* \\ \sigma_* \end{pmatrix} = \begin{pmatrix} -\frac{1}{c_-} \\ \rho_- \end{pmatrix} g \left( t + t_0 - \frac{x}{c_-} \right), \qquad (57) $$

or (ii) null initial conditions with a time-dependent forcing at a Dirac source point (Section 4.2). Both configurations involve the source function $g(t)$, chosen as the following combination of sinusoids with bounded support:

$$ g(t) = \begin{cases} G \sum_{m=1}^{4} a_m \sin(b_m \omega_c t) & \text{if } 0 < t < \frac{1}{f_c}, \\ 0 & \text{otherwise,} \end{cases} \qquad (58) $$

where $b_m = 2^{m-1}$, the coefficients $a_m$ are $a_1 = 1$, $a_2 = -21/32$, $a_3 = 63/768$, $a_4 = -1/512$, and the amplitude factor is $G \equiv 1$, unless otherwise stated. It entails that $g$ is a smooth function ($g \in C^6([0, +\infty])$), see the display in Figure 2(a).

Moreover, $g(t)$ is a wide-band signal with a central frequency $f_c = \omega_c/2\pi$. As announced in Remark 1, this frequency is used to define a dimensionless parameter characteristic of the “low-frequency” regime underlying the formal asymptotic expansions used throughout the paper:

$$ \epsilon_c = \ell f_c/c_0. $$

We denote $\hat{g}(\omega)$ the Fourier transform in time of $g(t)$, where $\omega = 2\pi f$ and $f$ is the frequency. Figure 2(b) displays the normalized modulus $|\hat{g}|(\omega)$ in terms of $\epsilon = \ell f/c_0 = \ell \omega/(2\pi c_0)$, for four values of $f_c$. Note that even for the smallest central frequencies, non-negligible values of $\epsilon$ are reached by the signal.

## 5.2. Numerical set-up

In the following experiments, comparisons are made between full-field simulations in the microstructured configurations and simulations involving their homogenized counterparts at different orders. All the numerical simulations are done on a uniform grid, with space and time discretization parameters $\Delta x$ and $\Delta t$. The
discretization is fine enough to guarantee that the numerical artifacts are negligible. A much coarser grid can be used for simulations in homogenised media than for simulations in microstructured media (this is one of the great interests of homogenisation!). Nevertheless, as this is not the object of the present work, we use here the same grid $Δx = 1$ m (to be compared with cell size $ℓ = 20$ m) for both types of configurations.

The reference simulations in the microstructured configurations are done using a fourth-order ADER scheme [29]. This explicit two-step finite-volume scheme has a stencil of width 2. In the homogenized configurations, we prefer to use the standard Lax-Wendroff scheme, of width 1. Doing so ensures that the stencil does not cross the whole enlarged interphase, and only one ghost value is required in the interphase [39]. Both ADER and Lax-Wendroff schemes are stable under the usual CFL condition $ζ = \max(c)Δt/Δx \le 1$; in practice, $ζ = 0.95$ is chosen.

In the case of second-order homogenization, the stress-gradient formulation yields a non-homogeneous first-order system (17). A splitting strategy is then followed, where one solves alternatively the homogeneous part by the Lax-Wendroff scheme, and the relaxation part exactly thanks to the exponential (19). The reader is referred to [7] for details.

The discretization of the transmission conditions relies on the Explicit Simplified Interface Method. Roughly speaking, the numerical scheme is modified at grid points surrounding the interfaces, based on the jump conditions (31) or (46). Once again, the reader is referred to [30] for details. Lastly, the numerical approximations of the fields are denoted as follows:

* $v_h$ for the microstructured velocity $v_ℓ$ (1);

* $v_0$ for the zeroth-order homogenized velocity (25) with the transmission conditions (26);

* $v_1$ for the first-order homogenized velocity (31) without $o(ℓ)$, endowed with the conditions (31);

* $v_2$ for the “total” homogenized velocity (44) without $o(ℓ)$, endowed with the conditions (46).

### 5.3. Transmission towards a microstructured half-space

A bounded domain $[-400, 1100]$ is considered. An interface at $x_0 = 350$ m separates the domain into two areas: (i) on the left, a homogeneous medium; (ii) on the right, the bilaminate described above. A right-going wave is initially supported in the homogeneous medium, with a time shift $t_0 = 0.18$ s in (57).

To highlight the importance of correctors for the transmission conditions, the physical parameters of the homogeneous medium are chosen as $E_- = E_0$ and $ρ_- = ρ_0$, to ensure impedance matching ($E_-ρ_- = E_0ρ_0$) between the left and right media. At leading order and first order, the whole domain has thus homogeneous physical parameters, and the (first-order) transmission condition is the only signature of the microstructure.

Choosing $d = d_{min}$ in (40), one finds $d = 0.017$, yielding an enlarged interphase width $2dℓ = 0.68$ m. In homogenized configurations, two interfaces are thus located at $x_0 \pm dℓ$, thus at 349.66 m and at 350.34 m. The interface parameters deduced from (38) are then $A_1 = 9.09 \times 10^{-11}$ Pa$^{-1}$ and $B_1 = 0$ Pa m$^{-2}$ s$^2$.

Figure 3 illustrates the low-frequency configuration: $f_c = 3$ Hz ($ε_c = 0.043$). The wave is shown at the initial time (a) and at $t = 0.5$ (b), where a reflected wave is observed in the Cauchy medium. For the model at leading order (Section 3.2), the medium appears as homogeneous: no reflected wave is generated (c), which clearly shows the inability of this naive approach to handle the boundary effects of a homogenized medium. Moreover, the small-scale gradients of the solution in the microstructured medium are not captured. The first-order model (Section 3.3) fixes these deficiencies: the reflected wave is captured, and the correctors handle the cell evolution of the transmitted wave (d) (note that at the scale of the figure, the two interfaces cannot be distinguished from each other). The same observations can be done for the total model (Section 3.4): at low frequency, the latter captures equally well the wave phenomena (e). A close-up on the wave in the transmitted medium shows that the mean field $w$ (13) captures the slow variations of the solution, whereas the correctors in (15) describe finely the gradients of the solution (f).

A higher value of the central frequency is considered in Figure 4: $f_c = 9$ Hz ($ε_c = 0.123$). Compared with the low-frequency case, the microstructured field $v_h$ exhibits a dispersive behavior. The model at leading order does not capture any phenomena: no reflected wave, no dispersion, no small-scale gradient (a). The first-order model captures two of the phenomena (reflected wave and gradients), but the dispersion in the right domain
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Wave in a homogeneous half-space (a-b, left) impacting a microstructured half-space (a-b, right) with impedance matching ($E_{-}\rho_{-} = E_{0}\rho_{0}$). The central frequency is $f_{c} = 3$ Hz ($\varepsilon_{c} = 0.043$).</figcaption></figure>
cannot be handled (b): this model is intrinsically non-dispersive. On the contrary, the total model captures all these features (c). A close-up in the right subdomain illustrates the quality of the correctors of the transmitted wave (d). The agreement is slightly worse on the main front of the reflected wave. We do not know how to explain this difference in behaviour, which is also observed in other simulations not presented here. A loss of accuracy is also observed behind the main front of the reflected and transmitted waves. These oscillations are associated with high-frequency components that are not captured by the dispersive model.

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Interface between a homogeneous half-space and a microstructured half-space. The central frequency is $f_c = 9$ Hz ($\varepsilon_c = 0.123$), and time is $t = 0.5$.</figcaption></figure>

Lastly, Figure 5 illustrates the high central frequency $f_c = 12$ Hz ($\varepsilon_c = 0.172$), where dispersion becomes the predominant feature. The limitations of models at the leading order and first order are exemplified (a,b). Even for this large $\varepsilon_c$, the total model captures finely the main fronts of the transmitted dispersive wave (c,d). The amplitude errors are larger on the reflected wave.
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Interface between a homogeneous half-space and a microstructured half-space. The central frequency is $f_c = 12$ Hz ($\epsilon_c = 0.172$), and time is $t = 0.5$.</figcaption></figure>

Other simulations have been done in cases without impedance matching (not shown here). Similar conclusions are obtained: (i) the model at leading order yields a large error in the reflected wave; (ii) the first-order model captures well the reflected wave and the small-scale variations of the solution, but it is restricted to low frequencies; (iii) the total model captures all the features, and maintains accuracy even for large values of $\epsilon_c$.

The accuracy of the effective models is then examined quantitatively. A source point and a receiver are respectively located at $x_s < x_0$ and $x_r > x_0$. We impose that these positions correspond to nodes of the grid to avoid discretization effects. Moreover, the relative position of the receiver with respect to the microstructure is kept constant. Finally, the distances $L$, the measurement time $T$ and the forcing amplitude $G$ must be kept constant in scaled variables (denoted by bars) $\bar{L} = kL$, $\bar{T} = \omega T$, and $\bar{G} = G/\omega^2$, with $k = \omega/c$. To satisfy
<figure><img src="image_6.png" /><figcaption>Figure 6. Interface between a homogeneous half-space and a microstructured half-space. Error measurements for the velocity (a) and stress (b) fields. Dashed lines indicate expected asymptotic slopes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7. Interfaces between a homogeneous half-space and a microstructured half-space. (a) Error measurements on velocity fields $v_1$ and $v_2$ for variable interface thickness and two frequencies for the configuration of the previous figures. The minimum thickness $2d_{min}l = 0.68m$ is indicated by the dashed vertical line. (b) Example of instability occurring for a left-going wave impacting a flat interface (dashed vertical line). In this example the homogeneous and homogenized media occupy the right and left half-spaces, respectively.</figcaption></figure>

these criteria, the indices $m$ and $n$ are introduced and one defines

$$x_s(m) = x_0 - \left(m + \frac{1}{2}\right) D,$$

$$x_r(m) = x_0 + \left(m + \frac{1}{2}\right) l,$$

$$f_c(m) = \frac{n + \frac{1}{2}}{m + \frac{1}{2}} f_c(n),$$

$$T(m) = \frac{m + \frac{1}{2}}{n + \frac{1}{2}} T(n),$$

$$G(m) = \left( \frac{n + \frac{1}{2}}{m + \frac{1}{2}} \right)^2 G(n).$$

The spatial scales are $D = 4$ m and $l = 20$ m. The index $n = 30$ is chosen, with $f_c(30) = 1$ Hz, $T(30) = 3$ s and $G(30) = 1$, which defines the measurement parameters according to the index $m$ alone. The following values of $m$ are considered: 1, 2, 3, 4, 5, 8, 10, 15 and 30. The relative error in norm $L^2(0, T)$ is measured at receiver
position $x_r$ by comparing simulations in microstructured and effective media, i.e. :

$$Err(v_n) = \frac{\|v_n(x_r, \cdot) - v_h(x_r, \cdot)\|_{L^2(0,T)}}{\|v_h(x_r, \cdot)\|_{L^2(0,T)}}.$$

Figure 6-(a) represents the error on velocity fields {$v_n$}$_{n=0,1,2}$ as a function of the asymptotic parameter $\epsilon_c$. The expected orders of convergence are asymptotically retrieved, in particular the “total” model remains of first order due to the choice of total field and associated transmission conditions. However, accounting for dispersion using the stress-gradient model brings a clear improvement: additionally to the qualitative gains observed in Figures 4 and 5, the prefactor in front of the error on $v_2$ is divided by approximately 5 compared to that of the error on $v_1$. With reference to Remark 2, the error on stress fields is also plotted in Figure 6-(b), and the same observations can be made, highlighting the identical treatment given to the velocity and stress fields.

We also examine the influence of larger values of the interface thickness on the accuracy for two central frequencies (in previous experiments, the minimal value $2d\ell = 2d_{min}\ell \approx 0.68$ m was chosen). To do so, we consider integer values of $d\ell$, from which are deduced the interface parameters $A_1$ et $A_2$ (39). Figure 7(a) shows the relative error as a function of the interphase thickness for $v_1$ and $v_2$. In both cases, the error grows strictly with $d\ell$. However, this growth is moderate, which allows to enlarge the interphase while keeping an almost constant accuracy. One can therefore choose a given thickness from numerical requirements (e.g. mesh size or needed “ghost” nodes inside the interface) without deteriorating the solution too much. Since we do not have such requirements with the chosen numerical methods, the minimum value $d = d_{min}$ is kept in the following numerical experiments.

Finally, to illustrate what may happen when the sufficient conditions of stability on the transmission condition coefficients are not satisfied, Figure 7(b) represents a left-going wave impacting a flat interface between a homogeneous and a homogenized media occupying the right and left half-spaces, respectively, i.e. the configuration described in Section 3.5. The interface is not thickened, i.e. $d' = 0$ in (47), which results in $A'_1 = -7.81 \times 10^{-11} < 0$ and $B'_1 = 46.875$. When the wave crosses the interface, an instability appears and quickly grows: here it reaches the same amplitude than the incident wave in a few milliseconds (while the characteristic times of propagation are in seconds in these numerical experiments).

## 5.4. Source point in a slab

This second example illustrates both the modelling of a Dirac source term and of a bounded laminated medium. A domain $[-300, 1700]$ is divided into three parts: a heterogeneous slab of bilaminate [350, 1090] is surrounded by two homogeneous Cauchy media with properties $\rho_{\pm} = 1000$ kg/m³ and $E_{\pm} = 2.25 \times 10^9$ Pa. Contrary to the previous example, there is no impedance matching between the homogeneous and homogenized media. A Dirac source point is located at the middle of the slab at $x_s = 720$ m. It is handled through the systems (55) for first-order homogenization, or (56) for second-order homogenization.

Taking again $d = d_{min}$ in (40) and $d' = d'_{min}$ in (48), the interface thicknesses are $2d\ell = 0.788$ m (left interface) and $2d'\ell = 3.8$ m (right interface). The interface parameters deduced from (38) and (47) are: $A_1 = 9.42 \times 10^{-11}$ Pa⁻¹ and $B_1 = 0$ Pa m⁻² s² (left interface) and $A'_1 = 0$ Pa⁻¹ and $B'_1 = 273.3$ Pa m⁻² s² (right interface).

Figure 8 shows snapshots of the microstructured solution $v_h$ at several times, for various central frequencies. For the sake of clarity, the internal interfaces of the bilaminate are not shown. As time increases, one observes the waves emitted by the source point and scattered by the external interfaces, leading to transmitted waves (in the homogeneous external medial) and reflected waves (in the internal microstructured medium). At low frequency $f_c = 3$ Hz (Figure 8-(a)), the dispersion is negligible. The only small scale variations are due to the internal microstructure. As the central frequency increases (b,c,d), the dispersive effects are becoming increasingly important. A slight asymmetry between left and right parts of the figures may be observed. This is due to the grid-induced asymmetry of the Dirac source point.

Figure 9 compare the microstructured velocity $v_h$ with the velocity $v_2$ of the total model (Section 3.4). Once again, the internal interfaces are not shown. The velocities of the models at the leading order and first order
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Waves emitted in a microstructured slab surrounded by two homogeneous half-spaces: microstructured velocity $v_h$ for several source frequency and times. For the clarity of the figure, the internal interfaces are not shown.</figcaption></figure>

are discarded, since they are insufficient to capture accurately dispersion. Figure 9-(a,b) illustrates the low-frequency forcing $f_c = 3$ Hz. The fields are shown at $t = 0.25$ s (a), where the emitted waves are still in the slab, and $t = 0.6$ s (b), where waves have been transmitted in the surrounding media. As in Figure 5, dispersive effects are not visible at this frequency. Figures 9-(c-d) illustrate the case of a higher frequency $f_c = 9$ Hz, where the dispersive effects become more visible in the microstructured model. The total model captures the first wavefronts and their transmission to the surrounding media. Figures 9-(e-f) display the high-frequency case $f_c = 12$ Hz. Some inaccuracies of the total model can be observed, especially inside the slab, but the qualitative agreement with the microstructured solution remains good, even for this large value of the parameter $\epsilon_c = 0.172$.

For completeness, the stress fields are also displayed in Figure 10 for $f_c = 9$ Hz corresponding to the intermediate configuration of Figure 9-(c,d), and for thickened interfaces: $d = d' = 4$ m, so that $A_1 = 2.28 \times 10^{-10}$
<figure><img src="image_9.png" /><figcaption>(c) $f_c = 9$ Hz ($\epsilon_c = 0.123$), $t = 0.25$ s</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) $f_c = 9$ Hz ($\epsilon_c = 0.123$), $t = 0.6$ s</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(e) $f_c = 12$ Hz ($\epsilon_c = 0.172$), $t = 0.25$ s</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(f) $f_c = 12$ Hz ($\epsilon_c = 0.172$), $t = 0.6$ s</figcaption></figure>

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

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

**Figure 9.** Wave emitted by a source point in a heterogeneous slab: comparisons between “total” velocity $v_2$ (Section 3.4) and the microstructured field $v_h$ at three central frequencies (rows), represented at $t = 0.25$ (left column) and $t = 0.6$ (right column).
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Wave emitted by a source point in a heterogeneous slab: comparisons between “total” stress $\sigma_2$ (Section 3.4) and the microstructured field $\sigma_h$ at central frequency $f_c = 9$ Hz (central row of Figure 9), represented at $t = 0.25$ (a) and $t = 0.6$ (b).</figcaption></figure>

$Pa^{-1}$, $B_1 = 503.12$ Pa (left interface), $A'_1 = 8.57 \times 10^{-11}$ Pa$^{-1}$ and $B'_1 = 521.87$ Pa (right interface). Once again, and despite these non-optimal values of $d$ and $d'$, excellent agreement is observed between the microstructured and “total” fields, especially for the first wavefronts of the fields transmitted into the left and right homogeneous domains.

Finally, to provide another representation of these approximations, a receiver put at $x_r = 340$ m records the field transmitted to the left homogeneous domain. Figure 11 shows the time history of the absolute errors $|v_h(x_r, t) - v_n(x_r, t)|$. The signal is null up to $t \approx 0.3$ s, which corresponds of the travel time from the source point to the receiver. A low frequency (a), the error of the model at leading order is much bigger than that of the first-order model and total model. It emphasizes the role of the first-order transmission conditions. At higher frequencies (b-c-d), the models at leading order and first order are unable to capture the dispersive effects, which yields similar errors. On the contrary, the benefit induced by the total model is clearly observed. Lastly, the maximal amplitude of the error for each model increases with the frequency.

## 6. CONCLUSION AND PERSPECTIVES

In this paper, a consistent “hybrid” model was proposed for transient waves in periodic media, combining second-order correctors of the wave equation to account for dispersion and first-order correctors for transmission conditions from or toward homogeneous domains. Based on a reformulation as a hyperbolic system, its well-posedness was proven, and its efficiency was established through numerical simulations.

Many follows-up come in mind. First, designing effective jump conditions of order 2 as discussed in Remark 7 would lead to a fully consistent model at second order. This could perhaps correct the observed discrepancy between errors on reflected waves and transmitted ones, whose reason is unknown. Moreover, a rigorous asymptotic error analysis of these “total” models, either hybrid or second-order, is still missing and is necessary to completely validate the proposed methodology.

Then, higher frequencies could be addressed with a similar formalism using the extension of double-scale approaches established by [16] and developed afterwards e.g. [24, 25]. For time-harmonic boundary problems in 1D, a similar approach relying on Bloch waves is found in [38], but a time-domain counterpart is yet to be proposed.
<figure><img src="image_16.png" /><figcaption>Figure 11. Waves emitted in a microstructured slab surrounded by two homogeneous half-spaces. Time evolution of the errors $|v_h(x_r, t) - v_n(x_r, t)|$ at receiver $x_r$, with the orders $n = 0, 1, 2$.</figcaption></figure>

Finally, going to higher dimensions to tackle many more real-life configurations is an important research direction. Relying on the work [8, 20] for boundary conditions, a hyperbolic reformulation of the wave equation could for instance enable effective simulations in homogenized media.

## APPENDIX A. CELL FUNCTIONS AND HOMOGENIZED COEFFICIENTS

This appendix gathers existing results on the cell functions, effective coefficients and particular closed-form solutions for bilaminates, and Section A.2 additionally provides the proof of Proposition 1. The notation $\langle f \rangle_Y$ is used for the mean value of a function $f$ on the periodicity cell $Y = ]0, 1[$:

$$ \langle f \rangle_Y = \int_0^1 f(y)dy $$
### A.1. Definitions of cell functions and homogenized coefficients

Cell functions $P_j$ and $Q_j$, associated respectively with displacement/velocity and stress correctors, see (3), satisfy the following static problems on $Y$:

$$ \partial_y Q_j = f_j, \quad Q_j = \frac{E}{E_0}(P_j + \partial_y P_{j+1}), \quad P_j \text{ is 1-periodic, } \langle P_j \rangle_Y = 0, \quad j \in \{0, 1, 2\}, $$

where $P_0 = 1$ and the source-terms $f_j$ are defined recursively as:

$$ f_0 = 0, \quad f_1 = \frac{\rho}{\rho_0} - Q_0 \quad \text{and} \quad f_2 = \frac{\rho}{\rho_0} P_1 - Q_1, $$

Homogenized coefficients $(\rho_j, E_j)$ are then defined as weighted means of their periodic counterparts, involving the cell functions:

$$ \rho_j = \langle \rho P_j \rangle_Y \quad \text{and} \quad E_j = \langle E(P_j + \partial_y P_{j+1}) \rangle_Y = E_0 \langle Q_j \rangle_Y, $$

and one can verify that these definitions lead to the expressions (5) for $(\rho_0, E_0)$.

**Remark 9.** With the definitions above, one has $P_0 = Q_0 = 1$ and therefore:

$$ \partial_y Q_1 = \frac{\rho}{\rho_0} - 1 \quad \text{and} \quad \partial_y P_1 = \frac{E_0}{E} - 1. \qquad (59) $$

Finally, the coefficient $\beta$ which is involved in the second-order models (7) aggregates the second-order contributions:

$$ \beta = \frac{E_2}{E_0} - \frac{\rho_2}{\rho_0}. \qquad (60) $$

### A.2. Reciprocity identities and properties of the coefficients

Additional relations between the effective coefficients can be obtained by using the weak form of the cell problems:

$$ \text{Find } P_j \in H_{\#}^1, \quad \int_0^1 E(\partial_y P_j)(\partial_y w)dy = F_j(w) \quad \forall w \in H_{\#}^1, \quad j \in \{1, 2, 3\}, \qquad (61) $$

with $H_{\#}^1 = \{w \in H^1([0, 1]), \langle w \rangle = 0, w \text{ is 1-periodic}\}$, and where:

$$ F_1(w) = \int_0^1 -E \partial_y w \ dy, $$

$$ F_2(w) = \int_0^1 \left( -EP_1 \partial_y w + E(1 + \partial_y P_1)w - E_0 \frac{\rho}{\rho_0} w \right) dy, $$

$$ F_3(w) = \int_0^1 \left( -EP_2 \partial_y w + E(P_1 + \partial_y P_2)w - E_0 \frac{\rho}{\rho_0} P_1 w \right) dy. $$

Then *reciprocity identities* are obtained by "testing" a problem with another cell solution, and exploiting the symmetry of the left-hand-side of (61).

**Relations between leading- and first-order coefficients.** Setting $w = P_2$ in (61) for $j = 1$, and $w = P_1$ for $j = 2$, one obtains $F_1(P_2) = F_2(P_1)$, which leads to $E_1 = E_0 \rho_1 / \rho_0$ as in [15, Lemma 1].

**Alternative expressions for the second-order coefficient $\beta$.** Similarly, one can consider the following combination:

$$ F_2(P_2) - F_1(P_3) + F_3(P_1) = E_2 - E_0 \frac{\rho_2}{\rho_0} + \int_0^1 \left( E - E_0 \frac{\rho}{\rho_0} \right) P_1^2 dy. $$
Starting from (60), and since $F_2(P_2) = \int_0^1 E(\partial_y P_2)^2 dy$ and $F_1(P_3) = F_3(P_1)$, one obtains:

$$ \beta = \frac{E_2}{E_0} - \frac{\rho_2}{\rho_0} = \int_0^1 \frac{E}{E_0} (\partial_y P_2)^2 + \left( \frac{\rho}{\rho_0} - \frac{E}{E_0} \right) P_1^2 dy = \int_0^1 \frac{E_0}{E} Q_1^2 - 2Q_1 P_1 + \frac{\rho}{\rho_0} P_1^2 dy. \quad (62) $$

Finally the Proposition 1 is proven by working a bit more on the expression (62). In particular, using the relations (59),

$$ \begin{aligned} \beta &= \int_0^1 (1 + \partial_y P_1) Q_1^2 - 2Q_1 P_1 + (1 + \partial_y Q_1) P_1^2 dy \\ &= \int_0^1 (P_1 - Q_1)^2 dy + \mathcal{I}, \end{aligned} \quad (63) $$

where the integral $\mathcal{I}$ can be expressed on the one hand using integration by part as:

$$ \mathcal{I} = \int_0^1 \partial_y P_1 Q_1^2 + \partial_y Q_1 P_1^2 dy = -2 \int_0^1 (\partial_y P_1 + \partial_y Q_1) P_1 Q_1 dy, \quad (64) $$

and, on the other hand, introducing the term $(P_1 - Q_1)^2$ and using (64):

$$ \begin{aligned} \mathcal{I} &= \int_0^1 (\partial_y P_1 + \partial_y Q_1) [(P_1 - Q_1)^2 + 2P_1 Q_1] + \partial_y P_1 P_1^2 + \partial_y Q_1 Q_1^2 dy, \\ &= \int_0^1 (\partial_y P_1 + \partial_y Q_1)(P_1 - Q_1)^2 dy - \mathcal{I}, \end{aligned} \quad (65) $$

because $\langle \partial_y P_1 P_1^2 \rangle_Y = \langle \partial_y P_1^3 \rangle_Y / 3 = 0$, and similarly $\langle \partial_y Q_1 Q_1^2 \rangle_Y = 0$. Finally, from (63), (65), and (59), a third expression of $\beta$, clearly positive, is obtained:

$$ \begin{aligned} \beta &= \frac{1}{2} \int_0^1 (P_1 - Q_1)^2 (2 + \partial_y P_1 + \partial_y Q_1) dy, \\ &= \frac{1}{2} \int_0^1 (P_1 - Q_1)^2 \left( \frac{E_0}{E} + \frac{\rho}{\rho_0} \right) dy \ge 0, \end{aligned} $$

### A.3. Closed-form formula for bilaminates

This section provides closed-form expressions, in our notation, for the effective coefficients, the parameters of an optimal (mt) model, and the cell functions $(P_1, Q_1, P_2, Q_2)$ assuming periodic bilaminate with the unit cell $Y = ]0, 1[ $ as depicted in figure 1. Additional detail and justifications are provided in [15] and the references therein, and the multi-laminate generalization is performed in [36].

**Effective coefficients.** The homogenization indexes at leading order are first defined as:

$$ n_E = \frac{\gamma_E}{(1 - \alpha) + \alpha \gamma_E} \quad \text{and} \quad n_\rho = \alpha + (1 - \alpha) \gamma_\rho, \quad (66) $$

in terms of the contrasts $\gamma_E = E_B/E_A$ and $\gamma_\rho = \rho_B/\rho_A$, so that the effective Young's modulus and density are:

$$ E_0 = n_E E_A, \quad \text{and} \quad \rho_0 = n_\rho \rho_A. $$

The effective wavespeed $c_0$ is such that $c_0^2 = c_A^2/n_0$ with $n_0 = n_\rho/n_E$ and $c_A = \sqrt{E_A/\rho_A}$. Finally, the second-order coefficient $\beta$ is:

$$ \beta = \frac{1}{12} \left[ \frac{\alpha(1 - \alpha)(1 - \gamma_E \gamma_\rho)}{n_0 \gamma_E} \right]^2. \quad (67) $$
**Optimized (mt) model.** The optimal two-parameter (mt) model from [15, Sect. 4.1.1] is given by:

$$ \beta_m = \frac{-1 - 4\beta + 4\bar{\beta}}{10} \quad \text{and} \quad \beta_t = \frac{1 - 6\beta - 4\bar{\beta}}{10} \qquad (68) $$

where

$$ \bar{\beta} = \frac{1}{12} \left[ \frac{\alpha^2 \gamma_E - (1 - \alpha)^2 \gamma_\rho}{n_0 \gamma_E} \right]^2 . $$

**First cell problem.** As already seen, $Q_0 = 1$ and the first cell function $P_1$ is given by:

$$ P_1(y) = \begin{cases} (n_E - 1)\left(y - \frac{\alpha}{2}\right), & y \in [0, \alpha], \\ -(n_E - 1)\frac{\alpha}{1 - \alpha}\left(y - \frac{1 + \alpha}{2}\right), & y \in [\alpha, 1], \end{cases} $$

in terms of the index $n_E$ given by (66).

**Second cell problem.** The cell function $Q_1$ is:

$$ Q_1(y) = \begin{cases} \left(\frac{1}{n_\rho} - 1\right)\left(y - \frac{\alpha}{2}\right), & y \in [0, \alpha], \\ -\left(\frac{1}{n_\rho} - 1\right)\frac{\alpha}{1 - \alpha}\left(y - \frac{1 + \alpha}{2}\right), & y \in [\alpha, 1], \end{cases} $$

in terms of the index $n_\rho$ given by (66). The affiliated “displacement” cell function $P_2$ is:

$$ P_2(y) = \begin{cases} P_2(0) + p_A y(y - \alpha), & y \in [0, \alpha], \\ P_2(0) + p_B [y(y - (1 + \alpha)) + \alpha], & y \in [\alpha, 1], \end{cases} $$

where

$$ P_2(0) = \frac{\alpha^2}{12} \left[ \frac{n_E}{n_\rho} + \frac{1 - 3\alpha}{\alpha} n_E - \frac{1 - \alpha}{\alpha n_\rho} + 1 \right] $$

and

$$ p_A = \frac{1}{2} \left[ \frac{n_E}{n_\rho} - 2n_E + 1 \right], \quad p_B = \frac{1}{2} \left( \frac{\alpha}{1 - \alpha} \right)^2 \left[ \frac{n_E}{n_\rho} + \frac{1 - 2\alpha}{\alpha} n_E - \frac{1}{\alpha n_\rho} + 1 \right]. $$

**Third cell problem.** Finally, the “stress” cell function $Q_2$ is:

$$ Q_2(y) = \begin{cases} Q_2(0) + q_A y(y - \alpha), & y \in [0, \alpha], \\ Q_2(0) + q_B [y(y - (1 + \alpha)) + \alpha], & y \in [\alpha, 1], \end{cases} $$

where

$$ Q_2(0) = \frac{n_E}{6} \left[ \alpha^3 q_A + \frac{(1 - \alpha)^3}{\gamma_E} q_B \right]. $$

and:

$$ q_A = \frac{1}{2} \left[ \frac{n_E}{n_\rho} - \frac{2}{n_\rho} + 1 \right], \quad q_B = \frac{1}{2} \left( \frac{\alpha}{1 - \alpha} \right)^2 \left[ \frac{n_E}{n_\rho} + \frac{1 - 2\alpha}{\alpha n_\rho} - \frac{n_E}{\alpha} + 1 \right]. $$

*Acknowledgements.* The authors wish to thank Sonia Fliss, Marie Touboul and Cédric Bellis for fruitful discussions. The numerous and relevant comments and corrections of the reviewers are also acknowledged.
REFERENCES

[1] A. ABDULLE AND T. POUCHON, *A priori error analysis of the finite element heterogeneous multiscale method for the wave equation over long time*, SIAM Journal on Numerical Analysis 54-3 (2016), 1507-1534.

[2] G. ALLAIRE, M. BRIANE AND M. VANNINATHAN, *A comparison between two-scale asymptotic expansions and Bloch wave expansions for the homogenization of periodic structures*, SeMA Journal 73-3 (2016), 237-259.

[3] G. ALLAIRE AND T. YAMADA, *Optimization of dispersive coefficients in the homogenization of the wave equation in periodic structures*, Numerische Mathematik 140-2 (2018), 265-326.

[4] G. ALLAIRE, A. LAMACZ-KEYMLING AND J. RAUCH, *Crime pays: homogenized wave equations for long times*, Asymptotic Analysis (2021), 1-42.

[5] I. ANDRIANOV, V. BOLSHAKOV, V. DANISHEVS AND D. WEICHERT, *Higher order asymptotic homogenization and wave propagation in periodic composite materials*, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 464 (2008), 1181-1201.

[6] S. ARMSTRONG, AND T. KUUSI, J.C. MOURRAT AND C. PRANGE, *Quantitative analysis of boundary layers in periodic homogenization*, Archive for Rational Mechanics and Analysis 226-2 (2017), 695-741.

[7] C. BELLIS AND B. LOMBARD, *Simulating transient wave phenomena in acoustic metamaterials using auxiliary fields*, Wave Motion 86 (2019), 175-194.

[8] C. BENETEAU, *Modèles homogénéisés enrichis en présence de bords : Analyse et traitement numérique*, PhD Thesis, Institut Polytechnique de Paris (2021),

[9] A. BENSOUSSAN, J.L. LIONS AND G. PAPANICOLAU, *Asymptotic Analysis for Periodic Structures*, North-Holland (1978).

[10] F. CAKONI, B. GUZINA AND S. MOSKOW, *On the homogenization of a scalar scattering problem for highly oscillating anisotropic media*, SIAM Journal on Mathematical Analysis 48-4 (2016), 2532-2560.

[11] F. CAKONI, B. GUZINA, S. MOSKOW AND T. PANGBURN, *Scattering by a bounded highly oscillating periodic medium and the effect of boundary correctors*, SIAM Journal on Applied Mathematics 79-4 (2019), 1448-1474.

[12] Y. CAPDEVILLE, L. GUILLOT AND J.J. MARIGO, *1-D non-periodic homogenization for the seismic wave equation*, Geophysical Journal International 181-2 (2010), 897-910.

[13] Y. CAPDEVILLE, *Homogenization of seismic point and extended sources*, Geophysical Journal International 226-2 (2021), 1390-1416.

[14] R. CORNAGGIA AND C. BELLIS, *Tuning effective dynamical properties of periodic media by FFT-accelerated topological optimization*, International Journal for Numerical Methods in Engineering 121-14 (2020), 3178-3205.

[15] R. CORNAGGIA AND B. GUZINA, *Second-order homogenization of boundary and transmission conditions for one-dimensional waves in periodic media*, International Journal of Solids and Structures, 188-189 (2020), 88-102.

[16] R.V. CRASTER, J. KAPLUNOV AND A.V. PICHugin, *High-frequency homogenization for periodic media*, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 466 (2010), 2341-2362.

[17] B. DELOURME, H. HADDAR AND P. JOLY, *Approximate models for wave propagation across thin periodic interfaces*, Journal de Mathématiques Pures et Appliquées 98-1 (2012), 28-71.

[18] J. FISH AND W. CHEN, *Higher-order homogenization of initial/boundary-value problem*, Journal of Engineering Mechanics 127-12 (2001), 1223-1230.

[19] J. FISH, W. CHEN AND G. NAGAI, *Non-local dispersive model for wave propagation in heterogeneous media: one-dimensional case*, International Journal for Numerical Methods in Engineering 54-3 (2002), 331-346.

[20] S. FLISS, *Wave propagation in periodic media : mathematical analysis and numerical simulation*, Habilitation à diriger des recherches, Université Paris Sud (Paris 11) (2019).

[21] S. FOREST AND K. SAB, *Finite-deformation second-order micromorphic theory and its relations to strain and stress gradient models*, Mathematics and Mechanics of Solids (2017), 1-21.

[22] D. GÉRARD-VARET AND N. MASMOUDI, *Homogenization and boundary layers*, Acta Mathematica 209-1 (2012), 133-178.

[23] E. GODLEWSKI AND P.A. RAVIART, *Numerical Approximation of Hyperbolic Systems of Conservation Laws*, Springer (1996).

[24] B. GUZINA, S. MENG AND O. 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 (2019), 20180547.

[25] D. HARUTYUNYAN, G.W. MILTON AND R.V. CRASTER, *High-frequency homogenization for travelling waves in periodic media*, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 472 (2016).

[26] M. JOSIEN, *Some quantitative homogenization results in a simple case of interface Communications in Partial Differential Equations*, Informa UK Limited (2019), 44, 907-939

[27] A. LAMACZ, *Dispersive effective models for waves in heterogeneous media*, Mathematical Models and Methods in Applied Sciences 21-9 (2011), 1871-1899.

[28] V. LAUDE, *Phononic crystals*, De Gruyter (2015).

[29] R.J. LEVEQUE, *Finite Volume Methods for Hyperbolic Problems*, Cambridge University Press (2002).

[30] B. LOMBARD AND J. PIRAUX, *How to incorporate the spring-mass conditions in finite-difference schemes*, SIAM Journal on Scientific Computing 24-4 (2003), 1379-1407.
[31] J.J. MARIGO AND A. MAUREL *Second Order Homogenization of Subwavelength Stratified Media Including Finite Size Effect SIAM Journal on Applied Mathematics* (2017) 77, 721-743.

[32] A. MAUREL AND J.J. MARIGO, *Sensitivity of a dielectric layered structure on a scale below the periodicity: A fully local homogenized model*, Physical Review B 98-2 (2018).

[33] A. MAUREL, K. PHAM AND J.J. MARIGO, *Scattering of gravity waves by a periodically structured ridge of finite extent*, Journal of Fluid Mechanics, Cambridge University Press (2019), 871, 350-376.

[34] S. MOSKOW AND M. VOGELIUS, *First-order corrections to the homogenised eigenvalues of a periodic composite medium. A convergence proof*, Proceedings of the Royal Society of Edinburgh: Section A Mathematics (1997) 127, 1263-1299.

[35] F. SANTOSA AND W.W. SYMES, *A dispersive effective medium for wave propagation in periodic composites*, SIAM Journal on Applied Mathematics 51-4 (1991), 984-1005.

[36] L. SCHWAN, N. FAVRIE, N. COTTEREAU AND B. LOMBARD, *Extended stress gradient elastodynamics: Wave dispersion and micro-macro identification of parameters*, International Journal of Solids and Structures 219-220 (2021), 34-50.

[37] D.P. SHAHRAKI AND B. GUZINA, *Homogenization of the wave equation with non-uniformly oscillating coefficients*, Mathematics and Mechanics of Solids (2022).

[38] D.P. SHAHRAKI AND B. GUZINA, *From d'Alembert to Bloch and back: A semi-analytical solution of 1D boundary value problems governed by the wave equation in periodic media*, International Journal of Solids and Structures 234-235 (2022), 111239.

[39] M. TOUBOUL, B. LOMBARD, C. BELLIS, *Time-domain simulation of wave propagation across resonant meta-interfaces*, Journal of Computational Physics 414-1 (2020), 109474.

[40] V. VINOLES, *Interface problems with metamaterials: modelling, analysis and simulations*, PhD Thesis, ENSTA Paris-Saclay (2016).

[41] A. WAUTIER AND B. GUZINA, *On the second-order homogenization of wave motion in periodic media and the sound of a chessboard*, Journal of the Mechanics and Physics of Solids 78 (2015), 382-414.