Complex Dispersion Relation Recovery from 2 D Periodic Resonant Systems of Finite Size

The complex dispersion relations along the main symmetry directions of two-dimensional finite size periodic arrangements of resonant or non-resonant scatterers are recovered by using an extension of the SLaTCoW (Spatial LAplace Transform for COmplex Wavenumber) method. This method relies on the analysis of the spatial Laplace transform instead of the usual spatial Fourier transform of the measured wavefield in the frequency domain. We apply this method to finite dimension square periodic arrangements of both rigid and resonant scatterers embedded in air, i.e., to finite size sonic crystals and finite size acoustic metamaterials, respectively. The main hypothesis considered in this work is the mirror symmetry of the finite structure with respect to its median axis along the analyzed direction. However, we show that the method is robust enough to provide excellent results even if this hypothesis is not fully satisfied. Effectively, a minor asymmetry could be considered as a side effect when the structure is large enough because Laplace transforming the field along the main symmetry directions also implies averaging the field in the perpendicular one. The calculated complex dispersion relations are in excellent agreement with those obtained by an already validated technique, like the Extended Plane Wave Expansion (EPWE). The methodology employed in this work is intended to be directly used for the experimental characterization of real 2D periodic and resonant systems.


Introduction
Periodic structures or locally resonant systems, also known as metamaterials, have revolutionized the control of waves these last years [1,2].The main physical properties of these structured media, as for example the opening of band gaps [3,4], the slow wave frequency band [5,6], or the spatial filtering [7], among others, can be directly derived from the complex dispersion relation [8][9][10][11][12].In addition to dispersive effects, wave attenuation can be caused by factors such as geometric attenuation or intrinsic material loss (e.g., heat or viscous dissipation).In all cases, wave attenuation can be interpreted in terms of complex wavenumbers, which are functions of the angular frequency ω via a complex dispersion relation k(ω).The real part of the wavenumber translates the propagative (dispersive) properties of the system while its imaginary part is related to the attenuation.This rich but sometimes hardly to interpret information is of crucial relevance for the design of metamaterial-based devices.Therefore, the recovery of the complex dispersion relation, i.e., the measurement or calculation of both the dispersion and the attenuation of waves in a finite system is of utmost importance in wave physics and material science, since it provides insightful information on the tested system also allowing the characterization of the effective properties of the material.
Several methods have been proposed during these last years to calculate the complex dispersion relation of waves in discrete or continuous systems.The extended plane wave expansion [9,10] or methods based on Multiple Scattering Theory (MST) [13,14] are efficient tools to evaluate the complex dispersion relation in non-resonant structures made of scatterers of simple geometry.Numerical approaches are more relevant [15] when the scatterer geometry is more complicated or the scatterers are resonant.
Recently, the SLaTCoW (Spatial LAplace Transform for COmplex Wavenumber) [11,12] method has been proposed as a general procedure for the recovery of the complex dispersion relation.It has been successfully applied to several systems of different kind, such as Lamb-like modes in lossy materials in the low kHz regime, elastic waves in a Duralumin plate in the MHz regime at the Zero-Group Velocity (ZGV) point, surface acoustic waves in a two-dimensional microscale granular crystal adhered to a substrate near 100 MHz [11] and also spoof surface acoustic waves at a lossy resonant metasurface [12].The method is based on the spatial Laplace transform (LT) instead of the usual spatial Fourier transform (FT) of the measured wavefield in the frequency domain.The LT, providing information on both the real and imaginary parts of the wavenumber, is analyzed by minimizing a correctly chosen cost function.The complex wavenumbers (as well as the complex amplitudes) of the modes are reconstructed even when they are interacting one with another or overlapping.Moreover, the SLaTCoW method is particularly suitable for the experimental characterization of real periodic systems of finite-size (resonant or non-resonant), which, while exhibiting similar properties as systems of infinite extent such as band gaps [16][17][18] or other dispersive effects [19,20], can not be experimentally characterized using methods considering infinite structures, such as EPWE, MST or solving an eigenvalue problem using the Finite Element Method (FEM).
In this work, we present an advanced version of the SLaTCoW method to characterize the complex dispersion relation of two-dimensional systems.We apply this methodology to both periodic media as Sonic Crystals (SC) and locally resonant acoustic materials.The article is organized as follows.The methodology and the extension of the method is presented Section 2. Section 3 depicts the different two-dimensional configurations the methodology is applied to.Finally, results are provided and validated in Section 4. Numerical reconstructions from finite size two-dimensional systems agree well with complex dispersion relation calculated from analytical techniques.A summary and the concluding remarks are given in Section 5.

SLatCoW
The SLatCow method relies on the analysis of the spatial LT of the acoustic field p(x, ω) in the frequency domain.The analysis is performed in the linear harmonic regime at the circular frequency ω with the implicit time dependence e −iωt .For two-dimensional problems, the pressure field is recorded on a surface map of dimensions x ∈ [0, L s ] and y ∈ [0, L s ], as depicted in Figure 1.The excitation is assumed to take the form of a plane incident wave along the positive x-axis.Therefore, the pressure field takes the form where pm is the complex amplitude of the m-th mode, K m = (K m x , K m y ) is the complex wavenumber of the m-th mode, M is the set of modes, and Π (z, L) is the gate function equal to 1 when z ∈ [0, L] and equal to 0 elsewhere.For the sake of explanation the ω dependence is further dropped.
The complex wavenumber is also defined as K m = k m r + ik m i , with k m r and k m i are two component vectors, one along the x-axis and the other along the y-axis.Please note that along the specific direction X, K m • X = k m r + k m i , with k m r ≥ 0 and k m i ≥ 0 ensuring that Equation (1) only involves forward propagating waves along positive X.Applying the usual spatial FT only enables the recovery of k m r .In order to recover both real and imaginary parts of K m , a spatial LT is applied to p(x) denoted x dx, where s is a complex valued vector s = s r + is i , with s r and s i being two component vectors, again one along the x-axis and the other along the y-axis.The spatial LT takes the form In what follows, we are only interested in specific directions, namely ΓX and ΓM, obtained by rotating the structure in the area [0, For a structure of infinite lateral extent, i.e., along the y-direction, the continuity of the real part of the K m y is ensured between the outer and the inner domain.In the present case of finite lateral dimensions, every modes are excited because of wave diffraction by the edges.Fixing s y = 0 ensures that we are only looking for modes propagating along the x-axis.Nevertheless, the LT reduces in the general case to . Considering a symmetric structure with respect to y = L s /2 excited by a plane wave along the x-axis ensures that K m y = 0, ∀m ∈ M [21].Therefore, the two dimensional LT reduces to Please note that the symmetry is ensured in the case depicted in Figure 1a along both directions ΓX and ΓM, while it is only ensured along the ΓX direction in the case depicted in Figure 1b.Along the ΓM direction, the symmetry condition is not satisfied because of the rotation of the quarter-wavelength resonator (QWR) openings.Nevertheless, this will be considered and proven as a side effect in what follows.Note also that Laplace transforming along the y-direction for s y = 0 consists in averaging the pressure field measured along this axis or alternatively Fourier transforming the field with k y = 0.The problem becoming also independent from y, we drop the x dependence on both s and K m , so as The spatial LT of the field therefore reduces to a form close to the one used in Ref. [11].The recovery of K m is performed in a very similar way for each frequency by minimizing the following cost function where P mes (s) = P mes (s x , s y = 0) is the spatial LT of the measured field p(x), and | pm | and φ m are, respectively, the theoretical amplitude and phase of the m-th mode.The minimization is performed under constraints with the Nelder-Mead simplex algorithm ( R Matlab function fminsearchbnd).

Spatial Sampling Criteria for Dispersion Relation Recovery
As stated previously, the SLaTCoW method relies on analyzing the pressure field in the complex wavenumber and frequency domains (p(s r , s i , ω)).Hence, a double transformation of the acoustic pressure is needed to apply the method.First, time to frequency ( p(x, t) → p(x, ω)) through the temporal FT, then space to wavenumber (p(x, ω) → P(s r , s i , ω)) through the spatial LT.While the temporal FT is not limited by or independent from the spatial periodicity of the structure, the spatial LT greatly depends on both the size of the structure, i.e., the number of periodic elements, and the spatial sampling of the acoustic field.To illustrate this, we consider a spatial signal as represented in Figure 2a and applied successively a temporal and a spatial FT.By analogy with temporal signals and following basic sampling theory, the resolution in s i -space is determined by the length of the signal, ∆s i = 1/L s , where L s = ∆xN with ∆x the spatial sampling and N the number of samples of the spatial signal.Lets consider now a finite size one-dimensional periodic structure with lattice constant a formed by alternating two different materials, as depicted in Figure 2b.From the Nyquist-Shannon sampling theorem [22] and considering periodicity in s i -space, the spatial sampling rate 1/∆x has to be equal or greater than twice the highest wavenumber component, i.e., s i /2 ≥ k max /2π, where k max = π/a.As a result, the following condition for the spatial sampling ∆x ≤ a should be satisfied to ensure that no information is lost when performing the FT.In summary, for a finite dimension periodic structure comprising a specific number of unit cells N uc and assuming that the wavefield is measured at a single point per periodic element, the number of these periodic elements (unit cells) determines the number of points and hence the resolution in the s i -space as follows Figure 2c-f illustrates the effects of the variation of the number of periodic elements, and therefore of the spatial sampling used on the FT of the pressure field, used to evaluate the dispersion relation of the multilayer structure shown in Figure 2b.A double FT is applied to numerically simulated spatio-temporal acoustic pressure values, p(x, t).Results are presented for positive and negative values of the wavenumber and are compared to the theoretical dispersion relation calculated using the Rytov formula [23].Figures 2c,d analyze the influence of the number of unit cells on the double FT using the critical sampling rate ∆x = a.The total number of wavevector points in the s i -space is not high enough to extract accurate information to represent the dispersion relation in the case of a structure containing only N uc = 10.Please note that the number of points in the s i -space is distributed along the positive and negative parts of the wavenumber.Figures 2e,f show the consequences of spatially under-sampling and over-sampling, respectively.The former causes the FT to be folded, while the latter creates replicates of the FT.However, no gain in information is obtained by increasing the spatial sampling above its critical value.Please note that the wavenumber resolution can be artificially increased by zero-padding the spatial signal in the space domain, although no actual information will be extracted from this process.

Numerical Set-Up
Two different 2D periodic structures of finite sizes are considered in this study, as shown in Figure 1.Both structures are composed of N × N rigid scatterers arranged in a square lattice with a lattice constant a = 0.085 m, embedded in air, where N = 21.The SC shown in Figure 1a consists of cylindrical rigid scatterers with radius r = 0.003 m.The second periodic system is a metamaterial consisting of square scatterers with QWRs, as shown in Figure 1b.The geometric parameters of the resonator are its external side length, L = 0.05 m, and internal dimensions, d qw = 0.035 m and l qw = 0.04 m.See Refs.[24,25] for specific details on the resonators.The host medium in both cases is air, with density and speed of sound ρ = 1.213 kg/m 3 and c = 343 m/s, respectively.
Numerical experiments are carried out using FEM.A scattering problem in the frequency domain is solved for the two structures described previously.The excitation is an incident plane wave traveling along the positive x-axis.Perfectly matched layers are considered surrounding the numerical domain to approximate the Sommerfeld radiation condition.Following the spatial sampling criteria for periodic structures introduced in Section 2.2, the acoustic pressure is recorded at multiple positions separated by the lattice constant, ∆x = a, giving rise to a total of N × N measurement points.Please note that the structures presented in Figure 1 are used to calculate the complex dispersion relation along the ΓX direction.To retrieve the dispersion relation along the ΓM direction both structures are rotated by 45 degrees while the incident plane wave is unchanged, i.e., remains traveling along the positive x-axis.Note also that the spatial sampling in this case is modified accordingly, ∆x = √ 2a.

Complex Dispersion Relation. SLaTCoW vs. EPWE
We start by considering a reference case which consists in a classical SC of finite size composed of 21 × 21 rigid cylinders arranged in a square lattice and embedded in air, where no losses are considered.In this situation, the imaginary part of the wavenumber is non-zero only because of the evanescent modes created within the band gap.To validate the results obtained using SLaTCoW, we make use of EPWE [9,10], an already established technique that allows us to retrieve both the real and imaginary parts of the wavenumber for non-resonant periodic structures of simple geometries.
Figure 3 illustrates the complex band structure of the finite SC presented in Figure 1 using SLaTCoW and the one obtained using EPWE for an infinite SC with the same geometrical and physical properties.The imaginary part of the wavenumber is shown in the right (left) panel in the ΓX (ΓM) direction, while the real part of the wavevector for both directions is shown in the central panel.Both the real and imaginary parts of the wavenumber show an excellent agreement with the EPWE calculations for both directions of symmetry.Moreover, the typical signature of the Bragg interference on the imaginary part of the wavevector is captured, with a maximum attenuation at the Bragg frequency and smooth variations from this point to the upper and lower limits of the band gap.Minor differences are noted in the band gap frequency range along the ΓM direction, where the imaginary part of the wavenumber using SLaTCoW is slightly lower than the one obtained using EPWE.Small differences are also noticed on the real part of the wavenumber in that frequency range.Both are therefore attributed to the finite size of the SC used for SLaTCoW calculations.Numerical results obtained for a N × N finite structure are compared to the ones obtained using EPWE for the same periodic structure having infinite size.

Complex Dispersion Relation of a 2D Periodic Resonant System
We then focus on the dispersion relation of the 2D acoustic metamaterial of finite size.Specifically, we consider the structure shown in Figure 1b, i.e., a periodic system composed of QWRs arranged in a square lattice embedded in air.Both lossless and lossy cases are considered with visco-thermal losses considered solely inside the QWR resonators [24,25].Figure 4 shows the complex dispersion relation results of the 2D acoustic metamaterial.The real part of the wavenumber is shown in the central panel and results obtained using SLaTCoW are compared to those obtained solving an eigenvalue problem using the FEM over an infinite periodic structure with the same geometric and physical characteristics.Please note that the dispersion relation of the host medium alone, i.e., the air in this case, is also plotted in gray dashed lines.Results show an excellent agreement for all the propagative bands between the results obtained using SLaTCoW and FEM for both ΓX and ΓM directions.However, differences are noticed in between the propagative bands, where the eigenvalue problem does not possess real solution.Effectively, no real solution satisfy the problem ω(k) for a periodic system of infinite extent.Moreover, acoustic modes corresponding to the air medium are also recovered via the SLaTCoW method, because the structure under analysis is of finite size and excited by a plane incident wave initially propagating in a homogeneous medium.Their magnitude are much smaller than those of periodic structure modes and decrease rapidly with the number of unit cells of the finite structure.In the band gap frequency range, air modes become more significant, because no propagative mode in the periodic system exists.They are detected by the algorithm as for example along the ΓM direction around 1500 Hz.This prompts a reduction of the imaginary part of the wavenumber at a small range of frequencies within the resonant gap.A closer look to ΓX results for both lossless and lossy cases, indicate the typical characteristics of resonant and Bragg gaps when the air modes are not detected.The attenuation peak of the resonant band gap (around 1600 Hz) is greatly reduced when the visco-thermal losses are accounted for in the QWRs.On the other hand, the attenuation is unaltered between lossless and lossy cases for the Bragg gap (around 2200 Hz), because the losses are only considered inside the resonators and not in the host medium.

Conclusions
We have presented a methodology to calculate the complex dispersion relation of 2D periodic resonant and non-resonant systems of finite-size via a slight modification of the SLaTCoW method presented in Ref. [11].By assuming the finite size structures to be mirror symmetric with respect to their median axes along the analyzed direction, the SLaTCoW method is extended to a two-dimensional problem allowing for the calculation of the dispersive and attenuation properties of the system at specific symmetry directions, namely ΓX and ΓM.To properly record the necessary acoustic pressure, which happens to be distributed in a surface map to account for 2D structures, a basic spatial sampling criteria applicable to periodic structures is provided and discussed.Two periodic structures with different properties are considered in the study.First, we validated the results provided by the present methodology by comparing them to those obtained via a well-established technique, the EPWE, exhibiting an excellent agreement for both the real and imaginary parts of the wavevector for the two symmetry directions.Then, we show the potential of the methodology to analyze a 2D acoustic metamaterial including QWRs and accounting for visco-thermal losses.The obtained results also show an excellent agreement with an already validated technique for the real part of the dispersion relation and clearly manifest the typical signature of resonant and Bragg gaps in the imaginary part of the wavenumber.Limitations related to the appearance of air modes in the dispersion results, due to the finite size of the structure and its excitation by plane incident waves are also discussed.As shown in Refs.[11,12] this method is applicable to systems of different kind, at different frequency ranges, and should also be applicable for other type of waves, such as TE and TM modes for electromagnetic waves, where there is a close similarity with the 2D acoustic wave equation.Finally, the methodology presented in this work is intended to directly serve for the experimental characterization of real 2D periodic and resonant systems.

Figure 1 .
Figure 1.Set-up of the numerical experiments used to calculate the complex dispersion relation: (a) 2D SC composed of N × N cylindrical scatterers arranged in a square lattice.(b) 2D Metamaterial composed of N × N QWRs.

TheoryFigure 2 .
Figure 2. (a) Example of real and imaginary parts of the spatial pressure signal along the multilayer structure at the normalized frequency ω/Ω 0 = 0.25.(b) Periodic multilayer structure composed of N uc = 50 periodic elements, combining two artificial materials with different acoustic properties.(c-f) Double FT of the pressure field and dispersion relation of the multilayer structure considering different number of periodic elements and different values of the spatial sampling ∆x.

Figure 3 .
Figure 3. Complex dispersion relation for a SC composed of rigid cylinders arranged in a square lattice.Numerical results obtained for a N × N finite structure are compared to the ones obtained using EPWE for the same periodic structure having infinite size.

Figure 4 .
Figure 4.Complex dispersion relation for a periodic resonant system formed by QWRs.Numerical simulations using an infinite structure and solving the corresponding eigenvalue problem in FEM are compared to the results obtained applying SLatCoW from numerical data using a finite structure.