A multiple scattering formulation for elastic wave propagation in space-time modulated metamaterials

Space-time modulation of material parameters offers new possibilities for manipulating elastic wave propagation by exploiting time-reversal symmetry breaking. Here we propose and validate a general framework based on the multiple scattering theory to model space-time modulated elastic metamaterials, namely elastic waveguides equipped with modulated resonators. The formulation allows to consider an arbitrary distribution of resonators with a generic space-time modulation profile and compute the wavefield within and outside the resonators' region. Additionally, under appropriate assumptions, the same framework can be exploited to predict the waveguide dispersion relation. We demonstrate the capabilities of our formulation by revisiting the dynamics of two representative space-time modulated systems, e.g. the non-reciprocal propagation of (i) flexural waves along a metabeam and (ii) surface acoustic waves along a metasurface. Given its flexibility, the proposed method can pave the way towards the design of novel devices able to realize unidirectional transport of elastic energy for vibration isolation, signal processing and energy harvesting purposes.


Introduction
In the last decade, the research on active (or activated) materials has fueled the discovery of novel dynamic functionalities to design devices for vibrations and waves control [1,2].Activated materials are often characterized by constitutive properties that are modulated in space and time according to an external energy source.The study of such space-time modulated materials was originally pioneered in optics [3] and, shortly afterward, extended to acoustics [4] and elasticity [1].Elastic waves propagating in these space-time varying media are of particular interest since the modulation can create a directional bias that breaks the time-reversal symmetry.Breaking reciprocity allows to realize rich and unconventional phenomena, including, but not limited to, unidirectional wave propagation, adiabatic energy pumping [5,6], frequency conversion [7].These effects can be leveraged to design novel devices such as acoustic rectifiers [8], circulators [9], and topological insulators [10], which can find applications in acoustic communication, signal processing, energy harvesting and vibration isolation [11,12,13].
In the context of elastodynamics, space-time modulation can be achieved following two strategies.The first one relies on a bias directly introduced in the waveguide, as a modulation of the elastic and/or mass properties, so to obtain a modulated phononic crystal [14,15,16].The second option utilizes space-time modulated mechanical oscillators attached to a non-modulated waveguide [17,18,19] to obtain a modulated elastic metamaterial.Both approaches proved to be technically feasible by a series of experimental works where programmable electric components were used to modulate the media/oscillators [20,21,22,7,23].Nonetheless, modulated metamaterials, compared to their phononic counterpart, are easier to realize, since only the resonant elements need to be modulated, and support non-reciprocal effects at sub-wavelength scales.
Besides the numerous examples of modulated waveguides [24,25], most of the conducted studies rely on the use of numerical simulations, typically developed via finite element (FE) or finite difference (FD) algorithms, to describe the expected non-reciprocal effects.Nonetheless, numerical simulations are always bounded by their computational cost which inherently limits the development of design and optimization studies.Computationally inexpensive analytical tools for modulated media are thus desirable, not only to reduce the computational burden but also to gain a deeper understanding of non-reciprocal effects.Currently, analytical formulations for time-modulated systems are mainly used to predict the dispersion relations of both discrete [26,27,28] and continuous media [15,17,20,21,22,18,19].
Although knowledge of the dispersion relations provides physical insights into the existence of directional band gaps, evidence of non-reciprocal phenomena can be found only by computing transient or steady-state responses across finite modulated systems.To the best of our knowledge, analytical methods for the computation of wavefields and transmission/reflection coefficients are currently limited to one-dimensional (1D) problems [29,30,20,16].Additionally, there is no unified framework that enables the computation of both the dispersion relation and the wavefield of a generic modulated system.
To fill this gap, we here propose a generalized multiple scattering formulation able to model the dynamic response of space-time modulated resonators coupled to a generic elastic waveguide.As observed in experiments, space-time modulated resonators can generate scattered fields at lower and higher harmonics with respect to the excitation frequency [20].To capture this response, we first describe the coupling between the vibrating resonators and the waveguide motion with an ad-hoc impedance operator able to account for the expected additional harmonics.Then, we compute the scattered fields in the waveguide by means of Green's functions.
Finally, we set our multiple scattering scheme to couple the incident and scattered fields and compute the related unknown amplitudes by imposing proper boundary conditions at each resonator base.The proposed formulation allows us to investigate the dynamic of an arbitrary number N of resonators with an arbitrary spatial-temporal modulation profile, since all the space-time varying oscillators can be described individually.
Additionally, by introducing appropriate assumptions, the same formulation can be used to derive the related dispersion equations.
The details of the methodology and its modeling capabilities are discussed in what follows.In particular, in Section 2 we describe the proposed general multiple scattering formulation for the computation of the wavefield and the dispersion relation of waveguides coupled with space-time-modulated resonators.In Section 3, we apply the formulation to model flexural waves in a beam and Rayleigh waves on a substrate, both coupled with an array of modulated surface resonators.For both scenarios we show the capability of the formulation to predict non-reciprocal guided waves.Finally, we derive conclusions and outlook of the work in Section 4.

Statement of the problem
We propose a general analytical framework to model a cluster of space-time-modulated oscillators attached to a given elastic waveguide (Fig. 1).The formulation includes the following three steps: (i) the definition of the elastic force exerted on the waveguide by a time-modulated resonator when excited by a base motion; (ii) the use of Green's functions to describe the scattered wavefields generated by the resonators in the waveguide; (iii) the construction of a multiple scattering formulation to couple the waveguide with an arbitrary number of time-modulated resonators.The approach allows computing the lower-and higher-order scattered harmonics, generated by the collective response of the time-modulated resonators, and responsible for the non-reciprocal wave motion in the waveguide.
First, we present the framework in its most general form, i.e. considering a finite array of time-modulated oscillators with mechanical properties obeying the same modulation period T m and arbitrarily arranged over the waveguide surface.Then, we show how to derive the waveguide dispersion relation by introducing appropriate assumptions, e.g., considering an infinite array of identical resonators regularly arranged along the elastic support.

Incidence Scattering
Fig. 1.Schematics of space-time modulated resonators laying over an elastic waveguide.

Elastic force of a time-modulated resonator
Let us recall the dynamics of the generic nth resonator attached to the waveguide surface at the location r n (see Fig. 1).The resonator has a mass m n , damping coefficient c n , and time-modulated spring stiffness k n (t): where T m is the modulation time period.The governing equation of the nth resonator motion reads: in which W n (t) = W (r n , t) denotes the mass vertical displacement while w n (t) = w(r n , t) is the vertical motion at the resonator base.Accordingly, the point force F n (t) = F (r n , t) exerted by the resonator onto the waveguide surface reads: Since the modulated stiffness k n (t) in Eq. ( 1) is time-periodic, we express it in Fourier series form as: in which i = √ −1 is the imaginary unit, ω m = 2π/T m is the modulation frequency, and where the Fourier coefficients are defined as: As we will see in the next section, the motion along the waveguide excited by a harmonic (e iωt ) incident field, contains several lower-and higher-order harmonics generated by the scattering of the time-modulated mechanical resonators.As a result, the vertical motion at the resonator base, namely the motion at the waveguide surface, can be written as [29]: so that the solution of Eq. ( 2) is sought in the form [21,20,16]: Substituting Eqs. ( 4), ( 6) and (7) into Eq.( 2), yields: Exploiting the orthogonality of harmonic functions, we simplify Eq. ( 8) by multiplying it for ω m e −ipωmt /(2π), and integrating it from −π/ω m to π/ω m , to obtain: By truncating the orders from −P to P , Eq. ( 9) can be reorganized in matrix form as: with: in which m(j) The vertical force at the base of the resonator can thus be obtained by substituting Eq. ( 7) into Eq.( 3): where the F (h) n coefficients from h = −P to h = P , collected in the vector Fn , read: with: In Eq. ( 13), the matrix Z n is the impedance operator which relates the resonator base motion to the resonator base force.It can be observed that the force exerted by each modulated resonator on the elastic substrate comprises multiple harmonics (ω + hω m ).In the next section, we discuss how these forces generate the related multiple scattered wavefields.

Elastic wave field of a finite cluster of modulated resonators
We now consider an arbitrary distribution of N space-time modulated resonators arranged on top of a given elastic waveguide.We assume that the resonators have an identical stiffness modulation frequency ω m .
When an incident wave field u 0 = [u 0 , v 0 , w 0 ] impinges the bases of such resonators, it triggers their vibrations which, in turn, generate scattered waves in the waveguide.Following a standard multiple scattering description [31,32,33], the total wave field u = (u, v, w) at the generic position r along the waveguide can be expressed as the summation of the incident and scattered wave fields of the N resonators: w(r, t) = w 0 (r, t) where G u , G v , G w are the related Green's functions in terms of displacements along x, y, z.As in the previous section, we express the displacements of Eqs.(15a), (15b), (15c) accounting for the multiple harmonics: Truncating the harmonic terms from h = −P to h = P , Eqs. (16a), (16b), (16c) can be rewritten as: with: and where φ0 has non zero components only for the incident field ϕ 0 = u 0 , v 0 , w 0 : Note that in Eqs.(17a), (17b), (17c) the total displacement components û, v, ŵ and the elastic force coefficients Fn are unknown.Nonetheless, following a classical multiple scattering approach, we can obtain the coefficients Fn by setting boundary conditions at resonator bases.In particular, we substitute Eq. ( 13) into Eq. (17c) and specify it at the resonator location r m : Eq. ( 19) leads to a system of m = N equations that we can recast in matrix form as: with: It follows that for a given incident wave field ŵ0 , the vector X of the force amplitudes Fn can be computed as X = A −1 B, and the total wave field in the waveguide evaluated by using Eqs.(15a), (15b), (15c).
In addition, we will show in the following subsection that Eq. ( 20), under appropriate assumptions, allows to derive the dispersion relation of time-modulated waveguides.

Dispersion relation
We here consider an infinite array of equally spaced resonators, arranged atop an elastic waveguide (see Fig. 2) at lattice distance a.We restrict our interest to oscillators with identical mass and with spring constant modulated in time and space with a wave-like modulation of period T m = 2π/ω m and wavelength λ m = 2π/κ m , whose general form reads: As before, we express k(x, t) in a Fourier series form: where the Fourier coefficients are computed as: As discussed in [29,19], a stable response of the modulated system requires each modulation amplitude k(j) (j = 0) to be sufficiently small with respect to the static stiffness k(0) .Under this assumption, for the assumed infinite (N → ∞) periodic array of identical resonators, the scattering Eqs.(19) are the same at any location x m and satisfy: where x m has been conveniently set equal to 0. Following the effective medium approach [34,35,36], namely considering the lattice spacing a much smaller than the characteristic wavelength, we approximate the discrete point force as an average line load.As a result, the total vertical displacement at the generic resonator base can be computed as: Due to the space-time modulation of the resonator properties, we can express the force in Eq. ( 12) in the form: Substituting Eq. ( 27) into Eq.( 26) and truncating the orders from h = −P to h = P we obtain: F (−P +1) . . .
Some minor algebra yields the following system of homogeneous equations: . . .
in which G(κ + P κ m , ω + P ω m ) is P th order Green's function in space-domain which is obtained as: Non-trivial solutions of Eq. ( 29) provide the dispersion relation:

Examples and applications
To show the potential of our formulation, we consider two space-time-modulated waveguides that have been thoroughly discussed in previous works [17,20,18,19].We begin our investigation by considering an Euler beam coupled with an array of modulated resonators.For this example, we validate our approach against the results of Transfer Matrix Method (TMM).For the sake of completeness we report in Appendix A the full derivation of the adopted TMM [20].As a second example, we consider a 2D elastic half-space coupled to modulated resonators.For this configuration, given the absence of closed-form formulations, we compare our findings vs. those obtained via finite element simulations, as in Ref. [19].

Modeling non-reciprocal flexural waves in a space-time modulated beam
We consider an Euler-Bernoulli beam equipped with an array of undamped resonators, see Fig. 2a, modulated in a wave-like fashion according to the relationship [27,17,20,26]: where k 0 denotes the static stiffness, k a the amplitude of the modulation, ω m the modulation angular frequency, κ m the modulation wavenumber.At any location x n , the modulated stiffness is time-periodic and its non-zero Fourier coefficients read: as k(j) n = 0 for |j| > 1.For the numerical example, we consider the mechanical parameters originally adopted in Ref. [17], i.e., a resonator mass m 0 = ρAa, where ρ is the mass density of the beam and A is the cross-section area of the beam.Similarly, the modulation frequency is set as ω m = 0.25ω 0 and modulation amplitude as k a = 0.2k 0 , in which ω 0 is the resonance frequency of resonators and k 0 = m 0 ω 2 0 is the unmodulated stiffness; the modulation wavenumber is κ m = 1.25κ 0 , where κ 0 = 4 k 0 /(aD), D the bending stiffness of the beam.

Dispersion relation
According to the Euler-Bernoulli beam theory, the P th order governing equation under the action of a harmonic point force can be written as: We Fourier transform Eq. ( 34) along the x direction, and obtain the transformed P th order Green's function in space-domain as: where the shifted frequency and wavenumber are defined as: First, by setting P = 0 we get the non-modulated impedance parameter Z from Eq. ( 13) as: Substituting Eqs. ( 35) and (37) into Eq.( 31) we obtain immediately the dispersion relation of a non-modulated metabeam: This dispersion equation is identical to the one obtained in Refs.[17,20] and matches the dispersion curve provided by FE simulations, see Appendix B for details.
In the presence of modulation, scattered waves are expected when the phase matching condition is satisfied, i.e., C(κ, ω) = C(κ P , ω P ) [27].As an example, in Fig. 3a we show the original (C) and the two shifted dispersion curves (C ±1 ) for P = ±1, respectively.The phase matching condition is met at the intersections between the original curve and the shifted ones, namely at six magenta points of the pairs (A), (B) and (C).The asymmetric distribution of these intersections suggests the breaking of time-reversal symmetry which, in turn, leads to direction-dependent phenomena within the metabeam [20].
We now predict the dispersion properties of the modulated metabeam.To do so, we substitute Eq. (35) into Eq.( 31) by truncating waves to the first order (P = 1), which yields: with the impedance operator: We remark that the coupled dispersion relation in Eq. (39) holds only near the above-mentioned intersections in Fig. 3a [28].Thus, we compute and plot the coupled dispersion in the range of ±0.1κ and ±0.1ω around each crossing point, as shown in Fig. 3b (red circular markers).For comparison, we also provide the unmodulated dispersion curve (solution of Eq. ( 38)) and its shifted analogs on the same figure.As discussed in Ref. [17], in the vicinity of pair B no directional band gap is generated, since both modes have positive group velocities.
Conversely, for contra-directional branches such as pairs A and C, the repulsion effect can lead to narrow band directional gaps, for instance, around A. Within these gaps, waves are hindered only when propagating along the specific direction (dictated by the sign of the related wavenumber); conversely, they are fully transmitted when propagating along the opposite direction [27].This directional wave-filtering is usually accompanied by the generation of lower/higher-order waves at the phase-matched frequencies, thus resulting in a reflection combined with a frequency conversion [17].Evidence of these effects is provided in the next section where the steady-state solution of waves propagating along a finite modulated metabeam is computed.

Steady-state solutions
To evidence the non-reciprocal behavior predicted by the dispersion analysis, we utilize Eq. (17c) to compute the steady-state response of a finite metabeam.In particular, we are interested in verifying the non-reciprocal reflection/transmission in the directional band gap at pair A in Fig. 3.As an example, an array of 50 resonators is considered for these investigations.The response is recorded at locations x r and x t , and later used to compute the reflection and transmission values, respectively.In both scenarios, the harmonic point source e iωt and the receiver are located at distances of d s = 600a and d r = 300a from the closest oscillator.
According to the formulation discussed in Section 2, the impedance operators Z 1 to Z N are obtained from Eq. ( 13) while the P th order Green's function in Eq. ( 20) is obtained by applying the inverse Fourier transform to Eq. ( 35) as: where the P th order wavenumber for flexural waves reads: Substituting Eq. (41) into Eq.( 20) we obtain the elastic force coefficients Fn , which are inserted into Eq.
(17c) for the calculation of the displacement components ŵ(x) in the beam.
We begin our investigation considering a right-propagating (κ > 0) flexural wave at frequency ω = 1.66ω 0 , i.e., the intersection at pair A in Fig. 3a.The reflection coefficient, normalized with respect to the incident wave, |w r /w 0 |, is displayed in Fig. 3c, considering the scattered waves truncated at P = ±5 order.As expected, right-propagating incident waves at ω = 1.66ω 0 (dashed line in Fig. 3c) undergo a strong reflection with different frequency contents, including the largest component at the first-order harmonic (1.66ω 0 − ω m ) and non-negligible components at the second (1.66ω 0 − 2ω m ) and third-order harmonic (1.66ω 0 − 3ω m ).
The amplitude of other higher-order harmonics is negligible.Conversely, the left-propagating wave (opposite to the modulation direction) with the same frequency ω = 1.66ω 0 can travel through the resonators almost undisturbed, as shown by the normalized transmission |w t /w 0 | in Fig. 3d.To verify the predictions provided by our approach, we compute the same transmission and reflection coefficients using the transfer matrix method.
The results, marked by solid lines in Figs.3c,d, are in excellent agreement with our analytical solutions (see more details on the transfer matrix method in Appendix A).

Modeling non-reciprocal Rayleigh wave propagation in a space-time modulated metasurface
We now consider the propagation of Rayleigh waves across a cluster of modulated resonators.Such a problem has been recently investigated with the aid of FE numerical simulations [18,19].Our purpose is to show the capability of the proposed analytical formulation to reproduce both the non-reciprocal dispersion and the reflection/transmission coefficients in this complex configuration.

Dispersion relation
Let us briefly recall the Green's function for a 2D isotropic elastic half-space actuated by a harmonic vertical load acting at the surface.For this configuration, the equilibrium equation can be formulated as a boundary value problem: in which c L and c T denote the longitudinal (L) and transverse (T) wave velocities, and τ zx , σ z represent the shear and normal stresses, respectively; u is the displacement field with components u and w; δ(x) is the Dirac delta function.
In analogy with the metabeam problem, we Fourier transform the equilibrium Eqs.(43a) and (43b) along the x direction, and obtain the transformed P th order Green's function at z = 0 as: where ρ is the density of the substrate, and the shifted frequency ω P and wavenumber κ P are defined in Eq. (36).Substituting Eqs. ( 37) and (44) into Eq.( 31) and setting P = 0, we obtain immediately the dispersion relation for Rayleigh waves existing in a non-modulated metasurface: This dispersion equation is identical to the one obtained in Refs.[35,36], and matches the numerical dispersion curve computed via FEM, see Appendix B.
As for the metabeam scenario, we first plot the unmodulated C(κ, ω) and the shifted C(κ P , ω P ) dispersion curves for P = ±1, Fig. 4a.Again, phase matching points (e.g., pairs A to E) are found when C(κ, ω) = C(κ P , ω P ) is met.We predict the dispersion properties of the modulated metasurface around these points using Eq.(31).To this purpose, we substitute Eq. (44) into Eq.( 31) and truncate the expansion to the first order, using the impedance operator Z computed according to Eq. (38).
We display the modulated dispersion relation in the range of ±0.1κ and ±0.1ω around each intersection in Fig. 4b (red circular markers).As an example, the intersection between contra-directional branches gives rise to the locking pair C which results in a directional band gap.Harmonic waves propagating with wavenumberfrequency falling within the directional gap (1.21κ r , 1.185ω 0 ) are reflected by the metasurface as a propagating mode at the phase-matched frequency-wavenumber pair (1.21κ r −κ m , 1.185ω 0 −ω m ).Conversely, such reflection by conversion does not occur for waves propagating along the opposite direction at the same frequency 1.185ω 0 , confirming the non-reciprocity due to the broken time-reversal symmetry [19].Again, clear evidence of these effects predicted by the dispersion curve is provided in the next section by computing the steady-state solutions of Rayleigh waves propagating along a finite modulated metasurface.

Steady-state solutions
To shown the non-reciprocal Rayleigh wave propagation discussed above, we use Eq.(17c) to compute the steady-state response of a finite metasurface.The impedance operators Z 1 to Z N are computed from Eq. ( 13), while the Green's function in Eq. ( 20) is obtained by applying the inverse Fourier transform to Eq. (44), yielding the P th order wave field (Green's function) at z = 0: Ĝ(P ) w (x, 0, We note that unlike the Green's function of an Euler beam, Eq. ( 41), the one for a 2D elastic substrate is divergent at the origin, Eq. (46).To avoid any convergence issue, we introduce a small footprint of length s for each resonator so that the associated Green's functions are [33]: for the vertical displacement components and for the horizontal ones, where: Substituting Eq. (47a) into Eq.( 20) we obtain the elastic force coefficients Fn , which are used in Eqs.
For our example, we compute the steady-state response at locations (x r , 0) and (x t , 0) on the substrate surface considering an array of 100 resonators with footprint width s = a/20, where the harmonic point source and the receiver are located at distances of d s = 600a and d r = 300a from the closest resonator.For a right-going (κ > 0) incident Rayleigh wave (dashed line) at ω = 1.185ω 0 the reflected field, shown in Fig. 4c, confirms a back-scattering at the coupled frequency ω = 1.185ω 0 −ω m .Conversely, the left-going (κ < 0) incident Rayleigh wave (dashed line) at the same frequency ω = 1.185ω 0 propagates without reflection or frequency conversion phenomena (see Fig. 4d).
We now resort to the FEM to verify our analytical solutions.To this purpose, we build a 2D plane-strain model in a commercial FE software (COMSOL Multiphysics).The reader can find the details of the numerical model in Appendix C. Specifically, we compute the transient response of the system actuated by a vertical tone-burst-shaped force having central frequency ω = 1.185ω 0 , and analyze the vertical displacement field w (see the insets of Fig. 4).The corresponding frequency spectra (solid line) computed through the Fourier transform (FFT) of the record time-domain data at the receiver are displayed in Figs.4c,d.The reader can appreciate how the numerical results match the analytical solutions.
Finally, we inspect the steady-state response at the "veering pair" (intersection between two co-directional branches) [19] E in Fig. 4, where we expect the Rayleigh wave to be transmitted and converted from one harmonic to another [18,19].To evidence such a conversion, we utilize the same model (see Fig. 5) excited by a right-going incident Rayleigh wave at frequency ω = 0.734ω 0 .The results obtained with both the analytical solutions and the FE model are collected in Fig. 5b.As expected, the modulated metasurface can convert the incident wave (ω = 0.734ω 0 ) into a transmitted wave with a different frequency content, e.g., the phase matched first-order harmonic at ω = 0.734ω 0 + ω m .
Both scattering fields exhibit a clear asymmetry, with the right-hand side having a greater amplitude than the left-hand side, a clear feature of the forward scattering behavior at veering pairs of the modulated metasurface.4).The total wave field, free field, fundamental scattered field (ω = 0.734ω 0 ), and the first-order scattered field (ω = 0.734ω 0 + ωm) excited by the source at ω = 0.734ω 0 are shown in (c), (d), (e) and (f), respectively.
Appendix A. Details on the transfer matrix method (TMM) In this Appendix, we provide the details of the transfer matrix method for the modulated beam (see Fig. A.1) [20].According to Euler beam theory, the pth order displacement in the nth cell can be expressed as: in which β (p) = 4 ρA(ω + pω m ) 2 /D is the pth order wavenumber.By truncating the orders from p = −P to p = P , the displacement can be expressed in matrix form ŵn (x) = L n (x)A n , with: Hence, the vertical force Fn in Eq. ( 13) can be written as: For an arbitrary pth order, continuities of the displacement, slope, moment and shear force at . . .

Fig. 3 .
Fig.3.(a) Dispersion curve of a non-modulated metabeam (black dashed lines) and its shifted analogs for P = −1 and P = 1, respectively.(b) Dispersion curves (circular red markers) of a modulated metabeam in proximity of the phase matching pairs.Normalized reflection and transmission coefficients for flexural waves propagating at ω = 1.66ω 0 inside the directional band gap (pair A) for (c) a right and (d) a left traveling incident wave (star marker), respectively.For comparison, results obtained by the transfer matrix method (TMM) are also provided (blue solid lines).

Fig. 4 . 5 .
Fig. 4. Analytical modeling of a metasurface.(a) Dispersion relation of a non-modulated metasurface (black solid lines) and its shifted analogs for P = −1 and P = 1, respectively.(b) Dispersion relation (circular markers) of a modulated metasurface in the vicinity of phase matching pairs.Steady-state solutions of Rayleigh wave propagation at ω = 1.185ω 0 inside the narrow directional band gap (pair C) for (c) a right and (d) a left traveling incident wave (star marker), respectively.

Fig. A. 1 .
Fig. A.1.Schematic of transfer matrix method: (a) the nth cell, (b) the global system.

22 M
A.14) can be further transformed to:ψ out = Sψ in , 11 − M 12 M −122 M 21 M 12 M (A.17) we can compute both the transmission and reflection coefficients directly.It is worth mentioning that, due to the presence of exponential amplification terms in α(p)n−1 in Eq. (A.5), the transfer matrix method may encounter numerical divergence in some occasions, e.g., when considering a large number of oscillators or large values of resonators spacing.Such a limitation can be well addressed by the multiple scattering formulation proposed in this work.

Fig. C. 1 .
Fig. C.1.Schematic of the FE model for transient simulations.