Transmission properties of time-dependent one-dimensional metamaterials

We solve the wave equation with periodically time-modulated material parameters in a one-dimensional high-contrast resonator structure in the subwavelength regime exactly, for which we compute the subwavelength quasifrequencies numerically using Muller's method. We prove a formula in the form of an ODE using a capacitance matrix approximation. Comparison of the exact results with the approximations reveals that the method of capacitance matrix approximation is accurate and significantly more efficient. We prove various transmission properties in the aforementioned structure and illustrate them with numerical simulations. In particular, we investigate the effect of time-modulated material parameters on the formation of degenerate points, band gaps and k-gaps.


Introduction
Numerous papers have tackled the problem of understanding and manipulating wave propagation in two-and three-dimensional systems with subwavelength resonant structures [5,24,31].Systems of subwavelength resonant structures are of particular interest due to their ability to manipulate waves at subwavelength scales in two-and three-dimensional materials [26][27][28].Such media are made up of a background medium and highly contrasting inclusions, which we call subwavelength resonators.The fact that these inclusions are highly contrasting leads to subwavelength resonances, frequencies at which the resonators interact with incident waves with wavelengths of possibly larger magnitudes [22].This kind of structure appears in various application areas.Subwavelength resonances in highly contrasted structures can be found, for instance, in elastic media [20,29], in plasmonic particles [15,[17][18][19], Helmholtz resonators [16,27] and in dielectric high-index particles [14,32].The plethora of applications of subwavelength resonances make this topic of more general scientific interest.
Wave propagation through a two-or three-dimensional structure with highly contrasting resonators is modelled by a high-contrast Helmholtz problem [7].It has been shown that the high material contrast within the structure is a key assumption for the existence of resonant behaviors at subwavelength scales [9,33].The way the aforementioned Helmholtz problem is approximately solved is to use single-layer potentials based on the fundamental solution of the Laplace problem ing equations exactly in order to find the subwavelength quasifrequencies, for which we make use of the Dirichlet-to-Neumann approach.In Section 4 we provide a brief explanation of Muller's method -the root-finding algorithm used to solve the Helmholtz equations.In Section 5 we shift our attention to a further novel contribution of this paper, which consists of the introduction of a capacitance matrix approximation of the subwavelength quasifrequencies.We prove that such a discrete approximation is a suitable replacement for the numerical scheme solving the wave problem exactly.Lastly, we move on to apply the capacitance matrix approximation to investigate the formation of band gaps, k-gaps and degeneracies and analyze the reciprocity of the wave propagation in Section 6.We summarize our results in Section 7.
2 Problem formulation and preliminary theory

Problem formulation
We seek to solve the one-dimensional wave equation on a domain composed of contrasting materials.In this section, we first introduce the setting which we shall consider in the remainder of this paper.Moreover, we define the material parameters to be time-dependent and assume quasi-periodic boundary conditions.
We consider the case of a one-dimensional system of periodically reoccurring chains of N disjoint subwavelength resonators D i := (x − i , x + i ), where (x ± i ) 1≤i≤N are the 2N boundary points of the resonators satisfying x + i < x − i+1 , for any 1 ≤ i ≤ N − 1.We denote by (x ± i ) i∈N the infinite sequence defined by x ± i+N := x ± i + L, where L ∈ R >0 is the period of an infinite chain of resonators.Furthermore, we denote the length of the i-th resonator D i by ℓ i := x + i − x − i , and the length of the gap between the i-th and the (i + 1)-th resonator by ℓ i(i+1 1 throughout this paper.We refer to Figure 1 for an illustration of the hereby introduced setting.In what follows, we denote by Y := (0, L) the periodic unit cell and by the union of the N resonators in the unit cell.With this notation, the region within R which is taken up by the resonators, is given by

Time-dependent material parameters
We assume that the material parameter distributions are periodic in x with period L and in t with period T := 2π/Ω and are given by Here, ρ and κ represent in acoustics the density and the bulk modulus of the material, respectively, and Ω is the frequency of the time-modulations of the material parameter distributions.We define the contrast parameter and the wave speeds by respectively.To achieve subwavelength resonance, we will assume that the contrast parameter δ is small: Typically, the most interesting regime of the frequency of modulations of ρ i (t) and κ i (t) is Ω = O(δ 1/2 ), i.e., of the same order as the static subwavelength resonances [13].This allows strong coupling between the time modulations and the response time of the structure.We aim at finding ω = O(δ 1/2 ) such that the wave equation has a non-trivial solution u(x, t) which is essentially supported in the low-frequency regime.
By substituting the time-harmonic wave field u(x, t) = ℜ v(x, t)e iωt into the wave equation ( 6), we obtain Due to the assumption that u(x, t)e −iωt is T -periodic with respect to time t, we write the Fourier series expansion Note that any L 2 -function v n (x) can be decomposed into a superposition of Bloch waves as follows: where α is the so-called momentum and vn (x, α) is L-periodic in x.The function vn is defined by Thus, we can write Inserting the expansion (11) into the differential equation ( 7), we conclude that for any n ∈ Z, vn must satisfy for x ∈ R and t ∈ R.
Recall that we have assumed the chain of N resonators to be repeated periodically with period L. Therefore, we study the one-dimensional spectral problem in the unit cell (0, L) for the quasi-periodic function v n (x, α) := vn (x, α)e iαx : where we use the notation The functions v * i,n (x, α) and v * * i,n (x, α) are defined in each resonator D i through the convolutions where r i,m and k i,m are the Fourier series coefficients of 1/ρ i (t) and 1/κ i (t), respectively.Furthermore, we define the wave number outside and inside the resonators corresponding to the n-th mode through respectively.We assume that the time-modulations of ρ i and κ i have finite Fourier series in each resonator D i , that is, for some M ∈ N satisfying M = O Ä δ −γ/2 ä , for some γ ∈ (0, 1) [13].Note that the solution to (13) is invariant under scaling.Hence, we can assume the solution to be normalized.As u is continuously differentiable in t, we have where ∥ • ∥ 2 denotes the L 2 -norm on (0, L).Due to folding (see Definition 5.3), we need to specify the subwavelength quasifrequencies in terms of the oscillations in their associated modes [13].As said before, the subwavelength quasifrequencies are those associated with Bloch modes essentially supported in the low-frequency regime as δ → 0. Therefore, we shall assume that there exists some as δ → 0, where the sequence of functions (v n ) n∈Z is a nontrivial solution to (13).
In order to perform some numerical and analytic analysis in this regime, we adapt the Dirichlet-to-Neumann approach of [21,22] to the one-dimensional, quasi-periodic and timemodulated case to solve (13).

Exact solution
In this section we seek to solve the coupled Helmholtz problem (13) exactly.We first present a characterization of the solution to the exterior problem and then to the interior problem.Lastly, we use the Dirichlet-to-Neumann map to derive a system of equations based on the boundary condition.

Exterior problem
In this section we seek to characterize the Dirichlet-to-Neumann map of the Helmholtz operator on the domain (0, L) with the quasi-periodic boundary condition.
We denote the Sobolev space of quasi-periodic complex-valued functions by H 1 per,α (R).We also denote by C 2N,α the set of quasi-periodic boundary data f ≡ (f ± i ) i∈Z such that where f + i (resp.f − i ) refers to the component associated with x + i (resp.with x − i ).The space of such quasi-periodic sequences is clearly finite-dimensional; specifically, it is of dimension 2N .
The following lemma from [22] provides an explicit expression for the solution to the exterior problem on R \ (D + LZ).Lemma 3.1.Assume that k n = (ω +nΩ)/v 0 , for some fixed n ∈ Z, is not of the form mπ/ℓ i(i+1) for some non-zero integer m ∈ Z\{0} and index 1 ≤ i ≤ N .Then, for any quasi-periodic sequence Furthermore, when k n ̸ = 0, the solution v α f,n reads explicitly where, for fixed n ∈ Z, α n i and β n i are given by the matrix-vector product Proof.Identical to [22, Lemma 2.1].
Definition 3.1.For any k n ∈ C, for fixed n ∈ Z, which is not of the form mπ/ℓ i(i+1) for some m ∈ Z\{0} and 1 ≤ i ≤ N − 1, the Dirichlet-to-Neumann map with wave number k n is the linear operator T k n ,α : C 2N,α → C 2N,α defined by where v α f,n is the unique solution to (21).
Using the exponential Ansatz presented in Lemma 3.1 to solve (21) gives rise to a closed form definition of the Dirichlet-to-Neumann map, which we introduce in the following proposition.Proposition 3.2.For fixed n ∈ Z, the Dirichlet-to-Neumann map T k n ,α admits the following explicit matrix representation: for any So far, we have found a way to solve the exterior problem explicitly for some given boundary data, as stated in Lemma 3.1.Moreover, we have provided an explicit matrix representation of the Dirichlet-to-Neumann map in Proposition 3.2, which we will make use of when dealing with the Neumann boundary condition of (13) in order to solve the interior problem.

Interior problem
Having dealt with the exterior problem in the previous section, we now focus on the solution of the interior problem.We can formulate the interior part of problem ( 13) using the Dirichlet-to-Neumann map, which leads to for n ∈ Z. Recall that v * i,n and v * * i,n are the convolutions defined by (15).We now recall the definition of a subwavelength quasifrequency [13].
Next, we state the following lemma, which provides us with the solution to the interior problem upon using an exponential Ansatz.Lemma 3.3.For each resonator D i , for i = 1, . . ., N , the interior problem (27) can be written as an infinitely-dimensional system of ODEs, ∆A i v i + B i v i = 0, with where we define the coefficients and v i,j = v j | D i .By the definition of A i , the matrix is invertible and we may write Let { λi j } j∈Z be the set of all eigenvalues of C i with corresponding eigenvectors {f j,i } j∈Z .Using the square-roots ±λ i j of the eigenvalues λi j , the solution to the interior problem (27) over D i takes the form Remark 3.3.Note that in order to use the Ansatz (31) to obtain a solution numerically, we introduce a truncation parameter K ∈ N such that K ≥ M and such that ( 19) is satisfied.Thus, we truncate the infinitely-dimensional vector and matrices to . Hence, we use the following Ansatz as a solution to the interior problem: Moreover, we denote the n-th entry of f j,i by f j,i n .We shall illustrate the influence of the parameter K in Figure 2 and Figure 3b.
Using the Ansatz (32) and the transmission conditions in (27), we construct a non-linear system of equations in ω, which characterizes the coefficients {(a i j , b i j )} −K≤j≤K ⊂ R 2 for each i = 1, . . ., N .Theorem 3.4.Let K be a fixed and sufficiently large truncation parameter.The subwavelength quasifrequencies ω to the wave problem (27) are approximately satisfying, as δ → 0, the following truncated system of non-linear equations: where the unknown vector is and the matrices G n,j = G n,j (ω) and V n,j = V n,j (ω) are given by Here, we use the convention f k,i K+1−n = 0, if |n| > K.We define the 2N × 2N (2K + 1) matrix such that (33) can equivalently be expressed as A n (ω, δ)w = 0, for all n = −K, . . ., K, where To this end, (33) can be written as Proof.Using the Ansatz (32), we express any v n on the boundary of a resonator D i as Inserting this into the boundary condition we obtain which can be written as (38) by evaluating it for each −K ≤ j ≤ K.The value of K directly impacts the runtime of solving the root-finding problem arising from (38) and the accuracy of the result.Figure 2 shows how the runtime of Muller's method solving the root-finding problem depends on the truncation parameter K. Due to calculation of the eigenvalues of (2K + 1)2N × (2K + 1)2N and (2K + 1) × (2K + 1) matrices, the runtime increases algebraically in K. Contrarily, we require K to be large in order for the Ansatz (32) to be accurate.Thus, we seek to introduce an alternative way of computing the subwavelength quasifrequencies ω α i , which does not involve the choice of a truncation parameter K.For this matter we shall introduce the formulation of an approximation formula based on the capacitance matrix in Section 5.1.

Explicit choice of parameters
We now assume that ρ i and κ i are specifically given by for all 1 ≤ i ≤ N , where ε ρ,i , ε κ,i are the amplitudes of the time-modulations and ϕ ρ,i , ϕ κ,i the phase shifts.This means that we can set M = 1 with the Fourier coefficients defined as follows: The remainder of this paper treats the case of material parameters given by (40).We will modulate the amplitudes ε ρ,i and ε κ,i in our numerical experiments to investigate the effect of time-modulated materials on propagating waves.

Muller's method
We solve problem (13) with the help of Muller's method.In particular, we use Lemma 3.4 to construct a 2N (2K + 1) × 2N (2K + 1) system of equations A * (ω, δ)w = 0, which provides the correct coefficients a n i and b n i of the n-th mode v n in each resonator D i .We seek to find the subwavelength quasifrequencies ω α (δ), which are those values of ω for which the interior problem ( 27) admits a non-trivial solution.Note that these are exactly the values of ω for which A * (ω, δ) is non-invertible, i.e., the matrix A * (ω, δ) has a zero-eigenvalue.Therefore, we define the function whose zeros we must find, for a fixed δ.Note that σ (A * (ω, δ)) is defined to be the spectrum of the matrix A * (ω, δ).In order to find the zeros of the non-linear function f (ω), we use Muller's method upon three initial guesses per root.For a detailed explanation of Muller's method we refer the reader to [10, Section 1.6].One of the reasons for using Muller's method to solve the root-finding problem is that, unlike other root-finding algorithms, Muller's method is well-suited for complex-valued problems [10, Section 1.6].Muller's method requires the definition of three initial guesses to find a zero of f (ω).In our numerical computations we make use of the already known definition of the capacitance matrix C α in the static case [1], as further explained in Appendix A. Namely, we compute the eigenvalues λ α i , 1 ≤ i ≤ N , of the static generalized capacitance matrix C α and employ the asymptotic approximation from [22] To initialize Muller's method we use (44) and two perturbations of this value.The use of Muller's method requires the definition of a tolerance, which we chose consistently to be 10 −12 .By the definition of Muller's method we need to supply the algorithm with three initial guesses in order to find a zero of f (ω), which is not trivial.Furthermore, for fixed K, the runtime of Muller's method grows exponentially in N , as illustrated in Figure 2. Therefore, we seek to introduce an alternative characterization of the subwavelength quasifrequencies ω α , for which we do not require the exact solution of (27).In view of this, we introduce a discrete approximation of (13) in Section 5.1.
Lemma 5.1.As δ → 0, the functions v * i,n (x, α) are approximately constant inside the resonator: For simplicity of notation, for any smooth function f : R → R, we define I ∂D j [f ] by Passing the convolution in the definition (15) of v * i,n (x) to the time domain, we obtain Recall that 1/ρ i (t) has a finite number of Fourier coefficients {r i,m } −M ≤m≤M .Thus, ρ i (t) has an infinite number of Fourier coefficients, which we denote by {r i,m } m∈Z .Therefore, it follows from the definition u(x, t) = ρ i (t)u * i (x, t), for x ∈ D i , that v n (x) can be expressed through Lemma 5.2.The expression (48) can be extended to the whole space given by where the functions V α i : R → R are solutions to the following equations: Proof.By Lemma 5.1, we have Between D i and D i+1 , v n (x) satisfies the following Helmholtz equation: The solution is hence given by for some coefficients A and B. We also have Without loss of generality, we can assume that x − i = 0 and the solution has the following expansion: This leads to a linear interpolation between the resonators.Hence, we can extend the characteristic function χ D i to V α i for all i = 1, ..., N .
Applying the operator I ∂D i to v n , we conclude that where the coefficients of the capacitance matrix C α are given by We give the explicit definition of the coefficients of C α in Appendix A.
On the other hand, we can use the transmission conditions in (13) to obtain an alternative characterization of I ∂D i [v n ] as follows: Equating ( 55) and (57) leads to Next, we define the following functions: We have Letting δ → 0, we obtain To this end, we obtain the following theorem.
Theorem 5.3.Assuming that the material parameters are given by (3), as δ → 0, the quasifrequencies in the subwavelength regime are, at leading order, given by the quasifrequencies of the system of ordinary differential equations where w i (t) := ρ i (t)c i (t).
Remark 5.1.We note here that the capacitance formulation given by (62) does not depend on the material parameter ρ i (t).This is a direct consequence of the wave equation ( 6) and the transmission condition cancelling out the ρ i (t)-dependency.However, (62) only holds true for small δ and we may assume that the quasifrequencies for δ = O(1) do depend on ρ i (t).
Remark 5.2.An equivalent way of formulating the ODE (62) is through the following system of ODEs: where M α (t) = δκr ρr W 1 (t)C α W 2 (t) + W 3 (t) with W 1 , W 2 and W 3 being diagonal matrices defined as for i = 1, . . ., N , with Our numerical results presented in Figure 3 corroborate our analytically proven claim that the capacitance matrix approximation is an efficient and effective alternative to the exact computation of the subwavelength quasifrequencies using Muller's method.Here, we define the absolute error to be given by where ω α i,muller are the quasifrequencies calculated by Muller's method and ω α i,cap are the quasifrequencies calculated through the capacitance approximation.
Based on our results shown in Figure 3b we conclude that one should choose K ≥ 3 in order to minimize the absolute error between the result obtained by the capacitance approximation and Muller's method.However, it is important to note that as K increases, the calculation of the eigenvalues of the (2K + 1)2N × (2K + 1)2N and (2K + 1) × (2K + 1) matrices becomes more complex and thus more error-prone.Therefore, the choice of the parameter K is a delicate matter and we shall replace the exact computation of ω α i by the herein introduced capacitance approximation, which is independent of K. Nevertheless, the capacitance matrix depends on the number N of resonators.But, the influence of N on the runtime of the capacitance approximation is less significant than the influence it has in the exact computation.We conclude this from Figure 4.This observation further motivates the introduction of the capacitance matrix approximation.

Asymptotic analysis
In order to analyze the reciprocity properties of the wave transmission we make use of asymptotic Floquet analysis developed in [12] and we closely follow [4].For simplicity, we consider the modulation amplitudes of ρ i and κ i to be the same over all resonators, i.e., ε ρ,i = ε κ,i = ε, for all 1 ≤ i ≤ N .Then we assume that the matrix M α (t) in ( 63) is analytic in ε, whence, we can expand M α (t) as follows [4]: for small ε > 0. If ρ i (t) and κ i (t) have finitely many Fourier coefficients, we can assume that the series (67) converges for any |ε| < ε 0 , for some ε 0 > 0 [3,12].Note that we omit the superscript α in the remainder of this section for the sake of convenience.Next, we rewrite the second order ODE (63) into the first order ODE [4] dy dt (t) = A(t)y(t), A(t) := where Id N is the N × N identity matrix.By Floquet's theorem, the fundamental solution of (68) can be written as X(t) = P (t)e F t , for some matrices P (t) and F [37].As a consequence of M (t) being analytic in ε, we can write [4] The coefficients A 0 and P 0 are not time-dependent, as they correspond to ε = 0, which represents exactly the static case.Due to the T -periodicity of the material parameters, A(t) is T -periodic and, thus, A j (t) are also T -periodic, for all j ≥ 1.Hence, we may write [4] A We now aim to derive asymptotic expansions of the eigenvalues f = f 0 + εf + . . . of F in ε.
Assume the first coefficient A 0 in the expansion of A(t) to be diagonal.Then, according to [4], we have with m i being the folding number of (A 0 ) ii , which is defined as follows.
We are specifically interested in investigating perturbations due to the modulations at a degenerate point f 0 of F 0 , which can be obtained through folding, for which we make use of the following lemma from [3]. numerically.We use the capacitance matrix approximation in order to conduct some numerical experiments under different conditions.We seek to analyze the so-called band structure of the material, which describes the quasifrequency-to-momentum relationship of the propagating waves [4].We are especially interested in the occurrence of band gaps and k-gaps as a consequence of time-modulated material parameters.Previous work [11] has proven the occurrence of subwavelength gaps in three-dimensional high-contrast materials in the static regime.
Moreover, we want to understand the effect of time-modulation on the reciprocity of wave transmission properties.The reciprocity of wave transmission is defined as follows.Definition 6.1.We say that a wave propagates reciprocally if for each α in the space Brillouin zone Y * , the quasifrequencies of the wave problem (6) at α coincide with the quasifrequencies at −α [4].It becomes apparent from Figure 5c that there is a degenerate point at α = 0 if the gap size between each resonator is equal, which can be treated equivalently to the case of N = 1 resonator in the unit cell.
Comparing Figure 5 with Figure 6, it becomes apparent that modulating κ in time turns degenerate points into k-gaps.Furthermore, measuring the size of the k-gaps shows that the gaps forming in the regime α < 0 do not have the same size as the gaps forming in the regime α > 0. This means that the wave transmission is non-reciprocal in the time-modulated case.The following theorem has been proven in higher dimensions in [4,Theorem 8], but can equivalently be proven in the one-dimensional case.

Conclusion
In this paper we have provided the mathematical foundation to solve the quasi-periodic Helmholtz equation in one dimension with periodically time-dependent material parameters.We presented a discretization of the problem (13), which led to a scheme solving the interior problem exactly up to a negligible numerical error induced by Muller's method.The solution of this scheme involves the calculation of eigenvalues and eigenvectors of large matrices which is very time-consuming.
Furthermore, we have introduced a novel capacitance matrix approximation to the subwavelength quasifrequencies in one dimension assuming the problem to be quasi-periodic and the material periodically time-modulated, which is equivalent to the approximation formula valid in higher dimensions [13].This approximation formula is advantageous because it recovers the quasifrequencies in the subwavelength range much more efficiently.In the static case, the subwavelength quasifrequencies can be approximated by the formula ω α i ≈ v r λ α i δ, for all 1 ≤ i ≤ N , where λ α i are the eigenvalues of the generalized capacitance matrix C α .Whereas in the time-dependent case, the quasifrequencies ω α i are obtained through the solution of an ODE which depends on the capacitance matrix and the material parameters.We have showed in Figure 4 that for increasing N , approximating the subwavelength quasifrequencies with the help of the capacitance matrix is indeed much faster than computing them with Muller's method.Moreover, the capacitance approximation does not depend on the truncation length K of the system of Helmholtz equations (13), as opposed to the exact solution.We showed numerically in Figure 2 that the runtime of the exact computation of ω α i depends heavily on the truncation parameter K.
Our numerical analysis led to the conclusion that under time-modulated material parameters, the wave transmission is non-reciprocal, which aligns with the asymptotic analysis in Section 5.Moreover, it became apparent that periodic time-modulations in the κ i 's lead to the formation of k-gaps.

A Code availability
The codes that were used to generate the results presented this paper are openly available under https://github.com/liorarueff/1D_quasiperiodic_timemod.

B Capacitance matrix approximation to the static problem
In this section, we recall results from [1] regarding the capacitance matrix approximation to the static problem.

Figure 1 :
Figure 1: An illustration of the one-dimensional setting for N = 3 resonators in the unit cell.

Figure 2 :
Figure 2: The time it takes to formulate the root-finding problem and for Muller's method to solve it depending on the truncation parameter K.The runtime depends algebraically on the parameter K.

Figure 3 :
Figure 3: We compare the quasifrequencies obtained by Muller's method with the quasifrequencies obtained through the capacitance matrix approximation.

Figure 5 :
Figure 5: Subwavelength quasifrequencies for three resonators repeated periodically in the static case.The figures on the right-hand side illustrate the setting corresponding to the numerical results shown in the left-hand side figures.