Marie Touboul
email: m.touboul@imperial.ac.uk

Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

Raphaël Assier
email: raphael.assier@manchester.ac.uk

Sébastien Guenneau
email: s.guenneau@imperial.ac.uk

Richard V Craster
email: r.craster@imperial.ac.uk

High-order homogenisation of the time-modulated wave equation: non-reciprocity for a single varying parameter

Laminated media with material properties modulated in space and time in the form of travelling waves have long been known to exhibit non-reciprocity. However, when using the method of low frequency homogenisation, it was so far only possible to obtain non-reciprocal effective media when both material properties are modulated in time, in the form of a Willis-coupling (or bi-anisotropy in electromagnetism) model. If only one of the two properties is modulated in time, while the other is kept constant, it was thought impossible for the method of homogenisation to recover the expected non-reciprocity since this Willis-coupling coefficient then vanishes. Contrary to this belief, we show that effective media with a single time-modulated parameter are non-reciprocal, provided homogenization is pushed to the second order. This is illustrated by numerical experiments (dispersion diagrams and time-domain simulations) for a bilayered modulated medium.

Introduction

Modulation of the electric permittivity in space and time was proposed 60 years ago in seminal papers by Simon, Oliner, Hessel and Cassedy [START_REF] Oliner | Wave propagation in a medium with a progressive sinusoidal disturbance[END_REF][START_REF] Simon | Action of a progressive disturbance on a guided electromagnetic wave[END_REF][START_REF] Cassedy | Dispersion relations in time-space periodic media: Part i-stable interactions[END_REF][START_REF] Cassedy | Dispersion relations in time-space periodic media part II-unstable interactions[END_REF]. It is observed that when modulating the properties in space and time, band diagrams are rotated in reciprocal space, and this is a general feature observed in all kinds of physical contexts, from optics [START_REF] Galiffi | Photonics of time-varying media[END_REF], to acoustics and elastic waves in solids and plates [START_REF] Nassar | Nonreciprocity in acoustic and elastic materials[END_REF], including surface Rayleigh waves [START_REF] Palermo | Surface wave non-reciprocity via time-modulated metamaterials[END_REF] and even flexural-gravity waves [START_REF] Farhat | Spacetime modulation in floating thin elastic plates[END_REF].

The physical interpretation of rotated, thus non-reciprocal, band diagrams is that directionality of space-time modulations such as travelling-wave modulations breaks time-reversal symmetry. The breaking of reciprocity induced by a time-harmonic modulation of the refractive index has been exploited in the realisation of photonic isolators and circulators [START_REF] Fang | Photonic aharonov-bohm effect based on dynamic modulation[END_REF] without the need of an applied external static magnetic field [START_REF] Śmigaj | Magneto-optical circulator designed for operation in a uniform external magnetic field[END_REF]. Similar effects have been explored in the context of acoustic and elastic waves [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF] or for acoustic modulated gratings [START_REF] Pham | Diffraction grating with space-time modulation[END_REF].

In-phase modulations of two parameters with the same strength within the acoustic, elastic and electromagnetic wave equations (e.g. the shear modulus and the compressional modulus, or the permittivity and the permeability) result in the closing of the high-frequency band gaps, at the same time as keeping the nonreciprocal character of these systems (induced by the rotated k-space). Interestingly, it has been shown in the electromagnetic context that giant linear non-reciprocity can be achieved at high-frequency within space-time modulated media [START_REF] Taravati | Giant linear nonreciprocity, zero reflection, and zero band gap in equilibrated space-time-varying media[END_REF].

A successful method to approximate wave propagation in heterogeneous media, for periodic media, is the two-scale asymptotic expansion method and the notion of slow and fast variables [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Bakhvalov | Homogenisation: Averaging Processes in Periodic Media[END_REF][START_REF] Cioranescu | An Introduction to Homogenization. An Introduction to Homogenization[END_REF][START_REF] Conca | Homogenization of periodic structures via bloch decomposition[END_REF]. Low-frequency homogenisation of laminates with in-phase modulations of two parameters leads to Willis-coupling [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF] in acoustics and magneto-electric coupling [START_REF] Huidobro | Fresnel drag in space-time-modulated metamaterials[END_REF][START_REF] Huidobro | Homogenization theory of spacetime metamaterials[END_REF] in electromagnetics. A leadingorder asymptotic analysis of the low-frequency regime for the time-modulated wave equation is also addressed in [START_REF] Lurie | An Introduction to the Mathematical Theory of Dynamic Materials[END_REF]. Homogenisation of a larger class of periodic multiscale space-time media has been studied with various techniques in the context of linear and nonlinear parabolic and hyperbolic equations [START_REF] Holmbom | Multiscale convergence and reiterated homogenization of parabolic problems[END_REF][START_REF] Flodén | Homogenization of some parabolic operators with several time scales[END_REF][START_REF] Ariel | Iterated averaging of three-scale oscillatory systems[END_REF][START_REF] Dehamnia | Multiscale homogenization of nonlinear hyperbolic-parabolic equations[END_REF]. This makes possible interesting features ranging from frequency-momentum transitions to compression and amplification of electromagnetic signals, one-way hyperbolic metasurfaces, and even non-Hermitian and topological phenomena. However, while non-reciprocity is expected with only one parameter modulated in time, the leading-order low-frequency homogenisation does not encapsulate this behaviour [START_REF] Lurie | An Introduction to the Mathematical Theory of Dynamic Materials[END_REF]. Therefore, leading-order homogenisation is able to describe nonreciprocity when two parameters are modulated in time, but not as soon as one of them is kept constant.

However, there is no clear physical reason for this necessity of in-phase modulations of two parameters to get non-reciprocity at low frequency. Moreover, while in elasticity the modulation of the elastic constants is possible [START_REF] Danas | Experiments and modeling of iron-particlefilled magnetorheological elastomers[END_REF][START_REF] Swinteck | Bulk elastic waves with unidirectional backscattering-immune topological states in a time-dependent superlattice[END_REF][START_REF] Deymier | Tailoring phonon band structures with broken symmetry by shaping spatiotemporal modulations of stiffness in a one-dimensional elastic waveguide[END_REF][START_REF] Bou Matar | Band gap tunability of magneto-elastic phononic crystal[END_REF][START_REF] Vasseur | Band structures tunability of bulk 2d phononic crystals made of magneto-elastic materials[END_REF], it seems much more complicated to modulate the mass density. In a similar fashion, modulating the permittivity in electromagnetism has been achieved [START_REF] Tirole | Saturable time-varying mirror based on an epsilon-near-zero material[END_REF][START_REF] Amr | Spatiotemporal light control with active metasurfaces[END_REF][START_REF] Caloz | Spacetime metamaterials-part i: General concepts[END_REF] but modulating the permeability while natural materials are non magnetic is much more demanding. This is also the case for acoustic waves propagating in a fluid, where a modulation of both physical parameters cannot be done without generating a background flow or breaking mass conservation [START_REF] Martin | Acoustics and dynamic materials[END_REF]. This motivates the present work that establishes the second-order homogenised model (see [START_REF] Boutin | Rayleigh scattering in elastic composite materials[END_REF][START_REF] Chen | A dispersive model for wave propagation in periodic heterogeneous media based on homogenization with multiple spatial and temporal scales[END_REF][START_REF] Fish | Higher-order homogenization of initial/boundary-value problem[END_REF][START_REF] Fish | Space-time multiscale model for wave propagation in heterogeneous media[END_REF][START_REF] Lamacz | Dispersive effective models for waves in heterogeneous media[END_REF][START_REF] Allaire | Optimization of dispersive coefficients in the homogenization of the wave equation in periodic structures[END_REF] in the non modulated case) of the time-modulated scalar wave equation and proves that non-reciprocity is achieved even if only one of the parameters is modulated.

In Section 2, the configuration of a time-modulated laminate is recalled in the framework of antiplane elasticity and homogenisation is performed up to the second order. For each order, we study the limit case of only one parameter being modulated and it is proved that non-reciprocity still happens in this case provided we consider the second-order model. In Section 3, dispersion diagrams obtained by Bloch-Floquet theory and by homogenisation are compared, to validate the second-order model and to quantify the non-reciprocal effects. Time-domain simulations are also performed in the microstructured medium to validate the conclusions obtained through homogenisation. Finally, Appendix B explains why these results presented in the framework of antiplane elasticity, i.e. of non-reciprocity even if only one of the parameters is modulated in time, still hold when considering transverse electromagnetic waves.

Homogenisation of a time-modulated laminate in antiplane elasticity

We first consider a 1D h-periodic elastic medium of shear modulus µ h and mass density ρ h in the antiplane elasticity framework. Therefore the equations for the elastic displacement U h (X, T ) and the stress Σ h (X, T ) read

∂ T (ρ h ∂ T U h ) = ∂ X Σ h + F with Σ h = µ h ∂ X U h , (1)
where X (in m) and T (in s) denote the space coordinate and the time, respectively, and F (X, T ) is a source term. The partial derivatives with respect to space and time are denoted by ∂ X and ∂ T , respectively. Moreover, the physical parameters ρ h (in kg m -3 ) and µ h (in Pa) are modulated in space and time in a wave-like fashion of modulation speed c m :

ρ h (X, T ) = ρ(X -c m T ) and µ h (X, T ) = µ(X -c m T ) (2)
with ρ and µ h-periodic functions assumed to be piecewise continuous, see Figure 1.

In order to satisfy some smoothness assumptions [START_REF] Cassedy | Dispersion relations in time-space periodic media: Part i-stable interactions[END_REF][START_REF] Cassedy | Dispersion relations in time-space periodic media part II-unstable interactions[END_REF], we further assume that c m satisfies:

c m < min ξ∈(0,h) c(ξ) or c m > max ξ∈(0,h) c(ξ) with c = µ ρ . ( 3
)
Figure 1: Bilaminate modulated in a wave-like fashion at T = 0 (top) and T = T 0 (bottom).

Moving frame

As in [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF], we reformulate the problem in the moving frame, which will be useful to determine the continuity conditions for the fields, and for time-domain numerical simulations in Section 3.

The momentum equation and the constitutive law (1) can be written

∂ X V h = ∂ T Σ h µ h with ∂ X Σ h = ∂ T (ρ h V h ) -F, (4)
where

V h (X, T ) = ∂ T U h (X, T ) is the velocity. Setting W h = V h Σ h , N h = 0 -1/µ h -ρ h 0 , F = 0 -F , ( 5
)
the equations (4) are written

∂ X W h (X, T ) + ∂ T (N h (X, T ) W h (X, T )) = F (X, T ). (6)
From [START_REF] Simon | Action of a progressive disturbance on a guided electromagnetic wave[END_REF], it follows that

N h (X, T ) = N (X -c m T ) = 0 -1/µ(X -c m T ) -ρ(X -c m T ) 0 . ( 7
)
One introduces a moving frame with coordinates (ξ, T ) = (X -c m T, T ); the fields in this frame are denoted by g(ξ, T ) = g(X, T ) and satisfy the chain rules

∂ ξ g = ∂ X g and ∂ T g = ∂ T g + c m ∂ X g. ( 8
)
One then defines

Ψ(ξ, T ) = (I -c m N (ξ)) Wh (ξ, T ) =   Ṽh (ξ, T ) + c m µ(ξ) Σh (ξ, T ) Σh (ξ, T ) + c m ρ(ξ) Ṽh (ξ, T )   ≡ Ψ 1 (ξ, T ) Ψ 2 (ξ, T ) , B(ξ) = N (ξ) (I -c m N (ξ)) -1 = 1 c(ξ) 2 -c 2 m c m -1/ρ(ξ) -µ(ξ) c m , ( 9
)
where I is the 2 × 2 identity matrix. In the limit-case c m = 0, Ψ = W h . Using the chain rules ( 8) on (6) yields

∂ ξ Ψ(ξ, T ) + B(ξ) ∂ T Ψ(ξ, T ) = F (ξ, T ). ( 10
)
Let us consider an interface at ξ n . Integrating [START_REF] Śmigaj | Magneto-optical circulator designed for operation in a uniform external magnetic field[END_REF] on [ξ n -ε, ξ n + ε] with ε → 0 gives Ψ ξn = 0, where g a = g(a + ) -g(a -) refers to the jump across a. Using (9), these jump conditions in the moving frame read Ṽh + c m µ

Σh ξn = 0, Σh + c m ρ Ṽh ξn = 0, (11)
which are the jump conditions one would obtain by Reynolds' transport theorem. 8) and (4), the first jump condition in [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF]

Remark 1 Since ∂ T Ũh (ξ, T ) = V h (X, T ) + cm µ(ξ) Σ h (X, T ) with (
leads to ∂ T Ũh ξn = ∂ T Ũh ξn = 0
in the moving frame where the interfaces are at fixed positions ξ n . Since the jump is zero at T = 0, then we get Ũh ξn =0. Equations [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF] are therefore equivalent to continuity for both Ũh and Σh + c m ρ Ṽh in the moving frame. Therefore, the same holds for both U h and Σ h + c m ρ h V h in the fixed frame.

Setting

A(ξ) = B -1 (ξ) = -c m -1/ρ(ξ) -µ(ξ) -c m , ( 12
)
and

Φ(ξ + c m T, T ) = B -1 (ξ) F (ξ, T ) = F (ξ + c m T, T ) c m F (ξ + c m T, T ) , ( 13
)
the system (10) leads to

∂ T Ψ(ξ, T ) + A(ξ) ∂ ξ Ψ(ξ, T ) = Φ(ξ + c m T, T ). ( 14
)
The source term is therefore no longer fixed in the moving frame. The eigenvalues of A are ±c 0 -c m . Lastly, an energy balance of ( 14) without source term leads to

dE dT = 0, with E(T ) = 1 2 ρ(ξ)Ψ 1 (ξ, T ) 2 + 1 µ(ξ) Ψ 2 (ξ, T ) 2 dξ. ( 15
)
This can be interpreted as the conservation of the sum of a kinetic and a potential energy in the moving frame. However, in the fixed frame, there is no conservation of energy due to the time modulation of the constitutive properties [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF].

Non-dimensionalisation

In order to simplify the mathematical notations, we start by non-dimensionalising the physical problem (1). To do so, we introduce a characteristic dimensional density ρ , Young's modulus µ , wavespeed c = µ /ρ , and wavenumber k . These can be used in order to define the following non-dimensional quantities

η = k h, x = k X, t = (k c )T, u η = k U h , σ η = 1 µ Σ h , α = ρ ρ , β = µ µ , c = c m c , f = F k µ .
(16) Using these quantities, (1) can be rewritten as the non-dimensional governing equation

∂ t (α∂ t u η ) = ∂ x (β∂ x u η ) + f ( 17
)
together with continuity for u η and σ η + c ∂ t u η at interfaces (see Remark 1).

Remark 2 Due to (3), β -c 2 α is either strictly positive or strictly negative everywhere.

Two-scale analysis

We now make the assumption that the characteristic wavelength λ = 2π k is much larger than the periodicity h so that η 1 and is the small positive parameter of the low-frequency setting. The material parameters α and β vary on a fine scale associated with the rescaled coordinate y = k ξ =

x-ct η . Following the two-scale asymptotic technique, we further assume that the field u η have small-scale variations that are described by y, and slow continuous features as well, which can be described by the variable x. Accordingly the field u η is expanded using the following ansatz

u η (x, t) = j≥0 η j u j (x, y, t) (18)
where x and y are assumed to be independent variables implying that

∂ x ↔ ∂ x + 1 η ∂ y and ∂ t ↔ ∂ t -c η ∂ y .
We introduce the stress and momentum

σ j = β(∂ x u j + ∂ y u j+1 ) and m j = α(∂ t u j -c∂ y u j+1 ). ( 19
)
To ensure the correct continuity and periodicity conditions for u η and σ η , the quantities u j , σ j and m j are assumed to have the following properties: u j continuous with respect to the first variable, u j 1-periodic with respect to the second variable, u j and σ j + cm j continuous with respect to the second variable (see Remark 1).

The non-dimensional problem can hence be rewritten as the governing equation [START_REF] Lurie | An Introduction to the Mathematical Theory of Dynamic Materials[END_REF] together with the properties given above.

∂ t - c η ∂ y   α(y) ∂ t - c η ∂ y j≥0 η j u j   = ∂ x + 1 η ∂ y   β(y) ∂ x + 1 η ∂ y j≥0 η j u j   + f (x, t),

Homogenised equation at leading order

The homogenised model at leading order has been obtained in various references, see [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF] for elasticity. We summarize the procedure here for sake of completeness. We start by introducing the mean value g = 1 0 g(y)dy for a function g(y). Identifying the terms of order η -2 in (20) we get

∂ y ((β -c 2 α)∂ y u 0 ) = 0, (21)
which leads to

∂ y u 0 (x, y, t) = 1 β(y) -c 2 α(y) q(x, t),
for some unknown function q(x, t). Integrating this equation on a unit cell leads to q(x, t) = 0, implying that u 0 (x, y, t) = U 0 (x, t). [START_REF] Flodén | Homogenization of some parabolic operators with several time scales[END_REF] Collecting now the terms of order η -1 in (20), we get

∂ y (β -c 2 α)∂ y u 1 = -c∂ y (α∂ t u 0 ) -∂ y (β∂ x u 0 ). ( 23
)
We can then write u 1 as

u 1 (x, y, t) = U 1 (x, t) + P (y)∂ x U 0 (x, t) + Q(y)∂ t U 0 (x, t), ( 24
)
with P and Q satisfying the following cell problems

   ∂ y (β -c 2 α)P + β = 0
1-periodicity of P , P = 0, and continuity of P and (β -c 2 α)P + β, [START_REF] Danas | Experiments and modeling of iron-particlefilled magnetorheological elastomers[END_REF] and

   ∂ y (β -c 2 α)Q + cα = 0 1-periodicity of Q, Q = 0, and continuity of Q and (β -c 2 α)Q + cα, ( 26
)
where hereafter g will denote the derivative of any function g of one variable. Then, collecting terms of order η 0 we get

∂ t m 0 -c∂ y m 1 = ∂ x σ 0 + ∂ y σ 1 + f, (27)
which, upon integration on a unit cell and due to continuity of σ 1 + cm 1 , leads to

∂ t m 0 = ∂ x σ 0 + f. ( 28
)
Using [START_REF] Dehamnia | Multiscale homogenization of nonlinear hyperbolic-parabolic equations[END_REF], we get

m 0 = α(1 -c Q ) ∂ t U 0 -c α P ∂ x U 0 and σ 0 = β Q ∂ t U 0 + β(1 + P ) ∂ x U 0 . ( 29
)
We can prove by integrations by parts (see reciprocity identities in Appendix A, and equation (120) in the present case) that the so-called Willis coupling coefficient satisfies

W 0 := c α P = β Q . ( 30
)
We also introduce the effective parameters α 0 = α(1 -cQ (y)) and β 0 = β(1 + P (y)) [START_REF] Amr | Spatiotemporal light control with active metasurfaces[END_REF] so that the leading-order equation [START_REF] Bou Matar | Band gap tunability of magneto-elastic phononic crystal[END_REF] can finally be written

α 0 ∂ 2 tt U 0 -2W 0 ∂ 2 tx U 0 -β 0 ∂ 2 xx U 0 = f. ( 32
)
Upon coming back to dimensionalised variables using [START_REF] Cioranescu | An Introduction to Homogenization. An Introduction to Homogenization[END_REF], setting f = 0, and looking for wave-like solutions, we obtain the dispersion relation

-α 0 ω 2 + 2W 0 c ωk + β 0 (c ) 2 k 2 = 0 ( 33
)
where ω and k are the dimensionalised angular frequency and wavenumber, respectively. Odd orders of space derivatives in the effective equation [START_REF] Caloz | Spacetime metamaterials-part i: General concepts[END_REF], or in an equivalent manner odd powers of k in the dispersion relation [START_REF] Martin | Acoustics and dynamic materials[END_REF], implies that reciprocity is broken since there is no longer symmetry between k and -k. However, this occurs only if the Willis coupling coefficient W 0 is non-zero.

Remark 3 What happens at leading order if only one parameter is modulated in time?

If α does not depend on y, we can choose ρ such that α = 1. Then W 0 = c αP = c P = 0, the Willis-coupling term vanishes and the effective equation is the usual one

α 0 ∂ 2 tt U 0 -β 0 ∂ 2 xx U 0 = f. ( 34
)
In preparation for the following orders, we also notice that Q = 0 and α 0 = 1 in this case. Similarly, if β = 1, then W 0 = βQ = Q = 0 and the Willis-coupling term vanishes. We also notice that P = 0 and β 0 = 1.

Consequently, the Willis coupling coefficient W 0 vanishes if either α or β does not depend on y, therefore non-reciprocity appears in this leading-order effective equation only if both physical parameters are modulated in space and time. This is consistent with results of the literature of low-frequency homogenization of time-modulated media [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF][START_REF] Huidobro | Homogenization theory of spacetime metamaterials[END_REF]. However, we expect non-reciprocity even if one parameter is modulated in time according to Bloch-Floquet analysis and dispersion relations obtained in the literature [START_REF] Cassedy | Dispersion relations in time-space periodic media: Part i-stable interactions[END_REF][START_REF] Cassedy | Dispersion relations in time-space periodic media part II-unstable interactions[END_REF]. In the forthcoming sections, we prove that the modulation of both parameters is not a necessary condition for non-reciprocity in the low-frequency setting, but rather a limitation of the leading-order homogenised model.

Homogenised equation at order 1

We start by inserting [START_REF] Huidobro | Homogenization theory of spacetime metamaterials[END_REF] for j = 0, 1 in ( 27):

∂ y (β -c 2 α)∂ y u 2 = α∂ 2 tt u 0 -β∂ 2 xx u 0 -αc∂ 2 ty u 1 -c∂ y (α∂ t u 1 ) -β∂ 2 xy u 1 -∂ y (β∂ x u 1 ) -f. ( 35
)
Using [START_REF] Dehamnia | Multiscale homogenization of nonlinear hyperbolic-parabolic equations[END_REF], and (32) to replace the ∂ 2 tt U 0 terms, we can write u 2 as

u 2 (x, y, t) = U 2 (x, t) + P (y)∂ x U 1 (x, t) + Q(y)∂ t U 1 (x, t) + R(y)∂ 2 xx U 0 (x, t) + S(y)∂ 2 xt U 0 (x, t) + B(y)f (x, t) (36)
where R, S and B are solutions of the following cell problems

       ∂ y (β -c 2 α)R + βP + β 0 α 0 cαQ = -β(1 + P ) + β 0 α 0 α(1 -cQ ) 1-periodicity of R, R = 0, and continuity of R and (β -c 2 α)R + βP + β 0 α 0 αcQ, (37)
       ∂ y (β -c 2 α)S + βQ + cαP + 2W 0 α 0 cαQ = -αcP -βQ + 2W 0 α 0 α(1 -cQ )
1-periodicity of S, S = 0, and continuity of S and (β -

c 2 α)S + βQ + cαP + 2W 0 α 0 cαQ, ( 38
)        ∂ y (β -c 2 α)B + cα α 0 Q = α(1 -cQ ) α 0 - 1
1-periodicity of B, B = 0, and continuity of B and (β -

c 2 α)B + cα α 0 Q, ( 39
)
with the continuity conditions coming from continuity of u 2 and σ 1 + cm 1 . Collecting the terms of order η in (20) yields:

∂ t m 1 -c∂ y m 2 = ∂ x σ 1 + ∂ y σ 2 (40)
which once averaged gives

∂ t α(∂ t u 1 -c∂ y u 2 ) = ∂ x β(∂ x u 1 + ∂ y u 2 ) . ( 41
)
Inserting ( 24) and ( 36) in ( 41) leads to the final effective equation for the first-order mean field:

α 0 ∂ 2 tt U 1 -2W 0 ∂ 2 tx U 1 -β 0 ∂ 2 xx U 1 = F(U 0 ) + A(f ) ( 42
)
where the first source term reads

F(g) = -αQ ∂ 3 ttt g + -αP + c αS ∂ 3 ttx g + c αR + βQ + βS ∂ 3 txx g + βP + βR ∂ 3 xxx g ( 43
)
for any function g(X, t) and the second one reads

A(f ) = βB ∂ x f + c αB ∂ t f. ( 44
)

Remark 4 What happens at order 1 if only one parameter is modulated in time?

Let us assume that mass density is constant (α = 1). Because the cell problem solutions are continuous and have zero mean-values, the source term simplifies to

F(U 0 ) = βS ∂ 3 txx U 0 + βP + βR ∂ 3 xxx U 0 . ( 45
)
The reciprocity identities (123), ( 121) and (125), given in Appendix A allow us to conclude that

F(U 0 ) = 0 and A(f ) = 0. ( 46
)
Similarly, if β = 1, the source terms simplify to

F(U 0 ) + A(f ) = -αQ ∂ 3 ttt U 0 + c αS ∂ 3 ttx U 0 + c αR ∂ 3 txx U 0 + c αB ∂ t f, ( 47
)
and the reciprocity identities (124), (122), and (126), together with [START_REF] Boutin | Rayleigh scattering in elastic composite materials[END_REF] and W 0 = 0 lead to

F(U 0 ) + A(f ) = 0. ( 48
)
Therefore, as soon as only one parameter is modulated in time, the effective equation reduces to

α 0 ∂ 2 tt U 1 -β 0 ∂ 2 xx U 1 = 0 (49)
and the dispersion relation is the same as at the leading order.

Consequently, if either α or β does not depend on y, both the Willis coupling coefficient W 0 and the source term F(U 0 ) vanish; then in the absence of source f the equation for the first-order mean field U 1 is the same as for the zeroth-order one. Therefore, up to the first order, non-reciprocity appears in the effective equation only if both physical parameters are modulated in space and time, but we will see in the next section that it is not the case at order 2.

Homogenised equation at order 2

We start by inserting [START_REF] Huidobro | Homogenization theory of spacetime metamaterials[END_REF] in [START_REF] Lombard | How to incorporate the spring-mass conditions in finite difference schemes[END_REF] to get an equation for u 3 . We then make use of the expressions obtained for u 1 (24) and u 2 [START_REF] Fish | Higher-order homogenization of initial/boundary-value problem[END_REF]. With a similar methodology as for the previous order, we replace the terms ∂ 2 tt U 1 , ∂ 3 ttt U 0 , and ∂ 3 xxx U 0 , using [START_REF] Vinoles | Problèmes d'interface en présence de métamatériaux : modélisation, analyse et simulations[END_REF], [START_REF] Caloz | Spacetime metamaterials-part i: General concepts[END_REF] differentiated with respect to time, and (32) differentiated with respect to space, respectively. The field u 3 can therefore be written

u 3 (x, y, t) = U 3 (x, t) + P (y)∂ x U 2 (x, t) + Q(y)∂ t U 2 (x, t) + R(y)∂ 2 xx U 1 (x, t) + S(y)∂ 2 tx U 1 (x, t) + M (y)∂ 3 ttx U 0 (x, t) + N (y)∂ 3 txx U 0 (x, t) + B x (y)∂ x f (x, t) + B t (y)∂ t f (x, t) (50)
with M and N solutions of the following cell problems

       ∂ y (β -c 2 α)M + cαS + α 0 β 0 βR = α P -cS + 2W 0 α 0 Q - α 0 β 0 (βP + βR )
1-periodicity of M , M = 0, and continuity of M and (β -

c 2 α)M + cαS + α 0 β 0 βR, ( 51
)          ∂ y (β -c 2 α)N + cαR + βS - 2W 0 β 0 βR = -β(Q + S ) + α β 0 α 0 Q -cR + 2W 0 β 0 (βP + βR )
1-periodicity of N , N = 0, and continuity of N and (β -

c 2 α)N + cαR + βS - 2W 0 β 0 βR, ( 52
)    ∂ y (β -c 2 α)B t + cαB = α α 0 Q -αcB 1-periodicity of B t , B t = 0, and continuity of B t and (β -c 2 α)B t + cαB, ( 53
)        ∂ y (β -c 2 α)B x - 1 β 0 βR + βB = β β 0 (P + R ) -βB 1-periodicity of B x , B x = 0, and continuity of B x and (β -c 2 α)B x - 1 β 0 βR + βB, (54)
where the continuity conditions come from the continuity of u 3 and σ 2 +cm 2 . Moreover, when replacing the high-order temporal derivatives, some mean values appear as source terms of the cell problems. They can be shown to cancel out using the reciprocity identitites (124) and ( 121) in ( 51), the identities (121), ( 122), ( 123) in ( 52), the identities (126) in (53), and finally ( 121) and ( 125) in (54).

Collecting the terms of order η 2 in [START_REF] Lurie | An Introduction to the Mathematical Theory of Dynamic Materials[END_REF] and averaging the equation so obtained yields:

∂ t α(∂ t u 2 -c∂ y u 3 ) = ∂ x β(∂ x u 2 + ∂ y u 3 ) . ( 55
)
Using ( 36) and ( 50) leads to the final effective equation for the mean field at the second order

α 0 ∂ 2 tt U 2 -2W 0 ∂ 2 tx U 2 -β 0 ∂ 2 xx U 2 = F(U 1 ) + E(U 0 ) + B(f ) (56)
with F defined in [START_REF] Maurel | Sensitivity of a dielectric layered structure on a scale below the periodicity: A fully local homogenized model[END_REF], the second source term given by

E(U 0 ) = -α(S -cM ) ∂ 4 tttx U 0 -α(R -cN ) ∂ 4 ttxx U 0 + βR ∂ 4 xxxx U 0 + β(S + N ) ∂ 4 txxx U 0 + βM ∂ 4 ttxx U 0 (57)
and the third one defined as

B(f ) = β(B + B x ) ∂ 2 xx f + βB t + cαB x ∂ 2 tx f + α(cB t -B) ∂ 2 tt f. ( 58
)
We already now that if either α or β is constant, W 0 vanishes in (56). We will show that the same holds for F(U 1 ). The question is therefore to know whether the terms with an odd number of spatial derivatives (associated with non-reciprocial behaviour) in E cancel out if only one of the parameters is modulated in time: This is the subject of the next two subsections.

If ρ is constant

Let us first consider the case where ρ and therefore α is constant. As a reminder, ρ is chosen such that α = 1. The reciprocity identities (123) and (121) given in Appendix A allow us to conclude that F(U 1 ) = 0. Similarly, in this case B = B t = 0 so that the source term in (58) reduces to

B(f ) = βB x ∂ 2 xx f (59)
and the effective equation for U 2 is

∂ 2 tt U 2 -β 0 ∂ 2 xx U 2 = E(U 0 ) + βB x ∂ 2 xx f. ( 60
)
Since S = R = 0, and M = N = 0 by continuity of M and N , the source term (57) then reduces to

E(U 0 ) = βR ∂ 4 xxxx U 0 + βM ∂ 4 ttxx U 0 + β(S + N ) ∂ 4 txxx U 0 . ( 61
)
Remark 5 If c = 0, the effective equation for the second-order term reduces to

∂ 2 tt U 2 -β 0 ∂ 2 xx U 2 -β 0 P 2 ∂ 4 xxxx U 0 = P 2 ∂ 2 xx f (62)
and we therefore recover the high-order homogenized equation of [START_REF] Allaire | Optimization of dispersive coefficients in the homogenization of the wave equation in periodic structures[END_REF].

Consequently, non-reciprocity will be achieved if the following coefficient is non zero:

N β = β(S + N ) . ( 63
)
Using the reciprocity identity (129) in the case α = 1 allows us to get

βN = c RP -c R P + βSP -βS P . ( 64
)
Similarly, (127) leads to βSP -βS P = -βS

+ c P R -P R . ( 65
)
Combining ( 64) and (65) together with the continuity of P and R allows us to simplify N β to

N β = -4c P R . ( 66
)
Since P and R appear in (66), let us first rewrite the cell problems ( 25) and ( 37) when α = 1:

   ∂ y (β -c 2 )P + β = 0
1-periodicity of P , P = 0, and continuity for P and (β -c 2 )P + β, (67)

   ∂ y (β -c 2 )R + βP = -β(1 + P ) + β 0
1-periodicity of R, R = 0, and continuity for R and (β -c 2 )R + βP .

(68)

Integration of (67) leads to (β

-c 2 )P + β = C ∂P , ( 69
)
where C ∂P is a constant of integration, which, after averaging this equation on a unit cell, can be shown to be C ∂P = β 0 . Consequently, we get from (69) that the first equation of (68) reduces to

∂ y (β -c 2 )R + βP = -c 2 P . ( 70
)
It follows that

R = - β + c 2 β -c 2 P + C ∂R β -c 2 , ( 71
)
for some integration constant C ∂R , which, using that R = 0 is given by

C ∂R = β + c 2 β -c 2 P 1 β -c 2 -1 . ( 72
)
Therefore N β can be written in terms of P only:

N β = -4c - β + c 2 β -c 2 P 2 + β + c 2 β -c 2 P 1 β -c 2 -1 1 β -c 2 P . ( 73
)
We know that β -c 2 is of constant sign depending on whether c m > max c h or c m < min c h , and we have now to treat both cases separately.

Case 1: β -c 2 > 0 In this subsonic regime, we introduce the following scalar product

(f, g) β,sub = 1 β -c 2 f g (74)
where sub denotes the subsonic regime, and we choose µ so that

(1, 1) β,sub = 1 β -c 2 = 1. ( 75
)
Then, N β simplifies to

N β = 4c P 2 + 2c 2 (P, P ) β,sub (1, 1) β,sub -(P, 1) 2 β,sub . ( 76
)
We know that P 2 > 0. The term in the bracket, combined with (75), is positive due to Cauchy-Schwarz inequality. Consequently N β > 0 and we get non-reciprocity.

Case 2: β -c 2 < 0 In this supersonic case, we introduce the following scalar product

(f, g) β,sup = β c 2 -β f g (77)
and we choose µ so that (1, 1) β,sup = β c 2 -β = 1. Using P = 0, the terms in N β can be written as:

           1 β -c 2 = 1 c 2 -1 + β β -c 2 = - 2 c 2 , - β + c 2 β -c 2 P 2 = P 2 + 2 β c 2 -β P 2 - β + c 2 β -c 2 P = 2 β c 2 -β P , and 1 β -c 2 P = - 1 c 2 β c 2 -β P (78)
so that N β simplifies to

N β = -4c P 2 + 2 (P, P ) β,sup (1, 1) β,sup - 1 2 (P, 1) 2 β,sup . ( 79
)
Consequently

N β ≤ -4c P 2 + 2 (P, P ) β,sup (1, 1) β,sup -(P, 1) 2 β,sup < 0 (80)
where the last inequality comes from the fact that P 2 > 0 and (P, P ) β,sup (1, 1) β,sup -(P, 1) 2 β,sup ≥ 0 with Cauchy-Schwarz inequality.

If µ is constant

Let us now consider the second case where µ, and consequently β, is constant. As a reminder, µ is chosen such that β = 1. In this case, the reciprocity identities (124), (122) together with (49) leads to F(U 1 ) = 0, while the source term (58) reduces to

B(f ) = c αB x ∂ 2 tx f -α(B -cW t ) ∂ 2 tt f (81)
so that the mean equation for U 2 reads

α 0 ∂ 2 tt U 2 -∂ 2 xx U 2 = E(U 0 ) + c αB x ∂ 2 tx f -α(B -cW t ) ∂ 2 tt f. ( 82
)
Eventually, the source term (57) reduces to

E(U 0 ) = -α(S -cM ) ∂ 4 tttx U 0 -α(R -cN ) ∂ 4 ttxx U 0 (83)
so that what matters is to prove that N α = -α(S -cM ) is non zero. Using (130) and (127) together with continuity of Q and R, we get N α = -4α 0 QR with Q (26) and R (37) solutions of the following simpler cell problems when β = 1

   ∂ y (1 -c 2 α)Q + cα = 0
1-periodicity of Q, Q = 0, and continuity for Q and (1 -

c 2 α)Q + cα, (84)
       ∂ y (1 -c 2 α)R + cα α 0 Q = -1 + α α 0 (1 -cQ )
1-periodicity of R, R = 0, and continuity for R and (1 -

c 2 α)R + αc α 0 Q. ( 85
)
Similarly to the case α = 1, we get (1 -

c 2 α)Q + cα = C ∂Q with C ∂Q = cα 0 and R = - 1 + αc 2 cα 0 (1 -αc 2 ) Q + C ∂R 1 -αc 2 with C ∂R = 1 cα 0 1 + αc 2 1 -αc 2 Q 1 1 -αc 2 -1 . ( 86
)
Consequently,

N α = 4 c 1 + αc 2 1 -αc 2 Q 2 - 1 1 -αc 2 Q 1 1 -αc 2 -1 1 + αc 2 1 -αc 2 Q . ( 87
)
We know that 1 -αc 2 has constant sign and we need to consider two cases again.

Case 1: αc 2 > 1 In this case, we introduce the scalar product

(f, g) α,sup = 1 αc 2 -1 f g (88)
and we choose ρ such that (1, 1) α,sup = 1. Consequently N α writes

N α = 4 c 2 -(Q, Q) α,sup (1, 1) α,sup + (Q, 1) 2 α,sup -Q 2 . ( 89
)
Cauchy-Schwarz inequality for the first bracket and Q 2 > 0 yields N α < 0.

Case 2: αc 2 < 1 In this last case, we introduce the scalar product

(f, g) α,sub = αc 2 1 -αc 2 f g (90)
and we choose ρ so that (1, 1) α,sub = 1. Consequently, we get the following simplifications, using Q = 0:

         1 1 -αc 2 = 2 , 1 + αc 2 1 -αc 2 Q 2 = Q 2 + 2 αc 2 1 -αc 2 Q 2 , 1 1 -αc 2 Q = αc 2 1 -αc 2 Q and 1 + αc 2 1 -αc 2 Q = 2 αc 2 1 -αc 2 Q (91) so that N α writes N α = 4 c Q 2 + 2 (Q, Q) α,sub (1, 1) α,sub - 1 2 (Q, 1) 2 α,sub . ( 92
)
Therefore,

N α ≥ 4 c Q 2 + 2 (Q, Q) α,sub (1, 1) α,sub -(Q, 1) 2 α,sub > 0 ( 93
)
where we used for the last inequality that Q 2 > 0 and Cauchy-Schwarz inequality. Consequently, in each case, we recover the fact that non-reciprocity is present through the source term that appears in the second-order homogenised model even if only one of the physical parameters is modulated in time. We will illustrate this property in the next section through the expression of the dispersion relation.

Total mean field and dispersion relation

The second-order mean field ū2 is defined by

ū2 (x, t) = U 0 (x, t) + ηU 1 (x, t) + η 2 U 2 (x, t). ( 94
)
Collecting ( 32), ( 42) and (56), and neglecting the higher-order terms, this second-order mean field satisfies

α 0 ∂ 2 tt ū2 -2W 0 ∂ 2 tx ū2 -β 0 ∂ 2 xx ū2 = ηF(ū 2 ) + η 2 E(ū 2 ) + f + ηA(f ) + η 2 B(f ). ( 95
)
Since our main interest is to investigate non-reciprocity for only one varying parameter, we will distinguish this case in the next two subsections. Our interest being mainly on the non-reciprocal behaviour, we will consider f = 0 and look at the associated dispersion relations.

Dispersion relation in the case of ρ being constant

Using ( 34), ( 49), ( 60) and ( 61) for f = 0, (95) reduces to

∂ 2 tt ū2 -β 0 ∂ 2 xx ū2 = η 2 β 1 β 0 R + M ∂ 4 ttxx ū2 + N β ∂ 4 txxx ū2 (96)
which, in dimensionalised coordinates, gives for the dimensionalized mean field Ū2 = ū2 /k

∂ 2 T T Ū2 -β 0 (c ) 2 ∂ 2 XX Ū2 = h 2 β 1 β 0 R + M ∂ 4 T T XX Ū2 + N β c ∂ 4 T XXX Ū2 . ( 97
)
Considering a harmonic solution Ū2 = exp(i(ωT -kX)) yields the dispersion relation

-ω 2 + β 0 (c ) 2 k 2 = β 1 β 0 R + M h 2 k 2 ω 2 -N β c h 2 ωk 3 . ( 98
)
Since N β = 0, the term in k 3 is non zero, which ensures non-reciprocity. For numerical examples later on, we write explicitly the different coefficients in the dispersion relation. Firstly, we know from (69) that β 0 = C ∂P . Moreover, we deduce from P = 0 that

β 0 = β β -c 2 1 β -c 2 -1 . ( 99
)
Secondly, integration of P and P = 0 leads to the following expression of P

P (y) = y 0 β 0 -β(z) β(z) -c 2 dz - 1 0 s 0 β 0 -β(z) β(z) -c 2 dz ds, (100)
and we want to write the right-hand-side terms of (98) in terms of P only. It has already been done for N β in (73), but we still have to compute β( 1 β 0 R + M ) . Integration of (71) together with R = 0 gives

R(y) = y 0 -(β(z) + c 2 )P (z) + C ∂R β(z) -c 2 dz + 1 0 s 0 (β(z) + c 2 )P (z) -C ∂R β(z) -c 2 dz ds, (101)
with C ∂R given by (72). Concerning βM , the reciprocity identity (131) together with (71) yields

βM = 1 - β β 0 P 2 - β β 0 - β + c 2 β -c 2 P + C ∂R β -c 2 P -2c S P + 1 β 0 βRP , ( 102
)
where S is solution of

   ∂ y (β -c 2 )S + cP = -cP
1-periodicity of S, S = 0, and continuity of S and (β -c 2 )S + cP .

(103)

We deduce from this system together with S = 0 that

S = - 2c β -c 2 P + C ∂S β -c 2 with C ∂S = 2c 1 β -c 2 P 1 β -c 2 -1 . ( 104
)
Consequently the dispersion relation ( 98) is given by

-ω 2 + β 0 (c ) 2 k 2 = L(P )h 2 k 2 ω 2 -N β c h 2 ωk 3 , ( 105
)
with β 0 given by (99), N β by (73) and

L(P ) = 1 + 2c 2 β -c 2 β β 0 + 2 P 2 - β β 0 C ∂R + 2cC ∂S 1 β -c 2 P + 1 β 0 βR(P + 1) , ( 106
)
with C ∂R given by (72), C ∂S by (104), R by (101), P by (69) and P by (100).

Dispersion relation in the case of µ being constant

The same approach is followed in this case. It results in the following equation in dimensionalised variables

α 0 ∂ 2 T T Ū2 -(c ) 2 ∂ 2 XX Ū2 = h 2 -α R -cN ∂ 4 T T XX Ū2 + N α c ∂ 4 T T T X Ū2 , (107)
with associated dispersion relation

-α 0 ω 2 + (c ) 2 k 2 = M(Q)h 2 k 2 ω 2 - N α c h 2 ω 3 k, ( 108
)
where

                                       M(Q) = -1 + 2α α 0 1 1-c 2 α + 4 1-c 2 α Q 2 + αR -1+c 2 α 0 1-c 2 α -(cC ∂R α + 2C ∂S ) 1 1-c 2 α Q Q(y) = c y 0 α 0 -α(z) 1-c 2 α(z) dz + C Q and R(y) = - y 0 1+c 2 α(z) cα 0 (1-c 2 α(z)) Q(z)dz + C ∂R y 0 1 1-c 2 α(z) dz + C R C ∂R = 1 cα 0 1 + c 2 α 1 -c 2 α Q 1 1 -c 2 α -1 and C ∂S = 2 1 -c 2 α Q 1 1 -c 2 α -1 α 0 = α 1 -c 2 α 1 1 -c 2 α -1 and C Q = c 1 0 s 0 α(z) -α 0 1 -c 2 α(z) dz ds C R = 1 0 s 0 1 + c 2 α(z) cα 0 (1 -c 2 α(z)) Q(z)dz ds -C ∂R 1 0 s 0 1 1 -c 2 α(z)
dz ds.

(109)

Numerical experiments

In this section, we illustrate on a bilayered example the results obtained through asymptotic developments in the previous section. We first focus on the dispersion diagrams associated with both the space-time modulated laminate and the homogenised models. We then illustrate the results expected through time-domain simulations in the microstructured medium.

Dispersion diagrams

To start with, we consider a bilayered medium of periodicity h = 20 with piecewise constant parameters given by

(ρ, µ)(y) = (ρ A , µ A ) if y ∈ (0, 0.25) (ρ B , µ B ) if y ∈ (0.25, 1), (110)
and a modulation velocity c m . We also choose the reference parameters such that

ρ = ρ and µ = 1/µ -1 , hence c = 1/ ρ 1/µ . ( 111
)
Table 1: Physical parameters for the bilayered medium (110) with both parameters modulated. Mass densities are given in kg • m -3 , shear moduli in Pa, and velocities in m • s -1 .

Fig. number ρ

A ρ B µ A µ B c A c B c B /c A c m Fig. 2a

Leading-order model

Fig. 2 displays the dispersion diagrams associated with this configuration for the leading-order model [START_REF] Martin | Acoustics and dynamic materials[END_REF] previously obtained in the literature [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF]. Different values of the parameters are chosen and summarized in Table 1. In each subfigure, the plain line denotes the dispersion diagram for the case described by Table 1 where both parameters are modulated, and shows asymmetry with respect to k = 0 underlying therefore non-reciprocity. Through the different examples, we can see that nonreciprocity happens in both supersonic and subsonic regimes, and increases as the contrast of the parameters increases. The dotted line denotes the limit case W 0 = 0 which happens as soon as one of the parameters is constant leading therefore to reciprocity. where both parameters are modulated described by Table 1, while the dotted line denotes the case where W 0 = 0. The y-axis is given in terms of the small parameter η = k h = ωh/c .

Second-order model

Now we consider the case where only µ or ρ are modulated and compare the exact dispersion diagrams, and the homogenised one at both leading and second order.

Firstly, the exact dispersion relation for a bilaminate alternating two phases of parameters (ρ i , µ i ) and length h i (i = A, B) and modulated in time in a wave-like fashion, such as in Fig. 1, is given by the zeros of the following dispersion function [START_REF] Lurie | An Introduction to the Mathematical Theory of Dynamic Materials[END_REF][START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF]:

Disp(ω, k) = cos kh -c m (ω -c m k) h A c 2 A -c 2 m + h B c 2 B -c 2 m -cos ω-cmk c 2 A -c 2 m c A h A cos ω-cmk c 2 B -c 2 m c B h B + 1 2 Z A Z B + Z B Z A sin ω-cmk c 2 A -c 2 m c A h A sin ω-cmk c 2 B -c 2 m c B h B , ( 112
)
with c i = µ i /ρ i and Z i = √ ρ i µ i for i = A, B.
Secondly, the homogenised dispersion relation at the leading (or first) order is computed solving for [START_REF] Martin | Acoustics and dynamic materials[END_REF] with W 0 = 0. Thirdly, the homogenised dispersion relation associated to the second-order model developped in this paper is plotted solving either (105) or (108) for the cases α constant and β constant, respectively. The case α = 1 is described by Table 2 and Fig. 3, while the case β = 1 is described by Table 4 and Fig. 6 for different choices of parameter contrasts or modulation velocities. In the figures, the colormap denotes the logarithm of the exact dispersion function (112). Consequently, the dark lines represent its zeros which are solutions of the exact dispersion relation for the time-modulated bilaminate. We recover the existence of asymetric band-gaps in the subsonic regime, and of stop-bands in wavenumber instead of frequency (the so called k-gaps) in the supersonic regime [START_REF] Cassedy | Dispersion relations in time-space periodic media: Part i-stable interactions[END_REF][START_REF] Cassedy | Dispersion relations in time-space periodic media part II-unstable interactions[END_REF][START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF]. The dispersion relation associated to the leading-order (or first-order) model and the second-order model are also plotted with a dashed black line, and a dotted blue line, respectively. The y-axis is given in terms of the small parameter η(ω) = ωh/c to investigate the quality of the approximation as we move away from the assumption η 1 of the asymptotic process. To investigate more closely, the agreement between the dispersion diagrams of the microstructured configuration and of the homogenised models, dispersion diagrams close to the origins are obtained by root finding and plotted in Fig. 4 and7. These figures underline that the second-order homogenised model investigated in this paper allows to get a better approximation of the exact dispersion diagram compared to the leading-order or first-order model, especially as we go further away from the origin. In practice, one notices that, as long as we do not enter a bang-gap or a k-gap, the agreement is even quite good as η approaches 1 while the assumption of the asymptotic process is η 1. For more quantitative measures, we introduce the following relative errors between the microstructured dispersion relation ω micro and the homogenised ones ω homog,i (i = 0, 2 for the leading/first-order and the second-order one, respectively):

i = I (ω micro (k) -ω homog,i (k)) 2 dk I ω micro (k) 2 dk ∈ (0, 1). ( 113
)
Integrations are performed on an interval I defined by

I = {k such that k ∈ [-k lim ; k lim ]}
where k lim is chosen so that both kh π/2 and the interval stops about the entrance of any band-gap or k-gap, for which both homogenised models are no longer accurate, as noticed in Fig. 3 and6. The values of the chosen k lim and of these relative errors are given in the first columns of Tables 3 and5, and confirm the better agreeement of the second-order homogenised model. To validate the developments of the high-order homogenization, we plot the absolute error [ω micro (k) 2 -ω homog,i (k) 2 | as a function of both positive and negative wavenumbers. Fig. 5 represents the results for the very first example of Table 2. We can see that we recover an error of order 4 for the leading-order model, and of order 6 for the second-order model. This validates the developments of the higher-order terms conducted in this paper. This second-order model captures dispersive effects that are missed at the previous orders, and more importantly non-reciprocal behaviour even if only one of the parameters is space-time dependent. To measure this non-reciprocal behaviour more quantitatively we introduce the following ratio:

R = I ω odd (k) 2 dk I ω homog,2 (k) 2 dk ∈ (0, 1) (114)
where ω odd (k) = (ω homog,2 (k) -ω homog,2 (-k))/2 is the odd part of the frequency as a function of the wavenumber. In the reciprocal case, ω(k) is an even function and R is consequently zero. On the contrary, the closer to 1 R is, the more asymmetric the dispersion curve is. These measures of non-reciprocity are given in the last columns of Tables 3 and5. It confirms that the secondorder homogenised model captures non-reciprocal effects: as soon as one parameter is modulated, R is non-zero, while at the leading order it vanishes unless both parameters are modulated. These measurements also give an insight on how reciprocity evolves together with the parameters: it increases as the material contrast increases, and as the modulation velocity approaches the velocity range of the material. R Fig. 3a π/2 7.4761% 2.7001% 3.2335% Fig. 3b and 4a π/2 9.4515% 2.4897% 4.6307% Fig. 3c and4b 0.69 2.7655% 1.9600% 0.6545% Fig. 3d π/2 0.7259% 0.1905% 0.3804% Table 4: Physical parameters for the bilayered medium (110) with only mass density modulated (shear modulus is constant µ A = µ B = 2.66 × 10 9 Pa). Mass densities are given in kg • m -3 .

= ρ B = 1.37 × 10 3 kg • m -3 ). Fig. number µ A (Pa) µ B (Pa) c A (m • s -1 ) c B (m • s -1 ) c B /c A c m (m • s -1 ) Fig. 3a 10

Fig. number ρ

A ρ B c A (m • s -1 ) c B (m • s -1 ) c A /c B c m (m • s -1 ) Fig. 6a & 7a 10

Numerical setup for time-domain simulations

We now focus on time-domain simulations in the microstructured medium. The system ( 14) is solved in the moving frame. For this purpose, a ADER (Arbitrary high order using Derivatives) scheme is 2 used on a uniform grid with mesh ∆x and time step ∆t. This two-time step finite-difference scheme is fourth-order accurate in space and time. It is stable under the Courant-Friedrichs-Lewy condition γ = max c∆t/∆ ≤ 1; in practice, one uses γ = 0.95. Outgoing conditions ensure that no waves are reflected by the edges of the computational domain. An immersed interface method [START_REF] Lombard | How to incorporate the spring-mass conditions in finite difference schemes[END_REF] is implemented to discretize the jump conditions (11) on a Cartesian grid. This method ensures a subcell resolution of interfaces inside the meshing and maintains the accuracy of the scheme despite the non-smoothness of the solution across interfaces. The computational domain [-2500, 1500] is discretized on 4000 grid nodes, hence ∆x = 1. A finite heterogeneous slab with 101 interfaces is considered on [-1000, 0]. It is a bilayered medium of periodicity h = 20 with piecewise constant parameters (110). The parameters in the two media surrounding the slab are the leading-order effective parameters ρ 0 = ρ α 0 and µ 0 = µ β 0 , with α 0 and β 0 defined in [START_REF] Amr | Spatiotemporal light control with active metasurfaces[END_REF], respectively. Doing so minimizes the reflection at the edges of the slab and puts the focus on the interaction with the bulk.

A pulse is emitted by a Dirac source point at X = X s and with a central frequency f c . The time evolution of the source is a smooth combination of sinusoids with bounded support:

s(t) = 4 m=1 a m sin(b m ω c t) if 0 < t < 1 f c , 0 otherwise, (115)
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, ensuring C 6 smoothness (differentiable 6 times and the 6-th derivative is continuous). The typical small parameter (16) for this type of signal is chosen to be the one associated with the central frequency i.e.

η c = 2πhf c /c , ( 116
)
with c given by (111).

In the moving frame, this source moves at speed -c m . Similarly, a receiver is considered at X = X r ; it moves at speed -c m and records the fields V h and Σ h at each time step. Numerical integration provides the displacement U h . Lastly, the simulations are stopped before the time where the source or receiver exits from the computational domain or enters into the slab.

Numerical results

First we consider the case where both parameters are modulated given by the following values of the parameters in (110): 10

(ρ A , µ A ) = (
and the results of the corresponding time-domain simulations are given by Fig. 11. In both figures, we observe a greater dispersive behaviour for the left-going waves (the ones emited by the source on the right), which is consistent with the dispersion diagrams of Fig. 3 and6.

Conclusion

In the present paper, we performed high-order homogenisation of the 1D wave equation when the physical parameters are modulated in space and time in a wave-like fashion. Approximating fields up to the second order allows to take into account dispersive effects, and most importantly, non-reciprocity that appears at this order does not vanish as we vary only one of the parameters. We therefore proved that it is enough to modulate a single parameter in order to achieve low-frequency non-reciprocity for the acoustic wave equation. Similar properties then hold for the corresponding equations in the case of transverse electric (138) and transverse magnetic (142) waves. The natural choice for SH waves is to modulate in space and time the shear modulus µ [START_REF] Danas | Experiments and modeling of iron-particlefilled magnetorheological elastomers[END_REF] and keep constant the mass density ρ. For transverse electromagnetic waves, it is to time-modulate the relative permittivity ε and assume the relative permeability µ = 1 (i.e. no magnetism) in ( 138) and (142). Thus, non-reciprocity does not require modulation of both parameters as was assumed in elasticity and electromagnetism [START_REF] Nassar | Modulated phononic crystals: Non-reciprocal wave propagation and willis materials[END_REF][START_REF] Huidobro | Homogenization theory of spacetime metamaterials[END_REF]. We hope this will foster efforts towards an experimental demonstration, e.g. of the Fresnel drag effect in time-modulated dielectric layers [START_REF] Huidobro | Fresnel drag in space-time-modulated metamaterials[END_REF]. Some theoretical aspects have been disregarded in this paper, such as the boundary effects in the case of a finite medium [START_REF] Beneteau | Modèles homogénéisés enrichis en présence de bords : Analyse et traitement numérique[END_REF][START_REF] Vinoles | Problèmes d'interface en présence de métamatériaux : modélisation, analyse et simulations[END_REF][START_REF] Maurel | Sensitivity of a dielectric layered structure on a scale below the periodicity: A fully local homogenized model[END_REF][START_REF] Cornaggia | A homogenized model accounting for dispersion, interfaces and source points for transient waves in 1d periodic media[END_REF], or the behaviour for long times [START_REF] Lamacz | Dispersive effective models for waves in heterogeneous media[END_REF][START_REF] Allaire | Crime pays: homogenization for long times[END_REF]; they could be investigated in future works. Field pattern space-time composites with PT-symmetry such as chessboards as introduced by Milton and Mattei [START_REF] Mattei | Field patterns: a new mathematical object[END_REF] also represent a further challenge. We could further explore optimal design of time-modulated media [START_REF] Cherkaev | Optimal structures for focusing and energy accumulation: mathematical models and intuition[END_REF]. As another perspective, our asymptotic analysis can be applied to time-modulated conductivity, hydrodynamic and elasticity equations (we note in passing the fascinating work [START_REF] Ab Movchan | Frontal waves and transmissions for temporal laminates and imperfect chiral interfaces[END_REF] on temporal elastic laminates with imperfect chiral interfaces), as well as other governing equations in physics described by linear partial differential equations, where we expect that similar non-reciprocal effects can be unveiled in the homogenisation regime. Another perspective concerns the case of locally-resonant media [START_REF] Felbacq | Theory of mesoscopic magnetism in photonic crystals[END_REF][START_REF] Auriault | Long wavelength inner-resonance cut-off frequencies in elastic composite materials[END_REF], to combine both low-frequency non-reciprocity and low-frequency resonances for more important macroscopic effects. Finally, high-frequency homogenisation [START_REF] Craster | High-frequency homogenization for periodic media[END_REF] could be applied to get effective properties around the opening of stop-bands either in frequency or in wavenumber.

A Reciprocity identities

The cell problems introduced to express u 1 , u 2 and u 3 lead to a family of reciprocity identities that are listed here and extensively used in the homogenisation process.

Integration by parts of (26

) × P -(25) × Q leads to βQ = c αP (120)
Integration by parts of (37

) × P -(25) × R leads to βR = -βP + β 0 α 0 α(1 -cQ )P + β 0 α 0 c αQP (121)
Integration by parts of (37

) × Q -(26) × R leads to c αR = -βQ + βP Q + βQP + β 0 α 0 αQ (122)
Integration by parts of (38) × P -(25) × S leads to

βS = -βQ P + βQP + 2W 0 α 0 α(1 -cQ )P + cαQP (123)
Integration by parts of (26

) × S -(38) × Q leads to c αS = -c αP Q + c αP Q + 2W 0 α 0 αQ (124)
Integration by parts of (25) × B -(39) × P leads to

βB = 1 α 0 αP + c α(QP -P Q ) (125)
Integration by parts of (26

) × B -(39) × Q leads to c αB = 1 α 0 αQ (126)
Integration by parts of (37) × S -(38) × R leads to

-β(1 + P )S + βP S + β 0 α 0 α(1 -cQ )S + β 0 α 0 c αQS = -(αcP + βQ )R + (cαP + βQ)R + 2W 0 α 0 α(1 -cQ )R + 2W 0 α 0 c αQR

B Transverse electromagnetic waves

We now consider the two-scale homogenisation of the vector Maxwell system of a time-modulated layered medium. Let us consider a fixed space-time Cartesian coordinate system (x, t) = (x 1 , x 2 , x 3 , t) and a time modulated layered periodic medium. The propagation is along the direction x 1 of stacking of layers. In what follows, the subscript denotes dependence of the field upon the periodicity η in the space-time variable x 1 -c 1 t, where c 1 is the modulation speed along x 1 .

We first consider the Maxwell-Faraday equation that relates the magnetic inductance field B η and the electric field

E η via curl E η + ∂ ∂t B η = 0. ( 132
)
This equation implies that ∇ • B η (x, t) = 0, from which it follows that there is a vector magnetic potential field such that B η (x, t) = curl A η (x, t).

(133)

We use this vector potential as an unknown rather than the electric and magnetic fields. Combining (132) and ( 133), it follows that E η + ∂ ∂t A η is curl free and thus there exists a scalar electric potential φ η such that

E η (x, t) = - ∂ ∂t A η (x, t) -∇φ η (x, t). ( 134
)
Besides, the magnetic field H η and the electric displacement field D η satisfy Maxwell-Ampère's equation

curl H η - ∂ ∂t D η = 0 . ( 135
)
Furthermore, the electric displacement D η (resp. magnetic inductance B η ) is related to the electric field E η (resp. magnetic field H η ) via

D η = ε 0 ε x 1 -c 1 t η E η , B η = µ 0 µ x 1 -c 1 t η H η , ( 136
)
where ε 0 µ 0 = 1/c 2 0 , c 0 being the speed of light in vacuum, and where ε and µ respectively denote the relative permittivity and permeability. These are 1 periodic functions of x 1 and 1/c 1 periodic functions of t.

From (133)-(136), we deduce that the magnetic potential A η and the electric potential φ η satisfy curl µ

-1 x 1 -c 1 t η curl A η = - 1 c 2 0 ∂ ∂t ε x 1 -c 1 t η ∂ ∂t A η - 1 c 2 0 ∂ ∂t ε x 1 -c 1 t η ∇φ η (137)
Let us further assume that the electric field is orthogonal to the plane (x 1 x 2 ) of the layering (transverse electric, TE, case). Eq. (134) implies that (0, 0, E η (x 1 , t)) T = -( ∂ ∂x 1 φ η (x 1 , t), ∂ ∂x 2 φ η (x 1 , t), ∂ ∂t A η (x 1 , t)+ ∂ ∂x 3 φ η (x 1 , t)) T = -(0, 0, ∂ ∂t A η (x 1 , t)) T , where A η is the only non-zero component of A η , and so ∇φ η must be a null vector field in (137).

We deduce that in TE polarisation and for the considered layering, (137) reduces to

∂ ∂x 1 µ -1 x 1 -c 1 t η ∂ ∂x 1 A η = 1 c 2 0 ∂ ∂t ε x 1 -c 1 t η ∂ ∂t A η . ( 138
)
We note that ( 17) and ( 138) are identical upon the replacement of α by µ -1 and of β by c -2 0 ε. Maxwell's equations are fully symmetric, so we can now interchange the roles of the vector magnetic potential A η and the scalar electric potential φ η with a vector electric potential V η and a scalar electric potential ψ. More precisely, (135) implies that ∇ • D η (x, t) = 0, from which it follows that there is a vector magnetic potential field such that D η (x, t) = curl V η (x, t).

(139)

Combining ( 135) and (139), we know there exists a scalar magnetic potential ψ η such that

H η (x, t) = - ∂ ∂t V η (x, t) -∇ψ η (x, t). ( 140
)
From ( 132), (136), ( 139) and (140), we deduce that V η and ψ η satisfy

curl ε -1 x 1 -c 1 t η curl V η = - 1 c 2 0 ∂ ∂t µ x 1 -c 1 t η ∂ ∂t V η - 1 c 2 0 ∂ ∂t µ x 1 -c 1 t η ∇ψ η ( 141
)
Let us now assume that the magnetic field is orthogonal to the plane (x 1 x 2 ) of the layering (transverse magnetic, TM, case). From Eq. (139) we deduce that ∇ψ η is a null vector field in (141).

Therefore, in TM polarisation and for the considered layering, (141) reduces to

∂ ∂x 1 ε -1 x 1 -c 1 t η ∂ ∂x 1 V η = 1 c 2 0 ∂ ∂t µ x 1 -c 1 t η ∂ ∂t V η (142)
We note that ( 17) and (142) are identical upon the replacement of α by ε -1 and of β by c -2 0 µ. Therefore, it is enough to modulate the permittivity ε in Eq. (138) and Eq. (142) to achieve some non-reciprocity in the homogenisation frequency regime in TE and TM polarizations. Our work thus complements leading-order homogenisation results in [START_REF] Huidobro | Fresnel drag in space-time-modulated metamaterials[END_REF][START_REF] Huidobro | Homogenization theory of spacetime metamaterials[END_REF] wherein space-time modulation of both permittivity and permeability was a prerequisite for non-reciprocity.

10 3 1 . 5 × 2 6 ×

 1526 Fig. number ρA ρ B µ A µ B c A c B c B /c A c m Fig. 2a10 3 1.5 × 10 3 10 9 6 × 10 9 10 3 2 × 10 3 2 5 × 10 2 Fig. 2b 10 3 1.5 × 10 3 10 9 6 × 10 9 10 3 2 × 10 3 2 3 × 10 3 Fig. 2c 10 3 10 2 10 9 6 × 10 9 10 3 7.74 × 10 3 7.74 5 × 10 2

Figure 2 :

 2 Figure 2: Dispersion diagrams for the leading-order model. The plain line denotes a configuration where both parameters are modulated described by Table1, while the dotted line denotes the case where W 0 = 0. The y-axis is given in terms of the small parameter η = k h = ωh/c .

Fig. number k lim 0 2

 2 Fig. numberk lim

Figure 3 :

 3 Figure 3: Dispersion diagrams when modulating only µ. The colormap denotes the logarithm of the dispersion function for the microstructured configuration, the dark lines therefore represent the exact dispersion diagram. The dotted blue line denotes the dispersion diagram obtained at the second order in (105). The dashed black line denotes the dispersion diagram obtained at the leading [11] or first order. Parameters are given in Table 2

Figure 4 :

 4 Figure 4: Comparison of dispersion diagrams when modulating only µ. Parameters are given in Table 2

Figure 5 :

 5 Figure 5: Errors (in a log-log scale) made on the dispersion curves for the physical parameters of Fig. 3a. Plain lines denote the error |ω 2 micro -ω 2 homog,i | between the frequency obtained by Bloch-Floquet analysis and the one obtained through homogenization at order i. Dashed lines denote reference orders of convergence.

Figure 6 :

 6 Figure 6: Dispersion diagrams when modulating only ρ. The colormap denotes the logarithm of the dispersion function for the microstructured configuration, the dark lines therefore represent the exact dispersion diagram. The dotted blue line denotes the dispersion diagram obtained at the second order in (108). The dashed black line denotes the dispersion diagram obtained at the leading [11] or first order. Parameters are given in Table 4

Figure 7 :

 7 Figure 7: Comparison of dispersion diagrams when modulating only ρ. Parameters are given in Table 4

  (a) Source on the left (t = 0.15 s) (b) Source on the right (t = 0.20 s) Source on the left (t = 0.76 s) (d) Source on the right (t = 1.00 s) Source on the left (t = 2.50 s) (f) Source on the right (t = 2.05 s)

Figure 8 :

 8 Figure 8: Snapshots of U h at different times, with source such as η c = 0.54 (116) on the left (a-c-e) or on the right (b-d-f) of the modulated slab (grey zone). Simulations are performed in the moving frame, with a left-going source (red dotted line) and a left-going receiver (green dotted line).

Fig. number k lim 0 2 RFig. 6a and 7a 1 7 .

 27 Fig. numberk lim

  3 kg • m -3 , 10 9 Pa) and (ρ B , µ B ) = (1.5 × 10 3 kg • m -3 , 6 × 10 9 Pa),(117) with the modulation velocity being c m = 5×10 2 m•s -1 < min(c A , c B ), which corresponds to a subsonic case. Fig.8displays snapshots of U h for f c = 6Hz, i.e. η c = 0.54 (116), at three different times. In the left column (a-c-e), the source is on the left of the slab, and the receiver is on the right. In the right column (b-d-f), the situation is reversed. Qualitatively, one observes the difference in scattered fields in both cases. Fig.9displays the time evolution of the recorded U h , for various central frequencies f c . The differences by exchanging source and receiver are clearly observed. Fig.10illustrates the case where only µ h is modulated with a velocity c m = 5 × 10 2 m • s -1 . The parameters of the bilaminate for this case are given by:(ρ A , µ A ) = (1.375 × 10 3 kg • m -3 , 1.375 × 10 9 Pa), (ρ B , µ B ) = (1.375 × 10 3 kg • m -3 ,5.5 × 10 9 Pa). (118) In this case, the differences obtained by exchanging source and receiver are still clearly observed, which is a signature of non-reciprocity. Similar observations have been done by modulating only ρ h following (ρ A , µ A ) = (10 3 kg • m -3 , 6 × 10 9 Pa) and (ρ B , µ B ) = (3 × 10 3 kg • m -3 , 6 × 10 9 Pa),

  η c = 0.27 (f c =3 Hz) η c = 0.54 (f c =6 Hz) c = 0.81 (f c =9 Hz) η c = 1.08 (f c =12 Hz)

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Modulation of both parameters. Time history of U h measured at the receiver, for different small parameters η c (116). Plain and dotted lines denote the case of a source on the left and on the right of the slab, respectively.

( 127 )

 127 Integration by parts of (52) × Q -(26) × N leads toc αN = ( β 0 α 0 α -β)Q 2 + c α(RQ -R Q) + β(SQ -S Q) + 2W 0 β 0 βP Q + β(R Q -RQ ) (128)Integration by parts of (52) × P -(25) × N leads to βN = ( β 0 α 0 α -β)P Q + c α(RP -R P ) + β(SP -S P ) + 2W 0 β 0 βP 2 + β(R P -RP ) (129)Integration by parts of (26) × M -(51) × Q leads to c αM = αP Q + c α(SQ -S Q) + 2W 0 α 0 αQ 2 -α 0 β 0 βP Q + β(R Q -RQ )(130)Integration by parts of (25) × M -(51) × P leads to βM = (αα 0 β 0 β)P 2 + c α(SP -S P ) + 2W 0 α 0 αP Q + α 0 β 0 β(RP -R P )(131)

Table 2 :

 2 Physical parameters for the bilayered medium (110) with only shear modulus modulated (mass density is constant ρ A

Table 3 :

 3 Case of a bilayered medium (110) with only shear modulus modulated (physical parameters are given in Table2). Relative errors (113) and measure of non-reciprocity (114).

		9	6 × 10 9	8.53 × 10 2	2.09 × 10 3	2.45	5 × 10 2
	Fig. 3b and 4a	10 9	1.5 × 10 10	8.53 × 10 2	3.31 × 10 3	3.87	5 × 10 2
	Fig. 3c and 4b	10 9	6 × 10 9	8.53 × 10 2	2.09 × 10 3	2.45	2.5 × 10 3
	Fig. 3d	10 9	6 × 10 9	8.53 × 10 2	2.09 × 10 3	2.45	3.5 × 10 3

Table 5 :

 5 Case of a bilayered medium (110) with only mass density modulated (physical parameters are given in Table4). Relative errors (113) and measure of non-reciprocity (114).

Acknowledgments

SG and RVC were funded by UK Research and Innovation (UKRI) under the UK government's Horizon Europe funding guarantee [grant number 10033143]. MT and RA would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Mathematical theory and applications of multiple wave scattering where work on this paper was undertaken.