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

HAL
open science

# Diffraction grating with space-time modulation

Kim Pham, Agnès Maurel

► To cite this version:

Kim Pham, Agnès Maurel. Diffraction grating with space-time modulation. Journal of Computational Physics, 2022, 469, pp.111528. 10.1016/j.jcp.2022.111528. hal-03848638

HAL Id: hal-03848638

https://hal.science/hal-03848638v1

Submitted on 13 Oct 2025

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.

Public Domain
# Diffraction grating with space-time modulation

Kim Pham

IMSIA, CNRS, EDF, CEA, ENSTA Paris, Institut Polytechnique de Paris,
828 Bd des Maréchaux, 91732 Palaiseau, France

Agnès Maurel

Institut Langevin, ESPCI Paris, Université PSL, CNRS,
1 rue Jussieu, 75005 Paris, France

## Abstract

We present a theoretical and numerical analysis of the diffraction of acoustic waves by space-time modulated gratings with rigid-type modulations. This is done by deriving the form of the modes which are exact, uncoupled, solutions of the problem in the unbounded regions, inside and outside the grating. The dispersion of the modes is studied as a function of the ratio of the modulation speed to the speed of sound which shows that each spatial diffraction order is associated with a single temporal diffraction order. For a grating of finite extent, the solution is obtained as a superposition of these modes, which couple at the grating interfaces. This provides a numerical, multimodal, method when considering a truncated version of the solution. We provide analysis of the solutions in the harmonic and in the transient regimes.

**Keywords:** Space-time modulated grating. Multimodal method.

## 1. Introduction

Starting with the pioneering works of Simon [1] and of Oliner and co-workers [2, 3], the interaction of waves with media whose properties are space-time dependent has aroused deep interest in the past years since they offer additional degrees of freedom for the control of sound propagation when compared to their conventional space modulated counterparts. This includes non-reciprocal propagation which is similar although non identical to the non-reciprocity in the presence of a flow [4, 5, 6, 7, 8, 9]. Another remarkable feature of space-time modulated media is their capacity to generate frequency conversions to almost arbitrarily high frequencies thanks to the appearance of space-time diffraction orders dictated by the speed of the modulation [10, 11, 12, 13]. Based on these

<footnote>*Fully documented templates are available in the elsarticle package on CTAN.</footnote>
properties, a plethora of applications have been proposed within few years,
including experimental realizations [5, 9, 11, 12, 14, 15].

Besides fabrication challenges, the numerical implementation of these problems involving several spatial and temporal scales which can be very different, is an issue as well. Indeed, the commercial solvers are not yet performant to do that and, if in principle standard methods such as the finite-difference time-domain FDTD can be employed, the simulations can become prohibitively expensive. An exception is the FDTD schemes developed with success in [16, 17, 18] owing to the simplifications which arise when the metasurface is treated as an idealized spatial discontinuity. To circumvent these computational limitations, multimodal methods, already used for classical patterned devices, have been extended to devices with space and time modulations [19, 20, 21, 22, 23]. Multimodal methods use the discrete set of modes which are the exact solutions in each unbounded region, say the unbounded grating and the unbounded free space. The solution of diffraction problems is then written as a superposition of these modes which propagate without coupling in each sub-region and couple only at their interface. For periodic structures, the reduction of the computational domain to a single period and the absence for requirement of any spatial meshing often enable multimodal methods to outperform other numerical solvers.

In this study, we consider a diffraction grating made of a periodic arrangement of sound-rigid and sound-soft layers modulated along its interfaces at velocity $v$ (figure 1). We aim to understand the dispersion of the fundamental solutions, or modes, inside and outside the grating and to provide semi-analytical solutions of scattering problems thanks to the aforementioned multimodal analysis. In §2, we derive the dispersion of the modes, being propagating or evanescent, and we show that our diffraction grating produces spatial diffraction orders

<figure><img src="image_2.png" /><figcaption>Figure 1: Space-time modulated grating of length $\ell$. The grating consists in a $h$-periodic lamination alternating air and a rigid parts modulated in time at constant speed $v$ in direction $e_y$.</figcaption></figure>
associated with up-converted frequencies. The numerical method based on this multimodal analysis is detailed in §3. Solutions of scattering problems for incident plane waves and incident beams in the harmonic regime are given and commented (a script is provided). Eventually in §4, we analyse the temporal and spatial separations of the diffraction orders (or modes) for an incident narrow Gaussian pulse in the transient regime. We notice that this study is closely related to that of [24] which provide a comparison between such a multimodal method and FDTD. In this reference, the authors consider a similar grating with sinusoidal variations of the permittivity, being the electromagnetic equivalent of the bulk modulus in acoustics. The sound-rigid layers can be thought as high values of the mass density and bulk modulus and we consider the limiting case of infinite values.

## 2. Analysis of the modes inside and outside the space-time modulated grating

We consider the linearized Euler equations for acoustic waves

$$ \rho \frac{\partial \mathbf{u}}{\partial t} = -\nabla p, \quad \mathrm{div} \mathbf{u} + \frac{1}{\rho c^2} \frac{\partial p}{\partial t} = 0, \qquad (1) $$

with $\mathbf{u}$ and $p$ the acoustic pressure and velocity, $\rho$ the mass density and $c$ the sound speed (the bulk modulus is $\rho c^2$). Following [25], we introduce the velocity potential $\varphi$ so that (1) can be written

$$ \left( \begin{array}{c} \mathbf{u} \\ \frac{p}{Z} \end{array} \right) = \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right) \left( \begin{array}{c} \nabla \varphi \\ -\frac{1}{c} \frac{\partial \varphi}{\partial t} \end{array} \right), \qquad (2) $$

with $I$ the 2x2 identity matrix and $Z = \rho c$, along with

$$ \mathrm{div} \mathbf{u} + \frac{1}{Zc} \frac{\partial p}{\partial t} = 0. \qquad (3) $$

Eventually, at the moving rigid boundaries (with velocity $v$), we have

$$ Zu_y = Mp, \quad y = \{vt, \xi h + vt\}, \qquad (4) $$

with $M = v/c$ (see Appendix A).

### 2.1. Modes in the moving frame

The analysis of the problem is performed in the moving frame owing to the change of variables $(x, y, t) \rightarrow (x, y, t)$,

$$ \begin{aligned} & x = x, \quad y = y - vt, \quad t = t, \\ & \phi(x, y, t) = \varphi(x, y, t), \end{aligned} \qquad (5) $$
(the same for $(\mathbf{u}, \mathbf{p}) \to (\mathbf{u}, \mathbf{p})$). Eqs. (2)-(3) become

$$ egin{cases} 
\mathbf{u} = \nabla \phi, & \frac{\mathbf{p}}{Z} = -\frac{1}{c} \left( \frac{\partial \phi}{\partial t} - v \frac{\partial \phi}{\partial y} \right), \\ 
\operatorname{div} \mathbf{u} + \frac{1}{Zc} \left( \frac{\partial \mathbf{p}}{\partial t} - v \frac{\partial \mathbf{p}}{\partial y} \right) = 0. 
\end{cases} 
\tag{6} $$

Accordingly, the wave equation takes the form

$$ \frac{\partial^2 \phi}{\partial x^2} + (1 - M^2) \frac{\partial^2 \phi}{\partial y^2} + 2M \frac{1}{c} \frac{\partial^2 \phi}{\partial y \partial t} - \frac{1}{c^2} \frac{\partial^2 \phi}{\partial t^2} = 0. \tag{7} $$

In the harmonic regime, we consider an incident wave at frequency $\omega_0$ in the laboratory frame with $(\alpha_0, \beta_0)$ wavevector and $k_0 = \sqrt{\alpha_0^2 + \beta_0^2} = \frac{\omega_0}{c}$, $\phi^{\text{inc}}(x, y, t) = e^{-i\omega_0 t} e^{i(\alpha_0 x + \beta_0 y)}$. In the moving frame, the incident wave reads

$$ \phi^{\text{inc}}(x, y, t) = e^{-i(\omega_0 - \beta_0 v)t} e^{i(\alpha_0 x + \beta_0 y)}, \tag{8} $$

and we look for solutions of (7) at the shifted frequency

$$ \phi(x, y, t) = e^{-i\hat{\omega}_0 t} \hat{\phi}(x, y), \quad \hat{\omega}_0 = \omega_0 - v\beta_0. \tag{9} $$

To that end, we look for the set of elementary solutions, or modes, labeled by $n$, namely $\phi(x, y, t) = \phi_n(x, y, t)$ with

$$ \phi_n(x, y, t) = e^{-i\hat{\omega}_0 t} \hat{\phi}_n(x, y), \tag{10} $$

that we shall determine outside and inside the modulated grating. We notice that they satisfy exactly (7) in the (shifted) harmonic regime, if they are solutions to

$$ \frac{\partial^2 \hat{\phi}_n}{\partial x^2} + (1 - M^2) \frac{\partial^2 \hat{\phi}_n}{\partial y^2} - 2iM \frac{\hat{\omega}_0}{c} \frac{\partial \hat{\phi}_n}{\partial y} + \frac{\hat{\omega}_0^2}{c^2} \hat{\phi}_n = 0. \tag{11} $$

### 2.1.1. Modes outside the grating

Outside the modulated grating, the pseudo-periodic condition, imposed by the incident wave, applies along $y$. Accordingly, the transverse dependence of $\hat{\phi}_n$ has to be sought in the form

$$ \Phi_n(y) = \frac{1}{\sqrt{h}} e^{i\beta_n y}, \quad \beta_n = \beta_0 + \frac{2n\pi}{h}, \tag{12} $$

with $n \in \mathbb{Z}$. The associated mode is determined by inserting $\hat{\phi}_n(x, y) = f_n(x)\Phi_n(y)$ in (11), from which $f_n$ is solution to $f_n'' + \alpha_n^2 f_n = 0$ with $\alpha_n^2 + \beta_n^2 = \left(\frac{\hat{\omega}_0}{c} + M\beta_n\right)^2$ and with $\hat{\omega}_0$ in (9). Hence the modes

$$ \left\{ \begin{array}{l} \hat{\phi}_n(x, y) = e^{i\alpha_n x} \Phi_n(y), \\ \alpha_n = \sqrt{k_n^2 - \beta_n^2}, \quad k_n = k_0 + M \frac{2n\pi}{h}, \end{array} \right. \tag{13} $$
are exact solutions of the problem outside the grating associated with the $\Phi_n$ forming a basis of pseudo-periodic functions satisfying the orthogonality relation

$$ \int_{0}^{h} \Phi_{n}(y) \Phi_{m}^{*}(y) dy = \delta_{mn}, \qquad (14) $$

where * denotes the complex conjugate. We notice that the modes with $k_n > \beta_n$ correspond to right-going propagating modes (from $-\infty$ to $+\infty$ along x) and left-going modes are deduced using $\alpha_n \to -\alpha_n$; modes with $k_n \le \beta_n$ are evanescent modes vanishing when $x \to +\infty$ with the convention on the root $\alpha_n = +i\tilde{\alpha}_n$ with $\tilde{\alpha}_n$ real positive. We shall use these properties to write the reflected and transmitted waves in the forthcoming equation (27).

### 2.1.2. Modes inside the grating

The boundary condition in (4) now reads $Zu_y = Mp$, at $y = \{0, \xi h\}$. With $(p, u_y)$ in (6) we look for modes $\hat{\phi}_n$ satisfying the boundary conditions

$$ (1 - M^2) \frac{\partial \hat{\phi}_n}{\partial y} = i \frac{\hat{\omega}_0}{c} M \hat{\phi}_n, \quad y = \{0, \xi h\}, \qquad (15) $$

from which they can be thought in the form $\hat{\phi}_n(x, y) = g_n(x)\Psi_n(y)$ with

$$ \Psi_n(y) = \sqrt{\frac{2 - \delta_{n0}}{\xi h}} e^{i \frac{M(k_0 - M\beta_0)}{1 - M^2} y} \cos\left(\frac{n\pi}{\xi h} y\right), \quad n \in \mathbb{N}, \qquad (16) $$

which satisfy exactly (15). Inserting the form of $\hat{\phi}_n$ in (11), we obtain that $g_n(x)$ satisfies $g_n'' + \gamma_n^2 g_n = 0$, hence the modes

$$ \left\{ \begin{array}{l} \hat{\phi}_n(x, y) = e^{i\gamma_n x} \Psi_n(y), \\ \gamma_n = \sqrt{\frac{(k_0 - M\beta_0)^2}{1 - M^2} - (1 - M^2)\left(\frac{n\pi}{\xi h}\right)^2}, \end{array} \right. \qquad (17) $$

are exact solutions inside the grating associated with the basis of transverse functions $\Psi_n(y)$ satisfying the orthogonality relation

$$ \int_{0}^{\xi h} \Psi_{n}(y) \Psi_{m}^{*}(y) dy = \delta_{mn}. \qquad (18) $$
## 2.2. The modes in the laboratory frame

Owing to the previous analyses, we now express the discrete sets of modes in the laboratory frame.

### 2.2.1. Dispersion outside the modulated grating

Outside the modulated grating, we use (12)-(13) along with (5) which results in the forms of the modes

$$ \varphi_n(x, y, t) = e^{-i\omega_n t} e^{i(\alpha_n x + \beta_n y)}, \qquad (19) $$

with

$$ \omega_n = \omega_0 + Mc\frac{2n\pi}{h}, \quad n \in \mathbb{Z}, \qquad (20) $$

$$ \beta_n = \beta_0 + \frac{2n\pi}{h}, \quad \alpha_n^2 + \beta_n^2 = k_n^2, \quad k_n = \frac{\omega_n}{c}, $$

(for simplicity, we have omitted the normalization factor $1/\sqrt{h}$ in (12)). We show in figure 2 the typical dispersion of the propagating modes ($\alpha_n$ real) outside the grating by means of ($\alpha_n, \beta_n$) against $\omega_0$ ($M = 0.8$ and $\xi = 0.8$). For $M = 0$, $\alpha_n^2 + \beta_n^2 = k_0^2 = (\omega_0/c)^2$, the dispersion of all the modes would fall on the same cone of radius $\omega_0/c$ (of apex at $(0,0)$) in the plane $(\alpha_n, \beta_n)$. In contrast for $M \neq 0$, the propagating modes describe only a part of the cones of radius $k_n > 0$ (and apex at $(0,0)$). We notice that in the subwavelength regime ($k_0h/2\pi < 1$), the condition of real $\alpha_n$ in (20) imposes that for $M > 0$, the mode $n$ is propagating if

$$ -\frac{1+s}{1+M} \frac{k_0 h}{2\pi} \le n \le \frac{1-s}{1-M} \frac{k_0 h}{2\pi}, \qquad (21) $$

with $s = \beta_0/k_0 < 1$. The modulation fosters higher diffraction orders $n > 0$ up to arbitrarily high values when $M \to 1$, that is associated with up-shifted

<figure><img src="image_3.png" /><figcaption>Figure 2: Dispersion of the propagating waves outside the modulated grating ($M = 0.8, \xi = 0.8$).</figcaption></figure>
frequencies $\omega_n > \omega_0$ from (20). In contrast, it prevents from the appearance of diffraction orders $n < 0$ associated with down-shifted frequencies $\omega_n < \omega_0$ (with $\omega_n > 0$ since, from (21), for $n < 0$, we have $\omega_n \ge -\frac{1-Ms}{1+s} \frac{2n\pi c}{h} > 0$). For instance, $n = -1$ becomes propagating for $k_0h/2\pi \ge (1+M)/(1+s)$, that is for a higher frequency than in the absence of modulation. Note that the result holds for $M < 0$ since the problem is invariant for $(M, n, \beta_0) \to -(M, n, \beta_0)$.

### 2.2.2. Dispersion inside the modulated grating

Inside the modulated grating, we use (16)-(17) along with (5) which results in the forms of the modes

$$ \varphi_n(x, y, t) = e^{-i\omega_0(1-Ms)t} e^{i\gamma_n x} \Psi_n(y-vt), \quad (22) $$

where $s = \beta_0/k_0$. To get insight of the dispersion in the grating, we write the mode $n$ as the superposition of two waves, namely

$$ \varphi_n(x, y, t) = (e^{-i\Omega_n t} e^{i\delta_n y} + e^{-i\Omega_{-n} t} e^{i\delta_{-n} y}) e^{i\gamma_n x}, \quad (23) $$

where, extending $n$ to $\mathbb{Z}$, we have

$$ \Omega_n = \Omega_0 + Mc\frac{n\pi}{\xi h}, \quad \Omega_0 = \omega_0\frac{1-Ms}{1-M^2}, \quad (24) $$

$$ \delta_n = M K_0 + \frac{n\pi}{\xi h}, \quad \gamma_n^2 + \delta_n^2 = K_n^2, \quad K_n = \frac{\Omega_n}{c}, $$

see figure 3. Note that the above expression of $\gamma_n$ is consistent with that given in (17), namely we have $\gamma_n^2 = (1-M^2)(K_0^2 - (n\pi/\xi h)^2)$, hence $\gamma_{-n} = \gamma_n$. As it should be, for $M=0$, we recover

$$ \varphi_n^{M=0}(x, y, t) = 2e^{-i\omega_0 t} e^{i\gamma_n x} \cos\left(\frac{n\pi}{\xi h} y\right), \quad (25) $$

with $\gamma_n = \sqrt{k_0^2 - \left(\frac{n\pi}{\xi h}\right)^2}$ providing cut-on frequencies for $\omega_0 \ge c\frac{n\pi}{\xi h}$. When $M \ne 0$, (24) is close, although not identical, to the inverse Doppler effect reported in [26] for $s \to 1$ which corresponds to a wave propagating along the modulation. Indeed the fundamental frequency $\Omega_0 = \omega_0/(1+M)$ is down-shifted for $M > 0$ corresponding to the usual Doppler effect for receding walls. In contrast higher harmonics $n$ with $n\pi(1+M) > \xi k_0h$ are up-shifted which correspond to inverse Doppler effects.

The propagating modes inside the modulated grating correspond to $\gamma_n$ real. At frequency $\omega_0$, $(\gamma_n, \delta_n)$ describe the hyperbola of semi major/minor axis $\frac{n\pi}{\xi h}(M, \sqrt{1-M^2})$ and center $(\gamma_n, \delta_n) = \frac{n\pi}{\xi h}(1, 0)$. (specifically, the hyperbola is $\frac{(X-1)^2}{M^2} - \frac{Y^2}{1-M^2} = 1$ with $\delta_n = \tau X, \gamma_n = \tau Y$ and $\tau = n\pi/(\xi h)$). If we account for the constraint $s \in (-1, 1)$, $\Omega_0$ is bounded, which imposes that $\delta_n$ describes only the upper part of the hyperbola. Namely for $M \ge 0$ we have

$$ (1+M)\frac{n\pi}{\xi h} \le \delta_n \le \frac{Mk_0}{1-M} + \frac{n\pi}{\xi h}. \quad (26) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Dispersion of the propagating waves inside the modulated grating.</figcaption></figure>

### 3. Numerical implementation using multimodal method

We work in the laboratory frame where $\varphi(x, y, t)$ can be expressed as the superposition of the modes (19) for $x \in (-\infty, 0)$ and for $x \in (\ell, +\infty)$ and a superposition of the mode (22) for $x \in (0, \ell)$. We shall consider wave sources which can be decomposed as a superposition of plane waves and we consider below a plane wave with wave-vector $(\alpha_0, \beta_0)$ at frequency $\omega_0$ satisfying (20). Accordingly, the solution inside and outside the grating reads

$$ \varphi = \begin{cases} e^{-i\omega_0 t} e^{i\alpha_0 x} \Phi_0(y) + \sum_{n \in \mathbb{Z}} r_n e^{-i\omega_n t} e^{-i\alpha_n x} \Phi_n(y), & x \in (-\infty, 0) \\ e^{-i\hat{\omega}_0 t} \sum_{n \in \mathbb{N}} (a_n^+ e^{i\gamma_n x} + \bar{a}_n e^{-i\gamma_n(x-\ell)}) \Psi_n(y-vt), & x \in (0, \ell), \\ \sum_{n \in \mathbb{Z}} t_n e^{-i\omega_n t} e^{i\alpha_n(x-\ell)} \Phi_n(y), & x \in (\ell, +\infty), \end{cases} \quad (27) $$

where $\hat{\omega}_0 = (\omega_0 - v\beta_0)$ and with $(\Phi_n, \alpha_n)$ in (12)-(13), $(\Psi_n, \gamma_n)$ in (16)-(17). It is worth noticing that (27) is consistent as we have

$$ e^{-i\omega_n t} \Phi_n(y) = e^{-i\hat{\omega}_0 t} \Phi_n(y - vt), \quad (28) $$

hence the matchings at $x=0$ and $x=\ell$ can be done once for all, being valid for any time $t$ in $y = (y - vt) \in (0, \xi h)$. The continuity of $\varphi$ on $y \in \mathcal{Y}_\xi = (0, \xi h)$ applies at $x=0$ and $x=\ell$, resulting in

$$ \begin{aligned} \int_{\mathcal{Y}_\xi} \varphi(0^{-}, y, t) \Psi_m^*(y) dy &= \int_{\mathcal{Y}_\xi} \varphi(0^{+}, y, t) \Psi_m^*(y) dy, \\ \int_{\mathcal{Y}_\xi} \varphi(\ell^{-}, y, t) \Psi_m^*(y) dy &= \int_{\mathcal{Y}_\xi} \varphi(\ell^{+}, y, t) \Psi_m^*(y) dy. \end{aligned} \quad (29) $$
Next, with $Y = (0, h)$, the continuity $\partial_x\varphi$ on $Y_\xi$ is used along with the Neumann condition on $Y\setminus Y_\xi$, namely $\partial_x\varphi(0^-, y, t) = \partial_x\varphi(\ell^+, y, t) = 0$ for $y \in (Y\setminus Y_\xi)$. Using (28), thus we have

$$ 
\begin{align} 
\int_Y \partial_x \varphi(0^-, y, t) \Phi_m^*(y) dy &= \int_{Y_\xi} \partial_x \varphi(0^+, y, t) \Phi_m^*(y) dy, \tag{30} \\ 
\int_{Y_\xi} \partial_x \varphi(\ell^-, y, t) \Phi_m^*(y) dy &= \int_Y \partial_x \varphi(\ell^+, y, t) \Phi_m^*(y) dy. \nonumber 
\end{align} 
$$

The system (29)-(30) is classical owing to the orthogonality relations (14)-(18). It involves the coupling matrices

$$ b_{mn} = \int_{Y_\xi} \Psi_m^*(y) \Phi_n(y) dy = \frac{\sqrt{(2-\delta_{m0})\xi}}{2} \left( e^{iq^+} \mathrm{sinc}(q^+) + e^{iq^-} \mathrm{sinc}(q^-) \right) \quad (31) $$

where $q^\pm = \frac{(\beta_0 - Mk_0)\xi h}{2(1-M^2)} + n\pi\xi \pm \frac{m\pi}{2}$. With the above definition of $b$, the matchings (29) provide for $m \in \mathbb{N}$

$$ b_{m0} + \sum_{n \in \mathbb{Z}} b_{mn} r_n = a_m^+ + a_m e^{i\gamma_m \ell}, \quad a_m^+ e^{i\gamma_m \ell} + a_m^- = \sum_{n \in \mathbb{Z}} b_{mn} t_n, \quad (32) $$

and the matchings (30) provide, for $m \in \mathbb{Z}$

$$ 
\begin{align} 
\alpha_0 \delta_{m0} - \alpha_m r_m &= \sum_{n \in \mathbb{N}} b_{nm}^* \gamma_n (a_n^+ - a_n^- e^{i\gamma_m \ell}), \tag{33} \\ 
\sum_{n \in \mathbb{N}} b_{nm}^* \gamma_n (a_n^+ e^{i\gamma_m \ell} - a_n^-) &= \alpha_m t_m. \nonumber 
\end{align} 
$$

The numerical implementation of the multimodal analysis is then also classical using a truncated version of (29)-(30). In doing so, we obtain the unknown scattering coefficients $(r_n, t_n)$, $n = -M, \dots, M$, outside the grating and the amplitudes $a_n^\pm$, $n = 0, \dots, N$, of the right-going and left-going waves within the grating, hence the whole solution (27). Specifically, the final system reads $M\mathbf{v} = \mathbf{s}$, with

$$ M = \begin{pmatrix} -b & l & d(e^{i\gamma\ell}) & 0 \\ 0 & d(e^{i\gamma\ell}) & l & -b \\ d(\alpha) & cd(\gamma) & -cd(\gamma e^{i\gamma\ell}) & 0 \\ 0 & -cd(\gamma e^{i\gamma\ell}) & cd(\gamma) & d(\alpha) \end{pmatrix} \quad (34) $$

where we have defined $d(e)$ the diagonal matrix with diagonal elements $e_n$, $c = (b^*)^\mathrm{T}$ the transpose complex conjugate of $b$, and

$$ \mathbf{v} = \begin{pmatrix} \mathbf{r} \\ \mathbf{a}^+ \\ \mathbf{a}^- \\ \mathbf{t} \end{pmatrix}, \quad \mathbf{s} = \begin{pmatrix} \mathbf{s}^1 \\ \mathbf{0} \\ \mathbf{s}^2 \\ \mathbf{0} \end{pmatrix}, \quad (35) $$
with $s_m^1 = b_{m0}$ and $s_m^2 = \alpha_0 \delta_{m0}$, which can be solved by simple matrix inversion. For completeness, we provide in the appendix Appendix B the corresponding Matlab script. It is worth noticing that the resolution holds for any time, and the representation of the solution at different times is a simple post-treatment.

175

### 3.1. Scattering of an incident plane wave

To begin with, we consider a plane wave at frequency $\omega_0 h/c = 1$ and at oblique incidence $\theta = -40^{\circ}$ (hence $(\alpha_0, \beta_0)h = (0.77, -0.64)$) on a space-time modulated grating with $\xi = 0.8$ and $\ell = 10h$. At this relatively low frequency, for $M = 0$, there is a single propagating mode (or diffraction order) outside and inside the grating and this holds true up to $M = 0.74$. We report in figure 4 the real part of the potential field for increasing modulation speeds producing the appearance of higher space-time diffraction orders: for $M = 0.8$, the mode $n = 1$ with $(\alpha_1, \beta_1)h = (2.12, 5.64)$ and $\omega_1 h/c = 6.03$; for $M = 0.9$, the mode $n = 1$ with $(\alpha_1, \beta_1)h = (3.53, 5.64)$ and $\omega_1 h/c = 6.65$, and the mode $n = 2$ with

180

185

<figure><img src="image_5.png" /><figcaption>Figure 4: Real part of the potential fields for a plane wave at incidence $\theta = -40^{\circ}$ on the space-time modulated grating ($\xi = 0.8$ and $\ell/h = 10$). The left panels show the diffraction patterns for $x/\ell \in (-1, 2)$. The right panels show magnifications for increasing times. Higher space-time diffraction orders appear for $M = 0.8$ ($n = 1$) and for $M = 0.9$ ($n = 1, 2$).</figcaption></figure>
($\alpha_2, \beta_2$)h = (3.06, 11.92) and $\omega_2h/c = 12.31$. It is noticeable that, although the sound-rigid parts of the grating are thin, the diffraction strength is significant; in particular strong evanescent fields are triggered at the junctions of the grating, particularly visible on the temporal sequences for $t \in (0, 0.6T)$ with $T = h/v$ the period of the modulation.

## 3.2. Case of a Gaussian beam

Following [24], we now illustrate the non reciprocal and asymmetric sound diffraction by considering an incident beam. This is done by computing the elementary contributions $\varphi(x, y, t; \theta) = \varphi(x, y, t)$ in (27) (we make explicit the dependence in the incidence $\theta$) for a range of incidence $\theta \in (\theta_0 - \Delta\theta, \theta_0 + \Delta\theta)$ and by weighting these contributions according to

$$ \varphi_{\text{beam}}(x, y, t) = \int_{\theta_0 - \Delta\theta}^{\theta_0 + \Delta\theta} f(\theta)\varphi(x, y, t; \theta)d\theta. \quad (36) $$

Various weighting functions $f(\theta)$ are possible and we choose

$$ f(\theta) = e^{-i(k \cos \theta x_0 + k \sin \theta y_0)} e^{-\sigma^2 (\cos \theta - \cos \theta_0)^2}, $$

with $\sigma = 10$ the width of the Gaussian beam and $(x_0, y_0) = (-50, 80)$ which fixes the position of the beam and $\Delta\theta = \theta_0/2$. We notice that the spatial separation of the different diffraction orders is attributable to two ingredients. On the one hand, from (19) (see also figure 2), the mode $n$ propagates with a direction $\theta_n$ given by

$$ \sin \theta_n = \frac{\sin \theta_0 + \frac{2n\pi}{k_0h}}{1 + M \frac{2n\pi}{k_0h}}, \quad n \in \mathbb{Z}, \quad (37) $$

with the limiting case $M=0$ being classical. It is also consistent with Eq. (9) in [24] with $K = 2\pi/h$ and $\Omega = M\frac{2\pi}{k_0}\omega_0$ the spatial and temporal periodicities of the modulated grating. However, note that, in our case, a spatial diffraction order is associated with a single temporal diffraction order. The striking feature is that for $Mn > 0$ new interference orders can appear in the subwavelength regime $k_0h \ll 1$. Next, within the grating, the modulation at relative speed $M$ produces a shift of the direction of propagation which, from (24) (see also figure 3), reads

$$ \sin \Theta_n = \frac{M}{\sqrt{1 - \frac{(1-M^2)^3}{(1-M \sin \theta_0)^2} \left( \frac{n\pi}{\xi k_0 h} \right)^2}}, \quad n \in \mathbb{N}. \quad (38) $$

In particular, $\sin \Theta_0 = M$, which is conform with a propagation at speed $c$ with a vertical component $v$ imposes by the grating modulation.

We report in figures 5 and 6 the diffraction patterns for a beam with centre incidence angle $\theta_0 = \pm 40^\circ$ and different values of the modulation speeds $M \ge 0$. We still consider $\omega_0h/c = 1$, $\xi = 0.8$ and to make the spatial separation more visible, we use $\ell = 30h$.
<figure><img src="image_6.png" /><figcaption>Figure 5: Scattering of an incident beam at incidence $\theta = 40^\circ$ and $k_0h = 1$ in the reference case $M = 0$, for $M = 0.65$ producing an almost straight propagation of the beam ($\Theta_0 \simeq \theta_0$) and $M = 0.9$. The directions $\theta_0$ and $\Theta_0 = \sin^{-1}(M)$ are shown with plain lines (black and red, respectively). The grating has length $l = 30h$ with $\xi = 0.8$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Scattering of an incident beam at main incidence $\theta = -40^\circ$ and $k_0h = 1$. Same representation as in figure 5. Outside and inside the grating, increasing $M$ produces the appearance of higher interference orders: $n = 0$ for $M = 0.7$, $n = 0, 1$ for $M = 0.8$ and $n = 0, 1, 2$ for $M = 0.9$.</figcaption></figure>
For $\theta_0 = 40^\circ$, there is a single diffraction order $n = 0$ and a single propagating mode within the grating. Accordingly, the beam is partly reflected and transmitted in the directions $\pm\theta_0$ and within the grating it propagates along $\Theta_0 = \sin^{-1}(M)$ (independent of $\theta_0$). In particular for $M = 0.65$, $\Theta_0 \simeq \theta_0$ resulting in a straight propagation of the beam.

For $\theta_0 = -40^\circ$, the most striking feature in the diffraction patterns is the apparent negative refraction within the grating (e.g. for $M = 0.7$), a feature that has been used for unmodulated grating to produce perfectly transmitting metamaterial [27]. Next, several diffraction orders take place when increasing $M > 0$ whose characteristics are collected in table 1.

<table><caption>Table 1: Space-time diffraction orders outside the grating from (20) and (37), and modes inside the grating from (38) for the modulation speeds considered in figure 6 ($\theta_0 = -40^\circ$, $\omega_0h/c = 1$).</caption><thead><tr><th></th><th>θ<sub>1</sub></th><th>ω<sub>1</sub>h/c</th><th>θ<sub>2</sub></th><th>ω<sub>2</sub>h/c</th><th>Θ<sub>0</sub></th><th>Θ<sub>1</sub></th><th>Θ<sub>2</sub></th></tr></thead><tbody><tr><td>M = 0.7</td><td>—</td><td>—</td><td>—</td><td>—</td><td>44.4°</td><td>—</td><td>—</td></tr><tr><td>M = 0.8</td><td>69.4°</td><td>6.03</td><td>—</td><td>—</td><td>53.1°</td><td>75.0°</td><td>—</td></tr><tr><td>M = 0.9</td><td>57.9°</td><td>6.65</td><td>75.6°</td><td>12.31</td><td>64.1°</td><td>66.9°</td><td>81.0°</td></tr></tbody></table>

As already commented in [24], the diffraction is asymmetrical and nonreciprocal. Both properties are due to the breaking of the invariance $(\theta_0, n) \rightarrow (-\theta_0, -n)$ which is replaced by the invariance $(M, \theta_0, n) \rightarrow (-M, -\theta_0, -n)$, see (19)-(20). In the presence case, this is equivalent to the effect of flows on sound propagation, a property which has been used to break the time invariance [28, 29] and non-reciprocal scattering [30].

## 4. Time-domain analysis

We now move to the time-domain and focus on the temporal separation of the harmonics, in addition to their spatial separation. We consider an incident pulse at normal incidence described by the real signal $s(t)$ coming from $x = -\infty$ and we denote $\hat{s}(\omega)$ its Fourier transform ($\hat{s}(\omega) = \frac{1}{2\pi} \int_{-\infty}^{\infty} s(t)e^{i\omega t}dt$ and $s(t) = 2\Re \int_{0}^{\infty} \hat{s}(\omega)e^{-i\omega t}d\omega$ where $\Re(z)$ denotes the real part of $z$). With an incident wave in the form

$$ \varphi^{\text{inc}}(x, y, t) = 2\Re \int_{0}^{\infty} e^{-i\omega t} e^{i\alpha x} \hat{s}(\omega) d\omega = s(t - x/c), \quad (39) $$
with $\alpha = \omega/c$, and using (27), the potential field $\varphi = \varphi(x, 0, t)$ at $y = 0$ reads

$$ \varphi = \begin{cases} \varphi^{\text{inc}}(x, 0, t) + 2\Re \sum_{n \in \mathbb{Z}} \int_0^\infty r_n(\omega) e^{-i\omega_n t} e^{-i\alpha_n x} \hat{s}(\omega) d\omega, & (40) \\ 2\Re \sum_{n \in \mathbb{N}} \int_0^\infty e^{-i\omega t} (a_n^+(\omega) e^{i\gamma_n x} + a_n^-(\omega) e^{-i\gamma_n(x-\ell)}) \Psi_n(-vt) \hat{s}(\omega) d\omega, & \\ 2\Re \sum_{n \in \mathbb{Z}} \int_0^\infty t_n(\omega) e^{-i\omega_n t} e^{i\alpha_n(x-\ell)} \hat{s}(\omega) d\omega, & \end{cases} $$

for $x \in (-\infty, 0)$, $x \in (0, \ell)$, $x \in (\ell, +\infty)$, respectively, and with $\omega_n = \omega + M2n\pi$ and $\alpha_n = \sqrt{\omega_n^2 - (2n\pi)^2}$, from (20) (for simplicity, we work in dimensionless form by setting $h = c = 1$). Within the grating, $\Psi_n(-vt)$ is given by (16) when $-vt \in (p, p + \xi)$ with $p$ an integer, otherwise $\Psi_n(-vt) = 0$ (the line $y = 0$ falls in a region of air or in a rigid part of the grating).

We consider a narrow band Gaussian pulse with center frequency $\omega_c$ and bandwidth $\tau$

$$ s(t) = \cos(\omega_c t) e^{-(t/\tau)^2}, \quad \hat{s}(\omega) = \frac{\tau}{4\sqrt{\pi}} \left( e^{-(\omega-\omega_c)^2 \tau^2 / 4} + e^{-(\omega+\omega_c)^2 \tau^2 / 4} \right), \quad (41) $$

and we choose $\tau = 30$, $\ell = 100$, $\xi = 0.8$ and $\omega_c/2\pi = 0.1, 0.2$ and $0.3$ and $M = 0.8$. In all cases, the potential field (40) is calculated using $10^4$ values of $\omega/2\pi \in (0, 1)$. Then we use two representations:

(1) Snapshots for which $\varphi(x, 0, t_j)$ is represented against $x \in (-200, 300)$ ($2.10^3$ values of $x$ are used) at discrete times $t_j \in (-100, 500)$, $j = 1, \dots, 24$.

(2) Time records $\varphi(-d, 0, t)$ at discrete values $x = -d$ for $d = 0, 20, 40$ against time $t \in (-100, 600)$ using $10^4$ values of times.

We report in figure 7 the snapshots and time records for $\omega_c/2\pi = 0.1$. At this centre frequency, a single mode $n = 0$ is propagating inside and outside the grating. Accordingly, we observe that the incident pulse is multiply reflected and transmitted as it would be for $M = 0$ (in panel (a) we observe two low reflected pulses and a single high transmitted pulse). Next, within the grating, the pulse propagates with an angle $\Theta_0 = \sin^{-1}(M)$. As a consequence we observe that the reflected pulses emerge from the grating with a time separation $2\ell_0/c$, $\ell_0 = \ell/\cos\Theta_0 = \ell/\sqrt{1-M^2}$ in panel (b), conform with the apparent lower speed of sound $c\cos\Theta_0$ in the panel (a). In this simple case, the signature of the temporal harmonics is the contribution of evanescent modes triggered at the grating interface, particularly visible on the first pulse of $\varphi(0, 0, t)$. As these harmonics have a time-dependence $e^{-i\omega_n t}$ from (20), the spectral content of the signal near the interface is much wider than that of the incident pulse.

In figure 8, the same representation is shown for a pulse with a center frequency $\omega_c/2\pi = 0.2$. As previously, a single mode is propagating outside and
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Case $\omega_c/2\pi = 0.1$ – (a) Snapshots of the diffraction of the incident pulse by a modulated grating at time $t_i \in (-100, 500)$. (b) Temporal records at $x = -d$, with $d = 0, 20$ and $40$. The two reflected pulses have a time separation $\ell_0 = \ell/\sqrt{1 - M^2}$.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: Case $\omega_c/2\pi = 0.2$ – Same representation as in figure 7; the dashed lines show the edge of the cone $\theta_1 \in (\theta_{1m}, 90^\circ)$ near the grazing angle.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 9: Case $\omega_c/2\pi = 0.3$ – Same representation as in figure 7.</figcaption></figure>
inside the grating; accordingly, in the reported time range, we observe two reflected pulses and one transmitted pulse (panel (a)). Next, the mode $n = 1$ is at its cut-off frequency with $\sin\theta_1 = 1$ (from (37) with $\omega_c/2\pi = 1 - M = 0.2$). This mode emerges at the grazing angle and it generates continuously the mode 0 by mode conversion; this is particularly visible on the time record $\varphi(-40,0,t)$ where we observe a dissymmetry of the first pulse for $t > 100$ (panel (b)). Eventually, due to the narrow but finite spectral content of the incident pulse, we observe a long wave packet with higher frequencies. A Fourier transform of this wave packet shows that it has been generated by mode conversion of the fundamental frequency range $\omega \in (\omega_c, \omega_m)$, $\omega_m/2\pi \approx 0.22$, corresponding to $\theta_1 \in (\theta_{1m}, 90^\circ)$ with $\theta_{1m} \approx 80^\circ$. This is conform with the observations in panel (a) where, for $t > 0$, part of the reflected signal remains close the interface, the same in transmission for $t > 120$ (see additional results in appendix Appendix C).

Eventually for $\omega_c/2\pi = 0.3$, reported in figure 9, two modes $n=0$ and $n=1$ are propagating inside and outside the grating. The two modes are now associated with well defined reflected and transmitted pulses (without distorsion), and the pulses, with different frequencies, are separated in space and time far enough from the interface. This is visible for instance on $\varphi(-40,0,t)$ for the first reflected pulses at $\omega_c$ and $\omega_{c1}/2\pi = \omega_c + M = 1.1$ around $t = 50$ and $t = 100$ (panel (b)). Besides, the mode conversions at each interfaces, mode 0 $\rightarrow$ mode 1 and mode 1 $\rightarrow$ mode 0, now produce additional pulses in reflection and transmission. For $t > 250$, 3 reflected pulses associated with the mode 0 and 3 reflected pulses associated with the mode 1 are superimposed (and they can be separated by Fourier transform, revealing each pulse without distorsion, see appendix Appendix C).

## 5. Concluding remarks

We have analyzed the dispersion of a space-time modulated gratings alternating sound-rigid and sound-soft layers. This has been done by deriving the exact dispersion relations of the modes inside and outside the grating. The use of these modes is the basis of the numerical (multimodal) method proposed to solve easily classical diffraction problems. Examples for an incident plane wave and for an incident beam have been given in the harmonic regime. Expectedly, the grating allows for the generation of space and time diffraction orders and the diffraction is asymmetrical and nonreciprocal. These properties are similar, although non identical, to the properties of sound scattering in the presence of flow. In particular, each spatial diffraction order is associated to a single temporal harmonic. We also provided the analysis of the diffraction of an incident pulse in the transient regime revealing a complex mixing of the space-time diffraction orders, notably at the threshold of appearance of a new diffraction order. The present analysis applies straightforwardly to the case of a space-time modulated grating alternating sound-soft layers with contrasts in the mass density and/or in the bulk modulus.
## Appendix A. Boundary condition at the moving rigid boundaries

The boundary condition at the moving rigid boundaries is deduced from the equation of mass conservation over a volume $\Omega$ in which a source of mass $\dot{m}$ located along $\Gamma$ is used to produce the space-time modulation of the mass density (figure A.10). Specifically, we have

$$ \frac{dM}{dt} = \frac{d}{dt} \int_{\Omega} \rho_t dr = \dot{m}. \qquad (A.1) $$

Because of the source term, the boundary $\Gamma \subset \partial\Omega$ moves at velocity $\mathbf{v}$ along its normal $\mathbf{n}$ and we denote $\dot{m} = \int_{\Gamma} j dx$. The application of the Reynolds transport theorem to the above equation provides

$$ \int_{\Omega} \left( \frac{\partial \rho_t}{\partial t} + \operatorname{div}(\rho_t \mathbf{u}_t) \right) d\mathbf{r} + \int_{\Gamma} \rho_t (\mathbf{v} - \mathbf{u}_t) \cdot \mathbf{n} \, dx = \int_{\Gamma} j \, dx, \qquad (A.2) $$

and at this stage, $(\rho_t, \mathbf{u}_t)$ are the total mass density and velocity of the air. From the above, we deduce the local equations

$$ \frac{\partial \rho_t}{\partial t} + \operatorname{div}(\rho_t \mathbf{u}_t) = 0, \quad \rho_t (\mathbf{v} - \mathbf{u}_t) \cdot \mathbf{n} = j. \qquad (A.3) $$

<figure><img src="image_14.png" /><figcaption>Figure A.10: Discontinuity $\Sigma$ moving at velocity $v$ in the $(x, y, t)$ space.</figcaption></figure>

Now, we linearize the problem for small acoustic perturbations (the equilibrium is the air at rest with mass density $\rho$); introducing $\varepsilon \ll 1$, we have $\rho_t = \rho + \varepsilon\rho'$, $\mathbf{u}_t = \varepsilon\mathbf{u}$. The left hand-side equation in (A.3) provides the usual (local) equation of mass conservation for the acoustic perturbation after introduction of the acoustic pressure $p = c^2\rho'$ which defines the acoustic velocity $c$. The right hand-side equation in (A.3) provides the boundary condition we are looking for. Namely we obtain after linearization

$$ \left( \rho + \varepsilon \frac{p}{c^2} \right) (\mathbf{v} - \varepsilon \mathbf{u}) \cdot \mathbf{n} = j, \quad \rightarrow \quad j = \rho \mathbf{v} \cdot \mathbf{n}, \quad \rho \mathbf{u} \cdot \mathbf{n} - v \frac{p}{c^2} = 0. \qquad (A.4) $$

The relation $j = \rho \mathbf{v} \cdot \mathbf{n}$ is expected as the mass is either created ($\mathbf{n} = \mathbf{e}_y$) or removed ($\mathbf{n} = -\mathbf{e}_y$) at the moving boundary; the resulting boundary condition links the acoustic pressure and velocity through the relation announced in (4), with $\mathbf{u} \cdot \mathbf{n} = u_y$, $Z = \rho c$ and $M = v/c$. We notice that such boundary condition is used in [25] from compatibility relations; it is also a relation similar to Rankine-Hugoniot condition between two fluids in the presence of a line of discontinuity.
Appendix B. Complement on the numerical implementation

We provide below a MATLAB M-file script used to compute the solutions presented in the present study. The first script computes the solution $\varphi(x, y, t)$ in (27) as described in §3 (specifically, the system (34)-(35)). The same procedure has been used to implement the problem of the scattering of an incident beam and of incident pulse thanks to the linearity of the solutions for the incident waves in (36) and (39). For the ease of use, we also provide an additional script that can be used to represent the corresponding two-dimensional field as shown in figures 4.

% entries

c=1;h=1;om0=1;xi=.8;V=-.5;L=10;theta=40*pi/180;

% constant

m=V/c; k0=om0/c; bet0=k0*sin theta;

%% pseudo periodic modes

M1=100; N2=100; n1=-M1:M1; N1 =length(n1); bet=bet0+2*n1*pi/h; alp=sqrt((k0+m*2*n1*pi/h).^2-bet.^2);

%% modes in the slab

n2=0:N2-1; no=sqrt((2-(n2==0))/(xi*h));
gam=sqrt( (k0-m*bet0)^2/(1-m^2)-(1-m^2)*(n2*pi/(xi*h)).^2);

% the system

matn=ones(N2,1)*n1; matm=n2' ones(1,N1);
temp=(bet0-m*k0)*(xi*h)/(2*(1-m^2))+matn*pi*xi;
qm= temp-matm*pi/2; qp= temp-matm*pi/2;
temp=1/2*sqrt((2-(matm==0))*xi);
mb= temp.* (exp(1i*qp).*sinc(qp/pi)+exp(1i*qm).*sinc(qm/pi));
mc=mb'; z11= zeros(N1,N1); z21= zeros(N2,N1);

mM=[-mb, eye(N2), diag(exp(1i*gam*L)), z21;
diag(alp), mc*diag(gam), -mc*diag(gam.*exp(1i*gam*L)), z11;
z21, diag(exp(1i*gam*L)), eye(N2), -mb;
z11,mc*diag(gam.*exp(1i*gam*L)), -mc*diag(gam), -diag(alp)];
s1=mb(:,M1+1); s2= zeros(N1,1); s2(M1+1)= alp(M1+1);
s=[s1;s2; zeros(N1+N2,1)];
V=mM's;

% the coefficients entering in the solution

R=V(1:N1); p=[R]; n =length(p);
ap=V(n+1:n+N2); p=[p;ap]; n =length(p);
am=V(n+1:n+N2); p=[p;am]; n =length(p);
T=V(n+1:end);
```matlab
% entries
Lplot=10; Nplot=20; Nx=100*Lplot; Ny=20;

% field in (-Lplot,0)
xm=linspace(-Lplot,0,Nx)*ones(1,Ny); ym=ones(Nx,1)*linspace(0,h,Ny);
Phi=exp(1i*bet0*ym); phim=exp(1i*alp(M1+1)*xm).*Phi;
for ii=1:N1,
    Phi=exp(1i*bet(ii)*ym); phim=phim+R(ii)*exp(-1i*alp(ii)*xm).*Phi;
end

% field in the slab (0,L)
xs=linspace(0,L,Nx)*ones(1,Ny); ys=ones(Nx,1)*linspace(0,xi*h,Ny);
phis=0;
for j=1:N2,
    Psi=no(j)*exp(1i*m*(k0-m*bet0)/(1-m^2)*ys).*cos(n2(ii)*pi/(xi*h)*ys);
    phis=phis+(Ap(j)*exp(1i*gam(j)*xs)+Am(j)*exp(-1i*gam(j)*(xs-L))).*Psi;
end

% transmitted field in (L,L+Lplot)
xp=linspace(L,L+Lplot,Nx)*ones(1,Ny); yp=ones(Nx,1)*linspace(0,h,Ny);
phip=0;
for j=1:N1,
    Phi=exp(1i*bet(j)*yp);
    phip=phip+T(j)*exp(1i*alp(j)*(xp-L)).*Phi;
end

% graphical representation
for j=1:Nplot,
    pcolor(xm,ym+(j-1)*h,real(phim.*exp(1i*(j-1)*bet0*h))), hold on
    pcolor(xs,ys+(j-1)*h,real(phis.*exp(1i*(j-1)*bet0*h))),
    pcolor(xp,yp+(j-1),real(phip.*exp(1i*(j-1)*bet0*h))),
end
shading flat, colormap jet, axis auto, axis equal
```

## Appendix C. Additional results

345 In this section, we provide additional results which complement the analysis of the figures 8 and 9. This is done by windowing the Fourier transform $\hat{\varphi}(-d, 0, \omega)$, $\omega \in (0, \omega_f)$ of $\varphi(-d, 0, t)$, and then by coming back in the time domain, namely we use

$$ \varphi = \varphi_0 + \varphi_1, \qquad (C.1) $$

with

$$ \left\{ \begin{array}{ll} \varphi_0(-d, 0, t) = \int_{W_0} \hat{\varphi}(-d, 0, \omega) e^{-i\omega t} \, d\omega, & W_0 = (0, \omega_c + \Delta\omega), \\ \\ \varphi_1(-d, 0, t) = \int_{W_1} \hat{\varphi}(-d, 0, \omega) e^{-i\omega t} \, d\omega, & W_1 = (\omega_{c1} - \Delta\omega, \omega_f), \end{array} \right. \qquad (C.2) $$

350 with $\omega_{c1} = \omega_c + 2\pi M$ the centre frequency of the mode 1, from (19). The two Fourier components are well separated, hence this basic filtering is sufficient, and we choose $\Delta\omega = \frac{1}{2}(\omega_{c1} - \omega_c)$.
The result is shown in figure C.11 for $\omega_c/2\pi = 0.2$. As previously said, each reflected pulses in $\varphi_0$ are asymmetrical as a long coda is generated by mode conversion of the grazing mode 1. The same codas are observed on $\varphi_1$ due to the grazing mode 1 without mode conversion. Next, increasing $d$ produces arrival times of highly distorted pulses alimented by the cone of propagating modes near the grazing angle. It is worth noticing that the amplitude of this group of harmonics is higher than that of the fundamental mode which is conform with the classical Wood's anomaly (as a diffracted order appears at grazing angle, the power in that order is removed from the zero order signal).

<figure><img src="image_15.png" /><figcaption>Figure C.11: Case $\omega_c/2\pi = 0.2$ – Snapshots of $\varphi_0$ (top panels, in blue) and $\varphi_1$ (bottom panels, in green) defined in (C.2) and complementing the figure 8.</figcaption></figure>

The result for $\omega_c/2\pi = 0.3$ is reported in figure C.12. We already notice an increase in the number of reflected pulses for the two propagating modes 0 and 1. The first pulses are due to direct reflection of the incoming mode 0 without $(0 \to 0)$ and with $(0 \to 1)$ mode conversion. The following 3 pulses emerge in
reflection after a round trip in the grating involving a first transmission of the incoming mode 0 at $x = 0$, a reflection at $x = \ell$ and a second transmission at $x = 0$ toward $x < 0$. For $\varphi^0$, the 3 pulses, in time order, are the result of: for the 2th pulse, zero conversion ($0 \to 0$ at each interface), for 3th pulse, two conversions $0 \to 1$ at $x = 0$ and $1 \to 0$ at $x = \ell$ (or the opposite), and for the 4th pulse, two conversions $0 \to 1$ and $1 \to 0$ at $x = 0$. We observe an apparent lower velocity $c \cos \theta_1$ of the mode 1 due to oblique propagation outside the grating, forth and back.

370

<figure><img src="image_16.png" /><figcaption>Figure C.12: Same representation as in figure C.11 for $\omega_c/2\pi = 0.3$ complementing figure 9.</figcaption></figure>

References

375

[1] J.-C. Simon, Action of a progressive disturbance on a guided electromagnetic wave, IRE Transactions on Microwave Theory and Techniques 8 (1)
(1960) 18–29.

[2] A. Oliner, A. Hessel, Wave propagation in a medium with a progressive sinusoidal disturbance, IRE Transactions on Microwave Theory and Techniques 9 (4) (1961) 337–343.

380

[3] E. S. Cassedy, A. A. Oliner, Dispersion relations in time-space periodic media: Part i—stable interactions, Proceedings of the IEEE 51 (10) (1963) 1342–1359.

[4] H. Nassar, X. Xu, A. Norris, G. Huang, Modulated phononic crystals: Non-reciprocal wave propagation and willis materials, Journal of the Mechanics and Physics of Solids 101 (2017) 10–29.

385

[5] L. Zhang, X. Q. Chen, R. W. Shao, J. Y. Dai, Q. Cheng, G. Castaldi, V. Galdi, T. J. Cui, Breaking reciprocity with space-time-coding digital metasurfaces, Advanced materials 31 (41) (2019) 1904069.

[6] J. Li, C. Shen, X. Zhu, Y. Xie, S. A. Cummer, Nonreciprocal sound propagation in space-time modulated media, Physical Review B 99 (14) (2019) 144311.

390

[7] P. A. Huidobro, E. Galiffi, S. Guenneau, R. V. Craster, J. B. Pendry, Fresnel drag in space-time-modulated metamaterials, Proceedings of the National Academy of Sciences 116 (50) (2019) 24943–24948.

395

[8] J. Marconi, E. Riva, M. Di Ronco, G. Cazzulani, F. Braghin, M. Ruzzene, Experimental observation of nonreciprocal band gaps in a space-time-modulated beam using a shunted piezoelectric array, Physical Review Applied 13 (3) (2020) 031001.

[9] A. E. Cardin, S. R. Silva, S. R. Vardeny, W. J. Padilla, A. Saxena, A. J. Taylor, W. J. Kort-Kamp, H.-T. Chen, D. A. Dalvit, A. K. Azad, Surface-wave-assisted nonreciprocity in spatio-temporally modulated metasurfaces, Nature communications 11 (1) (2020) 1–9.

400

[10] Z. Wu, A. Grbic, Serrodyne frequency translation using time-modulated metasurfaces, IEEE Transactions on Antennas and Propagation 68 (3) (2019) 1599–1606.

405

[11] S. Taravati, G. V. Eleftheriades, Full-duplex nonreciprocal beam steering by time-modulated phase-gradient metasurfaces, Physical Review Applied 14 (1) (2020) 014027.

[12] L. Zhang, M. Z. Chen, W. Tang, J. Y. Dai, L. Miao, X. Y. Zhou, S. Jin, Q. Cheng, T. J. Cui, A wireless communication scheme based on space- and frequency-division multiplexing using digital metasurfaces, Nature electronics 4 (3) (2021) 218–227.

410
415

[13] X. Wang, V. S. Asadchy, S. Fan, S. Tretyakov, Space-time metasurfaces for perfect power combining of waves, arXiv preprint arXiv:2105.14627 (2021).

420

[14] Y. Chen, X. Li, H. Nassar, A. N. Norris, C. Daraio, G. Huang, Nonreciprocal wave propagation in a continuum-based metamaterial with space-time modulated resonators, Physical Review Applied 11 (6) (2019) 064052.

[15] G. Castaldi, L. Zhang, M. Moccia, A. Y. Hathaway, W. X. Tang, T. J. Cui, V. Galdi, Joint multi-frequency beam shaping and steering via space-time-coding digital metasurfaces, Advanced Functional Materials 31 (6) (2021) 2007620.

425

[16] S. A. Stewart, T. J. Smy, S. Gupta, Finite-difference time-domain modeling of space-time-modulated metasurfaces, IEEE Transactions on Antennas and Propagation 66 (1) (2017) 281-292.

430

[17] S. A. Stewart, T. J. Smy, S. Gupta, Finite-difference time-domain modeling of space-time-modulated metasurfaces, IEEE Transactions on Antennas and Propagation 66 (1) (2017) 281-292.

[18] N. Chamanara, Y. Vahabzadeh, C. Caloz, Simultaneous control of the spatial and temporal spectra of light with space-time varying metasurfaces, IEEE Transactions on Antennas and Propagation 67 (4) (2019) 2430-2441.

435

[19] S. Inampudi, M. M. Salary, S. Jafar-Zanjani, H. Mosallaei, Rigorous space-time coupled-wave analysis for patterned surfaces with temporal permit-tivity modulation, Optical Materials Express 9 (1) (2019) 162-182.

440

[20] B. M. Goldsberry, S. P. Wallen, M. R. Haberman, Nonreciprocal vibrations of finite elastic structures with spatiotemporally modulated material properties, Physical Review B 102 (1) (2020) 014312.

[21] B. M. Goldsberry, S. P. Wallen, M. R. Haberman, Nonreciprocity and mode conversion in a spatiotemporally modulated elastic wave circulator, arXiv preprint arXiv:2110.01745 (2021).

445

[22] Y. Khorrami, D. Fathi, A. Khavasi, R. C. Rumpf, Dynamical control of multilayer spacetime structures using extended fourier modal method, IEEE Photonics Journal 13 (6) (2021) 1-10.

[23] C. Scarborough, Z. Wu, A. Grbic, Efficient computation of spatially discrete traveling-wave modulated structures, IEEE Transactions on Antennas and Propagation 69 (12) (2021) 8512-8525.

[24] S. Taravati, G. V. Eleftheriades, Generalized space-time-periodic diffraction gratings: Theory and applications, Physical Review Applied 12 (2) (2019) 024026.

450

[25] K. A. Lurie, Effective properties of smart elastic laminates and the screening phenomenon, International Journal of Solids and Structures 34 (13) (1997) 1633-1643.
455

[26] Z.-L. Deck-Léger, M. Skorobogatiy, C. Caloz, Diagrammatic explanation of the reverse doppler effect in space-time modulated photonic crystals, in: 2016 IEEE International Symposium on Antennas and Propagation (APSURSI), IEEE, 2016, pp. 2101–2102.

[27] R. Porter, Plate arrays as a perfectly-transmitting negative-refraction metamaterial, Wave Motion 100 (2021) 102673.

460

[28] P. Roux, M. Fink, Experimental evidence in acoustics of the violation of time-reversal invariance induced by vorticity, EPL (Europhysics Letters) 32 (1) (1995) 25.

[29] M. Fink, D. Cassereau, A. Derode, C. Prada, P. Roux, M. Tanter, J.-L. Thomas, F. Wu, Time-reversed acoustics, Reports on progress in Physics 63 (12) (2000) 1933.

465

[30] C. Saverna, Y. Aurégan, V. Pagneux, Non-reciprocal scattering in shear flow, The Journal of the Acoustical Society of America 146 (2) (2019) 1051–1060.