Acoustic metamaterials for sound focusing and confinement

We give a theoretical design for a locally resonant two-dimensional cylindrical structure involving a pair of C-shaped voids in an elastic medium which we term as double ‘C’ resonators (DCRs) and imbedded thin stiff bars, that displays the negative refraction effect in the low frequency regime. DCRs are responsible for a low frequency band gap which hybridizes with a tiny gap associated with the presence of the thin bars. Using an asymptotic analysis, typical working frequencies are given in closed form: DCRs behave as Helmholtz resonators modeled by masses connected to clamped walls by springs on either side, while thin bars behave as a periodic bi-atomic chain of masses connected by springs. The discrete models give an accurate description of the location and width of the stop band in the case of the DCR and the first two dispersion bands for the periodic thin bars. We then combine our asymptotic formulae for arrays of DCR and thin-bars to design a composite structure that displays a negative refraction effect and has a negative phase velocity in a frequency band, and thus behaves in many ways as a negative refractive acoustic medium (NRAM). Finite element computations show that at this frequency, a slab of such NRAM works as a phononic flat superlens whereas two corners of such NRAM sharing a vertex act as an open resonator and can be used to confine sound to a certain extent.


Optical metamaterials
In 1967, Veselago wrote a visionary paper in which materials with simultaneously negative permittivity (ε) and magnetic permeability (µ) were shown to have a negative refractive index [1]. It was shown by a ray analysis that a slab of such a negative refractive index material (NIM) can act as a flat lens that imaged a source on one side to a point on the other. But this result remained as an academic curiosity for almost thirty years, until Pendry et al [2,3] proposed designs of structured materials which would have effectively negative ε and µ. The experimental demonstration of a negative refractive index at GHz frequencies [4] provided a fillip to research in this area (see [5,6] for recent reviews). One should note that these materials are structured at subwavelength lengthscales (typically one tenth of the wavelength) and it is possible to regard them as homogeneous and describe their response by dispersive effective medium parameters (see [7] for a comprehensive survey on homogenization theory). Potential applications of negative refraction came when Pendry [8] demonstrated that the Veselago slab lens not only involves the propagation waves but also the evanescent near-field components of a source in the image formation. It was also shown that two corners of NIM combined in a checkerboard fashion can act as a unique resonator [9]- [12].

Acoustic metamaterials
In 2000, Liu et al [13] provided the first numerical and experimental evidence of localized resonant structures for elastic waves in three-dimensional (3D) arrays of thin coated spheres. This work paved the way towards acoustic analogous of electromagnetic metamaterials. Movchan and Guenneau [14] subsequently proposed to use arrays of cylinders with a split ring cross-section as building blocks for 2D localized resonant structures. Li and Chan [15] independently proposed a similar type of negative acoustic metamaterial. In a recent work, Fang et al [16] experimentally demonstrated a dynamic effective negative stiffness for a chain double 3 'c' resonator (DCR) for ultrasonic waves. Milton et al [17] provided a thorough mathematical frame for such effects including cloaking. Interestingly, similar localization effects have also been observed experimentally by Russell et al [18] at megahertz frequencies in sonic band gap crystal fiber preforms with defects.
In this paper, we improve the design of acoustic metamaterials we proposed in [14] by using a new design that combines the cylindrical DCR with thin stiff sheets. We further show that sound displays a negative refraction effect at the interface of this metamaterial with a homogeneous medium and also has oppositely oriented group and phase velocities in the sense of a negative refractive acoustic medium (NRAM). Further, a flat slab of this metamaterial can focus light in the manner of the Veselago lens [1] and the image of the source has sub-wavelength resolution [8]. We also show that sound can be confined by putting together two corners of the NRAM in a checkerboard fashion.

Set-up of the spectral problem
As a preamble to the phenomenological analysis of the problem, we shall recall briefly its finite element set-up. Thanks to the cylindrical geometry of the problem, elastodynamic waves (governed by the vector Navier equations), decouple into anti-plane shear waves and in-plane coupled pressure and shear waves [19]. In this paper, we focus on the analysis of anti-plane shear waves. For this, we invoke the harmonic acoustic wave equation with derivatives taken in generalized sense where ρ, µ and U are functions of spatial variables satisfying the usual conditions for selfadjointness, boundedness and ellipticity of the wave operator, i.e. there exist two real positive constants m and M such that M µ, ρ m > 0. Due to the periodicity of the lattice, it is customary to require that the eigenfunctions be of the Floquet-Bloch type, i.e.
where k x and k y are components of the so-called Bloch vector k which lies within the Brillouin zone Y * = [0, π] 2 . Solving this spectral problem amounts to looking for the countable set of real positive eigenvalues ω and associated non-zero eigensolutions U of finite energy within the periodic cell Y = [0; 1] 2 for a fixed Bloch vector k within reciprocal space.
The implementation within the finite element package is fairly straightforward. We first multiply equation (1) by a smooth test function V and using Green's formula, we obtain the so-called weak form of the acoustic equation: This expression is then discretized using test functions V taking values on nodes of a triangular mesh of the basic cell Y . From equation (3), we note that setting traction free boundary conditions on an inclusion amounts to assuming Neumann (natural) homogeneous data, whereas ideal contact conditions are such that the quantities U and µ∂U/∂n are preserved across the interfaces.

4
To enforce Floquet-Bloch conditions in equation (3), it is enough to link values of V on opposite sides of the basic cell Y [20]. The finite element formulation was implemented in FEMLAB (COMSOL Multiphysics package), which requires some changes in the readily available boundary conditions. We now wish to discuss this in conjunction with our main numerical results on band diagrams.

Main numerical results on the band structures of the metamaterials
The numerical results presented in this paper were obtained for a matrix material defined by normalized elastic parameters ρ = 1 (mass density) and µ = 1 (shear modulus or first Lamé coefficient). Let us further assume that the other Lamé coefficient for this material is λ = 2.3 (pressure modulus), even though it does not appear within the governing equation (3) for shear waves. In this case, the elastic material parameters used are those of fused silica, ρ = 2.2 × 10 3 kg m −3 , µ = 16.05 × 10 9 Pa and λ = 31.15 × 10 9 Pa. The elastic material considered is micro-structured. It displays two types of inclusions, some of which are C-shaped voids (which are freely vibrating) and the others are highly stiff thin elongated material of density ρ = µ = 140 (which are clamped). Reinforced fibers such as carbon nanotubes are a good candidate for such a stiff phase since such values for ρ and µ lie within the theoretical bounds for realizable effective elastic material parameters [7]. We note that the wave speed c = √ ρ/µ is conserved throughout interfaces (no impedance mismatch). Note that the Poisson ratio ν(µ, λ) of the corresponding material ought to be within reasonable bounds (therefore λ cannot take any possible value). In the case of fused silica, it is around ν = 0.3. Interestingly, Milton [7] has shown that a Poisson ratio ν can be as close as one wants to −1 (when ν < 0 the material shrinks/expands in one direction when compressed/stretched in the orthogonal direction). Let us now outline our main numerical results. We now wish to compute the band structure diagrams of this elastic meta-material. In order to find corresponding Bloch modes in FEMLAB, some changes have to be performed with respect to readily available boundary conditions (Dirichlet, Neumann or periodic ones). A scalar discrete field U (x, y) on the square cell Y with Floquet-Bloch conditions relates the left and the right sides. The set of nodes is separated in three subsets: the nodes on the left side, i.e. with x = 0, corresponding to the column array of unknown u l , the nodes on the right side i.e. with x = 1, corresponding to the column array of unknown u r , and the internal nodes i.e. with x ∈ (0, 1), corresponding to the column array of unknown u. One has the following structure for the matrix problem (corresponding in fact to natural boundary conditions in FEMLAB i.e. Neumann homogeneous boundary conditions): where A is the square Hermitian matrix of the system and b a column array. The solution to be approximated by FEMLAB is a discrete Bloch function satisfying (2). It is convenient to re-express this Floquet-Bloch condition as U (x, y) = U (x, y)e i(k x x+k y y) , where U is Yperiodic and in particular U (x + 1, y) = U (x, y). This way, one can write that U (1, y) = U (1, y)e i(k x +k y y) = U (0, y)e ik x , so that the relation between the left and the right sides is u r = u l e ik x . Consequently, the set of unknowns can be expressed in FEMLAB in function of 5 the reduced set u and u l thanks to: where I and 0 are identity and null matrices, respectively with suitable dimensions. The finite element equations related to the eliminated nodes have now to be taken into account. Thanks to the periodicity of the structure, an element located on the left of the right side of the basic cell Y corresponds to an element on the left of the left side of the basic cell. Therefore, their contributions (i.e. equations corresponding to u r ) must be added to the equations corresponding to u l with the correct phase factor e −ik x . This amounts to multiplying the matrix system by the Hermitian conjugateP of P.
Finally, the linear system to be solved in FEMLAB is: where it should be noted that the matrix system is still Hermitian. To sum up, the implementation of Floquet-Bloch conditions in FEMLAB amounts to transforming a generalized eigenvalue problem (with natural boundary conditions) Au = λBu into an eigenvalue problem with Floquet-Bloch conditions according toPAPu = λPBPu . We will consider our structures to have an axis of invariance so that effectively our models become 2D. Numerical results for the case of a pair of C-shaped voids (shown in figure 1) which we term as DCR, are reported as a band diagram on figure 2. One can see a band gap in the normalized frequency interval [1.1; 1.4], hence which is below the Bragg regime. The convex shape of the second curve is typical for a low frequency local resonance band that falls within the first Bragg zone (heavy photon bands in electromagnetism which arise usually due to Mie or plasmon resonances). For the case of a stiff inclusion with ideal contact conditions which is shaped as a thin elongated bar, we can see in the right hand panel of figure 4, a band gap within the interval [2.6; 3.3]. The concave curvature of the second dispersion curve means there is multiple scattering at play here. We were able to tune down the frequency of this acoustic gap by taking a high-contrast inclusion (yet preserving the wave number across the interface). If we now combine these two types of inclusions (DCR and a non-intersecting stiff thin bar) within a cell, we observe that the second propagating band on figure 5 falls neatly within the gap of figure 2: this is the result of hybridization of eigenfields associated with the two types of inclusions. Most importantly, the slope of the second curve is negative for directions of Bloch vector sitting within M. Hence, we are in possession of an acoustic metamaterial that should display negative refraction within the low frequency regime.

Phenomenological analysis
These numerical results presented above motivate the following phenomenological discussion.

Modeling of DCR using multi-structures
Firstly, let us try to gain some insight into the mechanics of the problem by using asymptotic methods. Such tools enable us to reduce the partial differential equation (1) to certain ordinary differential equation in a thin domain supplied with appropriate boundary conditions.
We model the DCR as a multi-structure, in a way similar to what was done in [14] (see also [21] for the full mathematical details in the analogous case of transverse electric waves propagating within thin-walled photonic crystal fibers).
We denote by the double-split rings as shown figure 1 to preserve the twofold symmetry. Formally, where a and b are functions of variables x and y unless the rings are circular and is a thin ligament of length l j between the 'ends of the letter C'. Here, εh j the thickness of the jth bridge, with ε a small positive non-dimensional parameter. In our case, we have two thin-bridges (1) ε and (2) ε . To derive the asymptotic expansions, we introduce the scaled variable ξ = y/ε so that 1 , the wave equation (1) takes the rescaled form where the derivatives are taken in classical sense (µ is constant in the thin-bridge). The field U is approximated in the form To leading order we obtain (see (3), (9) and (10)) Hence, Assuming that U (0) is given, we derive that the function U (1) satisfies the following model problem on the scaled cross-section of 1 The condition of solvability for the problem (12) has the form Hence, we have shown that to the leading order we can approximate the field U within the thin bridge ( j) ε by the function U (0) which satisfies the wave equation in one-space dimension. We now assume that the field is periodic over the macro-cell since it is localized. This shows that the average of the eigenfield over the macro-cell vanishes. Indeed, let χ 1 denote the value of the field in the large body = { x 2 + y 2 < l j } of the multi-structure (union of the two C-shaped voids) and let χ 2 (which we normalize to 1) denote the value of the field within the complementary area of the macro-cell Y \ excluding the ligaments. Taking V = 1 in (3), we deduce that since U is periodic on ∂Y (and the normal anti-periodic) and U is traction free on ∂ . This shows that the average of the displacement field U over Y vanishes, hence by neglecting the area of the thin bridges, we obtain Since, we have two thin bridges, we have two separate eigensolutions V j , j = 1, 2, corresponding to the vibrations of the thin domains ( j) where εh j and l j are the thickness and the length of the thin ligament ( j) ε , and M is the mass of the body . The bridges are both connected to , hence V 1 (l 1 ) = V 2 (l 2 ) = V , where V is the anti-plane displacement of the rigid body . We note that V j (0) is equal to a non-zero constant unlike in [14] (in that case we assumed that χ 2 = 0, i.e. U = 0 where the bridges meet the region outside ). Here, the constant is chosen in such a way that the average of the field over the basic cell vanishes, as should be expected for a localized (stationary) field. We can see a phononic band gap in the range of normalized frequencies ωd/c ∈ [1.14, 1.29], which is associated with the resonance of the split ring. In that case the eigenfield is axi-symmetric and the inclusion moves as a rigid body (as in figure 6). The lower and upper-edge eigenfrequencies of the gap ω 1 d/c and ω 2 d/c are both in excellent agreement with two asymptotic estimates (23) and (22). The associated eigenfield is displayed in figure 6(a).
The solution of the problem (16)-(18) has the form where c = √ µ/ρ and the frequency ω is given as the solution of the following equation where we invoked Newton's second law that states that total shear force is equal to mass times shear wave acceleration. Looking at a first low frequency, we deduce an explicit asymptotic approximation This estimate actually holds for the frequency of the upper edge of the first phononic band gap of figure 2. It matches that of a LC resonant circuit as depicted in figure 1(d).
We report in figure 2 finite element computations for a periodic cell of side length d with a double C-shaped void. The central disk has a radius of 0.3d and the two cuts have the same length 0.1d and a thickness 0.025d. Therefore, the frequency estimate is which is in excellent agreement with the finite element value ω 1 d/c = 1.1507 for the lower-edge of the first gap (occurring at the K point when k = (π, π)). We notice the presence of a second phononic band gap in the range of normalized frequencies [3.2, 4.1] which occurs due to multiple scattering between inclusions within the array (Bragg regime). We do not know of any asymptotic work showing how to estimate accurately the location or the width of this gap. Interestingly, this gap remains unchanged if one replaces the DCR by circular voids whose radius matches their outer boundary (r = 0.4d), hence one way forward might be to use the multipole method.

Modeling of a periodic array of thin rigid bars using a chain of masses connected by springs
Secondly, we analyze acoustic waves propagating within a doubly periodic array of elongated thin stiff bars (see figure 3(c)). We first implement the weak formulation (3) in the finite element package FEMLAB with normalized material parameters ρ = µ = 1 in the matrix material and ρ = µ = 140 in the thin elongated vertical stiff bar of length 0.95d and thickness 0.01d, where d is the size of the square periodic cell. The corresponding band diagram shown in figure 4 exhibits a phononic band gap in the range of normalized frequencies [1.9, 3.5]. For waves propagating in the M-direction (horizontal x-axis), we can actually approximate the acoustic and optical bands of figure 4. Before undertaking this phenomenological analysis, it is, however, important to note that this gap does overlap with the second dispersion curve and the second gap for an array of DCR inclusions, as depicted in figure 2. This is the principle underlying the effect of negative refraction when combining both types of inclusions.
For waves propagating in the M-direction, let us assume the bars to be infinite along the vertical direction. The problem therefore reduces to the 1D Helmholtz equation with µ j the shear modulus and ρ j the material density in layer S n j = n(a + j (b − a)). This equation is supplied with ideal contact conditions together with Floquet-Bloch conditions where k ∈ [ − π, π] is the Bloch parameter.
To preserve the wave velocity (hence the impedance) in every layer, we now assume that ρ j = µ j . For long wavelengths, if we further assume that ε = µ 1 /µ 2 1, we obtain the approximate dispersion equation (see also [19] for a derivation for thin layers)  where k j is the wave number k j = ω ρ j /µ j and 12 (28) and The above expression bears some resemblance with the dispersion relation accounting for oscillations of a bi-atomic chain of masses m 1 and m 2 connected by massless springs of elasticity constant (stiffness) C (see [22], chapter 4) Indeed, the equations (27) and (30) are identical when The dispersion curves deduced from (27) are displayed in figure 3. We used the geometric parameters a = 0.01, b = 0.99, d = 1 and the small parameter ε = 1/140.

Negative refraction: combining the DCR and the thin stiff elongated bars
Thirdly, we consider Floquet-Bloch waves propagating within a square array of thin vertical bars combined with DCR inclusions. Within a basic cell, the bar of thickness 0.01d (in fact it is a thin sheet, since the structure is invariant in the third direction) is located at a distance 0.04d away from the right-hand side of the cell and a distance 0.05d away from the DCR as shown in figure 7. The corresponding band diagram is shown in figure 5. The second dispersion curve plotted in red color is the result of hybridization of the gap for thin bars of figure 4 and the second dispersion curve in figure 2 (associated with a local resonance of the DCR). Its slope is negative in both K -and M-directions, indicating that the wave propagates with a negative group velocity. It is interesting to notice in figure 7 that the displacement field associated with the local resonance of the DCR at ωd/c = 1.19 is not affected by the presence of the thin bar (the corresponding eigenfrequency is nearly unchanged compared with a DCR taken on its own). The lower acoustic band of figure 6 looks actually identical to that of figure 2. On the other hand, displacement field at the frequency associated with the intersection of the free sound wave (blue line) and the red curve in figure 7 clearly shows some interaction with the DCR, although it is maximum near the location of the thin bar, and hence seems to be driven by the periodic array of stiff thin bars. This illustrates the hybridization of the stop band and the second dispersion curve in figure 4. We emphasize that the eigenfields displayed in figures 7(a) and (b) are of very different nature: one is localized whereas the other one is propagating. This is the cornerstone of our route to negative refraction, which can be actually transposed to the electromagnetic context. Importantly, the red curve of figure 5 (and 6(a)) is nearly flat in the M-direction and it is located within the low-frequency stop band for an array of DCR (see figure 2). Indeed, it has been shown in the electromagnetic context that the nearly flat dispersion on a PC band diagram leads to a better coupling between a wave propagating within the PC and a wave incident on to this composite structure [23]. Also, the fact that the negative slope occurs below the Bragg regime, unlike for PC, could enable some subwavelength resolution, as we shall check in the sequel. To obtain such a focusing effect, the first requirement is to meet the conditions for negative refraction. To check this, we considered a Gaussian beam incident from the left on a finite slab of our NRAM, as reported in figure 6(b). The Gaussian beam makes an angle of 35 • with the vertical axis, centered on the sixth stiff bar and is three cells in width. When it reaches the slab on NRAM, it gets reflected according to the Snell-Descartes' law (the incident angle being visibly equal to the reflected angle). But the beam is transmitted in the NRAM according to inverted Snell-Descartes' law (both incident and transmitted beams are located on the left-hand side of the vertical axis). On leaving the NRAM, the transmitted Gaussian beam is once again negatively refracted. Importantly, it should be noted in figure 6(b) that the barycenters of incident and doubly refracted beams are indeed parallel to each other. This suggests that the heterogeneous slab of NRAM behaves to a certain extent as if it were made of an effective material displaying a negative refractive index close to −1. Nevertheless, our NRAM also exhibits some heterogeneous features at the working wavelength, due to the high material contrast between silica and thin-stiff bars. It is actually hopeless to homogenize such a material in a classical sense. To see this, it is enough to divide by ρ all throughout in (1). We then notice that in the stiff phase the density ρ is a small parameter sitting in front of the higher order derivatives in the equation. Solutions to this type of problem typically involve large derivatives and possibly highly oscillating fields inside the thin stiff inclusions and in its immediate neighborhoods in the silica matrix. It is known that the occurrence of such large derivatives do not preclude a satisfactory homogenization treatment of such singularly perturbed problems [19]. We note here that the phase velocity for waves moving along the -M-direction in this metamaterial i.e.
√ ρ eff /µ eff is complex with opposite signs for the real and imaginary parts.
Since the constituent materials of the metamaterial are lossless, we will only have radiation losses [24] or dissipation that comes in here through the numerical methods used for computing (some dissipation is usually required in most numerical algorithms for stability by damping out fast growing spurious modes). Hence, we expect the imaginary parts of ρ eff and µ eff to contribute to a Im(k) > 0 (which will be reasonably small for propagating modes). Hence, we conclude that Re(k) < 0 unambiguously. On grounds of causality [25], we suggest that the acoustic metamaterial considered here would have Re(ρ eff ) < 0 and Re(µ eff ) < 0. A rigorous parameter retrieval procedure on the lines of those developed for electromagnetic metamaterials [26]- [28] will be required to be implemented on this system to obtain the ρ eff and µ eff . But since the metamaterial structures in our case are not far subwavelength in size at the operating frequency (λ ∼ 2.15d), such a homogenization of all properties via effective medium parameters is difficult. Many parameter retrieval procedures rely on proper choice of the branch for multivalued functions and are difficult to implement once homogenization begins to breakdown and the number of branches increases rapidly. However, we satisfy the requirements for the negative refraction effect to make a flat lens and an open resonator as we will describe now.
To meet the next condition for superlensing described by Pendry in [8], it is necessary to obtain a negative effective refractive index as close as −1 as possible (some imaginary part will nevertheless be present, which is a well known fact in the effective medium theory [24]). This will occur for the normalized frequency of 1.19 which corresponds to a point at the intersection between the acoustic band and the free wave sound line as shown in figure 7: such a point meets the required criteria of both a negative group velocity (negative slope) and a minimized impedance mismatch between the metamaterial and surrounding material (here silica). The inset also shows that the 1.19 isofrequency contour of the red line is nearly that of a free wave (blue circle surrounding the red ellipse as it should): the metamaterial indeed behaves as an effective medium with |n eff | 1 at this frequency (within acceptable margins for anisotropic behavior).
The fact that n eff close to −1 is confirmed by the numerical experiment reported in figure 8(a) where a plane wave of that frequency is incident upon a slab of this metamaterial: it appears that the structure is invisible at this frequency thanks to impedance matching between the surrounding medium (silica) and the metamaterial (up to the unavoidable absorptive nature of the composite slab). Such a phenomenon was experimentally observed in few instances in flat super lenses in electromagnetism [29].
We can now proceed to investigate the focusing effect by a flat slab of the metamaterial. For this, we look at the eigenfunctions associated with a finite structure consisting of 16 square cells (slab lens) combining DCR and thin bars, embedded in an infinite outer medium made of silica. We model the scattering problem using the weak form of the acoustic wave equation (3) where the spectral parameter (frequency term) is now replaced by a forcing term (a line source). Perfectly matched layers (PML) provide a reflectionless interface between the region of interest (large middle square containing the metamaterial in figure 8) and the PML (four rectangles and four small squares) at all incident angles. PML were originally introduced by Berenger [30] in electromagnetism, in which case the regions of PML consist of anisotropic absorptive media. In our case, we consider some orthotropic absorptive media defined by where T is a representation of the metric (rank two) tensor with non-zero entries in Cartesian coordinates given by Here, s x , s y and s z represent complex stretched coordinates defined by in the case of a wave propagating in the x-direction. The real parameters a and b are chosen empirically so that the amplitude of the outgoing acoustic wave vanishes smoothly within the PML region which was meshed with a mesh size of about a tenth of the wavelength. In our numerical simulations reported in figure 8, we took the normalized elastic parameters for silica ρ = µ = 1 and λ = 2.3 in (32). We further chose a = b = 1 for the working normalized frequency ωd √ ρ/µ = 1.19, with d = 0.2 being the cell's size of our metamaterial's building block. The thickness of the PML regions was taken as large as 5d to ensure convergence towards an accurate numerical solution. In figure 8(c), we plot the displacement field (shear wave) emitted by a line source at normalized frequency ωd/c = 1.19 which is located one cell below the slab lens. We can clearly see an image located one cell above the slab lens, in accordance with the ray analysis originally proposed by Veselago for a homogeneous flat lens whose refractive index is −1 [1]. An interference pattern due to reflected waves arising from the imperfect impedance matching for waves at larger oblique incidence can be see in the field map. In order to show the potential of this metamaterial for the super focusing effect as discussed by Pendry [8], we plot the modulus of the displacement field in the source and image planes, see figure 8(c). We can see in figure 8(e) that the full-width at half-maximum (FWHM) of the function is significantly smaller than half the wavelength with the working wavelength which is actually 2.15d where d is the cell's size. Note that we do not have an infinite image resolution. That is due to the presence of the unavoidable effective dissipation and dispersion that are inherent in any numerical calculation [31] and any such discrepancy from the ideal perfect lens conditions of n = −1 is well known to severely limit the extent of subwavelength image resolution [32]. Indeed, in figure 8(b), the domain for finite element computations excluding PML is a square of computational size 2 and the cell's size is d = 0.2. In figure 8(e), we plot the modulus of the profile of the displacement field along the x-direction. It shows transverse oscillations associated with the line source whose location is at a computational distance 1: this represents 1/0.43 = 2.33 wavelengths or 5 unit cells. The measure of the distance in terms of the wavelength are actually added on top of the horizontal axis for convenience (note that the wavelength corresponds to the distance between three successive maxima-minima, as we plot the modulus of the eigenfield). We note the asymmetry of the profile which we attribute to the anisotropic nature of the effective behavior of the metamaterial at the working frequency (the isofrequency contour in the inset of figure 7 is not perfectly circular). We checked that adding a thin stiff bar on the left-hand side of the basic cell in a symmetric fashion makes the profile more symmetric.
We conclude that sonic surface modes possibly enhance the evanescent modes and contribute to subwavelength image resolution similar to the surface plasmon polariton excitations in electromagnetism [16].
Lastly, we will demonstrate sound confinement through negative refraction in corner reflectors, following the original idea of Notomi [11] for PC in electromagnetism. As can be seen from the ray analysis in figure 8(b), sound trajectories circle around the point of intersection of the two corners in closed loops. We first checked that some resonances occurred for a corner reflector surrounded by an infinite outer medium (silica) which we modeled using the spectral formulation (3) together with PML. We observed that the larger the number of cells, the smaller the imaginary part of the resonant eigenfrequency: eigenvalues are complex with an imaginary part accounting for the leakage of the corresponding eigenmode in a non-dissipative medium.