Subwavelength negative refraction and flexural wave lens design via resonant double-negative piezoelectric metamaterial

We report the concept and demonstration of a double-negative, resonant metamaterial characterized by both dynamic negative mass and stiffness for negative refraction of flexural wave modes by means of a lens designed using this concept. The negative equivalent material properties are obtained in the subwavelength regime by concurrently exploiting both the effect of mechanical resonators (negative mass) and of piezoelectric patches with inductive resonant shunts (negative stiffness), leading to double-negative behavior. Following the theoretical foundations based on a modal framework, we analytically derive the frequency-dependent mass and stiffness properties as a function of the electromechanical parameters. The findings are corroborated by numerical computation of dispersion properties and simulations showing the focusing of a point source. As a case study, energy harvesting performance enhancement by exploiting the piezoelectric effect at the focal spot is also discussed.


Introduction: file preparation and submission
The early work of Veselago (1968) [1] is commonly recognized to be the first suggesting the pioneering idea that a material exhibiting at the same time both negative dielectric permittivity and negative magnetic permeability leads to a negative refractive index.The actual implementation of such an idea was demonstrated much later (2000) [2][3][4][5][6], with the advent of the concept of 'metamaterial', when Pendry and coworkers showed that the required double-negative effective properties could be achieved by combining split ring resonators and periodic arrangements of wires.In his seminal paper [3], Pendry also pointed out that a flat lens made by a slab of such material could overcome the diffraction limit, i.e. generate images with details smaller than the wavelength as a 'superlens'.
Negative refraction has also been sought for elastic waves in natural crystals exhibiting strong anisotropy in the slowness surface [7,8], such that group velocity and wavevectors are not aligned.In a similar fashion, artificial photonic [9][10][11][12][13] and phononic crystals [8,14] have been designed to exhibit negative refraction, exploiting negative group velocities arising in second branches near the Γ point in reciprocal space, or folding of the first branch near the boundaries of the first Brillouin zone.For example, lattices of steel rods [15,16] or coated tungsten bars [17] in water have been used to experimentally show negative refraction of ultrasound waves in water.Norris and coworkers [18] showed that a bi-mode material designed to be transparent in water in the long wavelength regime, actually exhibits an optical branch that can be used to obtain negative refraction too.
Other than acoustic waves in fluids, Morvan et al [19] showed that a phononic crystal made by cylindrical air cavities in a solid matrix can be used for transverse elastic waves, while steel rods arrangements in epoxy [20] were used to demonstrate negative refraction for longitudinal elastic waves.More than that, there have also been efforts to design phononic crystals whose dispersion relation unlocks the possibility to achieve flat lensing for elastic waves in waveguides where energy is confined between two free surfaces, i.e. for shear horizontal [21] and antisymmetric (A) [22][23][24][25] Lamb wave modes.
Much less work has been done on focusing phonons following the original idea of Veselago, i.e. implementing a metamaterial exhibiting both negative mass and modulus.For example, Oh et al [26] showed that a spring mass network can be finely designed such that isotropic negative refraction index can be achieved through double-negative equivalent properties for in-plane motion, to then use such scheme to design a unit cell made by beams and blocks in an aluminum plate to focus the S 0 Lamb mode.Zhu et al [27] used simultaneous in-plane rotational and translational resonances observed in a chiral metamaterial to show a similar result.
In this paper, we follow another path to obtain doublenegative properties.It is long known that the inertial contribution of mechanical resonators near the resonance can be understood as an equivalent negative dynamic mass [28].We exploit this effect on the lowest antisymmetric Lamb mode of a plate by covering its surface with a periodic arrangement of mechanical harmonic oscillators.We combine this effect with that of an equivalent negative stiffness which is obtained by bonding piezoelectric elements shunted to inductive resonant circuits, thereby exploiting an electromechanical resonance.The use of piezo patches with shunting circuits has a long history of applications starting from vibration control [29] and energy harvesting [30] to the more recent applications as in resonant metastructures for bandgap formation [31][32][33][34] and wave steering [35,36] or for tunable topologically protected edge waves [37] graded index focusing [38] time modulation of stiffness properties [39] and rainbow trapping [40].Sugino et al [33] first underlined that merging mechanical and electromechanical bandgaps in a 1D beam resulted in the annihilation of the two stop bands while later Gao and Wang [41] used mechanical resonators and piezoelectric patches shunted with a non-resonant synthetic negative capacitance circuit to obtain a doubly negative metamaterial.Nonetheless, the electromechanical resonators have not yet been exploited in conjunction with mechanical ones in a way that their resonance frequency is finely tuned such that both negative properties occur in the same frequency range, resulting in a negative slope pass-band in the dispersion relation.One advantage of this method of obtaining negative refraction flat lensing is that the frequency associated with the resonators can be moved to low values in order to work in the subwavelength regime, thus obtaining a very compact device.Moreover, the conversion between mechanical and electrical energy that naturally occurs in the piezoelectric elements, opens the chance to harvest vibration energy, especially exploiting the amplification that is observed around the image that forms inside the lens.
In the next section, a recent electromechanical framework [34] is here adopted to show the individual contribution of both types of resonators in creating the homogenized doublenegative material properties making use of Kirchhoff's theory of flexural waves in thin plates, to which the lowest antisymmetric Lamb mode converges when dealing with thin structures at low frequencies.The aforementioned method has been shown to provide reliable estimates of the frequency range where the double-negative properties are expected [32,34,42] and validated against both three-dimensional (3D) finite element models [34], the plane wave expansion method [33], and experiments [42] (for this reason, being the focus of the work restricted to thin structures, we will loosely use the terms 'flexural wave' both in the context of the thin plate theory and when dealing with 3D finite elements waveguides).The numerical dispersion relation is thus computed and associated modes are analyzed, before assessing the performance of the flat lens with a case study on energy harvesting.Finally, conclusions are drawn regarding this new class of negative refraction metamaterials.

Governing equations and equivalent double-negative properties
Figure 1 illustrates the layout of the metastructure adopted to obtain the double-negative effective material properties: a substrate metallic plate is sandwiched between two piezoelectric plates of the same shape and thickness h p = h s .Both piezoelectric elements are poled in the same direction, with the poling axis oriented along the thickness of the laminate.The aluminum substrate works as an electrode for the inner faces of the piezoelectric patches, which are grounded and shorted to the outer ones, except for the areas covered by a hexagonal pattern of thin circular metallic electrodes, each one used to create an electric parallel connection to an inductor, as shown in the inset in figure 1.At each lattice point, a mechanical resonator is also added.According to Kirchhoff's theory of in thin plates, the partial differential equation governing the outof-plane flexural motion w(P, t) of the sandwich reads: where v(P, t) is the voltage across the electrodes, J(P, t) is the density of current, and P is the position vector.The source distribution f(P, t) will be used to take into account the inertial effect of the resonators.The other coefficients appearing in the equations can be computed from the physical and geometrical properties of the plates adopted: 11 , ē31 and εS 33 are the density, the equivalent short circuit elastic constant, the equivalent piezoelectric coupling coefficient and the equivalent strain free dielectric permittivity of the piezo plates [34].The voltage distribution can be approximated as [34]: where d j (P) is the indicator function for the surface of the electrodes, i.e.
with D j indicating the surface of the jth electrode.Defining as C p,j the capacitance associated to the same jth electrode: where ∆D j is the area of the jth electrode, and defining ˆDj J (P, t) the current associated with each electrode (that depends upon v j through the operator Y j ), by integrating the electrical equation (1b) over the area of each electrode, one obtains an equation in the form: for each electrode.To include mechanical resonators, consider the force f(P, t) to be the sum of an external force f 1 (P, t) and the force exchanged with the resonators: where δ(P) is the Dirac delta distribution, P j the position of the jth resonator, k is the stiffness of the resonator, and u j (t) is the relative displacement of the resonator with respect to the laminate, governed by: Substituting equation ( 7) into (1a) and using ( 6) and ( 8), we obtain: Adopting the modal decomposition w(P, t) = ∑ N r=1 ϕ r (P)η r (t), with ϕ r (P) the eigenmodes of the plate computed for short circuit of each electrode pair (v j = 0): Multiplying the first equation by a given mode, integrating over the whole structure, and exploiting the orthogonality conditions: being δ rs the Kronecker delta, and ω r the natural circular frequency associated to the rth mode, one obtains N equations in the form: ϕr (P j ) mpϕ k (P j ) ηr (t) = ˆD f 1 ϕr (P) dD (12) where m j = µ∆D j m p has been used, with µ = ∑ s j =1 m j /(m p ∆D) i.e. the ratio between the total mass of the resonators over the total mass of the beam.Assuming a harmonic forcing q r = ´f1 ϕ r dD = Q r e jΩt , a harmonic motion of the type: is obtained.The system can thus be rewritten as: where Γ r,j = ˆdj (P) ∇ 2 ϕ r (P) dD (15) and Y j represents the electrical admittance of the shunt circuit.
Expressing the voltage from equation (14b) the following is found: In the same way, the equation ( 14c) is used to express the displacement of the resonator: where ω 2 m = k j /m j is the natural circular frequency of the mechanical resonator.Substituting equations ( 16) and ( 17) into (14a), we get: Using and, considering that for the average theorem there exists a point P rj such that: Equation ( 18) can be rewritten as: ϕr (P j ) mpϕ k (P j ) ∆D j H k = 0. ( Considering now to increase the number of resonators to infinite, while taking the size of the structure to be fixed, i.e. shrinking the size of each resonator to zero: the equation finally becomes (24) When a purely inductive shunt is used: with ω e the resonance frequency of the electromechanical resonator, we finally write: with It is thus shown that the mass becomes negative when and the stiffness becomes negative when The two different types of resonators (mechanical and electromechanical) can be thus adjusted to obtain double-negative effective material properties.

Hybrid electromechanical unit cell design
The commercial software COMSOL Multiphysics is used to compute via the finite element method the dispersion band diagram of an infinite repetition of unit cells as the one depicted in figure 1, in order to verify the behavior of an infinite structure designed by the use of the formulas derived in the previous.To do so, a full 3D model of the substrate and the piezoelectric patches composing a unit cell is adopted, while both the shunt circuit and the mechanical resonator are modeled by means of lumped parameters systems, with the spring-mass resonator attached to the mesh node at the center of the top surface of the unit cell.Opposite boundaries of the hexagonal cell are linked with periodic boundary conditions, according to Bloch theorem, looking for the dispersion relation of elastic waves confined between two free surfaces.The resulting eigenvalue problem is solved for wavenumbers along the irreducible Brillouin zone, which is depicted in figure 2(a).The material considered for the substrate plate is aluminum, while PZT-5 H is adopted as piezoelectric material.Physical and geometrical parameters adopted in the computations are listed in table 1.Note that practical implementations of the electrical circuit can benefit from synthetic impedance circuits [43,44], so that analog inductors are not required.The resulting band structure is reported in figure 2(b).A bandgap is created at the edges of the Brillouin zone due to Bragg scattering at approximately 1 kHz for the flexural mode, i.e. the one characterized by quadratic dispersion.Below that frequency, a band showing a negative slope can be inspected just above 300 Hz.Indeed, with the parameter's choice as reported in table 1, the dynamic mass is predicted to be negative in the range 300 < f < 410 Hz, and the dynamic stiffness in the range 294 < f < 350 Hz. Figure 3 shows a detail of the dispersion relation along the ΓM direction, along with the computed mode belonging to the first and second branches.Although the computational domain is restricted to one unit cell only, the mode can be plotted on an assembly of cells by exploiting the nature of the Bloch solution, i.e. the ratio of the displacement vector computed on points in space that are separated by a linear combination of lattice vectors a i b i , a i ∈ Z is equal to e κ•(aibi) , κ being the wavevector.Gray arrows in figure 3 represent the wavevector, while superimposed on the distribution of the out-of-plane displacements, black arrows indicate the component of the Poynting vector along the direction of the wavevector.As it can be seen from the figure, the first branch is characterized by the typical long wavelength flexural mode, with Poynting vector and wavevector κ pointing in the same direction.In contrast, the negative slope branch shows a backward antisymmetric mode with black arrows pointing in the opposite direction of κ, witnessing that energy flows in the opposite direction.The displacement of the center of each cell moving counter phase with respect to the others witnesses the out-of-phase relative motion of the mechanical resonators with respect to the plate.

Design of the flat lens with negative refraction
Next, we aim at designing a flat lens that exploits the doublenegative material properties as introduced in the previous sections.First of all, in figure 4(a) the dispersion relation of the metamaterial is superimposed to that of the same piezosubstrate-piezo laminate, but with removed mechanical resonators and removed electrodes and circuit.At 321 Hz the flexural mode of the homogeneous laminate crosses the second branch of the double-negative metamaterial: at that frequency, they both show the same phase speed but opposite group velocity, thus enabling negative refraction.To obtain all-angle negative refraction, the isofrequency contour at the coupling frequency should be as isotropic as possible.The hexagonal lattice has indeed been chosen to exploit the maximum available crystal symmetry.which almost overlaps with the circular one of the homogeneous laminate.
A finite-element (COMSOL) model of a flat lens surrounded by homogeneous regions is thus built as illustrated in the scheme in figure 5(a) and a point force is applied acting perpendicularly to the plate as a point source for flexural waves.Note that the overall structure is simply made by an aluminum plate covered uniformly above and below by two equally sized piezo patches, and the lens is just obtained by attaching in the central region a finite lattice of 8 × 20 mechanical resonators and shunted segmented electrodes.The number of unit cells in the propagation direction was chosen to be sufficiently high for the effects of periodicity to be observed as per prior work [31,42,45], where the effect of the number of mechanical/electromechanical resonators on bandgap formation in locally resonant metamaterials is explored.Based on such efforts, eight resonators were deemed sufficient in the propagation direction.As for the perpendicular direction, it governs the lens aperture, which can affect aspects like focal spot size, as further discussed below.The computational domain is surrounded on all sides by a perfectly matched layer, in order to suppress scattering thus mimicking an unbounded domain.The obtained out-of-plane displacement field is shown in figures 5(b)-(e), where both the image inside the lens and in the homogeneous regions can be observed.The maximum displacement in the image outside the lens is measured at a distance from the source equal to 1.86 times the size of the lens.The full width at half maximum evaluated along a vertical line passing across the maximum is equal to 0.59λ, where λ is the wavelength at 321 Hz, as computed from the dispersion relation.These discrepancies from the theoretically expected results (i.e. the distance between image and source being different from two times the size of the lens and the size of the image being greater than half the wavelength) can be attributed to (i) the isofrequency contours not being exactly circular, with waves traveling along the ΓK direction being faster than those traveling along the ΓM one, and thus not focalizing all in the same point, (ii) the finite size of the lens in the direction perpendicular to propagation.Another matter of concern regards the fact that there are modes in the lens other than the backward one at the selected frequency, though their coupling with the wave in the exterior plate is expected to be very low. Figure 5(c) shows the amplitude of the out-of-plane vibrations, while figure 5(d) shows the absolute value of the electric potential on the outer surface of the piezo patch.Note that in terms of out-of-plane displacement, the image inside the lens is less strong than the image outside it.This is related to the fact that part of the energy is converted to the electrical domain, as shown in figure 5(d), where the value of the electric potential shows a strong maximum near the image.Note that the distribution of electric potential shows high values inside the 'cone' whose aperture is defined by the set of rays that impinge on the right-hand side of the lens, i.e. the rays that are focalized in the image outside the lens.
Figure 6 shows the sensitivity of the solution to various levels of structural damping.Absorption in the simulation is introduced via the use of an isotropic loss factor, both in the aluminum substrate (η al ) and in the piezoelectric patches (η p ).It can be seen how damping in the substrate has little to no effect on the solution, since increasing by a factor of 10 η al has no appreciable impact on the focusing capabilities of the lens.Conversely, the solution is more sensitive to damping in the piezo patches.Nonetheless, it has been shown [43] that there are ways to reduce absorption in piezo patches by adding negative damping via synthetic shunt circuits.This means that while the behavior is sensitive to losses, this could be still compensated in cases where losses would be too high.

Case study on energy harvesting
The focusing of energy inside the lens exhibited in the electrical field in figure 5(d) suggests that part of that energy could be harvested as a demonstrative case study.To check the performance of the system with respect to energy harvesting, a resistor is placed in parallel to each inductor in the shunt circuit, as a representative load that is powered by the lens.The shunting circuit can now be characterized by the dimensionless number τ ω e where τ = RC p,j is the time constant.Figure 7 shows the out-of-plane displacements and the electric power P = |V| 2 /R collected by each electrode, for various values of τ ω e , obtained varying R from 275 Ω to 27.5 MΩ.As it can be seen, for low values of the resistor, the inductor is short-circuited, and the electrical resonance basically does not happen.As a result, the assembly of unit cells does not behave anymore as a double-negative metamaterial, but as an elastic metamaterial with a locally resonant bandgap, resulting from the negative mass associated with mechanical resonators.As a result, for low values of the resistance the lens behaves as a mirror instead, and the maximum power is collected near the interface with the background plate in the first row of unit cells, where the amplitude of the evanescent waves excited inside the mirror is at its maximum.For the value of τ ω e tending to infinite, the resistor behaves as an open circuit, and the behavior tends to the one depicted in figure 5. Note that when one changes the order of magnitude for τ ω e not only does the distribution of the power inside the lens change but also the maximum harvested power for each electrode changes the order of magnitude.Figure 8 shows the maximum power observed among the electrodes and the overall power harvested obtained as a sum of the power of every single electrode as a function of the τ ω e parameters.A peak maximum is observed for both curves, around a range of resistance in between the two extreme conditions of open and short circuit.

Conclusions
In this work, mechanical and electromechanical resonances have been combined in order to obtain effective doublenegative material properties in the sub-wavelength regime.An analytical method to compute the frequency range in which both the dynamic effective mass and stiffness take negative values has been introduced, and finite element computation of the dispersion relation has been used to validate such formulas.Indeed, a branch characterized by negative group velocity appears in the frequency range in which the bandgap generated by negative mass and negative stiffness overlap.The subwavelength nature of such unit cells is key for obtaining compact flat lenses that exploit negative refraction to focus flexural waves in plates.As a case study, the designed flat lens is used for vibration energy harvesting.

Figure 1 .
Figure 1.Schematic of the unit cell comprising the mechanical and electromechanical resonators.The inset shows the wiring scheme used to connect the piezos to the inductive shunt circuit.

Figure 2 .
Figure 2. (a) Diagram showing a small region with removed electrodes between the circular electrode and the grounded area used in the numerical computations to avoid discontinuities in the electrical quantities [34].The radial thickness is dr = 0.025a, with a the length of the lattice vectors.The inset shows the Brillouin zone and the associated irreducible Brillouin zone.(b) Dispersion diagram along the boundaries of the irreducible Brillouin zone.

Figure 3 .
Figure 3. Detail of the dispersion diagram along the ΓM direction, showing the band characterized by negative group velocity.The out-of-plane displacement associated with the first and second bands is also displayed for the same wavenumber κ (represented in the insets as a gray arrow).Black arrows represent the component of the Poynting vector along the wavenumber direction, showing the direction of the energy flux.
Figure 4(b) shows the second branch of the dispersion relation over the first Brillouin zone between 310 and 330 Hz, along with the isofrequency contour at 321 Hz

Figure 4 .
Figure 4. (a) Detail of the dispersion relation of the double-negative unit cell, superimposed to that of an open aluminum-piezo laminate with no electrodes and no mechanical resonators, showing the coupling between the flexural mode and the branch characterized by negative group velocity.(b) Isofrequency contour (blue) at 321 Hz of the second branch, which almost overlaps with the isotropic flexural mode of the bare aluminum-piezo laminate (red).

Figure 5 .
Figure 5. (a) Layout of the Veselago lens.(b) Out of plane displacement.(c) The absolute value of the out-of-plane displacement.(d) Voltage distribution of the piezo elements.(e) Out-of-plane displacement (3D).

Figure 6 .
Figure 6.Influence of various levels of damping on the solution: η al stands for the loss factor in the substrate, while ηp is the isotropic loss factor used for the piezo patches.

Figure 7 .
Figure 7. Out of plane displacement and electric power fields for various values of τ ωe.

Figure 8 .
Figure 8. Maximum and overall harvested electrical power as a function of the dimensionless load resistance τ ωe.