Beam paths of flexural Lamb waves at high frequency in the first band within phononic crystal-based acoustic lenses

This work deals with an analytical and numerical study of the focusing of the lowest order anti-symmetric Lamb wave in gradient index phononic crystals. Computing the ray trajectories of the elastic beam allowed us to analyze the lateral dimensions and shape of the focus, either in the inner or behind the phononic crystal-based acoustic lenses, for frequencies within a broad range in the first band. We analyzed and discussed the focusing behaviors inside the acoustic lenses where the focalization at sub-wavelength scale was achieved. The focalization behind the gradient index phononic crystal is shown to be efficient as well: we report on FMHM = 0.63λ at 11MHz.


I. INTRODUCTION
The past decade has witnessed the development of focusing of acoustic waves through phononic crystals (PCs). Analogously to heterostructures designed to achieve negative refraction of electromagnetic waves, 1,2 PCs have been conceived in order to focus acoustic waves in a frequency range where the negative refraction is possible. [3][4][5][6][7][8] It has been proved both theoretically and experimentally that such PCs can realize the subwavelength focusing or imaging. [3][4][5][6][7][8] To a certain extent, the resolution of the focusing was beyond the diffraction limit with the evanescent waves being gathered at the focusing through the bound modes. 9,10 These studies of acoustic waves in fluids have inspired a series of works dealing with other polarized waves both on the negative refraction phenomenon [11][12][13] and on the focusing of negative-refractive PCs. [14][15][16][17][18] On the other hand, there has been a great deal of work on systems based onto a positive refraction process to realize the focusing. In the relatively early studies, [19][20][21][22] the PCs were usually composed of an uniform structure with various geometrical shapes, and the same physical limitations as for conventional optical lenses were therefore observed. Later, analogously to the studies in optics, 23 gradient-index (GRIN) PCs were proposed, especially for waves at low frequency, on the first acoustic band. 24,25 To now the focusing achievements for waves with different polarizations have been investigated. [26][27][28][29][30][31][32][33] The flat GRIN PC features a graded sound velocity, making it feasible to control the propagation of waves through processes as focusing, collimating, cloaking and so on. Moreover, GRIN PCs have been designed to focus waves at a subwavelength resolution, generally with a highly gathered energy at the focus which is quite different from the focusing through negative refraction. Indeed, both numerical and experimental observed behaviors of the lowest-order flexural (A 0 ) Lamb wave promise important applications in many fields, as for example nondestructive testing, acoustical imaging, and wave guiding. [31][32][33] However, up to now there is only few works aiming at analyzing the shape, size, and location of the focus area. In this work, we analyze the consequences of the anisotropy along some lines of inclusions in the PCs and we propose an analytical model that allows predicting the achievements of a GRIN PC designed to focus Lamb waves in the MHz range.
a Author to whom correspondence should be addressed.

II. PHYSICS MODEL AND RAY TRAJECTORIES
The GRIN PC is an analogue of the optical GRIN lens however with discrete structures. 25 One of the simplest way to get a graded sound velocity along one direction of the 2D PC, is to modify the filling ratio (ff) along that direction. The GRIN PC we have elaborated for this work was obtained by drilling air holes throughout a silicon plate, in a square lattice. The diameter of the holes, and consequently the filling ratio, varied along the gradient direction. The first band of the A 0 mode shifts gradually as the filling ratio changes, as shown in Fig. 1 for a PC in square lattice featuring the periodic constant a=100µm and the plate thickness h=110µm. 34 To determine the size of air holes, a hyperbolic secant profile was used to define the transverse refraction index (i.e. y axis) of the GRIN PC while the refraction index keeps uniform along x axis as: where δ is the gradient coefficient, k ΓX stands for the wavenumber of the flexural Lamb mode propagating along ΓX in the PC, and k B is the wavenumber for the wave along the same direction in the silicon plate. Due to the anisotropy along certain lines of inclusions in the GIRN PC, the hyperbolic profile may not focus all ray trajectories at a single point, but rather on a focal spot. 34,35 The GRIN PC we designed for this study was optimized to focalize flexural Lamb waves at 5MHz. The refractive index along the row of inclusions at x=0 and the gradient coefficient were n 0 =1.166 and δ=0.059 respectively. The radii of the air holes were decreasing by steps from the central layer to the edges, according to Eq. (1). The calculated discrete filling ratios are given in the label in Fig. 1, together with the sequence of the unit cell along y axis. In this work, we took account of the crystallographic cubic symmetry of silicon, 34 with the x-axis and y-axis being parallel to the crystallographic directions <100>and <010>respectively.
A formalism which takes into account the anisotropy along some lines of inclusions in the PC is applied here to calculate the ray trajectories of A 0 mode, whatever is the frequency in the first band. Actually, while optimized for a flexural mode at 5 MHz, the geometry of the GRIN PC is also efficient at higher frequencies, 34 up to 11 MHz, close to the lower edge of the band gap of the PC with ff=0.503 (see Fig. 1). The formalism described below is based on the consideration of the exact shape of the equal-frequency contours (EFCs) of the PCs. The modulus of k vector is assumed to obey 34 along the y axis. k ΓM is the wavenumber along ΓM in the PC, α is a dimensionless coefficient accounting for the fact that the EFCs may depart from a circle, and θ is the angle between k vector in the first Brillouin zone and ΓX, running in the interval [0, 2π].  Figure 2 shows the EFCs for A 0 mode at different frequencies for a PC with ff=0.503. Solid lines indicate the EFCs calculated by the finite element method (FEM). The stars, crosses and circles are the fitting data deduced from Eq. (2) for 5MHz, 8MHz and 11MHz, respectively. Good agreement was found between the results of both groups. However, the fitting formula Eq. (2) losses its validness when EFC is no more centered at Γ in the first Brillouin zone. This situation is out of the scope of this work although similar procedures can be applied in that case. Details on the calculation of the ray trajectories may be found elsewhere 34 and are not repeated here.

III. NUMERICAL RESULTS
We used FEM to study the wave behaviors inside the GRIN PC featuring 55 identical columns as described above. A line source with monochromatic frequency was placed in front of the GRIN PC with large unstructured silicon areas on either side to avoid the effects of reflected waves from the boundaries and to better mimic a realistic experimental configuration. The out-of-plane displacement u z at any point on the surface of the GRIN PC was then computed as a function of time, and we recorded the maximum value at each measured point.

A. Focusing inside GRIN PC
We show in Figure 3 the result when the frequency is tuned to (a) 5MHz, (b) 8MHz, and (c) 11MHz. The dashed lines are for the corresponding ray trajectories derived from Eq. (2). 34 In Fig. 3(a), a focusing zone appears at x∼[27a, 42a] with the overall maximum of u z located at x=35a. Inside this focusing zone, the ray trajectories are dispersed, because of aberrations which are very similar to the spherical aberrations encountered in optical lenses. The focusing zone much more extends along x axis than it does along y axis, where the lateral dimensions of the spot is y∼[−1.5a, 1.5a]. Furthermore, the ray trajectories issued from the edges converge at a shorter distance than the ray trajectories issued from the zone close to the central axis of the GRIN PC do. For example, the rays impinging the lens at y=±8a converge at a distance x=26a, while the value is 40a for rays at y=±2a. This is related to the anisotropy coefficient α( y) which takes larger values for y close to the central axis than it does for larger values of y, where the filling ratios are small. When the anisotropy within each PC constituting the lens is ignored by setting k ΓM =k ΓX , all rays focalize to a single point, at a distance 27a (π/2δ) from the origin. However, when considering this anisotropy the directions of phase velocity (parallel to k) and group velocity (tangent to the ray trajectory) are not parallel anymore, as shown in the inset in Fig. 4. As a consequence of Snell's law that states the conservation of the transverse wave vector k x , for incident waves close to the central layer the angle θ between x axis and the phase velocity is limited to a very small region where the angle ϕ made by the group velocity is much smaller than θ, as shown by the ratio tan(ϕ)/tan(θ) for the rays incident at y=−2a. Thus the ray trajectory is prolonged along x coordinate, and the focal distance becomes larger than expected. For rays incident far from the central axis angle ϕ is first larger than θ, and then becomes smaller as x coordinate increases, as shown by the ratio tan(ϕ)/tan(θ) for the ray incident at y=−8a. This characteristic may compensate the anisotropy to a certain extent, resulting in a focal distance very close to the value derived from π/2δ. In Fig. 3(b), the overall maximum of u z at 8 MHz was found to be centered at x=30a, inside a focusing zone extending over x=[27a, 32a] and the focal distance becomes shorter. Actually, the  refraction index fits with a hyperbolic secant profile with n 0 =1.129 and δ=0.063. Meanwhile, the ratio tan(ϕ)/tan(θ) keeps a value close to 1, with smaller variations for rays impinging the lens at y=±8a than it was at 5MHz; even, this ratio deviates less from 1 for the ray trajectory with y=±2a. This behavior corresponds to the reduced influence caused by the anisotropy of the GRIN PC. In turn, the focus area is less elongated along x-axis, and the focal distance becomes smaller. Similar characteristics can be observed in Fig. 3(c) for focalization at 11MHz with even a shorter focal distance at x=25a: the refraction index can well be fitted by a hyperbolic secant profile with n 0 =1.151 and the δ=0.068; moreover, the anisotropy is expected to have less influence as can be seen from the almost constant value of the ratio tan(ϕ)/tan(θ) at 11MHz. Fig. 3(d) further shows the focal distance as a function of frequency. This quasi-linear behavior is attributed to the increase of the gradient coefficient δ as the frequency increases. Fig. 5 depicts the transverse profiles of the normalized maximum of u z at the peak value at 5MHz (dashed line), 8MHz (dot dashed line) and 11MHZ (solid line) respectively. The full width at half maximum (FWHM) of these profiles is 3.62a, 2.69a and 2.22a, respectively. Note that measuring all the lower FWHM than the frequency increases is not surprising since the incident waves get better focused for short curvatures of the trajectories. The FWHMs, when normalized to the wavelength in the central layer λ min , are 0.83λ min , 0.83λ min and 0.85λ min , respectively, i.e. less than one wavelength. Indeed, the pattern at the focus depends on several factors, including the anisotropy of the GRIN PC, the gradient coefficient δ, 35 and the ratio between the wavelength and the lattice constant a. These factors may compensate each other, making the normalized resolution more or less the same whatever the frequency we considered in this work was, as shown in the inset in Fig. 5.
The maximum field intensity (MFI), defined as the ratio between the overall maximum value of u z at the focusing and the amplitude of the incident wave, is 3.73, 4.97 and 4.63 at 5, 8, and 11 MHz respectively. We show in Fig. 6(a) the MFI (circle markers) which features an overall increasing trend as the frequency increases, with the maximum value at 9.5 MHz. The reduced focusing zone along y axis, evaluated through the FWHM (in units of a) in Fig. 6(a), results in the increasing MFI. Meanwhile, since the focus is formed in the air/silicon heterostructure, the size of the holes can also, to a certain extent, affect the MFI. This is explained below and by Fig. 7.
The relatively smaller MFI at 11 MHz relates to the large reflective coefficient at the interface between the background medium and the GRIN lens. To illustrate this, we show in Fig. 6(b) the reflective coefficient computed by FEM. 37 The geometry we used is depicted in Fig. 6(c): forces polarized along z direction, uniform in space and time-varying according to a pulsed function f(t), were applied in the form of fringes in front of a PC made of identical air holes in a silicon plate. The distance between two consecutive fringes was set to one wavelength. We defined the function f(t) in such a way that the fringes produced wave pulses of A 0 Lamb mode with a carrier frequency being exactly the wanted one. The normal displacement u z was recorded in front of the PC (dashed line) in a place where both the incident and the reflected waves appear well separated. We obtained that way, the overall reflective coefficient defined as the squared ratio of the amplitude of the reflected waves to the incident one. 38 We calculated also the normal displacement u z just behind the interface between the background and the PC. We show in Fig. 7 the field intensity (FI), i.e. the ratio between the amplitude of u z in the inner of the PC to that of the incident waves, at several frequencies with ff=0.503. The FI in Fig. 7 is always larger than one, meaning that, even in the absence of any focusing process, the rate of waves transmitted in the inhomogeneous PC gets enhanced in comparison to the incident waves. An increasing FI against frequency indicates that the displacement enhancement in the actual inhomogeneous PC becomes larger as the frequency increases. Such a behavior caused by the size of the air holes, together with the decreasing of the FWHM, leads to the increase of the MFI as the frequency increases. On the other hand, we show in Fig. 6(b) that, whatever the filling ratio is, the rate of reflection increases in a nonlinear way as the frequency increases. This behavior must be regarded against the slow decrease of the FMHM shown in Fig. 6(a) or the linear increase of FI depicted in Fig. 7. At some frequency, the influence of the reflectivity overcomes those combined of the FWHM and FI, and the MFI decreases as frequency increases. This is what happens around 9.5 MHz in Fig. 6(a). One should notice in Fig. 6(b), that the reflectivity at 11 MHz and ff=0.503 corresponds to an amplitude for the reflected u z about 22% of the incident waves. Following the above discussion, it is coherent to expect even stronger wave reflection at higher frequencies, and accordingly smaller MFI. Actually, at 13 MHz, we got 50% for the reflected amplitude with ff=0.503 and 4.35 for the MFI instead of 4.63 measured at 11 MHz.

B. Focusing behind the GRIN PC
Having analyzed the focalization of an elastic beam inside the GRIN PC, we present in this section the focalization out of the lens. To this end, we designed a GRIN PC featuring only 21 columns, instead of 55 as previously. In order to allow for comparisons, the geometrical parameters of the unit cell remain unchanged. The incident waves, as already described in the last section, enter the GRIN PC at normal incidence, and they are expected to focalize in the substrate behind the patterned area. This focusing scheme may be of great interest for applications, as for example exciting into resonance a resonator laying in the background or focusing elastic waves into a waveguide.
We show in Fig. 8 mappings of normalized u z at the moment when the displacement at the focus is maximum, recorded inside (x=[0, 21a]) as well as outside (x=[21a, 30a]) the GRIN PC, together with the calculated ray trajectories which are drawn as dashed lines. The analytical results suggest a focus zone at a position behind the GRIN PC, in good agreement with the numerical results. However, behind the GRIN PC, u z is much larger than the peak value measured inside the preceding lens at x=25a (see Sect. III A above). There are two main reasons for that: first, there is a strong mismatch between the impedances of the GRIN PC and background, especially at the central layer, where the filling ratio is maximum (see Fig. 6(b) and the following analysis). As a result, the elastic energy is partially reflected at the interface, which is not depicted by the ray trajectories. Secondly, the GRIN PC is actually an inhomogeneous structure: the filling ratio is much larger along the layer at the center than it is along the lateral layers. Accordingly, the displacement may be a little larger than expected. The position of the focusing zone is shifted to a shorter distance as compared with the focusing inside the GRIN PC, as shown in Fig. 3(a). Waves are refracted at the interface with the substrate, with a refractive angle larger than the incident angle, to form a focusing zone in the near field.
Almost similar conclusions can be drawn when the frequency is tuned to 8MHz, as shown in Fig. 8(b). In this case, the peak value appears at x=24a. This slight decrease of the focal distance must be ascribed to the enlarged focusing of incident wave through the GRIN PC, as it is the case for the focalization inside the structure. At 11 MHz, the peak value appears just behind the interface at x=22a for the same reason. The snapshots of u z shown in Figs. 8(a)-8(c), together with the corresponding ray trajectories, give clear evidence that the formalism applies in a broad band of working frequency. This provides a useful tool for studies dealing with the homogenization of PC's properties or of metamaterails. 36 We show in Fig. 8(d) the shift of the focal position against the frequency.
To further evaluate the focusing achievements of the lens, we determined the normalized transverse profiles of u z , as presented in Fig. 9 for 5 MHz (dashed line), 8MHz (dot dashed line) and 11MHz (solid line). The FWHMs are 4.49a, 3.07a and 1.97a respectively. This last value is smaller than its counterpart calculated in the case of inner lens focusing.
When normalized to the wavelength λ B in the background, the FWHMs are respectively 0.88λ B , 0.80λ B and 0.63λ B for 5, 8 and 11 MHz, with a down low trend against the frequency (see the inset in Fig. 9). In comparison to the focusing inside the GRIN PC (normalized to the wavelength in the PC), the FWHM at 5MHz is larger when focusing in the near field of the lens than on the inner. At 8MHz and 11MHz the focusing properties get even better than they were with the preceding system. It is of great interest to understand the physical reasons accounting for this result. First, the FWHM is normalized to the wavelength λ B instead of λ min for the reason that the focusing is formed in the background. The ratios between λ min and λ B are 0.86, 0.85 and 0.83 with increasing frequency, i.e. always less than one. As a result, the FWHM gets reduced. Secondly, the transverse wave number k y normalized to the maximum value k M is k M /k B = 0.66, 0.70, 0.73 at 5, 8, and 11 MHz respectively. An increasing ratio indicates enlarged Fourier components at the focus as the frequency increases, which can be quite helpful to enhance the resolution of focusing.
The MFI is 2.31, 2.31 and 2.23 at 5, 8 and 11 MHz respectively, as shown in Fig. 10(a). As compared to previous data for the inner focusing, the MFI here becomes smaller and are almost constant in the overall range of frequency. To illustrate this property, Fig. 10(a) represents the FWHM (triangular markers) whereas Fig. 10(b) shows the coefficient of transmission (defined as the squared ratio of the normal amplitude u z measured behind the PC to that of the incident waves 38 ) for three filling ratios, derived from the geometry drawn in Fig. 10(c): a PC composed of 21 air holes in an homogeneous silicon plate with all other conditions being the same as in Fig. 6(c). At 5 MHz, the focusing zone behind the GRIN lens features a FWHM (4.5a) larger than the one obtained for the inner focusing (3.6a). One should also notice that in this configuration the wave undergoes partial reflections on both the interfaces between the GRIN PC and the background and all things otherwise being equal, less elastic energy gets focalized. As a result, the normal amplitude u z is smaller for the focusing behind lens than it is in the inner of the lens. At 8 MHz and higher, the FWHM for the external focusing gets closer to the value obtained in the inner of the lens (about 3a in both cases). As a consequence of the lower transmission, the MFI for the external focusing (2.31) is much smaller than observed for the inner focusing (4.97). In this latter case, the transmission at 8 MHz, as high as 0.9, together with a FI of 1.66 (see Fig. 7), indicate that the effects induced by the patterning are at the origin of the differences between the MFI for external and for inner focusing.

IV. CONCLUSION
In this work, we have investigated the focalization of Lamb waves in GRIN PCs. The computation of the dispersion curves for different filling ratios allowed us to optimize the GRIN PC geometry and to demonstrate focalization of the elastic waves in the inner or behind the lens.
Although the index gradient was optimized for focusing at low frequency, the focalization was found to be even more efficient at higher frequencies (8 MHz and 11 MHz), close to the X point in the Brillouin zone. This is due to the fact that, for a given filling ratio, the refractive index strongly varies at higher frequencies. Then, we proposed a formalism allowing calculating the ray trajectory of Lamb wave in GRIN PCs with a square lattice. This analysis well explains the spreading of the elastic energy in the focal zone, resulting from spherical aberrations along the direction of propagation Nevertheless, the lens we have conceived leads to transverse profiles with FWHM close to 0.8λ min whatever the frequency is in the first band. The results show that the maximum amplitude inside the phononic lenses is about four times that of the source.
We have also investigated the possibility of focusing anti-symmetric Lamb waves behind the PC. We reduced the length of the lens to 21a which appeared to be the best compromise between efficient resolution and long focusing distance. Here, the focal distance calculated using a FEM method is in very good agreement with the ray trajectories analysis. The resolution at the focus was estimated to 0.8λ B at 8 MHz and ∼0.63λ B at 11MHz, which is closer to diffraction limit. Future developments of this work include the study of another lattice geometry like, for example, the honey-comb lattice, in order to increase the value of refractive index gradient and to compensate the spherical aberrations.

APPENDIX A
The tangent of angle ϕ is parallel to the slope of the ray trajectories, so that we used the ratio tan(ϕ)/tan(θ) in Fig. 4 to analyze the anisotropy along certain horizontal layers of the GRIN PCs. Using Eq. (2), the components of the wave vector k can be written as k x = k ΓX ( y) 1 − α ( y) cos 4θ 1 − α ( y) cos θ k y = k ΓX ( y) 1 − α ( y) cos 4θ 1 − α ( y) sin θ The partial derivatives of k x and k y with respect to θ are then derived from Eq. (A1) to be Eq. (A5) implies that tan(ϕ)/tan(θ) only depends on the anisotropy coefficient α at the input of lens.
To verify both the results in Fig. 4 and Eq. (A5), we found at 5 MHz the coefficient α along y=−8a and y=−2a, being −0.014 and 0.030 respectively. Eq. (A5) gives rise to 1.221 and 0.505 for the ratio tan(ϕ)/tan(θ) at the input of the lens, in excellent agreement with the results in Fig. 4. Similar comparisons can be made for other frequencies, and for each case we found excellent agreement.