Nonlocal effective medium theory for phononic temporal metamaterials

We have developed a nonlocal effective medium theory (EMT) for phononic temporal metamaterials using the multiscale technique. Our EMT yields closed-form expressions for effective constitutive parameters and reveals these materials as reciprocal media with symmetric band dispersion. Even without spatial symmetry breaking, nonzero Willis coupling coefficients emerge with time modulation and broken time-reversal symmetry, when the nonlocal effect is taken into account. Compared to the local EMT, our nonlocal version is more accurate for calculating the bulk band at high wavenumbers and essential for understanding nonlocal effects at temporal boundaries. This nonlocal EMT can be a valuable tool for studying and designing phononic temporal metamaterials beyond the long-wavelength limit.


Introduction
The bulk composite materials with meta-atoms spatially arranged in the subwavelength scale, which are well-known as metamaterials, possess the ability to manipulate classical waves in various novel ways due to their abnormal macroscopic properties [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] and thereby have drawn great attention in the past few decades.Based on the idea of effective medium theory (EMT), the macroscopic properties of metamaterials can be well described by a set of effective constitutive parameters when a homogenization process is imposed, which simplifies the theoretical and numerical study greatly.Although the EMT was initially applied to the regime where the long-wavelength limit holds, i.e. both the frequency (ω) and wavenumber (k) are approaching zero, more and more Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Recently, there has been a surge of interest in the timevarying media the material properties of which are modulated periodically in time [44][45][46][47].In addition to a new degree of freedom to tailoring the material properties, the time-modulation also breaks time-reversal symmetry and induces temporal discrete translational symmetry of the system when the modulation frequency is higher than or comparable with the operating frequency.In phononics, these time-modulations can be achieved much easier due to the low acoustic and elastic wave frequencies, and they have been used recently for realizing Floquet topological states [48,49], nonreciprocal wave transportations [50][51][52][53][54], frequency conversion [55,56] and selective wave filtering [57].
For the modulation frequency much higher than the operating frequency, the time-varying media can serve as the temporal analog of metamaterials, which are called the temporal metamaterials.Besides other various investigations [58][59][60][61], the EMT concept has also been applied to the temporal metamaterials recently [50,[62][63][64][65][66].The EMT with closed-from expressions have been introduced to the photonic temporal multi-stepped metamaterials in the long-wavelength limit [61] and later verified in acoustic wave experiment [63].Moreover, homogenization theories have also been applied to phononic and photonic space-time metamaterials with their constitutive parameters modulated in a travel wave form [51,64].These works obtain effective Willis coupling and effective bi-anisotropic coupling tensors which are directly related to the breaking of reciprocity.Recently, nonlocal EMTs have been proposed to study the spatial dispersions [65] and nonlocal effects at the temporal boundaries [66] in photonic temporal metamaterials.
In this study, we harnessed the general multiscale technique [37,50,67] to establish a nonlocal EMT for phononic temporal metamaterials.Our novel framework offers analytical expressions for effective constitutive parameters, applicable to materials with spatially homogeneous but periodically time-modulated stiffness and mass density.Remarkably, we unveiled that the emergence of Willis coupling is contingent upon the simultaneous temporal modulation of stiffness and mass density.A crucial condition for Willis coupling is the disruption of time-reversal symmetry, distinguishing it from traditional Willis metamaterials where spatial inversion symmetry breaking is the norm [32][33][34][35][36]68].The nonlocal origin of Willis coupling makes it beyond the predictive capability of conventional local EMT.Our expressions elucidate that effective stiffness and mass density are real, expressed in even power series of the wavenumber, while Willis coupling coefficients are purely imaginary and linear to the wavenumber.Consequently, despite the presence of nonzero Willis coupling, the band dispersion of phononic temporal metamaterials remains symmetric, a marked departure from space-time modulated Willis metamaterials that often exhibit strong asymmetry [50,51].In comparison to the local EMT, our nonlocal EMT excels not only in generating more precise bulk band dispersions but also in accurately forecasting nonzero reflections at impedance-matched temporal boundaries-an aspect that the local EMT cannot address.It is crucial to underscore that our theory is inherently applicable to phononic temporal metamaterials featuring diverse forms of periodic time modulation, and paves the way for theoretical exploration of temporal phononic metamaterials characterized by strong spatial dispersions and nonlocal effects.

Analysis based on multiscale technique
Let us consider that the stiffness κ and mass density ρ of an isotropic phononic medium are modulated periodically in time with the same time-modulation period t 0 , namely κ (t + t 0 ) = κ (t) , ρ (t + t 0 ) = ρ (t).Here, we assume that the medium is dispersiveless.Based on this assumption, κ (t) and ρ (t) are real and positive valued for any time.The phononic medium we considered is spatially homogenous.While a plane wave propagates inside the medium, the wavenumber k must be conserved due to the homogeneity.Suppose the wave propagates along the x direction, the time-dependent motion equations are expressed as where v (x, t) , σ (x, t) are the particle velocity and stress, ε (x, t) , µ (x, t) are strain and momentum density, and according to the constitutive relations they are related as For the sake of convenience, hereafter, we introduce the parameter ζ = ρ −1 as the inverse of mass density.Then the two equations in equation ( 2) are in the similar forms.We stress that, for dispersive medium, instead of the simple linear relations as shown in equation ( 2), the constitutive relations should be σ , where * denotes the convolution operator.
Here, we define a new parameter K = Ω/c 0 , where Ω = 2π /t 0 is the modulation frequency, c 0 is the phase velocity of acoustic wave in the static case, to characterize the modulation speed in the momentum space.c 0 can be chosen as √ κζ , where κ, ζ are the temporal averages of κ (t) and ζ (t) over a period t 0 , respectively.When we focus on the region where the wavenumber k is smaller than K, it is convenient to introduce the small parameter η = k K < 1 and the fast temporal coordinate T = t/η.Based on the idea of multiscale technique [62], the fast coordinate T can be regarded as an independent variable from the slow temporal coordinate t.Then, following the idea of [37], any time-dependent field can be expanded in Maclaurin series of η, and decomposed into the temporal averaged parts and the corresponding zero mean residuals, namely where A = σ, ε, v, p, κ, ζ, the temporal averages Ān (x, t) and zero mean residuals Ãn (x, t, T) are defined as We note that the temporal average of any function in the following is based on the first of equation ( 4).Since the constitutive parameters are time-periodic functions with the period t 0 , they can be expanded in Fourier series with e inΩηT .Therefore, they are functions of the fast coordinate T only: With the additional temporal coordinate T, the time partial derivative becomes As such, substituting equations (3) and ( 5) into (1), extracting each order of η and separating the temporal averaged fields and their zero mean residuals [37], we obtain for the temporal averaged fields and for the zero mean residuals.In equation ( 6), ζ μn and κε n represent the temporal averages of ζ μn and κε n , respectively.Multiplying η n to equation (6) and summing over all the terms, we obtain the motion of equations in the macroscopic scale as where with ∆v n = ζ μn , ∆σ n = κε n .We can first solve equation (7) for n = 0, and then obtain solutions of higher orders iteratively, see details in appendix A.

Willis coupling
According to equation (10), the constitutive relations for the macroscopic fields can be written as where with From equation ( 15), we see that the phononic temporal metamaterial exhibits the Willis coupling when k ̸ = 0, ψ ̸ = 0.As a consequence, there is no Willis coupling in the limit k → 0. Because the time-dependent stiffness and mass density considered here are always real, the expansion coefficients satisfy the relations κ where * denotes the complex conjugate.Then, according to (equation ( 13)), where Im (•) denotes the imaginary part.Therefore, ψ is purely real, making that the Willis coupling coefficient χ e = ikψ = χ * e is purely imaginary.This imaginary Willis coupling coefficients satisfy χ * e (k) = χ e (−k) , χ * e (k) = χ e (−k) which meets the requirement by reciprocity in Fourier space [36].According to equation (18), ψ vanishes when either κ n̸ =0 = 0 or ζ n̸ =0 = 0, namely there is no Wills coupling when a single constitutive parameter is temporally modulated.However, if both the stiffness and mass density are temporally modulated, ψ is zero only when Im (κ * n ζ n ) = 0 for all n, which is not easily fulfilled if no other constrains are implied.
If time-reversal symmetry is preserved, there must be temporal inversion centers Therefore, there will be no Willis coupling when timereversal symmetry is preserved for the phononic temporal metamaterials.
In the following, we consider that κ (t) and ζ (t) are even functions about t c1 = 2Nπ /Ω − ϕ 1 /Ω and t c2 = 2Nπ /Ω − ϕ 2 /Ω, respectively, where N is an arbitrary integer, ϕ 1 , ϕ 2 are the initial phases of time modulation functions for the stiffness and mass density, which of course depend on the origin of temporal coordinate chose.Then, for complex valued κ n and ζ n , they can be expressed as Inserting equation (31) into equation ( 28), we obtain Therefore, the Willis coupling coefficient is nonzero when time-reversal symmetry is broken since there will be no temporal inversion center t c that κ Therefore, the nonzero Willis coupling in this case just originates from the breaking of time-reversal symmetry, which is in sharp contrast to conventional Willis metamaterials [27][28][29][30][31]63] where the Willis coupling stems from the breaking of spatial inversion symmetry.However, it should be noted that the breaking of time-reversal symmetry is a necessary but not the sufficient condition for achieving nonzero Willis couplings.For example, when the stiffness κ is static (κ n̸ =0 = 0) and the time-dependent mass density ρ (t) is not symmetric about any temporal coordinate, time-reversal symmetry is broken but the Willis coupling coefficient ikψ is zero.

Spatial dispersion and reciprocity
Besides the Willis coupling, the nonlocal EMT based on equation ( 15) also takes the spatial dispersion into account, with the strengths characterized by the parameters 16) and χ * e (k) = χ e (−k), following the constrains of Fourier space reciprocity [36] and leading to ω (k) = ω (−k), see equation (34) in the following.Therefore, in despite of the nonzero Willis coupling, the phononic temporal metamaterials are reciprocal exhibiting symmetry band structure, which is different from space-time modulated materials [50,51].This is because asymmetric band structure requires the simultaneous breaking of inversion and time-reversal symmetry [69][70][71], while inversion symmetry of the spatial homogeneous phononic temporal metamaterial is preserved.As a consequence, no matter which order of η we truncate to in the analysis, the reciprocity guarantees that the effective constitutive parameters κ e , ζ e must be in even powers series of k.For example, if we truncate to the fourth order of η, the first and third order terms vanish in κ e , ζ e , ) , (23) where γ ′ 1 , γ ′ 2 are constants to be determined.When the stiffness and mass density are temporally modulated in phase (ϕ 1 = ϕ 2 ), γ 1 , γ 2 will be independent of the initial phase ϕ 1 , since the initial phases are eliminated in the terms This difference is due to the fact that the media with different ∆ϕ are corresponding to different phononic time-Floquet crystals.
In the small wavenumber limit k → 0, the spatial dispersion can be ignored and then the effective constitutive parameters are reduced to the well-known results based on the local effective medium description [62]

Bulk band dispersions
Substituting equation (15) into equation ( 1), the motion equation is rewritten as From equation (26), we obtain the dispersion relations of the phononic temporal metamaterial as and the corresponding plane wave modes as According to the group velocity ∂ω/∂k, the solution with plus (minus) sign in equation ( 27) is for the wave propagates forward (backward).In the following, we will consider two typical cases to show that the nonlocal EMT can well predict the bulk band dispersions even for large wavenumbers.

Temporal stratified media
The first typical case of phononic temporal metamaterial is the temporal stratified medium with the time-dependent stiffness and mass density in the A and B layers given by where N is an integer, t 1 < t 0 denotes the temporal thickness of layer A, (κ A , ρ A ) and (κ B , ρ B ) are the constitutive parameters of layer A and layer B, respectively, which are all real and positive.Based on equation ( 29), the stiffness and mass density are temporally modulated in phase, namely ϕ 1 = ϕ 2 = 0, and both κ (t) and ρ (t) are symmetric about the temporal centers t c = 2Nπ /Ω.Therefore, there is no Willis coupling in this kind of temporal metamaterial.The Fourier coefficients are then calculated as where When the temporal coordinate origin is shifted by δt, all the Fourier coefficients attain a phase term as: However, γ 1 , γ 2 are kept the same as we have discussed previously, making the effective constitutive parameters independent of the temporal origin.Substituting equation (30) into equations ( 13) and ( 14), we obtain where with B 4 (x) being the fourth order Bernoulli polynomial.
Because ϑ 2 ≪ 1, it is safe and convenient to ignore the term related to it.By doing so, the parameters γ 1 , γ 2 are approximated as , (33a) In figure 1(b), we showed the dispersion of temporal stratified medium calculated by different means.The dispersion can be rigorously obtained by using the temporal transfer matrix method (TMM), see details in appendix B, with the result shown by the black line.Since the band dispersion is symmetric about k = 0 as we have discussed previously, and the EMT can only be applied to lowest band, we only show the band of order 1 with k ⩾ 0 in the first Brillouin zone here.In the framework of EMT, the dispersion can be also approximated by equation (27).We use equation (25) based on the local EMT and equation (16) based on the nonlocal EMT to calculate the effective constitutive parameters, respectively.
The EMT totally fails inside the k-gap region, see figure 1(b).This is because the EMT works only for the single valued dispersion relation ω(k).However, inside the kgap, because two bands coalescence in their real parts, each wavenumber k corresponds to a complex conjugate frequency pair [47].Therefore, the nonlocal EMT is valid within the range −k edg ∼ k edg , where k edg is the lower edge of the k-gap.
In figure 1(a), the parameters γ 1 , γ 2 as functions of the ratio ρ B /ρ A for different κ B /κ A are plotted.It is shown that γ 1 , γ 2 approach zero when the two temporal layers are impedance matched (κ A ρ A = κ B ρ B ).To the contrary, γ 1 , γ 2 can be large when the impedance mismatch between the two temporal layers is significant.However, the larger the impedance mismatch is, the smaller k edg will be.
We calculated the dispersions based on both local and nonlocal EMTs, represented by the red and blue solid lines, respectively.Notably, the results derived from the nonlocal EMT align more closely with the rigorous results, especially as the wavenumber increases.However, it is worth noting that the disparity between the two results may not appear substantial.This is primarily because the parameters γ 1 and γ 2 employed here are small, limiting the strength of spatial dispersion at lower wavenumbers.Efforts to enhance the mismatch between the two temporal layers and subsequently increase γ 1 and γ 2 are constrained by the reduction of k edg .Consequently, the spatial dispersion remains weak due to the small wavenumber.On the other hand, it is essential to recognize that the results based on the nonlocal EMT may deviate from the rigorous results near the k-gap.Fortunately, we can significantly enhance the accuracy of band dispersion calculations by truncating the effective constitutive parameters to much higher orders.As demonstrated by the blue dashed line in figure 1(b), the results are notably improved when considering up to the fourth order (as detailed in equation ( 23)), and the discrepancy between the local and nonlocal EMTs is enlarged.In principle, the nonlocal EMT can accurately predict the entire band by incorporating a sufficient number of orders, in contrast to the local EMT which is only valid for very small wavenumber.
To quantitatively study the accuracy of the local and nonlocal EMTs, we plotted the relative errors of the calculated band dispersions by different means as functions of the wavenumber, see figure 2. These plots reveal that relative errors increase across all means as the wavenumber grows, with a particularly sharp increase near the k-gap.Due to the neglect of spatial dispersion, the local EMT consistently yields larger errors.In contrast, the accuracy of the nonlocal EMT depends on the order of truncation in the multiscale analysis.When truncated to the second order, the band can be calculated with very high accuracy (relative error less than 2%) for k/K < 0.3.However, when truncated to the fourth order, the range of accurate calculation (less than 2% relative error) can be extended to k/K < 0.35, which is very close to the edge of the k-gap (k edge /K = 0.385).Thus, the nonlocal EMT truncated to the fourth order proves to be accurate enough for nearly the entire band.

Sinusoidal time-modulation
Another typical case we consider here is that the stiffness and the inverse of mass density vary harmonically in time as where ∆ 1 , ∆ 2 are small dimensionless parameters, κ 0 , ρ 0 = ζ −1 0 are the stiffness and mass density in the static.Based equation (34), κ (t) and ρ (t) are symmetric about the temporal inversion centers t c1 = 2Nπ /Ω − ϕ 1 /Ω and t c2 = 2Nπ /Ω − ϕ 2 /Ω, respectively.As we have discussed previously, the Willis coupling is nonzero when ϕ 1 ̸ = ϕ 2 .Because the stiffness and mass density must be positive and real, ∆ 1 , ∆ 2 are both real and their magnitudes are less than 0.5.Then the Fourier series of the two parameters are where δ ij is the Kronecker delta function.Substituting equation (35) into equation ( 18), the parameter ψ is And combining (equations ( 13)-( 17), (35a) and ( 35b)), we obtain According to equations ( 36) and (37), ψ , γ 1 , γ 2 are all periodic functions of ∆ϕ .In figure 3(a), we plotted the effective constitutive parameters as functions of the initial phase difference ∆ϕ when ∆ 1 = 0.2, ∆ 2 = 0.3.It is shown that γ 1 , γ 2 reach maximum when ∆ϕ = (2N + 1) π , and reach minimum when ∆ϕ = 2Nπ .For both cases, the Willis coupling coefficient ψ vanishes.Figures 3(b)-(d) depict the dispersion relations of the phononic temporal metamaterials for different initial phase differences.The dispersion relations are calculated rigidly by using the plane wave expansion (PWE) method, see details in appendix C.And similarly, the dispersion relations are also calculated based on equation (27), where the effective constitutive parameters are obtained according to either the local EMT or nonlocal EMT.We can see that the results based on the nonlocal EMT (blue solid lines) are always more accurate than the results based on the local EMT (red dashed lines).For ∆ϕ = π because γ 1 , γ 2 are large, the difference between results based on the two EMTs are remarkable, see figure 3(d).However, for ∆ϕ = 0, there is no Willis coupling and the spatial dispersion is very weak due to the small γ 1 , γ 2 , leading to a vanishing small difference between the two results, see figure 3(b).Nevertheless, both the local and nonlocal EMTs can produce very accurate result in this case.

Temporal boundary
In this section, we will discuss the nonlocal effect induced by the temporal boundary.Consider the stiffness and mass density are given by where κ s , ρ s are time-independent stiffness and mass density, H (t) is the Heaviside step function.Therefore, there is a temporal boundary at t = 0, across which the medium varies from static to dynamic.Assume that the momentum density field before and after the temporal boundary are expressed as where k is the wavenumber, µ i is the complex amplitude of momentum field before t = 0, μt , μr are the complex amplitudes of averaged momentum density fields of the forward and backward propagating waves after t = 0, ω s (k) = √ κs ρs k and ω (k) are the frequencies before and after the temporal boundary.In the framework of nonlocal EMT, ω (k) is given by equation (27).To specify the forward propagating wave, ω (k) takes the positive solution in equation (27).In analogy to the spatial counterparts, µ i , μt , μr can be regarded as the complex amplitudes of incident, transmitted and reflected waves for the temporal boundary.
According to the motion equations, the amplitude of incident strain field is where Z s is the impedance of the static medium.Based on equation ( 28), the amplitudes of transmitted and reflected strain fields are where Z 1 , Z 2 are the effective impedances for the forward and backward propagating waves.According to equation (1), the strain and momentum density field should be continuous across the temporal boundary.However, it should be noted that μt and μr are not the realistic fields.If the wavenumber k is small, according to equations ( 3) and (A3), the realistic momentum density and strain at t > 0 truncated up to the first order of η are given as where Z 0 = √ κ 0 ρ 0 and Combining equations ( 39)-( 43) and imposing the temporal boundary conditions, we obtain the transmission and reflection coefficients with respect to the momentum density field as Equation ( 42) reveals that the transmission and reflection depend on the wavenumber of acoustic wave when k is not zero.Moreover, since a 1 , a 2 rely on the initial time-modulation phases ϕ 1 , ϕ 2 , see equation ( 41), the transmission and reflection also depend on how the temporal boundary is truncated.In the vanishing wavenumber limit k → 0, reduces to the local form [71] Therefore, the nonlocal and local EMTs produce different transmission and reflection coefficients for a temporal boundary even when the effective constitutive parameters based on the two EMTs are almost the same for small wavenumbers.Especially, when Z s = Z 0 , there will be no reflection based on the local EMT, which is in sharp contrast to the case that the reflection coefficient is nonzero and depends on the initial modulation phases ϕ 1 , ϕ 2 based on the nonlocal EMT.For the sake of simplicity, we consider that only the mass density is temporally modulated, namely ∆ 1 = 0, ∆ 2 = 0.4.Then a 2 = 0, and considering that Z 1 = Z 2 ≈ Z 0 for a small wavenumber, the reflectance predicted by the local and nonlocal EMTs are 0 and respectively.According to equation (46), the reflectance is almost proportional to the square of wavenumber and varies sinusoidally with 2ϕ 2 when the wavenumber is small.As a consequence, the smaller the wavenumber is, the tinier the reflectance will be.The reflectance of the acoustic wave incident to the temporal boundary can be numerically calculated by full wave  44) and ( 45) and the FDTD solutions are shown by the blue solid, red dashed and black symbol lines, respectively.In the simulations, κs = κ 0 , ρs = ρ 0 , ∆ 1 = 0, ∆ 2 = 0.4, k = 0.15K.simulations in conjunction with the Fourier transform process.We used the finite-difference time domain (FDTD) method to obtain the time-dependent momentum density field µ (t), then imposed the Fourier transform to acquire the spectra µ (ω).In figures 4(a) and (c), we plotted the normalized value |µ (ω) /µ i | versus the initial modulation phase ϕ 1 when the impedance is matched Z s = Z 0 for the initial modulation phase differences being ∆ϕ = 0 and ∆ϕ = π /2, respectively.In both cases, ϕ 1 = π /2 (namely ϕ 2 = π /2 and π , respectively).Here, the wavenumber is chosen as k = 0.15K.It is seen that there are three peaks located at ω = −0.8493Ω,− 0.1507Ω and 0.1507Ω, respectively.The second and third peak are corresponding to the reflected and transmitted modes, which are discussed in equation (39), with the peak values denoting the normalized amplitudes of the two modes.The frequencies of the two peaks ω = ±0.1507Ωaccord with the results ±0.1505Ω given by equation ( 27) very well and deviate from the results ±0.15Ω based on the local EMT, again confirming the accuracy of the nonlocal EMT.The mode corresponding to the first peak is not predicted in either local or nonlocal EMTs.Frequency of the first peak in figures 4(a) and (c) is outside the first Brillouin zone −Ω/2 ⩽ ω ⩽ Ω/2 belonging to the band of negative-first order [42] and is induced by the frequency down-conversion due to temporal discrete translational symmetry of the system.Therefore, it should not be excited when the medium can be assumed temporally homogeneous.However, we note that the amplitude of this mode is vanishing small compared with the transmitted and reflected modes, which guarantees the applicability of EMT.The reflectance of the acoustic wave for the temporal boundary is then well approximated by |µ (ω = −0.1507Ω)/µ i | 2 .We showed the reflectance as functions of the initial modulation phase ϕ 1 for ∆ϕ = 0 in figure 4(b).The results are calculated according to equations (44) and (45) and |µ (ω = −0.1507Ω)/µ i | 2 and shown by the blue solid, red dashed and black symbol lines, respectively.While equation (45) produces zero reflectance for all modulation phases, which is quite different from the full wave simulation results, equation ( 44) outputs nonzero sinusoidal reflectance which agree with the full wave simulation results well.Besides some numerical errors, the slight difference between the results given by equation ( 44) and the results based on full wave simulations are incurred by the ignorance of high order terms about the wavenumber in equation ( 44) and the existence of non-zeroth order band modes.If we change ∆ϕ = 0 to ∆ϕ = π /2, according to equation ( 46), the curve reflectance vs. ϕ 1 is also shifted by π /2, which is verified by the full wave simulation results as shown in figure 4(d).As such, the nonlocal EMT can well predict the tiny reflectance at the impedance-matched temporal boundary induced by nonlocality.

Conclusion
In conclusion, we have successfully established a nonlocal EMT with closed-form expressions for effective constitutive parameters within phononic temporal metamaterials.These metamaterials feature periodic time modulation of stiffness and mass density, employing a general multiscale technique.Notably, our research has unveiled that the application of pure temporal modulation can induce Willis coupling, leading to symmetric band structures in these materials, even in the presence of non-zero Willis coupling coefficients.This insight reveals that, in addition to breaking inversion symmetry, temporal modulation offers an alternative route to achieve reciprocal Willis metamaterials.The accuracy and computational efficiency of our nonlocal EMT have been rigorously validated through the calculation of bulk dispersions for two distinct types of phononic temporal metamaterials and the evaluation of minute reflectance at impedance-matched temporal boundaries.
Our nonlocal EMT stands as a powerful instrument for crafting phononic metamaterials through temporal modulations.This approach can effectively cater to diverse applications, including tailoring the flow of acoustic waves for communication purposes, amplifying Willis coupling and nonlocal effects for advanced sensing and diagnostics.
However, it is important to emphasize that our nonlocal EMT is specifically applicable to nondispersive and lossless materials, where the time-dependent stiffness and mass density are purely real and independent of frequency.In cases where material dispersion and loss are factors, more realistic models are necessitated, prompting necessary modifications to the wave equations.Furthermore, our nonlocal EMT is limited to uniform temporal modulations.Since real-world materials typically exhibit some level of temporal dispersion and loss, and space-time modulations could induce more novel phenomena [51,64,72], we acknowledge that the development of EMTs for phononic temporal and space-time metamaterials with strong material dispersion remains a direction for future research.

Appendix C. Plane wave expansion method
Based on the Floquet-Bloch theorem [42,73], the acoustic fields inside the phononic temporal metamaterial can be expressed as where Q is the quasi-energy, and ↔ H F is the Floquet matrix expressed as The band dispersion Q (k) can be obtained by solving the eigenvalue equation (C3).

Figure 1 .
Figure 1.(a) γ 1 (solid lines) and γ 2 (dashed lines) as functions of ρ B /ρ A for different κ B /κ A .(b) The dispersion of the phononic temporal metamaterial calculated by different means.Q is the quasi-energy, which is the effective frequency in the temporal metamaterial.The shaded region denotes the k-gap.The parameters used in (b) are κ B /κ A = 0.5, ρ B /ρ A = 0.2.For both (a) and (b), ∆t 1 = 0.5.

Figure 2 .
Figure 2. The relative errors of calculated band dispersions by different means as functions of the wavenumber.The parameters used are the same as those of figure 1(b).