%~Mouliné par MaN_auto v.0.34.2 (f7f49b77) 2025-07-21 14:36:04
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Darche_20250278.bib}

\usepackage{cases}
\usepackage{siunitx}
\usepackage{subcaption}
%\usepackage{mathrsfs}

\usepackage[scr=boondox,cal=euler]{mathalfa}
%\newcommand{\mathscr}{\mathcal}

\newcommand*{\dd}{\mathrm{d}}

\newcommand{\bfGamma}{\boldsymbol{\Gamma}}
\newcommand{\bfPsi}{\boldsymbol{\Psi}}
\newcommand{\bfPhi}{\boldsymbol{\Phi}}
\newcommand{\bfA}{\mathbf{A}}
\newcommand{\bfB}{\mathbf{B}}
\newcommand{\bfC}{\mathbf{C}}
\newcommand{\bfD}{\mathbf{D}}
\newcommand{\bfE}{\mathbf{E}}
\newcommand{\bfF}{\mathbf{F}}
\newcommand{\bfI}{\mathbf{I}}
\newcommand{\bfU}{\mathbf{U}}
\newcommand{\bfV}{\mathbf{V}}
\newcommand{\bfW}{\mathbf{W}}
\newcommand{\bfzero}{\mathbf{0}}

\newcommand{\incident}{\mathrm{in}}
\newcommand{\reflected}{\mathrm{ref}}
\newcommand{\total}{\mathrm{tot}}
\newcommand{\transmitted}{\mathrm{trans}}

\newcommand{\Na}{\mathfrak{a}}
\newcommand{\Nb}{\mathfrak{b}}
\newcommand{\Nh}{\mathfrak{h}}
\newcommand{\Ng}{\mathfrak{g}}

\makeatletter
\renewcommand\p@subfigure{}
\makeatother

%%%---------------------------------------------------------------------------

%%% BAR & OVERLINE

%%% Better widebar macro
\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
  }
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
  }%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\jump}{\llbracket}{\rrbracket}
\newcommand{\llangle}{{\langle\!\!\langle}}
\newcommand{\rrangle}{{\rangle\!\!\rangle}}
\DeclarePairedDelimiter{\dmean}{\llangle}{\rrangle}
\DeclarePairedDelimiter{\dmeanbig}{{\bigl\langle\!\!\!\bigl\langle}}{{\bigr\rangle\!\!\!\bigr\rangle}}
\DeclarePairedDelimiter{\floor}{\lfloor}{\rfloor}

%%%---------------------------------------------------------------------------

%%% VERTICAL BAR

%%% ''Evaluated at'' macro
\newcommand\restr[2]{{				
 \left.\kern-\nulldelimiterspace 
 #1 								
 \vphantom{\big\vert} 				
 \right\vert_{#2} 	
 }}

%%% ''Such that'' macro
\newcommand*{\st}[1][]{\;#1\vert\;}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}

%%%---------------------------------------------------------------------------

\makeatletter
\def\@setafterauthor{%
  \vglue3mm%
\begingroup\hsize=12cm\advance\hsize\abstractmarginL\raggedright
\noindent
   \leftskip\abstractmarginL
  \normalfont\Small
  \@afterauthor\par
\endgroup
\vskip2pt plus 3pt minus 1pt
}
\makeatother

%%%---------------------------------------------------------------------------

\title{Scattering of transient waves by an interface with time-modulated jump conditions}
\alttitle{Propagation d'ondes à travers une interface avec des conditions de saut modulées en temps}

\author{\firstname{Micha\"el} \lastname{Darche}\CDRorcid{0000-0001-5731-4551}}
\address{Aix Marseille Université, CNRS, Centrale Marseille, LMA UMR 7031, Marseille, France}
\email{darche@lma.cnrs-mrs.fr}

\author{\firstname{Rapha\"el} \lastname{Assier}}
\address{Department of Mathematics, University of Manchester, Oxford Road, Manchester M13 9PL, UK}
\email{raphael.assier@manchester.ac.uk}

\author{\firstname{S\'ebastien} \lastname{Guenneau}}
\address{UMI 2004 Abraham de Moivre-CNRS, Imperial College London, London SW7 2AZ, UK}
\address{Department of Physics, The Blackett Laboratory, Imperial College London, London SW7 2AZ, UK}
\email{s.guenneau@imperial.ac.uk}

\author{\firstname{Bruno} \lastname{Lombard}\IsCorresp}
\address[1]{Aix Marseille Université, CNRS, Centrale Marseille, LMA UMR 7031, Marseille, France}
\email{lombard@lma.cnrs-mrs.fr}

\author{\firstname{Marie} \lastname{Touboul}}
\address[3]{UMI 2004 Abraham de Moivre-CNRS, Imperial College London, London SW7 2AZ, UK}
\address{Department of Mathematics, Imperial College London, Huxley Building, Queen's Gate, London~SW7~2AZ, UK}
\address{POEMS, ENSTA, CNRS, INRIA, Institut Polytechnique de Paris, 91120 Palaiseau, France}
\email{marie.touboul@ensta.fr}

\thanks{M.D.\ was funded as a post-doctoral researcher by the Institut M\'ecanique et Ing\'enierie (Marseille, France). S.G.\ was funded by UK Research and Innvovation (UKRI) under the UK government's Horizon Europe funding guarantee (grant number 10033143)}
\CDRGrant[UK Research and Innovation (UKRI)]{10033143}

\keywords{\kwd{Elastic waves} \kwd{imperfect jump conditions} \kwd{time-varying media} \kwd{non-reciprocity} \kwd{numerical methods for hyperbolic equations}}
\altkeywords{\kwd{Ondes élastiques} \kwd{conditions de sauts imparfaites} \kwd{milieux modulés en temps} \kwd{non-réciprocité} \kwd{méthodes numériques pour les équations hyperboliques}}

\dateposted{2025-07-28}
\begin{document}

\maketitle

\begin{abstract}
Time modulation of the physical parameters offers interesting new possibilities for wave control. Examples include amplification of waves, harmonic generation and non-reciprocity, without resorting to non-linear mechanisms. Most of the recent studies focus on the time-modulation of the bulk physical properties. However, as the temporal modulation of these properties is difficult to achieve experimentally, we will concentrate here on the special case of an interface with time-varying jump conditions, which is simpler to implement. This work is focused on wave propagation in a one-dimensional medium containing one modulated interface. Properties of the scattered waves are investigated theoretically: energy balance, generation of harmonics, impedance matching and non-reciprocity. A fourth-order numerical method is also developed to simulate transient scattering. Numerical experiments are conducted to validate the numerical scheme and to illustrate the theoretical findings.
\end{abstract}

\begin{altabstract}
La modulation temporelle des propriétés d'un milieu de propagation offre de nouvelles possibilités intéressantes pour le contrôle des ondes. Les exemples incluent l'amplification des ondes, la génération d'harmoniques et la non-réciprocité, sans avoir recours à des mécanismes non linéaires. La plupart des études récentes se concentrent sur la modulation temporelle des propriétés physiques volumiques d'un milieu de propagation. Cependant, cette modulation étant difficile à réaliser expérimentalement, on se concentre ici sur le cas particulier d'une interface avec des conditions de saut variables en temps, plus simple à mettre en {\oe}uvre. Ce travail se concentre sur la propagation des ondes dans un milieu monodimensionnel contenant une interface modulée. On étudie théoriquement de nombreuses propriétés des ondes diffusées: bilan d'énergie, génération d'harmoniques, adaptation d'impédance et non-réciprocité. Une méthode numérique d'ordre~4 est développée pour simuler les ondes transitoires diffractées. Des expériences numériques sont menées pour valider le schéma numérique et illustrer les résultats théoriques attendus.
\end{altabstract}

\section{Introduction}\label{SecIntro}

Since the early 2000s, periodic modulation of physical properties in space has led to significant developments in wave control. Based on homogenization theory and Floquet--Bloch analysis, these so-called \emph{metamaterials} have made it possible to achieve cloaking, negative refraction or perfect lensing, to cite a few of the exotic effects. We refer the interested reader to~\cite{craster2013} for an overview of this abundant subject in the case of acoustics.

In recent years, the technical possibility of modulating physical properties in time has opened new perspectives in the field of metamaterials~\cite{caloz2019spacetime1,caloz2019spacetime2,galiffiA2022}. Time-varying metamaterials exhibit unusual phenomena such as time reflection and time refraction~\cite{wang_temporal_2025}, nontrivial topology~\cite{Swinteck2015,Lustig2018}, frequency conversion~\cite{Salehi2022,yi_frequency_2017} unidirectional and parametric amplification~\cite{torrentPRB2018,kimPRE2023,kiorpelidisPRB2024}. Numerous works have emerged since the late fifties, investigating this physics both theoretically~\cite{cullen1958travelling,morgenthaler1958velocity,tien1958parametric,Cassedy1963,Cassedy1967,fante1971transmission,chu1972wave,weekes2001numerical,maestre2007spatio,lurie2007introduction,maestre2008dynamic,jensen2009space,to2009homogenization,jensen2010optimization,sanguinet2011homogenized,fang2012realizing,yuan2016photonic,milton2017field,movchanPTRSMPES2022,ammariJoMP2023,Ammari2022,Nassar2017,Huidobro2019,farhat2021spacetime,Huidobro2021,Pham2022,li2023space,Touboul2024} and experimentally~\cite{ashkin1958parametric,couder2005walking,lira2012electrically,taravati2016mixer,fink2016loschmidt,mallejac2023scattering,tessierbrothelandeAPL2023,Tirole2022,Moussa2023,Lustig2023,Harwood2024,Goldsberry2025}. In particular, systems that are periodically modulated in both space and time have the ability to alleviate some of the constraints of static media, such as the breaking of reciprocity~\cite{Goldsberry2022}. As a result, exciting wave control possibilities have been identified, such as unidirectional amplification~\cite{wenCP2022}, coherent perfect absorption~\cite{Galiffi2024} and non-reciprocity~\cite{Croenne2019,Nassar2020,tessierbrothelandeAPL2023,yi_reflection_2018}.

Nevertheless, a current limitation relies on the difficulty to achieve time modulation of bulk parameters experimentally (e.g.\ density or Young's modulus of elasticity). It is easier to modify properties at discrete points, for example, by modifying the stiffness of a membrane or a surface impedance~\cite{zhuPRB2020,zhuAPL2020}, or by vertical oscillations of a submerged plate within a water tank~\cite{KOUKOURAKI2025103530}. Recent work has focused on the theoretical analysis of wave propagation in systems with time modulation at discrete points using multiple scattering~\cite{puJoSaV2024}, transmission-line theory~\cite{mallejac2023scattering}, and capacitance matrix~\cite{ammari_scattering_2024,ammari_spacetime_2025}.

The aim of this paper is then to study theoretically and numerically the propagation of waves across a time-modulated interface. For this purpose, the article is structured as follows. Section~\ref{SecContinu} describes the time-dependent jump conditions, which generalize the static jump conditions widely used to describe imperfect contacts~\cite{assierPRSMPES2020}. Interface dissipation is also incorporated to be closer to experimental devices. An energy balance is conducted. The generation of harmonics is studied through a harmonic balance analysis. The particular case of reflectionless modulated interface is discussed. Section~\ref{SecNum} describes the time-domain numerical methods used to simulate transient wave phenomena. Integration of the momentum equation and of the constitutive law is done by a fourth-order finite-difference ADER scheme. The time-varying jump conditions are discretized by the Explicit Simplified Interface Method (ESIM) \cite{lombardSJSC2003,lombardSJSC2006}, requiring new developments of this method. Validation of the numerical results is done by comparisons with a semi-analytical solution. Section~\ref{SecSimus} presents different numerical experiments, illustrating the theoretical findings: amplification of waves, generation of harmonics, impedance matching, and non-reciprocity. Lastly, Section~\ref{SecConclu} draws future lines of research such as homogenization of a network of such modulated interfaces.

\section{Physical modeling}\label{SecContinu}

\subsection{Problem statement}\label{SecContinutPb}

Let us consider a one-dimensional linear elastic medium with density $\rho(x)$ and Young's modulus $E(x)=\rho(x)\,c^2(x)$, where $c$ is the sound speed (the analysis is however relevant to other wave physics problems such as acoustics or s-polarization and p-polarization in electromagnetism, see Remark~\ref{remark_1}). These parameters may be discontinuous at $x=x_0$. The conservation of momentum and the Hooke's law are formulated as follows:
\begin{subnumcases}{\label{EDP-Hooke}}
\partial_t\parens[\big]{\rho(x)\,\partial_{t}u(x,t)} =\partial_x\sigma(x,t)+F(x,t),\label{EDP-Hooke-a} \\
\sigma(x,t) = E(x)\,\partial_x u(x,t),\label{EDP-Hooke-b}
\end{subnumcases}
with $u$ being the displacement field, $\sigma$ the stress field, and $F$ a body force with compact support that does not contain $x_0$. The velocity field is denoted by $v=\partial_t u$. For any function $g(x)$, we define the jump and mean operators $\jump{\,\cdot\,}_{x_0}$ and $\dmean{\,\cdot\,}_{x_0}$ as
\begin{equation}\label{jump_mean}
\jump{g}_{x_0}=g^+(x_0)-g^-(x_0),\qquad \dmean{g}_{x_0}=\frac{1}{2}\parens[\big]{g^+(x_0)+g^-(x_0)},
\end{equation}
where
\begin{equation}\label{Gpm}
g^\pm(x_0)=\lim_{\eta \rightarrow 0^+}g(x_0\pm \eta).
\end{equation}
To simplify the notations, the indices and arguments denoting space and time will be omitted when no risk of ambiguity occurs; in particular, $x_0$ is omitted in the jump and mean operators from now on. As can be seen from direct computations, for any two functions $g$ and $h$, these operators satisfy the property
\begin{equation}\label{PropJumpMean}
\jump{g\,h}=\jump{g}\,\dmean{h}+\dmean{g}\,\jump{h}.
\end{equation}
Introducing interface parameters of stiffness $\mathscr{K}(t)>0$, compliance $\mathscr{C}(t)=1/\mathscr{K}(t)$, inertia $\mathscr{M}(t)\geq 0$ and dissipation $\mathscr{Q}_{C,M}(t)\geq 0$, the imperfect interface is modelled by the time-dependent jump conditions:
\begin{subnumcases}{\label{JC1interf}}
\jump[\big]{v(\,\cdot\,,t)}=\partial_t\parens[\big]{\mathscr{C}(t)\,\dmeanbig{\sigma(\,\cdot\,,t)}}+\mathscr{Q}_C(t)\,\dmeanbig{\sigma(\,\cdot\,,t)},\label{JC1interf-K}
\\
\jump[\big]{\sigma(\,\cdot\,,t)}=\partial_t\parens[\big]{\mathscr{M}(t)\,\dmeanbig{v(\,\cdot\,,t)}}+\mathscr{Q}_M(t)\,\dmeanbig{v(\,\cdot\,,t)}.\label{JC1interf-M}
\end{subnumcases} Five remarks arise from~\eqref{JC1interf}:
\begin{itemize}
\item the modulated jump conditions are assumed to result from an external stimulus; such a modulation of $\jump{v}$ and $\jump{\sigma}$ can be obtained by piezoelectric components~\cite{tessierbrothelandeAPL2023} and time-modulated membranes~\cite{zhuAPL2020}, respectively;
\item as seen further in Section~\ref{SecContinuNRJ}, the compliance $\mathscr{C}(t)$ and the inertia $\mathscr{M}(t)$ must be \emph{inside} the time derivative to ensure a sound energy balance;
\item the terms in~\eqref{JC1interf-K} are reminiscent of the Maxwell model of viscoelasticity in the case of a static dissipation parameter $\mathscr{Q}_C(t)=\mathscr{Q}_C$;
\item when $\mathscr{C}=0$, $\mathscr{M}=0$, $\mathscr{Q}_C=0$ and $\mathscr{Q}_M=0$, the conditions of perfect bonding are recovered;
\item if $\mathscr{Q}_C=0$, then~\eqref{JC1interf-K} can be recast in the more usual form with displacement:
\[
\jump[\big]{u(\,\cdot\,,t)}=\mathscr{C}(t)\,\dmeanbig{\sigma(\,\cdot\,,t)}.
\]
\end{itemize}

\begin{rema}\label{remark_1}
The problem studied in this article is generic to various wave physics problems. For instance, the case of acoustics is obtained by changing $u$ into the acoustic pressure, $E$ into the inverse of mass density and $\rho$ into the compressibility. Similarly, the case of s-polarization (resp.\ p-polarization) in electromagnetism can be obtained by changing $u$ into the transverse electric (resp.\ magnetic) field, $\rho$ into the permittivity (resp.\ permeability) and $E$ into the inverse of the permeability (resp.\ inverse of the permittivity).
\end{rema}

\subsection{Energy balance}\label{SecContinuNRJ}

We take the product of the momentum equation~\eqref{EDP-Hooke-a} by $v=\partial_t u$ and integrate over $(-\infty,x_0)$. Integrating by parts and using~\eqref{EDP-Hooke-b} give
\begin{equation}\label{IPP1}
\begin{split}
\int_{-\infty}^{x_0}\rho\,v\,\partial_t v \, \dd x
	& = \int_{-\infty}^{x_0}v\,\partial_x\sigma \, \dd x+\int_{-\infty}^{x_0}F\,v \, \dd x
\\	& = -\int_{-\infty}^{x_0}\sigma\,\partial_xv \, \dd x+[v\,\sigma]_{-\infty}^{x_0}+\int_{-\infty}^{x_0}F\,v \, \dd x
\\	& = -\int_{-\infty}^{x_0}\frac{1}{E}\,\sigma\,\partial_t\sigma \, \dd x+v^-(x_0,t)\,\sigma^-(x_0,t)+\int_{-\infty}^{x_0}F\,v \, \dd x.
\end{split}
\end{equation}
Similarly, integrating over $(x_0,+\infty)$ yields
\begin{equation}\label{IPP2}
\int_{x_0}^{+\infty}\rho\,v\,\partial_t v \, \dd x= -\int_{x_0}^{+\infty}\frac{1}{E}\,\sigma\,\partial_t\sigma \, \dd x-v^+(x_0,t)\,\sigma^+(x_0,t)+\int^{+\infty}_{x_0}F\,v \, \dd x,
\end{equation}
where we have used the fact that $F$ is compactly supported and that the fields vanish at infinity. Summing~\eqref{IPP1} and~\eqref{IPP2} leads to
\begin{equation}\label{IPP3}
\frac{\dd}{\dd t} \parens[\Bigg]{\frac{1}{2}\int_{\mathbb{R}}\parens[\Bigg]{\rho\,v^2+\frac{1}{E}\,\sigma^2} \, \dd x}=-\jump{v\,\sigma}+\int_{\mathbb{R}}F\,v \, \dd x.
\end{equation}
The property~\eqref{PropJumpMean} and the jump conditions~\eqref{JC1interf} lead to
\[
\begin{split}
 \jump{v\,\sigma}
& = \jump{v}\dmean{\sigma}+\dmean{v}\jump{\sigma} 
\\	& = \parens[\big]{\mathscr{C}'(t)\,\dmean{\sigma}+\mathscr{C}(t)\,\dmean{\partial_t \sigma}+\mathscr{Q}_C(t)\,\dmean{\sigma}}\dmean{\sigma}+\parens[\big]{\mathscr{M}'(t)\,\dmean{v}+\mathscr{M}(t)\,\dmean{\partial_t v}+\mathscr{Q}_M(t)\,\dmean{v}}\dmean{v}
\\ 	& = \frac{\dd}{\dd t} \parens[\Bigg]{\frac{1}{2}\parens[\big]{\mathscr{C}(t)\,\dmean{\sigma}^2 + \mathscr{M}(t)\,\dmean{v}^2}}+\frac{1}{2}\parens[\big]{\mathscr{C}'(t)\dmean{\sigma}^2+\mathscr{M}'(t)\dmean{v}^2}+\mathscr{Q}_C(t)\,\dmean{\sigma}^2+\mathscr{Q}_M(t)\,\dmean{v}^2.
\end{split}
\]
It leads to the following energy balance.

\begin{Proposition}\label{PropNRJ}
Let ${\mathcal E}_m={\mathcal E}_b+{\mathcal E}_i$ be the total energy, with the bulk energy
\begin{equation}\label{NRJv}
{\mathcal E}_b(t)=\frac{1}{2}\int_{\mathbb{R}}\parens[\Bigg]{\rho\,v^2+\frac{1}{E}\,\sigma^2} \, \dd x
\end{equation}
and the interface energy
\begin{equation}\label{NRJi}
{\mathcal E}_i(t)=\frac{1}{2}\mathscr{M}\,\dmean{v}^2+\frac{1}{2}\mathscr{C}\,\dmean{\sigma}^2.
\end{equation}
Introducing the power of external forces due to body forcing
\begin{equation}
{\mathcal P}(t)=\int_{\mathbb{R}}F\,v \, \dd x,
\end{equation}
one has
\begin{equation}\label{dEdT}
\frac{\dd}{\dd t}{\mathcal E}_m(t)={\mathcal P}(t)-\frac{1}{2}\parens[\big]{\mathscr{M}'(t)\,\dmean{v}^2+ \mathscr{C}'(t)\,\dmean{\sigma}^2}-\parens[\big]{\mathscr{Q}_C(t)\,\dmean{\sigma}^2+\mathscr{Q}_M(t)\,\dmean{v}^2}.
\end{equation}
\end{Proposition}

Five remarks arise from Proposition~\ref{PropNRJ} and positivity of $\mathscr{M}$, $\mathscr{C}$, $\mathscr{Q}_C$ and $\mathscr{Q}_M$:
\begin{itemize}
\item both ${\mathcal E}_b$ and ${\mathcal E}_i$ are positive, hence ${\mathcal E}_m \geq 0$;
\item in ${\mathcal E}_i$ one recognizes the kinetic energy and the potential energy of a spring-mass system when static parameters are considered;
\item if there are no body force ($F=0$), no time-modulation (both $\mathscr{M}$ and $\mathscr{C}$ are constant in time) and no dissipation term in the jump conditions (${\mathscr{Q}_C}=0$ and ${\mathscr{Q}_M}=0$), then the total mechanical energy is conserved;
\item if no body force is applied ($F=0$) and both $\mathscr{M}$ and $\mathscr{C}$ are constant in time, but ${\mathscr{Q}_C}$ or ${\mathscr{Q}_M}$ is non-zero, then the total mechanical energy is dissipated;
\item if $F=0$ but either $\mathscr{M}$ or $\mathscr{C}$ varies with time, then the total mechanical energy may vary, due to the power of external forces required to modify $\mathscr{M}$ and $\mathscr{C}$; the sign of the right-hand-side in~\eqref{dEdT} is arbitrary, so that ${\mathcal E}_m$ may increase or decrease.
\end{itemize}

\subsection{Some additional properties} The latter remark about the amplification of energy raises the question of possible parametric amplification. Under suitable assumptions, the next proposition states that such an amplification is not possible, and that the solution remains bounded. The proof is given in Appendix~\ref{AppAmpli}.

\begin{Proposition}\label{PropAmpli}
Let us assume:
\begin{enumerate}\romanenumi
\item a sinusoidal modulation of the interface compliance $\mathscr{C}(t)$ and the dissipation term $\mathscr{Q}_C(t)$;
\item no jump of stress (i.e.\ $\mathscr{M}(t)=0$ and $\mathscr{Q}_M(t)=0$);
\item a body force in~\eqref{EDP-Hooke-a} in the form $F(x,t)=\delta(x-x_s)\,S(t)$, where $S$ is bounded.
\end{enumerate}
Then the scattered fields $v$ and $\sigma$ remain bounded for all $t$.
\end{Proposition}

This result is counter-intuitive, since modulated systems often exhibit resonance (e.g.\ a modulated spring-mass system). An interpretation is that the energy may be evacuated on both sides of the domain, preventing from an unbounded increase. A similar conclusion holds by interchanging the roles of $\mathscr{C}$ and $\mathscr{M}$. The case where both interface compliance and interface inertia are modulated remains open but this limitation seems, \emph{a priori}, purely technical. Indeed, our numerical experiments suggest that the same result will hold in this case.

\begin{Proposition}\label{PropEnvelop}
When only $\mathscr{Q}_C(t)\neq 0$, i.e.\ $\mathscr{M}(t)=0$, $\mathscr{C}(t)=0$ and $\mathscr{Q}_M(t)=0$, the limits between which the solutions evolve are determined analytically (Appendix~\ref{AppAnalytic}) and correspond to the static cases with the extreme values of $\mathscr{Q}_C(t)$. The limits for the case with only $\mathscr{Q}_M(t)\neq0$ can be obtained using the same reasoning.
\end{Proposition}

\subsection{Generation of harmonics}\label{SecContinuHarm}

The interaction of a single-frequency wave with a modulated interface generates an infinite number of harmonics. Here we shortly describe how to compute them. For simplicity, constant physical parameters are considered around the interface at $x_0=0$. In this part, a sinusoidal modulation for $\mathscr{C}$, $\mathscr{M}$, $\mathscr{Q}_C$ and $\mathscr{Q}_M$ is chosen:
\begin{equation}\label{ModulSinus}
\begin{aligned}
& \mathscr{C}(t)=\mathscr{C}_0\parens[\big]{1+\varepsilon_C\,\sin(\Omega t)},\\
& \mathscr{M}(t)=\mathscr{M}_0\parens[\big]{1+\varepsilon_M\,\sin(\Omega t)},\\
& \mathscr{Q}_C(t)=\mathscr{Q}_{C_0}\parens[\big]{1+\varepsilon_{Q_C}\,\sin(\Omega t)},\\
& \mathscr{Q}_M(t)=\mathscr{Q}_{M_0}\parens[\big]{1+\varepsilon_{Q_M}\,\sin(\Omega t)},
\end{aligned}
\end{equation}
with $\Omega=2\,\pi\,f_m$, where $f_m$ is the modulation frequency, and $\mathscr{C}_0\geq0$, $\mathscr{M}_0\geq0$ and $-1<\varepsilon_X<1$, with $X=C,M,Q_C,Q_M$. A harmonic incident wave $v^{\incident}$ and $\sigma^{\incident}$ of the form
\[
v^{\incident} (x, t) = e^{i \omega (t - x / c)},\qquad \sigma^{\incident} (x, t) = -\rho c \,e^{i \omega (t - x / c)},
\]
impacts the interface from the left. Thus, we can write
\begin{align*}
v^{\total} (x, t) & =
\begin{dcases}
v^{\incident} (x, t) + v^{\reflected} (x, t) 	& \text{for $x < 0$},\\
v^{\transmitted} (x, t) 						& \text{for $x > 0$},
\end{dcases}
\intertext{and}
\sigma^{\total} (x, t) & =
\begin{dcases}
\sigma^{\incident} (x, t) + \sigma^{\reflected} (x, t) 	& \text{for $x < 0$},\\
\sigma^{\transmitted} (x, t) 							& \text{for $x > 0$},
\end{dcases}
\end{align*}
for some reflected and transmitted wave fields $v^{\reflected}$, $\sigma^{\reflected}$, $v^{\transmitted}$ and $\sigma^{\transmitted}$ that also satisfy the governing wave equation. Floquet theorem states that these fields write
%\begin{alignat}{3}
%v^{\reflected} (x, t) & = \sum_{k \in \mathbb{Z}} R_k e^{ i \omega_k (t + x / c)}, 	& \qquad \hfill \sigma^{\reflected} (x, t) & = \rho c\sum_{k \in \mathbb{Z}} R_k e^{ i \omega_k (t + x / c)}, \label{eq:ansatz-ref}
%\\
%v^{\transmitted} (x, t) & = \sum_{k \in \mathbb{Z}} T_k e^{ i \omega_k (t - x / c)},	& \qquad \sigma^{\transmitted} (x, t) & = - \rho c \sum_{k \in \mathbb{Z}} T_k e^{ i \omega_k (t - x / c)}, \label{eq:ansatz-trans}
%\end{alignat}
\begin{gather}
v^{\reflected} (x, t) = \sum_{k \in \mathbb{Z}} R_k e^{ i \omega_k (t + x / c)}, 	\qquad \sigma^{\reflected} (x, t) = \rho c\sum_{k \in \mathbb{Z}} R_k e^{ i \omega_k (t + x / c)}, \label{eq:ansatz-ref}
\\
v^{\transmitted} (x, t) = \sum_{k \in \mathbb{Z}} T_k e^{ i \omega_k (t - x / c)},	\qquad \sigma^{\transmitted} (x, t) = - \rho c \sum_{k \in \mathbb{Z}} T_k e^{ i \omega_k (t - x / c)}, \label{eq:ansatz-trans}
\end{gather}
for some reflection and transmission coefficients $R_k$ and $T_k$ to be determined, and where we denote $\omega_k \equiv \omega + k \Omega$. One introduces the auxiliary quantities $\Psi_k = T_k - R_k$ and $\Phi_k = T_k + R_k$.

\subsubsection{Dealing with the first jump equation}

Using~\eqref{eq:ansatz-ref} and~\eqref{eq:ansatz-trans}, together with the definition of the jump operator, we find that
\[
\llbracket v^{\total} \rrbracket = e^{ i \omega t} \parens[\Bigg]{- 1 +
\sum_{k \in \mathbb{Z}} \Psi_k e^{ik \Omega t} } \qquad \text{and} \qquad \llangle \sigma^{\total} \rrangle = - \frac{\rho\,c}{2} e^{ i \omega t} \parens[\Bigg]{1+\sum_{k \in \mathbb{Z}} \Psi_k e^{ik \Omega t} }.
\]
Similarly, though it is slightly longer, we find that
\begin{multline*}
\partial_t \dmeanbig{\mathscr{C}(t)\, \sigma^{\total}} = e^{i
\omega t} \parens[\Bigg]{\frac{\mathcal{C} \varepsilon_C}{4}\omega_{-1} e^{- i \Omega t} -
\frac{i\mathcal{C}}{2}\omega - \frac{\mathcal{C} \varepsilon_C}{4} \omega_1e^{+ i
\Omega t}} \\
+ e^{i \omega t}\sum_{k \in \mathbb{Z}} \parens[\Bigg]{- \frac{\mathcal{C}
\varepsilon_C}{4} \omega_{k} \Psi_{k - 1} - \frac{i\mathcal{C}}{2}
\omega_{k} \Psi_k + \frac{\mathcal{C} \varepsilon_C}{4} \omega_{k} \Psi_{k + 1} } e^{ ik \Omega t},
\end{multline*}
where we have introduced the reduced quantity $\mathcal{C} = \rho\,c\,\mathscr{C}_0$. Therefore the first jump condition~\eqref{JC1interf-K}, after dividing through by $e^{i \omega t}$ and rearranging slightly becomes
\begin{multline}
 \sum_{k \in \mathbb{Z}} \parens[\Bigg]{\frac{\mathcal{C} \varepsilon_C\omega_{k}-i\mathcal{Q}_C\varepsilon_{Q_C}}{4} \Psi_{k - 1} + \parens[\Bigg]{1 +\frac{\mathcal{Q}_C}{2}+ i\frac{\mathcal{C}}{2} \omega_k } \Psi_k -\frac{\mathcal{C} \varepsilon_C\omega_{k}-i\mathcal{Q}_C\varepsilon_{Q_C}}{4} \Psi_{k +1} } e^{ik\Omega t}
\\ = \frac{\mathcal{C} \varepsilon_C\omega_{-1}-i\mathcal{Q}_C\varepsilon_{Q_C}}{4}\,e^{- i \Omega t} + \parens[\Bigg]{1-\frac{\mathcal{Q}_C}{2} - i\frac{\mathcal{C} \omega}{2} } - \frac{\mathcal{C} \varepsilon_C\omega_1-i\mathcal{Q}_C\varepsilon_{Q_C}}{4}\,e^{+ i \Omega t},
\end{multline}
where $\mathcal{Q}_C = \rho\,c\,\mathscr{Q}_{C_0}$. Since the $e^{ ik \Omega t}$ functions are linearly independent, we obtain this system in matrix form:
\begin{equation}\label{SysPsi}
\mathbb{A} (\mathcal{C}, \varepsilon_C, \mathcal{Q}_C, \varepsilon_{Q_C}, \omega, \Omega) \,\bfPsi = \bfV (\mathcal{C}, \varepsilon_C, \mathcal{Q}_C, \varepsilon_{Q_C}, \omega, \Omega),
\end{equation}
where
\begin{align*}
	& \bfPsi \equiv \parens{\dotsc, \Psi_{-N}, \dotsc, \Psi_{-1}, \Psi_0, \Psi_{ 1}, \dotsc, \Psi_{ N}, \dotsc}^\top,
\\	& \bfV (\mathcal{C}, \varepsilon_C, \mathcal{Q}_C, \varepsilon_{Q_C},
\omega,\Omega) =
\\	& \textstyle \hspace{2cm} \parens[\big]{\dotsc, \, 0, \dotsc, \frac{\mathcal{C} \varepsilon_C}{4}\omega_{-1}-i\frac{\mathcal{Q}_C\varepsilon_{Q_C}}{4},\, 1-\frac{\mathcal{Q}_C}{2} - i\frac{\mathcal{C}}{2}\omega,\, -\frac{\mathcal{C} \varepsilon_C}{4}\omega_1+i\frac{\mathcal{Q}_C\varepsilon_{Q_C}}{4}, \dotsc, \, 0, \dotsc }^\top,
\end{align*}
and $\mathbb{A}$ is a tri-banded antidiagonal matrix. Now we can truncate this system so that $\mathbb{A}$ becomes a $(2 N + 1) \times (2 N + 1)$ matrix, and the truncated vector $\bfPsi $ can be recovered by simple inversion.

\subsubsection{Dealing with the second jump equation}

One obtains directly
\[
\llbracket \sigma^{\total} \rrbracket = \rho c\,e^{i \omega t} \parens[\Bigg]{1 - \sum_{k \in \mathbb{Z}} \Phi_k e^{ ik \Omega t} } \qquad \text{and} \qquad \llangle v^{\total} \rrangle=\frac{1}{2} e^{i\omega t} \parens[\Bigg]{1+\sum_{k\in\mathbb{Z}} \Phi_k e^{ ik \Omega t} }.
\]
Similarly, we find that
\begin{multline*}
\partial_t \dmeanbig{\mathscr{M} (t) {v^{\total}}} = e^{i
\omega t} \parens[\Bigg]{- \frac{\mathscr{M}_0 \,\varepsilon_M}{4}\omega_{-1} e^{- i \Omega t} +
\frac{i\mathscr{M}_0}{2}\omega + \frac{\mathscr{M}_0 \, \varepsilon_M}{4}\omega_1 e^{+ i
\Omega t}}
\\	+ e^{i \omega t}\sum_{k \in \mathbb{Z}} \parens[\Bigg]{\frac{\mathscr{M}_0\,
\varepsilon_M}{4} \omega_{k} \Phi_{k - 1} + \frac{i\mathscr{M}_0}{2}
\omega_{k} \Phi_k - \frac{\mathscr{M}_0\,\varepsilon_M}{4} \omega_{k} \Phi_{k + 1} } e^{ ik \Omega t}.
\end{multline*}
Upon introducing $\mathcal{M} = \mathscr{M}_0/(\rho c)$ and $\mathcal{Q}_M=\mathscr{Q}_{M_0}/(\rho c)$, the second jump condition~\eqref{JC1interf-M} leads to
\begin{multline}
\sum_{k \in \mathbb{Z}} \parens[\Bigg]{\frac{\mathcal{M}\,\varepsilon_M\omega_{k}-i\mathcal{Q}_M\varepsilon_{Q_M}}{4} \Phi_{k - 1} + \parens[\Bigg]{1 +\frac{\mathcal{Q}_M}{2}+ i\frac{\mathcal{M}}{2} \omega_k } \Phi_k -\frac{\mathcal{M}\,\varepsilon_M\omega_{k}-i\mathcal{Q}_M\varepsilon_{Q_M}}{4} \Phi_{k +1} } e^{ik\Omega t}
\\ = \frac{\mathcal{M}\, \varepsilon_M\omega_{-1}-i\mathcal{Q}_M\varepsilon_{Q_M}}{4}\,e^{- i \Omega t} + \parens[\Bigg]{1 -\frac{\mathcal{Q}_M}{2} - i\frac{\mathcal{M} \omega}{2}} - \frac{\mathcal{M}\,\varepsilon_M\omega_1-i\mathcal{Q}_M\varepsilon_{Q_M}}{4}\,e^{+ i \Omega t}.
\end{multline}
We end up with a system of the form
\begin{equation}\label{SysPhi}
\mathbb{B} (\mathcal{M}, \varepsilon_{M},\mathcal{Q}_M, \varepsilon_{Q_M}, \omega, \Omega) \,\bfPhi = \bfW (\mathcal{M}, \varepsilon_M,\mathcal{Q}_M, \varepsilon_{Q_M},\omega,\Omega).
\end{equation}
The matrix $\mathbb{B}$ and the vector $\bfW $ are obtained by replacing $(\mathcal{C}, \varepsilon_C, \mathcal{Q}_C, \varepsilon_{Q_C}, \omega, \Omega)$ by $(\mathcal{M}, \varepsilon_M, \mathcal{Q}_M, \varepsilon_{Q_M},\omega,\Omega)$ in $\mathbb{A}$ and $\bfV $. The vector $\bfPhi $ is
\[
(\dotsc, \Phi_{-N}, \dotsc, \Phi_{-1}, \Phi_0, \Phi_{ 1}, \dotsc, \Phi_{ N},\dotsc)^\top.
\]
As before, this system can be solved for any truncation by simple inversion.

\subsubsection{Final reflection and transmission coefficients}

Once $\Psi_k$ and $\Phi_k$ are computed, the reflection and transmission coefficients are given by:
\begin{equation}\label{RkTk}
R_k = \frac{1}{2} (\Phi_k - \Psi_k),\qquad T_k = \frac{1}{2} (\Psi_k + \Phi_k).
\end{equation}
Two limit-cases are highlighted.
\begin{itemize}
\item Without modulation, the only non-zero coefficients are $R_0$ and $T_0$. These coefficients are given in Appendix~\ref{AppRT}.
\item If $\mathcal{M}=\mathcal{C}$, $\varepsilon_M = \varepsilon_C$, $\mathcal{Q}_C=\mathcal{Q}_M$ and $\varepsilon_{Q_C} = \varepsilon_{Q_M}$, then the two matrix systems~\eqref{SysPsi} and~\eqref{SysPhi} are exactly the same, and so $\bfPhi =\bfPsi $, which leads to $R_k = 0$ for all $k$. This property of \emph{impedance matching} is summarized in the following proposition.
\end{itemize}

\begin{Proposition}\label{PropImpedance}
Let us consider a homogeneous medium of impedance $Z=\rho\,c$. Under the following equalities
\begin{equation}\label{Z_KMT}
\mathscr{M}_0=Z^2\,\mathscr{C}_0, \qquad \mathscr{Q}_{M_0}=Z^2\,\mathscr{Q}_{C_0}, \qquad
\varepsilon_C=\varepsilon_M, \qquad
\varepsilon_{\mathcal{Q}_C} = \varepsilon_{\mathcal{Q}_M},
\end{equation}
then no reflected wave is generated at the modulated interface.
\end{Proposition}

Proposition~\ref{PropImpedance} generalizes the impedance matching condition already obtained with static jump conditions and no dissipation in~\cite{assierPRSMPES2020}.

\begin{rema}\label{RqueImpedance}
Proposition~\ref{PropImpedance} can be extended to any periodic modulation, while respecting $\mathscr{M}(t)=Z^2\mathscr{C}(t)$ and $\mathscr{Q}_M(t)=Z^2\mathscr{Q}_C(t)$, using the linearity of the equations and the fact that periodic functions can be decomposed into Fourier series.
\end{rema}

\section{Numerical modeling}\label{SecNum}

\subsection{Numerical scheme}\label{SecNumScheme}

\begin{enonce*}[remark]{First-order formulation}
A velocity-stress formulation of the evolution equations is used. We introduce the compact notations
\begin{equation}\label{SysOrdre1}
\begin{alignedat}{3}
	& \bfU (x,t)=
\begin{pmatrix}
v(x,t)\\
\sigma(x,t)
\end{pmatrix}, \qquad
	&& \bfF (x,t)=
\begin{pmatrix}
\frac{F(x,t)}{\rho(x)}\\
0
\end{pmatrix}, \qquad
	&
\\	& \bfA (x)=
\begin{pmatrix}
0 & -\frac{1}{\rho(x)}\\
-E(x) & 0
\end{pmatrix}, \qquad
	&& \bfB (t)=
\begin{pmatrix}
0 & \mathscr{C}(t)\\
\mathscr{M}(t) & 0
\end{pmatrix}, \qquad
	& \bfE (t)=
\begin{pmatrix}
0 & \mathscr{Q}_C(t)\\
\mathscr{Q}_M(t) & 0
\end{pmatrix},
\end{alignedat}
\end{equation}
where $\bfA $ and $\bfE $ can be discontinuous at the interface. From~\eqref{EDP-Hooke} and from~\eqref{JC1interf}, one obtains the boundary-value problem
\begin{subnumcases}{\label{BVP}}
\partial_t \bfU +\bfA \,\partial_x\bfU =\bfF ,	& $x\neq x_0$,\label{BVP-EDP} \\
\jump{\bfU (\,\cdot\,,t)}=\partial_t\parens[\big]{\bfB (t)\,\dmean{\bfU (\,\cdot\,,t)}}+\bfE (t)\,\dmean{\bfU (\,\cdot\,,t)}. & \label{BVP-JC}
\end{subnumcases}
It is assumed that the body force is zero at the interface: $\bfF (x_0,t)=\bfzero $, which is not restrictive.
\end{enonce*}

\begin{enonce*}[remark]{ADER-4 scheme}
A uniform grid with space step $\Delta x$ and time step $\Delta t$ is introduced, to determine an approximation $\bfU _j^n$ of $\bfU (x_j=j\,\Delta x, t_n=n\,\Delta t)$. Integration of~\eqref{BVP-EDP} is done by a fourth-order ADER scheme. This explicit two-step finite-difference scheme is fourth-order accurate in space and time, using a stencil of 2 neighboring points on the right and on the left. It is stable under the usual CFL condition
\begin{equation}\label{CFL}
\zeta=\max\parens[\big]{c(x)}\,\frac{\Delta t}{\Delta x}\leq 1,
\end{equation}
where $c(x)=\sqrt{\frac{E(x)}{\rho(x)}}$ is the sound speed. The reader is referred to~\cite{benjaziaWM2014} for implementation details.
\end{enonce*}

\subsection{Immersed interface method for time-modulated jump conditions}\label{SecNumIIM}

Discretization of the modulated jump conditions~\eqref{BVP-JC} is done by implementing the Explicit Simplified Interface Method (ESIM). The philosophy of this numerical method is as follows: based on the jump conditions and the interface geometry, smooth extensions of the solution~$\bfU $ on either side of the interface are constructed at each time step (Figure~\ref{fig:IIM}). Evaluation of these extensions leads to modified values $\bfU ^\star_j$ of the solution. At a point $x_j$ on one side of the interface, the integration scheme (ADER 4 here) then uses the usual values $\bfU _j^n$ if $x_j$ is on the same side, and $\bfU _j^\star$ if $x_j$ is on the other side of the interface. This method has been applied to various interface problems, notably to constant imperfect interfaces in 1D~\cite{lombardSJSC2003} and 2D~\cite{lombardSJSC2006}.

\begin{figure}[h!]
\includegraphics[width=0.65\textwidth]{IIM}
\caption{Scheme of the principle of ESIM at an interface with a discontinuous solution~$\bfU (x,t)$ (blue line) and its smooth extensions $\bfU ^\star(x,t)$ (red dotted lines).}\label{fig:IIM}
\end{figure}

An important ingredient of ESIM is to determine the jump conditions satisfied by the solution and its successive spatial derivatives. As the jump conditions here vary as a function of time, an important modification has to be made to the method, as detailed now. From~\eqref{jump_mean} and~\eqref{BVP-JC}, we obtain
\begin{equation}\label{Esim1}
\bfU ^+-\frac{1}{2}\,\partial_t(\bfB\,\bfU^+)-\frac{1}{2}\,\bfE \,\bfU ^+=\bfU ^-+\frac{1}{2}\,\partial_t(\bfB\,\bfU^-)+\frac{1}{2}\,\bfE \,\bfU ^-,
\end{equation}
where the superscripts $\pm$ refer to the limit values, as in~\eqref{Gpm}. Using~\eqref{BVP-EDP}, we then have
\begin{equation}\label{Esim2}
\parens[\Bigg]{\bfI -\frac{1}{2}\parens{\partial_t\bfB +\bfE }}\,\bfU ^++\frac{1}{2}\bfB\,\bfA^+\,\partial_x\bfU ^+=\parens[\Bigg]{\bfI +\frac{1}{2}\parens{\partial_t\bfB +\bfE }}\,\bfU ^--{\frac{1}{2}}\bfB\,\bfA^-\,\partial_x\bfU ^-,
\end{equation}
where $\bfI $ is the $2\times 2$ identity matrix. Successive time-derivatives of~\eqref{Esim2} up to an order $k$ yields the matrix relation
\begin{equation}\label{Esim3}
\bfC_k^+(t)\,\bfU _k^++\bfGamma_{k+1}^+=\bfC_k^-(t)\,\bfU _k^-+\bfGamma_{k+1}^-,
\end{equation}
with the $2\,(k+1)$ vectors of limit values
\begin{equation}
\begin{aligned}
	& \bfU _k^{\pm}=\parens{\bfU ^{\pm},\,\partial_x\bfU ^{\pm},\dotsc,\,\partial^k_x\bfU ^{\pm}}^\top,
\\	& \bfGamma_{k+1}^\pm=\parens[\Big]{\bfzero ,\,\bfzero ,\dotsc,\,\bfzero ,\,\mp\tfrac{(-1)^{k+1}}{2}\,\bfB \parens{\bfA ^\pm}^{k+1}\,\partial_x^{k+1}\bfU ^\pm}^\top,
\end{aligned}
\end{equation}
where $\bfC_k^\pm$ are $\parens[\big]{2(k+1)}\times \parens[\big]{2(k+1)}$ time-dependent matrices, $\bfzero $ is the $2 \times 2$ zero matrix, and $\bfGamma_{k+1}^\pm$ are $2(k+1)$ vectors. The matrices $\bfC_k^\pm$ depend on the bulk parameters (through the matrices $\bfA ^\pm$) and the successive time derivatives of the jump conditions (through the matrix $\bfB (t)$ and $\bfE (t)$ in~\eqref{SysOrdre1}). Moreover, the matrices $\bfC^\pm$ have a structure of block lower diagonal matrix; the $2 \times 2$ blocks below the diagonal depend on the successive time derivatives of $\mathscr{C}(t)$, $\mathscr{M}(t)$, $\mathscr{Q}_C(t)$ and $\mathscr{Q}_M(t)$.

Neglecting the $(k+1)$-th spatial derivatives ($\bfGamma_{k+1}^\pm=\bfzero _{2(k+1)}$) and inverting~\eqref{Esim3} yields
\begin{equation}\label{Esim4}
\begin{split}
\bfU _k^+ 
	& = \parens[\big]{\bfC_k^+(t)}^{-1}\,\bfC_k^-(t)\,\bfU _k^-
\\	& \coloneqq \bfD_k(t)\,\bfU _k^-.
\end{split}
\end{equation}
The transfer matrices $\bfD_k(t)$ of order $k$~\eqref{Esim4} are the building-block of the ESIM. These are used to calculate extrapolation matrices, which are used to determine the modified values $\bfU ^\star_i$ at each time step. Once computed, the rest of the method follows the usual lines. The reader is referred e.g.\ to~\cite{touboulJoCP2020} for practical details. Four remarks can be made concerning the implementation and the properties of the ESIM.
\begin{itemize}
\item The transfer matrices $\bfD_k$ in~\eqref{Esim4} depend on time. Consequently, the extrapolation matrices must be recalculated at each time step. This is an important difference with the cases previously dealt with by ESIM, where the extrapolation matrices were calculated once during a preprocessing step.
\item The calculation of $\bfD_k$ is tedious, especially for large values of $k$. A naive idea to simplify this calculation is to \emph{freeze} the jump conditions at the instant under consideration, so that $\partial_t \bfB =\bfzero $ in~\eqref{Esim2}, and $\bfD_k(t)$ depends only on $\bfB $ and $\bfE $ and not on their successive derivatives. This is a bad idea: the time evolution of the modulated interface conditions is then poorly described, leading to inaccurate numerical results.
\item The problem under study has similarities with other time-varying interface problems. But here we are in a \emph{kinematic} case, where the jump conditions are known and imposed, unlike \emph{dynamic} cases such as~\cite{assierJFM2014}, where the jump conditions are unknown and depend on the solution in a non-trivial way.
\item The optimal choice of the derivation order $k$ is important. In the case of static jump conditions (see e.g.~\cite{lombardSJSC2003}), a numerical analysis of the ESIM has been carried out, and we will reuse these results here. In the limit case where the interface disappears ($\mathscr{C}_0= 0$, $\mathscr{M}_0=0$, $\mathscr{Q}_{C_0}=0$ and $\mathscr{Q}_{M_0}=0$), $\bfU _j^\star=\bfU _j^n$ is needed to recover the scheme in homogeneous medium. This consistency property is obtained if $k$ is odd and $k \geq 2s-1$, where $s$ is the width of the stencil. Additionally, the local truncation error of a $r$-th order accurate scheme is maintained at the interface if $k\geq r$. Putting together these two conditions yields
\begin{equation}\label{Esim-AN}
k=\max\parens{r,2s-1} \quad \text{and} \quad k \ \text{odd}.
\end{equation}

Here $r=4$ and $s=2$, so that $k=5$ is the optimal value, which will be used in the numerical experiments. We will examine numerically the convergence properties using this value in Section~\ref{ValidNum} since we have no proof that it still holds in the current modulated case.
\end{itemize}

\subsection{Numerical set-up}\label{SecSimusSetup}

We consider a medium of length $L=400$~\unit{m}, discretized on $N_x=400$ grid nodes. The bulk parameters are constant and correspond to those of a bar of Plexiglass: $\rho=1200$~\unit{kg/m^3} and $c=2800$~\unit{m/s}. The CFL number in~\eqref{CFL} is $\zeta=0.95$.

A modulated interface lies at $x_0=200$ m. The values of the interface parameters used further have no physical origin, and they are only chosen to illustrate clearly the expected phenomena. Typically, modulation magnitudes of the order of 0.75 or 0.9 are used, whereas real systems lead to magnitudes of the order of 0.2~\cite{mallejac2023scattering}.

Two types of forcing are considered: (i) a Cauchy problem with an initial data associated with a right-going wave
\begin{equation}\label{Cauchy}
\bfU (x,0) =
\begin{pmatrix}
\rho\\
-1/c
\end{pmatrix}
\,
S\parens[\Bigg]{t_0-\frac{x}{c}},
\end{equation}
where $t_0$ allows to position the initial wave on the left of the interface, or (ii) zero initial conditions with a time-dependent forcing at a Dirac source point
\begin{equation}\label{Dirac}
F(x,t)=S(t)\,\delta(x-x_s).
\end{equation}
Both forcings involve the source function $S$, chosen as a combination of truncated sinusoids:
\begin{equation}\label{Source}
S(\xi)=
\begin{dcases}
\sum_{m=1}^{4}a_m\sin(b_m\,\omega_c\,\xi) 	& \text{if $-0<\xi <1/f_c$}, \\
0 											& \text{otherwise},
\end{dcases}
\end{equation}
where $f_c=\omega_c/(2\pi)$ is the central frequency, $b_m=2^{m-1}$, the coefficients $a_m$ are $a_1=1$, $a_2=-{21}/{32}$, $a_3 = {63}/{768}$, and $a_4=-{1}/{512}$. This source entails $C^6\parens[\big]{[0,+\infty\mathclose{[}}$ smoothness. Figure~\ref{FigSource} illustrates \eqref{FigSource_a} a Cauchy problem and \eqref{FigSource_b} the spectral content of $S$.

\begin{figure}[h!]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_t0}
\caption{}\label{FigSource_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Source-FFT}
\caption{}\label{FigSource_b}
\end{subfigure}
\caption{Forcing of simulations. \eqref{FigSource_a} Cauchy problem~\eqref{Cauchy} with $f_c=45$~\unit{Hz}. \eqref{FigSource_b} Spectrum of the source~\eqref{Source}; the vertical dotted line denotes the central frequency $f_c$.}\label{FigSource}
\end{figure}

The time dependence of $\mathscr{C}$, $\mathscr{M}$, $\mathscr{Q}_C$ and $\mathscr{Q}_M$ writes
\begin{equation}\label{KM-T}
\begin{aligned}
\mathscr{C}(t)		& = \mathscr{C}_0\parens[\big]{1+\varepsilon_C\,\phi_C(t)},\\
\mathscr{M}(t)		& = \mathscr{M}_0\parens[\big]{1+\varepsilon_M\,\phi_M(t)},\\
\mathscr{Q}_C(t)		& = \mathscr{Q}_{C_0}\parens[\big]{1+\varepsilon_{Q_C}\,\phi_{Q_C}(t)},\\
\mathscr{Q}_M(t)		& = \mathscr{Q}_{M_0}\parens[\big]{1+\varepsilon_{Q_M}\,\phi_{Q_M}(t)},
\end{aligned}
\end{equation}
where $\mathscr{C}_0\geq0$, $\mathscr{M}_0\geq0$, $\mathscr{Q}_{C_0}\geq0$, $\mathscr{Q}_{M_0}\geq0$, and $\varepsilon_X\,\phi_X >-1$, with $X=C,M,Q_C,Q_M$. The functions $\phi_C$, $\phi_M$, $\phi_{Q_C}$ and $\phi_{Q_M}$ may differ from each other. If $\phi_X=0$ or $\varepsilon_X=0$, then $\mathscr{C}$, $\mathscr{M}$, $\mathscr{Q}_C$ and $\mathscr{Q}_M$ do not depend on $t$, and one recovers the static conditions of imperfect contact investigated in~\cite{lombardSJSC2003}. In addition to the sinusoidal modulation~\eqref{ModulSinus}, two modulation functions are introduced for the numerical simulations: a quasi-periodic modulation and a rectangular one. It follows
\begin{equation}\label{FuncModul}
\phi_{X}(t)=
\begin{dcases}
\sin(\Omega t) & \text{(sinusoidal)},\\
\sin(\Omega t)+\sin\parens{\sqrt{2}\,\Omega t} & \text{(quasi-periodic)},\\
(-1)^n, \ \text{with} \ n= \floor[\Big]{\parens[\Big]{\tfrac{\Omega \,t}{2\pi} \mod 1} - \nu} \quad & \text{(rectangular)},
\end{dcases}
\end{equation}
where $\lfloor g \rfloor$ denotes the floor function of $g$, $0>\nu>1$, and $\Omega=2\pi\,f_m$, where $f_m$ is the \emph{modulation frequency}. In this work, we assume $\phi_C=\phi_M=\phi_{Q_C}=\phi_{Q_M}$ (which implies that the modulation frequencies are the same), but this is only for simplicity. The numerical scheme presented above does not depend on this assumption. Therefore, this scheme and our implementation of it remain valid for functions that differ from one another.

\subsection{Validation}\label{ValidNum}

In this part, a validation of the numerical scheme is proposed.

\begin{figure}[htpb]
\begin{subfigure}{0.49\linewidth}
\includegraphics[width=\textwidth]{F0}
\caption{}\label{FigModulK_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{F0-eta2m7-mod.pdf}
\caption{}\label{FigModulK_b}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{F100}
\caption{}\label{FigModulK_c}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{F100-eta2m7-mod.pdf}
\caption{}\label{FigModulK_d}
\end{subfigure}
\caption{Scattering by an interface with modulated stiffness, for various frequencies of sinusoidal modulation $f_m$. Snapshot of the solution at $t=0.045$~\unit{s} with: \eqref{FigModulK_a}--\eqref{FigModulK_b} no modulation; \eqref{FigModulK_c}--\eqref{FigModulK_d} $f_m=100$~\unit{Hz}. Left: cases without dissipation. Right: $\mathscr{Q}_{C_0}=2\,\times 10^{-7}$~\unit{m^2\cdot s/kg} and $\mathscr{Q}_{M_0}=0$.}\label{FigModulK}
\end{figure}

\begin{figure}[htpb]
\setcounter{figure}{2}
\renewcommand{\thefigure}{\arabic{figure} (cont.)}
\begin{subfigure}{0.49\textwidth}
\setcounter{subfigure}{4}
\includegraphics[width=\textwidth]{F500}
\caption{}\label{FigModulK_e}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{F500-eta2m7-mod.pdf}
\caption{}\label{FigModulK_f}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{F500-zoom}
\caption{}\label{FigModulK_g}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{F500-eta2m7-mod-zoom.pdf}
\caption{}\label{FigModulK_h}
\end{subfigure}
\caption{Scattering by an interface with modulated stiffness, for various frequencies of sinusoidal modulation $f_m$. \eqref{FigModulK_e}--\eqref{FigModulK_f} Snapshot of the solution at $t=0.045$~\unit{s} with $f_m=500$~\unit{Hz}. \eqref{FigModulK_g}--\eqref{FigModulK_h} Zoom on the reflected wave when $f_m=500$~\unit{Hz}. Left: cases without dissipation. Right: $\mathscr{Q}_{C_0}=2\,\times 10^{-7}$~\unit{m^2\cdot s/kg} and $\mathscr{Q}_{M_0}=0$.}\label{FigModulKbis}
\end{figure}

A pulse is initially located on the left of the interface~\eqref{Cauchy}. The interface parameters are $\mathscr{K}_0=2.45$~\unit{GPa/m}, $\mathscr{M}_0=0$. The central frequency is $f_c=30$~\unit{Hz}. A semi-analytic solution of~\eqref{BVP} can be obtained by the method of characteristics when $\mathscr{M}_0=0$ or $\mathscr{C}_0=0$. Technical details are given in Appendix~\ref{AppAnalytic} in the case where only the stiffness is modulated.

Figures~\ref{FigModulK}\eqref{FigModulK-NH_a} and~\ref{FigModulK}\eqref{FigModulK_b} illustrate the case without modulation: $\varepsilon_C=0$. The other subfigures are computed using $\varepsilon_C=\varepsilon_{Q_C}= 0.9$ and various values of the modulation frequency $f_m=\frac{\Omega}{2\pi}$. Cases without dissipation in the interface conditions ($\mathscr{Q}_{C_0}=0$ and $\mathscr{Q}_{M_0}=0$) are presented on the left column and cases with modulated dissipation ($\mathscr{Q}_{C_0}=2\,\times 10^{-7}$~\unit{m^2\cdot s/kg} and $\mathscr{Q}_{M_0}=0$) are on the right side. The value of $\mathscr{Q}_{C_0}$ have been chosen to have an effect of the same order as the compliance effect. One observes an excellent agreement between the exact and numerical values, which validates the numerical method detailed in Section~\ref{SecNum}. There is also an enrichment of peaks in the scattered waves when $\abs{f_m-f_c}$ increases.

The particular case where only $\mathscr{Q}_{C}$ is non-zero and modulated is presented on Figure~\ref{FigModulQ_AN}. Here again, the agreement between the analytical solution and the result of the numerical model is excellent. On the left side, the envelope curves has been added to highlight the Proposition~\ref{PropEnvelop}. In this case with a sinusoidal modulation, the extreme values of $\mathscr{Q}_{C}(t)$, used to determine the dotted curves, are equal to $\mathscr{Q}_{C_{0}}(1\pm\varepsilon_{Q_C})$.

\begin{figure}[htpb]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Q_ananumZoom}
\caption{}\label{FigModulQ_AN_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Q_Envelop2}
\caption{}\label{FigModulQ_AN_b}
\end{subfigure}
\caption{Scattering by an interface with modulated dissipation ($\mathscr{Q}_{C_0}=2\times10^{-7}$~\unit{m^2\cdot s/kg} and the other parameters are zero) at the frequency $f_m=500$~\unit{Hz} at $t=0.045$~\unit{s}. \eqref{FigModulQ_AN_a} Snapshot of the solution with a zoom on the reflected wave. \eqref{FigModulQ_AN_b} Analytic solution with envelope curves (dashed lines).}\label{FigModulQ_AN}
\end{figure}

In the rest of this work, we only consider cases with $\mathscr{Q}_{C_0}=0$ and $\mathscr{Q}_{M_0}=0$.

\begin{figure}[htpb]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{FQP-100}
\caption{}\label{FigModulK-NH_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Fsquare-100}
\caption{}\label{FigModulK-NH_b}
\end{subfigure}
\caption{Scattering by an interface with modulated stiffness, for non-harmonic modulation at the frequency $f_m=100$~\unit{Hz} at $t=0.045$~\unit{s}. \eqref{FigModulK-NH_a} Quasi-periodic modulation. \eqref{FigModulK-NH_b} Square modulation.}\label{FigModulK-NH}
\end{figure}

Figure~\ref{FigModulK-NH} illustrates the two other modulation types given by~\eqref{FuncModul} at the frequency $f_m=100$~\unit{Hz}. Figure~\ref{FigModulK-NH}\eqref{FigModulK-NH_a} displays the results obtained in the quasi-periodic modulation, while Figure~\ref{FigModulK-NH}\eqref{FigModulK-NH_b} investigates the rectangular modulation ($\nu=0.65$ in~\eqref{FuncModul}). Here again, there is an excellent agreement between the exact and numerical solution. In the case of a square modulation, the terms corresponding to the derivatives of the modulation in~\eqref{Esim4} are set to zero.

Measures of convergence have been performed to verify the conclusions obtained in~\cite{lombardSJSC2003}. The Figure~\ref{figConvergence} presents the evolution of the error $\varepsilon_v$ committed on the velocity field $v$ at $t=0.045$~\unit{s}, and defined by:
\begin{equation}
\varepsilon_v=\sqrt{\Delta x\sum_{i=0}^{N_x}\, \parens{v(x_i,t)-v_i^n}^2}.
\end{equation}
On the left side, convergence measurement are shown for a sinusoidal modulation of the compliance $\mathscr{C}$ using $\mathscr{K}_0=2.45$~\unit{GPa/m} and $\mathscr{M}=0$. On the right side, the analog case is presented (i.e.\ $\mathscr{M}_0=2\,\times 10^4$~\unit{kg/m^2} and $\mathscr{C}=0$). The other numerical values of the parameters are $\varepsilon_{C,M}=0.75$, $f_c=45$~\unit{Hz} and two different frequencies of a sine modulation are tested: $f_m=100$~\unit{Hz} (plain blue line) and $f_m=10$~\unit{Hz} (orange dashed line).

\begin{figure}[htpb]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Convergence-C.pdf}
\caption{}\label{figConvergence_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Convergence-M}
\caption{}\label{figConvergence_b}
\end{subfigure}
\caption{Convergence measurements for $f_m =100$~\unit{Hz} in blue line and $f_m=10$~\unit{Hz} in orange dashed line in the case where \eqref{figConvergence_a} $\mathscr{C}(t)\neq 0$ and $\mathscr{M}(t)=0$, and \eqref{figConvergence_b} $\mathscr{C}(t)=0$ and $\mathscr{M}(t)\neq 0$. The green dotted line shows a slope of order 4 (in log-log scale).}\label{figConvergence}
\end{figure}

These measurements show a convergence of order $4$ for both modulations and both frequencies. The two curves corresponding to $f_m=100$~\unit{Hz} present a slope discontinuity for the higher value of $\Delta x$ but we only consider the limit for small values of $\Delta x$. Moreover, due to the generation of harmonics with larger frequencies, the value of the error is more significant for a modulation frequency at $f_m=100$~\unit{Hz} than at $f_m=10$~\unit{Hz}.

\section{Numerical experiments}\label{SecSimus}

Here we illustrate the various properties analyzed in Section~\ref{SecContinu}. We conclude with an illustration of non-reciprocity.

\subsection{Modulation of energy}\label{SecSimusNRJ}

\begin{figure}[htpb]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_t1_inc}
\caption{}\label{FigNinterfAmpli_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_energy_inc}
\caption{}\label{FigNinterfAmpli_b}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_t1_dec}
\caption{}\label{FigNinterfAmpli_c}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_energy_dec}
\caption{}\label{FigNinterfAmpli_d}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_t1_squaremod}
\caption{}\label{FigNinterfAmpli_e}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_energy_squaremod}
\caption{}\label{FigNinterfAmpli_f}
\end{subfigure}
\caption{Interaction of a pulse centered on $f_c=45$~\unit{Hz} with a modulated interface. Left: snapshot of the scattered velocity $v$ at $t=0.045$~\unit{s}. Right: time evolution of the energy, where the gray shaded areas denote the evolution of the interface parameters with modulation frequency $f_m=100$~\unit{Hz}. Sinusoidal modulation using \eqref{FigNinterfAmpli_a}--\eqref{FigNinterfAmpli_b} $\varepsilon_{C,M}=0.75$, \eqref{FigNinterfAmpli_c}--\eqref{FigNinterfAmpli_d} $\varepsilon_{C,M}=-0.75$, and a rectangular modulation \eqref{FigNinterfAmpli_e}--\eqref{FigNinterfAmpli_f} with $\varepsilon_{C,M}=0.75$ and $\nu=0.65$.}\label{FigNinterfAmpli}
\end{figure}

A pulse of the form~\eqref{Cauchy} with central frequency $f_c=45$~\unit{Hz} is initially located on the left of the interface. The interface parameters are $\mathscr{K}_0=2.45$~\unit{GPa/m} and $\mathscr{M}_0=2\,\times 10^4$~\unit{kg/m^2} with a modulation frequency $f_m=100$~\unit{Hz}. Figure~\ref{FigNinterfAmpli} investigates the scattered waves for different modulations. The left column shows a snapshot of the velocity $v$. The right column shows the temporal evolution of the energies~\eqref{NRJv}--\eqref{NRJi} as a function of time. The background gray scale represents the time evolution of the interface parameters.

The total energy ${\mathcal E}_m$ is constant until $t\approx 0.01$~\unit{s}, when the incident wave impacts the interface. The initially zero interface energy ${\mathcal E}_i$ then increases and remains non-zero as long as the wave/interface interaction lasts. As a consequence, ${\mathcal E}_m$ also varies. After the wave has passed through the interface, the interface energy is radiated and reconverted into bulk energy, which is then conserved. Depending on the sign of $\varepsilon_{C,M}$, different evolutions of energy are observed.
\begin{itemize}
\item In the first case $\varepsilon_{C,M}=0.75$ (\eqref{FigNinterfAmpli_a}--\eqref{FigNinterfAmpli_b}), the final total energy value ($=0.0013$) has almost doubled compared with the initial energy ($=0.007$).
\item In the second case, the modulation is inverted by choosing $\varepsilon_{C,M}=-0.75$ (\eqref{FigNinterfAmpli_c}--\eqref{FigNinterfAmpli_d}). The total amount of energy decreases from $=0.007$ to $=0.006$.
\item In the third case, a rectangular modulation is applied with $\varepsilon_{C,M}=0.75$ and $\nu=0.65$ (\eqref{FigNinterfAmpli_e}--\eqref{FigNinterfAmpli_f}). The interface energy (and consequently the total energy) is then discontinuous.
\end{itemize}
In the first case, the interface is reached when the value of $\varepsilon_{C,M} \sin(\Omega t)$ is at its maximum value, whereas these values are negative in the second case. In conclusion, the work used to modulate the interface parameters can logically increase or decrease the total energy, as expected and predicted in our remarks to Proposition~\ref{PropNRJ}.

\begin{figure}
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Deltaener}
\caption{}\label{fig:deltaener_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{DeltaenerFMZoomRat}
\caption{}\label{fig:deltaener_b}
\end{subfigure}
\caption{Influence of the modulation parameters on the variation of energy after crossing the interface. \eqref{fig:deltaener_a} Influence on the amplitude of the modulation $\varepsilon_{C,M}$. \eqref{fig:deltaener_b} Influence of the ratio $f_m/f_c$.}\label{fig:deltaener}
\end{figure}

Figure~\ref{fig:deltaener}\eqref{fig:deltaener_a} presents the variation of the energy $\Delta{\mathcal E}=\mathcal{E}_m(t)-\mathcal{E}_m(0)$ between $t=0$~s and $t=0.04$~s (time after which the wavelet has completely crossed the interface) for different values of the amplitude of modulation $\varepsilon_{C,M}\in [-0.95,0.95]$, and for four different modulation frequencies ($f_m=50$~\unit{Hz}, $f_m=100$~\unit{Hz}, $f_m=150$~\unit{Hz} and $f_m=750$~\unit{Hz}). This figure shows that the gain of energy is not linear with respect to $\varepsilon_{C,M}$. Moreover, for low values of $f_m$, there exists a range of modulation for which energy is lost, but this range depends also on the choice of the modulation frequency. For the highest frequencies, the variation of energy is always positive and is symmetric because the phase shift due to the modification of sign of $\varepsilon_{C,M}$ does not have a significant effect. Figure~\ref{fig:deltaener}\eqref{fig:deltaener_b} shows the variation of the energy $\Delta{\mathcal E}$ as a function of $\frac{f_m}{f_c}$ for three values of $\varepsilon_{C,M}$. At lower frequencies (until $f_m=400$~\unit{Hz}, i.e.\ $\frac{f_m}{f_c}\approx 9$, not presented here), the variation of energy~$\Delta{\mathcal E}$ depends strongly on the modulation frequency $f_m$. This is because, at such low modulation frequency, the state of the interface when the pulse crosses it varies a lot with $f_m$. On the contrary, at higher frequencies (presented here), $\Delta{\mathcal E}$ appears to tend monotonously to a limit value, that could be interpreted as an effective (averaged in time) behavior.

\subsection{Generation of harmonics}\label{SecSimusHarmonics}

Here we illustrate the calculations carried out in Section~\ref{SecContinuHarm}. The reflection and transmission coefficients $R_k$ and $T_k$ in~\eqref{eq:ansatz-ref} are determined by solving~\eqref{RkTk} with $N=12$ Fourier modes. This number has been chosen experimentally to ensure the convergence. The interface parameters are $\mathscr{K}_0=2.45$~\unit{GPa/m} and $\mathscr{M}_0=2\,\times 10^4$~\unit{kg/m^2}, with modulation amplitudes $\varepsilon_C=\varepsilon_M=0.75$.

\begin{figure}[htbp]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{HBM_R_fc100_fm30}
\caption{}\label{FigHarmonics_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{HBM_T_fc100_fm30}
\caption{}\label{FigHarmonics_b}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{HBM_R_fc30_fm30}
\caption{}\label{FigHarmonics_c}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{HBM_T_fc30_fm30}
\caption{}\label{FigHarmonics_d}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{HBM_R_fc30_fm100}
\caption{}\label{FigHarmonics_e}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{HBM_T_fc30_fm100}
\caption{}\label{FigHarmonics_f}
\end{subfigure}
\caption{Generation of harmonics. Coefficients of reflection $R_k$ (left) and transmission $T_k$ (right) in~\eqref{eq:ansatz-ref}. The vertical dotted line denotes the frequency $f_c$ of the incident wave. \eqref{FigHarmonics_a}--\eqref{FigHarmonics_b} $f_c=100$ and $f_m=30$. \eqref{FigHarmonics_c}--\eqref{FigHarmonics_d} $f_c=f_m=30$. \eqref{FigHarmonics_e}--\eqref{FigHarmonics_f} $f_c=30$ and $f_m=100$. Blue lines denote the FFT of direct simulations, whereas plain circles and plain stars denote the results of harmonic balance respectively for $\omega_k<0$ and $\omega_k\geq0$.}\label{FigHarmonics}
\end{figure}

The forcing is a Dirac source point~\eqref{Dirac} on the left of the interface $(x_s<x_0$). The velocity field is recorded at a receiver on the right of the interface ($x_r>x_0$), and then a FFT is performed. Three cases are studied considering different frequencies of the source ($f_c=\omega/(2\,\pi)$) and of the modulation ($f_m=\Omega/(2\,\pi)$). Their values are (i) $f_c=100$~\unit{Hz} and $f_m=30$~\unit{Hz} in the first case, (ii) $f_c=f_m=30$~\unit{Hz} in the second one, and (iii) $f_c=30$~\unit{Hz} and $f_m=100$~\unit{Hz} in the third one.

Figure~\ref{FigHarmonics} compares the Fourier spectra of the reflected and transmitted signals obtained by FFT of direct simulations (blue lines) with the theoretical values obtained by harmonic balance (plain circles and stars). We represent here only the positive part of the frequencies as the spectrum is symmetric. The theoretical coefficient corresponding to negative (resp.\ positive) values of $\omega_k$ are then plotted in green circles (resp.\ purple stars). A vertical dotted line is set at the frequency $f_c$. In the three cases, the higher coefficient corresponds to the fundamental mode ($k=0$) and the coefficient values decrease when $\abs{k}$ increase. In each case, an excellent agreement is obtained between results of simulation and of harmonic balance, even for a small number of Fourier modes. The following observations are done.
\begin{itemize}
\item When $f_c>f_m$, the influence of negative frequencies is weak (almost zero) because it corresponds to higher modes. Moreover, the values of the coefficients decrease rapidly. The coefficients of the first harmonics are smaller than the half of the value at the frequency $\omega$.
\item When $f_c=f_m$, then $\omega_{-1}=0$ and the calculations show that $R_{-1}= T_{-1}=0$. In this case, green circles and purple stars are located at the same frequency because $\omega_{-k}=\omega_{k-2}$ and in practice, the coefficients computed for $\omega_k\leq 0$ are strictly zero and the same observation applies if $f_c$ is a multiple of $f_m$. This is an expected behavior in the light of~\cite[Section~4]{KOUKOURAKI2025103530}.
\item When $f_c<f_m$, there is a stronger influence of the modes at negative frequencies. It can be explained by the fact that, as visible on the Figure~\ref{FigModulK}, the signal is strongly modified with a high-frequency modulation. The first harmonic for the reflection coefficient is for example almost equal to the fundamental frequency ($\approx0.29$), and there are twelve modes that can be noticed while only five or six were visible for a smaller modulation frequency.
\end{itemize}

\subsection{Impedance matching}\label{SecSimusImpedance}

The interface is modulated sinusoidally at $f_m=100$~\unit{Hz} and its parameters are $\mathscr{M}_0=12\,\times 10^3$~\unit{kg/m^2}, $\mathscr{K}_0=941$~\unit{MPa/m}, so that the properties~\eqref{Z_KMT} are verified. Both parameters are modulated at the same frequency and with the same amplitude $\varepsilon_{C,M}=0.75$. Bulk parameters are given in Section~\ref{SecSimusSetup}.

\begin{figure}[htbp]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_t1_RL}
\caption{}\label{FigNinterfRL_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{FFTSigs}
\caption{}\label{FigNinterfRL_b}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_energy_RL}
\caption{}\label{FigNinterfRL_c}
\end{subfigure}
\caption{Reflectionless interaction of a pulse centered on $f_c=45$~\unit{Hz} with a modulated interface due to impedance adaptation. \eqref{FigNinterfRL_a} Snapshot of the scattered velocity $v$ at $t=0.045$~\unit{s}. \eqref{FigNinterfRL_b} Spectra of the incident wave and of the transmitted wave with harmonic generation; the vertical dotted line denotes $f_c$. \eqref{FigNinterfRL_c} Time evolution of the energy, where the gray shaded areas denote the evolution of the interface parameters. The modulation function is sinusoidal.}\label{FigNinterfRL}
\end{figure}

Figure~\ref{FigNinterfRL} illustrates the behavior of a wave at this modulated interface. On the bottom left side, a snapshot of $v$ is presented. At the left of the interface, no reflected wave propagates because of the impedance matching. The shape of the transmitted signal is similar to that obtained with the first configuration in Figure~\ref{FigNinterfAmpli}\eqref{FigNinterfAmpli_a}. The spectrum of the transmitted wave is presented in Figure~\ref{FigNinterfRL}\eqref{FigNinterfRL_b} and compared to the spectrum of the source. The harmonic generation is then visible. As previously, the total amount of energy is modified when the pulse crosses the interface (Figure~\ref{FigNinterfRL}\eqref{FigNinterfRL_c}).

Proposition~\ref{PropImpedance} about the impedance matching has only been proven for the case of sinusoidal modulation. Nevertheless, Figure~\ref{FigNinterfRLNP} shows that the results can be extended to other functions defined by~\eqref{FuncModul} if $\phi_C(t)=\phi_M(t)$. Two cases are tested: \eqref{FigNinterfRLNP_a} a quasi-periodic modulation with $f_m=100$~\unit{Hz}, and \eqref{FigNinterfRLNP_b} a rectangular modulation with $f_m=100$~\unit{Hz} and $\nu=0.65$. The latter case corresponds to the framework of Remark~\ref{RqueImpedance}. In both cases, no reflection is observed.

\begin{figure}[htbp]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_t1_RL_QP}
\caption{}\label{FigNinterfRLNP_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{1Interface_t1_RL_rect}
\caption{}\label{FigNinterfRLNP_b}
\end{subfigure}
\caption{Snapshots of reflectionless interaction of a pulse centered on $f_c=45$~\unit{Hz} with a modulated interface at $t=0.045$~\unit{s}. Impedance adaptation \eqref{FigNinterfRLNP_a} with a quasi-periodic modulation and \eqref{FigNinterfRLNP_b} with a rectangular modulation.}\label{FigNinterfRLNP}
\end{figure}

\subsection{Non-reciprocity}\label{SecQSimusNR}

Lastly, we investigate numerically the non-reciprocity induced by time-modulating the parameters of the interface. The configuration is the same as in the previous Section. Figure~\ref{FigNonReciprocity} shows a comparison of temporal signals measured when we swap source and sensor (Figure~\ref{FigNonReciprocity}\eqref{FigNonReciprocity_a}) for different values of $f_c$ and $f_m$. The positions of the sources and sensors are $x_s=150$~\unit{m} and $x_r=225$~\unit{m}, respectively at a distance $\ell_s$ and $\ell_r$ from the interface.

\begin{figure}[htpb]
\begin{subfigure}{0.49\textwidth}
\centering
\includegraphics[width=.65\textwidth]{ConfigNR}
\caption{}\label{FigNonReciprocity_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{NR_fc30_fm30_Zoom}
\caption{}\label{FigNonReciprocity_b}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{NR_fc100_fm30_Zoom}
\caption{}\label{FigNonReciprocity_c}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{NR_fc30_fm100_Zoom}
\caption{}\label{FigNonReciprocity_d}
\end{subfigure}
\caption{Illustration of non-reciprocity. \eqref{FigNonReciprocity_a} Scheme of the numerical experiment where the positions of the sensor and the source are inverted. The signals measured by the sensor are represented for \eqref{FigNonReciprocity_b} $f_m=f_c=30$~\unit{Hz}, \eqref{FigNonReciprocity_c} $f_c=100$~\unit{Hz} and $f_m=30$~\unit{Hz}, and \eqref{FigNonReciprocity_d} $f_c=30$~\unit{Hz} and $f_m=100$~\unit{Hz}.}\label{FigNonReciprocity}
\end{figure}

The transmission of waves through the interface is altered because of the change of state of the interface, when the wave reaches it. It is due to the difference in travel time from the source to the interface. In both cases, the measured signals are different, which illustrates the fact that the reciprocity theorem is no longer verified.

We propose in this section to quantify the non-reciprocity by computing the difference between signals. Figure~\ref{FigQuantNonReciprocity} presents a measure of the difference between the two signals $v_s$, measured by the sensor at $x=x_r$ due to the source in $x=x_s$, and $v_r$, measured by the sensor at $x=x_s$ due to the source in $x=x_r$ at $t=0.07$~\unit{s}. This quantity is defined by
\begin{equation}
\vartheta_v=f_c\sqrt{\sum_{i=0}^{N_x}(v_s-v_r)^2}.
\end{equation}
Two types of modulation are studied, a sinusoidal one~\eqref{ModulSinus} and a quasi-periodic one~\eqref{FuncModul} for 100 values of the modulation frequency $f_m$ from $0$ to $400$~\unit{Hz}. For both the periodic and quasi-periodic modulations, two different frequencies of the pulse are tested: $f_c=30$~\unit{Hz} represented with the plain blue line and $f_c=50$~\unit{Hz} with orange dashed line. As the wavelength differs in the two cases, $\vartheta_v$ is normalized using the source frequency.

\begin{figure}[htpb]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{NRP2}
\caption{}\label{FigQuantNonReciprocity_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{NRQP2}
\caption{}\label{FigQuantNonReciprocity_b}
\end{subfigure}
\caption{Quantification of non-reciprocity for various modulation frequencies in~\eqref{FuncModul}. \eqref{FigQuantNonReciprocity_a}~Pe\-ri\-od\-ic modulation with two different values of $f_c$: blue line for 30~\unit{Hz} and orange dashed line for 50~\unit{Hz}. \eqref{FigQuantNonReciprocity_b} Quasi-periodic modulation with $f_c=30$~\unit{Hz} and $f_c=50$~\unit{Hz}.}\label{FigQuantNonReciprocity}
\end{figure}

Measurement with a periodic modulation, as the sinusoidal one treated on Figure~\ref{FigQuantNonReciprocity}\eqref{FigQuantNonReciprocity_a}, shows that some frequencies exist, for a couple source/sensor positions, for which $\vartheta_v=0$. The link between these values is $\ell_s-\ell_r=p \frac{c}{f_m}$, with $p\in\mathbb{Z}$. For these frequencies, the interface is in the same state (due to being phase-shifted by a temporal period) when the incident wave crosses the interface. In particular, the case $p=0$ corresponds to a symmetric case and the measures of non-reciprocity are always zero. It is important to notice that the modulation frequencies with zero non-reciprocity do not depend on the source frequency. One can notice that the non-reciprocity tends to increase with the ratio $\frac{f_m}{f_c}$.

Figure~\ref{FigQuantNonReciprocity}\eqref{FigQuantNonReciprocity_b} corresponds to a quasi-periodic modulation. Except for the case $\ell_s=\ell_r$, there are no longer values of $f_m\neq 0$ for which $\vartheta_v=0$, so the measures of non-reciprocity are always strictly positive and the signals are always different. Similarly to the case of a periodic modulation in Figure~\ref{FigQuantNonReciprocity}\eqref{FigQuantNonReciprocity_a}, the non-reciprocity tends to increase with the ratio $\frac{f_m}{f_c}$, albeit at a slower rate and with a more complex pattern.

\section{Conclusion}\label{SecConclu}

In this study, we investigated the interaction of waves with an interface whose jump conditions are modulated in time. These jump conditions can be experimentally implemented by considering e.g.\ membranes, spring-mass systems or vibrating plates whose thickness is small compared to the wavelength. Numerical methods for time-domain simulations have been developed. Numerical experiments have demonstrated and validated theoretically predicted properties: amplification of the energy, harmonic generation, impedance matching, and non-reciprocity.

Under the assumptions of Proposition~\ref{PropAmpli} (only one sinusoidally modulated interface parameter, the other parameter being zero), we proved that no parametric amplification resonance can occur. That said, we were not able to conclude on this when these assumptions were relaxed. A more complete theoretical study of the energy evolution remains to be carried out, especially in order to better understand the high-frequency regime.

While mentioning resonance, it is worth noting here that we have only considered non-resonant jump conditions. A relaxation of this assumption could be envisaged, and would allow, for example, to study the case of Helmholtz resonators whose physical or geometrical properties vary with time~\cite{ammariJoMP2023,mallejac2023scattering}.

Another natural follow-up to this work is to consider a set of modulated interfaces, in phase or not, and periodically distributed in space. In the case where the wavelength is large compared with the spacing between interfaces, asymptotic homogenization can be applied. The work carried out in~\cite{bellisJotMaPoS2021} for the case of non-linear but static interfaces is currently generalized to the case of modulated interfaces. This work, currently under finalization, will be the subject of a forthcoming publication. It could also be interesting to study the effective properties induced by a network of time-modulated resonators, as considered in~\cite{ammari_effective_2024,tachet2024effective}.

\appendix

\section{Scattering coefficients in the static case}\label{AppRT}

In the case of static jump conditions ($\mathscr{C}(t)=\mathscr{C}_0$, $\mathscr{M}(t)=\mathscr{M}_0$, $\mathscr{Q}_{C}(t)=\mathscr{Q}_{C_0}$ and $\mathscr{Q}_{M}(t)=\mathscr{Q}_{M_0}$), the scattering coefficients can be calculated analytically. The physical parameters are assumed to be piecewise constant:
\[
(\rho,\,c)=
\begin{dcases}
(\rho_0,\,c_0)	& \text{in $\Omega_0 \ \ (x<x_0)$},\\
(\rho_1,\,c_1)	& \text{in $\Omega_1 \ \ (x>x_0)$}.
\end{dcases}
\]
Upon introducing the following quantities
\begin{equation}\label{NotaImpedance}
\begin{gathered}
Z_0=\rho_0\,c_0,\qquad Z_1=\rho_1\,c_1,\qquad Z_2=\rho_0\,c_1,\\
Y_0=\mathscr{C}_0\,Z_0\,Z_1-\mathscr{M}_0,\qquad Y_1=\mathscr{C}_0\,Z_0\,Z_1+\mathscr{M}_0,\qquad \omega_\sharp^2=\frac{\mathscr{C}_0\mathscr{M}_0}{4},
\end{gathered}
\end{equation}
the reflection coefficient $R$ and the transmission coefficient $T$ at frequency $\omega$ are given by
\[
\begin{aligned}
& R(\omega)=\frac{(Z_1 \! - \! Z_0)\,\parens[\Big]{1 \! - \! \parens[\big]{\frac{\omega}{\omega_\sharp}}^2 \! + \! \frac{\mathscr{Q}_{C_0}\,\mathscr{Q}_{M_0}}{4}} \! - \! \mathscr{Q}_{C_0}\,Z_0\,Z_1 \! + \! \mathscr{Q}_{M_0} \! - \! i\,\omega\parens[\big]{Y_0 \! + \! \frac{1}{4}\parens{\mathscr{Q}_{C_0}\,\mathscr{M}_0 \! + \! \mathscr{Q}_{M_0}\,\mathscr{C}_0}\,(Z_1 \! - \! Z_0)}}{(Z_1 \! + \! Z_0)\,\parens[\Big]{1 \! - \! \parens[\big]{\frac{\omega}{\omega_\sharp}}^2 \! + \! \frac{\mathscr{Q}_{C_0}\,\mathscr{Q}_M}{4}} \! + \! \mathscr{Q}_{C_0}\,Z_0\,Z_1 \! + \! \mathscr{Q}_{M_0} \! + \! i\,\omega\parens[\big]{Y_1 \! + \! \frac{1}{4}\parens{\mathscr{Q}_{C_0}\,\mathscr{M}_0 \! + \! \mathscr{Q}_{M_0}\,\mathscr{C}_0}\,(Z_1 \! + \! Z_0)}},\\ 
& T(\omega)=\frac{2\,Z_2\,\parens[\Big]{1 \! + \! \parens[\big]{\frac{\omega}{\omega_\sharp}}^2 \! - \! \frac{\mathscr{Q}_{C_0}\,\mathscr{Q}_{M_0}}{4} \! - \! i\,\omega\frac{1}{4}\parens{\mathscr{Q}_{C_0}\,\mathscr{M}_0 \! + \! \mathscr{Q}_{M_0}\,\mathscr{C}_0}}}{(Z_1 \! + \! Z_0)\,\parens[\Big]{1 \! - \! \parens[\big]{\frac{\omega}{\omega_\sharp}}^2 \! + \! \frac{\mathscr{Q}_{C_0}\,\mathscr{Q}_{M_0}}{4}} \! + \! \mathscr{Q}_{C_0}\,Z_0\,Z_1 \! + \! \mathscr{Q}_{M_0} \! + \! i\,\omega\parens[\big]{Y_1 \! + \! \frac{1}{4}\parens{\mathscr{Q}_{C_0}\,\mathscr{M}_0 \! + \! \mathscr{Q}_{M_0}\,\mathscr{C}_0}\,(Z_1 \! + \! Z_0)}}.
\end{aligned}
\]

In the particular case where the density and the Young's modulus are constant, the impedance $Z=\rho\,c=\sqrt{\rho\,E}$ is constant. Additionally, let us assume that the two following conditions are satisfied
\begin{equation}
\mathscr{M}_0=Z^2\,\mathscr{C}_0, \qquad \mathscr{Q}_{M_0}=Z^2\,\mathscr{Q}_{C_0}.
\end{equation}
Then $R(\omega)=0$: no reflected wave is generated at a static imperfect interface, at any frequency. One recovers the impedance matching condition given in Proposition~\ref{PropImpedance}. Moreover, these coefficients can also be obtained by calculating $R_0$ and $T_0$ using~\eqref{RkTk} when the interface parameters are static.

\section{Analytical solution of scattered fields}\label{AppAnalytic}

Here we determine the exact solution of~\eqref{BVP} with initial data~\eqref{Cauchy}. As in Appendix~\ref{AppRT}, the physical parameters are assumed to be piecewise constant. We consider the case where only $\mathscr{C}(t)$ and $\mathscr{Q}_C(t)$ are modulated, whereas the interface inertia and the corresponding dissipation term is taken to be zero ($\mathscr{M}(t)=0$ and $\mathscr{Q}_M(t)=0$). The reciprocal case $\mathscr{M}(t)\neq 0$, $\mathscr{Q}_C(t)=0$ and $\mathscr{C}(t)=0$ can be treated in a similar manner. In both cases, the calculation relies on the method of characteristics, illustrated on Figure~\ref{fig:Riemann}. The proof follows the lines of~\cite{lombardJoCaAM2007}.

\begin{figure}[htbp]
\includegraphics[width=0.5\textwidth]{Caracteristics.pdf}
\caption{Sketch of the method of characteristics for computing semi-analytical solution using Riemann invariants.}\label{fig:Riemann}
\end{figure}

To express $\bfU (x,\,t)$ in terms of limit-values of the fields at $x_0$, we use the \emph{Riemann invariants}~$J^{R,L}$ that are constant along the \emph{characteristics} $\gamma_{R,L}$. The invariants for linear Partial Differential Equations with constant coefficients satisfy
\begin{equation}\label{Jrl}
\begin{dcases}
\gamma_R:\,\frac{\dd x}{\dd t}=+c \; \Longrightarrow \; \restr{\frac{\dd J^R}{\dd t}}{\gamma_R}=0,	& \text{with} \ J^R(x,\,t)=\frac{1}{2}\parens[\Bigg]{v-\frac{1}{\rho\,c}\sigma}(x,\,t),\\
\gamma_L:\,\frac{\dd x}{\dd t}=-c \; \Longrightarrow \; \restr{\frac{\dd J^L}{\dd t}}{\gamma_L}=0,	& \text{with} \ J^L(x,\,t)=\frac{1}{2}\parens[\Bigg]{v+\frac{1}{\rho\,c}\sigma}(x,\,t).
\end{dcases}
\end{equation}
Using the continuity of stress~\eqref{JC1interf-M} and the invariant~\eqref{Jrl}, along with the initial data condition of compact support in $\Omega_0$, one obtains for $t\geq t_0$
\begin{equation}\label{SUpm}
\begin{aligned}
& \sigma^\pm(t)=-\rho_1\,c_1\,v^+(t),\\
& v^-(t)=-\frac{\rho_1\,c_1}{\rho_0\,c_0}\,v^+(t)+2\,J_0^R\parens{x_0-c_0\,t,\,0},
\end{aligned}
\end{equation}
where the subscript $i$ on $J^{R,L}_i$ refers to $\Omega_i$. The solution $\bfU (x,\,t)$ can be expressed in terms of $v^+(s)$, with $t_0\leq s \leq t$, and of the initial values of the Riemann invariants. Introducing the travel times
\[
t_A =t-\frac{1}{c_0}(x_0-x), \qquad t_B=t-\frac{1}{c_1}(x-x_0),
\]
the solution $\bfU (x,\,t)$ is given in $\Omega_0$ by
\begin{equation}
\begin{gathered}
\bfU (x,\,t)=
\begin{pmatrix}
1 & 1\\
-\rho_0\,c_0 & \rho_0\,c_0
\end{pmatrix}
\,
\begin{pmatrix}
J_0^R\parens{x-c_0\,t,\,0}\\
\Delta_A(x,\,t)
\end{pmatrix},\\
\text{with} \quad \Delta_A(x,\,t)=
\begin{cases}
-\frac{\rho_1\,c_1}{\rho_0\,c_0}\,v^+(t_A)+J_0^R(x_0-c_0\,t_A,\,0)	& \text{if $t_A \geq 0$},\\
J_0^L(x+c_0\,t,\,0) 													& \text{otherwise},
\end{cases}
\end{gathered}
\end{equation}
and in $\Omega_1$ by
\begin{equation}\label{U(x,t)}
\begin{gathered}
\bfU (x,\,t)=
\begin{pmatrix}
1\\
-\rho_1\,c_1
\end{pmatrix}
\Delta_B(x,\,t),\\
\text{with} \quad \Delta_B(x,\,t)=
\begin{dcases}
v^+(t_B)		& \text{if $t_B \geq t_0$},\\
0			& \text{otherwise}.
\end{dcases}
\end{gathered}
\end{equation}
To complete the analytical solution, it remains to determine $v^+(t)$. For this purpose, the jump of displacement~\eqref{JC1interf-K} is differentiated, leading to
\begin{equation}
v^+(t)-v^-(t)=\parens[\big]{\mathscr{C}(t)}'\sigma^+(t)+\mathscr{C}(t)\,\partial_t\sigma^+(t)+\mathscr{Q}_C(t)\,\sigma^+(t).
\end{equation}
From~\eqref{SUpm}, it follows
\begin{equation}\label{eq:eqdify}
\parens[\Bigg]{1+\frac{\rho_1\,c_1}{\rho_0\,c_0}+\rho_1c_1\parens[\big]{\mathscr{C}'(t)+\mathscr{Q}_C(t)}}\,v^+(t)-2\,J_0^R\parens{x_0-c_0\,t,\,0}=-\rho_1\,c_1\,\mathscr{C}(t)\,\partial_tv^+(t).
\end{equation}
Setting
\begin{equation}
y(t)=v^+(t),\qquad h(t)=2\,J_0^R\parens{x_0-c_0\,t,\,0},
\end{equation}
one obtains the non-autonomous Ordinary Differential Equation if $\mathscr{C}(t)\neq0$:
\begin{equation}\label{ODE}
y'(t)=\frac{1}{\rho_1\,c_1\,\mathscr{C}(t)}\parens[\Bigg]{h(t)-\parens[\Bigg]{1+\frac{\rho_1\,c_1}{\rho_0\,c_0}+\rho_1\,c_1\parens[\big]{\mathscr{C}'(t)+\mathscr{Q}_C(t)}}\,y(t)},
\end{equation}
which is integrated by the usual Runge--Kutta 4 method. If $\mathscr{C}(t)=0$ and $\mathscr{Q}_C\neq 0$, \eqref{eq:eqdify} gives:
\begin{equation}
v^+(t)=\frac{2\,J_0^R\parens{x_0-c_0\,t,\,0}}{1+\frac{ \rho_1\,c_1}{\rho_0\,c_0}+\rho_1c_1\mathscr{Q}_C(t)}.
\end{equation}
This solution lies within the bounds defined by the solutions associated with the extreme (minimum and maximum) values of $\mathscr{Q}_C(t)$.

The reciprocal case $\mathscr{M}(t)\neq 0$, $\mathscr{Q}_M\neq 0$ and $\mathscr{C}(t)=0$, $\mathscr{Q}_C(t)=0$ can be treated using the following expressions:
\begin{gather}
\Delta_A(x,\,t)=
\begin{cases}
-\frac{1}{\rho_1\,c_1}\,\sigma^+(t_A)-J_0^R(x_0-c_0\,t_A,\,0)	& \text{if $t_A \geq t_0$},\\
J_0^L(x+c_0\,t,\,0) 											& \text{otherwise},
\end{cases}
\\
\Delta_B(x,\,t)=
\begin{dcases}
-\frac{\sigma^+(t_B)}{\rho_1\,c_1}	& \text{if $t_B \geq t_0$},\\
0 									& \text{otherwise},
\end{dcases}
\end{gather}
and by solving
\begin{equation}\label{ODEM}
z'(t)=-\frac{\rho_1\,c_1}{\mathscr{M}(t)}\parens[\Bigg]{g(t)+\parens[\Bigg]{1+\frac{\rho_0\,c_0}{\rho_1\,c_1}+\frac{\parens[\big]{\mathscr{M}'(t) + \mathscr{Q}_M(t)}}{\rho_1\,c_1}}\,z(t)},
\end{equation}
where
\begin{equation}
z(t)=\sigma^+(t),\qquad g(t)=2\,\rho_0c_0J_0^R\parens{x_0-c_0\,t,\,0}.
\end{equation}
As previously, the particular case where $\mathscr{M}(t)=0$ and $\mathscr{Q}_M(t)\neq 0$ leads to
\begin{equation}
\sigma^+(t)=-\frac{2\,\rho_0c_0J_0^R\parens{x_0-c_0\,t,\,0}}{1+\frac{ \rho_0\,c_0}{\rho_1\,c_1}+\frac{ \mathscr{Q}_M(t)}{ \rho_1\,c_1}},
\end{equation}
which is bounded by the static solutions obtained with the extreme values of $\mathscr{Q}_M(t)$.

\section{Bounded amplification}\label{AppAmpli}

Here we prove Proposition~\ref{PropAmpli}. The assumptions are:
\begin{itemize}
\item homogeneous medium, with $\rho_0\,c_0=\rho_1\,c_1=Z$;
\item source point at $x_s<x_0$, with a bounded source function $S$;
\item $T$-periodic modulation of the compliance and of the corresponding dissipation parameter~\eqref{KM-T}, with $\phi_X(t+T)=\phi_X(t)=\varphi_X(\Omega t)$, $T=2\,\pi/\Omega$, and $X=C,Q_C$;
\item continuity of stress ($\mathscr{M}_0=0$ and $\mathscr{Q}_M=0$); similar conclusions hold with a modulated inertia and continuity of velocity ($\mathscr{C}_0=0$ and $\mathscr{Q}_C=0$).
\end{itemize}

The mean of a function $w(t)$ over $T$ is denoted by
\begin{equation}
\bar{w}=\frac{1}{T}\int_0^Tw(t) \, \dd t.
\end{equation}
We need the following lemmas.

\begin{lemm}\label{LemmaDuhamel}
Consider the Ordinary Differential Equation
\begin{equation}
y'(t)=\alpha(t)\,y+\beta(t),\label{ODE_Duhamel}
\end{equation}
where $\alpha(t)$ is a $T$-periodic function, and $\beta(t)$ is bounded. If $\bar{\alpha}\leq 0$, then $y(t)$ is bounded.
\end{lemm}

\begin{proof}
The integrating factor method yields
\begin{equation}\label{Duhamel}
y(t)=e^{A(t)}\parens[\Bigg]{y_0+\int_0^t e^{-A(s)}\,\beta(s) \, \dd s},\quad \text{with} \ A(t)=\int_0^t\alpha(s) \, \dd s.
\end{equation}
If $\bar{\alpha}\leq 0$, then $e^{A(t)}$ does not grow exponentially and $y(t)$ is bounded.
\end{proof}

\begin{lemm}\label{LemmaG}
Let the function
\begin{equation}\label{FuncG}
\Ng(\Na,\Nb)=1+\frac{\Na}{\Nb}\parens[\big]{\sqrt{1-\Nb^2}-1}
\end{equation}
defined on the domain ${\mathcal{D}}=\braces[\big]{(\Na,\Nb)\in\mathbb{R}^2 \st[\big] \abs{\Na}<1, \: \abs{\Nb}<1, \: \Nb\neq 0}$. Then $\Ng(\Na,\Nb)>0$ for all $(\Na,\Nb)\in {\mathcal{D}}$.
\end{lemm}

\begin{proof}
One introduces
\[
\Nh(\Na,\Nb)=-\frac{\Na}{\Nb}\parens[\big]{1-\sqrt{1-\Nb^2}}.
\]
For $(\Na,\Nb)\in{\mathcal{D}}$, one has
\[
1-\sqrt{1-\Nb^2}\in(0,1),\quad \abs[\Bigg]{\frac{\Na}{\Nb}}<\frac{1}{\abs{\Nb}},
\]
so that
\[
\abs[\big]{\Nh(\Na,\Nb)}<\frac{1}{\abs{\Nb}}\parens[\big]{1-\sqrt{1-\Nb^2} }\coloneqq H(\Nb).
\]
The real function $H$ satisfies
\[
H(\Nb)=\frac{\abs{\Nb}}{2}+o(\Nb)\xrightarrow[\Nb\rightarrow 0]{\hphantom{\Nb\rightarrow \pm 1}}0,\qquad H(\Nb)\xrightarrow[\Nb\rightarrow \pm 1]{}1.
\]
For $\Nb>0$, its derivative writes
\[
H'(\Nb)=\frac{1-\sqrt{1-\Nb^2}}{\Nb^2 \sqrt{1-\Nb^2}}>0,
\]
so that $H$ is increasing from 0 (at $\Nb=0$) to 1 (at $\Nb=1$). Lastly, $H$ is a continuous even function. For all $\Nb\in(-1,1)$, one deduces that
\[
0<H(\Nb)<1.
\]

It follows $\abs[\big]{\Nh(\Na,\Nb)}<1$, and hence
\[
\Ng(\Na,\Nb)=1+\Nh(\Na,\Nb)>1-1=0,
\]
which concludes the proof.
\end{proof}

Following the same lines as in Appendix~\ref{AppAnalytic}, it can be proven that $y(t)=v^+(t)$ satisfies the Ordinary Differential Equation
\begin{equation}\label{ODE-bounded}
y'=a(t)\,y+b(t),
\end{equation}
with
\begin{equation}\label{ABt}
a(t)=-\frac{1}{Z}\,\frac{2+Z\,\parens[\big]{\mathscr{C}'(t)+\mathscr{Q}_C(t)}}{\mathscr{C}(t)},\qquad b(t)=-\frac{1}{Z}\,\frac{S(t)}{\mathscr{C}(t)}.
\end{equation}
The function $a(t)$ is $T$-periodic, and $b(t)$ is bounded. If $\bar{a}\leq 0$, then the conclusion of Lemma~\ref{LemmaDuhamel} holds and $y(t)$ is bounded. Let us consider again the modulation $\varphi_X(t)=\sin(t)$, with $X=C,Q_C$. From~\eqref{ABt}, one obtains
\begin{equation}\label{FuncA}
a(t)=-\frac{1}{Z\,\mathscr{C}_0}\,f(t),
\end{equation}
with
\begin{equation}
f(t)=\frac{2+Z\mathscr{Q}_{C_0}+\gamma\,\cos(\Omega t)+Z\mathscr{Q}_{C_0}\varepsilon_{Q_C}\sin(\Omega t)}{1+\varepsilon_C\,\sin(\Omega t)},\quad \gamma=Z\,\mathscr{C}_0\,\varepsilon_C\,\Omega.
\end{equation}

The mean value of $f$ writes
\[
\begin{split}
\bar{f}
	& = \frac{1}{T}\int_0^T\frac{2+Z\mathscr{Q}_{C_0}+\gamma\,\cos(\Omega t)+Z\mathscr{Q}_{C_0}\varepsilon_{Q_C}\sin(\Omega t)}{1+\varepsilon_C\,\sin(\Omega t)} \, \dd t,
\\	& = \frac{1}{2\,\pi}\parens*{(2+Z\mathscr{Q}_{C_0})\int_0^{2\pi}\frac{\dd s}{1+\varepsilon_C\,\sin(s)}+\gamma\int_0^{2\pi}\frac{\cos(s)}{1+\varepsilon_C\,\sin(s)} \, \dd s\,+\,Z\mathscr{Q}_{C_0}\varepsilon_{Q_C}\int_0^{2\pi}\frac{\sin(s)}{1+\varepsilon_C\,\sin(s)} \, \dd s }.
\end{split}
\]
Classical change of variables, and integration of a periodic function over one period, gives
\begin{equation}
\begin{gathered}
	\int_0^{2\pi}\frac{\dd s}{1+\varepsilon_C\,\sin(s)}=\frac{2\,\pi}{\sqrt{1-\varepsilon_C^2}},\qquad \int_0^{2\pi}\frac{\cos(s)}{1+\varepsilon_C\,\sin(s)} \, \dd s=0, 
\\	\int_0^{2\pi}\frac{\sin(s)}{1+\varepsilon_C\,\sin(s)} \, \dd s=\frac{2\pi}{\varepsilon_C}\parens[\Bigg]{1-\frac{1}{\sqrt{1-\varepsilon_C^2}}}.
\end{gathered}
\end{equation}

Straightforward calculations give
\begin{equation}
\bar{f}=\frac{1}{\sqrt{1-\varepsilon_C^2}}\parens[\big]{2+Z\mathscr{Q}_{C_0}\,\Ng(\varepsilon_C,\varepsilon_{Q_C})},
\end{equation}
where $\Ng(\,\cdot\,,\cdot\,)$ has been introduced in~\eqref{FuncG}. Equation~\eqref{FuncA} and Lemma~\ref{LemmaG} yield
\begin{equation}
\bar{a}=-\frac{1}{Z\,\mathscr{C}_0}\,\frac{2+Z\mathscr{Q}_{C_0}\,\Ng(\varepsilon_C,\varepsilon_{Q_C})}{\sqrt{1-\varepsilon_C^2}}<0.
\end{equation}
Lemma~\ref{LemmaDuhamel} implies that the solution of~\eqref{ODE-bounded} is bounded, which concludes the proof of Proposition~\ref{PropAmpli}.

\section*{Acknowledgments}

The manuscript was written through contributions of all authors.

We would like to thank the reviewers for their useful comments, in particular for the suggestion to study the evolution of energy. Lastly, our thanks go to Vincent Pagneux for the valuable discussions that contributed to this work.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography

\end{document}